[QUESTION] How to archive emission variable in GEOS-Chem?

[debatoshpartha]

Overview:

I'm trying to calculate and parameterize a specific species' emission from a particular fraction of global land. I have written the source code for it which uses dry and wet nitrate deposition. I have two questions,

1. How to properly call and use the stored data of the dry and wet deposition of nitrate from the archived 'DryDep' and 'WetLossLS' from State_Diag_Mod.F90 in any newly created module?

   2 How to properly archive the final emission variable of the equation in geoschem so that I can get valid values in netcdf output files? I'm following the below guidelines of archiving variable in state_diag_mod.F90 and in HISTORY.rc,

   

However, I'm not sure how to correctly populate the state_diag archived array with the values I want to output. I'm passing the archived variable in my own emission module like this,

    ! dep_nit is the sum of dry and wet nitrogen deposition and I'm planning to call them from state_diag   
    emis_x(I,J) = emis_x(I,J) +  (dep_nit(I,J) * constants * global land fraction(I,J) ! x is the species for example
        !X emissions [kg/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = State_Diag%EmisX(I,J) + emis_x(I,J) ! EmisX is archived in state_diag following the guideline mentioned above and in HISTORY.rc collection as well. 
       ENDIF

Any general or specific suggestion on these issues would be really great. Thanks!

[yantosca]

Thanks for writing @debatoshpartha. The way you have described it above will work just fine. Then you need to add EmisX to one of the collections in the HISTORY.rc file.

We are in the process of updating the GEOS-Chem and HEMCO diagnostics for the upcoming 14.0.0 release. Some of the docs on the wiki may be a bit out of date. Please see:

• https://geos-chem.readthedocs.io
• https://hemco.readthedocs.io

These are about 95% complete but we are adding/editing a few pages.

[debatoshpartha]

Thanks Mr. @yantosca for your reply. I already had a newly added collection named 'XEmis' and added EmisX to that collection like shown below, 'X' being the species of interest from a specific source here,

COLLECTIONS: 'Restart',
             #'Metrics',
             #'SpeciesConc',
             'XEmis',  !** newly created collection **
             #'AdvFluxVert',
             #'AerosolMass',
             #'Aerosols',
             #'Budget',
             #'CloudConvFlux',
             ##'ConcAboveSfc',
             #'ConcAfterChem',
             #'DryDep',
             #'JValues',
             #'KppDiags',
             #'LevelEdgeDiags',
             #'ProdLoss',
             ##'RRTMG',
             ##'RxnRates',
             #'StateChm',
             #'StateMet',
            ##'StratBM',
             #'WetLossConv',
             #'WetLossLS',
             ##'BoundaryConditions',

#==============================================================================
  XEmis.template:       '%y4%m2%d2_%h2%n2z.nc4',
  XEmis.format:         'CFIO',
  XEmis.frequency:      00000001 000000
  XEmis.duration:       00000001 000000
  XEmis.mode:           'time-averaged'
  XEmis.fields:         'EmisX                   ', 'GIGCchem',                        
::

But I got zero values of the variable 'EmisX' under the newly created daily output file named 'GEOSChem.XEmis.20150101_0000z.nc4'. Have I missed something here?

Question 2: Initially, I'm using this 'get_nit_dep_mod.F90' to get the dry and wet deposited nitrate on soil but I know the dry and wet deposition/loss of nitrate is already archived in 'DryDep' and 'WetLossLS' collections for all species. Could you please tell me how can I use this archived data into my x_mod.F90 so that I don't need to use another module to calculate the nitrate deposition?

[yantosca]

Thanks for the feedback @debatoshpartha. One thing: I think I may have erred in one of my replies to you earlier. You have:

    ! Snow NOx emissions [kg/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = State_Diag%EmisX(I,J) + emis_x(I,J)
       ENDIF

but it is sufficient to do

    ! Snow NOx emissions [kg/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = emis_x(I,J)
       ENDIF

because the HISTORY component will handle summing the data (if it is a time-averaged collection).

Also: in state_diag_mod.F90, have you registered the variable X_Emis? If you haven't then the HISTORY diagnostic won't know what to do with it. You can take a look at the other calls to Init_and_Register in Init_State_Diag to see how to set that up.

Also in your get_ndep_mod.F90, you don't need to make e.g. State_Diag%WetDepNIT and State_Diag%DryDepNIT. All of the fields in State_Chm can be requested in HISTORY.rc by prefixing the field name with Chem_. So if you added Chem_WetDepNIT and Chem_DryDepNIT to your HISTORY.rc file (you can even put these in a new collection if you wish) then you will get output.

• CAVEAT: If you added State_Chm%DryDepNIT and State_Chm%WetDepNIT yourself, you will need to make sure that you have registered these variables with calls to Init_and_Register in state_chm_mod.F90. Otherwise the HISTORY diagnostics won't know that they are archivable.

Also -- in general it is not advisable to use fields of State_Diag for anything other than diagnostic output. If you need a field to be carried from one module to another, you should place it into State_Chm, because pretty much all routines will take State_Chm as an argument. Also, if you are using e.g. State_Diag%WetDepNIT in another module and that diagnostic has been turned off, your GEOS-Chem run will stop with an error because the State_Diag%WetDepNIT array isn't allocated.

Hope this helps!

[debatoshpartha]

Thank you Mr. @yantosca for correcting the IF loop. I have revised that.

I'm not sure if you're trying to mention the variables EmisXor emis_x as X_Emis cause I don't have any variable named X_Emis in my modules, this is just the name of the HISTORY.rc collection. But to answer your question, yes, right now I'm using the same variable name as EmisX both in my main equation (see below) and in state_diag_mod.F90 and I have registered this variable in Init_and_Register in Init_State_Diag and in other places of state_diag_mod.F90 as well.

! dep_nit is the sum of dry and wet nitrogen deposition
! I'm calculating the dry deposition of nitrate below (followed aero_drydep.F90), 

!=================================================================    
! DEP_NIT_FLUX begins here!
!=================================================================

    ! Assume success
    RC        = GC_SUCCESS

    ! Number of dry-deposited species
    nDryDep   = State_Chm%nDryDep

    ! DTCHEM is the chemistry timestep in seconds
    DTCHEM    = GET_TS_CHEM()

    ! Initialize pointers
    Spc                  => State_Chm%Species  ! Chemistry species [kg]

     ! First-time setup
    IF ( FIRST ) THEN
! Define species ID flags
       id_NIT = Ind_('NIT')

    !$OMP PARALLEL DO                                      &
    !$OMP PRIVATE( L, J, I, AREA_CM2, BIN ) &
    !$OMP PRIVATE( ID, X, Y0, Y )                      &
    !$OMP DEFAULT( SHARED )                                &
    !$OMP SCHEDULE( DYNAMIC )
    DO I = 1, State_Grid%NX
    DO J = 1, State_Grid%NY
    DO L = 1, State_Grid%MaxChemLev

    ! Initialize for safety's sake
       AREA_CM2 = 0d0

       ! Dry deposit NIT (DBP, 7/28/22)
       Y0 = Spc(I,J,L,id_NIT)

       !==============================================================
       ! Drydep flux of NIT [molec/cm2/s]
       !==============================================================

          ! Surface area [cm2]
          AREA_CM2 = State_Grid%Area_M2(I,J) * 1e+4_fp

          ! Convert from [kg/timestep] to [molec/cm2/s]
          DRYNIT_FLUX = Y0 
          DRYNIT_FLUX = DRYNIT_FLUX / (State_Chm%SpcData(id_NIT)%Info%MW_g * &
                 1.e-3_fp) / AREA_CM2 / DTCHEM * AVO

      ! ENDIF

    ENDDO
    ENDDO
    ENDDO

    !$OMP END PARALLEL DO
    ! Free pointers
    Spc      => NULL()

    ! Check that species units are still in [kg] (ewl, 8/13/15)
    IF ( TRIM( State_Chm%Spc_Units ) /= 'kg' ) THEN
       MSG = 'Incorrect species units at end of routine: ' &
             // TRIM(State_Chm%Spc_Units)
       LOC = 'Routine EmisX_FLUX in emisx_mod.F90'
       CALL GC_Error( MSG, RC, LOC )
    ENDIF

    !=================================================================
    ! EmisX_FLUX begins here!
    !=================================================================

 ! Initialize
   EmisX(I,J) = 0e+0_fp

    ! Emission timestep [s]
    !DTSRCE = GET_TS_EMIS()

    !$OMP PARALLEL DO       &
    !$OMP DEFAULT( SHARED ) &
    !$OMP PRIVATE( I,              J,    L,            NN                 ) &
    !$OMP PRIVATE( FRAC_SNOW_OR_ICE ) &
    !$OMP PRIVATE( IS_SNOW_OR_ICE                                       )

    DO I  = 1, State_Grid%NX
    DO J  = 1, State_Grid%NY

    ! Don't let fraction be greater than 1
    FRAC_SNOW_OR_ICE = MIN( State_Met%FRSNO(I,J)     + &
                            State_Met%FRSEAICE(I,J)  + &
                            State_Met%FRLANDIC(I,J), 1e+0_fp)
    IS_SNOW_OR_ICE   = ( FRAC_SNOW_OR_ICE > 0e+0_fp )

    ! Check if there is snow on the ground, or if this is sea ice
    IF ( IS_SNOW_OR_ICE ) THEN

        ! Mass of emitted EmisX, molec/cm2/s
          EmisX(I,J) = EmisX(I,J) + ( DRYNIT_FLUX * constants * FRAC_SNOW_OR_ICE)

          ELSE

          ! No flux from non-land surfaces 
          EmisX(I,J) = 0e+0_fp

       ENDIF

    ! Archive EmisX flux [molec/cm/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = EmisX(I,J)
       ENDIF

    ENDDO
    ENDDO
    !$OMP END PARALLEL DO

I'm also attaching the state_diag_mod.F90 with my addition of EmisX. You can search with my initials 'DBP' to see those additions/modifications.

state_diag_mod.txt

Also I've tried running 1-day simulation with Chem_DryDepNIT and Chem_WetDepNIT under my own HISTORY.rc collection XEmis but all the variables i.e. EmisX,Chem_DryDepNIT and Chem_WetDepNIT in the output came out as zero. I've also deleted DryDepNIT and WetDepNIT from state_diag_mod.F90 and only kept them in state_chm_mod.F90.

I'm still not sure what's going wrong here. Any suggestions/insights here?

Question:

How can I specify/tell the GEOSChem model that I am trying to calculate a specific chemical species emission in my defined variable EmisX? I know there is a way to add the species of interest in HEMCO_Diagn.rc but how can I do it in GEOSChem?

[debatoshpartha]

Hi Mr. @yantosca, haven't heard from you or any other GCST member lately. Could anyone please share their opinion on how to solve this kind of issue?

[yantosca]

Thanks for your patience @debatoshpartha. We have had several things happening on our end this month.

I see that State_Diag%Emis_X is the array name. The name of the collection (which you said is XEmis) is what you need to specify in Init_and_Register via the diagId field (this can be different than the name of the array in State_Diag if you wish). That should all be OK.

Also a couple of things I noticed about your code above:

1. It looks as is if the DO loop where DRYNIT_FLUX is computed is in the IF ( FIRST ) THEN block. Code in this block will only execute the first time the simulation is called. I am not sure if this is what you intended but this may explain why the diagnostic is not updating. The only thing that should be in this IF block is the setting of id_NIT and then resetting FIRST = .FALSE. so that the IF block won't execute again.

   IF ( FIRST ) THEN 
      id_NIT = Ind_('NIT')
      FIRST = .FALSE.
   ENDIF

2. For better execution efficiency, the order of the DO loops should reversed, as this will step through arrays in the same order that they are laid out in memory. So use this ordering instead:

   DO L = 1, State_Grid%MaxChemLev
   DO J = 1, State_Grid%NY
   DO I = 1, State_Grid%NX

3. Try commenting out the !$OMP statements by adding an extra ! in front of them (i.e. !!$OMP) while you are developing. This will prevent any parallelization errors from confusing your output. When you have everything working again you can remove the extra comment (which will restore the parallelization commands). For more information about Parallelization in GEOS-Chem, see our wiki page: http://wiki.geos-chem.org/Parallelizing_GEOS-Chem.

4. You can also avoid the ELSE statement (which may lead to more optimal execution) by rewriting the code where you add into State_Diag%EmisX as:

     ! No flux from non-land surfaces 
     EmisX(I,J) = 0e+0_fp

     ! Check if there is snow on the ground, or if this is sea ice
     IF ( IS_SNOW_OR_ICE ) THEN

        ! Mass of emitted EmisX, molec/cm2/s
        EmisX(I,J) = EmisX(I,J) + ( DRYNIT_FLUX * constants * FRAC_SNOW_OR_ICE)

     ENDIF

     ! Archive EmisX flux [molec/cm/s]
     IF ( State_Diag%Archive_EmisX ) THEN
         State_Diag%EmisX(I,J) = EmisX(I,J)
     ENDIF

So these are a few things to try. Let us know how it goes!

[yantosca]

Oops I did not mean to close the issue, I just hit the wrong button. Now reopened.

Also tagging @Jourdan-He

[debatoshpartha]

Thank you Mr. @yantosca for your reply. I have done all the things you've suggested yet the EmisX variable value is coming out as zero. I have renamed the HISTROY.rc collection name as EmisX from XEmis as it seemed more straightforward and is already defined in state_diag_mod.f90. The HISTORY.rc file is copied below,

COLLECTIONS: 'Restart',
             #'Metrics',
             #'SpeciesConc',
             'EmisX,
#=======================================================================
  EmisX.template:        '%y4%m2%d2_%h2%n2z.nc4',
  EmisX.format:           'CFIO',
  EmisX.frequency:     00000000 010000,
  EmisX.duration:        00000000 010000,
  EmisX.mode:            'time-averaged'
  EmisX.fields:             'EmisX.       ', 'GIGCchem',
::

I have also revised the calculation of EmisX variable as below,

   !=================================================================
    ! EmisX begins here!
    !=================================================================

     ! Assume success
    RC        = GC_SUCCESS

    ! DTCHEM is the chemistry timestep in seconds
   DTCHEM    = GET_TS_CHEM()

    ! Initialize pointers
    Spc                  => State_Chm%Species  ! Chemistry species [kg]

     ! First-time setup
      IF ( FIRST ) THEN 
      id_NIT = Ind_('NIT')
      FIRST = .FALSE.
   ENDIF

 !!$OMP PARALLEL DO       &
    !!$OMP DEFAULT( SHARED ) &
    !!$OMP PRIVATE( I,              J,    L,            NN                 ) &
    !!$OMP PRIVATE( FRAC_SNOW_OR_ICE ) &
    !!$OMP PRIVATE( IS_SNOW_OR_ICE                                       )                                    
    !!$OMP PRIVATE( AREA_CM2 ) &
    !!$OMP PRIVATE( ID, X )                      &                              &
    !!$OMP SCHEDULE( DYNAMIC )

    DO L = 1, State_Grid%MaxChemLev
    DO J = 1, State_Grid%NY
    DO I = 1, State_Grid%NX

    ! Initialize for safety's sake
       AREA_CM2 = 0d0
       !==============================================================
       ! Drydep flux of NIT [molec/cm2/s]
       !==============================================================

          ! Surface area [cm2]
          AREA_CM2 = State_Grid%Area_M2(I,J) * 1e+4_fp

          ! Convert from [kg/timestep] to [molec/cm2/s]

          DRYNIT_FLUX = Spc(I,J,1,id_NIT)/ (State_Chm%SpcData(id_NIT)%Info%MW_g * &
                 1.e-3_fp) / AREA_CM2 / DTCHEM * AVO

    ! Don't let fraction be greater than 1
    FRAC_SNOW_OR_ICE = MIN( State_Met%FRSNO(I,J)     + &
                            State_Met%FRSEAICE(I,J)  + &
                            State_Met%FRLANDIC(I,J), 1e+0_fp)
    IS_SNOW_OR_ICE   = ( FRAC_SNOW_OR_ICE > 0e+0_fp )

     ! No flux from non-land surfaces (snow & ice)
         EmisX(I,J) = 0e+0_fp

    ! Check if there is snow on the ground, or if this is sea ice
    IF ( IS_SNOW_OR_ICE ) THEN

        ! Mass of EmisX, kg !
          EmisX(I,J) = EmisX(I,J) +   DRYNIT_FLUX(I,J) * constants * FRAC_SNOW_OR_ICE

       ENDIF

    !EmisX [kg/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = EmisX(I,J)
       ENDIF

   ! Print out ** Just FYI, this portion is not printing out while the simulation is running **
          WRITE( 6, 100 ) REPEAT( '=', 79 )
          WRITE( 6, 110 ) SUM( EmisX )
          WRITE( 6, 100 ) REPEAT( '=', 79 )

100    FORMAT( a                                            )
110    FORMAT( 'EmisX   E M I S S I O N S : ', f7.3, ' [molec/cm2/s]' )

    ENDDO
    ENDDO
    ENDDO
    !!$OMP END PARALLEL DO

I have also tried to remove DRYNIT_FLUX(I,J) from the main equation to check if I can get any values in the EmisX variable of the NetCDF output but I got no values. I am not sure why my EmisX array is not populating with values. Could you please take a look? I am attaching the state_diag_mod.F90 below,

state_diag_mod.txt

Also if you could mention how can I direct/specify GEOS-Chem that I am trying to calculate a specific species emission and add it to the already existing array of the species emission? For example, there is dust emission already in GEOSChem and suppose I am trying to calculate dust emission from a specific source that's not yet considered in GEOSChem, so after calculating the amount of dust from that specific source, how can I add that emission to the already existing dust emission?

[yantosca]

@debatoshpartha: Most emission diagnostics are done via HEMCO and not through the GEOS-Chem History diagnostics. You would have to request a new diagnostic in HEMCO_Diagn.rc for the new emissions source as a separate container. Then you can combine the data in post-processing to get totals, etc. We have some new documentation on the HEMCO diagnostics at https://hemco.readthedocs.io.

[yantosca]

@debatoshpartha Your state_diag_mod.F90 looks good.

Have you checked to see that the subroutine in which state_diag%EmisX is being updated is being called? If for some reason the subroutine isn't being called then there wouldn't be any emissions being added to the diagnostic.

[debatoshpartha]

Hi Mr. @yantosca, I have couple of updates, I have added a new diagnostics inHEMCO_DIagn.rc called EmisXsourceunder NO emission as below,

###############################################################################
#####  NO emissions                                                       #####
#####                                                                     #####
##### - Separate fertilizer NOx emissions are only available when the     #####
#####   SoilNOx extension is enabled                                      #####
##############
#################################################################
EmisNO_Total       NO     -1     -1  -1   3   kg/m2/s  NO_emission_flux_from_all_sectors
EmisNO_Aircraft    NO     0  20  -1   3   kg/m2/s  NO_emission_flux_from_aircraft
EmisNO_Anthro      NO     0  1   -1   3   kg/m2/s  NO_emission_flux_from_anthropogenic
EmisNO_BioBurn     NO     111    -1  -1   2   kg/m2/s  NO_emission_flux_from_biomass_burning
EmisNO_Lightning   NO     103    -1  -1   3   kg/m2/s  NO_emission_flux_from_lightning
EmisNO_Ship        NO     102    -1  -1   2   kg/m2/s  NO_emission_flux_from_ships
EmisNO_Soil        NO     0      3   -1   2   kg/m2/s  NO_emission_flux_from_soil_including_fertilizer
#EmisNO_Fert        -1     104    -1  -1   2   kg/m2/s  NO_emission_flux_from_fertilizer_only
EmisXsource        NO     109    -1   -1   2   mole/cm2/s  X_emission_flux_from_source

but I am not sure how to pass our GESO-Chem/GeosCore code variable EmisX to HEMCO. Previously in one of my github issues #315, you suggested to use the bridge module GeosCore/hcoi_gc_main_mod.F90 and pass data to HEMCO via State_Chm. In the v13 of GEOSChem, I don't see the module hcoi_gc_main_mod.F90 anymore. Could you please direct me on the way on how to pass data of GEOSCHEM variable EmisX to EmisXsource in HEMCO?

Also my GEOS-Chem/GeosCore module that I'm using to calculate EmisX is being compiled. Still, when I try to print out something in the module while running the model simulation, it doesn't print out anything and finishes the simulation. I am attaching my module here,

emisx_mod.txt

[yantosca]

Thanks for your patience @debatoshpartha. We have had a lot of things going on here lately.

In GC 13 I believe hcoi_gc_main_mod.F90 has been renamed to hco_interface_gc_mod.F90. Sorry for that.

Also in your emiss_mod.txt I noticed that you are using AD32. That is not used unless you build GEOS-Chem with -DBPCH_DIAG=y. That array is from the older emissions diagnostics which are being phased out. In GC 13 I think we have since removed that array. You can just do:

       ! source X emissions [kg/s]
       IF ( State_Diag%Archive_EmisX ) THEN
          State_Diag%EmisX(I,J) = EmisX(I,J)
       ENDIF

Also, the History diagnostics are separate from the HEMCO diagnostics. So if you add something in State_Diag it won't be printed out via HEMCO (and can't be scheduled via HEMCO_Diagn.rc). We do have some documentation on HEMCO diagnostics online: see this link.

Most emission diagnostics are handled via HEMCO. But since you are doing something in the GEOS-Chem code you might want to implement that as a History diagnostic (i.e. via State_Diag as you have done, and scheduled in HISTORY.rc).

Also I think I see another problem. In your SourcePack_X_Flux routine you are passing EmisX as an argument. But at the top of your module you have:

MODULE EmisX_MOD
!
! !USES:
!
  USE HCO_ERROR_MOD    ! For real precisions (hp)
  USE PRECISION_MOD    ! For GEOS-Chem Precision (fp, f4, f8)

  IMPLICIT NONE
  !PRIVATE

! !PUBLIC DATA MEMBERS:
!
  REAL(fp),  ALLOCATABLE :: EmisX(:,:)  <---(1) This is declaring EmisX as a global module variable; 
                                        <--- it is visible to everything below the CONTAINS statement.   
! Pointers
    REAL(fp), POINTER :: Spc(:,:,:,:)

  !============================================================================
  !PUBLIC 
  !============================================================================
  PUBLIC :: SOURCEPACK_X_FLUX
  PUBLIC :: INIT_EmisX
  PUBLIC :: CLEANUP_EmisX

..etc ...

CONTAINS

 SUBROUTINE SOURCEPACK_X_FLUX( EmisX, State_Grid, State_Met, State_Diag, State_Chm, RC )

!
! !USES:
!
... etc ...

!
! !INPUT PARAMETERS:
!
    TYPE(GrdState), INTENT(IN)  :: State_Grid  ! Grid State object
    TYPE(MetState), INTENT(IN)  :: State_Met   ! Meteorology State object
!
! !OUTPUT PARAMETERS:
!
!X flux [kg/s]

    REAL(fp),       INTENT(INOUT) :: EmisX(State_Grid%NX,State_Grid%NY)    
    ^^^^^ (2) This is declaring EmisX as an argument, which will override the global module variable

    TYPE(DgnState), INTENT(INOUT) :: State_Diag  ! Diagnostics State object
    TYPE(ChmState), INTENT(INOUT) :: State_Chm   ! Chemistry State object
    INTEGER,        INTENT(OUT)   :: RC          ! Success or failure?   

so in your subroutine this statement:

    ! Check if there is source on the ground, or if this is sea ice
    IF ( IS_SNOW_OR_ICE ) THEN

        ! Mass of emitted EmisX, kg !
          EmisX(I,J) = EmisX(I,J) +   DRYNIT_FLUX(I,J) * constants * FRAC_SNOW_OR_ICE

       ENDIF

will update the array that is being passed into SourcePack_X_Flux but not the module variable.

Long story short: you can't have both constructs (1) and (2) in your code. If you don't need to use Emis_X as a global variable I would get rid of (1),

If you are new to Fortran I would recommend taking a look at at some of these links on our Fortran language resources wiki page. Also see this link: http://www.cs.rpi.edu/~szymansk/OOF90/bugs.html.

[yantosca]

@debatoshpartha: if you get rid of (1), then you can also get rid of routines Init_EmisX and Cleanup_EmisX.

[yantosca]

Also, @debatoshpartha, in which routine are you calling SourcePack_X_Flux from? You should have a CALL SourcePack_X_Flux(...) statement somewhere in the code.

[debatoshpartha]

Hi Mr. @yantosca, your last couple of replies helped me out and now I am getting values out of EmisX variable, finally! Thank you very much for taking the time to reply. I have couple of follow-up questions,

1. Longer simulation run-time: Until this point, I was running 2x2.5 resolution, 1-hour simulations to test my code but right now, I am trying to run a 2-year simulation (2014-2015, 2014 is for spin-up simulation) but it seems with the updated version of GEOS-Chem v13, the simulation runtime is much more longer (it's taking 1 day to complete 6-days of simulation where I have changed the transport and chemical timesteps as 900 and 1800 seconds according to this doc http://wiki.seas.harvard.edu/geos-chem/index.php/Speeding_up_GEOS-Chem. Could you please suggest how to increase the simulation speed so that I can run 1-month of simulation in 1 day?

2. Adding newly calculated flux to existing flux: In my newly developed code, I'm calculating NOx emission from a specific source, how can we add this new emission flux to the already existing total NOx flux in geoschem model so that it can be added up with the total chemistry scheme? I have tried adding EmisX to the total NOx flux the way copied below but the simulation doesn't go through in this way,

   ! Add EmisX to total NOx flux and archive it
   
     IF (State_Diag%Archive_NOFLUX) THEN
   
   !calculate the already existing NOx flux
   NOFLUX = Spc(I,J,1,id_NO)/ (State_Chm%SpcData(id_NO)%Info%MW_g * 1.e-3_fp) / AREA_CM2 / DTCHEM * AVO
   
    NOFLUX(I,J) = NOFLUX(I,J) + EmisX(I,J) 
   
    State_Diag%NOFLUX(I,J) = State_Diag%NOFLUX(I,J) + NOFLUX(I,J)
   
     ENDIF

   Please note that I have defined and registered the NOFLUX in both state_diag and state_chm and added this into my EmisX HISTORY.rc collection.

[yantosca]

@debatoshpartha: Due to travel, I will be unable to respond to your issue for the next couple of weeks. Tagging @Jourdan-He @lizziel @msulprizio. Thanks for your understanding.

[debatoshpartha]

Hi @Jourdan-He @lizziel @msulprizio, could you please take some time and reply to this comment https://github.com/geoschem/geos-chem/issues/1322#issuecomment-1231960509?

[lizziel]

Hi @debatoshpartha, one of us will take a look at this issue next week.

[debatoshpartha]

Hi @lizziel, thanks for replying. Do you or anyone have any thoughts on https://github.com/geoschem/geos-chem/issues/1322#issuecomment-1231960509?

[debatoshpartha]

Hi Mr. @yantosca, I hope you're doing good and back from travel. During this time, I have tried reaching out to GCST but unfortunately haven't received any replies from them. I have tried different ways to solve my issue of adding new emission to total NOx emissions but couldn't do so. Right now, the main issue is to add the new sector's NOx emission to the total NOx emission (which is archived in HEMCO as EmisNO_Total) so that the GEOS-Chem model can address the new emission flux and update its chemistry mechanism accordingly. If you have some insights, please provide some.

I have tried different ways of calling IDTNO which is the instance for total NOx emission (saw it in hcox_soilnox_mod.F90and hcox_lightnox_mod.F90) and added my EmisX toInst%IDTNO like this (also added subroutine 'InstGet' in the same module.

! Get Instance
    CALL InstGet ( ExtState%EmisX, Inst, RC )
    IF ( RC /= HCO_SUCCESS ) THEN
       WRITE(MSG,*) 'Cannot find X NOx instance Nr. ', ExtState%EmisX
       CALL HCO_ERROR(HcoState%Config%Err,MSG,RC)
       RETURN
    ENDIF

    ! Add flux to emission array
    CALL HCO_EmisAdd( HcoState, EmisX, Inst%IDTNO, &
                      RC,       ExtNr=Inst%ExtNr )
    IF ( RC /= HCO_SUCCESS ) THEN
       CALL HCO_ERROR( HcoState%Config%Err, 'HCO_EmisAdd error', RC )
       RETURN
    ENDIF

but whenever I do it (i.e. call and write anything to HEMCO), the simulation just stops there after this line,

********************************************
* B e g i n   T i m e   S t e p p i n g !! *
********************************************

---> DATE: 2014/01/03  UTC: 00:00  X-HRS:      0.000000
 HEMCO already called for this timestep. Returning.
===============================================================================
TPCORE_FVDAS (based on GMI) Tracer Transport Module successfully initialized
===============================================================================
HEMCO (VOLCANO): Opening /wsu/home/groups/yhgroup/ExtData/HEMCO/VOLCANO/v2019-08/2014/01/so2_volcanic_emissions_Carns.20140103$
--- Initialize surface boundary conditions from input file ---

I am not sure if I can use GEOS-Chem to calculate emissions anymore. So I have shifted my code to HEMCO/Extensions (attached), added a new module of nitrate deposition (attached) cause my emission module uses nitrate deposition and I passed DryDepNIT and WetDepNIT from GEOS-Chem to HEMCO via hco_interface_gc_mod.F90 and made necessary changes/registration of arrays to hcox_state_mod.F90 and state_chm_mod.F90.

get_nitdep_mod.txt

hcox_emisx_mod.txt

For HEMCO diagnostics, I have added EmisX in HEMCO_Config.rc (attached) at L158. I have also added the diagnostic contained EmisNO_X in HEMCO_Diagn.rc (attached) at L344.

HEMCO_Diagn.txt

HEMCO_Config.txt

But I am not getting values in HEMCO_diagnostics.YYYYMMDDMMSS.nc, I have checked the hcox_emisx_mod.F90 to see if there's any issue of variable overriding just like we saw in our previous GEOS-Chem emisx_mod.F90 but I haven't encountered any.

Now the question is,

1. Do I have to shift to HEMCO/Extension code to calculate and add this emission to the total NOx emission or GEOS-Chem code will work just fine? If GEOS-Chem code works, how can I add this new emission to the total NOx emission?
2. If calculating and adding new emission in HEMCO/Extension makes things straightforward, how to archive EmisNO_X diagnostic and get values out of it and add this emission to the total NOx emission flux?

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. If there are no updates within 7 days it will be closed. You can add the "never stale" tag to prevent the Stale bot from closing this issue.

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. If there are no updates within 7 days it will be closed. You can add the "never stale" tag to prevent the Stale bot from closing this issue.

[msulprizio]

Hi @debatoshpartha. We recommend doing this via a HEMCO extension. Please see my response in https://github.com/geoschem/geos-chem/issues/1532.