Merge pull request #1993 from jimmielin/hplin/14.1.1+cesmfix

Updates for GEOS-Chem within CESM: convective scavenging correction and MPI fixes

CHANGELOG.md

@@ -5,6 +5,12 @@ This file documents all notable changes to the GEOS-Chem repository starting in
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased 14.1.2]
+### Added
+- CESM-only update: Added option for correctConvUTLS for correcting buildup of soluble tracers in the UT/LS to match CAM-chem behavior
+
+### Changed
+- CESM-only update: extend existing KppError, KppStop to CESM for model stability
+- CESM-only update: Removed mpi_bcast in ucx_mod NOXCOEFF_INIT to be handled at coupler level to support spectral-element dynamical core
 
 ## [14.1.1] - 2023-03-03
 ### Added

GeosCore/fullchem_mod.F90

@@ -49,6 +49,13 @@ MODULE FullChem_Mod
   INTEGER               :: id_ISOPNOO1, id_ISOPNOO2, id_INO2B, id_INO2D
   INTEGER               :: id_IDHNBOO, id_IDHNDOO1, id_IDHNDOO2
   INTEGER               :: id_IHPNBOO, id_IHPNDOO, id_ICNOO, id_IDNOO
+#endif
+#ifdef MODEL_CESM
+  INTEGER               :: id_TSOA0, id_TSOA1, id_TSOA2, id_TSOA3
+  INTEGER               :: id_ASOA1, id_ASOA2, id_ASOA3, id_ASOAN
+  INTEGER               :: id_TSOG0, id_TSOG1, id_TSOG2, id_TSOG3
+  INTEGER               :: id_ASOG1, id_ASOG2, id_ASOG3
+  INTEGER               :: id_NIT, id_SO4s, id_NITs, id_HNO3
 #endif
   INTEGER               :: id_SALAAL, id_SALCAL, id_SO4, id_SALC ! MSL
   LOGICAL               :: ok_OH, ok_HO2, ok_O1D, ok_O3P
@@ -213,10 +220,15 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     REAL(f4)               :: TROPv_NOx_mass(State_Grid%NX,State_Grid%NY)
     REAL(dp)               :: localC(NSPEC)
 #endif
-#ifdef MODEL_WRF
+#if defined( MODEL_WRF ) || defined( MODEL_CESM )
     REAL(dp)               :: localC(NSPEC)
 #endif
 
+#if defined( MODEL_CESM )
+    ! Sink rate for artificial UT/LS sink
+    REAL(dp)               :: ScaleCESMLossRate
+#endif
+
     ! Grid box integration time diagnostic
     REAL(fp)               :: TimeStart, TimeEnd
 
@@ -281,7 +293,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
                                       State_Diag%KppAutoReducerNVAR   = 0.0_f4
        IF (State_Diag%Archive_KppcNONZERO)  State_Diag%KppcNONZERO    = 0.0_f4
     ENDIF
-    
+
     ! Also zero satellite diagnostic archival arrays
     IF ( State_Diag%Archive_SatDiagnLoss ) State_Diag%SatDiagnLoss    = 0.0_f4
     IF ( State_Diag%Archive_SatDiagnProd ) State_Diag%SatDiagnProd    = 0.0_f4
@@ -397,7 +409,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
        CALL Timer_Start( "=> FlexChem",           RC ) ! ended in Do_Chemistry
     ENDIF
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
     ! Init diagnostics
     IF ( ASSOCIATED(State_Diag%KppError) ) THEN
        State_Diag%KppError(:,:,:) = 0.0
@@ -513,7 +525,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
     !$OMP PRIVATE( NOxTau,     NOxConc, localC                              )&
     !$OMP PRIVATE( NOx_weight, NOx_tau_weighted                             )&
 #endif
-#ifdef MODEL_WRF
+#if defined( MODEL_WRF ) || defined( MODEL_CESM )
     !$OMP PRIVATE( localC                                                   )&
 #endif
     !$OMP COLLAPSE( 3                                                       )&
@@ -727,6 +739,52 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
           ENDIF
        ENDIF
 
+#if defined( MODEL_CESM )
+       !=====================================================================
+       ! Unphysical fix: Photolyze soluble aerosol tracers
+       ! This removes unphysical values of soluble tracers in the UT/LS due
+       ! to decoupling of convection and wet scavenging in CESM dynamics.
+       !
+       ! This process has to be done before InChemGrid as it is supposed to
+       ! be active everywhere, especially the stratosphere.
+       ! (hplin, 5/30/23)
+       !=====================================================================
+
+       IF ( Input_Opt%correctConvUTLS .and. L .ge. State_Met%PBL_TOP_L(I,J) ) THEN
+
+           ! We operate directly on [molec/cm3] species concentrations in State_Chm,
+           ! because they have not been copied to C() in KPP yet. But, we can use
+           ! PHOTOL(11) which is J-NO2, and scale to create the artificial sink.
+           ! This is a consistent handling based off the MOZART-TS1 mechanism
+           ! in Emmons et al., 2020 JAMES.
+           ScaleCESMLossRate = MAX(0.0_dp, 1 - PHOTOL(11) * .0004_dp * DT)
+
+           State_Chm%Species(id_TSOA0)%Conc(I,J,L) = State_Chm%Species(id_TSOA0)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOA1)%Conc(I,J,L) = State_Chm%Species(id_TSOA1)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOA2)%Conc(I,J,L) = State_Chm%Species(id_TSOA2)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOA3)%Conc(I,J,L) = State_Chm%Species(id_TSOA3)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOA1)%Conc(I,J,L) = State_Chm%Species(id_ASOA1)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOA2)%Conc(I,J,L) = State_Chm%Species(id_ASOA2)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOA3)%Conc(I,J,L) = State_Chm%Species(id_ASOA3)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOAN)%Conc(I,J,L) = State_Chm%Species(id_ASOAN)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOG0)%Conc(I,J,L) = State_Chm%Species(id_TSOG0)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOG1)%Conc(I,J,L) = State_Chm%Species(id_TSOG1)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOG2)%Conc(I,J,L) = State_Chm%Species(id_TSOG2)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_TSOG3)%Conc(I,J,L) = State_Chm%Species(id_TSOG3)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOG1)%Conc(I,J,L) = State_Chm%Species(id_ASOG1)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOG2)%Conc(I,J,L) = State_Chm%Species(id_ASOG2)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_ASOG3)%Conc(I,J,L) = State_Chm%Species(id_ASOG3)%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_NIT  )%Conc(I,J,L) = State_Chm%Species(id_NIT  )%Conc(I,J,L) * ScaleCESMLossRate
+           State_Chm%Species(id_NITs )%Conc(I,J,L) = State_Chm%Species(id_NITs )%Conc(I,J,L) * ScaleCESMLossRate
+
+           ! Don't apply this to sulfate as it is not applied in CAM-chem either and will affect the SO4 budget.
+           !State_Chm%Species(id_SO4  )%Conc(I,J,L) = State_Chm%Species(id_SO4  )%Conc(I,J,L) * ScaleCESMLossRate
+           !State_Chm%Species(id_SO4s )%Conc(I,J,L) = State_Chm%Species(id_SO4s )%Conc(I,J,L) * ScaleCESMLossRate
+
+        ENDIF
+
+#endif
+
        !=====================================================================
        ! Test if we need to do the chemistry for box (I,J,L),
        ! otherwise move onto the next box.
@@ -1049,7 +1107,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
           WRITE(6,*) '### INTEGRATE RETURNED ERROR AT: ', I, J, L
        ENDIF
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
        ! Print grid box indices to screen if integrate failed
        IF ( IERR < 0 ) THEN
           WRITE(6,*) '### INTEGRATE RETURNED ERROR AT: ', I, J, L
@@ -1116,7 +1174,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
        !=====================================================================
        IF ( IERR < 0 ) THEN
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
           ! Save a copy of the C vector (GEOS and WRF only)
           localC    = C
 #endif
@@ -1215,7 +1273,7 @@ SUBROUTINE Do_FullChem( Input_Opt,  State_Chm, State_Diag,                 &
           IF ( IERR < 0 ) THEN
              WRITE(6,     '(a   )' ) '## INTEGRATE FAILED TWICE !!! '
              WRITE(ERRMSG,'(a,i3)' ) 'Integrator error code :', IERR
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
              IF ( Input_Opt%KppStop ) THEN
                 CALL ERROR_STOP(ERRMSG, 'INTEGRATE_KPP')
              ! Revert to start values
@@ -2648,6 +2706,29 @@ SUBROUTINE Init_FullChem( Input_Opt, State_Chm, State_Diag, RC )
     id_IDNOO    = Ind_( 'IDNOO'        )
 #endif
 
+#ifdef MODEL_CESM
+    ! hplin
+    id_TSOA0    = Ind_('TSOA0')
+    id_TSOA1    = Ind_('TSOA1')
+    id_TSOA2    = Ind_('TSOA2')
+    id_TSOA3    = Ind_('TSOA3')
+    id_ASOA1    = Ind_('ASOA1')
+    id_ASOA2    = Ind_('ASOA2')
+    id_ASOA3    = Ind_('ASOA3')
+    id_ASOAN    = Ind_('ASOAN')
+    id_TSOG0    = Ind_('TSOG0')
+    id_TSOG1    = Ind_('TSOG1')
+    id_TSOG2    = Ind_('TSOG2')
+    id_TSOG3    = Ind_('TSOG3')
+    id_ASOG1    = Ind_('ASOG1')
+    id_ASOG2    = Ind_('ASOG2')
+    id_ASOG3    = Ind_('ASOG3')
+    id_NIT      = Ind_('NIT')
+    id_SO4s     = Ind_('SO4s')
+    id_NITs     = Ind_('NITs')
+    id_HNO3     = Ind_('HNO3')
+#endif
+
     ! Set flags to denote if each species is defined
     ok_HO2      = ( id_HO2 > 0         )
     ok_O1D      = ( id_O1D > 0         )

GeosCore/ucx_mod.F90

@@ -3901,9 +3901,6 @@ SUBROUTINE NOXCOEFF_INIT( Input_Opt, State_Grid, State_Chm )
     USE State_Chm_Mod,      ONLY : ChmState
 #if defined( MODEL_CESM )
     USE UNITS,              ONLY : freeUnit
-#if defined( SPMD )
-    USE MPISHORTHAND
-#endif
 #endif
 !
 ! !INPUT PARAMETERS:
@@ -3929,9 +3926,6 @@ SUBROUTINE NOXCOEFF_INIT( Input_Opt, State_Grid, State_Chm )
     INTEGER            :: I, AS, IOS
     INTEGER            :: IMON, ITRAC, ILEV
     INTEGER            :: IU_FILE
-#if defined( MODEL_CESM ) && defined( SPMD )
-    INTEGER            :: nSize ! Number of elements in State_Chm%NOXCOEFF
-#endif
 
     ! Strings
     CHARACTER(LEN=255) :: NOX_FILE
@@ -4055,7 +4049,6 @@ SUBROUTINE NOXCOEFF_INIT( Input_Opt, State_Grid, State_Chm )
     State_Chm%NOXCOEFF = 0.0e+0_fp
 
 #if defined( MODEL_CESM )
-    nSize = State_Chm%JJNOXCOEFF * UCX_NLEVS * 6 * 12
     IF ( Input_Opt%amIRoot ) THEN
 #endif
     ! Fill array
@@ -4136,9 +4129,6 @@ SUBROUTINE NOXCOEFF_INIT( Input_Opt, State_Grid, State_Chm )
     ENDDO !IMON
 #if defined( MODEL_CESM )
     ENDIF
-#if defined( SPMD )
-    CALL MPIBCAST( State_Chm%NOXCOEFF, nSize, MPIR8, 0, MPICOM )
-#endif
 #endif
 
   END SUBROUTINE NOXCOEFF_INIT

Headers/input_opt_mod.F90

@@ -422,15 +422,21 @@ MODULE Input_Opt_Mod
      LOGICAL                     :: TurnOffHetRates
 #endif
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
      LOGICAL                     :: KppStop            = .TRUE. ! Stop KPP if integration fails twice
 #endif
 
 #if defined( MODEL_CESM )
-     LOGICAL                     :: onlineAlbedo       = .TRUE. ! Use albedo from land model
-     LOGICAL                     :: onlineLandTypes    = .TRUE. ! Use land types from land model
-     LOGICAL                     :: ddVel_CLM          = .TRUE. ! Use dry deposition velocities as computed by the Community Land Model
-     LOGICAL                     :: applyQtend         = .TRUE. ! Apply water vapor tendency to specific humidity
+     ! Use albedo from land model
+     LOGICAL                     :: onlineAlbedo       = .TRUE.
+     ! Use land types from land model
+     LOGICAL                     :: onlineLandTypes    = .TRUE.
+     ! Use dry deposition velocities as computed by the Community Land Model
+     LOGICAL                     :: ddVel_CLM          = .TRUE.
+     ! Apply water vapor tendency to specific humidity
+     LOGICAL                     :: applyQtend         = .TRUE.
+     ! Apply photolytic correction for convective scavenging of soluble tracers?
+     LOGICAL                     :: correctConvUTLS    = .TRUE.
 #endif
 
 #ifdef ADJOINT

Headers/state_diag_mod.F90

@@ -1181,9 +1181,6 @@ MODULE State_Diag_Mod
 
      !%%%%% Chemistry diagnostics %%%%%
 
-     REAL(f4),           POINTER :: KppError(:,:,:)
-     LOGICAL                     :: Archive_KppError
-
      REAL(f4),           POINTER :: O3concAfterChem(:,:,:)
      LOGICAL                     :: Archive_O3concAfterChem
 
@@ -1234,10 +1231,10 @@ MODULE State_Diag_Mod
      LOGICAL                     :: Archive_CO2photrate
 #endif
 
-#ifdef MODEL_WRF
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
      !----------------------------------------------------------------------
      ! The following diagnostics are only used when
-     ! GEOS-Chem is interfaced into WRF (as WRF-GC)
+     ! GEOS-Chem is interfaced into WRF (as WRF-GC) or CESM
      !----------------------------------------------------------------------
      REAL(f4),           POINTER :: KppError(:,:,:)
      LOGICAL                     :: Archive_KppError
@@ -2385,10 +2382,10 @@ SUBROUTINE Zero_State_Diag( State_Diag, RC )
     State_Diag%Archive_CO2photrate                 = .FALSE.
 #endif
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
     !=======================================================================
     ! These diagnostics are only activated when running GC
-    ! either in NASA/GEOS or in WRF
+    ! either in NASA/GEOS, WRF, or CESM
     !=======================================================================
     State_Diag%KppError                            => NULL()
     State_Diag%Archive_KppError                    = .FALSE.
@@ -6051,7 +6048,7 @@ SUBROUTINE Init_State_Diag( Input_Opt, State_Chm,       State_Grid,        &
           RETURN
        ENDIF
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
        !--------------------------------------------------------------------
        ! KPP error flag
        !--------------------------------------------------------------------
@@ -12041,10 +12038,10 @@ SUBROUTINE Cleanup_State_Diag( State_Diag, RC )
     IF ( RC /= GC_SUCCESS ) RETURN
 #endif
 
-#if defined(MODEL_GEOS) || defined(MODEL_WRF)
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
     !=======================================================================
     ! These fields are only used when GEOS-Chem
-    ! is interfaced to NASA/GEOS or to WRF (as WRF-GC)
+    ! is interfaced to NASA/GEOS, WRF (as WRF-GC), or CESM
     !=======================================================================
     CALL Finalize( diagId   = 'KppError',                                    &
                    Ptr2Data = State_Diag%KppError,                           &
@@ -12955,7 +12952,7 @@ SUBROUTINE Get_Metadata_State_Diag( am_I_Root,  metadataID, Found,         &
        IF ( isUnits   ) Units = 'kg m-2 s-1'
        IF ( isRank    ) Rank  = 2
 
-#if defined( MODEL_GEOS ) || defined( MODEL_WRF )
+#if defined( MODEL_GEOS ) || defined( MODEL_WRF ) || defined( MODEL_CESM )
     ELSE IF ( TRIM( Name_AllCaps ) == 'KPPERROR' ) THEN
        IF ( isDesc    ) Desc  = 'KppError'
        IF ( isUnits   ) Units = '1'