diff --git a/src/idl/dustem_plugin_ff_hii.pro b/src/idl/dustem_plugin_ff_hii.pro
new file mode 100755
index 0000000..255ae25
--- /dev/null
+++ b/src/idl/dustem_plugin_ff_hii.pro
@@ -0,0 +1,137 @@
+FUNCTION dustem_plugin_ff_hii,key=key,val=val,scope=scope,paramtag=paramtag,paramdefault=paramdefault,help=help
+
+;+
+; NAME:
+;    dustem_plugin_ff_hii
+; PURPOSE:
+;    DustEMWrap plugin to compute free-free emission
+; CATEGORY:
+;    DustEM, Distributed, Mid-Level, Plugin
+; CALLING SEQUENCE:
+;    freefree=dustem_plugin_freefree([,key=][,val=])
+; INPUTS:
+;    None
+; OPTIONAL INPUT PARAMETERS:
+;    key  = input parameter number
+;    val  = input parameter value
+; OUTPUTS:
+;    freefree = free-free spectrum (on dustem wavelengths)
+; OPTIONAL OUTPUT PARAMETERS:
+;    scope = scope of the plugin
+;    paramdefault = default values of parameters
+;    paramtag = plugin parameter names as strings
+; ACCEPTED KEY-WORDS:
+;    help                  = if set, print this help
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The DustEMWrap IDL code must be installed
+; PROCEDURE:
+;    This is a DustEMWrap plugin
+; EXAMPLES
+;    dustem_init
+;    vec=dustem_plugin_freefree(scope=scope)
+; MODIFICATION HISTORY:
+;    Written by JPB 2022 
+;    Evolution details on the DustEMWrap gitlab.
+;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
+;-
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_plugin_ff_hi'
+  output=0.
+  goto,the_end
+ENDIF
+
+;default values of input parameters
+Te=1.0E4             ;default electron temperature (K)
+;Adding an amplitude just for the sake of it (still not done)
+;Amplitude=       ;Amplitude for Te=1.E4 and nu=  GHz (arb) and lambda=1E4 microns
+;Not so sure about this emission measure
+EM = 100.             ;default emission measure  (cm^-6 pc)
+smallp=0.0            ;default polarization fraction
+psi=0.                ;default polarization angle
+scope='ADD_SED'
+paramtag=[textoidl('T_{e}')+' [K]','EM','p','Psi [deg]']
+paramdefault=[Te,EM,smallp,psi]
+IF keyword_set(key) THEN BEGIN 
+  a=where(key EQ 1,count1)
+  b=where(key EQ 2,count2)
+  c=where(key EQ 3,count3) ;default polarization fraction -newly added
+  d=where(key EQ 4,count4) ;default polarization angle -newly added
+  IF count1 NE 0 then Te=(val(a))(0)
+  IF count2 NE 0 then EM=(val(b))(0)
+  IF count3 NE 0 then smallp=val[c[0]]   
+  IF count4 NE 0 then psi=val[d[0]]    
+ENDIF
+
+lambir=dustem_get_wavelengths()
+Nwavs=n_elements(lambir)
+cmic=3.e14
+nu=cmic/lambir   ;Hz
+k_B = 1.380649E-23 ;J.K^-1 
+g_ff = alog(4.955*1D-2/nu*1D9)+1.5*alog(Te) ;gaunt factor
+tau_ff = 3.014*1D-2*Te^(-1.5)*(nu*1D-9)^(-2)*EM*g_ff ;ff optical depth 
+T_ff = Te*(1-exp(-tau_ff))
+
+out=fltarr(Nwavs,3)
+out[*,0] = 2*k_B*T_ff*nu^2/(cmic*1E-6)^2*1E20 ;convert to MJy by multiplying by 1E20
+
+polar_ippsi2iqu,out[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
+
+out[*,1]=Q
+out[*,2]=U
+
+; ;Not so sure about these two: 
+; mjy=1           ;output is in MJy/sr
+; lambir_ref=10000.
+; 
+; 
+; output=fltarr(Nwavs,3) ;newly added
+; 
+; use_method='Dickinson2003_norm'
+; 
+; CASE use_method OF
+;     'WallsGabaud1998':BEGIN
+; 		em=1.   ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
+; 		output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
+; 		fact=Amplitude/I_halpha_R
+;         ;=== normalize to the requested Halpha value.
+; 		output[*,0]=output[*,0]*fact
+; 	END
+;     'WallsGabaud1998_scaled':BEGIN
+; 		em=1.   ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
+; 		output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
+; 		norm=interpol(output[*,0],lambir,lambir_ref)
+; 		output[*,0]=output[*,0]/norm*Amplitude
+; 	END
+; 	'Deschenes2008':BEGIN
+; 		beta_freefree=2.+1./(10.48+1.5*alog(Tgas/8.e3)-alog(nu/1.e9))
+; 		output[*,0]=nu^(-1.*beta_freefree)
+; 		norm=interpol(output[*,0],lambir,lambir_ref)
+; 		output[*,0]=output[*,0]/norm*Amplitude
+; 	END
+; 	'Dickinson2003_norm':BEGIN
+; 		;stop
+; 		em=1.
+; 		T4=Tgas/1.e4
+; 		nu_ghz=nu/1.e9
+; 		a=0.366*(nu_ghz/1.e9)^(-0.15)*(alog(4.995*1e-2*(nu_ghz)^(-1.))+1.5*alog(Tgas))
+; 		Tb=8.369e3*a*(nu_ghz)^(-2.)*T4^0.667*10^(0.029*T4)*(1.+0.08)  ;in mK for 1 Rayleigh
+; 		convert_mk_mjy, lambir, Tb, I_Mjy, /RJ
+; 		norm=interpol(I_Mjy,lambir,lambir_ref)
+; 		output[*,0]=I_Mjy/norm*Amplitude
+; 	END
+; ENDCASE
+
+; polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
+; 
+; output[*,1]=Q
+; output[*,2]=U
+
+the_end:
+RETURN,out
+  
+END
--
libgit2 0.21.2