Skip to content

Commit

Permalink
Add additional error traps to prevent seg faults for aerosol-only sims
Browse files Browse the repository at this point in the history 
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 <yantosca@seas.harvard.edu>
  • Loading branch information
@yantosca
yantosca committed Sep 22, 2021
1 parent 5167c68 commit c744c29
Show file tree
Hide file tree
Showing 2 changed files with 74 additions and 32 deletions.
69 changes: 46 additions & 23 deletions GeosCore/aerosol_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -467,47 +467,70 @@ 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)
IF ( SO4_NH4_NIT(I,J,L) > 0e+0_fp ) THEN

! 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

Expand Down Expand Up @@ -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) + &
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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' )
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
!--------------------------------------
Expand Down
37 changes: 28 additions & 9 deletions GeosCore/isorropiaII_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit c744c29

Please sign in to comment.