User:Jmandel/chem

See Coupling with WRF-Chem for the users' guide.

Tracers

Tracers are compiled in only when WRF_CHEM is active. See dyn_em/start_em.F#L1173 dyn_em/solve_em.F#L2421

They are initialized in chem/module_input_tracer.F see particularly line 58 for the different initializations

Jon has discovered that tracer_opt=1 creates one tracer called smoke. I was unable to find any documentation for that or in fact any documentation about tracer_opt in WRF users guide whatsoever. I think this is TRACER_SMOKE as mentioned in chem/module_input_tracer.F line 15. Note that chem/module_input_tracer.F has #ifdef WRF_CHEM which seem to be useless because it is never compiled or accessed otherwise. Also note if tracer_opt == TRACER_SMOKE then the whole array chem is initialized to a nonzero constant (which seems to indicate it is useless)

I had changed Jon's code for inserting the tracer to account for the density of the atmosphere and the height of the first layer. It contains heat flux as a proxy for the mass burned (in the code, the heat/mass burned is constant). In future I can treat it similarly as the emissions in the chem array if you like (chem/module_input_tracer.F says it should be CO). I think to use it you run chem_opt=0 and tracer_opt=1.

Tracer accessed as grid%tracer(i,k,j,p_smoke) real,intent(inout),dimension( ims:ime,kms:kme,jms:jme,num_tracer ) :: tracer

from use module_state_description , only: num_tracer, p_smoke

which has INTEGER, PARAMETER :: tracer_smoke = 1

similar code in chem/module_add_emiss_burn.F: chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke)+ebu(i,k,j,p_ebu_co)*conv_rho

defined in ./inc/scalar_indices.inc: P_smoke = 1 ; F_smoke = .FALSE.

Tracers without WRF-Chem

The tracer_opt=2 seems to work with WRF-Chem only and mimics biomass burning emissions. The option that seems to work without WRF-Chem is tracer_opt=1 which turns on 8 tracers: tr17_1 ... tr17_8 We could use these tracers to passively resolve concentrations of the basic fire-emitted species, using existing emission factors for RADM2. For instance:

tr17_1 = co  
tr17_2 = no  
tr17_3 = no2  
tr17_4 = pm25 (= emissions factor of pm25i + pm25j or just p25 from mozart)
tr17_5 = oc  (= emission factor of oc1 + oc2)
tr17_6 = bc  (= emission factor of bc1 + bc2)
tr17_7 = NMHC
tr17_8 = NMOC

That would allow us to carry these variables for use by an external chemical model like CMAQ, which is used with BlueSky, or just to use a tracer field directly, to define plume properties, without engaging any component of WRF-Chem.

I don't know though, if these tracers get initialized in a similar way as chem variables used now and how tough would it be to do something like that.

Chem variables and conversion

from registry.chem: package radm2sorg chem_opt==2 - chem:so2,sulf,no2,no,o3,hno3,h2o2,ald,hcho,op1,op2,paa,ora1,ora2,nh3,n2o5,no3,pan,hc3,hc5,hc8,eth,co,ol2,olt,oli,tol,xyl,aco3,tpan,hono,hno4,ket,gly,mgly,dcb,onit,csl,iso,hcl,ho,ho2,so4aj,so4ai,nh4aj,nh4ai,no3aj,no3ai,naaj,naai,claj,clai,orgaro1j,orgaro1i,orgaro2j,orgaro2i,orgalk1j,orgalk1i,orgole1j,orgole1i,orgba1j,orgba1i,orgba2j,orgba2i,orgba3j,orgba3i,orgba4j,orgba4i,orgpaj,orgpai,ecj,eci,p25j,p25i,antha,seas,soila,nu0,ac0,corn

state real so2 ikjftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "so2" "SO2 concentration" "ppmv"

Call chain: chem_driver.F -> emissions_driver.F -> add_emiss_burn.F, plumerise_driver -> plumerise module_emissions_anthropogenics

gas emissions are converted to the concentration changes in module_emissions_anthropogenics.F located in WRFV3/chem. Basically, since each gas occupies the same volume under the same P and T, all gass species are converted from mole/km2h to delta ppmv (which is 1E6 times mixing ratio in mole/mole) in the same way.

The conversion looks like that:

                    E    is emission [mole/(km^2 hr)] 
     conv_rho    is converted emission in ppmv [-]
    rho_phy(i,jk)    is air density [kg/m3}
   dz8w(i,j,k)    is vertical grid spacing for the lowest model layer [m] 
         dtstep     is model time step [s]

8.047e-6 is conversion factor – km^2 hr to m^2 s (2.7778e-10) * molecular mass of air 0.02897 * 1e6 (mol/mol -> part per million) = 8.047e-6

conv_rho= E*dtstep*8.047E-6/(rho_phy(I,k,J)*dz8w(i,k,j))

Both rho_phy(I,k,J) and dz8w(i,k,j) are WRF variables responding to the change in the pressure and temperature so the plume expansion is taken care of automatically.

So, the code should be

    ! on the fire mesh
    !
    ! emission (g/m^2) = initial fuel load (kg/m^2) * fraction burned (1) * emission of the chemical species from the given fuel (g/kg) 
    emis_fgrid(i,j,p_species) = fgip(i,j) * fuel_frac_burnt(i,j) * emis_finn(ifuel(i,j),p_species)  
    ! 
    ! average to the atmospheric mesh to get emis_grid(i,j,p_species) (g/m^2)
    !
    ! on the atmospheric mesh
    !
    ! mols of air in the 1st layer (mol/m^2) =  (air density (kg/m^3) * layer height (m)) /  molecular mass of air (kg/mol)
    air_mols = rho_phy(i,kts,j)*dz8w(i,kts,j) / 28.97e-3
    !
    ! emissions (mols/m^2) = emissions (g/m^2) / molecular mass of the species (g/mol) 
    emis_mols = emis_grid(i,j,p_species) / mol(p_species)  ! the constants mols(p_species) are probably in chem somewhere??
    !
    !  add the concentration in ppmv to chem
    chem(i,kts,j,p_species) = chem(i,kts,j,p_species) + 1e-6 * emis_mols / air_mols

or composed:

    inv_mol = 1/mol(p_species)
    conv_fact = 28.97e-9 / (rho_phy(i,kts,j)*dz8w(i,kts,j))
    chem(i,kts,j,p_species) = chem(i,kts,j,p_species) + conv_fact * inv_mol