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:

  1. 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 the chem/ folder and make 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)

  1. 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)

! ... 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 Ptr3D is in GEOS-Chem indices. Use (II, JJ, k) for indexing. Inside the loops already provided in wrfgc_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.rc to 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