From 2b420f6fa0925595b86b9b7820d5e1e99a8aed93 Mon Sep 17 00:00:00 2001 From: Jonathan Moch Date: Tue, 22 Sep 2020 00:01:35 -0400 Subject: [PATCH 01/16] HMS chemistry (ref: Moch et al., 2020, JGR) --- GeosCore/aerosol_mod.F90 | 37 ++- GeosCore/input_mod.F90 | 1 + GeosCore/isorropiaII_mod.F90 | 16 +- GeosCore/sulfate_mod.F90 | 491 +++++++++++++++++++++++++++++++++-- Headers/species_database.yml | 15 ++ Headers/state_diag_mod.F90 | 162 +++++++++++- 6 files changed, 695 insertions(+), 27 deletions(-) diff --git a/GeosCore/aerosol_mod.F90 b/GeosCore/aerosol_mod.F90 index af0bbd5dd..709c076ab 100644 --- a/GeosCore/aerosol_mod.F90 +++ b/GeosCore/aerosol_mod.F90 @@ -43,6 +43,7 @@ MODULE AEROSOL_MOD ! SALC : Coarse mode seasalt aerosol [kg/m3] ! SO4_NH4_NIT : Lumped SO4-NH4-NIT aerosol [kg/m3] ! SO4 : Sulfate aerosol [kg/m3] + ! HMS : Hydroxymethane sulfonate aerosol [kg/m3] ! (jmm, 06/29/18) ! NH4 : Ammonium aerosol [kg/m3] ! NIT : Inorganic nitrate aerosol [kg/m3] ! SOILDUST : Mineral dust aerosol from soils [kg/m3] @@ -69,6 +70,7 @@ MODULE AEROSOL_MOD REAL(fp), ALLOCATABLE, PUBLIC :: SALC(:,:,:) REAL(fp), ALLOCATABLE, PUBLIC :: SO4_NH4_NIT(:,:,:) REAL(fp), ALLOCATABLE, PUBLIC :: SO4(:,:,:) + REAL(fp), ALLOCATABLE, PUBLIC :: HMS(:,:,:) ! (jmm, 06/29/18) REAL(fp), ALLOCATABLE, PUBLIC :: NH4(:,:,:) REAL(fp), ALLOCATABLE, PUBLIC :: NIT(:,:,:) REAL(fp), ALLOCATABLE, PUBLIC :: FRAC_SNA(:,:,:,:) @@ -130,7 +132,7 @@ MODULE AEROSOL_MOD INTEGER :: id_POA1, id_POA2, id_OPOA1, id_OPOA2 INTEGER :: id_TSOA1, id_TSOA2, id_TSOA3, id_TSOA0 INTEGER :: id_ASOAN, id_ASOA1, id_ASOA2, id_ASOA3 - INTEGER :: id_DUST1, id_SOAS, id_SALACL + INTEGER :: id_DUST1, id_SOAS, id_SALACL, id_HMS ! (jmm, 06/29/18) INTEGER :: id_SOAGX, id_SOAIE INTEGER :: id_INDIOL,id_LVOCOA @@ -233,7 +235,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & LOGICAL :: LUCX LOGICAL :: IS_OCPO, IS_OCPI, IS_BC LOGICAL :: IS_SO4, IS_NH4, IS_NIT - LOGICAL :: IS_SAL, IS_DST + LOGICAL :: IS_SAL, IS_DST, IS_HMS ! (jmm, 06/29/18) LOGICAL :: IS_TSOA, IS_ASOA LOGICAL :: IS_POA, IS_OPOA LOGICAL :: IS_SOAGX @@ -285,6 +287,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & IS_OCPO = ( id_OCPO > 0 ) IS_BC = ( id_BCPI > 0 .AND. id_BCPO > 0 ) IS_SO4 = ( id_SO4 > 0 ) + IS_HMS = ( id_HMS > 0 ) !(jmm, 06/29/18) IS_NH4 = ( id_NH4 > 0 ) IS_NIT = ( id_NIT > 0 ) IS_DST = ( id_DST1 > 0 .AND. id_DST2 > 0 ) @@ -465,12 +468,14 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! way (size and optics) as all other sulfate aerosol (DAR ! 2013) SO4_NH4_NIT(I,J,L) = ( Spc(I,J,L,id_SO4) + & + Spc(I,J,L,id_HMS) + & ! (jmm, 07/2/18) Spc(I,J,L,id_NH4) + & Spc(I,J,L,id_NIT)) / & !Spc(I,J,L,id_SO4s) + & !Spc(I,J,L,id_NITs) ) / & AIRVOL(I,J,L) SO4(I,J,L) = Spc(I,J,L,id_SO4) / AIRVOL(I,J,L) + HMS(I,J,L) = Spc(I,J,L,id_HMS) / AIRVOL(I,J,L) ! (jmm, 06/30/18) NH4(I,J,L) = Spc(I,J,L,id_NH4) / AIRVOL(I,J,L) NIT(I,J,L) = Spc(I,J,L,id_NIT) / AIRVOL(I,J,L) SLA(I,J,L) = 0e+0_fp @@ -481,6 +486,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! Tropospheric sulfate is zero in stratosphere SO4_NH4_NIT(I,J,L) = 0e+0_fp SO4(I,J,L) = 0e+0_fp + HMS(I,J,L) = 0e+0_fp ! (jmm, 06/30/18) NH4(I,J,L) = 0e+0_fp NIT(I,J,L) = 0e+0_fp SLA(I,J,L) = KG_STRAT_AER(I,J,L,1) / AIRVOL(I,J,L) @@ -505,10 +511,12 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! should be considered as a part of sea-salt (XNW Dec 7 2017) SO4_NH4_NIT(I,J,L) = ( Spc(I,J,L,id_SO4) + & + Spc(I,J,L,id_HMS) + & ! (jmm, 07/02/18) Spc(I,J,L,id_NH4) + & Spc(I,J,L,id_NIT)) / & AIRVOL(I,J,L) SO4(I,J,L) = Spc(I,J,L,id_SO4) / AIRVOL(I,J,L) + HMS(I,J,L) = Spc(I,J,L,id_HMS) / AIRVOL(I,J,L) !(jmm, 06/30/18) NH4(I,J,L) = Spc(I,J,L,id_NH4) / AIRVOL(I,J,L) NIT(I,J,L) = Spc(I,J,L,id_NIT) / AIRVOL(I,J,L) @@ -519,7 +527,8 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! Save these fractions for partitioning of optics ! until later when these may be treated independently - FRAC_SNA(I,J,L,1) = ( ( Spc(I,J,L,id_SO4 ) ) & + FRAC_SNA(I,J,L,1) = ( ( Spc(I,J,L,id_SO4 ) + & + Spc(I,J,L,id_HMS ) ) & ! (jmm, 07/02/18) ! & + Spc(I,J,L,id_SO4s) ) / AIRVOL(I,J,L) ) & / SO4_NH4_NIT(I,J,L) @@ -854,7 +863,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! ! Compute PM2.5 concentration [kg/m3] ! - ! PM25 = 1.33 (NH4 + NIT + SO4) + BCPI + BCPO + + ! PM25 = 1.33 (NH4 + NIT + SO4 + HMS) + BCPI + BCPO + ! 2.10 (OCPO + 1.16 OCPI) + 1.16 SOA* + ! DST1 + 0.38 DST2 + 1.86 SALA ! @@ -864,6 +873,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! NOTES: ! - We apply growth factors at 35% RH (computed above): ! 1.33 for SO4, NIT, and NH4 + ! 1.33 for HMS (jmm, 06/30/18) ! 1.16 for OCPI and SOA ! 1.86 for SALA ! - Ratio of OM/OC = 2.1 is applied to OCPI and OCPO above @@ -882,6 +892,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & PM25(I,J,L) = NH4(I,J,L) * SIA_GROWTH + & NIT(I,J,L) * SIA_GROWTH + & SO4(I,J,L) * SIA_GROWTH + & + HMS(I,J,L) * SIA_GROWTH + & ! (jmm, 06/30/18) BCPI(I,J,L) + & BCPO(I,J,L) + & OCPO(I,J,L) + & @@ -2293,6 +2304,7 @@ SUBROUTINE Init_Aerosol( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) id_SALACL = Ind_( 'SALACL' ) id_SO4 = Ind_( 'SO4' ) id_SO4s = Ind_( 'SO4s' ) + id_HMS = Ind_( 'HMS' ) ! (jmm, 06/30/18) id_NITs = Ind_( 'NITs' ) id_POA1 = Ind_( 'POA1' ) id_POA2 = Ind_( 'POA2' ) @@ -2364,6 +2376,12 @@ SUBROUTINE Init_Aerosol( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) IF ( RC /= GC_SUCCESS ) RETURN SO4 = 0.0_fp + ! (jmm, 06/30/18) + ALLOCATE( HMS( State_Grid%NX, State_Grid%NY, State_Grid%NZ ), STAT=RC ) + CALL GC_CheckVar( 'aerosol_mod.F:HMS', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + HMS = 0.0_fp + ALLOCATE( NH4( State_Grid%NX, State_Grid%NY, State_Grid%NZ ), STAT=RC ) CALL GC_CheckVar( 'aerosol_mod.F:NH4', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN @@ -2493,6 +2511,7 @@ SUBROUTINE CLEANUP_AEROSOL IF ( ALLOCATED( ACL ) ) DEALLOCATE( ACL ) IF ( ALLOCATED( SO4_NH4_NIT ) ) DEALLOCATE( SO4_NH4_NIT ) IF ( ALLOCATED( SO4 ) ) DEALLOCATE( SO4 ) + IF ( ALLOCATED( HMS ) ) DEALLOCATE( HMS ) ! (jmm, 06/30/18) IF ( ALLOCATED( NH4 ) ) DEALLOCATE( NH4 ) IF ( ALLOCATED( NIT ) ) DEALLOCATE( NIT ) IF ( ALLOCATED( FRAC_SNA ) ) DEALLOCATE( FRAC_SNA ) @@ -2794,6 +2813,16 @@ SUBROUTINE Set_AerMass_Diagnostic( Input_Opt, State_Chm, State_Diag, & State_Diag%AerMassSO4(I,J,L) = SO4(I,J,L) * kgm3_to_ugm3 ENDIF + !-------------------------------------- + ! AerMassHMS [ug/m3] + ! jmm 3/6/19 + !-------------------------------------- + IF ( State_Diag%Archive_AerMassHMS ) THEN + State_Diag%AerMassHMS(I,J,L) = HMS(I,J,L) * & + kgm3_to_ugm3 + ENDIF + + !-------------------------------------- ! AerMassSOAGX [ug/m3] !-------------------------------------- diff --git a/GeosCore/input_mod.F90 b/GeosCore/input_mod.F90 index 1dbd20cf1..845f6de35 100644 --- a/GeosCore/input_mod.F90 +++ b/GeosCore/input_mod.F90 @@ -5950,6 +5950,7 @@ SUBROUTINE Do_Error_Checks( Input_Opt, RC ) MAX( Ind_('SO2' ,'A'), 0 ) + & MAX( Ind_('SO4' ,'A'), 0 ) + & MAX( Ind_('SO4s','A'), 0 ) + & + MAX( Ind_('HMS' ,'A'), 0 ) + &! (jmm, 07/2/18) MAX( Ind_('MSA' ,'A'), 0 ) + & MAX( Ind_('NH3' ,'A'), 0 ) + & MAX( Ind_('NH4' ,'A'), 0 ) + & diff --git a/GeosCore/isorropiaII_mod.F90 b/GeosCore/isorropiaII_mod.F90 index 6ac9bc3e7..09e775e01 100644 --- a/GeosCore/isorropiaII_mod.F90 +++ b/GeosCore/isorropiaII_mod.F90 @@ -180,6 +180,7 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & LOGICAL, SAVE :: FIRST = .TRUE. INTEGER, SAVE :: id_HNO3, id_NH3, id_NH4 INTEGER, SAVE :: id_NIT, id_SALA, id_SO4 + INTEGER, SAVE :: id_HMS ! jmm 12/5/18 INTEGER, SAVE :: id_SALACL, id_HCL, id_SALCCL INTEGER, SAVE :: id_SO4s, id_NITs, id_SALC INTEGER, SAVE :: id_SALAAL, id_SALCAL @@ -279,6 +280,7 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & id_NIT = Ind_('NIT' ) id_SALA = Ind_('SALA') id_SO4 = Ind_('SO4' ) + id_HMS = Ind_('HMS' ) id_SALACL = Ind_('SALACL') id_HCL = Ind_('HCl') id_SALC = Ind_('SALC') @@ -295,6 +297,11 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & CALL GC_Error( ErrMsg, RC, ThisLoc ) RETURN ENDIF + IF ( id_HMS <= 0 ) THEN + ErrMsg = 'HMS is an undefined species!' + CALL GC_Error( ErrMsg, RC, ThisLoc ) + RETURN + ENDIF IF ( id_NH3 <= 0 ) THEN ErrMsg = 'NH3 is an undefined species!' CALL GC_Error( ErrMsg, RC, ThisLoc ) @@ -545,7 +552,9 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & ! Total SO4 [mole/m3], also consider SO4s in SALA TSO4 = (Spc(I,J,L,id_SO4)+Spc(I,J,L,id_SALA)*0.08*AlkR) & - * 1.e+3_fp / ( 96.e+0_fp * VOL) + * 1.e+3_fp / ( 96.e+0_fp * VOL) + & + Spc(I,J,L,id_HMS) * 0.5e+3_fp / ( 111.e+0_fp * VOL) + ! Total NH3 [mole/m3] TNH3 = Spc(I,J,L,id_NH4)*1.e+3_fp/(18.e+0_fp*VOL) + & @@ -557,7 +566,7 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & TSO4 = Spc(I,J,L,id_SO4s) & * 1.e+3_fp * AlkR / (31.4e+0_fp * VOL) + & Spc(I,J,L,id_SALC)*0.08_fp & - * 1.e+3_fp * AlkR / (96.e+0_fp * VOL) + * 1.e+3_fp * AlkR / (96.e+0_fp * VOL) ! Total NH3 [mole/m3] !TNH3 = Spc(I,J,L,id_NH4s)*1.e+3_fp*AlkR/(31.4e+0_fp*VOL)+ & @@ -896,7 +905,8 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & Spc(I,J,L,id_NH4) /18e+0_fp + & Spc(I,J,L,id_SALA)*0.3061e+0_fp/23.0e+0_fp ) - DEN_SAV = ( Spc(I,J,L,id_SO4) / 96e+0_fp*2e+0_fp + & + DEN_SAV = ( Spc(I,J,L,id_SO4) / 96e+0_fp*2e+0_fp + & + Spc(I,J,L,id_HMS) / 111e+0_fp + & Spc(I,J,L,id_NIT) / 62e+0_fp + & HNO3_DEN / 63e+0_fp + & Spc(I,J,L,id_SALA) *0.55e+0_fp / 35.45e+0_fp) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index b3afc9b36..4211682f8 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -106,6 +106,9 @@ MODULE SULFATE_MOD ! PMSA_DMS : P(MSA) from DMS [v/v/timestep] ! PSO2_DMS : P(SO2) from DMS [v/v/timestep] ! PSO4_SO2 : P(SO4) from SO2 [v/v/timestep] + ! PHMS_SO2 : P(HMS) from SO2 [v/v/timestep] (jmm, 6/15/18) + ! PSO2_HMS : P(SO2) from HMS [v/v/timestep] (jmm, 6/15/18) + ! PSO4_HMS : P(SO4) from HMS & radical chem [v/v/timestep] (jmm, 6/26/18) ! SSTEMP : Sea surface temperatures [K] ! Eev : SO2 em. from eruptive volcanoes [kg SO2/box/s] ! Env : SO2 em. from non-erup volcanoes [kg SO2/box/s] @@ -125,6 +128,9 @@ MODULE SULFATE_MOD REAL(fp), ALLOCATABLE :: PSO2_DMS(:,:,:) REAL(fp), ALLOCATABLE :: PSO4_SO2(:,:,:) REAL(fp), ALLOCATABLE :: PSO4_SS(:,:,:) + REAL(fp), ALLOCATABLE :: PHMS_SO2(:,:,:) ! jmm 06/13/2018 + REAL(fp), ALLOCATABLE :: PSO2_HMS(:,:,:) ! jmm 06/13/2018 + REAL(fp), ALLOCATABLE :: PSO4_HMS(:,:,:) ! jmm 06/26/2018 REAL(fp), ALLOCATABLE :: PNITs(:,:,:) REAL(fp), ALLOCATABLE :: PNIT_dust(:,:,:,:) ! tdf REAL(fp), ALLOCATABLE :: PSO4_dust(:,:,:,:) ! tdf @@ -181,6 +187,7 @@ MODULE SULFATE_MOD INTEGER :: id_HOBr, id_SO4H1, id_SO4H2 INTEGER :: id_HOCl, id_SO4H3, id_SO4H4 INTEGER :: id_HCOOH, id_ACTA !jmm 1/31/19 + INTEGER :: id_HMS, id_CH2O ! jmm 06/13/2018 ! Species drydep ID flags @@ -389,6 +396,9 @@ SUBROUTINE CHEMSULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, & PSO2_DMS = 0e+0_fp PMSA_DMS = 0e+0_fp PSO4_SO2 = 0e+0_fp + PHMS_SO2 = 0e+0_fp ! jmm 06/13/2018 + PSO2_HMS = 0e+0_fp ! jmm 06/13/2018 + PSO4_HMS = 0e+0_fp ! jmm 06/26/2018 PSO4_SS = 0e+0_fp PNITs = 0e+0_fp PSO4_dust = 0e+0_fp ! tdf 04/17/08 @@ -2332,12 +2342,15 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Reaction List (by Rokjin Park, rjp@io.harvard.edu) ! ============================================================================ ! (1 ) SO2 production: -! DMS + OH, DMS + NO3 (saved in CHEM_DMS) +! (a) DMS + OH, DMS + NO3 (saved in CHEM_DMS) +! (b) HMS -> SO2 + HCHO (aq) ! . ! (2 ) SO2 loss: ! (a) SO2 + OH -> SO4 ! (b) SO2 -> drydep ! (c) SO2 + H2O2 or O3 (aq) -> SO4 +! (d) SO2 + HCHO (aq)-> HMS +! (d) SO2 + HMS -> 2 SO4 ! . ! (3 ) SO2 = SO2_0 * exp(-bt) + PSO2_DMS/bt * [1-exp(-bt)] ! . @@ -2388,6 +2401,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fp) :: LSTOT, ALKdst, ALKss, ALKds, NH3, CL, TNA REAL(fp) :: SSCvv, aSO4, SO2_sr, SR, TANIT REAL(fp) :: TFA, TAA, TDCA ! (jmm, 12/03/2018) + REAL(fp) :: HCHO0, HMSc, KaqHCHO, KaqHMS ! (jmm, 06/07/2018) + REAL(fp) :: L7, L7S, L7_b, L7S_b, HMS0 ! (jmm, 06/07/2018) + REAL(fp) :: L8, L8S, OH0, KaqHMS2 ! (jmm, 06/26/2018) REAL(fp) :: PSO4d_tot, PNITd_tot REAL(fp) :: SO2_gas, PH2SO4d_tot REAL(fp) :: H2SO4_cd, H2SO4_gas @@ -2399,7 +2415,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fp) :: SO4H1_vv, SO4H2_vv, LSTOT0 REAL(fp) :: SO2_ss0, rSIV, fupdateHOBr_0 REAL(fp) :: HCO3, HCHOBr, KO3, KHOBr, f_srhobr, HOBr0 - REAL(fp) :: TMP + REAL(fp) :: TMP, LSTOT_HMS ! (jmm, 06/15/18) REAL(fp) :: KaqO2, L4, L4S, MnII, FeIII REAL(fp) :: DUST, Mn_ant, Mn_nat @@ -2557,6 +2573,8 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PRIVATE( Fe_d_ant, Fe_d_nat ) & !$OMP PRIVATE( HCHOCl, KHOCl, f_srhocl, HOCl0, L6, L6S, L6S_1 ) &!XW !$OMP PRIVATE( SRhocl, L6_1, SO4H3_vv, SO4H4_vv, fupdateHOCl_0 ) &!xw + !$OMP PRIVATE( HCHO0, HMSc, HMS0, OH0, KaqHCHO, KaqHMS, KaqHMS2 ) &!jmm + !$OMP PRIVATE( L7, L7S, L7_b, L7S_b, L8, L8S, LSTOT_HMS ) &!jmm !$OMP SCHEDULE( DYNAMIC ) DO L = 1, State_Grid%NZ DO J = 1, State_Grid%NY @@ -2566,6 +2584,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & Ld = 0.0_fp LSTOT0 = 0.0_fp LSTOT = 0.0_fp + LSTOT_HMS = 0.0_fp ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE @@ -2574,12 +2593,18 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & SO20 = Spc(I,J,L,id_SO2) H2O20 = Spc(I,J,L,id_H2O2) IF ( IS_FULLCHEM ) THEN - ! HOBr is only a defined species in fullchem simulations + ! HOBr, OH, and HCHO are only a defined species in fullchem simulations HOBr0 = Spc(I,J,L,id_HOBr) !(qjc, 01/30/17) HOCl0 = Spc(I,J,L,id_HOCl) !XW + HCHO0 = Spc(I,J,L,id_CH2O) !(jmm, 06/13/2018) + HMS0 = Spc(I,J,L,id_HMS) !(jmm, 06/15/2018) + OH0 = Spc(I,J,L,id_OH) !(jmm, 06/26/18) ELSE HOBr0 = 0.e+0_fp HOCl0 = 0.e+0_fp + HCHO0 = 0.e+0_fp !(jmm 06/13/2018) + HMS0 = 0.e+0_fp !(jmm, 06/15/2018) + OH0 = 0.e+0_fp !(jmm, 06/26/18) ENDIF ! Calculate O3, defined only in the chemistry grid @@ -2857,6 +2882,15 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L4S = 0.e+0_fp L6S = 0.e+0_fp !XW L6S_1 = 0.e+0_fp !XW + KaqHCHO = 0.e+0_fp !(jmm, 06/07/18) + KaqHMS = 0.e+0_fp !(jmm, 06/07/18) + KaqHMS2 = 0.e+0_fp !(jmm, 06/26/18) + L7 = 0.e+0_fp !(jmm, 06/13/18) + L7_b = 0.e+0_fp !(jmm, 06/13/18) + L7S = 0.e+0_fp !(jmm, 06/13/18) + L7S_b = 0.e+0_fp !(jmm, 06/13/18) + L8 = 0.e+0_fp !(jmm, 06/26/18) + L8S = 0.e+0_fp !(jmm, 06/26/18) ! If (1) there is cloud, (2) there is SO2 present, and ! (3) the T > -15 C, then compute aqueous SO2 chemistry @@ -2948,6 +2982,17 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & Spc(I,J,L,id_SO4s) * State_Met%AIRDEN(I,J,L) & / ( AIRMW * LWC) + IF ( IS_FULLCHEM ) THEN + ! Get HMS cloud concentration and convert from [v/v] to + ! [moles/liter] (jmm, 06/13/2018) + ! Use a cloud scavenging ratio of 0.7 + ! assume nonvolatile like sulfate for realistic cloud pH + HMSc = Spc(I,J,L,id_HMS) * State_Met%AIRDEN(I,J,L) * & + 0.7e+0_fp / ( AIRMW * LWC ) + ELSE + HMSc = 0e+0_fp + ENDIF + ! Get total ammonia (NH3 + NH4+) concentration [v/v] ! Use a cloud scavenging ratio of 0.7 for NH4+ TNH3 = ( Spc(I,J,L,id_NH4) * 0.7e+0_fp ) + Spc(I,J,L,id_NH3) @@ -3023,8 +3068,8 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ENDIF ! Calculate cloud pH - CALL GET_HPLUS( SO4nss, TNH3, TNO3, SO2_ss, CL, TNA, TDCA, & - TFA, TAA, TK, & + CALL GET_HPLUS( SO4nss, HMSc, TNH3, TNO3, SO2_ss, CL, TNA, & + TDCA, TFA, TAA, TK, & PATM, LWC, HPLUS_45, HPLUS ) ! Store the cloud pH @@ -3135,7 +3180,8 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & CALL AQCHEM_SO2( LWC, TK, PATM, SO2_ss, H2O20, & O3, HPLUS, MnII, FeIII, & KaqH2O2, KaqO3, KaqO3_1, KaqO2, & - HSO3aq, SO3aq ) + HSO3aq, SO3aq, HCHO0, KaqHCHO, & + KaqHMS, KaqHMS2 ) !---------------------------------------------------------- ! Compute loss by H2O2. Prevent floating-point exception @@ -3316,6 +3362,163 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & fupdateHOCl_0 = 1.e+0_fp ENDIF + !---------------------------------------------------------- + ! Compute loss by HCHO. Prevent floating-point exception + ! by not allowing the exponential to go to infinity if + ! the argument is too large. (jmm, 06/13/18) + !---------------------------------------------------------- + IF (IS_FULLCHEM ) THEN + ! Argument of the exponential + XX = ( SO2_ss - HCHO0 ) * KaqHCHO * DTCHEM + + ! Test if EXP(XX) can be computed w/o numerical exception + IF ( IS_SAFE_EXP( XX ) .and. ABS( XX ) > 0e+0_fp ) THEN + + ! Aqueous phase SO2 loss rate w/ HCHO [v/v/timestep] + L7 = EXP( XX ) + + ! Loss by HCHO + L7S = SO2_ss * HCHO0 * ( L7 - 1.e+0_fp ) / & + ( (SO2_ss * L7) - HCHO0 ) + ELSE + + ! NOTE from Jintai Lin (4/28/10): + ! However, in the case of a negative XX, L7S should be + ! approximated as SO2_ss, instead of HCHO0. In other words, + ! L7S = SO2_ss * HCHO0 * ( L7 - 1.D0 ) / ( (SO2_ss*L7) - HCHO0 ) + ! reaches different limits when XX reaches positive infinity + ! and negative infinity. + IF ( XX > 0.e+0_fp ) THEN + L7S = HCHO0 + ELSE IF ( XX < 0.e+0_fp) THEN + L7S = SO2_ss + ELSE + !(qjc, 04/10/16) different solution when SO2_ss = H2O20 + L7S = SO2_ss - 1/(KaqHCHO*DTCHEM+1/SO2_ss) + ENDIF + + ENDIF + + IF ( L7S > 1.0e+15_fp .or. L7S < 0 ) THEN + RC = GC_FAILURE + PRINT *,'Loc:',I,J,L + PRINT *,'L7S:',L7S + PRINT *,'HCHO :',HCHO0 + PRINT *,'SO2_ss',SO2_ss + PRINT *,'in exp:',XX + PRINT *,'L7:',L7 + PRINT *,'KaqHCHO',KaqHCHO + ENDIF + + + !---------------------------------------------------------- + ! Compute aqueous SO2 from HMS decomposition. Prevent floating-point + ! exception by not allowing the exponential to go to infinity if + ! the argument is too large. This reaction also produces HCHO. + ! (jmm, 06/13/18) + ! + ! This is treated as a pseudo first order reaction as it depends + ! on [HMS] and [OH-]. The [OH-] is incorporated into the rate + ! constant (KaqHMS) with HPLUS determined seperately. Updates to + ! [OH-] are incorporated in the next call of GET_HPLUS, which + ! includes [HMS] in the pH calculation + ! + ! Also note that [HMS] is in mol/L, and the conversion for the + ! rate from [mol/L/s] to [v/v/s] is done here + ! + ! Followed proceedure as done with oxidation by O2 (metal catalyzed) + !---------------------------------------------------------- + + ! Argument of the exponential + XX = -KaqHMS * DTCHEM + + ! Test if EXP(XX) can be computed w/o numerical exception + IF ( IS_SAFE_EXP( XX ) ) THEN + + ! Aqueous phase SO2 production rate w/ HMS [v/v/timestep] + L7_b = EXP( XX ) + + ! Production by HMS + L7S_b = HMSc * ( 1.e+0_fp - L7_b ) + ELSE + L7S_b = HMSc + ENDIF + + ! Convert HMS produced from [mol/L] to [v/v] (jmm, 06/15/18) + L7S_b = L7S_b *LWC * 0.08205e+0_fp * TK / PATM + + + + IF ( L7S_b > 1.0e+15_fp .or. L7S_b < 0 ) THEN + RC = GC_FAILURE + PRINT *,'Loc: ',I,J,L + PRINT *,'L7S_b: ',L7S_b + PRINT *,'HMSc: ',HMSc + PRINT *,'in exp: ',XX + PRINT *,'L7_b: ',L7_b + PRINT *,'KaqHMS: ',KaqHMS + PRINT *,'TK: ',TK + PRINT *,'LWC: ',LWC + PRINT *,'PATM :',PATM + ! CALL ERROR_STOP( 'L6s_b >1e15 or <0, point 3', LOC ) + ENDIF + + + + !---------------------------------------------------------- + ! Compute HMS loss by aqueous OH. Prevent floating-point exception + ! by not allowing the exponential to go to infinity if + ! the argument is too large. (jmm, 06/27/18) + ! + ! HMS oxidation by aqueous OH is converted to v/v/s in the KaqHMS2 + ! term. See aqchem_so2 subroutine for reaction details. + ! + ! This is a simplified version of a radical chain. + ! L8S consumes 1 HMS, 1 OH, and 1 SO2 and produces 2 SO4-- + !---------------------------------------------------------- + + ! Convert KaqHMS2 from [m^6 kg^-2 s^-1] to [v/v/s] + KaqHMS2 = KaqHMS2 * State_Met%AIRDEN(I,J,L) * & + State_Met%AIRDEN(I,J,L) + ! Argument of the exponential + XX = ( HMS0 - OH0 ) * KaqHMS2 * DTCHEM + + ! Test if EXP(XX) can be computed w/o numerical exception + IF ( IS_SAFE_EXP( XX ) .and. ABS( XX ) > 0e+0_fp ) THEN + + ! Aqueous phase HMS loss rate w/ OH(aq) [v/v/timestep] + L8 = EXP( XX ) + + ! Loss by OH(aq) + L8S = HMS0 * OH0 * ( L8 - 1.e+0_fp ) / & + ( (HMS0 * L8) - OH0 ) + ELSE + + ! NOTE from Jintai Lin (4/28/10): + ! However, in the case of a negative XX, L8S should be + ! approximated as , instead of HCHO0. In other words, + ! L8S = HMS0 * OH0 * ( L7 - 1.D0 ) / ( (HMS0*L8) - OH0 ) + ! reaches different limits when XX reaches positive infinity + ! and negative infinity. + IF ( XX > 0.e+0_fp ) THEN + L8S = OH0 + ELSE IF ( XX < 0.e+0_fp) THEN + L8S = HMS0 + ELSE + !(qjc, 04/10/16) different solution when HMS0 = OH0 + L8S = HMS0 - 1/(KaqHMS2*DTCHEM+1/HMS0) + ENDIF + + ENDIF + + ELSE + L7S = 0.e+0_fp + L7S_b = 0.e+0_fp + L8S = 0.e+0_fp + ENDIF + + + L2S = L2S * FC L3S = L3S * FC L3S_1 = L3S_1 * FC !(qjc, 06/20/16) @@ -3325,8 +3528,32 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L5S_1 = L5S_1 !(qjc, 06/20/16) L6S = L6S L6S_1 = L6S_1 + + L7S = L7S * FC ! Note: consumes SO2 but for HMS (jmm, 06/13/18) + L7S_b = L7S_b * FC ! Note: releases SO2 from HMS (jmm, 06/13/18) + L8S = L8S * FC ! Notes: releases 2 sulfate and 1 HCHO for 1 HMS and 1 SO2 (jmm, 06/29/18) + + + ! make sure HMS produced is less than HCHO available + ! (jmm, 06/28/18) + LSTOT_HMS = L7S - L7S_b - L8S + IF (LSTOT_HMS > HCHO0) THEN + L7S = HCHO0 + L7S_b + L8S + ELSE + L7S = L7S + ENDIF - LSTOT0 = L2S + L3S + L4S + L5S + L6S ! (qjc, 11/04/16) + ! make sure HMS decomposed to HCHO and SO2 is less than HMS available + ! (jmm, 06/28/18) + IF (-LSTOT_HMS > HMS0) THEN + L7S_b = ( HMS0 + L7S ) * L7S_b / ( L7S_b + L8S ) + L7S = ( HMS0 + L7S ) * L8S / ( L7S_b + L8S ) + ELSE + L7S_b = L7S_b + L8S = L8S + ENDIF + + LSTOT0 = L2S + L3S + L4S + L5S + L6S + L7S - L7S_b + L8S ! (qjc, 11/04/16), (jmm, 06/13/18) ! make sure sulfate produced is less than SO2 available ! (qjc, 06/20/16) @@ -3339,6 +3566,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L3S_1 = SO2_ss * L3S_1 / LSTOT0 L5S_1 = SO2_ss * L5S_1 / LSTOT0 L6S_1 = SO2_ss * L6S_1 / LSTOT0 + L7S = SO2_SS * L7S / LSTOT0 ! (jmm, 06/13/18) + L8S = SO2_SS * L8S / LSTOT0 ! (jmm, 06/29/18) + ! This is the ratio used to calculate the actual removal of ! HOBr by SO2 for use in gckpp_HetRates.F90 (qjc, 06/20/16) @@ -3355,6 +3585,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L3S_1 = L3S_1 L5S_1 = L5S_1 L6S_1 = L6S_1 + L7S = L7S ! (jmm, 06/13/18) + L8S = L8S ! (jmm, 06/29/18) + ! This is the ratio used to calculate the actual removal of ! HOBr by SO2 (qjc, 06/20/16) @@ -3525,7 +3758,8 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Get total in-cloud sulfate production based on bulk cloud pH ! calculations for use in HET_DROP_CHEM - LSTOT = (L2S + L3S + L4S + L5S + L6S) / FC !(qjc, 06/20/16) + ! added in SO4 from HMS (jmm, 06/29/18) + LSTOT = (L2S + L3S + L4S + L5S + L6S + L8S + L8S) / FC !(qjc, 06/20/16) ! Get coarse-mode sea-salt concentration for use in ! HET_DROP_CHEM [v/v] @@ -3539,7 +3773,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! This is to make sure HET_DROP_CHEM does not compute more ! sulfate then there is SO2 ! Add L5S (qjc, 11/04/16) - SO2_sr = MAX( SO2_ss - ( L2S + L3S + L4S + L5S + L6S ), MINDAT) + ! Add L7S, L7S_b, and L8S (jmm, 06/29/18) + SO2_sr = MAX( SO2_ss - ( L2S + L3S + L4S + L5S + L6S + & + L7S - L7S_b + L8S), MINDAT) CALL HET_DROP_CHEM( I, J, L, LSTOT, SSCvv, & aSO4, NH3, SO2_sr, H2O20, GNO3, SR, & @@ -3595,8 +3831,17 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Make sure SO2_ss and H2O20 are in the proper range ! Add L5S (qjc, 11/04/16) - SO2_ss = MAX( SO2_ss - ( L2S+L3S+L4S+L5S+L6S+SR ), MINDAT ) + ! Add L7S and L7S_b (jmm, 06/13/18) + ! Add loss and production of HCHO (jmm, 06/13/18) + ! Add loss and production of HMS (jmm, 06/15/18) + SO2_ss = MAX( SO2_ss - ( L2S+L3S+L4S+L5S+L6S+ & + L7S-L7S_b+L8S+SR ), MINDAT ) H2O20 = MAX( H2O20 - L2S, MINDAT ) + HCHO0 = MAX( HCHO0 - L7S + L7S_b + L8S, MINDAT ) + HMS0 = MAX( HMS0 + L7S - L7S_b - L8S, MINDAT ) + O3 = MAX( O3 - L3S, MINDAT ) + OH0 = MAX( OH0 - L8S, MINDAT ) + ! Factor to calculate effective SO3 and HSO3 in aqchem_so2 to ! be used in HOBr+HSO3/SO3 (qjc, 06/20/16) @@ -3613,8 +3858,13 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & H2O2s(I,J,L) = H2O20 ! SO2 chemical loss rate = SO4 production rate [v/v/timestep] + ! Add HMS (jmm, 06/13/18) PSO4_SO2(I,J,L) = LSTOT + PSO4E PSO4_ss (I,J,L) = PSO4F + PHMS_SO2(I,J,L) = L7S + PSO2_HMS(I,J,L) = L7S_b + PSO4_HMS(I,J,L) = L8S + L8S + ELSE @@ -3635,6 +3885,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & SRhocl = 0.e+0_fp SRo3 = 0.e+0_fp HPLUS = 0.e+0_fp + L7S = 0.e+0_fp !(jmm, 06/13/18) + L7S_b = 0.e+0_fp !(jmm, 06/13/18) + L8S = 0.e+0_fp !(jmm, 06/13/18) ! Store HSO3aq, SO3aq for use in gckpp_HetRates.F90 ! Avoid divide-by-zero errors @@ -3647,15 +3900,25 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & State_Chm%fupdateHOCl(I,J,L) = 0.e+0_fp ! SO2 chemical loss rate = SO4 production rate [v/v/timestep] + ! Add HMS (jmm, 06/13/18) PSO4_SO2(I,J,L) = PSO4E PSO4_ss (I,J,L) = PSO4F + PHMS_SO2(I,J,L) = 0.0e+0_fp + PSO2_HMS(I,J,L) = 0.0e+0_fp + PSO4_HMS(I,J,L) = 0.0e+0_fp + ENDIF ! Store updated SO2, H2O2 back to the tracer arrays + ! Add HCHO, OH, O3, and HMS (jmm, 06/13/18) If (FullRun) Then Spc(I,J,L,id_SO2) = SO2s( I,J,L) Spc(I,J,L,id_H2O2) = H2O2s(I,J,L) + Spc(I,J,L,id_CH2O) = HCHO0 + Spc(I,J,L,id_HMS) = HMS0 + Spc(I,J,L,id_OH) = OH0 + Spc(I,J,L,id_O3) = O3 ! Set SO4H1 and SO4H2 to zero at end of each timestep IF ( id_SO4H1 > 0 ) Spc(I,J,L,id_SO4H1) = 0.0e+0_fp @@ -3666,7 +3929,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! SO2 chemical loss rate = SO4 production rate [v/v/timestep] ! Add L5S (qjc, 11/04/16) - PSO4_SO2(I,J,L) = L1 + L2S + L3S + L4S + L5S + PSO4E + SR + L6S + ! Add L8S (jmm, 06/29/18) + PSO4_SO2(I,J,L) = L1 + L2S + L3S + L4S + L5S + PSO4E + SR + & + L6S + L8S + L8S #ifdef LUO_WETDEP ! Luo et al scheme: archive PSO4s @@ -3771,6 +4036,28 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ( L3S_1 * State_Met%AD(I,J,L) / TCVV_S ) / DTCHEM ENDIF + ! P(HMS) by SO2 and HCHO [kg S/s] + ! jmm (06/13/18) + IF ( State_Diag%Archive_ProdHMSfromSO2andHCHOinCloud ) THEN + State_Diag%ProdHMSfromSO2andHCHOinCloud(I,J,L) = & + ( L7S * State_Met%AD(I,J,L) / TCVV_S ) / DTCHEM + ENDIF + + ! P(SO2 and HCHO) from HMS decomp [kg S/s] + ! jmm (06/13/18) + IF ( State_Diag%Archive_ProdSO2andHCHOfromHMSinCloud ) THEN + State_Diag%ProdSO2andHCHOfromHMSinCloud(I,J,L) = & + ( L7S_b * State_Met%AD(I,J,L) / TCVV_S ) / DTCHEM + ENDIF + + ! P(SO4) from aqueous-phase oxidation of HMS [kg S/s] + ! jmm (07/06/18) + IF ( State_Diag%Archive_ProdSO4fromHMSinCloud ) THEN + State_Diag%ProdSO4fromHMSinCloud(I,J,L) = & + ( 2.e+0_fp * L8S * State_Met%AD(I,J,L) / & + TCVV_S ) / DTCHEM + ENDIF + !----------------------------------------------------------- ! Diagnostics for acid uptake on dust aerosol simulations !----------------------------------------------------------- @@ -4860,7 +5147,7 @@ END SUBROUTINE DUST_CHEM !\\ ! !INTERFACE: ! - SUBROUTINE GET_HPLUS( SO4nss, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & + SUBROUTINE GET_HPLUS( SO4nss, HMsc, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & TAA, T, PRES, LWC, iHPLUS, HPLUS ) ! ! !USES: @@ -4870,6 +5157,7 @@ SUBROUTINE GET_HPLUS( SO4nss, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & ! !INPUT PARAMETERS: ! REAL(fp), INTENT(IN) :: SO4nss ! Total nss sulfate mixing ratio [M] + REAL(fp), INTENT(IN) :: HMSc ! Total HMS mixing ratio [M] REAL(fp), INTENT(IN) :: TNO3 ! Total nitrate (gas+particulate) mixing ! ratio [v/v] REAL(fp), INTENT(IN) :: TNH3 ! NH3 mixing ratio [v/v] @@ -4892,7 +5180,7 @@ SUBROUTINE GET_HPLUS( SO4nss, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & ! ============================================================================ ! Solve the following electroneutrality equation: ! [H+] = 2[SO4--] + [Cl-] + [OH-] + [HCO3-] + 2[CO3--] + [HSO3-] + 2[SO3--] +! -! [NO3-] + [HCOO-] + [CH3COO-] - [Na] - 2[Ca] - [NH4] +! [NO3-] + [HCOO-] + [CH3COO-] +[HMS] - [Na] - 2[Ca] - [NH4] ! Uses Newton's method to solve the equation: ! x_1 = x_0 -f(x_0)/f'(x_0) ! iterate until converge @@ -4948,7 +5236,7 @@ SUBROUTINE GET_HPLUS( SO4nss, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & ! Non-volatile aerosol concentration [M] ! For now sulfate is the only non-volatile species - D = (2.e+0_fp*SO4nss) - TNA - (2.e+0_fp*TDCA) + D = (2.e+0_fp*SO4nss) - TNA - (2.e+0_fp*TDCA) + (1.e+0_fp * HMSc) ! Temperature dependent water equilibrium constant Kw_T = Kw*exp(DhrKw*((1.e+0_fp/T)-(1.e+0_fp/298.e+0_fp))) @@ -6400,7 +6688,8 @@ END SUBROUTINE CaCO3_PRECIP SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & O3, Hplus, MnII, FeIII, & KaqH2O2, KaqO3, KaqO3_1, KaqO2, & - HSO3aq, SO3aq ) + HSO3aq, SO3aq, HCHO, KaqHCHO, & + KaqHMS, KaqHMS2 ) ! ! !INPUT PARAMETERS: ! @@ -6413,6 +6702,8 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp), INTENT(IN) :: HPLUS ! Concentration of H+ ion (i.e. pH) [v/v] REAL(fp), INTENT(IN) :: MnII ! Concentration of MnII [mole/l] REAL(fp), INTENT(IN) :: FeIII ! Concentration of FeIII [mole/l] + REAL(fp), INTENT(IN) :: HCHO ! HCHO mixing ratio [v/v] (jmm, 06/13/18) + ! ! !OUTPUT PARAMETERS: ! @@ -6420,6 +6711,9 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp), INTENT(OUT) :: KaqO3 ! Reaction rate for O3 REAL(fp), INTENT(OUT) :: KaqO3_1 ! only the SO3-- oxidation, (qjc, 04/10/16) REAL(fp), INTENT(OUT) :: KaqO2 ! Reaction rate for O2 (metal cat) + REAL(fp), INTENT(OUT) :: KaqHCHO ! Reaction rate for SO2 and HCHO (jmm, 06/13/18) + REAL(fp), INTENT(OUT) :: KaqHMS ! Reaction rate for HMS and OH- (jmm, 06/13/18) + REAL(fp), INTENT(OUT) :: KaqHMS2 ! Reaction rate for HMS and OH(aq) (jmm, 06/28/18) REAL(fp), INTENT(OUT) :: HSO3aq ! Cloud bisulfite [mol/l] (qjc, 06/10/16) REAL(fp), INTENT(OUT) :: SO3aq ! Cloud sulfite [mol/l] (qjc, 06/10/16) ! @@ -6439,6 +6733,28 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! d[S(VI)]/dt = (k0[SO2(aq)] + k1[HSO3-] + K2[SO3--])[O3(aq)] ! [Seinfeld and Pandis, 1998, page 363] ! . +! (R3) HSO3- + HCHO(aq) => HMS +! SO3-- + HCHO(aq) => HMS + OH- +! [Moch et al., 2018; Olson and Hoffman, 1986] +! . +! d[S(HMS)]/dt = (k1[HSO3-] + k2[SO3--])[HCHO(aq)] +! [Seinfeld and Pandis, 2016, 309] +! +! (R4) HMS + OH- => HCHO(aq) + SO3-- +! [Moch et al., 2018; Deister et al., 1986] +! (note treated as 1st order in contrast to other reactions here) +! +! (R5) HMS + OH(aq) =(SO2,HO2,O2)=> HCHO + 2SO4-- + O2 + 3H+ + 2H2O +! [Jacob et al, 1986, Olson and Fessenden, 1992; +! Seinfeld and Pandis, 2016, Table 7A.7] +! Net reaction (R5): +! HMS + OH(aq) =(O2)=> SO5- + HCHO + H2O +! HO2 <=> H+ + O2- +! SO5- + O2- =(H2O)=> HSO5- + OH- + O2 +! SO2(aq) <=> HSO3- + H+ +! H+ + OH- <=> H2O +! HSO5- + HSO3- => 2SO4-- + 2H+ +! ! Reaction rates can be given as ! Ra = k [H2O2(ag)] [S(IV)] [mole/liter*s] OR ! Krate = Ra LWC R T / P [1/s] @@ -6464,7 +6780,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! (c ) Then Calculate Aquous phase, S[IV] concentrations ! S[IV] = Kheff * P(so2(g) in atm) [M] ! . -! (d ) The exact same procedure is applied to calculate H2O2(aq) +! (d ) The exact same procedure is applied to calculate H2O2(aq) and HCHO(aq) ! ! !REVISION HISTORY: ! (1 ) Updated by Rokjin Park (rjp, bmy, 12/12/02) @@ -6476,6 +6792,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! !DEFINED PARAMETERS: ! REAL(fp), PARAMETER :: R = 0.08205e+0_fp + REAL(fp), PARAMETER :: dOH = 1.0e-19_fp ! [M cm^3 molec^-1] (jmm, 06/28/18) ! ! !LOCAL VARIABLES: ! @@ -6483,7 +6800,10 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp) :: FHCSO2, XSO2G, SIV, HSO3, XSO2AQ REAL(fp) :: XHSO3, XSO3, KH1, HCH2O2, FHCH2O2 REAL(fp) :: XH2O2G, H2O2aq, KO0, KO1, KO2 - REAL(fp) :: HCO3, XO3g, O3aq + REAL(fp) :: HCO3, XO3g, O3aq, XHCHOg, HCHCHO ! (jmm, 06/13/18) + REAL(fp) :: FHCHCHO,KHCHO1, KHCHO2, KHMS ! (jmm, 06/13/18) + REAL(fp) :: KW1, KHC1, KHMS2 ! (jmm, 06/15/18) + !================================================================= ! AQCHEM_SO2 begins here! @@ -6633,7 +6953,116 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! Conversion rate from SO2 to SO4 via reaction with O2 (met cat) KaqO2 = FHCSO2 * XSO2g * ( (750e+0_fp * MnII ) + & ( 2600e+0_fp * FeIII ) + (1e+10_fp * MnII * FeIII ) ) * & - LWC * R * T ! [s-1] + LWC * R * T ! [s-1] + + !================================================================= + ! Aqueous reactions of SO2 with HCHO + ! HSO3- + HCHO(aq) => HMS + OH- (1) + ! SO3-- + HCHO(aq) => HMS (2) + ! + ! NOTE: + ! [Boyce and Hoffman, 1984] + ! KHCHO1 = 7.9E2 * EXP( -16.435 * ((298.15/T) - 1.)) + ! KHCHO2 = 2.5E7 * EXP( -6.037 * ((298.15/T) - 1.)) + ! + ! + ! Aqueous reaction of HMS with OH- + ! HMS + OH- => HCHO(aq) + SO3-- (3) + ! + ! NOTE: unclear where B (E/R) value in Seinfeld and Pandis from, + ! but close to Deister. Using Seinfeld and Pandis value for now + ! [Deister et al., 1986] + ! KHMS = 3.6E3 * EXP( -22.027 * ((298.15/T) - 1.)) + ! [Seinfeld and Pandis, 2016; Munger et al., 1986] + ! KHMS = 3.6E3 * EXP( -15.09 * ((298.15/T) - 1.)) + ! + ! + ! Aqueous reaction of HMS with OH(aq) + ! HMS + OH(aq) =(SO2,O2,HO2)=> 2SO4-- + HCHO + O2 + 3H+ + 2H2O (4) + ! + ! NOTE: O2, SO2, and HO2 particpate in the stoichiometry but not kinetics. + ! Assume steady state for sulfur radicals and the following reaction chain: + ! HMS + OH(aq) =(O2)=> SO5- + HCHO + H2O [Olsen and Fessenden, 1992] + ! HO2 <=> H+ + O2- [Jacob, 1986] + ! SO5- + O2- =(H2O)=> HSO5- + OH- + O2 [Jacob, 1986] + ! SO2(aq) <=> HSO3- + H+ + ! H+ + OH- <=> H2O + ! HSO5- + HSO3- => 2SO4-- + 2H+ [Jacob, 1986] + ! Instead of assuming Henry's law for OH, use the parameter from + ! Jacob et al, 2005 that relates gas phase OH to aqueous phase OH + ! accounting for the HO2(aq)/O2- cylcing in cloud droplets: + ! dOH = 1E-19 [M cm^3 molec^-1] + ! [Olson and Fessenden, 1992] + ! KHMS2 = 6.2E8 * EXP( -5.03 * ((298.15/T) -1.)) + ! + ! + ! (jmm, 06/28/18) + !================================================================= + KHCHO1 = 7.9e+2_fp * EXP( -16.44e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + KHCHO2 = 2.5e+7_fp * EXP( -6.04e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + KHMS = 3.6e+3_fp * EXP( -15.09e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + KHMS2 = 6.2e+8_fp * EXP( -5.03e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + + !================================================================= + ! HCHO equilibrium reaction + ! HCHO(aq) + H2O = HCH(OH)2 + ! + ! Reaction rate dependent on Temperature is given + ! H = A exp ( B (T./T - 1) ) + ! + ! For equilibrium reactions of HCHO + ! Ah1 Bh1 + ! Sienfeld and Pandis 2.53E3 13.48 [2016] + ! + ! (jmm, 06/15/18) + !================================================================= + Khc1 = 2.53e+3_fp * EXP( 13.48e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + + !================================================================= + ! H2O equilibrium reaction + ! H2O = H+ + OH- + ! + ! Reaction rate dependent on Temperature is given + ! H = A exp ( B (T./T - 1) ) + ! + ! For equilibrium reactions of HCHO + ! Ah1 Bh1 + ! Sienfeld and Pandis 1E-14 -22.51 [2016] + ! + ! (jmm, 06/15/18) + !================================================================= + Kw1 = 1e-14_fp * EXP( -22.51e+0_fp & + * ( 298.15e+0_fp / T - 1.e+0_fp ) ) + + ! Henry's constant [mol/l-atm] and Effective Henry's constant for HCHO + ! [Seinfeld and Pandis, 2016] + ! HCHCHO = 2.5 * EXP( 21.6e+0_fp * ( 298.15e+0_fp / T - 1.e+0_fp) ) + ! (jmm, -6/15/18) + HCHCHO = 2.5e+0_fp * EXP( 21.6e+0_fp & + * (298.15e+0_fp / T - 1.e+0_fp) ) + FHCHCHO = HCHCHO * (1.e+0_fp + Khc1 ) + + XHCHOg = 1.e+0_fp / ( 1.e+0_fp + ( FHCHCHO * R * T * LWC ) ) + + + ! Conversion rate from SO2 to HMS via reaction with HCHO + ! (jmm, 06/15/18) + KaqHCHO = (KHCHO1*XHSO3 + KHCHO2*XSO3) * FHCSO2 * XSO2G & + * P * HCHCHO * XHCHOg * LWC * R * T ! [v/v/s] + + ! Conversion rate from HMS to SO2 via reaction with OH- + ! (jmm, 06/15/18) + KaqHMS = KHMS * ( Kw1 / Hplus ) ! [mol/L/s] + + ! Conversion rate from HMS to SO42- & HCHO via reaction with OH(aq) + ! (jmm, 06/28/18) + KaqHMS2 = KHMS2 * dOH / AIRMW / AIRMW * 7.e-4_fp * AVO & + * T * R / P ! [m^6 kg^-2 s^-1] END SUBROUTINE AQCHEM_SO2 !EOC @@ -7019,7 +7448,7 @@ SUBROUTINE CHEM_SO4( Input_Opt, State_Chm, State_Diag, State_Grid, & ! SO4 chemistry !============================================================== - ! SO4 production from SO2 [v/v/timestep] + ! SO4 production from SO2 and HMS [v/v/timestep] SO4 = SO40 + PSO4_SO2(I,J,L) !============================================================== @@ -8293,6 +8722,24 @@ SUBROUTINE INIT_SULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) IF ( RC /= GC_SUCCESS ) RETURN PSO4_SO2 = 0e+0_fp + ALLOCATE( PHMS_SO2( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & + STAT=RC ) + CALL GC_CheckVar( 'sulfate_mod.F90:PHMS_SO2', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + PHMS_SO2 = 0e+0_fp + + ALLOCATE( PSO2_HMS( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & + STAT=RC ) + CALL GC_CheckVar( 'sulfate_mod.F90:PSO2_HMS', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + PSO2_HMS = 0e+0_fp + + ALLOCATE( PSO4_HMS( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & + STAT=RC ) + CALL GC_CheckVar( 'sulfate_mod.F90:PSO4_HMS', 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + PSO4_HMS = 0e+0_fp + ALLOCATE( PSO4_ss( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & STAT=RC ) CALL GC_CheckVar( 'sulfate_mod.F90:PSO4_ss', 0, RC ) @@ -8463,6 +8910,8 @@ SUBROUTINE INIT_SULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) id_SO4H4 = Ind_('SO4H4' ) id_HCOOH = Ind_('HCOOH' ) ! (jmm, 12/3/18) id_ACTA = Ind_('ACTA' ) ! (jmm, 12/3/18) + id_CH2O = Ind_('CH2O' ) ! (jmm, 06/15/18) + id_HMS = Ind_('HMS' ) ! (jmm, 06/15/18) ! Define flags for species drydep indices DRYNITs = Ind_('NITs', 'D') @@ -8521,6 +8970,10 @@ SUBROUTINE CLEANUP_SULFATE() IF ( ALLOCATED( PCCL ) ) DEALLOCATE( PCCL ) !xnw IF ( ALLOCATED( PSO2_DMS ) ) DEALLOCATE( PSO2_DMS ) IF ( ALLOCATED( PSO4_SO2 ) ) DEALLOCATE( PSO4_SO2 ) + IF ( ALLOCATED( PHMS_SO2 ) ) DEALLOCATE( PHMS_SO2 ) !(jmm, 06/15/18) + IF ( ALLOCATED( PSO2_HMS ) ) DEALLOCATE( PSO2_HMS ) !(jmm, 06/15/18) + IF ( ALLOCATED( PSO4_HMS ) ) DEALLOCATE( PSO4_HMS ) !(jmm, 06/26/18) + #ifdef APM IF ( ALLOCATED( PSO4_SO2APM ) ) DEALLOCATE( PSO4_SO2APM ) IF ( ALLOCATED( PSO4_SO2SEA ) ) DEALLOCATE( PSO4_SO2SEA ) diff --git a/Headers/species_database.yml b/Headers/species_database.yml index 22788c762..018d64ae8 100644 --- a/Headers/species_database.yml +++ b/Headers/species_database.yml @@ -3708,6 +3708,21 @@ SO4: WD_AerScavEff: 1.0 WD_KcScaleFac: [1.0, 0.5, 1.0] WD_RainoutEff: [1.0, 0.0, 1.0] +HMS: + DD_DvzAerSnow: 0.03 + DD_DvzMinVal: [0.01, 0.01] + DD_F0: 0.0 + DD_Hstar: 0.0 + Formula: HOCH2SO3 + FullName: Hydroxymethanesulfonate + Is_Advected: true + Is_Aerosol: true + Is_DryDep: true + Is_WetDep: true + MW_g: 111.10 + WD_AerScavEff: 1.0 + WD_KcScaleFac: [1.0, 0.5, 1.0] + WD_RainoutEff: [1.0, 0.0, 1.0] SO4D1: << : *DST1properties FullName: Sulfate on dust, Reff = 0.7 microns diff --git a/Headers/state_diag_mod.F90 b/Headers/state_diag_mod.F90 index d820981b6..0a2b8505f 100644 --- a/Headers/state_diag_mod.F90 +++ b/Headers/state_diag_mod.F90 @@ -260,6 +260,7 @@ MODULE State_Diag_Mod REAL(f4), POINTER :: AerMassPOA (:,:,: ) ! POA [ug/m3] REAL(f4), POINTER :: AerMassSAL (:,:,: ) ! Total seasalt [ug/m3] REAL(f4), POINTER :: AerMassSO4 (:,:,: ) ! Sulfate [ug/m3] + REAL(f4), POINTER :: AerMassHMS (:,:,: ) ! HMS [ug/m3] !(jmm, 06/28/19) REAL(f4), POINTER :: AerMassSOAGX (:,:,: ) ! SOAGX [ug/m3] REAL(f4), POINTER :: AerMassSOAIE (:,:,: ) ! SOAIE [ug/m3] REAL(f4), POINTER :: AerMassTSOA (:,:,: ) ! Terpene SOA [ug/m3] @@ -280,6 +281,7 @@ MODULE State_Diag_Mod LOGICAL :: Archive_AerMassPOA LOGICAL :: Archive_AerMassSAL LOGICAL :: Archive_AerMassSO4 + LOGICAL :: Archive_AerMassHMS ! jmm 06/29/18 LOGICAL :: Archive_AerMassSOAGX LOGICAL :: Archive_AerMassSOAIE LOGICAL :: Archive_AerMassTSOA @@ -365,6 +367,9 @@ MODULE State_Diag_Mod REAL(f4), POINTER :: ProdSO4fromSRHOBr (:,:,:) REAL(f4), POINTER :: ProdSO4fromO3s (:,:,:) REAL(f4), POINTER :: LossHNO3onSeaSalt (:,:,:) + REAL(f4), POINTER :: ProdHMSfromSO2andHCHOinCloud(:,:,:) ! (jmm, 06/29/18) + REAL(f4), POINTER :: ProdSO2andHCHOfromHMSinCloud(:,:,:) ! (jmm, 01/18/20) + REAL(f4), POINTER :: ProdSO4fromHMSinCloud (:,:,:) ! (jmm, 06/29/18) LOGICAL :: Archive_ProdSO2fromDMSandOH LOGICAL :: Archive_ProdSO2fromDMSandNO3 LOGICAL :: Archive_ProdSO2fromDMS @@ -382,6 +387,9 @@ MODULE State_Diag_Mod LOGICAL :: Archive_ProdSO4fromSRHOBr LOGICAL :: Archive_ProdSO4fromO3s LOGICAL :: Archive_LossHNO3onSeaSalt + LOGICAL :: Archive_ProdHMSfromSO2andHCHOinCloud ! (jmm, 06/29/18) + LOGICAL :: Archive_ProdSO2andHCHOfromHMSinCloud ! (jmm, 01/18/20) + LOGICAL :: Archive_ProdSO4fromHMSinCloud ! (jmm, 06/29/18) ! O3 and HNO3 at a given height above the surface REAL(f4), POINTER :: DryDepRaALT1 (:,: ) @@ -949,6 +957,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%AerMassPOA => NULL() State_Diag%AerMassSAL => NULL() State_Diag%AerMassSO4 => NULL() + State_Diag%AerMassHMS => NULL() ! jmm 06/28/19 State_Diag%AerMassSOAGX => NULL() State_Diag%AerMassSOAIE => NULL() State_Diag%AerMassTSOA => NULL() @@ -969,6 +978,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%Archive_AerMassPOA = .FALSE. State_Diag%Archive_AerMassSAL = .FALSE. State_Diag%Archive_AerMassSO4 = .FALSE. + State_Diag%Archive_AerMassHMS = .FALSE. ! jmm 06/28/19 State_Diag%Archive_AerMassSOAGX = .FALSE. State_Diag%Archive_AerMassSOAIE = .FALSE. State_Diag%Archive_AerMassTSOA = .FALSE. @@ -1051,6 +1061,9 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%ProdSO4fromSRHOBr => NULL() State_Diag%ProdSO4fromO3s => NULL() State_Diag%LossHNO3onSeaSalt => NULL() + State_Diag%ProdSO4fromHMSinCloud => NULL() ! (jmm, 06/29/18) + State_Diag%ProdHMSfromSO2andHCHOinCloud => NULL() ! (jmm, 06/29/18) + State_Diag%ProdSO2andHCHOfromHMSinCloud => NULL() ! (jmm, 01/18/20) State_Diag%Archive_ProdSO2fromDMSandOH = .FALSE. State_Diag%Archive_ProdSO2fromDMSandNO3 = .FALSE. State_Diag%Archive_ProdSO2fromDMS = .FALSE. @@ -1068,6 +1081,9 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%Archive_ProdSO4fromSRHOBr = .FALSE. State_Diag%Archive_ProdSO4fromO3s = .FALSE. State_Diag%Archive_LossHNO3onSeaSalt = .FALSE. + State_Diag%Archive_ProdSO4fromHMSinCloud = .FALSE. ! (jmm, 06/29/18) + State_Diag%Archive_ProdHMSfromSO2andHCHOinCloud= .FALSE. ! (jmm, 06/29/18) + State_Diag%Archive_ProdSO2andHCHOfromHMSinCloud= .FALSE. ! (jmm, 01/18/20) ! O3 and HNO3 at a given height above the surface State_Diag%DryDepRaALT1 => NULL() @@ -4016,6 +4032,66 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & IF ( RC /= GC_SUCCESS ) RETURN ENDIF + !-------------------------------------------------------------------- + ! Production of HMS from aqueous reaction of SO2 in cloud + ! (jmm, 06/29/18) + !-------------------------------------------------------------------- + arrayID = 'State_Diag%ProdHMSfromSO2andHCHOinCloud' + diagID = 'ProdHMSfromSO2andHCHOinCloud' + CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) + IF ( Found ) THEN + IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%ProdHMSfromSO2andHCHOinCloud( IM, JM, LM ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%ProdHMSfromSO2andHCHOinCloud = 0.0_f4 + State_Diag%Archive_ProdHMSfromSO2andHCHOinCloud = .TRUE. + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%ProdHMSfromSO2andHCHOinCloud, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + !-------------------------------------------------------------------- + ! Production of SO2 and HCHO from aqueous reaction of HMS in cloud + ! (jmm, 06/29/18) + !-------------------------------------------------------------------- + arrayID = 'State_Diag%ProdSO2andHCHOfromHMSinCloud' + diagID = 'ProdSO2andHCHOfromHMSinCloud' + CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) + IF ( Found ) THEN + IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%ProdSO2andHCHOfromHMSinCloud( IM, JM, LM ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%ProdSO2andHCHOfromHMSinCloud = 0.0_f4 + State_Diag%Archive_ProdSO2andHCHOfromHMSinCloud = .TRUE. + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%ProdSO2andHCHOfromHMSinCloud, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + !-------------------------------------------------------------------- + ! Production of SO4 from aqueous oxidation of HMS in cloud + ! (jmm, 06/29/18) + !-------------------------------------------------------------------- + arrayID = 'State_Diag%ProdSO4fromHMSinCloud' + diagID = 'ProdSO4fromHMSinCloud' + CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) + IF ( Found ) THEN + IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%ProdSO4fromHMSinCloud( IM, JM, LM ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%ProdSO4fromHMSinCloud = 0.0_f4 + State_Diag%Archive_ProdSO4fromHMSinCloud = .TRUE. + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%ProdSO4fromHMSinCloud, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + !-------------------------------------------------------------------- ! Loss of HNO3 on sea salt !-------------------------------------------------------------------- @@ -4130,6 +4206,26 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & IF ( RC /= GC_SUCCESS ) RETURN ENDIF + !------------------------------------------------------------------- + ! Aerosol mass of HMS [ug/m3] + ! (jmm, 06/29/18) + !------------------------------------------------------------------- + arrayID = 'State_Diag%AerMassHMS' + diagID = 'AerMassHMS' + CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) + IF ( Found ) THEN + IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) + ALLOCATE( State_Diag%AerMassHMS( IM, JM, LM ), STAT=RC ) + CALL GC_CheckVar( arrayID, 0, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + State_Diag%AerMassHMS = 0.0_f4 + State_Diag%Archive_AerMassHMS = .TRUE. + CALL Register_DiagField( Input_Opt, diagID, & + State_Diag%AerMassHMS, & + State_Chm, State_Diag, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + !------------------------------------------------------------------- ! PM2.5, aka prticulate matter with (r < 2.5 um) [ug/m3] !------------------------------------------------------------------- @@ -4326,7 +4422,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & ! being requested as diagnostic output when the corresponding ! array has not been allocated. !------------------------------------------------------------------- - DO N = 1, 21 + DO N = 1, 25 ! Select the diagnostic ID SELECT CASE( N ) @@ -4375,6 +4471,14 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & diagID = 'TotalOA' CASE( 21 ) diagID = 'TotalOC' + CASE( 22 ) ! (jmm, 06/29/18) + diagID = 'ProdSO4fromHMSinCloud' + CASE( 23 ) ! (jmm, 06/29/18) + diagID = 'ProdHMSfromSO2andHCHOinCloud' + CASE( 24 ) ! (jmm, 06/29/18) + diagID = 'AerMassHMS' + CASE( 25 ) ! (jmm, 06/29/18) + diagID = 'ProdSO2andHCHOfromHMSinCloud' END SELECT ! Exit if any of the above are in the diagnostic list @@ -6244,6 +6348,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, & State_Diag%Archive_AerMassPOA .or. & State_Diag%Archive_AerMassSAL .or. & State_Diag%Archive_AerMassSO4 .or. & + State_Diag%Archive_AerMassHMS .or. & !(jmm, 06/29/18) State_Diag%Archive_AerMassSOAGX .or. & State_Diag%Archive_AerMassSOAIE .or. & State_Diag%Archive_AerMassTSOA .or. & @@ -7097,6 +7202,27 @@ SUBROUTINE Cleanup_State_Diag( State_Diag, RC ) State_Diag%ProdSO4fromH2O2inCloud => NULL() ENDIF + ! (jmm, 06/29/18) + IF ( ASSOCIATED( State_Diag%ProdSO4fromHMSinCloud ) ) THEN + DEALLOCATE( State_Diag%ProdSO4fromHMSinCloud, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%ProdSO4fromHMSinCloud', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + ! (jmm, 06/29/18) + IF ( ASSOCIATED( State_Diag%ProdHMSfromSO2andHCHOinCloud ) ) THEN + DEALLOCATE( State_Diag%ProdHMSfromSO2andHCHOinCloud, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%ProdHMSfromSO2andHCHOinCloud', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + + ! (jmm, 06/29/18) + IF ( ASSOCIATED( State_Diag%ProdSO2andHCHOfromHMSinCloud ) ) THEN + DEALLOCATE( State_Diag%ProdSO2andHCHOfromHMSinCloud, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%ProdSO2andHCHOfromHMSinCloud', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + IF ( ASSOCIATED( State_Diag%ProdSO4fromO3inCloud ) ) THEN DEALLOCATE( State_Diag%ProdSO4fromO3inCloud, STAT=RC ) CALL GC_CheckVar( 'State_Diag%ProdSO4fromO3inCloud', 2, RC ) @@ -7262,6 +7388,13 @@ SUBROUTINE Cleanup_State_Diag( State_Diag, RC ) State_Diag%AerMassSO4 => NULL() ENDIF + ! (jmm, 06/28/18) + IF ( ASSOCIATED( State_Diag%AerMassHMS ) ) THEN + DEALLOCATE( State_Diag%AerMassHMS, STAT=RC ) + CALL GC_CheckVar( 'State_Diag%AerMassHMS', 2, RC ) + IF ( RC /= GC_SUCCESS ) RETURN + ENDIF + IF ( ASSOCIATED( State_Diag%AerMassSOAGX ) ) THEN DEALLOCATE( State_Diag%AerMassSOAGX, STAT=RC ) CALL GC_CheckVar( 'State_Diag%AerMassSOAGX', 2, RC ) @@ -8760,6 +8893,27 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isUnits ) Units = 'kg S s-1' IF ( isRank ) Rank = 3 + !(jmm, 06/29/18) + ELSE IF ( TRIM( Name_AllCaps ) == 'PRODSO4FROMHMSINCLOUD' ) THEN + IF ( isDesc ) Desc = 'Production of SO4 from aqueous ' // & + 'oxidation of HMS in clouds' + IF ( isUnits ) Units = 'kg S s-1' + IF ( isRank ) Rank = 3 + + !(jmm, 06/29/18) + ELSE IF ( TRIM( Name_AllCaps ) == 'PRODHMSFROMSO2ANDHCHOINCLOUD' ) THEN + IF ( isDesc ) Desc = 'Production of HMS from aqueous ' // & + 'reaction of SO2 and HCHO in clouds' + IF ( isUnits ) Units = 'kg S s-1' + IF ( isRank ) Rank = 3 + + !(jmm, 01/18/20) + ELSE IF ( TRIM( Name_AllCaps ) == 'PRODSO2ANDHCHOFROMHMSINCLOUD' ) THEN + IF ( isDesc ) Desc = 'Production of SO2 and HCHO from ' // & + 'aqueous reaction of HS and OH- in clouds' + IF ( isUnits ) Units = 'kg S s-1' + IF ( isRank ) Rank = 3 + ELSE IF ( TRIM( Name_AllCaps ) == 'PRODSO4FROMO3INCLOUD' ) THEN IF ( isDesc ) Desc = 'Production of SO4 from aqueous ' // & 'oxidation of O3 in clouds' @@ -8875,6 +9029,12 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isUnits ) Units = 'ug m-3' IF ( isRank ) Rank = 3 + ! (jmm, 06/29/18) + ELSE IF ( TRIM( Name_AllCaps ) == 'AERMASSHMS' ) THEN + IF ( isDesc ) Desc = 'Mass of hydroxymethanesulfonate aerosol' + IF ( isUnits ) Units = 'ug m-3' + IF ( isRank ) Rank = 3 + ELSE IF ( TRIM( Name_AllCaps ) == 'AERMASSSOAGX' ) THEN IF ( isDesc ) Desc = 'Mass of aerosol-phase glyoxal' IF ( isUnits ) Units = 'ug m-3' From 2ab1e2d34d490dc1645a9e4bde797c8f10d740c6 Mon Sep 17 00:00:00 2001 From: Jonathan Moch Date: Wed, 23 Sep 2020 23:08:34 -0400 Subject: [PATCH 02/16] Add limitiation by cloud pH and ionic strength for transition metal catalyzed SO2 oxidation (ref: Shao et al., 2019, ACP) --- GeosCore/sulfate_mod.F90 | 40 ++++++++++++++++++++++++++++++++++------ 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index 4211682f8..aa5a4374b 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -2422,6 +2422,8 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fp) :: Mn_tot, Mn_d, Fe_d REAL(fp) :: Fe_ant, Fe_nat, Fe_tot REAL(fp) :: Fe_d_ant, Fe_d_nat + REAL(fp) :: IONIC + REAL(fp) :: L6,L6S,SRhocl,L6_1,L6S_1 !XW REAL(fp) :: SO4H3_vv, SO4H4_vv !XW @@ -2570,7 +2572,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PRIVATE( HCO3, HCHOBr, KO3, KHOBr, f_srhobr, HOBr0, TMP ) & !$OMP PRIVATE( L4, L4S, KaqO2, DUST, Mn_ant, Mn_nat, Mn_tot ) & !$OMP PRIVATE( Fe_ant, Fe_nat, Fe_tot, Fe_d, Mn_d, FeIII, MnII ) & - !$OMP PRIVATE( Fe_d_ant, Fe_d_nat ) & + !$OMP PRIVATE( Fe_d_ant, Fe_d_nat, IONIC ) & !$OMP PRIVATE( HCHOCl, KHOCl, f_srhocl, HOCl0, L6, L6S, L6S_1 ) &!XW !$OMP PRIVATE( SRhocl, L6_1, SO4H3_vv, SO4H4_vv, fupdateHOCl_0 ) &!xw !$OMP PRIVATE( HCHO0, HMSc, HMS0, OH0, KaqHCHO, KaqHMS, KaqHMS2 ) &!jmm @@ -3176,9 +3178,13 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ENDIF + IONIC = 0.5e+0_fp * ( 4.0e+0_fp * SO4nss + & + TNH3*State_Met%AIRDEN(I,J,L)/(28.97e+0_fp*LWC) + & + TNO3*State_Met%AIRDEN(I,J,L)/(28.97e+0_fp*LWC) ) + ! Compute aqueous rxn rates for SO2 CALL AQCHEM_SO2( LWC, TK, PATM, SO2_ss, H2O20, & - O3, HPLUS, MnII, FeIII, & + O3, HPLUS, MnII, FeIII, IONIC, & KaqH2O2, KaqO3, KaqO3_1, KaqO2, & HSO3aq, SO3aq, HCHO0, KaqHCHO, & KaqHMS, KaqHMS2 ) @@ -6686,7 +6692,7 @@ END SUBROUTINE CaCO3_PRECIP ! !INTERFACE: ! SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & - O3, Hplus, MnII, FeIII, & + O3, Hplus, MnII, FeIII, IONIC, & KaqH2O2, KaqO3, KaqO3_1, KaqO2, & HSO3aq, SO3aq, HCHO, KaqHCHO, & KaqHMS, KaqHMS2 ) @@ -6702,6 +6708,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp), INTENT(IN) :: HPLUS ! Concentration of H+ ion (i.e. pH) [v/v] REAL(fp), INTENT(IN) :: MnII ! Concentration of MnII [mole/l] REAL(fp), INTENT(IN) :: FeIII ! Concentration of FeIII [mole/l] + REAL(fp), INTENT(IN) :: IONIC ! Ionic strength [mole/l]? REAL(fp), INTENT(IN) :: HCHO ! HCHO mixing ratio [v/v] (jmm, 06/13/18) ! @@ -6802,7 +6809,9 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp) :: XH2O2G, H2O2aq, KO0, KO1, KO2 REAL(fp) :: HCO3, XO3g, O3aq, XHCHOg, HCHCHO ! (jmm, 06/13/18) REAL(fp) :: FHCHCHO,KHCHO1, KHCHO2, KHMS ! (jmm, 06/13/18) - REAL(fp) :: KW1, KHC1, KHMS2 ! (jmm, 06/15/18) + REAL(fp) :: KW1, KHC1, KHMS2 ! (jmm, 06/15/18) + REAL(fp) :: Eff_mn, Eff_fe !jys + !================================================================= @@ -6951,9 +6960,28 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! =================================================================== ! Conversion rate from SO2 to SO4 via reaction with O2 (met cat) - KaqO2 = FHCSO2 * XSO2g * ( (750e+0_fp * MnII ) + & - ( 2600e+0_fp * FeIII ) + (1e+10_fp * MnII * FeIII ) ) * & + ! Commented out becasue using ionic strength pH modifiers version + !KaqO2 = FHCSO2 * XSO2g * ( (750e+0_fp * MnII ) + & + ! ( 2600e+0_fp * FeIII ) + (1e+10_fp * MnII * FeIII ) ) * & + ! LWC * R * T ! [s-1] + + ! Conversion rate from SO2 to SO4 via reaction with O2 (met cat) + ! added by shaojy16 10/13/2017 + ! takes into account pH and ionic strength + Eff_mn = 10.0**(-4.0*(SQRT(IONIC)/(1.0+SQRT(IONIC)))) + Eff_fe = 10.0**(-2.0*(SQRT(IONIC)/(1.0+SQRT(IONIC)))) + + IF ( Hplus > 10.0**(-4.2) ) THEN + KaqO2 = FHCSO2 * XSO2g * & + (3.72e+7_fp*Hplus**(-0.74)* & + (MnII*FeIII*Eff_fe*Eff_mn)) * & + LWC * R * T ! [s-1] + ELSE + KaqO2 = FHCSO2 * XSO2g * & + (2.51e+13_fp*Hplus**(0.67) * & + (MnII*FeIII*Eff_fe*Eff_mn)) * & LWC * R * T ! [s-1] + ENDIF !================================================================= ! Aqueous reactions of SO2 with HCHO From 0fbc2ce4e74ba143b2217720c3e9256986a53cb0 Mon Sep 17 00:00:00 2001 From: Jonathan Moch Date: Wed, 23 Sep 2020 23:24:46 -0400 Subject: [PATCH 03/16] add HMS species and diagnostics to input.geos and HISTORY.rc templates (ref: Moch et al., 2020, JGR) --- run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.benchmark | 4 ++++ run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.standard | 8 ++++++-- run/GCHPctm/input.geos.templates/input.geos.benchmark | 1 + run/GCHPctm/input.geos.templates/input.geos.standard | 1 + run/GEOS/input.geos.rc | 1 + 5 files changed, 13 insertions(+), 2 deletions(-) diff --git a/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.benchmark b/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.benchmark index a68e3cfcc..2672086cb 100644 --- a/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.benchmark +++ b/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.benchmark @@ -580,6 +580,7 @@ COLLECTIONS: 'AerosolMass', AerosolMass.mode: 'time-averaged' AerosolMass.fields: 'AerMassASOA ', 'GIGCchem', 'AerMassBC ', 'GIGCchem', + 'AerMassHMS ', 'GIGCchem', 'AerMassINDIOL ', 'GIGCchem', 'AerMassLVOCOA ', 'GIGCchem', 'AerMassNH4 ', 'GIGCchem', @@ -4574,6 +4575,9 @@ COLLECTIONS: 'AerosolMass', 'Prod_H2O2 ', 'GIGCchem', 'ProdBCPIfromBCPO ', 'GIGCchem', 'ProdOCPIfromOCPO ', 'GIGCchem', + 'ProdHMSfromSO2andHCHOinCloud ', 'GIGCchem', + 'ProdSO2andHCHOfromHMSinCloud ', 'GIGCchem', + 'ProdSO4fromHMSinCloud ', 'GIGCchem', 'ProdSO4fromH2O2inCloud ', 'GIGCchem', 'ProdSO4fromO2inCloudMetal ', 'GIGCchem', 'ProdSO4fromO3inCloud ', 'GIGCchem', diff --git a/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.standard b/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.standard index 557e5cf8e..221d2c6b7 100644 --- a/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.standard +++ b/run/GCHPctm/HISTORY.rc.templates/HISTORY.rc.standard @@ -316,6 +316,7 @@ COLLECTIONS: #'AerosolMass', AerosolMass.duration: 010000 AerosolMass.mode: 'time-averaged' AerosolMass.fields: 'AerMassBC ', 'GIGCchem', + 'AerMassHMS ', 'GIGCchem', 'AerMassINDIOL ', 'GIGCchem', 'AerMassLVOCOA ', 'GIGCchem', 'AerMassNH4 ', 'GIGCchem', @@ -1229,11 +1230,14 @@ COLLECTIONS: #'AerosolMass', 'Prod_H2O2 ', 'GIGCchem', 'ProdBCPIfromBCPO ', 'GIGCchem', 'ProdOCPIfromOCPO ', 'GIGCchem', - 'ProdSO4fromH2O2inCloud ', 'GIGCchem', + 'ProdHMSfromSO2andHCHOinCloud ', 'GIGCchem', + 'ProdSO2andHCHOfromHMSinCloud ', 'GIGCchem', + 'ProdSO4fromHMSinCloud ', 'GIGCchem', + 'ProdSO4fromH2O2inCloud ', 'GIGCchem', 'ProdSO4fromO2inCloudMetal ', 'GIGCchem', 'ProdSO4fromO3inCloud ', 'GIGCchem', 'ProdSO4fromO3inSeaSalt ', 'GIGCchem', - 'ProdSO4fromHOBrInCloud ', 'GIGCchem', + 'ProdSO4fromHOBrInCloud ', 'GIGCchem', 'ProdSO4fromSRO3 ', 'GIGCchem', 'ProdSO4fromSRHOBr ', 'GIGCchem', 'ProdSO4fromO3s ', 'GIGCchem', diff --git a/run/GCHPctm/input.geos.templates/input.geos.benchmark b/run/GCHPctm/input.geos.templates/input.geos.benchmark index d685cbea3..fd72a1bed 100644 --- a/run/GCHPctm/input.geos.templates/input.geos.benchmark +++ b/run/GCHPctm/input.geos.templates/input.geos.benchmark @@ -96,6 +96,7 @@ Species name : HCOOH Species name : HI Species name : HMHP Species name : HMML +Species name : HMS Species name : HNO2 Species name : HNO3 Species name : HNO4 diff --git a/run/GCHPctm/input.geos.templates/input.geos.standard b/run/GCHPctm/input.geos.templates/input.geos.standard index ff4237b49..2fa1e9820 100644 --- a/run/GCHPctm/input.geos.templates/input.geos.standard +++ b/run/GCHPctm/input.geos.templates/input.geos.standard @@ -89,6 +89,7 @@ Species name : HCOOH Species name : HI Species name : HMHP Species name : HMML +Species name : HMS Species name : HNO2 Species name : HNO3 Species name : HNO4 diff --git a/run/GEOS/input.geos.rc b/run/GEOS/input.geos.rc index 3eafb2c12..81256e478 100644 --- a/run/GEOS/input.geos.rc +++ b/run/GEOS/input.geos.rc @@ -96,6 +96,7 @@ Species name : HCOOH Species name : HI Species name : HMHP Species name : HMML +Species name : HMS Species name : HNO2 Species name : HNO3 Species name : HNO4 From e7d00d045a2e9457cd7d5987e8136d8bf3d5ec58 Mon Sep 17 00:00:00 2001 From: Melissa Sulprizio Date: Mon, 20 Sep 2021 12:53:33 -0400 Subject: [PATCH 04/16] Add bug fixes for creating integration test run directories Several of the GEOS-FP run directories were using MERRA-2 met fields by mistake. This was due to a mix-up of options in intTestCreate.sh following the removal of the trop-only chemistry option. Signed-off-by: Melissa Sulprizio --- test/GCClassic/intTestCreate.sh | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/test/GCClassic/intTestCreate.sh b/test/GCClassic/intTestCreate.sh index 72adbc52e..0780603cb 100755 --- a/test/GCClassic/intTestCreate.sh +++ b/test/GCClassic/intTestCreate.sh @@ -171,10 +171,10 @@ create_rundir "1\n7\n1\n2\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} #create_rundir "1\n1\n2\n1\n2\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_2x25_fullchem_complexSOA_merra2" -create_rundir "1\n3\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n3\n1\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_2x25_fullchem_complexSOA_SVPOA_merra2" -create_rundir "1\n3\n2\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n3\n2\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_2x25_fullchem_marinePOA_merra2" create_rundir "1\n4\n1\n2\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} @@ -343,31 +343,31 @@ dir="gc_4x5_aerosol_geosfp" create_rundir "2\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_geosfp" -create_rundir "1\n1\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_LuoWd_geosfp" -create_rundir "1\n1\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_aciduptake_geosfp" -create_rundir "1\n5\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n5\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_APM_geosfp" -create_rundir "1\n7\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n7\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_benchmark_geosfp" -create_rundir "1\n2\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n2\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_complexSOA_geosfp" -create_rundir "1\n3\n1\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n3\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_complexSOA_SVPOA_geosfp" -create_rundir "1\n3\n2\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n3\n2\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_marinePOA_geosfp" -create_rundir "1\n4\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n4\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_RRTMG_geosfp" -create_rundir "1\n8\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n8\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_TOMAS15_geosfp" create_rundir "1\n6\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} @@ -379,7 +379,7 @@ dir="gc_4x5_Hg_geosfp" create_rundir "5\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_POPs_BaP_geosfp" -create_rundir "6\n2\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "6\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_tagCH4_geosfp" create_rundir "7\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} From 67eb365c9372a562a39762928bba27cad454e186 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 20 Sep 2021 13:09:31 -0400 Subject: [PATCH 05/16] Update criterion to prevent div-by-zero in rate-law function kIIR1LTd In routine KPP/fullchem/rateLawUtilFuncs.F90, we removed the error trap on concEduct, as it was thought that any errors would be trapped by the Is_SafeDiv function. However, that was not the case, and Is_SafeDiv allowed excessively large reaction rates (which depleted the concentration). We now have set the error trap criterion to be 1e-20 molec/cm3. We should never have a concentration that low; anything lower is a denormal value. Signed-off-by: Bob Yantosca --- KPP/fullchem/rateLawUtilFuncs.F90 | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/KPP/fullchem/rateLawUtilFuncs.F90 b/KPP/fullchem/rateLawUtilFuncs.F90 index b5591cc07..0eb31b4b5 100644 --- a/KPP/fullchem/rateLawUtilFuncs.F90 +++ b/KPP/fullchem/rateLawUtilFuncs.F90 @@ -111,14 +111,11 @@ FUNCTION kIIR1Ltd( concGas, concEduct, kISource ) RESULT( kII ) kIEduct = 0.0_dp kII = 0.0_dp ! - ! Prevent div by zero - !--------------------------------------------------------------------- - ! NOTE: This line was presumably to prevent div-by-zero, but is - ! not really needed, as we now use SafeDiv to trap these errors. - ! This can prevent reaction rates from being computed if the - ! concEduct is very low. Comment out for now. (bmy, 9/15/21) - !IF ( concEduct < 100.0_dp ) RETURN - !--------------------------------------------------------------------- + ! Prevent div by zero. NOTE: Now use 1.0e-20 as the error trap for + ! concEduct. 100 (the previous criterion) was too large. We should + ! never have have a concentration as low as 1e-20 molec/cm3, as that + ! will definitely blow up the division. (bmy, 9/20/21) + IF ( concEduct < 1.0e-20_dp ) RETURN IF ( .not. Is_SafeDiv( concGas*kISource, concEduct ) ) RETURN ! ! Compute rates From ed6a82e0c3d7bad533b732cedc80805706f582a7 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 20 Sep 2021 14:50:14 -0400 Subject: [PATCH 06/16] Minor bug fixes for HMS chemistry GeosCore/sulfate_mod.F90: - Add missing comma in !$OMP PRIVATE statement in CHEM_SO2 Headers/state_diag_mod.F90 - Add HMS-related fields to TYPE(DgnState) - Add Init_and_Register calls for HMS-related fields of TYPE(DgnState) - Add Finalize calls for HMS-related fields of TYPE(DgnState) - Trimmed trailing whitespace run/GCClassic/input.geos.templates/input.geos.fullchem - Added HMS to list of advected species run/shared/species_database.yml - Added Background_VV = 1.0e-15 for HMS. This should be a better initial concentration than 1e-20 or 1e-30, and should not cause chemistry to exit with error. Signed-off-by: Bob Yantosca --- GeosCore/sulfate_mod.F90 | 2 +- Headers/state_diag_mod.F90 | 191 ++++++++++-------- .../input.geos.templates/input.geos.fullchem | 1 + run/shared/species_database.yml | 1 + 4 files changed, 106 insertions(+), 89 deletions(-) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index 1d21c4b50..69f8305a7 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -2636,7 +2636,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !$OMP PRIVATE( L6S, fupdateHOBr_0, SRhocl, L6_1, SO4H3_vv )& !$OMP PRIVATE( SO4H4_vv, fupdateHOCl_0, KaqO2, TNA, one_m_KRATE)& !$OMP PRIVATE( HCHO0, HMSc, HMS0, OH0, KaqHCHO )& - !$OMP PRIVATE( KaqHMS, KaqHMS2 L7, L7S, L7_b )& + !$OMP PRIVATE( KaqHMS, KaqHMS2, L7, L7S, L7_b )& !$OMP PRIVATE( L7S_b, L8, L8S, LSTOT_HMS )& !$OMP SCHEDULE( DYNAMIC, 1 ) DO L = 1, State_Grid%NZ diff --git a/Headers/state_diag_mod.F90 b/Headers/state_diag_mod.F90 index 318510b65..78682bd91 100644 --- a/Headers/state_diag_mod.F90 +++ b/Headers/state_diag_mod.F90 @@ -381,6 +381,9 @@ MODULE State_Diag_Mod REAL(f4), POINTER :: AerMassBC(:,:,:) LOGICAL :: Archive_AerMassBC + REAL(f4), POINTER :: AerMassHMS(:,:,:) + LOGICAL :: Archive_AerMassHMS + REAL(f4), POINTER :: AerMassINDIOL(:,:,:) LOGICAL :: Archive_AerMassINDIOL @@ -492,61 +495,61 @@ MODULE State_Diag_Mod !%%%%% Sulfur aerosols prod & loss %%%%% REAL(f4), POINTER :: ProdSO2fromDMSandOH(:,:,:) LOGICAL :: Archive_ProdSO2fromDMSandOH - + REAL(f4), POINTER :: ProdSO2fromDMSandNO3(:,:,:) LOGICAL :: Archive_ProdSO2fromDMSandNO3 - + REAL(f4), POINTER :: ProdSO2fromDMS (:,:,:) LOGICAL :: Archive_ProdSO2fromDMS - + REAL(f4), POINTER :: ProdMSAfromDMS(:,:,:) LOGICAL :: Archive_ProdMSAfromDMS - + REAL(f4), POINTER :: ProdNITfromHNO3uptakeOnDust(:,:,:) LOGICAL :: Archive_ProdNITfromHNO3uptakeOnDust - + REAL(f4), POINTER :: ProdSO4fromGasPhase(:,:,:) LOGICAL :: Archive_ProdSO4fromGasPhase - + REAL(f4), POINTER :: ProdSO4fromH2O2inCloud(:,:,:) LOGICAL :: Archive_ProdSO4fromH2O2inCloud - + REAL(f4), POINTER :: ProdSO4fromO3inCloud(:,:,:) LOGICAL :: Archive_ProdSO4fromO3inCloud - + REAL(f4), POINTER :: ProdSO4fromO2inCloudMetal(:,:,:) LOGICAL :: Archive_ProdSO4fromO2inCloudMetal - + REAL(f4), POINTER :: ProdSO4fromO3inSeaSalt(:,:,:) LOGICAL :: Archive_ProdSO4fromO3inSeaSalt - + REAL(f4), POINTER :: ProdSO4fromOxidationOnDust(:,:,:) LOGICAL :: Archive_ProdSO4fromOxidationOnDust - + REAL(f4), POINTER :: ProdSO4fromUptakeOfH2SO4g(:,:,:) LOGICAL :: Archive_ProdSO4fromUptakeOfH2SO4g - + REAL(f4), POINTER :: ProdSO4fromHOBrInCloud(:,:,:) LOGICAL :: Archive_ProdSO4fromHOBrInCloud - + REAL(f4), POINTER :: ProdSO4fromSRO3(:,:,:) LOGICAL :: Archive_ProdSO4fromSRO3 - + REAL(f4), POINTER :: ProdSO4fromSRHOBr(:,:,:) LOGICAL :: Archive_ProdSO4fromSRHOBr - + REAL(f4), POINTER :: ProdSO4fromO3s(:,:,:) LOGICAL :: Archive_ProdSO4fromO3s - + REAL(f4), POINTER :: LossHNO3onSeaSalt(:,:,:) LOGICAL :: Archive_LossHNO3onSeaSalt - + REAL(f4), POINTER :: ProdHMSfromSO2andHCHOinCloud(:,:,:) LOGICAL :: Archive_ProdHMSfromSO2andHCHOinCloud - + REAL(f4), POINTER :: ProdSO2andHCHOfromHMSinCloud(:,:,:) LOGICAL :: Archive_ProdSO2andHCHOfromHMSinCloud - + REAL(f4), POINTER :: ProdSO4fromHMSinCloud(:,:,:) LOGICAL :: Archive_ProdSO4fromHMSinCloud @@ -1409,7 +1412,7 @@ SUBROUTINE Zero_State_Diag( State_Diag, RC ) State_Diag%Archive_AerMassBC = .FALSE. State_Diag%AerMassHMS => NULL() - State_Diag_Archive_AerMassHMS = .FALSE. + State_Diag%Archive_AerMassHMS = .FALSE. State_Diag%AerMassINDIOL => NULL() State_Diag%Archive_AerMassINDIOL = .FALSE. @@ -2970,7 +2973,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! O3_MASS + ! O3_MASS !----------------------------------------------------------------------- diagID = 'O3_MASS' CALL Init_and_Register( & @@ -2992,7 +2995,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! GCCTO3 + ! GCCTO3 !----------------------------------------------------------------------- diagID = 'GCCTO3' CALL Init_and_Register( & @@ -3014,7 +3017,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! GCCTTO3 + ! GCCTTO3 !----------------------------------------------------------------------- diagID = 'GCCTTO3' CALL Init_and_Register( & @@ -3036,7 +3039,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! CHEMTOP + ! CHEMTOP !----------------------------------------------------------------------- diagID = 'CHEMTOP' CALL Init_and_Register( & @@ -3058,7 +3061,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! CHEMTROPP + ! CHEMTROPP !----------------------------------------------------------------------- diagID = 'CHEMTROPP' CALL Init_and_Register( & @@ -3080,7 +3083,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ENDIF !----------------------------------------------------------------------- - ! CONVCLDTOP + ! CONVCLDTOP !----------------------------------------------------------------------- diagID = 'CONVCLDTOP' CALL Init_and_Register( & @@ -6028,62 +6031,71 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ! Production of HMS from aqueous reaction of SO2 in cloud ! (jmm, 06/29/18) !-------------------------------------------------------------------- - arrayID = 'State_Diag%ProdHMSfromSO2andHCHOinCloud' diagID = 'ProdHMSfromSO2andHCHOinCloud' - CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) - IF ( Found ) THEN - IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) - ALLOCATE( State_Diag%ProdHMSfromSO2andHCHOinCloud( IM, JM, LM ), STAT=RC ) - CALL GC_CheckVar( arrayID, 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Diag%ProdHMSfromSO2andHCHOinCloud = 0.0_f4 - State_Diag%Archive_ProdHMSfromSO2andHCHOinCloud = .TRUE. - CALL Register_DiagField( Input_Opt, diagID, & - State_Diag%ProdHMSfromSO2andHCHOinCloud, & - State_Chm, State_Diag, RC ) - IF ( RC /= GC_SUCCESS ) RETURN + CALL Init_and_Register( & + Input_Opt = Input_Opt, & + State_Chm = State_Chm, & + State_Diag = State_Diag, & + State_Grid = State_Grid, & + DiagList = Diag_List, & + TaggedDiagList = TaggedDiag_List, & + Ptr2Data = State_Diag%ProdHMSfromSO2andHCHOinCloud, & + archiveData = State_Diag%Archive_ProdHMSfromSO2andHCHOinCloud,& + diagId = diagId, & + RC = RC ) + + IF( RC /= GC_SUCCESS ) THEN + errMsg = TRIM( errMsg_ir ) // TRIM( diagId ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN ENDIF !-------------------------------------------------------------------- ! Production of SO2 and HCHO from aqueous reaction of HMS in cloud ! (jmm, 06/29/18) !-------------------------------------------------------------------- - arrayID = 'State_Diag%ProdSO2andHCHOfromHMSinCloud' diagID = 'ProdSO2andHCHOfromHMSinCloud' - CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) - IF ( Found ) THEN - IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) - ALLOCATE( State_Diag%ProdSO2andHCHOfromHMSinCloud( IM, JM, LM ), STAT=RC ) - CALL GC_CheckVar( arrayID, 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Diag%ProdSO2andHCHOfromHMSinCloud = 0.0_f4 - State_Diag%Archive_ProdSO2andHCHOfromHMSinCloud = .TRUE. - CALL Register_DiagField( Input_Opt, diagID, & - State_Diag%ProdSO2andHCHOfromHMSinCloud, & - State_Chm, State_Diag, RC ) - IF ( RC /= GC_SUCCESS ) RETURN + CALL Init_and_Register( & + Input_Opt = Input_Opt, & + State_Chm = State_Chm, & + State_Diag = State_Diag, & + State_Grid = State_Grid, & + DiagList = Diag_List, & + TaggedDiagList = TaggedDiag_List, & + Ptr2Data = State_Diag%ProdSO2andHCHOfromHMSinCloud, & + archiveData = State_Diag%Archive_ProdSO2andHCHOfromHMSinCloud,& + diagId = diagId, & + RC = RC ) + + IF( RC /= GC_SUCCESS ) THEN + errMsg = TRIM( errMsg_ir ) // TRIM( diagId ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN ENDIF !-------------------------------------------------------------------- ! Production of SO4 from aqueous oxidation of HMS in cloud ! (jmm, 06/29/18) !-------------------------------------------------------------------- - arrayID = 'State_Diag%ProdSO4fromHMSinCloud' diagID = 'ProdSO4fromHMSinCloud' - CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) - IF ( Found ) THEN - IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) - ALLOCATE( State_Diag%ProdSO4fromHMSinCloud( IM, JM, LM ), STAT=RC ) - CALL GC_CheckVar( arrayID, 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Diag%ProdSO4fromHMSinCloud = 0.0_f4 - State_Diag%Archive_ProdSO4fromHMSinCloud = .TRUE. - CALL Register_DiagField( Input_Opt, diagID, & - State_Diag%ProdSO4fromHMSinCloud, & - State_Chm, State_Diag, RC ) - IF ( RC /= GC_SUCCESS ) RETURN + CALL Init_and_Register( & + Input_Opt = Input_Opt, & + State_Chm = State_Chm, & + State_Diag = State_Diag, & + State_Grid = State_Grid, & + DiagList = Diag_List, & + TaggedDiagList = TaggedDiag_List, & + Ptr2Data = State_Diag%ProdSO4fromHMSinCloud, & + archiveData = State_Diag%Archive_ProdSO4fromHMSinCloud, & + diagId = diagId, & + RC = RC ) + + IF( RC /= GC_SUCCESS ) THEN + errMsg = TRIM( errMsg_ir ) // TRIM( diagId ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN ENDIF - + !-------------------------------------------------------------------- ! Loss of HNO3 on sea salt !-------------------------------------------------------------------- @@ -6220,22 +6232,25 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & ! Aerosol mass of HMS [ug/m3] ! (jmm, 06/29/18) !------------------------------------------------------------------- - arrayID = 'State_Diag%AerMassHMS' diagID = 'AerMassHMS' - CALL Check_DiagList( am_I_Root, Diag_List, diagID, Found, RC ) - IF ( Found ) THEN - IF ( am_I_Root ) WRITE( 6, 20 ) ADJUSTL( arrayID ), TRIM( diagID ) - ALLOCATE( State_Diag%AerMassHMS( IM, JM, LM ), STAT=RC ) - CALL GC_CheckVar( arrayID, 0, RC ) - IF ( RC /= GC_SUCCESS ) RETURN - State_Diag%AerMassHMS = 0.0_f4 - State_Diag%Archive_AerMassHMS = .TRUE. - CALL Register_DiagField( Input_Opt, diagID, & - State_Diag%AerMassHMS, & - State_Chm, State_Diag, RC ) - IF ( RC /= GC_SUCCESS ) RETURN + CALL Init_and_Register( & + Input_Opt = Input_Opt, & + State_Chm = State_Chm, & + State_Diag = State_Diag, & + State_Grid = State_Grid, & + DiagList = Diag_List, & + TaggedDiagList = TaggedDiag_List, & + Ptr2Data = State_Diag%AerMassHMS, & + archiveData = State_Diag%Archive_AerMassHMS, & + diagId = diagId, & + RC = RC ) + + IF ( RC /= GC_SUCCESS ) THEN + errMsg = TRIM( errMsg_ir ) // TRIM( diagId ) + CALL GC_Error( errMsg, RC, thisLoc ) + RETURN ENDIF - + !------------------------------------------------------------------- ! PM2.5, aka prticulate matter with (r < 2.5 um) [ug/m3] !------------------------------------------------------------------- @@ -6426,7 +6441,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & mapData = State_Diag%Map_TotCol, & diagId = diagId, & diagFlag = 'S', & - RC = RC ) + RC = RC ) IF ( RC /= GC_SUCCESS ) THEN errMsg = TRIM( errMsg_ir ) // TRIM( diagId ) CALL GC_Error( errMsg, RC, thisLoc ) @@ -6564,7 +6579,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & CASE( 24 ) ! (jmm, 06/29/18) diagID = 'AerMassHMS' CASE( 25 ) ! (jmm, 06/29/18) - diagID = 'ProdSO2andHCHOfromHMSinCloud' + diagID = 'ProdSO2andHCHOfromHMSinCloud' END SELECT ! Exit if any of the above are in the diagnostic list @@ -8661,7 +8676,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm, State_Grid, & State_Diag%Archive_AerMassPOA .or. & State_Diag%Archive_AerMassSAL .or. & State_Diag%Archive_AerMassSO4 .or. & - State_Diag%Archive_AerMassHMS .or. & !(jmm, 06/29/18) + State_Diag%Archive_AerMassHMS .or. & !(jmm, 06/29/18) State_Diag%Archive_AerMassSOAGX .or. & State_Diag%Archive_AerMassSOAIE .or. & State_Diag%Archive_AerMassTSOA .or. & @@ -10457,7 +10472,7 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & ELSE IF ( TRIM( Name_AllCaps ) == 'NOY' ) THEN IF ( isDesc ) Desc = & - 'Reactive_nitrogen_=_NO_NO2_HNO3_HNO4_HONO_2xN2O5_PAN_OrganicNitrates_AerosolNitrates' + 'Reactive_nitrogen_=_NO_NO2_HNO3_HNO4_HONO_2xN2O5_PAN_OrganicNitrates_AerosolNitrates' IF ( isUnits ) Units = 'mol mol-1' IF ( isRank ) Rank = 3 @@ -10474,7 +10489,7 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isRank ) Rank = 3 ELSE IF ( TRIM( Name_AllCaps ) == 'O3_MASS' ) THEN - IF ( isDesc ) Desc = 'O3_grid_cell_mass_per_area' + IF ( isDesc ) Desc = 'O3_grid_cell_mass_per_area' IF ( isUnits ) Units = 'kg m-2' IF ( isRank ) Rank = 3 @@ -11802,7 +11817,7 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & IF ( isDesc ) Desc = 'Production of SO2 and HCHO from ' // & 'aqueous reaction of HS and OH- in clouds' IF ( isUnits ) Units = 'kg S s-1' - IF ( isRank ) Rank = 3 + IF ( isRank ) Rank = 3 ELSE IF ( TRIM( Name_AllCaps ) == 'PRODSO4FROMO3INCLOUD' ) THEN IF ( isDesc ) Desc = 'Production of SO4 from aqueous ' // & @@ -11813,7 +11828,7 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root, metadataID, Found, & ELSE IF ( TRIM( Name_AllCaps ) == 'AERMASSHMS' ) THEN IF ( isDesc ) Desc = 'Mass of hydroxymethanesulfonate aerosol' IF ( isUnits ) Units = 'ug m-3' - IF ( isRank ) Rank = 3 + IF ( isRank ) Rank = 3 ELSE IF ( TRIM( Name_AllCaps ) == 'AERMASSSOAGX' ) THEN IF ( isDesc ) Desc = 'Mass of aerosol-phase glyoxal' @@ -14124,7 +14139,7 @@ SUBROUTINE Get_Mapping( Input_Opt, State_Chm, TaggedDiagList, & mapData%slot2id(TagItem%index) = index ELSE - + ! Otherwise, this is a defined species. ! Call Ind_() to get the proper index mapData%slot2id(TagItem%index) = Ind_( TagItem%name, indFlag ) diff --git a/run/GCClassic/input.geos.templates/input.geos.fullchem b/run/GCClassic/input.geos.templates/input.geos.fullchem index be168a396..effecb20c 100644 --- a/run/GCClassic/input.geos.templates/input.geos.fullchem +++ b/run/GCClassic/input.geos.templates/input.geos.fullchem @@ -100,6 +100,7 @@ Species name : HCOOH Species name : HI Species name : HMHP Species name : HMML +Species name : HMS Species name : HNO2 Species name : HNO3 Species name : HNO4 diff --git a/run/shared/species_database.yml b/run/shared/species_database.yml index 7099ffc7f..2d221086f 100644 --- a/run/shared/species_database.yml +++ b/run/shared/species_database.yml @@ -1523,6 +1523,7 @@ HMML: MW_g: 102.10 WD_RetFactor: 2.0e-2 HMS: + Background_VV: 1.0e-15 DD_DvzAerSnow: 0.03 DD_DvzMinVal: [0.01, 0.01] DD_F0: 0.0 From 776e1c1ec32ba0ed32e3e0269a4c7f4d44472d98 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 20 Sep 2021 16:10:52 -0400 Subject: [PATCH 07/16] Add README.md for KPP/aciduptake folder The README.md file says that KPP/aciduptake will be used in future GC development, but is being ignored for now. Signed-off-by: Bob Yantosca --- KPP/aciduptake/README.md | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 KPP/aciduptake/README.md diff --git a/KPP/aciduptake/README.md b/KPP/aciduptake/README.md new file mode 100644 index 000000000..c316689db --- /dev/null +++ b/KPP/aciduptake/README.md @@ -0,0 +1,6 @@ +# README -- KPP/aciduptake folder + +NOTE: This folder contains files that will be used in future GEOS-Chem +development. These will be ignored for now. + +-- Bob Yantosca (20 Sep 2021) From 1bd56a3f153a843444510f8b0c97f4584733bd7a Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Mon, 20 Sep 2021 17:55:51 -0400 Subject: [PATCH 08/16] Update memory & time limits in test/GCClassic/intTestExecute_slurm.sh Change requested time from 16h to 5h. Most integration tests finish in about 3 to 3.5 hours. Change requested memory from 60000 to 75000. This is needed for the nested-grid NA simulation at 0.25 degree resolution. Signed-off-by: Bob Yantosca --- test/GCClassic/intTestExecute_slurm.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/GCClassic/intTestExecute_slurm.sh b/test/GCClassic/intTestExecute_slurm.sh index 7542424b9..1d663144e 100755 --- a/test/GCClassic/intTestExecute_slurm.sh +++ b/test/GCClassic/intTestExecute_slurm.sh @@ -2,9 +2,9 @@ #SBATCH -c 24 #SBATCH -N 1 -#SBATCH -t 0-16:00 +#SBATCH -t 0-05:00 #SBATCH -p huce_cascade -#SBATCH --mem=60000 +#SBATCH --mem=80000 #SBATCH --mail-type=END #------------------------------------------------------------------------------ From 885a198c4ca79d7cb80b7cfd44e3eccd922e6108 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 21 Sep 2021 13:32:57 -0400 Subject: [PATCH 09/16] Add formula for HMS to species_database.yml; point to new restart file run/GCClassic/createRunDir.sh - For fullchem and aerosol simulations, point to the restart files in ExtData/GEOSCHEM_RESTARTS/v2021-09. This has the same concentrations as the restart files in GC_13.0.0, plus HMS, SO3mm, and HSO3m. (will add C2H2 and C2H4 later) run/shared/species_database.yml - Add chemical formula for HMS Signed-off-by: Bob Yantosca --- run/GCClassic/createRunDir.sh | 12 +++++++----- run/shared/species_database.yml | 1 + 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/run/GCClassic/createRunDir.sh b/run/GCClassic/createRunDir.sh index 0cbe2528c..ca29075cc 100755 --- a/run/GCClassic/createRunDir.sh +++ b/run/GCClassic/createRunDir.sh @@ -1071,11 +1071,13 @@ if [[ ${met_name} = "MERRA2" ]] || [[ ${met_name} = "GEOSFP" ]]; then sample_rst=${rst_root}/v2020-02/GEOSChem.Restart.TOMAS40.${startdate}_0000z.nc4 else - # NOTE: When we add HSO3- and SO3-- as species (probably in - # 13.4.0), we will need to use the restart file in v2021-09. - # Leave this commented out for the time being. (bmy, 9/17/21). - #sample_rst=${rst_root}/v2021-09/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 - sample_rst=${rst_root}/GC_13.0.0/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 + # NOTE: We need to read the restart file in v2021-09 which + # contains the same species concentrations as in the GC_13.0.0 + # folder, with HSO3-, SO3--, HMS, C2H2, and C2H4 appended to + # it. Because HEMCO_Config.rc uses flag EFYO, the run will + # halt if not all species are found. (bmy, 9/21/21) + #sample_rst=${rst_root}/GC_13.0.0/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 + sample_rst=${rst_root}/v2021-09/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 fi elif [[ "x${sim_name}" == "xTransportTracers" ]]; then diff --git a/run/shared/species_database.yml b/run/shared/species_database.yml index 2d221086f..21b1ec05d 100644 --- a/run/shared/species_database.yml +++ b/run/shared/species_database.yml @@ -1528,6 +1528,7 @@ HMS: DD_DvzMinVal: [0.01, 0.01] DD_F0: 0.0 DD_Hstar: 0.0 + Formula: HOCH2SO3− FullName: Hydroxymethanesulfonate Is_Advected: true Is_Aerosol: true From 2bca7ea4d9a0d04ac7ae71eb8c88c3e1ff79442e Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Tue, 21 Sep 2021 17:49:30 -0400 Subject: [PATCH 10/16] Add HMS to the input.geos.aerosol template file for aerosol-only sims The HMS chemistry modifications in sulfate_mod.F90 and aerosol_mod.F90 assume HMS is an advected species. We have now added HMS to the run/GCClassic/input.geos.templates/input.geos.aerosol template file that is used in run directory creation. This should allow integration tests to pass. Signed-off-by: Bob Yantosca --- run/GCClassic/input.geos.templates/input.geos.aerosol | 1 + 1 file changed, 1 insertion(+) diff --git a/run/GCClassic/input.geos.templates/input.geos.aerosol b/run/GCClassic/input.geos.templates/input.geos.aerosol index 93a22dc15..ef2faf50c 100644 --- a/run/GCClassic/input.geos.templates/input.geos.aerosol +++ b/run/GCClassic/input.geos.templates/input.geos.aerosol @@ -33,6 +33,7 @@ Species name : DST2 Species name : DST3 Species name : DST4 Species name : H2O2 +Species name : HMS Species name : MSA Species name : NH3 Species name : NH4 From e0935cceb345764bdd4191b933fc4aea4898567a Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 11:54:17 -0400 Subject: [PATCH 11/16] In createRunDir.sh, TOMAS and fullchem restarts get copied from v2021-09/ We have added extra species (such as HMS) for new chemistry in 13.3.0 and later versions into the restart files in the ExtData/GEOSCHEM_RESTARTS/v2021-09 folder. When creating GEOS-Chem "Classic" run directories for fullchem or TOMAS simulations, these restart files must be placed in the run directory or else there will be a species not found error due to the HEMCO_Config.rc "EFYO" time cycle flag. Signed-off-by: Bob Yantosca --- run/GCClassic/createRunDir.sh | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/run/GCClassic/createRunDir.sh b/run/GCClassic/createRunDir.sh index ca29075cc..815e330ed 100755 --- a/run/GCClassic/createRunDir.sh +++ b/run/GCClassic/createRunDir.sh @@ -1061,22 +1061,20 @@ if [[ ${met_name} = "MERRA2" ]] || [[ ${met_name} = "GEOSFP" ]]; then if [[ "x${sim_name}" == "xfullchem" || "x${sim_name}" == "xaerosol" ]]; then - # For TOMAS simulations, use restarts provided by the TOMAS team - # For other fullchem simulations, use restart from latest benchmark + # NOTE: We need to read the fullchem and TOMAS restart files from + # the v2021-09/ folder. These contain extra species (e.g HMS), + # for chemistry updates that were added in 13.3.0. This is necessary + # to avoid GEOS-Chem Classic simulations from halting if these + # species are not found in the restart file (time cycle flag "EFYO"). + # -- Bob Yantosca (22 Sep 2021) + # # Aerosol-only simulations can use the fullchem restart since all of the - # aerosol species are included + # aerosol species are included. if [[ "x${sim_extra_option}" == "xTOMAS15" ]]; then - sample_rst=${rst_root}/v2020-02/GEOSChem.Restart.TOMAS15.${startdate}_0000z.nc4 + sample_rst=${rst_root}/v2021-09/GEOSChem.Restart.TOMAS15.${startdate}_0000z.nc4 elif [[ "x${sim_extra_option}" == "xTOMAS40" ]]; then - sample_rst=${rst_root}/v2020-02/GEOSChem.Restart.TOMAS40.${startdate}_0000z.nc4 + sample_rst=${rst_root}/v2021-09/GEOSChem.Restart.TOMAS40.${startdate}_0000z.nc4 else - - # NOTE: We need to read the restart file in v2021-09 which - # contains the same species concentrations as in the GC_13.0.0 - # folder, with HSO3-, SO3--, HMS, C2H2, and C2H4 appended to - # it. Because HEMCO_Config.rc uses flag EFYO, the run will - # halt if not all species are found. (bmy, 9/21/21) - #sample_rst=${rst_root}/GC_13.0.0/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 sample_rst=${rst_root}/v2021-09/GEOSChem.Restart.fullchem.${startdate}_0000z.nc4 fi From d4aebf2b07d117c27a13ec3e81d5d4859970fdd8 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 12:12:17 -0400 Subject: [PATCH 12/16] Bug fix: Make sure HMS0 is initialized for aerosol-only simulations Prior to this commit, HMS0 was set to zero for aerosol-only simulations. Now it is set to Spc(I,J,L,id_HMS). This had caused aerosol-only integration tests to fail. Signed-off-by: Bob Yantosca --- GeosCore/sulfate_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index 69f8305a7..26d3a88f1 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -2669,10 +2669,10 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & HMS0 = Spc(I,J,L,id_HMS) !(jmm, 06/15/2018) OH0 = Spc(I,J,L,id_OH) !(jmm, 06/26/18) ELSE + HMS0 = Spc(I,J,L,id_HMS) !(jmm, 06/15/2018) HOBr0 = 0.e+0_fp HOCl0 = 0.e+0_fp HCHO0 = 0.e+0_fp !(jmm 06/13/2018) - HMS0 = 0.e+0_fp !(jmm, 06/15/2018) OH0 = 0.e+0_fp !(jmm, 06/26/18) ENDIF From 8bf24a15351d26eaeae9dffad82470c6d06f4fcf Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 14:07:16 -0400 Subject: [PATCH 13/16] Additional bug fixes for aerosol-only simulation in sulfate_mod.F90 (1) Make sure that HOBr, HOCl, HCHO, HMS, and OH are only referenced from state_chm%species if it is a fullchem simulation. (2) Removed some ELSE statements (i.e. zero variables before IF). (3) Expanded IS_FULLCHEM to Input_Opt%ITS_A_FULLCHEM_SIM (4) Trimmed trailing whitespace. Signed-off-by: Bob Yantosca --- GeosCore/sulfate_mod.F90 | 253 ++++++++++++++++++++------------------- 1 file changed, 127 insertions(+), 126 deletions(-) diff --git a/GeosCore/sulfate_mod.F90 b/GeosCore/sulfate_mod.F90 index 26d3a88f1..e0110581d 100644 --- a/GeosCore/sulfate_mod.F90 +++ b/GeosCore/sulfate_mod.F90 @@ -109,7 +109,7 @@ MODULE SULFATE_MOD ! PSO4_SO2 : P(SO4) from SO2 [v/v/timestep] ! PHMS_SO2 : P(HMS) from SO2 [v/v/timestep] (jmm, 6/15/18) ! PSO2_HMS : P(SO2) from HMS [v/v/timestep] (jmm, 6/15/18) - ! PSO4_HMS : P(SO4) from HMS & radical chem [v/v/timestep] (jmm, 6/26/18) + ! PSO4_HMS : P(SO4) from HMS & radical chem [v/v/timestep] (jmm, 6/26/18) ! SSTEMP : Sea surface temperatures [K] ! Eev : SO2 em. from eruptive volcanoes [kg SO2/box/s] ! Env : SO2 em. from non-erup volcanoes [kg SO2/box/s] @@ -130,8 +130,8 @@ MODULE SULFATE_MOD REAL(fp), ALLOCATABLE :: PSO4_SO2(:,:,:) REAL(fp), ALLOCATABLE :: PSO4_SS(:,:,:) REAL(fp), ALLOCATABLE :: PHMS_SO2(:,:,:) ! jmm 06/13/2018 - REAL(fp), ALLOCATABLE :: PSO2_HMS(:,:,:) ! jmm 06/13/2018 - REAL(fp), ALLOCATABLE :: PSO4_HMS(:,:,:) ! jmm 06/26/2018 + REAL(fp), ALLOCATABLE :: PSO2_HMS(:,:,:) ! jmm 06/13/2018 + REAL(fp), ALLOCATABLE :: PSO4_HMS(:,:,:) ! jmm 06/26/2018 REAL(fp), ALLOCATABLE :: PNITs(:,:,:) REAL(fp), ALLOCATABLE :: PNIT_dust(:,:,:,:) ! tdf REAL(fp), ALLOCATABLE :: PSO4_dust(:,:,:,:) ! tdf @@ -414,7 +414,7 @@ SUBROUTINE CHEMSULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, & PSO4_SO2 = 0e+0_fp PHMS_SO2 = 0e+0_fp ! jmm 06/13/2018 PSO2_HMS = 0e+0_fp ! jmm 06/13/2018 - PSO4_HMS = 0e+0_fp ! jmm 06/26/2018 + PSO4_HMS = 0e+0_fp ! jmm 06/26/2018 PSO4_SS = 0e+0_fp PNITs = 0e+0_fp PSO4_dust = 0e+0_fp ! tdf 04/17/08 @@ -2369,14 +2369,14 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! ============================================================================ ! (1 ) SO2 production: ! (a) DMS + OH, DMS + NO3 (saved in CHEM_DMS) -! (b) HMS -> SO2 + HCHO (aq) +! (b) HMS -> SO2 + HCHO (aq) ! . ! (2 ) SO2 loss: ! (a) SO2 + OH -> SO4 ! (b) SO2 -> drydep ! (c) SO2 + H2O2 or O3 (aq) -> SO4 -! (d) SO2 + HCHO (aq)-> HMS -! (d) SO2 + HMS -> 2 SO4 +! (d) SO2 + HCHO (aq)-> HMS +! (d) SO2 + HMS -> 2 SO4 ! . ! (3 ) SO2 = SO2_0 * exp(-bt) + PSO2_DMS/bt * [1-exp(-bt)] ! . @@ -2431,7 +2431,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fp) :: TFA, TAA, TDCA ! (jmm, 12/03/2018) REAL(fp) :: HCHO0, HMSc, KaqHCHO, KaqHMS ! (jmm, 06/07/2018) REAL(fp) :: L7, L7S, L7_b, L7S_b, HMS0 ! (jmm, 06/07/2018) - REAL(fp) :: L8, L8S, OH0, KaqHMS2 ! (jmm, 06/26/2018) + REAL(fp) :: L8, L8S, OH0, KaqHMS2 ! (jmm, 06/26/2018) REAL(fp) :: PSO4d_tot, PNITd_tot REAL(fp) :: SO2_gas, PH2SO4d_tot REAL(fp) :: H2SO4_cd, H2SO4_gas @@ -2451,7 +2451,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & REAL(fp) :: Fe_ant, Fe_nat, Fe_tot REAL(fp) :: Fe_d_ant, Fe_d_nat REAL(fp) :: IONIC - + REAL(fp) :: L6,L6S,SRhocl,L6_1,L6S_1 !XW REAL(fp) :: SO4H3_vv, SO4H4_vv !XW @@ -2647,7 +2647,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & Ld = 0.0_fp LSTOT0 = 0.0_fp LSTOT = 0.0_fp - LSTOT_HMS = 0.0_fp + LSTOT_HMS = 0.0_fp ! Skip non-chemistry boxes IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE @@ -2658,38 +2658,36 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & one_m_KRATE = 1.0_fp - State_Chm%KRATE(I,J,L) #endif - ! Initial SO2, H2O2, HOBr, and O3 [v/v] + ! Initialize [v/v] SO20 = Spc(I,J,L,id_SO2) H2O20 = Spc(I,J,L,id_H2O2) - IF ( IS_FULLCHEM ) THEN - ! HOBr, OH, and HCHO are only a defined species in fullchem simulations - HOBr0 = Spc(I,J,L,id_HOBr) !(qjc, 01/30/17) - HOCl0 = Spc(I,J,L,id_HOCl) !XW - HCHO0 = Spc(I,J,L,id_CH2O) !(jmm, 06/13/2018) - HMS0 = Spc(I,J,L,id_HMS) !(jmm, 06/15/2018) - OH0 = Spc(I,J,L,id_OH) !(jmm, 06/26/18) - ELSE - HMS0 = Spc(I,J,L,id_HMS) !(jmm, 06/15/2018) - HOBr0 = 0.e+0_fp - HOCl0 = 0.e+0_fp - HCHO0 = 0.e+0_fp !(jmm 06/13/2018) - OH0 = 0.e+0_fp !(jmm, 06/26/18) + + ! These species are only needed for fullchem simulations + HOBr0 = 0.0_fp + HOCl0 = 0.0_fp + HCHO0 = 0.0_fp + HMS0 = 0.0_fp + OH0 = 0.0_fp + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN + HOBr0 = Spc(I,J,L,id_HOBr) + HOCl0 = Spc(I,J,L,id_HOCl) + HCHO0 = Spc(I,J,L,id_CH2O) + HMS0 = Spc(I,J,L,id_HMS) + OH0 = Spc(I,J,L,id_OH) ENDIF ! Calculate O3, defined only in the chemistry grid IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN ! Get O3 from State_Chm%Species [v/v] + O3 = 0.0_fp IF ( State_Met%InChemGrid(I,J,L) ) THEN O3 = State_Chm%Species(I,J,L,id_O3) - ELSE - O3 = 0e+0_fp ENDIF ELSE IF ( Input_Opt%ITS_AN_AEROSOL_SIM ) THEN ! Get offline mean O3 [v/v] for this gridbox and month + O3 = 0.0_fp IF ( L <= State_Grid%MaxChemLev ) THEN O3 = O3m(I,J,L) - ELSE - O3 = 0e+0_fp ENDIF ENDIF @@ -2786,7 +2784,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & CALL GET_ALK( I, J, L, ALK1, & ALK2, Kt1, Kt2, Kt1N, & Kt2N, Kt1L, Kt2L, Input_Opt, & - State_Chm, State_Grid, State_Met, RC ) + State_Chm, State_Grid, State_Met, RC ) ! Total alkalinity [kg] ALK = ALK1 + ALK2 @@ -2843,7 +2841,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! We also need to update the SALAAL and SALCAL species here ! when we are using KPP to perform sulfur chemistry, otherwise ! SALAAL and SALCAL will not get set properly. Think about - ! a better way to abstract this code later. + ! a better way to abstract this code later. ! -- Bob Yantosca (09 Sep 2021) IF ( .not. State_Chm%Do_SulfateMod_SeaSalt ) THEN Spc(I,J,L,id_SALAAL) = AlkA @@ -2998,7 +2996,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L7S = 0.e+0_fp !(jmm, 06/13/18) L7S_b = 0.e+0_fp !(jmm, 06/13/18) L8 = 0.e+0_fp !(jmm, 06/26/18) - L8S = 0.e+0_fp !(jmm, 06/26/18) + L8S = 0.e+0_fp !(jmm, 06/26/18) ! If we are not using KPP to compute reaction rates, then ! compute SO2 cloud chemistry (1) if there is cloud, (2) if there @@ -3100,17 +3098,17 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ( AIRMW * LWC ) ) #endif - IF ( IS_FULLCHEM ) THEN - ! Get HMS cloud concentration and convert from [v/v] to - ! [moles/liter] (jmm, 06/13/2018) - ! Use a cloud scavenging ratio of 0.7 - ! assume nonvolatile like sulfate for realistic cloud pH + ! Get HMS cloud concentration and convert from [v/v] to + ! [moles/liter] (jmm, 06/13/2018) + ! Use a cloud scavenging ratio of 0.7 + ! assume nonvolatile like sulfate for realistic cloud pH + ! NOTE: Only needed for fullchem sims, otherwise it's zero + HMSc = 0.0_fp + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN HMSc = Spc(I,J,L,id_HMS) * State_Met%AIRDEN(I,J,L) * & 0.7e+0_fp / ( AIRMW * LWC ) - ELSE - HMSc = 0e+0_fp ENDIF - + ! Get total ammonia (NH3 + NH4+) concentration [v/v] ! Use a cloud scavenging ratio of 0.7 for NH4+ #ifdef LUO_WETDEP @@ -3205,7 +3203,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Get total nitrate (HNO3 + NIT) concentrations [v/v] ! Use a cloud scavenging ratio of 0.7 for NIT - IF ( IS_FULLCHEM ) THEN + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN #ifdef LUO_WETDEP TNO3 = ( Spc(I,J,L,id_HNO3) & + Spc(I,J,L,id_NIT ) & @@ -3334,13 +3332,13 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & IONIC = 0.5e+0_fp * ( 4.0e+0_fp * SO4nss + & TNH3*State_Met%AIRDEN(I,J,L)/(28.97e+0_fp*LWC) + & TNO3*State_Met%AIRDEN(I,J,L)/(28.97e+0_fp*LWC) ) - + ! Compute aqueous rxn rates for SO2 CALL AQCHEM_SO2( LWC, TK, PATM, SO2_ss, H2O20, & O3, HPLUS, MnII, FeIII, IONIC, & KaqH2O2, KaqO3, KaqO3_1, KaqO2, & HSO3aq, SO3aq, HCHO0, KaqHCHO, & - KaqHMS, KaqHMS2 ) + KaqHMS, KaqHMS2 ) !---------------------------------------------------------- ! Compute loss by H2O2. Prevent floating-point exception @@ -3482,16 +3480,15 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !---------------------------------------------------------- ! Get SO4H (sulfate produced via HOBr) from the Spc array [v/v dry] - IF ( IS_FULLCHEM ) THEN + SO4H1_vv = 0.e+0_fp + SO4H2_vv = 0.e+0_fp + SO4H3_vv = 0.e+0_fp + SO4H4_vv = 0.e+0_fp + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN SO4H1_vv = Spc(I,J,L,id_SO4H1) SO4H2_vv = Spc(I,J,L,id_SO4H2) SO4H3_vv = Spc(I,J,L,id_SO4H3) SO4H4_vv = Spc(I,J,L,id_SO4H4) - ELSE - SO4H1_vv = 0.e+0_fp - SO4H2_vv = 0.e+0_fp - SO4H3_vv = 0.e+0_fp - SO4H4_vv = 0.e+0_fp ENDIF L5S = SO4H1_vv + SO4H2_vv @@ -3526,21 +3523,21 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! by not allowing the exponential to go to infinity if ! the argument is too large. (jmm, 06/13/18) !---------------------------------------------------------- - IF (IS_FULLCHEM ) THEN + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN ! Argument of the exponential XX = ( SO2_ss - HCHO0 ) * KaqHCHO * DTCHEM - + ! Test if EXP(XX) can be computed w/o numerical exception IF ( IS_SAFE_EXP( XX ) .and. ABS( XX ) > 0e+0_fp ) THEN - + ! Aqueous phase SO2 loss rate w/ HCHO [v/v/timestep] L7 = EXP( XX ) - + ! Loss by HCHO L7S = SO2_ss * HCHO0 * ( L7 - 1.e+0_fp ) / & ( (SO2_ss * L7) - HCHO0 ) ELSE - + ! NOTE from Jintai Lin (4/28/10): ! However, in the case of a negative XX, L7S should be ! approximated as SO2_ss, instead of HCHO0. In other words, @@ -3555,9 +3552,9 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & !(qjc, 04/10/16) different solution when SO2_ss = H2O20 L7S = SO2_ss - 1/(KaqHCHO*DTCHEM+1/SO2_ss) ENDIF - + ENDIF - + IF ( L7S > 1.0e+15_fp .or. L7S < 0 ) THEN RC = GC_FAILURE PRINT *,'Loc:',I,J,L @@ -3568,46 +3565,46 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & PRINT *,'L7:',L7 PRINT *,'KaqHCHO',KaqHCHO ENDIF - - + + !---------------------------------------------------------- - ! Compute aqueous SO2 from HMS decomposition. Prevent floating-point + ! Compute aqueous SO2 from HMS decomposition. Prevent floating-point ! exception by not allowing the exponential to go to infinity if ! the argument is too large. This reaction also produces HCHO. ! (jmm, 06/13/18) ! - ! This is treated as a pseudo first order reaction as it depends - ! on [HMS] and [OH-]. The [OH-] is incorporated into the rate + ! This is treated as a pseudo first order reaction as it depends + ! on [HMS] and [OH-]. The [OH-] is incorporated into the rate ! constant (KaqHMS) with HPLUS determined seperately. Updates to - ! [OH-] are incorporated in the next call of GET_HPLUS, which + ! [OH-] are incorporated in the next call of GET_HPLUS, which ! includes [HMS] in the pH calculation - ! - ! Also note that [HMS] is in mol/L, and the conversion for the - ! rate from [mol/L/s] to [v/v/s] is done here + ! + ! Also note that [HMS] is in mol/L, and the conversion for the + ! rate from [mol/L/s] to [v/v/s] is done here ! ! Followed proceedure as done with oxidation by O2 (metal catalyzed) !---------------------------------------------------------- - + ! Argument of the exponential XX = -KaqHMS * DTCHEM - + ! Test if EXP(XX) can be computed w/o numerical exception IF ( IS_SAFE_EXP( XX ) ) THEN - + ! Aqueous phase SO2 production rate w/ HMS [v/v/timestep] L7_b = EXP( XX ) - + ! Production by HMS L7S_b = HMSc * ( 1.e+0_fp - L7_b ) ELSE L7S_b = HMSc ENDIF - + ! Convert HMS produced from [mol/L] to [v/v] (jmm, 06/15/18) L7S_b = L7S_b *LWC * 0.08205e+0_fp * TK / PATM - - - + + + IF ( L7S_b > 1.0e+15_fp .or. L7S_b < 0 ) THEN RC = GC_FAILURE PRINT *,'Loc: ',I,J,L @@ -3621,9 +3618,11 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & PRINT *,'PATM :',PATM ! CALL ERROR_STOP( 'L6s_b >1e15 or <0, point 3', LOC ) ENDIF - - - + + L7S = 0.e+0_fp + L7S_b = 0.e+0_fp + L8S = 0.e+0_fp + !---------------------------------------------------------- ! Compute HMS loss by aqueous OH. Prevent floating-point exception ! by not allowing the exponential to go to infinity if @@ -3631,28 +3630,28 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! ! HMS oxidation by aqueous OH is converted to v/v/s in the KaqHMS2 ! term. See aqchem_so2 subroutine for reaction details. - ! + ! ! This is a simplified version of a radical chain. ! L8S consumes 1 HMS, 1 OH, and 1 SO2 and produces 2 SO4-- !---------------------------------------------------------- - + ! Convert KaqHMS2 from [m^6 kg^-2 s^-1] to [v/v/s] KaqHMS2 = KaqHMS2 * State_Met%AIRDEN(I,J,L) * & State_Met%AIRDEN(I,J,L) ! Argument of the exponential XX = ( HMS0 - OH0 ) * KaqHMS2 * DTCHEM - + ! Test if EXP(XX) can be computed w/o numerical exception IF ( IS_SAFE_EXP( XX ) .and. ABS( XX ) > 0e+0_fp ) THEN - + ! Aqueous phase HMS loss rate w/ OH(aq) [v/v/timestep] L8 = EXP( XX ) - + ! Loss by OH(aq) L8S = HMS0 * OH0 * ( L8 - 1.e+0_fp ) / & ( (HMS0 * L8) - OH0 ) ELSE - + ! NOTE from Jintai Lin (4/28/10): ! However, in the case of a negative XX, L8S should be ! approximated as , instead of HCHO0. In other words, @@ -3669,15 +3668,14 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ENDIF ENDIF - + ELSE L7S = 0.e+0_fp L7S_b = 0.e+0_fp L8S = 0.e+0_fp ENDIF - - - + + L2S = L2S * FC L3S = L3S * FC L3S_1 = L3S_1 * FC !(qjc, 06/20/16) @@ -3687,7 +3685,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L5S_1 = L5S_1 !(qjc, 06/20/16) L6S = L6S L6S_1 = L6S_1 - + L7S = L7S * FC ! Note: consumes SO2 but for HMS (jmm, 06/13/18) L7S_b = L7S_b * FC ! Note: releases SO2 from HMS (jmm, 06/13/18) L8S = L8S * FC ! Notes: releases 2 sulfate and 1 HCHO for 1 HMS and 1 SO2 (jmm, 06/29/18) @@ -3711,7 +3709,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L7S_b = L7S_b L8S = L8S ENDIF - + LSTOT0 = L2S + L3S + L4S + L5S + L6S + L7S - L7S_b + L8S ! (qjc, 11/04/16), (jmm, 06/13/18) ! make sure sulfate produced is less than SO2 available @@ -3727,7 +3725,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L6S_1 = SO2_ss * L6S_1 / LSTOT0 L7S = SO2_SS * L7S / LSTOT0 ! (jmm, 06/13/18) L8S = SO2_SS * L8S / LSTOT0 ! (jmm, 06/29/18) - + ! This is the ratio used to calculate the actual removal of ! HOBr by SO2 for use in gckpp_HetRates.F90 (qjc, 06/20/16) @@ -3746,7 +3744,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & L6S_1 = L6S_1 L7S = L7S ! (jmm, 06/13/18) L8S = L8S ! (jmm, 06/29/18) - + ! This is the ratio used to calculate the actual removal of ! HOBr by SO2 (qjc, 06/20/16) @@ -3917,7 +3915,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Get total in-cloud sulfate production based on bulk cloud pH ! calculations for use in HET_DROP_CHEM - ! added in SO4 from HMS (jmm, 06/29/18) + ! added in SO4 from HMS (jmm, 06/29/18) LSTOT = (L2S + L3S + L4S + L5S + L6S + L8S + L8S) / FC !(qjc, 06/20/16) ! Get coarse-mode sea-salt concentration for use in @@ -3932,7 +3930,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! This is to make sure HET_DROP_CHEM does not compute more ! sulfate then there is SO2 ! Add L5S (qjc, 11/04/16) - ! Add L7S, L7S_b, and L8S (jmm, 06/29/18) + ! Add L7S, L7S_b, and L8S (jmm, 06/29/18) SO2_sr = MAX( SO2_ss - ( L2S + L3S + L4S + L5S + L6S + & L7S - L7S_b + L8S), MINDAT) @@ -3992,7 +3990,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ! Add L5S (qjc, 11/04/16) ! Add L7S and L7S_b (jmm, 06/13/18) ! Add loss and production of HCHO (jmm, 06/13/18) - ! Add loss and production of HMS (jmm, 06/15/18) + ! Add loss and production of HMS (jmm, 06/15/18) SO2_ss = MAX( SO2_ss - ( L2S+L3S+L4S+L5S+L6S+ & L7S-L7S_b+L8S+SR ), MINDAT ) H2O20 = MAX( H2O20 - L2S, MINDAT ) @@ -4000,7 +3998,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & HMS0 = MAX( HMS0 + L7S - L7S_b - L8S, MINDAT ) O3 = MAX( O3 - L3S, MINDAT ) OH0 = MAX( OH0 - L8S, MINDAT ) - + ! Factor to calculate effective SO3 and HSO3 in aqchem_so2 to ! be used in HOBr+HSO3/SO3 (qjc, 06/20/16) @@ -4017,7 +4015,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & H2O2s(I,J,L) = H2O20 ! SO2 chemical loss rate = SO4 production rate [v/v/timestep] - ! Add HMS (jmm, 06/13/18) + ! Add HMS (jmm, 06/13/18) PSO4_SO2(I,J,L) = LSTOT + PSO4E PSO4_ss (I,J,L) = PSO4F PHMS_SO2(I,J,L) = L7S @@ -4046,7 +4044,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & HPLUS = 0.e+0_fp L7S = 0.e+0_fp !(jmm, 06/13/18) L7S_b = 0.e+0_fp !(jmm, 06/13/18) - L8S = 0.e+0_fp !(jmm, 06/13/18) + L8S = 0.e+0_fp !(jmm, 06/13/18) ! Store HSO3aq, SO3aq for use in gckpp_HetRates.F90 ! Avoid divide-by-zero errors @@ -4059,25 +4057,28 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & State_Chm%fupdateHOCl(I,J,L) = 0.e+0_fp ! SO2 chemical loss rate = SO4 production rate [v/v/timestep] - ! Add HMS (jmm, 06/13/18) + ! Add HMS (jmm, 06/13/18) PSO4_SO2(I,J,L) = PSO4E PSO4_ss (I,J,L) = PSO4F PHMS_SO2(I,J,L) = 0.0e+0_fp PSO2_HMS(I,J,L) = 0.0e+0_fp PSO4_HMS(I,J,L) = 0.0e+0_fp - + ENDIF ! Store updated SO2, H2O2 back to the tracer arrays - ! Add HCHO, OH, O3, and HMS (jmm, 06/13/18) + ! Add HCHO, OH, O3, and HMS (jmm, 06/13/18) If (FullRun) Then Spc(I,J,L,id_SO2) = SO2s( I,J,L) Spc(I,J,L,id_H2O2) = H2O2s(I,J,L) - Spc(I,J,L,id_CH2O) = HCHO0 - Spc(I,J,L,id_HMS) = HMS0 - Spc(I,J,L,id_OH) = OH0 - Spc(I,J,L,id_O3) = O3 + + IF ( Input_Opt%ITS_A_FULLCHEM_SIM ) THEN + Spc(I,J,L,id_OH) = OH0 + Spc(I,J,L,id_O3) = O3 + Spc(I,J,L,id_CH2O) = HCHO0 + Spc(I,J,L,id_HMS) = HMS0 + ENDIF ! Set SO4H1 and SO4H2 to zero at end of each timestep IF ( id_SO4H1 > 0 ) Spc(I,J,L,id_SO4H1) = 0.0e+0_fp @@ -4206,7 +4207,7 @@ SUBROUTINE CHEM_SO2( Input_Opt, State_Chm, State_Diag, State_Grid, & ( 2.e+0_fp * L8S * State_Met%AD(I,J,L) / & TCVV_S ) / DTCHEM ENDIF - + !----------------------------------------------------------- ! Diagnostics for acid uptake on dust aerosol simulations !----------------------------------------------------------- @@ -5334,7 +5335,7 @@ SUBROUTINE GET_HPLUS( SO4nss, HMsc, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & ! !INPUT PARAMETERS: ! REAL(fp), INTENT(IN) :: SO4nss ! Total nss sulfate mixing ratio [M] - REAL(fp), INTENT(IN) :: HMSc ! Total HMS mixing ratio [M] + REAL(fp), INTENT(IN) :: HMSc ! Total HMS mixing ratio [M] REAL(fp), INTENT(IN) :: TNO3 ! Total nitrate (gas+particulate) mixing ! ratio [v/v] REAL(fp), INTENT(IN) :: TNH3 ! NH3 mixing ratio [v/v] @@ -5357,7 +5358,7 @@ SUBROUTINE GET_HPLUS( SO4nss, HMsc, TNH3, TNO3, SO2, CL, TNA, TDCA, TFA, & ! ============================================================================ ! Solve the following electroneutrality equation: ! [H+] = 2[SO4--] + [Cl-] + [OH-] + [HCO3-] + 2[CO3--] + [HSO3-] + 2[SO3--] +! -! [NO3-] + [HCOO-] + [CH3COO-] +[HMS] - [Na] - 2[Ca] - [NH4] +! [NO3-] + [HCOO-] + [CH3COO-] +[HMS] - [Na] - 2[Ca] - [NH4] ! Uses Newton's method to solve the equation: ! x_1 = x_0 -f(x_0)/f'(x_0) ! iterate until converge @@ -6899,9 +6900,9 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp), INTENT(IN) :: HPLUS ! Concentration of H+ ion (i.e. pH) [v/v] REAL(fp), INTENT(IN) :: MnII ! Concentration of MnII [mole/l] REAL(fp), INTENT(IN) :: FeIII ! Concentration of FeIII [mole/l] - REAL(fp), INTENT(IN) :: IONIC ! Ionic strength [mole/l]? + REAL(fp), INTENT(IN) :: IONIC ! Ionic strength [mole/l]? REAL(fp), INTENT(IN) :: HCHO ! HCHO mixing ratio [v/v] (jmm, 06/13/18) - + ! ! !OUTPUT PARAMETERS: ! @@ -6931,7 +6932,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! d[S(VI)]/dt = (k0[SO2(aq)] + k1[HSO3-] + K2[SO3--])[O3(aq)] ! [Seinfeld and Pandis, 1998, page 363] ! . -! (R3) HSO3- + HCHO(aq) => HMS +! (R3) HSO3- + HCHO(aq) => HMS ! SO3-- + HCHO(aq) => HMS + OH- ! [Moch et al., 2018; Olson and Hoffman, 1986] ! . @@ -6943,7 +6944,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! (note treated as 1st order in contrast to other reactions here) ! ! (R5) HMS + OH(aq) =(SO2,HO2,O2)=> HCHO + 2SO4-- + O2 + 3H+ + 2H2O -! [Jacob et al, 1986, Olson and Fessenden, 1992; +! [Jacob et al, 1986, Olson and Fessenden, 1992; ! Seinfeld and Pandis, 2016, Table 7A.7] ! Net reaction (R5): ! HMS + OH(aq) =(O2)=> SO5- + HCHO + H2O @@ -6990,7 +6991,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! !DEFINED PARAMETERS: ! REAL(fp), PARAMETER :: R = 0.08205e+0_fp - REAL(fp), PARAMETER :: dOH = 1.0e-19_fp ! [M cm^3 molec^-1] (jmm, 06/28/18) + REAL(fp), PARAMETER :: dOH = 1.0e-19_fp ! [M cm^3 molec^-1] (jmm, 06/28/18) ! ! !LOCAL VARIABLES: ! @@ -7001,9 +7002,9 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & REAL(fp) :: HCO3, XO3g, O3aq, XHCHOg, HCHCHO ! (jmm, 06/13/18) REAL(fp) :: FHCHCHO,KHCHO1, KHCHO2, KHMS ! (jmm, 06/13/18) REAL(fp) :: KW1, KHC1, KHMS2 ! (jmm, 06/15/18) - REAL(fp) :: Eff_mn, Eff_fe !jys - - + REAL(fp) :: Eff_mn, Eff_fe !jys + + !================================================================= ! AQCHEM_SO2 begins here! @@ -7172,10 +7173,10 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & (2.51e+13_fp*Hplus**(0.67) * & (MnII*FeIII*Eff_fe*Eff_mn)) * & LWC * R * T ! [s-1] - ENDIF + ENDIF !================================================================= - ! Aqueous reactions of SO2 with HCHO + ! Aqueous reactions of SO2 with HCHO ! HSO3- + HCHO(aq) => HMS + OH- (1) ! SO3-- + HCHO(aq) => HMS (2) ! @@ -7184,11 +7185,11 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! KHCHO1 = 7.9E2 * EXP( -16.435 * ((298.15/T) - 1.)) ! KHCHO2 = 2.5E7 * EXP( -6.037 * ((298.15/T) - 1.)) ! - ! + ! ! Aqueous reaction of HMS with OH- ! HMS + OH- => HCHO(aq) + SO3-- (3) ! - ! NOTE: unclear where B (E/R) value in Seinfeld and Pandis from, + ! NOTE: unclear where B (E/R) value in Seinfeld and Pandis from, ! but close to Deister. Using Seinfeld and Pandis value for now ! [Deister et al., 1986] ! KHMS = 3.6E3 * EXP( -22.027 * ((298.15/T) - 1.)) @@ -7198,7 +7199,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! ! Aqueous reaction of HMS with OH(aq) ! HMS + OH(aq) =(SO2,O2,HO2)=> 2SO4-- + HCHO + O2 + 3H+ + 2H2O (4) - ! + ! ! NOTE: O2, SO2, and HO2 particpate in the stoichiometry but not kinetics. ! Assume steady state for sulfur radicals and the following reaction chain: ! HMS + OH(aq) =(O2)=> SO5- + HCHO + H2O [Olsen and Fessenden, 1992] @@ -7207,13 +7208,13 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & ! SO2(aq) <=> HSO3- + H+ ! H+ + OH- <=> H2O ! HSO5- + HSO3- => 2SO4-- + 2H+ [Jacob, 1986] - ! Instead of assuming Henry's law for OH, use the parameter from + ! Instead of assuming Henry's law for OH, use the parameter from ! Jacob et al, 2005 that relates gas phase OH to aqueous phase OH ! accounting for the HO2(aq)/O2- cylcing in cloud droplets: ! dOH = 1E-19 [M cm^3 molec^-1] ! [Olson and Fessenden, 1992] - ! KHMS2 = 6.2E8 * EXP( -5.03 * ((298.15/T) -1.)) - ! + ! KHMS2 = 6.2E8 * EXP( -5.03 * ((298.15/T) -1.)) + ! ! ! (jmm, 06/28/18) !================================================================= @@ -7228,7 +7229,7 @@ SUBROUTINE AQCHEM_SO2( LWC, T, P, SO2, H2O2, & !================================================================= ! HCHO equilibrium reaction - ! HCHO(aq) + H2O = HCH(OH)2 + ! HCHO(aq) + H2O = HCH(OH)2 ! ! Reaction rate dependent on Temperature is given ! H = A exp ( B (T./T - 1) ) @@ -8950,13 +8951,13 @@ SUBROUTINE INIT_SULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) CALL GC_CheckVar( 'sulfate_mod.F90:PSO2_HMS', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN PSO2_HMS = 0e+0_fp - + ALLOCATE( PSO4_HMS( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & STAT=RC ) CALL GC_CheckVar( 'sulfate_mod.F90:PSO4_HMS', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN PSO4_HMS = 0e+0_fp - + ALLOCATE( PSO4_ss( State_Grid%NX, State_Grid%NY, State_Grid%MaxChemLev ), & STAT=RC ) CALL GC_CheckVar( 'sulfate_mod.F90:PSO4_ss', 0, RC ) @@ -9142,8 +9143,8 @@ SUBROUTINE INIT_SULFATE( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) id_SO4H4 = Ind_('SO4H4' ) id_HCOOH = Ind_('HCOOH' ) ! (jmm, 12/3/18) id_ACTA = Ind_('ACTA' ) ! (jmm, 12/3/18) - id_CH2O = Ind_('CH2O' ) ! (jmm, 06/15/18) - id_HMS = Ind_('HMS' ) ! (jmm, 06/15/18) + id_CH2O = Ind_('CH2O' ) ! (jmm, 06/15/18) + id_HMS = Ind_('HMS' ) ! (jmm, 06/15/18) ! Define flags for species drydep indices DRYNITs = Ind_('NITs', 'D') From 6f5db0c7d556d42828487d423d0e829028e5a962 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 14:23:15 -0400 Subject: [PATCH 14/16] Remove HMS from input.geos.aerosol template file HMS is only needed for the fullchem simulation. Signed-off-by: Bob Yantosca --- run/GCClassic/input.geos.templates/input.geos.aerosol | 1 - 1 file changed, 1 deletion(-) diff --git a/run/GCClassic/input.geos.templates/input.geos.aerosol b/run/GCClassic/input.geos.templates/input.geos.aerosol index ef2faf50c..93a22dc15 100644 --- a/run/GCClassic/input.geos.templates/input.geos.aerosol +++ b/run/GCClassic/input.geos.templates/input.geos.aerosol @@ -33,7 +33,6 @@ Species name : DST2 Species name : DST3 Species name : DST4 Species name : H2O2 -Species name : HMS Species name : MSA Species name : NH3 Species name : NH4 From 5167c68d17345cb5ddb43d8870eca6b5c8a9c700 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 15:09:55 -0400 Subject: [PATCH 15/16] Integration test fix: Only request 47L for TOMAS simulations Modified the sequence fed to createRunDir.sh in script test/GCClassic/intTestCreate.sh to create TOMAS15 and TOMAS40 integration test rundirs with 47 levels instead of 72. Signed-off-by: Bob Yantosca --- test/GCClassic/intTestCreate.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/GCClassic/intTestCreate.sh b/test/GCClassic/intTestCreate.sh index 0780603cb..114f1f48d 100755 --- a/test/GCClassic/intTestCreate.sh +++ b/test/GCClassic/intTestCreate.sh @@ -303,10 +303,10 @@ dir="gc_4x5_fullchem_RRTMG_merra2" create_rundir "1\n8\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_TOMAS15_merra2" -create_rundir "1\n6\n1\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n6\n1\n1\n1\n2\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_TOMAS40_merra2" -create_rundir "1\n6\n2\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n6\n2\n1\n1\n2\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_Hg_merra2" create_rundir "5\n1\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} @@ -370,10 +370,10 @@ dir="gc_4x5_fullchem_RRTMG_geosfp" create_rundir "1\n8\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_TOMAS15_geosfp" -create_rundir "1\n6\n1\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n6\n1\n2\n1\n2\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_fullchem_TOMAS40_geosfp" -create_rundir "1\n6\n2\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} +create_rundir "1\n6\n2\n2\n1\n2\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} dir="gc_4x5_Hg_geosfp" create_rundir "5\n2\n1\n1\n${root}\n${dir}\nn\n" ${root} ${dir} ${log} From c744c29681d8da8a2fb6832461f7359a5bf53803 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 16:06:06 -0400 Subject: [PATCH 16/16] Add additional error traps to prevent seg faults for aerosol-only sims GeosCore/aerosol_mod.F90: - Only include HMS in SO4_NH4_NIT if id_HMS > 0 - Set the HMS array to 0 if id_HMS < 0 GeosCore/isorropia_mod.F90 - Add an IS_HMS flag to denote when id_HMS > 0 - If id_HMS = 0 and it's an aerosol-only sim, let run proceed - Only add HMS to TSO4 if IS_HMS = TRUE - Only add HMS to DEN_SAV if IS_HMS = TRUE Signed-off-by: Bob Yantosca --- GeosCore/aerosol_mod.F90 | 69 ++++++++++++++++++++++++------------ GeosCore/isorropiaII_mod.F90 | 37 ++++++++++++++----- 2 files changed, 74 insertions(+), 32 deletions(-) diff --git a/GeosCore/aerosol_mod.F90 b/GeosCore/aerosol_mod.F90 index 272d14c6a..54927bf4c 100644 --- a/GeosCore/aerosol_mod.F90 +++ b/GeosCore/aerosol_mod.F90 @@ -467,29 +467,46 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! distribution but are currently simply treated in the same ! way (size and optics) as all other sulfate aerosol (DAR ! 2013) - SO4_NH4_NIT(I,J,L) = ( Spc(I,J,L,id_SO4) & - + Spc(I,J,L,id_HMS) & - + Spc(I,J,L,id_NH4) & - + Spc(I,J,L,id_NIT) ) & - / AIRVOL(I,J,L) + + IF ( IS_HMS ) THEN + + !%%%%% Fullchem simulations: add contribution from HMS + SO4_NH4_NIT(I,J,L) = ( Spc(I,J,L,id_SO4) & + + Spc(I,J,L,id_HMS) & + + Spc(I,J,L,id_NH4) & + + Spc(I,J,L,id_NIT) ) & + / AIRVOL(I,J,L) + + HMS(I,J,L) = Spc(I,J,L,id_HMS) / AIRVOL(I,J,L) + + ELSE + + !%%%%% Aerosol-only simulations: Skip contribution from HMS + SO4_NH4_NIT(I,J,L) = ( Spc(I,J,L,id_SO4) & + + Spc(I,J,L,id_NH4) & + + Spc(I,J,L,id_NIT) ) & + / AIRVOL(I,J,L) + + HMS(I,J,L) = 0.0_fp + ENDIF + SO4(I,J,L) = Spc(I,J,L,id_SO4) / AIRVOL(I,J,L) - HMS(I,J,L) = Spc(I,J,L,id_HMS) / AIRVOL(I,J,L) NH4(I,J,L) = Spc(I,J,L,id_NH4) / AIRVOL(I,J,L) NIT(I,J,L) = Spc(I,J,L,id_NIT) / AIRVOL(I,J,L) SLA(I,J,L) = 0.0_fp SPA(I,J,L) = 0.0_fp + ELSE ! Tropospheric sulfate is zero in stratosphere SO4_NH4_NIT(I,J,L) = 0.0_fp SO4(I,J,L) = 0.0_fp - HMS(I,J,L) = 0.0_fp ! (jmm, 06/30/18) + HMS(I,J,L) = 0.0_fp ! (jmm, 06/30/18) NH4(I,J,L) = 0.0_fp NIT(I,J,L) = 0.0_fp SLA(I,J,L) = KG_STRAT_AER(I,J,L,1) / AIRVOL(I,J,L) SPA(I,J,L) = KG_STRAT_AER(I,J,L,2) / AIRVOL(I,J,L) - ENDIF ! Add error check for safe division (bmy, 4/7/15) @@ -497,17 +514,23 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & ! Save these fractions for partitioning of optics ! until later when these may be treated independently - FRAC_SNA(I,J,L,1) = ( ( Spc(I,J,L,id_SO4 ) + & - Spc(I,J,L,id_HMS ) ) & - / AIRVOL(I,J,L) ) & - / SO4_NH4_NIT(I,J,L) + ! Only use HMS if it is defined (for fullchem sims) + IF ( IS_HMS ) THEN + FRAC_SNA(I,J,L,1) = ( ( Spc(I,J,L,id_SO4 ) + & + Spc(I,J,L,id_HMS ) ) & + / AIRVOL(I,J,L) ) & + / SO4_NH4_NIT(I,J,L) + ELSE + FRAC_SNA(I,J,L,1) = ( Spc(I,J,L,id_SO4 ) / AIRVOL(I,J,L) ) & + / SO4_NH4_NIT(I,J,L) + ENDIF - FRAC_SNA(I,J,L,2) = ( ( Spc(I,J,L,id_NIT ) ) & - / AIRVOL(I,J,L) ) & - / SO4_NH4_NIT(I,J,L) - FRAC_SNA(I,J,L,3) = ( Spc(I,J,L,id_NH4) & - / AIRVOL(I,J,L) ) / SO4_NH4_NIT(I,J,L) + FRAC_SNA(I,J,L,2) = ( Spc(I,J,L,id_NIT) / AIRVOL(I,J,L) ) & + / SO4_NH4_NIT(I,J,L) + + FRAC_SNA(I,J,L,3) = ( Spc(I,J,L,id_NH4) / AIRVOL(I,J,L) ) & + / SO4_NH4_NIT(I,J,L) ELSE @@ -860,7 +883,7 @@ SUBROUTINE AEROSOL_CONC( Input_Opt, State_Chm, State_Diag, & PM25(I,J,L) = NH4(I,J,L) * SIA_GROWTH + & NIT(I,J,L) * SIA_GROWTH + & SO4(I,J,L) * SIA_GROWTH + & - HMS(I,J,L) * SIA_GROWTH + & ! (jmm, 06/30/18) + HMS(I,J,L) * SIA_GROWTH + & ! (jmm, 06/30/18) BCPI(I,J,L) + & BCPO(I,J,L) + & OCPO(I,J,L) + & @@ -1813,7 +1836,7 @@ SUBROUTINE RDAER( Input_Opt, State_Chm, State_Diag, State_Grid, State_Met, & ! Volume of wet aerosol is also: VWet = 4/3*pi * RWet**3 * n ! So RWet = ( 3*VWet / (4 pi n) )**(1/3) ! RWet = RDry * ( 1 + VH2O/Vdry )**(1/3) - + ! Wet effective radius, um ! Here assume the dry radius of the mixture = SNA REFF = RW(1) * min( 3d0, & @@ -2285,7 +2308,7 @@ SUBROUTINE Init_Aerosol( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) id_SALACL = Ind_( 'SALACL' ) id_SO4 = Ind_( 'SO4' ) id_SO4s = Ind_( 'SO4s' ) - id_HMS = Ind_( 'HMS' ) ! (jmm, 06/30/18) + id_HMS = Ind_( 'HMS' ) ! (jmm, 06/30/18) id_NITs = Ind_( 'NITs' ) id_POA1 = Ind_( 'POA1' ) id_POA2 = Ind_( 'POA2' ) @@ -2362,7 +2385,7 @@ SUBROUTINE Init_Aerosol( Input_Opt, State_Chm, State_Diag, State_Grid, RC ) CALL GC_CheckVar( 'aerosol_mod.F:HMS', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN HMS = 0.0_fp - + ALLOCATE( NH4( State_Grid%NX, State_Grid%NY, State_Grid%NZ ), STAT=RC ) CALL GC_CheckVar( 'aerosol_mod.F:NH4', 0, RC ) IF ( RC /= GC_SUCCESS ) RETURN @@ -2797,11 +2820,11 @@ SUBROUTINE Set_AerMass_Diagnostic( Input_Opt, State_Chm, State_Diag, & ! jmm 3/6/19 !-------------------------------------- IF ( State_Diag%Archive_AerMassHMS ) THEN - State_Diag%AerMassHMS(I,J,L) = HMS(I,J,L) * & + State_Diag%AerMassHMS(I,J,L) = HMS(I,J,L) * & kgm3_to_ugm3 ENDIF - + !-------------------------------------- ! AerMassSOAGX [ug/m3] !-------------------------------------- diff --git a/GeosCore/isorropiaII_mod.F90 b/GeosCore/isorropiaII_mod.F90 index eda2037a3..6824b7817 100644 --- a/GeosCore/isorropiaII_mod.F90 +++ b/GeosCore/isorropiaII_mod.F90 @@ -181,6 +181,7 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & ! ! SAVEd scalars LOGICAL, SAVE :: FIRST = .TRUE. + LOGICAL, SAVE :: IS_HMS = .FALSE. INTEGER, SAVE :: id_HNO3, id_NH3, id_NH4 INTEGER, SAVE :: id_NIT, id_SALA, id_SO4 INTEGER, SAVE :: id_HMS ! jmm 12/5/18 @@ -293,13 +294,16 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & id_SALAAL = Ind_('SALAAL') id_SALCAL = Ind_('SALCAL') + ! Set a flag if HMS is defined + IS_HMS = ( id_HMS > 0 ) + ! Make sure certain tracers are defined IF ( id_SO4 <= 0 ) THEN ErrMsg = 'SO4 is an undefined species!' CALL GC_Error( ErrMsg, RC, ThisLoc ) RETURN ENDIF - IF ( id_HMS <= 0 ) THEN + IF ( id_HMS <= 0 .and. Input_Opt%ITS_A_FULLCHEM_SIM ) THEN ErrMsg = 'HMS is an undefined species!' CALL GC_Error( ErrMsg, RC, ThisLoc ) RETURN @@ -590,9 +594,14 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & IF ( N == 1 ) THEN ! Total SO4 [mole/m3], also consider SO4s in SALA - TSO4 = (Spc(I,J,L,id_SO4)+Spc(I,J,L,id_SALA)*0.08_fp*AlkR) * & - 1.e+3_fp / ( 96.0_fp * VOL) & - + Spc(I,J,L,id_HMS) * 0.5e+3_fp / ( 111.0_fp * VOL ) + IF ( IS_HMS ) THEN + TSO4 = (Spc(I,J,L,id_SO4)+Spc(I,J,L,id_SALA)*0.08_fp*AlkR) * & + 1.e+3_fp / ( 96.0_fp * VOL) & + + Spc(I,J,L,id_HMS) * 0.5e+3_fp / ( 111.0_fp * VOL ) + ELSE + TSO4 = (Spc(I,J,L,id_SO4)+Spc(I,J,L,id_SALA)*0.08_fp*AlkR) * & + 1.e+3_fp / ( 96.0_fp * VOL) + ENDIF ! Total NH3 [mole/m3] @@ -946,11 +955,21 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, & + Spc(I,J,L,id_NH4) / 18.0_fp & + Spc(I,J,L,id_SALA) * 0.3061_fp / 23.0_fp ) - DEN_SAV = ( Spc(I,J,L,id_SO4) / 96.0_fp * 2.0_fp & - + Spc(I,J,L,id_HMS) / 111.0_fp & - + Spc(I,J,L,id_NIT) / 62.0_fp & - + HNO3_DEN / 63.0_fp & - + Spc(I,J,L,id_SALA) * 0.55_fp / 35.45_fp ) + ! HMS is only defined for fullchem is simulations, + ! so skip it if it is not a defined species + IF ( IS_HMS ) THEN + DEN_SAV = ( Spc(I,J,L,id_SO4) / 96.0_fp * 2.0_fp & + + Spc(I,J,L,id_HMS) / 111.0_fp & + + Spc(I,J,L,id_NIT) / 62.0_fp & + + HNO3_DEN / 63.0_fp & + + Spc(I,J,L,id_SALA) * 0.55_fp / 35.45_fp ) + ELSE + DEN_SAV = ( Spc(I,J,L,id_SO4) / 96.0_fp * 2.0_fp & + + Spc(I,J,L,id_NIT) / 62.0_fp & + + HNO3_DEN / 63.0_fp & + + Spc(I,J,L,id_SALA) * 0.55_fp / 35.45_fp ) + ENDIF + ENDIF ENDDO