Adding new dummy species with loss rate dependent on meteorology

[cschooling]

Name and Institution (Required)

Name: Chlöe Schooling Institution: University of Edinburgh

Confirm you have reviewed the following documentation

• [x] Support guidelines
• [x] User manuals
• [x] Debugging GEOS-Chem and HEMCO errors
• [x] Current and past Github issues

Description of your issue or question

I would like to introduce a species (probably for the CO2 simulation run). This species will be NOx, however I am not wanting to include the NOx chemistry. Rather I have parameterised the loss and product rates of NOx as a function of some meteorological features (e.g. temperature, solar zenith angle, lat/long/level coordinates, surface incident shortwave flux, NOx concentration, O3 concentration, CO concentration). Would it be possible to include this as a dummy species, where the loss and product rates are calculated based on these features at each time point and grid point.

If so could you please give me an idea of how it would be best to go about introducing this into the model.

Thanks so much!

[yantosca]

Thanks for writing @cschooling. It might be easier to do this in the carbon simulation (which is KPP-based) than in the CO2 simulation (which is not). The carbon simulation has CH4-CO2-CO-OCS species, but you can also run with a single species if you wish. You could add a new species to the simulation and add a new rate-law function in GEOS-Chem/KPP/carbon/carbon_Funcs.F90. You would then probably have to add NOx and O3 as dummy species that are read in via HEMCO (see GeosCore/carbon_gases_mod.F90 for how the existing code is set up). You can get the pointers from HEMCO via calls to HCO_EvalFld.

[cschooling]

Hi @yantosca, thank you for this advice! I am initially attempting to run the carbon simulation for v 14.2.1 (without any changes). However, I have some issues, where I am receiving the following error message:

GEOS-Chem ERROR: Cannot get pointer to HEMCO field CH4_OIL
 -> at Emiss_Carbon_Gases (in module GeosCore/carbon_gases_mod.F90)
===============================================================================
===============================================================================
GEOS-Chem ERROR: Error encountered in "Emiss_Carbon_Gases"!
 -> at Emissions_Run (in module GeosCore/emissions_mod.F90)
===============================================================================
===============================================================================
GEOS-CHEM ERROR: Error encountered in "Emissions_Run"! after drydep!
STOP at  -> at GEOS-Chem (in GeosCore/main.F90)
===============================================================================

Additionally, I do not have the CH4_LOSS datafile ($ROOT/CH4/v2023-04/GC_CH4_LOSS/GCC14_72LM.ch4loss.4x5.nc4), which I can not find available on the database http://geoschemdata.wustl.edu. Perhaps this is contributing to these issues, but for now I have just turned off this inventory in HEMCO.

My run directory is attached, the GC source code has not been changed.

carbon_CO2.zip

[yantosca]

Hi @cschooling, thanks for the feedback. You should be able to get the missing data from the http://ftp.as.harvard.edu/gcgrid/data/ExtData/HEMCO/CH4/v2023-04/GC_CH4_LOSS/ folder.

There should have been a fix in 14.2.1 to prevent the error that you describe. Did you take the 14.2.1 release, or was it a pre-release version? See PR #1963

[cschooling]

Thanks @yantosca! I took the 14.2.1 release (downloaded last week), and also have tried with 14.2.2. However, I still get this error when running with a single species - fine when run with all 4.

In regards to introducing the NOx as an additional species here. I have introduced NOx, and two dummy placeholders to represent the NOx reactants/products within carbon.eqn:

NOx                 = IGNORE;  { Active NOx species                                    }
DummytoNOx   = IGNORE;  { Dummy placeholder for NOx reactant       }
DummyfrNOx    = IGNORE;  { Dummy placeholder for NOx product       }

With reactions defined as:

DummytoNOx    = NOx                 : k_Trop(4);
NOx           =     DummyfrNOx          : k_Trop(5);

I have initially set these reaction rates of these reactions as constants, within carbon_Funcs.F90, to get it running: k_Trop(4) = 1E4 k_Trop(5) = 1E7

What I would ideally like to do is calculate the loss rate of NOx using a decision tree model that I have trained in python using sklearn (I have the model saved as .pkl file), where model parameters are then these meteorological inputs. I've read a little about implementing such a python model in Fortran and it seems possible. Do you think this would work to incorporate into the GeosChem code?

I am no longer needing to use O3, but still figuring out how I can incorporate NOx as a dummy species read in by HEMCO. Should I add to GeosCore/carbon_gases_mod.F90, or add a new module. Also for incorporating NO/NO2 emission data, should I directly amend HEMCO_Config.rc file to read in the files including NOx emission in ExtData/HEMCO, or is this something that needs to be incorporated into the source code.

Thank you! Chloe

[cschooling]

Hi @yantosca I am still working on implementing this. I have added NOx into the model as per my previous comment (with the NOx product rate and loss rate set as constant values). I am not reading emissions or background concentrations through HEMCO yet, but so far my NOx mixing ratios are coming out as nan so I assume I have missed something or have an error somewhere. Attached are the files that I have amended in the source code. (This is v 14.2.2)

Archive.zip

[yantosca]

Hi @cschooling, thanks for your patience. I was away (business, vacation, Thanksgiving) for much of November.

I think you are missing an initialization for DummyNOx in the carbon_convertKgtoMolecCm3 routine in carbon_funcs.F90:

    ! Initialize placeholder species to 1 molec/cm3
    C(ind_DummyCH4)   = 1.0_dp
    C(ind_DummyNMVOC) = 1.0_dp

    ! Initialize fixed species to 1 molec/cm3
    ! These will later be set to values read via HEMCO
    C(ind_FixedCl) = 1.0_dp
    C(ind_FixedOH) = 1.0_dp

You'll need to add:

C(ind_DummyNOx) = 1.0_dp

and I think that might fix the NaN values.

[cschooling]

Hi @yantosca no worries about the delay.

Initialising DummyNOx has helped things. However, I am now seeing very quick decays of NOx until the concentrations become negative and then become -infinity.

I have tried to prevent this by setting the rate of decay as proportional to the concentration of NOx. And forcing the rate to be zero if NOx concentration is below zero. However, this hasn't fixed it. Why would the model decay NOx when there is no NOx left to decay? How do I prevent the concentration from becoming negative?

This is what I have for setting k_Trop(4) (NOx decay rate) within the carbon_ComputeRateConstants routine:

   IF (ConcNOxMnd > 0) THEN
      k_Trop(4) = ConcNOxMnd * airvol_cm3 * 1E-3 
   ELSE
      k_Trop(4) = 0
   ENDIF

Within carbon_gases_mod.F90 the concentration of NOx, ConcNOxmnd, is read in as Global_NOx(I,J,L) (equivalent to GLOBAL_Cl) from full chemistry species concentration outputs. I think this is why my condition to set NOx to zero isn't working.

I actually want the rate of change to be dependent on the local NOx concentration. How do I pull this concentration?

Thanks!

[yantosca]

Hi @cschooling. You can use C(ind_NOx) to get the concentration of NOx in molec/cm3. This will be the NOx species concentration rather than an offline field that is read into GlobalNOx.

[cschooling]

Thanks @yantosca.

I am continuing to have issues with the concentration of NOx becoming negative.

The rate of the NOx decay is determined by k_Trop(4) - my understanding is the units of this is [1/s], so the number of reactions that occur per second within the grid.

I have the condition set so that we only decay the NOx if the concentration is greater than zero. I.e.

IF (C(ind_NOx) > 0) THEN
          k_Trop(4) = 0.1
ELSE
          k_Trop(4) = 0

I assume that a rate of 0.1, would then give 6 reactions per minute / 60 reactions per 10 minutes in a single grid cell. I have printed out the NOx concentrations for a given grid point, and we start off with a concentration of C(ind_NOx) ~ 0.7E8 [molecule/cm^3], but after 10 minutes this concentration becomes negative. I assume that the volume of the grid would mean that it should take much longer for such a concentration of NOx to decay to nothing. However, I print out the air volume, airvol_cm3, and get values on the order of E-311 which seem far too small! I think this is the reason we lose the NOx so quickly as there is actually nothing within this tiny volume.

1- Do you know why this volume would be coming out so small?

2- How would you best recommend to prevent the NOx concentrations from becoming negative?

Note this is all run with the emissions turned off to just track how the NOx decays to start with.

[yantosca]

Hi @cschooling. 1e-311 is a denormal number, which means that it's an underflow. Have you compiled with -DCMAKE_RELEASE_TYPE=Debug? That should let you know if you are encountering e.g. array out-of-bounds errors, or numerical issues like NaN/overflow/underflow.

[yantosca]

TL;DR: https://cs.stackexchange.com/questions/101632/understanding-denormalized-numbers-in-floating-point-representation

[cschooling]

Hello,

I am still working on developing this model. I have introduced the NOx tracer and am able to set the product/loss rates based on some meteorological variables.

However, I want to add NOx emissions into the simulation and I am struggling to do this.

In HEMCO_Config.rc I add pointers to the CESv2 NO emissions, which is turned on and already reading CO emission data.

0 CEDS_NO_AGR     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_agr            1750-2019/1-12/1/0 C xy   kg/m2/s NOx    25        1 5
0 CEDS_NO_ENE     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_ene            1750-2019/1-12/1/0 C xyL* kg/m2/s NOx    25/315    1 5
0 CEDS_NO_IND     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_ind            1750-2019/1-12/1/0 C xyL* kg/m2/s NOx    25/316    1 5
0 CEDS_NO_TRA     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_tra            1750-2019/1-12/1/0 C xy   kg/m2/s NOx    25        1 5
0 CEDS_NO_RCO     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_rco            1750-2019/1-12/1/0 C xy   kg/m2/s NOx    25        1 5
0 CEDS_NO_SLV     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_slv            1750-2019/1-12/1/0 C xy   kg/m2/s NOx    25        1 5
0 CEDS_NO_WST     $ROOT/CEDS/v2021-06/$YYYY/NO-em-anthro_CMIP_CEDS_$YYYY.nc            NO_wst            1750-2019/1-12/1/0 C xy   kg/m2/s NOx    25        1 5

I have also added to the HEMCO_Diagn.rc:

###############################################################################
#####  NOx sources                                                         #####
###############################################################################
EmisNOx_Total        NOx     -1     -1  -1   3   kg/m2/s  NO_emission_flux_from_all_sectors

However, when I run the model HEMCO never opens the NO-em-anthro_CMIPCEDS$YYYY.nc file, and there is not a variable output of EmisNOx_Total inside the HEMCO_diagnostics file.

Is there a step that I'm missing to introduce NOx emissions into the carbon run? My aim for this model is to be able to run inversions by perturbing the NOx emissions, so I would like to set it up in a way that will make this possible.

Thanks!

[yantosca]

Hi @cschooling, thanks for the update. Could you post the entire HEMCO_Config.rc file here? It could be that the NOx emissions are within a ((( ))) bracket that isn't being activated.

Also your entries look OK in HEMCO_Config.rc and HEMCO_Diagn.rc.

[cschooling]

Attached

HEMCO_Config.rc.zip

[yantosca]

Hi @cschooling. I think that HEMCO by default only adds emissions into the advected species. If your NOx species is not advected (but is only in the KPP mechanism), then it might not receive the emissions.

There is this code in GeosCore/hco_interface_gc_mod.F90:

    !=======================================================================
    ! Initialize HEMCO state object and populate it
    !=======================================================================

    !-----------------------------------------------------------------------
    ! Extract species to use in HEMCO. nHcoSpc denotes the number of
    ! species that shall be used in HEMCO. The species properties are
    ! defined in the Register_Species call below.
    ! Typically, nHcoSpc is just the number of species defined in both
    ! the HEMCO configuration file and GEOS-Chem. However, additional
    ! species can be defined, e.g. those not transported in GEOS-Chem
    ! (e.g. SESQ) or tagged species (e.g. specialty simulations).
    !-----------------------------------------------------------------------
    CALL SetHcoSpecies ( Input_Opt, State_Chm, HcoState,  nHcoSpc, 1, HMRC )

    ! Trap potential errors
    IF ( HMRC /= HCO_SUCCESS ) THEN
       RC     = HMRC
       ErrMsg = &
         'Error encountered in ""SetHcoSpecies" (first call, to get species)!'
       CALL GC_Error( ErrMsg, RC, ThisLoc, Instr )
       CALL Flush( HcoState%Config%Err%Lun )
       RETURN
    ENDIF

and then in SetHcoSpecies:

!EOC
!------------------------------------------------------------------------------
!                  GEOS-Chem Global Chemical Transport Model                  !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: SetHcoSpecies
!
! !DESCRIPTION: Subroutine SetHcoSpecies defines the HEMCO species. These
! are typically just the GEOS-Chem species. Some additional species may be
! manually added, e.g. SESQ (which is not a species) or individual CO2 species
! per emission source (for CO2 specialty sim).
!\\
!\\
! This routine has two phases: phase 1 simply returns the number of species
! to be used by HEMCO. This is useful as this number needs to be passed to
! the HEMCO initialization call.
! Phase 2 sets the HEMCO species information in the HEMCO state object. This
! needs to be done after initialization of the HEMCO state object.
! !INTERFACE:
!
  SUBROUTINE SetHcoSpecies( Input_Opt, State_Chm, HcoState, nSpec, Phase, RC )
!
! !USES:
!
    USE ErrCode_Mod
    USE HCO_LogFile_Mod, ONLY : HCO_SPEC2LOG
    USE Input_Opt_Mod,   ONLY : OptInput
    USE Species_Mod,     ONLY : Species
    USE HCO_Types_Mod,   ONLY : ConfigObj
    USE State_Chm_Mod,   ONLY : ChmState
!
! !INPUT PARAMETERS:
!
    INTEGER,          INTENT(IN   )   :: Phase      ! 1=Init, 2=Run
    TYPE(ChmState),   INTENT(IN   )   :: State_Chm  ! Chemistry State object
!
! !INPUT/OUTPUT PARAMETERS:
!
    TYPE(OptInput),   INTENT(INOUT)   :: Input_Opt  ! Input Options object
    TYPE(Hco_State),  POINTER         :: HcoState   ! HEMCO state
    INTEGER,          INTENT(INOUT)   :: nSpec      ! # of species for HEMCO
    INTEGER,          INTENT(INOUT)   :: RC         ! Success or failure?
!
! !REMARKS:
!  (1) We now get physical parameters for species from the species database,
!       which is part of the State_Chm object.
!  (2) In the future, it will be easier to specify non-advected species
!       like SESQ and the CO2 regional species from the species database.
!       The species database flags if a species is advected or not.
!
! !REVISION HISTORY:
!  06 Mar 2015 - C. Keller   - Initial Version
!  See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! LOCAL VARIABLES:
!
    ! Scalars
    INTEGER                :: nSpc, HMRC
    INTEGER                :: N,    L,    M
    REAL(dp)               :: K0,   CR,   pKa

    ! Strings
    CHARACTER(LEN= 31)     :: ThisName
    CHARACTER(LEN=255)     :: ThisLoc, Instr
    CHARACTER(LEN=512)     :: ErrMsg,  Msg

    ! Pointers
    TYPE(Species), POINTER :: SpcInfo

    !=================================================================
    ! SetHcoSpecies begins here
    !=================================================================

    ! Initialize
    RC       = GC_SUCCESS
    HMRC     = HCO_SUCCESS
    ErrMsg   = ''
    ThisLoc  = &
       ' -> at SetHcoSpecies (in module GeosCore/hco_interface_gc_mod.F90)'
    Instr    = 'THIS ERROR ORIGINATED IN HEMCO!  Please check the '       // &
               'HEMCO log file for additional error messages!'

    !-----------------------------------------------------------------
    ! For most simulations (e.g. full-chem simulation, most of the
    ! specialty sims), just use the GEOS-Chem species definitions.
    !-----------------------------------------------------------------
    IF ( Input_Opt%ITS_A_FULLCHEM_SIM                                   .or. &
         Input_Opt%ITS_AN_AEROSOL_SIM                                   .or. &
         Input_Opt%ITS_A_CARBON_SIM                                     .or. &
         Input_Opt%ITS_A_CO2_SIM                                        .or. &
         Input_Opt%ITS_A_CH4_SIM                                        .or. &
         Input_Opt%ITS_A_MERCURY_SIM                                    .or. &
         Input_Opt%ITS_A_POPS_SIM                                       .or. &
         Input_Opt%ITS_A_TAGCH4_SIM                                     .or. &
         Input_Opt%ITS_A_TAGCO_SIM                                      .or. &
         Input_Opt%ITS_A_TAGO3_SIM                                      .or. &
         Input_Opt%ITS_A_TRACER_SIM                                     .or. &
         Input_Opt%ITS_A_TRACEMETAL_SIM ) THEN

       ! Get number of model species
       nSpc = State_Chm%nAdvect

       !%%%%% FOR SOA SIMULATIONS %%%%%
       ! Check for SESQ: SESQ is not transported due to its short lifetime,
       ! but emissions are still calculated (in MEGAN). SESQ is only used
       ! in the SOA simulation, i.e. if LIMO is defined. Thus, add one more
       ! species here if LIMO is a model species and calculate SESQ emissions
       ! along with LIMO!
       IF ( id_LIMO > 0 ) THEN
          nSpc = nSpc + 1
       ENDIF

       !%%%%% FOR THE CARBON OR TAGGED CO SIMULATIONS %%%%%
       ! Add 5 extra species (ISOP, ACET, MTPA, LIMO, MTPO) for tagged CO
       IF ( Input_Opt%ITS_A_TAGCO_SIM ) THEN
          nSpc = nSpc + 5
       ENDIF

       ! Assign species variables
       IF ( PHASE == 2 ) THEN

          ! Verbose (only written if debug printout is requested)
          IF ( Input_Opt%Verbose ) THEN
             Msg = 'Registering HEMCO species:'
             CALL HCO_MSG( HcoState%Config%Err, Msg, SEP1='-' )
          ENDIF

          ! Sanity check: number of input species should agree with nSpc
          IF ( nSpec /= nSpc ) THEN
             WRITE(ErrMsg,*) 'Input species /= expected species: ', nSpec, nSpc
             CALL HCO_ERROR( HcoState%Config%Err, ErrMsg, RC, ThisLoc )
             RETURN
          ENDIF

          DO N = 1, State_Chm%nAdvect

             ! Get info for this species from the species database
             SpcInfo => State_Chm%SpcData(N)%Info

             ! Model ID and species name
             HcoState%Spc(N)%ModID      = SpcInfo%ModelID
             HcoState%Spc(N)%SpcName    = TRIM( SpcInfo%Name )

             ! Actual molecular weight of species [g/mol]
             HcoState%Spc(N)%MW_g       = SpcInfo%MW_g

             ! Set Henry's law coefficients
             HcoState%Spc(N)%HenryK0    = SpcInfo%Henry_K0   ! [M/atm]
             HcoState%Spc(N)%HenryCR    = SpcInfo%Henry_CR   ! [K    ]
             HcoState%Spc(N)%HenryPKA   = SpcInfo%Henry_pKa  ! [1    ]

             ! Logfile output (only written if debug printout is requested)
             IF ( Input_Opt%Verbose ) THEN
                CALL HCO_SPEC2LOG( HcoState, N )
             ENDIF

             ! Free pointer memory
             SpcInfo => NULL()
          ENDDO

          !------------------------------------------------------------------
          ! %%%%% FOR SOA SIMULATIONS %%%%%
          !
          ! Add the non-advected species SESQ in the last species slot
          !------------------------------------------------------------------
          IF ( id_LIMO > 0 ) THEN
             N                           = nSpec
             HcoState%Spc(N)%ModID       = N
             HcoState%Spc(N)%SpcName     = 'SESQ'
             HcoState%Spc(N)%MW_g        = 150.0_hp
             HcoState%Spc(N)%HenryK0     = 0.0_hp
             HcoState%Spc(N)%HenryCR     = 0.0_hp
             HcoState%Spc(N)%HenryPKa    = 0.0_hp

             ! Logfile output (only written if debug output is requested)
             IF ( Input_Opt%Verbose ) THEN 
                CALL HCO_SPEC2LOG(  HcoState, N )
             ENDIF
          ENDIF

          !------------------------------------------------------------------
          ! %%%%% FOR THE TAGGED CO SIMULATION %%%%%
          !
          ! Add the non-advected species ISOP, ACET, MTPA, LIMO, MTPO
          ! in the last 5 species slots (bmy, ckeller, 6/1/16)
          !------------------------------------------------------------------
          IF ( Input_Opt%ITS_A_TAGCO_SIM ) THEN

             ! Add 5 additional species
             DO L = 1, 5

                ! ISOP, ACET, MONX follow the regular tagged CO species
                M = State_Chm%nAdvect + L

                ! Get the species name
                SELECT CASE( L )
                   CASE( 1 )
                      ThisName = 'ISOP'
                   CASE( 2 )
                      ThisName = 'ACET'
                   CASE( 3 )
                      ThisName = 'MTPA'
                   CASE( 4 )
                      ThisName = 'LIMO'
                   CASE( 5 )
                      ThisName = 'MTPO'
                END SELECT

                ! Add physical properties to the HEMCO state
                HcoState%Spc(M)%ModID      = M
                HcoState%Spc(M)%SpcName    = TRIM( ThisName )
                HcoState%Spc(M)%MW_g       = 12.0_hp
                HcoState%Spc(M)%HenryK0    = 0.0_hp
                HcoState%Spc(M)%HenryCR    = 0.0_hp
                HcoState%Spc(M)%HenryPKa   = 0.0_hp

                ! Logfile output (only written if debug printout is requested)
                IF ( Input_Opt%Verbose ) THEN
                   CALL HCO_SPEC2LOG( HcoState, M )
                ENDIF
             ENDDO
          ENDIF

          ! Add line to log-file
          IF ( Input_Opt%Verbose ) THEN
             CALL HCO_MSG( HcoState%Config%Err, SEP1='-' )
          ENDIF
       ENDIF ! Phase = 2

    !-----------------------------------------------------------------
    ! DEFAULT (RETURN W/ ERROR)
    !-----------------------------------------------------------------
    ELSE
       ErrMsg = 'Invalid simulation type - cannot define model species'
       CALL HCO_ERROR( HcoState%Config%Err, ErrMsg, RC, ThisLoc )
       RETURN
    ENDIF

    ! For phase 1, pass species to output
    nSpec = nSpc

    ! Return w/ success
    RC = HCO_SUCCESS

  END SUBROUTINE SetHcoSpecies

By default, for the carbon simulation, the nHcoSpc is set to the number of advected species. You could add some code there to add to add NOx as an extra species for the carbon simulation.

[cschooling]

@yantosca Hello thanks for this! I have added to this to add NOx as an extra species, and am now able to run with the NOx emissions.

However, I would actually like for the NOx to also be an advected species. I see that in carbon_gases_mod.F90 the advected species IDs are defined here:

            N = State_Chm%Map_Advect(NA)

How can I adjust the mapping for advected species to add NOx?

[cschooling]

I figured out I could advect the species by just adding NOx into the list of transported species in geoschem_config.yml. So I have solved the above issue.

I am wanting to test the model by reading in product/loss rates of NOx that were produced from a full chemistry run - and then hopefully see analogous outputs to the full chemistry run. What I am a little confused by is what units I should be using for rates of product and loss. For the loss of NOx I have:

NOx             = Dummy             : k_Trop(4);

And for the product I have:

DummytoNOx      = NOx               : k_Trop(5);

Looking at similar equations in the carbon model it seems k_Trop(4) should have units of [1/s] and k_Trop(5) should have units of [molec/cm3/s].

I am reading two variables LNOx_chem and PNOx_chem from the full chemistry ProdLoss files, which have units of [molec/cm3/s]. Therefore I don't adjust units of k_Trop(5), but do adjust k_Trop(4) as follows:

k_Trop(4) = LNOx_chem * ( AIRMW * 1.0e-3_dp )      / ( AVO    * State_Met%AirDen(I,J,L) * 1.0e-6_dp )   !! [1/s]

k_Trop(5) = PNOx_chem    !! [molec/cm3/s]

However, in doing this I get a large net gain of NOx as k_Trop(5) is significantly larger than k_Trop(4). Have I understood this correctly or should they both have units of [1/s]?

[cschooling]

@yantosca Hi, here with another query.

I am wanting to read the NOx product and loss rates from full chemistry outputs, I want this to be read and updated every 10 minutes. I have these pointers in my HEMCO_Config.rc file:

((PROD_NOx_CHEM
* Prod_NOx    $ROOT/GCClassic_Output/14.0.0/$YYYY/GEOSChem.ProdLoss.$YYYY$MM$DD_$HH$MNz.nc4         Prod_NOx     2019-2020/1-12/1/0 C xyz 1        * - 1 1
))PROD_NOx_CHEM
((LOSS_NOx_CHEM
* Loss_NOx    $ROOT/GCClassic_Output/14.0.0/$YYYY/GEOSChem.ProdLoss.$YYYY$MM$DD_$HH$MNz.nc4         Loss_NOx     2019-2020/1-12/1/0 C xyz 1        * - 1 1
))LOSS_NOx_CHEM

However, when I run the model HEMCO only opens the ProdLoss file for the starting time point. Is there a way to adjust this so it opens and reads from a new file in each timestep. Or would the only way to do it be to calculate relative scaling factors to apply?

In regards to units of k_Trop(4)/k_Trop(5) it seems that they need to do be in the same units or the balance between the loss/product becomes disproportionate. Keeping molecules/cm3/s the rate is far too high and the NOx all decays away at a much faster rate than it is emitted. But using units 1/s the NOx appears to decay too slowly in comparison to what I see when running full chemistry. Is there something I'm missing here?

Thanks so much!

[zion0210]

Hello, @yantosca and @cschooling!

Did you fix the following problem: "GEOS-Chem ERROR: Cannot get pointer to HEMCO field CH4_OIL" ?

I run difference cases of GEOS-Chem 14.3.0. Also, I set the horizontal grid with 4.0x5.0 deg resolution and the Z grid with 47 dimensions. CO2 - everything good, CH4 - everything good, Aerosols - everethyng good, Carbon - the following error: GEOS-Chem ERROR: Cannot get pointer to HEMCO field CH4_OIL -> at Emiss_Carbon_Gases (in module GeosCore/carbon_gases_mod.F90)

Do I need to use an older version of the GEOS-Chem for carbon case? Or can you suggest another solution to my problem?

Best regards, zion0210

[yantosca]

@zion0210, would you be able to open a new issue for your question? Thanks.

[cschooling]

Hi @yantosca.

I have now successfully set up the carbon model to include the NOx species, where I read the product/loss rates from full chemistry run outputs. When I turn off transport the difference between the full chemistry and tagged NOx run are in good agreement, with a small deviation of maximum 1.5% difference between the NOx columns appearing after one day of running. Would such a difference be expected, even when our product loss rates are taken from full chemistry and the restart is the same?

The thing I am more concerned about is that when the transport is turned on the difference in NOx columns between the two models increases quite significantly (up to around 25% deviation after one day). Is this because of the fact that NOx is treated as two individual species NO+NO2 in full chemistry but just one species in the tagged model so it is transported differently? Do you have any advice on the source of this error and if it is possible to reduce it?

[yantosca]

Hi @cschooling, thanks for your patience. It could be that the transport of the separate NO and NO2 species is causing the issue. I wouldn't know much more on what to do about it. You might try contacting the Transport WG and/or the Chemistry WG for advice.

[yantosca]

Closing out this issue