Skip to content

Issue with units of O3 loss in tagged simulation #2094

Open
@ncolombi

Description

Hello,

@xpye98 and I have previously reported issues with the tagged O3 run in GEOS-Chem version 13 and 14. The fix in #1109 corrected the issues involved with unit conversion from archived Prod/Loss of O3 (units of molec/cm3/s). P(O3) and L(O3) are now both converted from molec/cm3/s to units of kilograms (kg) in tagged_o3_mod.F90 and concentration of O3 is calculated in lines 397 and 461 as shown below:

   ! Add P(O3) [kg] to the total O3 species
   Spc(id_O3)%Conc(I,J,L) = Spc(id_O3)%Conc(I,J,L) + PO3_kg

   ! Apply chemical loss [kg]
   Spc(N)%Conc(I,J,L) = Spc(N)%Conc(I,J,L) - LO3_kg

In the second equation, where each term is in units of kilograms, as long as LO3_kg is greater than PO3_kg at a given timestep, the SpcConc will have a value of zero, which is what I see in my 14.2.0 simulation for all tagged species (zero everywhere). @xpye98 has corrected this issue in his v13.2.1 run by modifying the tagged_o3_mod.F90 file as such:

   ! Apply chemical loss [kg], LO3 is [1], Spc is [kg]
   Spc(I,J,L,N) = Spc(I,J,L,N) - LO3 * Spc(I,J,L,N)

In this case, LO3 has been changed to units of [1] and multiplied by SpcConc [kg], which fixes the issue. Since the SpcConc has an initial value of 0, this means that P(O3) will always be greater than LO3*SpcConc, and SpcConc will slowly reach a steady state concentration after some spin up. Outputting LO3 in units of [1] rather than [kg] requires dividing L(O3) by Ox mass in the fullchem_mod.F90 file. Is there anyway we can incorporate this fix in the next release of version 14?

Activity

yantosca

yantosca commented on Jan 2, 2024

@yantosca
Contributor

Thanks for the reminder @ncolombi! I'll work on getting this into the no-diff-to-benchmark branch, so that it can go into version 14.

self-assigned this
on Jan 2, 2024
ncolombi

ncolombi commented on Jan 2, 2024

@ncolombi
Author

Thanks Bob!

xpye98

xpye98 commented on Jan 2, 2024

@xpye98
Contributor

Thanks for tagging me. I believe previous GC (before 12.0.0 or so) used this kind of different units for Ox_prod and Ox_loss, which could produce the correct result. #973 also reported similar bug and fixed it by outputting different unit of Ox_loss.

yantosca

yantosca commented on Jan 2, 2024

@yantosca
Contributor

Question for you @xpye98 @ncolombi: In GeosCore/fullchem_mod.F90 the loss is computed from KPP as:

       !====================================================================
       ! HISTORY (aka netCDF diagnostics)
       !
       ! Prod and loss of families or species [molec/cm3/s]
       !
       ! NOTE: KppId is the KPP ID # for each of the prod and loss
       ! diagnostic species.  This is the value used to index the
       ! KPP "C" array (in module gckpp_Global.F90).
       !====================================================================

       ! Chemical loss of species or families [molec/cm3/s]
       IF ( State_Diag%Archive_Loss ) THEN
          DO S = 1, State_Diag%Map_Loss%nSlots
             KppId = State_Diag%Map_Loss%slot2Id(S)
             State_Diag%Loss(I,J,L,S) = C(KppID) / DT
          ENDDO
       ENDIF

for each of the loss families that are defined in KPP/fullchem/fullchem.kpp:

https://github.com/geoschem/geos-chem/blob/3cd802eb4a9b7a226254fc3938778c8db8c07cf2/KPP/fullchem/fullchem.kpp#L13C1-L20C14

that is LOx, LCO, LCH4. These are in molec/cm3, so we divide by DT to get molec/cm3/s.

Because LOx is a family species, is it appropriate to divide it by the concentration of O3? That would probably be approximate but am wondering if we'd have to also take into account the other species in the LOx family into the computation.

We could also add a LO3 family to KPP, which would just be the loss of the species O3, and then divide by the concentration of O3 to get the loss rate in the units that you want. Again, I'm not sure if that would be scientifically appropriate for the purposes of the tagged O3 simulation.

Also I'm thinking to add this to a different diagnostic, as the Prod_LOx might be used as input to other simulations, and thus if we change the units we would break those simulations.

ncolombi

ncolombi commented on Jan 2, 2024

@ncolombi
Author

Hi Bob, I believe you would need to divide LOx by the mass of the Ox family (not just the mass of O3). @xpye98 does that in his flexchem_mod.F90 (from version 13) in lines 1111-1289:

   !====================================================================
   ! HISTORY (aka netCDF diagnostics)
   !
   ! Prod and loss of families or species [molec/cm3/s]
   !
   ! NOTE: KppId is the KPP ID # for each of the prod and loss
   ! diagnostic species.  This is the value used to index the
   ! KPP "VAR" array (in module gckpp_Global.F90).
   !====================================================================

   ! Chemical loss of species or families [molec/cm3/s]
   ! Change unit to [1/m3/s] for TagO3 simulation
   ! xpye, 2022/02/12
   IF ( State_Diag%Archive_Loss ) THEN
      id_O3 = Ind_('O3')
      id_NO2 = Ind_('NO2')
      id_NO3 = Ind_('NO3')
      id_PAN = Ind_('PAN')
      id_PPN = Ind_('PPN')
      id_MPAN = Ind_('MPAN')
      id_HNO4 = Ind_('HNO4')
      id_N2O5 = Ind_('N2O5')
      id_HNO3 = Ind_('HNO3')
      id_BrO = Ind_('BrO')
      id_HOBr = Ind_('HOBr')
      id_BrNO2 = Ind_('BrNO2')
      id_BrNO3 = Ind_('BrNO3')
      id_MPN = Ind_('MPN')
      id_ETHLN = Ind_('ETHLN')
      id_MVKN = Ind_('MVKN')
      id_MCRHN = Ind_('MCRHN')
      id_MCRHNB = Ind_('MCRHNB')
      id_PROPNN = Ind_('PROPNN')
      id_R4N2 = Ind_('R4N2')
      id_PRN1 = Ind_('PRN1')
      id_PRPN = Ind_('PRPN')
      id_R4N1 = Ind_('R4N1')
      id_HONIT = Ind_('HONIT')
      id_MONITS = Ind_('MONITS')
      id_MONITU = Ind_('MONITU')
      id_OLND = Ind_('OLND')
      id_OLNN = Ind_('OLNN')
      id_IHN1 = Ind_('IHN1')
      id_IHN2 = Ind_('IHN2')
      id_IHN3 = Ind_('IHN3')
      id_IHN4 = Ind_('IHN4')
      id_INPB = Ind_('INPB')
      id_INPD = Ind_('INPD')
      id_ICN = Ind_('ICN')
      id_IDN = Ind_('IDN')
      id_ITCN = Ind_('ITCN')
      id_ITHN = Ind_('ITHN')
      id_ISOPNOO1 = Ind_('ISOPNOO1')
      id_ISOPNOO2 = Ind_('ISOPNOO2')
      id_INO2B = Ind_('INO2B')
      id_INO2D = Ind_('INO2D')
      id_INA = Ind_('INA')
      id_IDHNBOO = Ind_('IDHNBOO')
      id_IDHNDOO1 = Ind_('IDHNDOO1')
      id_IDHNDOO2 = Ind_('IDHNDOO2')
      id_IHPNBOO = Ind_('IHPNBOO')
      id_IHPNDOO = Ind_('IHPNDOO')
      id_ICNOO = Ind_('ICNOO')
      id_IDNOO = Ind_('IDNOO')
      id_MACRNO2 = Ind_('MACRNO2')
      id_ClO = Ind_('ClO')
      id_HOCl = Ind_('HOCl')
      id_ClNO2 = Ind_('ClNO2')
      id_ClNO3 = Ind_('ClNO3')
      id_Cl2O2 = Ind_('Cl2O2')
      id_OClO = Ind_('OClO')
      id_O = Ind_('O')
      id_O1D = Ind_('O1D')
      id_IO = Ind_('IO')
      id_HOI = Ind_('HOI')
      id_IONO = Ind_('IONO')
      id_IONO2 = Ind_('IONO2')
      id_OIO = Ind_('OIO')
      id_I2O2 = Ind_('I2O2')
      id_I2O3 = Ind_('I2O3')
      id_I2O4 = Ind_('I2O4')

     ! Error check
      IF ( id_O3 <= 0 ) THEN
        PRINT *, "Error"
     ENDIF

     ! Get ozone molecular weight from species database
      O3_MW_g  = State_Chm%SpcData(id_O3)%Info%MW_g  ! g/mol

      DO S = 1, State_Diag%Map_Loss%nSlots
         KppId = State_Diag%Map_Loss%slot2Id(S)
         State_Diag%Loss(I,J,L,S) = VAR(KppID) / DT
         ! Calculate Ox mass, according to Ox family in KPP 
         Ox_Mass        = State_Chm%Species(I, J, L, id_O3     ) &  
                        + State_Chm%Species(I, J, L, id_NO2    ) &
                      + State_Chm%Species(I, J, L, id_NO3    )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_PAN    ) &
                      + State_Chm%Species(I, J, L, id_PPN    ) &
                      + State_Chm%Species(I, J, L, id_MPAN    ) &
                      + State_Chm%Species(I, J, L, id_HNO4   ) &
                      + State_Chm%Species(I, J, L, id_N2O5   )*3.0_fp &
                      + State_Chm%Species(I, J, L, id_HNO3   ) &
                      + State_Chm%Species(I, J, L, id_BrO    ) &
                      + State_Chm%Species(I, J, L, id_HOBr   ) &
                      + State_Chm%Species(I, J, L, id_BrNO2  ) &
                      + State_Chm%Species(I, J, L, id_BrNO3  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_MPN    ) &
                      + State_Chm%Species(I, J, L, id_ETHLN  ) &
                      + State_Chm%Species(I, J, L, id_MVKN  ) &
                      + State_Chm%Species(I, J, L, id_MCRHN  ) &
                      + State_Chm%Species(I, J, L, id_MCRHNB  ) &
                      + State_Chm%Species(I, J, L, id_PROPNN  ) &
                      + State_Chm%Species(I, J, L, id_R4N2   ) &
                      + State_Chm%Species(I, J, L, id_PRN1   ) &
                      + State_Chm%Species(I, J, L, id_PRPN   ) &
                      + State_Chm%Species(I, J, L, id_R4N1   ) &
                      + State_Chm%Species(I, J, L, id_HONIT   ) &
                      + State_Chm%Species(I, J, L, id_MONITS   ) &
                      + State_Chm%Species(I, J, L, id_MONITU   ) &
                      + State_Chm%Species(I, J, L, id_OLND   ) &
                      + State_Chm%Species(I, J, L, id_OLNN   ) &
                      + State_Chm%Species(I, J, L, id_IHN1   ) &
                      + State_Chm%Species(I, J, L, id_IHN2   ) &
                      + State_Chm%Species(I, J, L, id_IHN3   ) &
                      + State_Chm%Species(I, J, L, id_IHN4   ) &
                      + State_Chm%Species(I, J, L, id_INPB   ) &
                      + State_Chm%Species(I, J, L, id_INPD   ) &
                      + State_Chm%Species(I, J, L, id_ICN  ) &
                      + State_Chm%Species(I, J, L, id_IDN  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_ITCN  ) &
                      + State_Chm%Species(I, J, L, id_ITHN  ) &
                      + State_Chm%Species(I, J, L, id_ISOPNOO1  ) &
                      + State_Chm%Species(I, J, L, id_ISOPNOO2  ) &
                      + State_Chm%Species(I, J, L, id_INO2B  ) &
                      + State_Chm%Species(I, J, L, id_INO2D  ) &
                      + State_Chm%Species(I, J, L, id_INA  ) &
                      + State_Chm%Species(I, J, L, id_IDHNBOO  ) &
                      + State_Chm%Species(I, J, L, id_IDHNDOO1  ) &
                      + State_Chm%Species(I, J, L, id_IDHNDOO2 ) &  
                      + State_Chm%Species(I, J, L, id_IHPNBOO ) &  
                      + State_Chm%Species(I, J, L, id_IHPNDOO ) &  
                      + State_Chm%Species(I, J, L, id_ICNOO   ) &  
                      + State_Chm%Species(I, J, L, id_IDNOO   ) *2.0_fp&  
                      + State_Chm%Species(I, J, L, id_MACRNO2 ) &  
                      + State_Chm%Species(I, J, L, id_ClO     ) &  
                      + State_Chm%Species(I, J, L, id_HOCl     ) &  
                      + State_Chm%Species(I, J, L, id_ClNO2     ) &  
                      + State_Chm%Species(I, J, L, id_ClNO3  )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_Cl2O2  )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_OClO   )*2.0_fp&
                      + State_Chm%Species(I, J, L, id_O      ) &
                      + State_Chm%Species(I, J, L, id_O1D    ) &
                      + State_Chm%Species(I, J, L, id_IO     ) &
                      + State_Chm%Species(I, J, L, id_HOI    ) &
                      + State_Chm%Species(I, J, L, id_IONO   ) &
                      + State_Chm%Species(I, J, L, id_IONO2  )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_OIO    )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O2   )*2.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O3   )*3.0_fp &
                      + State_Chm%Species(I, J, L, id_I2O4   )*4.0_fp 
         ! Grid box volume [cm3]
         BOXVL          = State_Met%AIRVOL(I, J, L) * 1e+6_fp
         ! Convert Ox_Mass from [molec/cm3] to [molec]
         Ox_Mass        = Ox_Mass*BOXVL
         ! Divide L(Ox) [molec/cm3/s] by Ox mass [molec] 
         ! and then convert to [1/m3/s] for HEMCO
         IF ( Ox_Mass > 0.0_fp ) THEN
                 State_Diag%Loss(I, J, L, S)  = &
           ( State_Diag%Loss(I, J, L, S) / Ox_Mass ) *  1.0e+6_fp
         ENDIF
      ENDDO
   ENDIF

   ! Chemical production of species or families [molec/cm3/s]
   ! Change unit to [kg/m3/s], for TagO3 simulation
   ! xpye, 2022/02/12
   IF ( State_Diag%Archive_Prod ) THEN
      DO S = 1, State_Diag%Map_Prod%nSlots
         KppID = State_Diag%Map_Prod%slot2Id(S)
         State_Diag%Prod(I,J,L,S) = VAR(KppID) / DT

         ! Convert P(Ox) from [molec/cm3/s] to [kg/cm3/s]
         State_Diag%Prod(I, J, L, S) = State_Diag%Prod(I, J, L, S) &
                                       / AVO*1.e-3_fp * O3_MW_g
         ! Convert to [kg/m3/s], for HEMCO 
         State_Diag%Prod(I, J, L, S)     = State_Diag%Prod(I, J, L, S)* 1e+6_fp
      ENDDO
   ENDIF

@xpye98, let me know if I misinterpreted anything here.

yantosca

yantosca commented on Jan 2, 2024

@yantosca
Contributor

@xpye98 e98 @ncolombi: Thanks. I'll see how I can implement this. A few things have changed since the 13.1 code, especially the update to KPP 3.0.0. Also, I think we can also use C(ind_O3) instead of Spc(id_O3)%Conc(I,J,L), etc. to simplify the code somewhat. Stay tuned.

xpye98

xpye98 commented on Jan 2, 2024

@xpye98
Contributor

@yantosca @ncolombi
Yes, in my practice (where I got the correct tag results), I divide it by the mass of Ox family, not just O3, but note that the definitions of Ox family are different for different GC versions.
I agree with you that changing the unit would affect other analysis, maybe outputting Ox_loss with both units (one for consistency with Ox_pord, and the other for tagged O3) would be good.

Outputting just PO3 and LO3 (not Ox family) to KPP won't work. Actually we've tried this before and got very strange results. The primary formation pathway of O3 is NO2 -> O3 + NO and the primary loss is O3 photolysis. Considering only O3 will get very high PO3 and LO3, and somehow will break the tagO3 simulation. Scientifically, it is also meaningful to consider Ox family to avoid fast reactions between O3 and NOx. But this is just my opinion and might need further test.

xpye98

xpye98 commented on Jan 2, 2024

@xpye98
Contributor

@yantosca Also, I'm more curious of why at the very beginning we used different units for prod and loss in previous GC versions. Is this the only way we can tag correctly?

added and removed on Jan 3, 2024

17 remaining items

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Assignees

Labels

category: BugSomething isn't workingnever staleNever label this issue as staletopic: DiagnosticsRelated to output diagnostic datatopic: Math and UnitsIssues pertaining to math expressions and/or unit conversionstopic: TagO3 SimulationsRelated to the GEOS-Chem Tagged O3 simulation

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions

    Issue with units of O3 loss in tagged simulation · Issue #2094 · geoschem/geos-chem