Revised Lopez lightning flash rate, and made it the default

CHANGELOG.md

@@ -13,9 +13,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - The file path was changed for anthropogenic CO emissions that are used by achem. Note that the previous version of the emissions have an incorrect seasonal cycle.
 - Update ESMF CMake target to `ESMF::ESMF`
+- Overhauled the Lopez lightning scheme, and made it the default scheme; note that lightning is used by GMI for computing NOx emissions; PCHEM does not use lightning
 
 ### Fixed
+
 - Updated GAAS_Gridcomp_Extdata.yaml in AMIP/ to avoid the model to crash when GAAS is turned on and AMIP emissions chosen.
+
 ### Deprecated
 
 

ChemEnv.rc

@@ -1,7 +1,7 @@
 #-----------------------
 # Settings for Lightning
 #-----------------------
-flashSource: MOIST       # MOIST (default), FIT, HEMCO, LOPEZ
+flashSource: LOPEZ       # MOIST, FIT, HEMCO, LOPEZ (default)
 
 # for FIT only:
 ratioGlobalFile: ExtData/g5chem/x/lightning/RatioGlobal.asc
@@ -20,9 +20,11 @@ MOIST_flashFactor_resvec_REPLAY: 2.0       2.0       2.0       1.84      1.84
   FIT_flashFactor_resvec_REPLAY: 1.0       1.0       1.0       1.0       1.0       1.0 
 
 #  originally called 'alpha'
-LOPEZ_flashFactor_resvec_CTM:    37.5      37.5      37.5      37.5    165300      37.5
-LOPEZ_flashFactor_resvec_FREE:   37.5      37.5      37.5      37.5    165300      37.5
-LOPEZ_flashFactor_resvec_REPLAY: 37.5      37.5      37.5      37.5    165300      37.5
+#  10.29.24 Manyin provided c90 and c180, FREE and REPLAY
+#           Assume CTM is close to REPLAY
+LOPEZ_flashFactor_resvec_CTM:    6705.50   6705.50   6705.50   6941.09   6941.09   6941.09
+LOPEZ_flashFactor_resvec_FREE:   5672.30   5672.30   5672.30   6495.07   6495.07   6495.07
+LOPEZ_flashFactor_resvec_REPLAY: 6705.50   6705.50   6705.50   6941.09   6941.09   6941.09
 
 #  originally called 'otdLisScale'
 #  NOTE: testing (Nov2020) indicated c-90 value of 9.48e-3 might be better

GEOS_ChemEnvGridComp.F90

@@ -251,7 +251,7 @@ subroutine SetServices ( GC, RC )
 
     call MAPL_AddImportSpec(GC,                                         &
          SHORT_NAME = 'PFI_CN',                                         &
-         LONG_NAME  = 'ice_convective_precipitation',                   &
+         LONG_NAME  = '3D_flux_of_ice_convective_precipitation',                   &
          UNITS      = 'kg m-2 s-1',                                     &
          DIMS       = MAPL_DimsHorzVert,                                &
          VLOCATION  = MAPL_VLocationEdge,                        __RC__ )
@@ -395,6 +395,24 @@ subroutine SetServices ( GC, RC )
          VLOCATION  = MAPL_VLocationNone,                         __RC__)
 
 
+    call MAPL_AddImportSpec(GC,                               &
+         SHORT_NAME='ZLCL',                                        &
+         LONG_NAME ='lifting_condensation_level',                  &
+         UNITS     ='m'  ,                                         &
+         DIMS      = MAPL_DimsHorzOnly,                            &
+         VLOCATION = MAPL_VLocationNone,                RC=STATUS  )
+    VERIFY_(STATUS)
+
+    call MAPL_AddImportSpec(GC,                               &
+         SHORT_NAME='ZLFC',                                        &
+         LONG_NAME ='level_of_free_convection',                    &
+         UNITS     ='m'  ,                                         &
+         DIMS      = MAPL_DimsHorzOnly,                            &
+         VLOCATION = MAPL_VLocationNone,                RC=STATUS  )
+    VERIFY_(STATUS)
+
+
+
 ! !EXPORT STATE:
 
 !    AIRDENS: Provided for Children
@@ -763,10 +781,10 @@ subroutine Initialize_ ( GC, impChem, expChem, clock, RC )
                                  __RC__ )
 
     IF(MAPL_AM_I_ROOT()) THEN
-      if ( flash_source_enum == FLASH_SOURCE_MOIST ) PRINT*,'MOIST_flashFactor is ',MOIST_flashFactor
-      if ( flash_source_enum == FLASH_SOURCE_FIT   ) PRINT*,'  FIT_flashFactor is ',  FIT_flashFactor
-      if ( flash_source_enum == FLASH_SOURCE_HEMCO ) PRINT*,'HEMCO_flashFactor is ',HEMCO_flashFactor
-      if ( flash_source_enum == FLASH_SOURCE_LOPEZ ) PRINT*,'LOPEZ_flashFactor is ',LOPEZ_flashFactor
+      if ( flash_source_enum == FLASH_SOURCE_MOIST ) PRINT*,'MOIST_flashFactor is ', MOIST_flashFactor
+      if ( flash_source_enum == FLASH_SOURCE_FIT   ) PRINT*,'  FIT_flashFactor is ',   FIT_flashFactor
+      if ( flash_source_enum == FLASH_SOURCE_HEMCO ) PRINT*,'HEMCO_flashFactor is ', HEMCO_flashFactor
+      if ( flash_source_enum == FLASH_SOURCE_LOPEZ ) PRINT*,'LOPEZ_flashFactor is ', LOPEZ_flashFactor
 
                                                      PRINT*,'usePreconCape = ', usePreconCape
     ENDIF
@@ -825,6 +843,9 @@ subroutine Run1 ( GC, IMPORT, EXPORT, CLOCK, RC )
   real, pointer, dimension(:,:)   ::         LWI => null()
   real, pointer, dimension(:,:)   ::        PBLH => null()
 
+  real, pointer, dimension(:,:)   ::          ZLFC => null()
+  real, pointer, dimension(:,:)   ::          ZLCL => null()
+
   real, pointer, dimension(:,:)   ::          TS => null()
   real, pointer, dimension(:,:)   ::     FROCEAN => null()
   real, pointer, dimension(:,:)   ::      FRLAND => null()
@@ -936,11 +957,11 @@ subroutine Run1 ( GC, IMPORT, EXPORT, CLOCK, RC )
       call MAPL_GetPointer ( IMPORT, MIDLAT_ADJ,  'MIDLAT_ADJ',   __RC__ )
     end if
 
-    call MAPL_GetPointer ( IMPORT, CNV_MFC,     'CNV_MFC', __RC__ )
-    call MAPL_GetPointer ( IMPORT, CNV_MFD,     'CNV_MFD', __RC__ )
-    call MAPL_GetPointer ( IMPORT, CN_PRCP,     'CN_PRCP', __RC__ )
-    call MAPL_GetPointer ( IMPORT, PHIS,        'PHIS',    __RC__ )
-    call MAPL_GetPointer ( IMPORT, PFI_CN,      'PFI_CN', ALLOC=.TRUE., __RC__  )     !  ??
+    call MAPL_GetPointer ( IMPORT, CNV_MFC,     'CNV_MFC',  __RC__ )
+    call MAPL_GetPointer ( IMPORT, CNV_MFD,     'CNV_MFD',  __RC__ )
+    call MAPL_GetPointer ( IMPORT, CN_PRCP,     'CN_PRCP',  __RC__ )
+    call MAPL_GetPointer ( IMPORT, PHIS,        'PHIS',     __RC__ )
+    call MAPL_GetPointer ( IMPORT, PFI_CN,      'PFI_CN',   __RC__ )
 
     call MAPL_GetPointer ( IMPORT, T,           'T',        __RC__ )
     call MAPL_GetPointer ( IMPORT, TH,          'TH',       __RC__ )
@@ -958,6 +979,10 @@ subroutine Run1 ( GC, IMPORT, EXPORT, CLOCK, RC )
     call MAPL_GetPointer ( IMPORT,  INHB_PRECON, 'INHB',    __RC__ )
     call MAPL_GetPointer ( IMPORT, BYNCY_PRECON, 'BYNCY',   __RC__ )
 
+    call MAPL_GetPointer ( IMPORT, ZLFC, 'ZLFC',   __RC__ )
+    call MAPL_GetPointer ( IMPORT, ZLCL, 'ZLCL',   __RC__ )
+
+
               BYNCY(:,:,:) = real(0)
                CAPE(:,:)   = real(0)
                 LFR(:,:)   = real(0)
@@ -986,7 +1011,7 @@ subroutine Run1 ( GC, IMPORT, EXPORT, CLOCK, RC )
                 minDeepCloudTop, lightNOampFactor, numberNOperFlash, &
                 MOIST_flashFactor, FIT_flashFactor, HEMCO_flashFactor, LOPEZ_flashFactor, &
                 CNV_MFD, usePreconCape, CAPE_PRECON, INHB_PRECON, BYNCY_PRECON, &
-                CAPE, BYNCY, LFR, LIGHT_NO_PROD, PHIS, &
+                CAPE, BYNCY, ZLFC, ZLCL, LFR, LIGHT_NO_PROD, PHIS, &
                 __RC__)
 
 !   call pmaxmin('LFR', LFR, flashRateMin, flashRateMax, nc, 1, 1.)

Shared/Chem_Shared/Lightning_mod.F90

@@ -133,7 +133,7 @@ subroutine getLightning (GC, ggState, CLOCK, &
      minDeepCloudTop, lightNOampFactor, numberNOperFlash, &
      MOIST_flashFactor, FIT_flashFactor, HEMCO_flashFactor, LOPEZ_flashFactor, &
      CNV_MFD, usePreconCape, CAPE_PRECON, INHB_PRECON, BYNCY_PRECON, &
-     CAPE, BYNCY, flashRate, light_NO_prod, PHIS,  &
+     CAPE, BYNCY, ZLFC, ZLCL, flashRate, light_NO_prod, PHIS,  &
      RC) 
 
   type(ESMF_GridComp),           intent(inout)    :: GC     ! Gridded component
@@ -178,7 +178,7 @@ subroutine getLightning (GC, ggState, CLOCK, &
   real, dimension(:,:,:),   intent(in)  :: PLE      ! edge pressures (Pa)
   real, dimension(:,:,:),   intent(in)  :: ZLE      ! geopotential height (m)             [top-down]
 
-  real, dimension(:,:,:),   intent(in)  :: PFICU    ! flux of ice convective precip (kg m-2 s-1)
+  real, dimension(:,:,:),   intent(in)  :: PFICU    ! flux of ice in convective updrafts (kg m-2 s-1)
   real, dimension(:,:,:),   intent(in)  :: T        ! air temperature (K)                 [top-down]
   real, dimension(:,:,:),   intent(in)  :: TH       ! potential temperature (K)           [top-down]
   real, dimension(:,:,:),   intent(in)  :: Q        ! specific humidity (kg kg-1)         [top-down]
@@ -189,14 +189,13 @@ subroutine getLightning (GC, ggState, CLOCK, &
   real, dimension(:,:,:),   intent(in)  :: CNV_MFD            ! detraining_mass_flux  (kg m-2 s-1)     [top-down]
 
   logical,                  intent(in)  :: usePreconCape   ! Whether to use CAPE, INHB and BYNCY from MOIST
-                                                           ! (only affects LOPEZ and MOIST)
-                                                           ! LOPEZ only uses CAPE, and should benefit from PRECON
-                                                           ! MOIST prefers CAPE computed using the old MERRA2 approach
-                                                           !       which needs BYNCY, so may benefit from PRECON
   real, dimension(:,:),     intent(in)  ::  CAPE_PRECON
   real, dimension(:,:),     intent(in)  ::  INHB_PRECON
   real, dimension(:,:,:),   intent(in)  :: BYNCY_PRECON
 
+  REAL, dimension(:,:),     INTENT(IN)  :: ZLFC            ! Level of Free Convection height [m]
+  REAL, dimension(:,:),     INTENT(IN)  :: ZLCL            ! Lifted Condensation Level height [m]
+
   real, dimension(:,:),     intent(out)  :: CAPE           ! can remove these eventually
   real, dimension(:,:,:),   intent(out)  :: BYNCY          ! can remove these eventually               [top-down]
   real, dimension(:,:),     intent(out)  :: flashRate      ! (flashes/km2/s)
@@ -356,7 +355,6 @@ subroutine getLightning (GC, ggState, CLOCK, &
 
      ! callLopezCalcuations() put above end subroutine getLightning
 
-     ! MEM- Recommend PRECON, but leaving the other option available
      if (usePreconCape) then
        CAPE  =  CAPE_PRECON
        INHB  =  INHB_PRECON
@@ -369,8 +367,8 @@ subroutine getLightning (GC, ggState, CLOCK, &
      flashRateLopez = real(0)
 !    flashRateLopez(:,:) = 0.01
 
-      call LOPEZ_FlashRate(IM, JM, LM, FRLAND, PBLH, CAPE, ZLE2, PFICU, &
-           CNV_QC, T, PLO, LOPEZ_flashFactor, flashRateLopez)
+     call LOPEZ_FlashRate(ggState, IM, JM, LM, FROCEAN, PBLH, CAPE, ZLE2, PFICU, &
+           CNV_QC, T, PLO, ZLFC, ZLCL, LOPEZ_flashFactor, flashRateLopez, __RC__ )
 
 !print*, "min/max LFR LOPEZ: ", minval(flashRateLopez), maxval(flashRateLopez)
 
@@ -419,10 +417,6 @@ subroutine getLightning (GC, ggState, CLOCK, &
 
      DM(:,:,1:LM) = ( PLE(:,:,K0+1:KM)-PLE(:,:,K0:KM-1) ) * (1./MAPL_GRAV)  ! DELP / g
 
-     ! (Don't use CAPE_PRECON since the units are different from CAPE_MERRA2)
-     ! The MOIST approach wants CAPE computed in the old MERRA2 way
-     ! which uses BYNCY, so usePreconCape really means usePreconBuoyancy in this case
-     ! Probably best to use PRECON
      if (usePreconCape) then
 
        CAPE  =  CAPE_PRECON
@@ -615,7 +609,7 @@ subroutine read_flash_source ( rcfilen, flash_source_enum, RC )
   ! How to calculate flashrate
   call ESMF_ConfigGetAttribute( esmfConfig, flashSource,    &
                                 Label    = "flashSource:",  &
-                                Default  = 'MOIST',  __RC__ )
+                                Default  = 'LOPEZ',  __RC__ )
 
   call identify_flash_source( flashSource, flash_source_enum )
   IF ( flash_source_enum == FLASH_SOURCE_UNDEFINED ) THEN
@@ -1534,44 +1528,76 @@ end subroutine HEMCO_FlashRate
 !-----------------------------------------------------------------------
 
 
-    subroutine LOPEZ_FlashRate(IM, JM, LM, FRLAND, ZKBCON, CAPE,  &
-                               ZLE, PFI_CN_GF, CNV_QC, TE, PLO, LOPEZ_flashFactor, &
-                               LFR_Lopez)
+    subroutine LOPEZ_FlashRate(STATE, IM, JM, LM, FROCEAN, ZKBCON, CAPE,  &
+                               ZLE, PFI_CN_GF, CNV_QC, TE, PLO, ZLFC, ZLCL, LOPEZ_flashFactor, &
+                               LFR_Lopez, RC)
      implicit none 
+     TYPE(MAPL_MetaComp), POINTER :: STATE ! Internal MAPL_Generic state
      integer ,intent(in)  :: im,jm,lm
-     real    ,intent(in), dimension(im,jm,0:lm) ::  &
-                              ZLE        & ! layer depths [m]
-                             ,PFI_CN_GF       ! 3d_ice_precipitation_flux_GF [kg/m2/s]
+     real,    intent(in), dimension(im,jm,0:lm) ::  ZLE          ! layer depths [m]
+     real,    intent(in), dimension(im,jm,0:lm) ::  PFI_CN_GF    ! 3d_ice_precipitation_flux_GF_CONVECTIVE [kg/m2/s]  in updrafts
+
+     real,    intent(in), dimension(im,jm,lm) ::  TE             ! air temp [K]
+     real,    intent(in), dimension(im,jm,lm) ::  PLO            ! press (hPa)
+     real,    intent(in), dimension(im,jm,lm) ::  CNV_QC         ! grid_mean_convective_condensate [kg/kg]
 
-     real    ,intent(in), dimension(im,jm,lm) ::  &
-                             TE             &  ! air temp [K]
-                             ,PLO           &  ! press (hPa)
-                             ,CNV_QC        ! grid_mean_convective_condensate [kg/kg]
 
+     real,    intent(in), dimension(im,jm) :: FROCEAN    ! fraction of ocean
+     real,    intent(in), dimension(im,jm) :: cape       ! Convective available potential energy [J/kg]
+     real,    intent(in), dimension(im,jm) :: zkbcon     ! cloud_base_height_deep_GF [m]
 
-     real    ,intent(in), dimension(im,jm) :: FRLAND    & ! fraction of land
-          , cape   &  ! Convective available potential energy [J/kg]
-          , zkbcon    ! cloud_base_height_deep_GF [m]
+     real,    intent(in), dimension(im,jm) :: ZLFC       ! Level of Free Convection height [m]
+     real,    intent(in), dimension(im,jm) :: ZLCL       ! Lifted Condensation Level height [m]
 
-     real    ,intent(in)                   :: LOPEZ_flashFactor   ! global scaling term
+     real,    intent(in)                   :: LOPEZ_flashFactor   ! global scaling term
                                                                   ! originally called 'alpha'
 
-     real    ,intent(out), dimension(im,jm) ::  &
-                              LFR_Lopez           ! lightning flash density rate (units: flashes/km2/day)
+     real,    intent(out), dimension(im,jm):: LFR_Lopez           ! lightning flash density rate (units: flashes/km2/day)
 
+     integer, intent(out), optional        :: RC
+
     !-- locals
+
+    integer :: STATUS
+    character(len=*), parameter :: Iam = "LOPEZ_FlashRate"
+
     real, parameter :: V_graup     = 3.0  ! m/s
     real, parameter :: V_snow      = 0.5  ! m/s
-    real, parameter :: beta_land   = 0.70 ! 1
-    real, parameter :: beta_ocean  = 0.45 ! 1
+!   real, parameter :: beta_land   = 0.70 ! 1
+!   real, parameter :: beta_ocean  = 0.45 ! 1
+!   real, parameter :: beta_land   = 0.55 ! 1
+!   real, parameter :: beta_ocean  = 0.15 ! 1
     real, parameter :: t_initial   =  0.0 + 273.15 ! K
     real, parameter :: t_final     = -25. + 273.15 ! K
+
+!   REAL :: a_lgt         ! Parameters for Lightning Flash Rate 
+!   REAL :: v_graup
+!   REAL :: v_snow
+    REAL :: b_land
+    REAL :: b_marine
+!   REAL :: cb_exp
+!   REAL :: cb_max
+!   REAL :: cb_ref
+!   REAL :: t_min_ec
+    REAL :: f_oc_thresh
+    REAL :: cb_min
+!   REAL :: c_lfc_land
+!   REAL :: c_lfc_marine
 
     integer :: i,j,k, k_initial, k_final                          
     real    :: tdif, td2
     real    :: Q_R, z_base,beta,prec_flx_fr,dz
     real,    dimension(1:lm) :: q_graup,q_snow,rho
 
+!   REAL :: f_oc_thresh
+!   REAL :: cb_min
+    REAL :: c_lfc_land
+    REAL :: c_lfc_marine
+
+!   f_oc_thresh  = 0.50
+!   cb_min       = 0.4
+    c_lfc_land   = 0.4
+    c_lfc_marine = 0.1
 
 !!! DEBUGGING
 !    print*, ""
@@ -1589,24 +1615,59 @@ subroutine LOPEZ_FlashRate(IM, JM, LM, FRLAND, ZKBCON, CAPE,  &
 !    print*, "zkbcon: ", minval(zkbcon), maxval(zkbcon)
 !!!
 
+    if (PRESENT(RC)) RC = ESMF_SUCCESS
+
+!   ! Parameters for Lightning Flash Rate from Lopez et al.
+!   ! ------------------------------------------------
+!   CALL MAPL_GetResource(STATE,       a_lgt,'A_LGT:',         DEFAULT= 10000.0,   __RC__ )
+!   CALL MAPL_GetResource(STATE,     v_graup,'V_GRAUP:',       DEFAULT= 3.0,       __RC__ )
+!   CALL MAPL_GetResource(STATE,      v_snow,'V_SNOW:',        DEFAULT= 0.5,       __RC__ )
+    CALL MAPL_GetResource(STATE,      b_land,'B_LAND:',        DEFAULT= 0.7,       __RC__ )
+    CALL MAPL_GetResource(STATE,    b_marine,'B_MARINE:',      DEFAULT= 0.45,      __RC__ )
+!   CALL MAPL_GetResource(STATE,      cb_exp,'CB_EXP:',        DEFAULT= 2.0,       __RC__ )
+!   CALL MAPL_GetResource(STATE,      cb_max,'CB_MAX:',        DEFAULT= 1.8,       __RC__ )
+!   CALL MAPL_GetResource(STATE,      cb_ref,'CB_REF:',        DEFAULT= 1.0,       __RC__ )
+!   CALL MAPL_GetResource(STATE,    t_min_ec,'T_MIN_EC:',      DEFAULT= 248.15,    __RC__ )
+    CALL MAPL_GetResource(STATE, f_oc_thresh,'F_OC_THRESH:',   DEFAULT= 0.50,      __RC__ )
+    CALL MAPL_GetResource(STATE,      cb_min,'CB_MIN:',        DEFAULT= 0.35,      __RC__ )  ! modified from MINDS, for better land/ocean ratio
+!   CALL MAPL_GetResource(STATE,  c_lfc_land,'C_LFC_LAND:',    DEFAULT= 0.4,       __RC__ )
+!   CALL MAPL_GetResource(STATE,c_lfc_marine,'C_LFC_MARINE:',  DEFAULT= 0.1,       __RC__ )
+
+
     DO j = 1, jm
         DO i = 1, im
            LFR_Lopez(i,j) = 0.0
 
-           if(ZKBCON(i,j) <= 0. .or. CAPE(i,j) == MAPL_UNDEF) cycle !-> no convection 
 
-           beta= frland(i,j)*beta_land + (1.-frland(i,j))*beta_ocean
+           !!!  z_base - as computed in MINDS:
+           IF (FROCEAN(i,j) >= f_oc_thresh) THEN
+             z_base = c_lfc_marine * ZLFC(i,j) + (1.0 - c_lfc_marine) * ZLCL(i,j)
+           ELSE
+             z_base = c_lfc_land   * ZLFC(i,j) + (1.0 - c_lfc_land  ) * ZLCL(i,j)
+           END IF
+
+           z_base = z_base * 0.001 - cb_min
+           z_base = MAX(z_base, 0.0)
+
+
+!          beta= frland(i,j)*b_land + (1.-frland(i,j))*b_marine
+           beta= frocean(i,j)*b_marine + (1.-frocean(i,j))*b_land
+
            q_graup(:) = 0. 
            q_snow (:) = 0. 
 
-           do k=lm,1,-1
+           do k=1,lm
                 rho(k) =  100.*PLO(i,j,k) / (MAPL_RGAS*TE(i,j,k) )
-                prec_flx_fr = PFI_CN_GF(i,j,k) / rho(k)
+                prec_flx_fr = (( PFI_CN_GF(i,j,k  )                      +    &
+                                 PFI_CN_GF(i,j,k-1)                      ) * 0.5 ) / rho(k)
 
                 q_graup(k) =      beta *prec_flx_fr/V_graup ! - graupel mixing ratio (kg/kg)
                 q_snow (k) =  (1.-beta)*prec_flx_fr/V_snow  ! - snow    mixing ratio (kg/kg)
 
            enddo
+
+           if(z_base      <= 0. .or. CAPE(i,j) == MAPL_UNDEF) cycle !-> no convection 
+
            k_initial = minloc(abs(te(i,j,:)-t_initial),1)  
            k_final   = minloc(abs(te(i,j,:)-t_final  ),1)  
 
@@ -1641,7 +1702,6 @@ subroutine LOPEZ_FlashRate(IM, JM, LM, FRLAND, ZKBCON, CAPE,  &
              Q_R = Q_R + dz*rho(k)*(q_graup(k)*(CNV_QC(i,j,k)+q_snow(k)))
          enddo
 
-         z_base = ZKBCON(i,j)/1000. !- convert to [km]
 
         !--- lightning flash density (units: number of flashes/km2/day) - equation 5
         !--- (to compare with Lopez 2016's results, convert to per year: LFR_Lopez*365)