This page details technical information useful for developers who wish to couple HEMCO
(the "Harmonized" Emissions Component) emissions component to other models.
The description of HEMCO
coupling to other models is available in Lin_et_al._2021
, which describes coupling to GEOS-Chem Classic, GCHP, WRF-GC, CESM2-GC
, and future NOAA models.
This work is made possible by a restructuring of HEMCO
, named HEMCO 3.0. HEMCO 3.0 separates model-specific components such as I/O, Regridding and the model speciation interface, into modular components, and isolate the HEMCO emissions Core.
This work is currently being actively worked on by the GEOS-Chem Support Team and Haipeng Lin (Harvard) as part of coupling GEOS-Chem with the CESM model.
- HEMCO Repository: geoschem/HEMCO on GitHub.
- Original description paper:
Keller_et_al._2014
. - Coupling and HEMCO 3.0 description paper:
Lin_et_al._2021
. - The HEMCO User's Guide
- HEMCO versions
As part of the HEMCO 3.0
restructuring, "HEMCO" is now divided into three pieces depending on their function:
- The HEMCO Core. Emissions calculations logic, containers, data types, etc.
- Data Input Layer. I/O (previously
HCOIO_Read/Write_*_Mod
), Regridding (HCO_MESSY_REGRID
,HCO_INTERP_MOD
), ... This will be rearranged intoRegrid/
andIO/
folders in a future version. Right now due to dependencies, some of these files still live in theCore/
folder. - Model Interface Layer. Code that couples
HEMCO
with other models. There are common utilities available atInterfaces/HCO_Interface_Common.F90
.
Note
Note that not all code pertinent to model coupling actually lives inside of HEMCO
; this is by design, as data types that are external to HEMCO
(i.e. GEOS-Chem types such as State_Met
, CESM types such as physics_state
, WRF types such as domain
) must be maintained with the model and not inside HEMCO. Some code lives in Interfaces/
, and some will live inside the model.
TBD
In order to interface HEMCO
with the target model, there are a few primary tasks that need to be performed as outlined below.
Data/code that needs to be provided to HEMCO
based on the target model's data structures include:
- The clock and time-step of the target model
- List of species and physical properties (molecular weight required; other properties such as Henry's law constants are optional, only for extensions such as SeaFlux)
- Grid information (
I
,J
,L
atmospheric '0-D box' dimensions required; if using HEMCO built-in regrid, then specifics are needed. See below)
Data/code that needs to be retrieved from HEMCO into the target model's data structures (i.e. state object for constituent flux/concentrations) include:
- Emissions fluxes (kg/m2/s format) retrieved from HEMCO, aggregated per species ID, for current time step
- Other data retrieved from HEMCO (using
HCO_GetPtr
orHCO_EvalFld
)
Important
Avoid calling HEMCO functions directly from outside of a specific module designed to interface HEMCO with the model. This is so the interface can be updated more easily if subroutines within HEMCO such as HCO_GetPtr
change, and the HEMCO state (:code`HcoState`) doesn't need to be passed to everywhere in your model that needs to retrieve data from HEMCO. It is also useful so regridding to/from HEMCO can be performed in a centralized location, if so needed by the model. For example, GEOS-Chem wraps HCO_GetPtr
and HCO_EvalFld
into its own interface, HCO_GC_GetPtr
, HCO_GC_EvalFld
, which will auto-magically add the HcoState
argument, in addition to handling regridding if necessary.
Things that come out-of-the-box and generally do not require customization to a specific model:
- Reading configuration file (
HEMCO_Config.rc
), although the path needs to be specified - HEMCO "driver" (run) routines
- Managing HEMCO memory (initializing HEMCO state in
HcoState
, extensions state inExtState
, etc.)
This is a three-step process. First initialize the configuration object (HcoConfig
):
call ConfigInit(HcoConfig, HMRC, nModelSpecies=nSpc)
You have to register the species first in addition to some other HcoConfig properties:
HcoConfig%amIRoot = masterproc
HcoConfig%MetField = 'MERRA2'
HcoConfig%GridRes = ''
HcoConfig%nModelSpc = nHcoSpc
HcoConfig%nModelAdv = nHcoSpc ! # of adv spc?
do N = 1, nHcoSpc
HcoConfig%ModelSpc(N)%ModID = N ! model id
HcoConfig%ModelSpc(N)%SpcName = trim(solsym(N))
enddo
Then open the configuration file in two phases; after phase 1, initialize the log file on the MPI root process:
call Config_ReadFile(HcoConfig%amIRoot, HcoConfig, HcoConfigFile, 1, HMRC, IsDryRun=.false.)
! Open the log file
if(masterproc) then
call HCO_LOGFILE_OPEN(HcoConfig%Err, RC=HMRC)
endif
call Config_ReadFile(HcoConfig%amIRoot, HcoConfig, HcoConfigFile, 2, HMRC, IsDryRun=.false.)
Warning
Note that the species count has to be populated three times. Once above at ConfigInit
, and twice inside the initialized HEMCO Config object.
Some species physical properties need to be defined for HEMCO
extensions, such as molecular weight and henry's law constants:
!-----------------------------------------------------------------------
! Register HEMCO species information (HEMCO state object)
!-----------------------------------------------------------------------
do N = 1, nHcoSpc
HcoState%Spc(N)%ModID = N ! model id
HcoState%Spc(N)%SpcName = trim(solsym(N)) ! species name
HcoState%Spc(N)%MW_g = adv_mass(N) ! mol. weight [g/mol]
! HcoState%Spc(N)%HenryK0 ! [M/atm]
! HcoState%Spc(N)%HenryCR ! [K]
! HcoState%Spc(N)%HenryPKA ! [1]
enddo
Note
If you are not using HEMCO extensions, only ModID
, SpcName
and MW_g
need to be defined.
HcoState%NX = my_IM
HcoState%NY = my_JM
HcoState%NZ = LM
There are many ways of defining the vertical discretization. Check HCO_VertGrid_Define
.
! Pass Ap, Bp values, units [Pa], [unitless]
call HCO_VertGrid_Define(HcoState%Config, &
zGrid = HcoState%Grid%zGrid, &
nz = HcoState%NZ, &
Ap = Ap, &
Bp = Bp, &
RC = HMRC)
Note
HEMCO requires HORIZONTAL grid information only if it is using internal regridding routines, i.e. MAP_A2A
or MESSy. Otherwise, this can be filled with dummy information.
Warning
If HEMCO
internal regridding (MAP_A2A
) regridding routines are used, only rectilinear grids are supported.
This is because XMid
, YMid
, ... arrays are 1-dimensional and thus curvilinear coordinates cannot be stored. The underlying MAP_A2A
algorithm can handle curvilinear; it is just due to the data structure. This will be fixed in a future HEMCO version.
! Point to grid variables
HcoState%Grid%XMID%Val => XMid (my_IS:my_IE , my_JS:my_JE )
HcoState%Grid%YMID%Val => YMid (my_IS:my_IE , my_JS:my_JE )
HcoState%Grid%XEdge%Val => XEdge (my_IS:my_IE+1, my_JS:my_JE )
HcoState%Grid%YEdge%Val => YEdge (my_IS:my_IE , my_JS:my_JE+1)
HcoState%Grid%YSin%Val => YSin (my_IS:my_IE , my_JS:my_JE+1)
HcoState%Grid%AREA_M2%Val => AREA_M2(my_IS:my_IE , my_JS:my_JE )
Here we point HEMCO
's variables to structures we have created in the model. Examples in how to create these structures are available in the HEMCO-CESM interface.
An example to translate and define meteorological quantities such as temperature, humidity, etc. is available in the HEMCO-CESM interface.
Prerequisites:
! HEMCO
use HCO_Interface_Common, only: GetHcoVal, GetHcoDiagn
use HCO_Clock_Mod, only: HcoClock_Set, HcoClock_Get
use HCO_Clock_Mod, only: HcoClock_EmissionsDone
use HCO_Diagn_Mod, only: HcoDiagn_AutoUpdate
use HCO_Driver_Mod, only: HCO_Run
use HCO_EmisList_Mod, only: Hco_GetPtr
use HCO_FluxArr_Mod, only: HCO_FluxArrReset
use HCO_GeoTools_Mod, only: HCO_CalcVertGrid, HCO_SetPBLm
Also make sure the time steps are set correctly. Use from the common utilities:
call HCOClock_Set(HcoState, year, month, day, &
hour, minute, second, IsEmisTime=.true., RC=HMRC)
call HCO_FluxArrReset(HcoState, HMRC)
HEMCO
needs an updated vertical grid at each time step. Data passed into HCO_CalcVertGrid
can vary and the definition can be checked for acceptable parameters.
call HCO_CalcVertGrid(HcoState, PSFC, ZSFC, TK, BXHEIGHT, PEDGE, HMRC)
call HCO_SetPBLm(HcoState, PBLM=State_HCO_PBLH, &
DefVal=1000.0_hp, & ! default value
RC=HMRC)
To document.
! Range of species and emission categories.
! Set Extension number ExtNr to 0, indicating that the core
! module shall be executed.
HcoState%Options%SpcMin = 1
HcoState%Options%SpcMax = -1
HcoState%Options%CatMin = 1
HcoState%Options%CatMax = -1
HcoState%Options%ExtNr = 0
! Use temporary array?
HcoState%Options%FillBuffer = .FALSE.
call HCO_Run( HcoState, 1, HMRC, IsEndStep=.false. )
call HCO_Run( HcoState, 2, HMRC, IsEndStep=.false. )
Necessary only if you are using HEMCO
extensions.
call HCOX_Run(HcoState, ExtState, HMRC)
!-----------------------------------------------------------------------
! Update "autofill" diagnostics.
! Update all 'AutoFill' diagnostics. This makes sure that all
! diagnostics fields with the 'AutoFill' flag are up-to-date. The
! AutoFill flag is specified when creating a diagnostics container
! (Diagn_Create).
!-----------------------------------------------------------------------
call HcoDiagn_AutoUpdate(HcoState, HMRC)
!-----------------------------------------------------------------------
! Tell HEMCO we are done for this timestep...
!-----------------------------------------------------------------------
call HcoClock_EmissionsDone(HcoState%Clock, HMRC)
You can either use the common utilities, where data is retrieved using GetHcoValEmis
, or tap into the arrays directly.
For generic data containers, pass the container name like so:
! For grabbing data from HEMCO Ptrs (uses HEMCO single-precision)
real(sp), pointer :: Ptr2D(:,:)
real(sp), pointer :: Ptr3D(:,:,:)
logical :: FND
call HCO_GetPtr(HcoState, 'CONTAINER_NAME', Ptr2D, HMRC, FOUND=FND)
Warning
Important: Note that deposition (sink terms) fluxes are handled separately from emissions in HEMCO. This is particularly important if you use HEMCO to calculate deposition terms, e.g. the sink term in SeaFlux
(sea-air exchange). The standard in HEMCO is that the sink terms are stored as deposition velocities (depv
, unit 1/s
) so HEMCO generally does not need to be aware of concentrations.
A thorough discussion of this is in the HEMCO GitHub issue tracker. The code to handle deposition velocities from HEMCO is generally as follows:
!------------------------------------------------------------------
! Also add drydep frequencies calculated by HEMCO (e.g. from the
! air-sea exchange module) to DFLX. These values are stored
! in 1/s. They are added in the same manner as the drydep freq values
! from drydep_mod.F90. DFLX will be converted to kg/m2/s later.
! (ckeller, 04/01/2014)
!------------------------------------------------------------------
CALL GetHcoValDep( NA, I, J, L, found, dep )
IF ( found ) THEN
dflx(I,J,NA) = dflx(I,J,NA) &
+ ( dep * spc(I,J,NA) / (AIRMW / ThisSpc%MW_g) )
ENDIF