Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bug fix for TOMAS sulfate production rates #1569

Merged
merged 4 commits into from
Apr 7, 2023
Merged

Bug fix for TOMAS sulfate production rates #1569

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d1c1af0
save initial updates for passing sulfate production to TOMAS
BettyCroft Dec 16, 2022
133816e
Updates for passing sulfate production rates from KPP to TOMAS - sepa…
BettyCroft Dec 19, 2022
a67a851
Code comment added for sulfate prod diag
BettyCroft Dec 20, 2022
baa98d0
Revert to orig createRunDir.sh for pull request
BettyCroft Dec 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
206 changes: 197 additions & 9 deletions GeosCore/fullchem_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@ MODULE FullChem_Mod
! Species ID flags (and logicals to denote if species are present)
INTEGER :: id_OH, id_HO2, id_O3P, id_O1D, id_CH4
INTEGER :: id_PCO, id_LCH4, id_NH3
#ifdef TOMAS
INTEGER :: id_NK5, id_NK8, id_NK10, id_NK20
#endif
#ifdef MODEL_GEOS
INTEGER :: id_O3
INTEGER :: id_A3O2, id_ATO2, id_B3O2, id_BRO2
Expand Down Expand Up @@ -131,6 +134,7 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
#ifdef TOMAS
#ifdef BPCH_DIAG
USE TOMAS_MOD, ONLY : H2SO4_RATE
USE TOMAS_MOD, ONLY : PSO4AQ_RATE
#endif
#endif
!
Expand Down Expand Up @@ -702,6 +706,8 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
C(N) = 0.0_dp
IF ( SpcId > 0 ) C(N) = State_Chm%Species(SpcID)%Conc(I,J,L)
ENDDO
C(ind_PH2SO4) = 0.0e+0_fp !for sulfate gas prod diag, BC,12/20/22
C(ind_PSO4AQ) = 0.0e+0_fp !for sulfate cloud prod diag, BC 12/20/22

!=====================================================================
! Start KPP main timer
Expand Down Expand Up @@ -964,7 +970,7 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
RCNTRL = 0.0_fp

! Initialize Hstart (the starting value of the integration step
! size with the value of Hnew (the last predicted but not yet
! size with the value of Hnew (the last predicted but not yet
! taken timestep) saved to the the restart file.
RCNTRL(3) = State_Chm%KPPHvalue(I,J,L)

Expand Down Expand Up @@ -1248,10 +1254,10 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
#ifdef BPCH_DIAG
#ifdef TOMAS
!always calculate rate for TOMAS
DO F = 1, NFAM
!! DO F = 1, NFAM

! Determine dummy species index in KPP
KppID = PL_Kpp_Id(F)
!! KppID = PL_Kpp_Id(F)

!-----------------------------------------------------------------
! FOR TOMAS MICROPHYSICS:
Expand All @@ -1263,19 +1269,42 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &

! Calculate H2SO4 production rate [kg s-1] in each
! time step (win, 8/4/09)
IF ( TRIM(FAM_NAMES(F)) == 'PSO4' ) THEN
!! IF ( TRIM(FAM_NAMES(F)) == 'PSO4' ) THEN
! Hard-coded MW
H2SO4_RATE(I,J,L) = C(KppID) / AVO * 98.e-3_fp * &
!! H2SO4_RATE(I,J,L) = C(KppID) / AVO * 98.e-3_fp * &
!! State_Met%AIRVOL(I,J,L) * &
!! 1.0e+6_fp / DT

!! IF ( H2SO4_RATE(I,J,L) < 0.0d0) THEN
!! write(*,*) "H2SO4_RATE negative in fullchem_mod.F90!!", &
!! I, J, L, "was:", H2SO4_RATE(I,J,L), " setting to 0.0d0"
!! H2SO4_RATE(I,J,L) = 0.0d0
!! ENDIF
!! ENDIF
!! ENDDO

!!print*,'Here is index' , ind_PH2SO4, ind_PSO4AQ

H2SO4_RATE(I,J,L) = C(ind_PH2SO4) / AVO * 98.e-3_fp * &
State_Met%AIRVOL(I,J,L) * &
1.0e+6_fp / DT
1.0e+6_fp / DT ! kg s-1 box-1

IF ( H2SO4_RATE(I,J,L) < 0.0d0) THEN
write(*,*) "H2SO4_RATE negative in fullchem_mod.F90!!", &
I, J, L, "was:", H2SO4_RATE(I,J,L), " setting to 0.0d0"
H2SO4_RATE(I,J,L) = 0.0d0
H2SO4_RATE(I,J,L) = 0.0d0
ENDIF


PSO4AQ_RATE(I,J,L) = C(ind_PSO4AQ) / AVO * 98.e-3_fp * &
State_Met%AIRVOL(I,J,L) * &
1.0e+6_fp ! kg per timestep box-1

IF ( PSO4AQ_RATE(I,J,L) < 0.0d0) THEN
write(*,*) "PSO4AQ_RATE negative in fullchem_mod.F90", &
I, J, L, "was:", PSO4AQ_RATE(I,J,L), " setting to 0.0d0"
PSO4AQ_RATE(I,J,L) = 0.0d0
ENDIF
ENDIF
ENDDO

#endif
#endif
Expand Down Expand Up @@ -1458,6 +1487,21 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
ENDIF
#endif


#ifdef TOMAS
!-----------------------------------------------------------------
! For TOMAS microphysics:
!
! SO4 from aqueous chemistry of SO2 (in-cloud oxidation)
!
! SO4 produced via aqueous chemistry is distributed onto 30-bin
! aerosol by TOMAS subroutine AQOXID.
!-----------------------------------------------------------------

CALL TOMAS_SO4_AQ( Input_Opt, State_Chm, State_Grid, State_Met, RC )
IF ( prtDebug ) CALL DEBUG_MSG( '### CHEMSULFATE: a TOMAS_SO4_AQ' )
#endif

!=======================================================================
! Archive OH, HO2, O1D, O3P concentrations after solver
!=======================================================================
Expand Down Expand Up @@ -1558,6 +1602,144 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &

END SUBROUTINE Do_FullChem
!EOC
#ifdef TOMAS
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: tomas_so4_aq
!
! !DESCRIPTION: Subroutine CHEM\_SO4\_AQ takes the SO4 produced via aqueous
! chemistry of SO2 and distribute onto the size-resolved aerosol number and
! sulfate mass as a part of the TOMAS aerosol microphysics module
! (win, 1/25/10)
!\\
!\\
! !INTERFACE:
!
SUBROUTINE TOMAS_SO4_AQ( Input_Opt, State_Chm, State_Grid, State_Met, RC )
!
! !USES:
!
USE ErrCode_Mod
USE ERROR_MOD
USE Input_Opt_Mod, ONLY : OptINput
USE State_Chm_Mod, ONLY : ChmState
USE State_Grid_Mod, ONLY : GrdState
USE State_Met_Mod, ONLY : MetState
USE TOMAS_MOD, ONLY : AQOXID, GETACTBIN,PSO4AQ_RATE
USE UnitConv_Mod, ONLY : Convert_Spc_Units
!
! !INPUT PARAMETERS:
!
TYPE(OptInput), INTENT(IN) :: Input_Opt ! Input Options object
TYPE(GrdState), INTENT(IN) :: State_Grid ! Grid State object
TYPE(MetState), INTENT(IN) :: State_Met ! Meteorology State object

! !INPUT/OUTPUT PARAMETERS:
!
TYPE(ChmState), INTENT(INOUT) :: State_Chm ! Chemistry State object
!
! !OUTPUT PARAMETERS:
!
INTEGER, INTENT(OUT) :: RC ! Success or failure?
!
! !REMARKS:
! NOTE: This subroutine is ignored unless we compile for TOMAS microphysics.
!
! !REVISION HISTORY:
! See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
INTEGER :: I, J, L
INTEGER :: k, binact1, binact2
INTEGER :: KMIN
REAL(fp) :: SO4OXID
CHARACTER(LEN=63) :: OrigUnit

!=================================================================
! TOMAS_SO4_AQ begins here!
!=================================================================

! Assume success
RC = GC_SUCCESS

! Convert species from to [kg]
CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, &
'kg', RC, OrigUnit=OrigUnit )
IF ( RC /= GC_SUCCESS ) THEN
CALL GC_Error('Unit conversion error', RC, &
'Start of TOMAS_SO4_AQ in sulfate_mod.F90')
RETURN
ENDIF

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( I, J, L ) &
!$OMP PRIVATE( KMIN, SO4OXID, BINACT1, BINACT2 ) &
!$OMP SCHEDULE( DYNAMIC )
DO L = 1, State_Grid%NZ
DO J = 1, State_Grid%NY
DO I = 1, State_Grid%NX

! Skip non-chemistry boxes
IF ( .not. State_Met%InChemGrid(I,J,L) ) CYCLE

SO4OXID = PSO4AQ_RATE(I,J,L)

! print*, 'Here is SO4OXID ',SO4OXID, I,J,L

! SO4OXID = PSO4_SO2AQ(I,J,L) * State_Met%AD(I,J,L) &
! / ( AIRMW / State_Chm%SpcData(id_SO4)%Info%MW_g ) ! convert v/v to kg/box

IF ( SO4OXID > 0e+0_fp ) THEN
! JKodros (6/2/15 - Set activating bin based on which TOMAS bin
!length being used)
#if defined( TOMAS12 )
CALL GETACTBIN( I, J, L, id_NK5, .TRUE. , BINACT1, State_Chm, RC )

CALL GETACTBIN( I, J, L, id_NK5, .FALSE., BINACT2, State_Chm, RC )
#elif defined( TOMAS15 )
CALL GETACTBIN( I, J, L, id_NK8, .TRUE. , BINACT1, State_Chm, RC )

CALL GETACTBIN( I, J, L, id_NK8, .FALSE., BINACT2, State_Chm, RC )
#elif defined( TOMAS30 )
CALL GETACTBIN( I, J, L, id_NK10, .TRUE. , BINACT1, State_Chm, RC )

CALL GETACTBIN( I, J, L, id_NK10, .FALSE., BINACT2, State_Chm, RC )
#else
CALL GETACTBIN( I, J, L, id_NK20, .TRUE. , BINACT1, State_Chm, RC )

CALL GETACTBIN( I, J, L, id_NK20, .FALSE., BINACT2, State_Chm, RC )
#endif

KMIN = ( BINACT1 + BINACT2 )/ 2.

CALL AQOXID( SO4OXID, KMIN, I, J, L, Input_Opt, &
State_Chm, State_Grid, State_Met, RC )
ENDIF
ENDDO
ENDDO
ENDDO
!$OMP END PARALLEL DO

! Convert species back to original units
CALL Convert_Spc_Units( Input_Opt, State_Chm, State_Grid, State_Met, &
OrigUnit, RC )
IF ( RC /= GC_SUCCESS ) THEN
CALL GC_Error('Unit conversion error', RC, &
'End of TOMAS_SO4_AQ in sulfate_mod.F90')
RETURN
ENDIF

END SUBROUTINE TOMAS_SO4_AQ
!EOC
#endif
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
Expand Down Expand Up @@ -2488,6 +2670,12 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
id_SALCAL = Ind_( 'SALCAL' )
id_SO4 = Ind_( 'SO4' )
id_SALC = Ind_( 'SALC' )
#ifdef TOMAS
id_NK5 = Ind_( 'NK5' )
id_NK8 = Ind_( 'NK8' )
id_NK10 = Ind_( 'NK10' )
id_NK20 = Ind_( 'NK20' )
#endif

#ifdef MODEL_GEOS
! ckeller
Expand Down
6 changes: 6 additions & 0 deletions GeosCore/tomas_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,7 @@ MODULE TOMAS_MOD
INTEGER :: lowRH = 1 !This is to match AW more with AERONET (JKODROS 6/15)

REAL(fp), ALLOCATABLE :: H2SO4_RATE(:,:,:) ! H2SO4 prod rate [kg s-1]
REAL(fp), ALLOCATABLE :: PSO4AQ_RATE(:,:,:) ! Cld chem sulfate prod rate [kg s-1]

#if defined( TOMAS12 ) || defined( TOMAS15 )
!tomas12 or tomas15
Expand Down Expand Up @@ -6344,6 +6345,10 @@ SUBROUTINE INIT_TOMAS( Input_Opt, State_Chm, State_Grid, RC )
IF ( AS /= 0 ) CALL ALLOC_ERR( 'H2SO4_RATE' )
H2SO4_RATE = 0.0e+0_fp

ALLOCATE( PSO4AQ_RATE(State_Grid%NX,State_Grid%NY,State_Grid%NZ), STAT=AS )
IF ( AS /= 0 ) CALL ALLOC_ERR( 'PSO4AQ_RATE' )
PSO4AQ_RATE = 0.0e+0_fp

!=================================================================
! Calculate aerosol size bin boundaries (dry mass / particle)
!=================================================================
Expand Down Expand Up @@ -10459,6 +10464,7 @@ SUBROUTINE CLEANUP_TOMAS
!=================================================================
IF ( ALLOCATED( MOLWT ) ) DEALLOCATE( MOLWT )
IF ( ALLOCATED( H2SO4_RATE ) ) DEALLOCATE( H2SO4_RATE )
IF ( ALLOCATED( PSO4AQ_RATE ) ) DEALLOCATE( PSO4AQ_RATE )

END SUBROUTINE CLEANUP_TOMAS
!EOC
Expand Down
16 changes: 9 additions & 7 deletions KPP/custom/custom.eqn
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -463,6 +463,8 @@ TOLU = IGNORE; {C7H8; Toluene}
TRO2 = IGNORE; {Peroxy radical from TOLU oxidation}
XYLE = IGNORE; {C8H10; Xylene}
XRO2 = IGNORE; {Peroxy radical from XYLE oxidation}
PH2SO4 = IGNORE; {SO4 from gas-phase chemistry}
PSO4AQ = IGNORE; {SO4 from cloud chemistry}

#DEFFIX

Expand All @@ -486,14 +488,14 @@ HNO3 + SALCAL = NITs : K_MT(6);
//
// Cloud
// S(IV) --> S(VI)
SO2 + H2O2 = SO4 : K_CLD(1);
SO2 + O3 = SO4 : K_CLD(2);
SO2 {+O2} = SO4 : K_CLD(3); {Mn & Fe catalysis + HET_DROP_CHEM()}
SO2 + H2O2 = SO4 + PSO4AQ : K_CLD(1);
SO2 + O3 = SO4 + PSO4AQ : K_CLD(2);
SO2 {+O2} = SO4 + PSO4AQ : K_CLD(3); {Mn & Fe catalysis + HET_DROP_CHEM()}
//
// HMS
CH2O + SO2 = HMS : K_CLD(4); {Sep 2021; Moch2020; MSL}
HMS = SO2 + CH2O : K_CLD(5); {Sep 2021; Moch2020; MSL}
HMS + OH = 2SO4 + CH2O - SO2 : K_CLD(6); {Sep 2021; Moch2020; MSL}
HMS + OH = 2SO4 + CH2O - SO2 + 2PSO4AQ : K_CLD(6); {Sep 2021; Moch2020; MSL}
//
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// %%%%% Gas-phase chemistry reactions %%%%%
Expand Down Expand Up @@ -721,7 +723,7 @@ MPN {+M} = MO2 + NO2 : GCJPLPR_abcabc(1.05d-02, 4.8d+00, -
DMS + OH = SO2 + MO2 + CH2O : GCARR_ac(1.20d-11, -280.0d0);
DMS + OH = 0.750SO2 + 0.250MSA + MO2 : GC_DMSOH_acac(8.20d-39, 5376.0d0, 1.05d-5, 3644.0d0);
DMS + NO3 = SO2 + HNO3 + MO2 + CH2O : GCARR_ac(1.90d-13, 530.0d0);
SO2 + OH {+M} = SO4 + HO2 : GCJPLPR_aba(3.30d-31, 4.3d+00, 1.6d-12, 0.6d0);
SO2 + OH {+M} = SO4 + HO2 + PH2SO4 : GCJPLPR_aba(3.30d-31, 4.3d+00, 1.6d-12, 0.6d0);
Br + O3 = BrO + O2 : GCARR_ac(1.60d-11, -780.0d0); {2012/06/07; Parrella2012; JPP}
BrO + HO2 = HOBr + O2 : GCARR_ac(4.50d-12, 460.0d0); {2012/06/07; Parrella2012; JPP}
Br + HO2 = HBr + O2 : GCARR_ac(4.80d-12, -310.0d0); {2012/06/07; Parrella2012; JPP}
Expand Down Expand Up @@ -918,11 +920,11 @@ CH2OO + H2O = 0.730HMHP + 0.210HCOOH +
CH2OO + H2O + H2O = 0.400HMHP +
0.540HCOOH + 0.060CH2O + 0.060H2O2 : GCARR_ac(2.88d-35, 1391.0d0); {2019/11/06; Bates2019; KHB}
CH2OO + O3 = CH2O : 1.40d-12; {2019/11/06; Bates2019; KHB}
CH2OO + SO2 = CH2O + SO4 : 3.70d-11; {2019/11/06; Bates2019; KHB}
CH2OO + SO2 = CH2O + SO4 + PH2SO4: 3.70d-11; {2019/11/06; Bates2019; KHB}
CH3CHOO + CO = ALD2 : 1.20d-15; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + NO = ALD2 + NO2 : 1.00d-14; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + NO2 = ALD2 + NO3 : 1.00d-15; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + SO2 = ALD2 + SO4 : 7.00d-14; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + SO2 = ALD2 + SO4 + PH2SO4 : 7.00d-14; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + H2O = ALD2 + H2O2 : 6.00d-18; {2015/09/25; Millet2015; DBM,EAM}
CH3CHOO + H2O = ACTA : 1.00d-17; {2015/09/25; Millet2015; DBM,EAM}
MTPA + OH = PIO2 : GCARR_ac(1.21d-11, 440.0d0); {2017/03/23; IUPAC2010; EVF}
Expand Down
Loading