From c744c29681d8da8a2fb6832461f7359a5bf53803 Mon Sep 17 00:00:00 2001 From: Bob Yantosca Date: Wed, 22 Sep 2021 16:06:06 -0400 Subject: [PATCH] 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