Skip to content

Commit

Permalink
Merge PR #1569 (branch 'feature/TOMASwork' of https://github.com/Bett…
Browse files Browse the repository at this point in the history 
…yCroft/geos-chem) into dev/14.2.0

Resolved conflicts in:
	GeosCore/fullchem_mod.F90
	KPP/fullchem/fullchem.eqn
	KPP/fullchem/gckpp_Function.F90
	KPP/fullchem/gckpp_Integrator.F90
	KPP/fullchem/gckpp_Jacobian.F90
	KPP/fullchem/gckpp_JacobianSP.F90
	KPP/fullchem/gckpp_LinearAlgebra.F90
	KPP/fullchem/gckpp_Monitor.F90
	KPP/fullchem/gckpp_Parameters.F90
	KPP/fullchem/gckpp_Rates.F90
	KPP/fullchem/gckpp_Util.F90

Signed-off-by: Melissa Sulprizio <mpayer@seas.harvard.edu>
  • Loading branch information
@msulprizio
msulprizio committed Apr 5, 2023
2 parents db22280 + baa98d0 commit 1cf8098
Show file tree
Hide file tree
Showing 12 changed files with 25,153 additions and 25,112 deletions.
241 changes: 202 additions & 39 deletions GeosCore/fullchem_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,10 @@ 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
INTEGER :: id_PSO4 !bc, 19/01/2022
INTEGER :: id_PSO4
#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 @@ -134,6 +137,7 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
USE UCX_MOD, ONLY : UCX_H2SO4PHOT
#ifdef TOMAS
USE TOMAS_MOD, ONLY : H2SO4_RATE
USE TOMAS_MOD, ONLY : PSO4AQ_RATE
#endif
!
! !INPUT PARAMETERS:
Expand Down Expand Up @@ -751,6 +755,9 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
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 @@ -1330,49 +1337,49 @@ SUBROUTINE Do_FullChem( Input_Opt, State_Chm, State_Diag, &
!-----------------------------------------------------------------
! FOR TOMAS MICROPHYSICS:
!
! Obtain P/L with a unit [kg S] at each timestep for tracing
! Obtain P/L with a unit [kg S] for tracing
! gas-phase sulfur species production (SO2, SO4, MSA)
! (win, 8/4/09),
! Note: VAR (id_PSO4) units are [molec/cm3] (bc, 19/01/2022)
! (win, 8/4/09)
!-----------------------------------------------------------------

! Calculate H2SO4 production rate [kg H2SO4 box-1 s-1] in each
! time step (win, 8/4/09, bc, 19/01/2022)
H2SO4_RATE(I,J,L) = VAR (id_PSO4)/ AVO * 98.e-3_fp * &
State_Met%AIRVOL(I,J,L) * &
1.0e+6_fp / DT
! Calculate H2SO4 production rate [kg s-1] in each
! time step (win, 8/4/09)
H2SO4_RATE(I,J,L) = C(ind_PH2SO4) / AVO * 98.e-3_fp * &
State_Met%AIRVOL(I,J,L) * &
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
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

#ifdef MODEL_CESM
!---------------------------------------------------------------------
! Calculate H2SO4 production rate for coupling to CESM
! (interface to MAM4 nucleation)
!---------------------------------------------------------------------
DO F = 1, NFAM

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

! Calculate H2SO4 production rate [mol mol-1] in this timestep
! (hplin, 1/25/23)
IF ( TRIM(FAM_NAMES(F)) == 'PSO4' ) THEN
! mol/mol = molec cm-3 * g * mol(Air)-1 * kg g-1 * m-3 cm3 / (molec mol-1 * kg m-3)
! = mol/molAir
State_Chm%H2SO4_PRDR(I,J,L) = C(KppID) * AIRMW * 1e-3_fp * 1.0e+6_fp / &
(AVO * State_Met%AIRDEN(I,J,L))

IF ( State_Chm%H2SO4_PRDR(I,J,L) < 0.0d0) THEN
write(*,*) "H2SO4_PRDR negative in fullchem_mod.F90!!", &
I, J, L, "was:", State_Chm%H2SO4_PRDR(I,J,L), " setting to 0.0d0"
State_Chm%H2SO4_PRDR(I,J,L) = 0.0d0
ENDIF
ENDIF
ENDDO
! mol/mol = molec cm-3 * g * mol(Air)-1 * kg g-1 * m-3 cm3 / (molec mol-1 * kg m-3) = mol/molAir
State_Chm%H2SO4_PRDR(I,J,L) = C(id_PSO4) * AIRMW * 1e-3_fp * 1.0e+6_fp /&
(AVO * State_Met%AIRDEN(I,J,L))

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

#ifdef MODEL_GEOS
Expand Down Expand Up @@ -1578,6 +1585,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 @@ -1678,6 +1700,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 @@ -2607,6 +2767,12 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
id_SALC = Ind_( 'SALC' )
id_SALCAL = Ind_( 'SALCAL' )
id_SO4 = Ind_( 'SO4' )
#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 Expand Up @@ -2731,21 +2897,18 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )

ENDIF

#ifdef TOMAS
#ifdef MODEL_CESM
!--------------------------------------------------------------------
! If we are finding H2SO4_RATE from a fullchem
! simulation for the TOMAS simulation, throw an error if we cannot find
! simulation for the CESM, throw an error if we cannot find
! the PSO4 prod family in this KPP mechanism.
!--------------------------------------------------------------------
IF ( State_Diag%Archive_Prod ) THEN

IF ( id_PSO4 < 1 ) THEN
ErrMsg = 'Could not find PSO4 in list of KPP families! This ' // &
'is needed for TOMAS H2SO4_RATE!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
ENDIF
IF ( id_PSO4 < 1 ) THEN
ErrMsg = 'Could not find PSO4 in list of KPP families! This ' // &
'is needed for State_Chm%H2SO4_PRDR and coupling to CESM!'
CALL GC_Error( ErrMsg, RC, ThisLoc )
RETURN
ENDIF
#endif

!--------------------------------------------------------------------
Expand Down
6 changes: 6 additions & 0 deletions GeosCore/tomas_mod.F90
View file
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 @@ -6518,6 +6519,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 @@ -10653,6 +10658,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
Loading

0 comments on commit 1cf8098

Please sign in to comment.