Merge pull request #2523 from yuanjianz/bugfix/Grell_Friedas_convection

Fix convection issues when meteorology is Grell Freitas scheme

CHANGELOG.md

@@ -15,6 +15,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 - Added KPP standalone interface (archives model state to selected locations)
 - Added `https://github/geoschem/KPP-Standalone` as a Git submodule
 - Added comments in `./run/sharedcleanRunDir.sh` describing the `--force` option (i.e. remove files w/o user confirmation)
+- Specified meteorology source in GCHP geoschem_config.yml
+- Added Input_Opt logical for whether to reconstruct convective precipitation fluxes rather than use met-fields
+- Added to run directory creation a warning about convection discontinuity and bug if GEOS-FP meteorology is chosen
+- Added surface precipitation flux fields as inputs to GCHP
 
 ### Changed
 - Renamed `Emiss_Carbon_Gases` to `CO2_Production` in `carbon_gases_mod.F90`
@@ -32,6 +36,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 - Fixed the satellite diagnostics counters from being inadvertently being reset
 - Fixed segmentation fault in qfyaml when running with certain compilers without debug flags on
 - Fixed errors in adjoint-only code preventing successful adjoint build
+- Fixed zero convective precipitation and high cloud base in runs using GEOS-FP (>=01Jun2020) or GEOS-IT
 
 ### Removed
 - Removed duplicate `WD_RetFactor` tag for HgClHO2 in `species_database.yml`

GeosCore/convection_mod.F90

@@ -342,7 +342,7 @@ SUBROUTINE DO_CONVECTION( Input_Opt,  State_Chm, State_Diag, &
     !$OMP END PARALLEL DO
 #endif
 
-    ! Return if COMPUTE_F returned an error
+    ! Return if Do_Cloud_Convection returned an error
     IF ( RC /= GC_SUCCESS ) THEN
        ErrMsg = 'Error encountered in "Do_Cloud_Convection"!'
        CALL GC_Error( ErrMsg, RC, ThisLoc )
@@ -526,21 +526,25 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     ! Arrays
     REAL(fp)               :: BMASS    (State_Grid%NZ)
     REAL(fp)               :: PDOWN    (State_Grid%NZ)
+    REAL(fp)               :: REEVAPCN (State_Grid%NZ)
+    REAL(fp)               :: DQRCU    (State_Grid%NZ)
 
     ! Pointers
-    REAL(fp),      POINTER :: BXHEIGHT (:)
-    REAL(fp),      POINTER :: CMFMC    (:)
-    REAL(fp),      POINTER :: DQRCU    (:)
-    REAL(fp),      POINTER :: DTRAIN   (:)
-    REAL(fp),      POINTER :: PFICU    (:)
-    REAL(fp),      POINTER :: PFLCU    (:)
-    REAL(fp),      POINTER :: REEVAPCN (:)
-    REAL(fp),      POINTER :: DELP_DRY (:)
-    REAL(fp),      POINTER :: T        (:)
-    REAL(fp),      POINTER :: H2O2s    (:)
-    REAL(fp),      POINTER :: SO2s     (:)
-    REAL(fp),      POINTER :: Q        (:)
-    TYPE(SpcConc), POINTER :: Spc      (:)
+    REAL(fp),      POINTER :: BXHEIGHT     (:)
+    REAL(fp),      POINTER :: CMFMC        (:)
+    REAL(fp),      POINTER :: DQRCU_MET    (:)
+    REAL(fp),      POINTER :: DTRAIN       (:)
+    REAL(fp),      POINTER :: PFICU        (:)
+    REAL(fp),      POINTER :: PFLCU        (:)
+    REAL(fp),      POINTER :: REEVAPCN_MET (:)
+    REAL(fp),      POINTER :: DELP_DRY     (:)
+    REAL(fp),      POINTER :: DELP         (:)
+    REAL(fp),      POINTER :: T            (:)
+    REAL(fp),      POINTER :: H2O2s        (:)
+    REAL(fp),      POINTER :: SO2s         (:)
+    REAL(fp),      POINTER :: Q            (:)
+    REAL(fp),      POINTER :: PRECCON
+    TYPE(SpcConc), POINTER :: Spc          (:)
     TYPE(Species), POINTER :: SpcInfo
 
     !========================================================================
@@ -553,18 +557,20 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     ThisLoc  = ' -> at Do_Cloud_Convection (in convection_mod.F90)'
 
     ! Point to columns of derived-type object fields
-    BXHEIGHT => State_Met%BXHEIGHT(I,J,:        ) ! Box height [m]
-    CMFMC    => State_Met%CMFMC   (I,J,2:State_Grid%NZ+1) ! Cloud mass flux
-                                                          ! [kg/m2/s]
-    DQRCU    => State_Met%DQRCU   (I,J,:        ) ! Precip production rate:
-    DTRAIN   => State_Met%DTRAIN  (I,J,:        ) ! Detrainment flux [kg/m2/s]
-    REEVAPCN => State_Met%REEVAPCN(I,J,:        ) ! Evap of precip'ing conv.
-    DELP_DRY => State_Met%DELP_DRY(I,J,:        ) ! Edge dry P diff [hPa]
-    T        => State_Met%T       (I,J,:        ) ! Air temperature [K]
-    H2O2s    => State_Chm%H2O2AfterChem(I,J,:   ) ! H2O2s from sulfate_mod
-    SO2s     => State_Chm%SO2AfterChem (I,J,:   ) ! SO2s from sulfate_mod
-    Spc      => State_Chm%Species              ! Chemical species vector
-    SpcInfo  => NULL()                            ! Species database entry
+    BXHEIGHT     => State_Met%BXHEIGHT(I,J,:        ) ! Box height [m]
+    CMFMC        => State_Met%CMFMC   (I,J,2:State_Grid%NZ+1) ! Cloud mass flux
+                                                              ! [kg/m2/s] [upper edge]
+    DQRCU_MET    => State_Met%DQRCU   (I,J,:        ) ! Precip production rate:
+    DTRAIN       => State_Met%DTRAIN  (I,J,:        ) ! Detrainment flux [kg/m2/s]
+    REEVAPCN_MET => State_Met%REEVAPCN(I,J,:        ) ! Evap of precip'ing conv.
+    DELP_DRY     => State_Met%DELP_DRY(I,J,:        ) ! Edge dry P diff [hPa]
+    DELP         => State_Met%DELP    (I,J,:        ) ! Edge P diff [hPa]
+    T            => State_Met%T       (I,J,:        ) ! Air temperature [K]
+    PRECCON      => State_Met%PRECCON (I,J          ) ! Surface precipitation flux [mm/day]
+    H2O2s        => State_Chm%H2O2AfterChem(I,J,:   ) ! H2O2s from sulfate_mod
+    SO2s         => State_Chm%SO2AfterChem (I,J,:   ) ! SO2s from sulfate_mod
+    Spc          => State_Chm%Species              ! Chemical species vector
+    SpcInfo      => NULL()                            ! Species database entry
 
     ! PFICU and PFLCU are on level edges
     PFICU    => State_Met%PFICU   (I,J,2:State_Grid%NZ+1) ! Dwnwd flx of conv
@@ -584,6 +590,9 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     ! Convection timestep [s]
     NDT      = TS_DYN
 
+    ! Is this a Hg simulation?
+    IS_Hg = Input_Opt%ITS_A_MERCURY_SIM
+
     IF ( State_Grid%NZ > 72 .or. &
          Input_Opt%MetField == "MODELE2.1" ) THEN 
        ! Higher vertical resolution runs need shorter convective timestep
@@ -597,6 +606,31 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     DNS      = DBLE( NS )               ! Num internal timesteps (real)
     SDT      = DBLE( NDT ) / DBLE( NS ) ! seconds in internal timestep
 
+    IF ( .NOT. Input_Opt%Reconstruct_Conv_Precip_Flux ) THEN
+       ! RAS scheme
+       REEVAPCN(:) = REEVAPCN_MET(:)
+       DQRCU(:) = DQRCU_MET(:)
+    ELSE
+       ! GF scheme
+       ! REEVAPCN_MET is cumulative re-evaporation
+       ! DQRCU_MET is net precipitation formation (need to add re-evaporation back)
+       REEVAPCN(NLAY) = REEVAPCN_MET(NLAY)
+       DO K = 2, NLAY-1
+         ! REEVAPCN (kg/kg/s) to REEVAPCN_FLUX (kg/m2/s)
+         ! subtraction between fluxes instead
+          REEVAPCN(K) = (REEVAPCN_MET(K) * DELP(K) &
+                     - REEVAPCN_MET(K+1) * DELP(K+1) ) &
+                     / DELP(K)
+       ENDDO
+       ! Zero DQRCU and some negative REEVAPCN at surface in GF scheme
+       ! Set as zero to avoid cloud base mistakenly at surface level
+       REEVAPCN(1) = 0.0_fp
+       DQRCU(1) = 0.0_fp
+       DO K = 2, NLAY
+          DQRCU(K) = DQRCU_MET(K) + REEVAPCN(K)
+       ENDDO
+    ENDIF
+
     !-----------------------------------------------------------------
     ! Determine location of the cloud base, which is the level where
     ! we start to have non-zero convective precipitation formation
@@ -616,24 +650,29 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     !-----------------------------------------------------------------
     ! Compute PDOWN and BMASS
     !-----------------------------------------------------------------
-    DO K = 1, NLAY
-
-       ! PDOWN is the convective precipitation leaving each
-       ! box [cm3 H2O/cm2 air/s]. This accounts for the
-       ! contribution from both liquid & ice precip.
-       ! PFLCU and PFICU are converted from kg/m2/s to m3/m2/s
-       ! using water and ice densities, respectively.
-       ! m3/m2/s becomes cm3/cm2/s using a factor of 100.
-       PDOWN(K) = ( ( PFLCU(K) / 1000e+0_fp ) &
-                +   ( PFICU(K) /  917e+0_fp ) ) * 100e+0_fp
-
-       ! BMASS is the dry air mass per unit area for the grid box
-       ! bounded by level K and K+1 [kg/m2]
-       ! BMASS is equivalent to deltaP (dry) * 100 / g
-       ! This is done to keep BMASS in the same units as CMFMC * SDT
-       BMASS(K) = DELP_DRY(K) * G0_100
+    ! PDOWN is the convective precipitation leaving each
+    ! box [cm3 H2O/cm2 air/s]. This accounts for the
+    ! contribution from both liquid & ice precip.
+    ! PFLCU and PFICU are converted from kg/m2/s to m3/m2/s
+    ! using water and ice densities, respectively.
+    ! m3/m2/s becomes cm3/cm2/s using a factor of 100.
+    IF ( .NOT. Input_Opt%Reconstruct_Conv_Precip_Flux ) THEN
+       PDOWN(:) = ( ( PFLCU(:) / 1000e+0_fp ) &
+                +   ( PFICU(:) /  917e+0_fp ) ) * 100e+0_fp
+    ELSE
+       PDOWN(NLAY) = DQRCU_MET(NLAY) * DELP(NLAY) * G0_100 * 100e+0_fp
+       DO K = NLAY-1, 1, -1
+          PDOWN(K) = PDOWN(K+1) + DQRCU_MET(K) & 
+                    * DELP(K) * G0_100 * 100e+0_fp
+                    ! kg/kg wet air/s to kg/m2/s then to cm3 H20/cm2 air/s
+       ENDDO
+    ENDIF
 
-    ENDDO
+    ! BMASS is the dry air mass per unit area for the grid box
+    ! bounded by level K and K+1 [kg/m2]
+    ! BMASS is equivalent to deltaP (dry) * 100 / g
+    ! This is done to keep BMASS in the same units as CMFMC * SDT
+    BMASS(:) = DELP_DRY(:) * G0_100
 
     !-----------------------------------------------------------------
     ! Compute MB, the mass per unit area of dry air below the cloud
@@ -645,9 +684,6 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
        MB = MB + BMASS(K)
     ENDDO
 
-    ! Is this a Hg simulation?
-    IS_Hg = Input_Opt%ITS_A_MERCURY_SIM
-
     !========================================================================
     ! (1)  A d v e c t e d   S p e c i e s   L o o p
     !========================================================================
@@ -1374,18 +1410,19 @@ SUBROUTINE DO_CLOUD_CONVECTION( Input_Opt,  &
     !================================================================
 
     ! Nullify pointers
-    NULLIFY( BXHEIGHT )
-    NULLIFY( CMFMC    )
-    NULLIFY( DQRCU    )
-    NULLIFY( DTRAIN   )
-    NULLIFY( PFICU    )
-    NULLIFY( PFLCU    )
-    NULLIFY( REEVAPCN )
-    NULLIFY( DELP_DRY )
-    NULLIFY( T        )
-    NULLIFY( H2O2s    )
-    NULLIFY( SO2s     )
-    NULLIFY( Spc      )
+    NULLIFY( BXHEIGHT     )
+    NULLIFY( CMFMC        )
+    NULLIFY( DQRCU_MET    )
+    NULLIFY( DTRAIN       )
+    NULLIFY( PFICU        )
+    NULLIFY( PFLCU        )
+    NULLIFY( REEVAPCN_MET )
+    NULLIFY( DELP_DRY     )
+    NULLIFY( DELP         )
+    NULLIFY( T            )
+    NULLIFY( H2O2s        )
+    NULLIFY( SO2s         )
+    NULLIFY( Spc          )
 
     ! Set error code to success
     RC                      = GC_SUCCESS

GeosCore/input_mod.F90

@@ -3202,7 +3202,9 @@ END SUBROUTINE Config_Photolysis
 ! !IROUTINE: config_convection_mixing
 !
 ! !DESCRIPTION: Copies convection & PBL mixing information from the Config
-!  object to Input_Opt, and does necessary checks.
+!  object to Input_Opt, and does necessary checks. Also sets whether
+!  to reconstruct convective precipitation flux based on meteorology
+!  source and simulation start date.
 !\\
 !\\
 ! !INTERFACE:
@@ -3286,13 +3288,40 @@ SUBROUTINE Config_Convection_Mixing( Config, Input_Opt, RC )
     ENDIF
     Input_Opt%LNLPBL = v_bool
 
+    !------------------------------------------------------------------------
+    ! Other settings based on inputs
+    !------------------------------------------------------------------------
+
     ! Set the PBL drydep flag. This determines if dry deposition is
     ! applied (and drydep frequencies are calculated) over the entire
     ! PBL or the first model layer only. For now, set this value
     ! automatically based upon the selected PBL scheme: 1st model layer
     ! for the non-local PBL scheme, full PBL for the full-mixing scheme.
     Input_Opt%PBL_DRYDEP = ( .not. Input_Opt%LNLPBL )
 
+    ! Set whether to reconstruct convective precipitation flux based on
+    ! meteorology source and simulation start date. This is important for
+    ! avoiding a bug where convective precipitation met-fields are all zero
+    ! in GEOS-IT for all years and in GEOS-FP following June 1, 2020.
+    !
+    ! IMPORTANT NOTE: The logic for GEOS-FP assumes (1) meteorology year
+    ! is the same as simulation year and (2) the simulation does not
+    ! run across June 1, 2020. Use the following rules to ensure your
+    ! simulation is correct:
+    !   (1) Manually update code below if GEOS-FP data year is
+    !       different than simulation year:
+    !          - Set to .FALSE. if data is prior to June 1, 2020
+    !          - Set to .TRUE. if data is on or after June 1, 2020
+    !   (2) Do not run a GEOS-FP simulation across June 1, 2020. Split
+    !       up the run in time to avoid this.
+    IF ( Input_Opt%MetField == 'GEOSIT' ) THEN
+       Input_Opt%Reconstruct_Conv_Precip_Flux = .TRUE.
+    ELSEIF ( Input_Opt%MetField == 'GEOSFP' .AND. Input_Opt%NYMDb >= 20200601 ) THEN
+       Input_Opt%Reconstruct_Conv_Precip_Flux = .TRUE.
+    ELSE
+       Input_Opt%Reconstruct_Conv_Precip_Flux = .FALSE.
+    ENDIF
+
     ! Return success
     RC = GC_SUCCESS
 
@@ -3303,6 +3332,16 @@ SUBROUTINE Config_Convection_Mixing( Config, Input_Opt, RC )
        WRITE( 6, 90  ) 'CONVECTION SETTINGS'
        WRITE( 6, 95  ) '-------------------'
        WRITE( 6, 100 ) 'Turn on cloud convection?   : ', Input_Opt%LCONV
+       WRITE( 6, 100 ) 'Reconstruct convective precipitation flux?   : ', &
+            Input_Opt%Reconstruct_Conv_Precip_Flux
+
+       IF ( Input_Opt%MetField == 'GEOSFP' ) THEN
+          IF ( Input_Opt%Reconstruct_Conv_Precip_Flux ) THEN
+             WRITE( 6, 90 ) 'WARNING: Convection will assume met data is on or after 01Jun2020!'
+          ELSE
+             WRITE( 6, 90 ) 'WARNING: Convection will assume met data is prior to 01Jun2020!'
+          ENDIF
+       ENDIF
 
        WRITE( 6, 90  ) 'PBL MIXING SETTINGS'
        WRITE( 6, 95  ) '-------------------'

Headers/input_opt_mod.F90

@@ -239,6 +239,7 @@ MODULE Input_Opt_Mod
      LOGICAL                     :: LCONV
      LOGICAL                     :: LTURB
      LOGICAL                     :: LNLPBL
+     LOGICAL                     :: Reconstruct_Conv_Precip_Flux
      INTEGER                     :: TS_CONV
 
      !----------------------------------------
@@ -759,6 +760,7 @@ SUBROUTINE Set_Input_Opt( am_I_Root, Input_Opt, RC )
     Input_Opt%LCONV                  = .FALSE.
     Input_Opt%LTURB                  = .FALSE.
     Input_Opt%LNLPBL                 = .FALSE.
+    Input_Opt%Reconstruct_Conv_Precip_Flux = .FALSE.
     Input_Opt%TS_CONV                = 0
 
     !----------------------------------------

Interfaces/GCHP/Includes_Before_Run.H

@@ -154,6 +154,10 @@
   ! transport is correctly implemented.
   State_Met%SPHU1               = SPHU1       *1.0d3         ! kg/kg -> g/kg
   State_Met%SPHU2               = SPHU2       *1.0d3         ! kg/kg -> g/kg
+  State_Met%PRECCON             = PRECCON * 86400d0          ! kg/m2/s to mm/day
+  State_Met%PRECLSC             = PRECLSC * 86400d0          ! kg/m2/s to mm/day
+  State_Met%PRECANV             = PRECANV * 86400d0          ! kg/m2/s to mm/day
+  State_Met%PRECTOT             = PRECTOT * 86400d0          ! kg/m2/s to mm/day
 
     ! If meteorology vertical index is top down, flip imports coming from ExtData.
   if (meteorology_vertical_index_is_top_down) then

Interfaces/GCHP/Registry/Chem_Registry.rc

@@ -25,7 +25,7 @@
 #  Name 	     |   Units     | Dim |Loc|Type| R | A |Tiles| ault |	Name
 # -------------------|-------------|-----|---|----|---|---|-----|------|--------------------------
 # AIRDENS            |  kg m-3     | xyz | C |    |   |   |     |      | air_density
-  AREA               |     m2      | xy  |   |    |   |   |	|      | 
+  AREA               |     m2    | xy  |   |    |   |   |	|      | 
   FRLAKE             |     1	   | xy  |   |    |   |   |	|      | 
   FRLAND             |     1	   | xy  |   |    |   |   |	|      | 
   FRLANDIC           |     1	   | xy  |   |    |   |   |	|      | 
@@ -96,9 +96,13 @@
   TMPU1              |     1	   | xyz | C |    |   |   |	|      | 
   SPHU2              |     1	   | xyz | C |    |   |   |	|      | 
   TMPU2              |     1	   | xyz | C |    |   |   |	|      | 
-  FLASH_DENS         |     1       | xy  |   |    |   |   |     |      | lightning_flash_density
-  CONV_DEPTH         |     1       | xy  |   |    |   |   |     |      | convective_cloud_depth
-  XLAIMULTI          | cm2_cm-2    | xyz | E |    |   |   |	|      | LAI_by_type 
+  FLASH_DENS         |     1     | xy  |   |    |   |   | |      | lightning_flash_density
+  CONV_DEPTH         |     1     | xy  |   |    |   |   | |      | convective_cloud_depth
+  XLAIMULTI          | cm2_cm-2  | xyz | E |    |   |   |	|      | LAI_by_type 
+  PRECCON            |     1	   | xy  |   |    |   |   |	|      | surface convective precipitation flux
+  PRECLSC            |     1	   | xy  |   |    |   |   |	|      | surface large-scale precipitation flux
+  PRECANV            |     1	   | xy  |   |    |   |   |	|      | surface anvil precipitation flux
+  PRECTOT            |     1	   | xy  |   |    |   |   |	|      | surface total precipitation flux
 # -------------------|-------------|-----|---|----|---|---|-----|------|--------------------------
 </ImportSpec>