Coupling HEMCO with other models

From Geos-chem
Jump to: navigation, search

This page details technical information useful for developers who wish to couple the HEMCO ("Harmonized" Emissions Component) emissions component to other models.

Overview

This work will be made possible by a restructuring of HEMCO, named HEMCO 3.0. We will separate 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.

Useful Resources

Terminology

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 into Regrid/ and IO/ folders in a future version. Right now due to dependencies, some of these files still live in the Core/.
  • Model Interface Layer. Code that couples HEMCO with other models. There are common utilities available at Interfaces/HCO_Interface_Common.F90. **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.

Technical Notes (Data Input Layer)

TBD

Technical Notes (Model Interface Layer)

HEMCO 3.0 Model Interface Layer Overview

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 or HCO_EvalFld)

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 in ExtState, etc.)

Reading the HEMCO configuration file and defining species list

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

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

If you are not using HEMCO extensions, only ModID, SpcName and MW_g need to be defined.

Defining Grid

Define atmospheric column numbers

HcoState%NX = my_IM
HcoState%NY = my_JM
HcoState%NZ = LM

Define the vertical grid

There are many ways of defining the vertical discretization. Check HCO_VertGrid_Define.

! Pass Ap, Bp values, units [Pa], [unitless]
! later remove masterproc
call HCO_VertGrid_Define(HcoState%Config,                &
                         zGrid = HcoState%Grid%zGrid,    &
                         nz    = HcoState%NZ,            &
                         Ap    = Ap,                     &
                         Bp    = Bp,                     &
                         RC    = HMRC)

Define horizontal grid parameters

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 (3.1+)

! 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  )

Defining Met Fields for HEMCO Extensions

Running HEMCO

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

Update the HEMCO clock

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)


Reset fluxes for new timestep

call HCO_FluxArrReset(HcoState, HMRC)

Update vertical grid parameters

HEMCO needs an updated vertical grid at each time step. Data passed into HCO_CalcVertGrid can vary. Check the code.

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)

Some dummy setup (advanced)

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.

Run HEMCO driver

call HCO_Run( HcoState, 1, HMRC, IsEndStep=.false. )
call HCO_Run( HcoState, 2, HMRC, IsEndStep=.false. )

Run HEMCO extensions driver

call HCOX_Run(HcoState, ExtState, HMRC)

Close timestep

!-----------------------------------------------------------------------
! 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)

Retrieving emissions data from HEMCO

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)

Retrieving deposition velocities (depv) from HEMCO

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 https://github.com/geoschem/HEMCO/issues/72#issuecomment-789409266. 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