We’ve been receiving lots of questions from WRF-GC users regarding whether it’s possible for WRF-GC to output custom quantities, such as internal diagnostic quantities in GEOS-Chem or HEMCO.
The answer is, maybe. It is important to remember that, as described in the model description paper WRF-GC does not output in GEOS-Chem HISTORY or HEMCO diagnostic format. So HEMCO_Diagn.rc’s output options, or the DIAGxx legacy bpch diagnostics in input.geos are all ignored. But not all hope is lost.
Quick rundown
WRF-GC history is controlled by a few different components:
- The WRF registry. In order for a field to be recognizable by the WRF framework, it needs to be declared in the WRF registry. The easiest way to do this is through editing
chem/registry.chem, then installing it by going to thechem/folder andmake install_registry.
Examples are built-in:
# WRF-GC dummy diagnostics
state real diagikj0 ikj misc 1 - rh "DiagIKJ-0" "Dummy 3-D diagnostic ikj,0" "TBD"
state real diagikj1 ikj misc 1 - rh "DiagIKJ-1" "Dummy 3-D diagnostic ikj,1" "TBD"
# ...
state real diagij0 ij misc 1 - rh "DiagIJ-0" "Dummy 2-D diagnostic ikj,0" "TBD"
state real diagij1 ij misc 1 - rh "DiagIJ-1" "Dummy 2-D diagnostic ikj,1" "TBD"
# ...
Read more about the WRF registry here (external link). In short, the diagikjX are samples for 3-D output variables, and diagijX are samples for 2-D output variables.
If you want to output something extra into the wrfout files, you need to add entries (or reuse the old ones - see below)
Note, if you modify the registry, you have to ./clean -a and rebuild the model (./configure -hyb, ./compile em_real)
- The WRF-GC coupler. Only the WRF-GC coupler can bridge data between WRF and GEOS-Chem. We actually use this to debug some “diagnostic” quantities within WRF-GC. Go to
wrfgc_convert_state_mod.F:
! The following code simply updates WRF diagnostics data for testing...
grid%diagikj0(i, k, j) = State_Met%CMFMC(II, JJ, k)
grid%diagikj1(i, k, j) = State_Met%PEDGE(II, JJ, k)
grid%diagikj2(i, k, j) = State_Met%PFLCU(II, JJ, k)
grid%diagikj3(i, k, j) = State_Met%TAUCLW(II, JJ, k)
grid%diagikj4(i, k, j) = State_Met%TAUCLI(II, JJ, k)
grid%diagikj5(i, k, j) = State_Met%AD(II, JJ, k)
grid%diagikj6(i, k, j) = State_Met%AIRVOL(II, JJ, k)
grid%diagij0(i, j) = State_Met%PARDR(II, JJ)
grid%diagij1(i, j) = State_Met%PARDF(II, JJ)
grid%diagij2(i, j) = State_Met%CONV_DEPTH(II, JJ)
grid%diagij3(i, j) = State_Met%FLASH_DENS(II, JJ)
If you want to output a 2-D data structure in GEOS-Chem, then put it on the right hand side of diagij variables; for 3-D data, use the diagikj containers. If there aren’t enough, add more through the registry. But you’re free to reuse the existing 7+4 containers for convenience.
How to access GEOS-Chem internal data structures?
HEMCO example
Let’s take a closer look at this line of code:
grid%diagikj0(i, k, j) = State_Met%CMFMC(II, JJ, k)
The left-hand side is in the WRF world, and right hand side is the GEOS-Chem world. Always use (i, j, k) indices in the WRF world, and (II, JJ, k) indices in the GEOS-Chem world. Relative or absolute indices that you specify using numbers in the horizontal dimensions I, J will not work and crash the model. WRF-GC runs in a distributed MPI parallelization environment, and the indices in WRF do not correspond to the GEOS-Chem indices, and each CPU can only see a chunk of the entire domain. Trust me, I wrote the code.
If you want to output data from State_Met or State_Chm, the example above can be easily edited. If you want to output data from HEMCO, you have to call data from HEMCO like so (the below is pseudo-code for conceptual purposes; you have to edit it a bit)
! ... add use statement
USE hco_interface_mod, ONLY: GetHcoDiagn
! ... create an array temporary
REAL(f4), POINTER :: Ptr3D(:,:,:)
! get the data from hemco
CALL GetHcoDiagn( 'BCPI_ANTH', .FALSE., ERR, Ptr3D=Ptr3D )
! ... inside the II, JJ loop
do j = jts, jte
do i = its, ite
! ...
! add this code:
if ( associated( Ptr3D ) ) then
grid%diagikj1(i, k, j) = Ptr3D(II, JJ, k)
endif
! ...
enddo
enddo
! add this for cleanup:
Ptr3D => NULL()
A few things to remember:
- HEMCO will pass you a pointer to the GEOS-Chem world, thus
Ptr3Dis in GEOS-Chem indices. Use(II, JJ, k)for indexing. Inside the loops already provided inwrfgc_convert_state_mod.F, there is loop index translation provided for you:
do j = jts, jte
do i = its, ite
II = i - its + 1
JJ = j - jts + 1
! ...
For the 3-D, use the 3-D loop:
do k = kts, kte
do j = jts, jte
do i = its, ite
II = i - its + 1
JJ = j - jts + 1
- Always clean-up after your pointers, and check if pointers are present before reading from them. HEMCO diagnostics, by design, do not always have to exist. You have to modify
HEMCO_Diagn.rcto add them.
GEOS-Chem “netCDF diagnostics” example
Refer to the example code in WRFGC_Set_WRF:
if(config_flags%gc_diagn_spc_n0 .ne. 0 .and. config_flags%gc_diagn_spc_n0 .le. State_Chm%nSpecies) then
grid%cldcnvflx_n0(i, k, j) = State_Diag%CloudConvFlux(II, JJ, k, config_flags%gc_diagn_spc_n0)
! Budget Diagnostics
! EmisDryDep uses nAdvect map
grid%gcemisdrydep_full_n0(i, j) = State_Diag%BudgetEmisDryDepFull(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcemisdrydep_trop_n0(i, j) = State_Diag%BudgetEmisDryDepTrop(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcemisdrydep_pbl_n0(i, j) = State_Diag%BudgetEmisDryDepPBL(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
! Mixing uses nAdvect map
grid%gcmixing_full_n0(i, j) = State_Diag%BudgetMixingFull(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcmixing_trop_n0(i, j) = State_Diag%BudgetMixingTrop(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcmixing_pbl_n0(i, j) = State_Diag%BudgetMixingPBL(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
! Convection uses nAdvect map
grid%gcconv_full_n0(i, j) = State_Diag%BudgetConvectionFull(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcconv_trop_n0(i, j) = State_Diag%BudgetConvectionTrop(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcconv_pbl_n0(i, j) = State_Diag%BudgetConvectionPBL(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
! Chemistry uses nAdvect map
grid%gcchem_full_n0(i, j) = State_Diag%BudgetChemistryFull(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcchem_trop_n0(i, j) = State_Diag%BudgetChemistryTrop(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
grid%gcchem_pbl_n0(i, j) = State_Diag%BudgetChemistryPBL(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%AdvectID)
! Wet Deposition uses nWetDep map
if(State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId .ne. MISSING_INT) then
grid%gcwetdep_full_n0(i, j) = State_Diag%BudgetWetDepFull(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
grid%gcwetdep_trop_n0(i, j) = State_Diag%BudgetWetDepTrop(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
grid%gcwetdep_pbl_n0(i, j) = State_Diag%BudgetWetDepPBL(II, JJ, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
grid%wetlscnvfrc_n0(i, k, j) = State_Diag%WetLossConvFrac(II, JJ, k, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
grid%wetlscnv_n0(i, k, j) = State_Diag%WetLossConv(II, JJ, k, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
grid%wetlossls_n0(i, k, j) = State_Diag%WetLossLS(II, JJ, k, &
State_Chm%SpcData(config_flags%gc_diagn_spc_n0)%Info%WetDepId)
endif
endif