search

geoschem / geos-chem

GEOS-Chem "Science Codebase" repository. Contains GEOS-Chem science routines, run directory generation scripts, and interface code. This repository is used as a submodule within the GCClassic and GCHP wrappers, as well as in other modeling contexts (external ESMs).
http://geos-chem.org
Other
164 stars 156 forks source link

[ISSUE] Wrong unit in the RxnRates Collection in GEOS-Chem version 13.4.1? #1459

Closed xin-tie closed 1 year ago

xin-tie commented 1 year ago

Ask a question about GEOS-Chem:

Hi there,

We have been running GEOS-Chem version 13.4.1 and archived some reaction rates from the RxnRates Collection.

Screen Shot 2022-10-25 at 1 44 26 PM

The unit is per second (/s) which is different from what is shown on the GEOS-Chem website which is molec/cm3/s (http://wiki.seas.harvard.edu/geos-chem/index.php/History_collections_for_chemistry_and_photolysis#The_RxnRates_Collection). I am confused. Could anyone help me with this?

Thanks!!

yantosca commented 1 year ago

Thanks for writing @xin-tie. This should be fixed now: https://github.com/geoschem/geos-chem/issues/965

However, note that simulations using KPP (fullchem, Hg) will archive Prod/Loss diagnostics in molec/cm3/s as that is the units that come out of the solver. However for other specialty simulations (e.g. tagO3 etc) the units are kg/s.

xin-tie commented 1 year ago

Thanks for the reply, @yantosca. The issue fixed in #965 is for Prod/Loss diagnostics, not the RxnRates diagnostic. I used the freshly downloaded code and ran the model out of the box for one day (20190701). I used the most recent stable version of GEOS-Chem (14.0.0) by using command,

git checkout main

The unit in the RxnRates collection is still in per second (please see attached file). The name of the run directory is gc_4x5_47L_merra2_fullchem GEOSChem.RxnRates.20190701_0000z.nc4.zip

Correct me if wrong, but I believe this unit definition is in GC14.0.0/src/GEOS-Chem/Headers/state_diag_mod.F90, line 10972-10976 (https://github.com/geoschem/geos-chem/blob/53cac3b937064bdceed21acec721a157b6fb8b1e/Headers/state_diag_mod.F90), which still appears to be "s-1",

 ELSE IF ( TRIM( Name_AllCaps ) == 'RXNRATE' ) THEN
       IF ( isDesc    ) Desc  = 'KPP equation reaction rates'
       IF ( isUnits   ) Units = 's-1'
       IF ( isRank    ) Rank  = 3
       IF ( isTagged  ) TagId = 'RXN'

Is it possible that this is a typo introduced into the diagnostics code?

Thanks!

yantosca commented 1 year ago

Thanks for replying @xin-tie. In KPP, the reaction rates should be s-1 for the fullchem simulation. The wiki is wrong and I can fix it.

xin-tie commented 1 year ago

Thanks for the reply @yantosca. I think the unit of the Rxnrate collection is very likely molec / cm3 / s. This is how I came to this conclusion.

Take CH4 + OH as an example. For the same simulation mentioned above, I archived the reaction rate for CH4 + OH in RxnRates file (I have attached the file and the related variable is RxnRate_EQ025) and calculate its mean surface value which is 2.5e5.

GEOSChem.RxnRates.20190701_0000z.nc4.zip

Next, I calculated the reaction rate of CH4 + OH using the code below (all based on the present-day values), and the result is about 3e5 molec cm-3 s-1. The value from RxnRate file has the same magnitude as the reaction rate I calculated. So it appears the unit should be in molec cm-3 s-1 instead of s-1.

k_298K=6.3e-15 #second-order reaction rate (cm3 molec-1 s-1) of reaction CH4 + OH at temperature of 298K
                             #from JPL handbook
OH=1e6 #mole / cm3
CH4_mix_ratio=2000e-9 #ppb, present-day value 
Temp=298
kB=1.38e-23 #Boltzmann’s constant
P=1013e2 #Pa
n_air=P / (kB*Temp) #molec m-3 #ideal gas law
n_air_2=n_air / 1e6 #convert to molec cm-3
CH4=CH4_mix_ratio*n_air_2 #molec cm-3

R=k_298K*OH*CH4 #reaction rate, molec cm-3 s-1

 R
[1] 310373.5

This also makes sense based on how Rxnrate is calculated in KPP. Correct me if I am wrong, but I think it is calculated following the code in KPP/fullchem/gckpp_Function.F90. Here A should be in the unit of molec/cm3/s. 

SUBROUTINE Fun ( V, F, RCT, Vdot, Aout, Vdotout )

! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)
! F - Concentrations of fixed species (local)
  REAL(kind=dp) :: F(NFIX)
! RCT - Rate constants (local)
  REAL(kind=dp) :: RCT(NREACT)
! Vdot - Time derivative of variable species concentrations
  REAL(kind=dp) :: Vdot(NVAR)
! Aout - Optional argument to return equation rate constants
  REAL(kind=dp), OPTIONAL :: Aout(NREACT)
! Vdotout - Optional argument to return time derivative of variable species
  REAL(kind=dp), OPTIONAL :: Vdotout(NVAR)

! Computation of equation rates
  A(1) = RCT(1)*V(144)*V(269)*V(286)
  A(2) = RCT(2)*V(144)*V(271)
  A(3) = RCT(3)*V(144)*V(260)
  A(4) = RCT(4)*V(135)*V(269)*V(286)
  A(5) = RCT(5)*V(135)*V(271)
  A(6) = RCT(6)*V(135)*V(260)
  A(7) = RCT(7)*V(208)*V(286)
  A(8) = RCT(8)*V(269)*V(286)
  A(9) = RCT(9)*V(286)
  A(10) = RCT(10)*V(262)*V(286)
  A(11) = RCT(11)*V(85)
  A(12) = RCT(12)*V(85)*V(279)
  A(13) = RCT(13)*V(269)*V(280)
  A(14) = RCT(14)*V(269)*V(279)
  A(15) = RCT(15)*V(269)*V(281)
  A(16) = RCT(16)*V(269)*V(270)
  A(17) = RCT(17)*V(263)*V(269)
  A(18) = RCT(18)*V(279)*V(279)
  A(19) = RCT(19)*V(279)*V(279)
  A(20) = RCT(20)*V(279)*V(281)
  A(21) = RCT(21)*V(208)*V(279)
  A(22) = RCT(22)*V(280)*V(281)
  A(23) = RCT(23)*V(281)*V(281)
  A(24) = RCT(24)*V(258)*V(279)
  A(25) = RCT(25)*V(186)*V(279)
  .......

And these A values are then passed to RxnRate diagnostics in GeosCore/fullchem_mod.F90 

      !=====================================================================
       ! HISTORY (aka netCDF diagnostics)
       !
       ! Archive KPP reaction rates [s-1]
       ! See gckpp_Monitor.F90 for a list of chemical reactions
       !
       ! NOTE: In KPP 2.5.0+, VAR and FIX are now private to the integrator
       ! and point to C.  Therefore, pass C(1:NVAR) instead of VAR and
       ! C(NVAR+1:NSPEC) instead of FIX to routine FUN.
       !=====================================================================
       IF ( State_Diag%Archive_RxnRate ) THEN
          !---------------------------------------------------
          ! Get equation rates (Aout)
          !---------------------------------------------------
          CALL Fun( V       = C(1:NVAR),                                     &
                    F       = C(NVAR+1:NSPEC),                               &
                    RCT     = RCONST,                                        &
                    Vdot    = Vloc,                                          &
                    Aout    = Aout                                          )

          DO S = 1, State_Diag%Map_RxnRate%nSlots
             N = State_Diag%Map_RxnRate%slot2Id(S)
             State_Diag%RxnRate(I,J,L,S) = Aout(N)
          ENDDO
       ENDIF

I hope this is helpful.

yantosca commented 1 year ago

Thanks for the feedback @xin-tie. I am also tagging @msl3v who may know more about this.

Also tagging @Jourdan-He and @SaptSinha.

msl3v commented 1 year ago

Hey all. I believe all the rates coming out of KPP to be mcl/cm3/s

yantosca commented 1 year ago

Thanks @msl3v. I'll change the units in the code for 14.0.2.

xin-tie commented 1 year ago

Thanks @yantosca @msl3v !!!

yantosca commented 1 year ago

This will be merged into the 14.0.2 development stream (see PR #1483)

yantosca commented 1 year ago

PR #1483 has now been merged and this update will ship with GEOS-Chem 14.0.2. We can now close out this issue.