Physics (PBL-SATMEDF) Update for Explicit Subcanopy Effects on Chemical Tracers

ccpp/config/ccpp_prebuild_config_fv3.py

@@ -189,6 +189,7 @@
     'physics/physics/PBL/MYNN_EDMF/mynnedmf_wrapper.F90',
     'physics/physics/PBL/SATMEDMF/satmedmfvdif.F',
     'physics/physics/PBL/SATMEDMF/satmedmfvdifq.F',
+    'physics/physics/PBL/SATMEDMF/canopy_driver.F',
     'physics/physics/PBL/YSU/ysuvdif.F90',
     'physics/physics/PBL/saYSU/shinhongvdif.F90',
     'physics/physics/Radiation/RRTMG/radsw_main.F90',

ccpp/data/GFS_typedefs.F90

@@ -40,7 +40,7 @@ module GFS_typedefs
    integer, parameter :: dfi_radar_max_intervals = 4 !< Number of radar-derived temperature tendency and/or convection suppression intervals. Do not change.
 
    real(kind=kind_phys), parameter :: limit_unspecified = 1e12 !< special constant for "namelist value was not provided" in radar-derived temperature tendency limit range
-   
+
    integer, parameter :: physics_no_tracer = -99
 
 !> \section arg_table_GFS_typedefs
@@ -687,7 +687,7 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: dqdt_qmicro(:,:) => null()  !< instantanious microphysics tendency to be passed from MP to convection
 
     !-- lake surface temperature from cdeps inline
-    real (kind=kind_phys), pointer :: mask_dat   (:) => null()   !< land-sea mask from cdeps inline  
+    real (kind=kind_phys), pointer :: mask_dat   (:) => null()   !< land-sea mask from cdeps inline
     real (kind=kind_phys), pointer :: tsfco_dat  (:) => null()   !< sfc temperature from cdeps inline
     real (kind=kind_phys), pointer :: tice_dat   (:) => null()   !< sfc temperature over ice from cdeps inline
     real (kind=kind_phys), pointer :: hice_dat   (:) => null()   !< sfc ice thickness from cdeps inline
@@ -792,7 +792,7 @@ module GFS_typedefs
     integer              :: dycore_active   !< Choice of dynamical core
     integer              :: dycore_fv3  = 1 !< Choice of FV3 dynamical core
     integer              :: dycore_mpas = 2 !< Choice of MPAS dynamical core
-    
+
 !--- coupling parameters
     logical              :: cplflx          !< default no cplflx collection
     logical              :: cplice          !< default no cplice collection (used together with cplflx)
@@ -1086,7 +1086,6 @@ module GFS_typedefs
     real(kind=kind_phys) :: fs_fac_rain     !< adjustment for rain fall speed
     real(kind=kind_phys) :: fs_fac_snow     !< adjustment for snow fall speed
 
-
     !--- GFDL microphysical paramters
     logical              :: lgfdlmprad      !< flag for GFDL mp scheme and radiation consistency
     logical              :: phys_hydrostatic
@@ -1573,6 +1572,11 @@ module GFS_typedefs
     integer              :: nchem           !< number of prognostic chemical species (vertically mixied)
     integer              :: ndvel           !< number of prognostic chemical species (which are deposited, usually =nchem)
     integer              :: ntchm           !< number of prognostic chemical tracers (advected)
+! "cplaqm" tracers
+    integer              :: nto3            !< tracer index for Ozone chemical species CMAQ
+    integer              :: ntno            !< tracer index for NO    chemical species CMAQ
+    integer              :: ntno2           !< tracer index for NO2   chemical species CMAQ
+
     integer              :: ntchs           !< tracer index for first prognostic chemical tracer
     integer              :: ntche           !< tracer index for last prognostic chemical tracer
     integer              :: ntdu1           !< tracer index for dust bin1
@@ -2156,6 +2160,11 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: dkt(:,:)       => null()  !< Eddy diffusitivity for heat
     real (kind=kind_phys), pointer :: dku(:,:)       => null()  !< Eddy diffusitivity for momentum
 
+!3-LAYER CANOPY
+    !--- Extra PBL diagnostics in canopy
+    real (kind=kind_phys), pointer :: dkt_can(:,:)   => null()  !< Eddy diffusitivity for heat
+    real (kind=kind_phys), pointer :: dku_can(:,:)   => null()  !< Eddy diffusitivity for momentum
+
 !
 !---vay-2018 UGWP-diagnostics instantaneous
 !
@@ -2260,7 +2269,6 @@ module GFS_typedefs
     ! Diagnostics for coupled air quality model
     real (kind=kind_phys), pointer :: aod   (:)   => null()    !< instantaneous aerosol optical depth ( n/a )
 
-!IVAI
     ! Diagnostics for coupled air quality model
     real (kind=kind_phys), pointer :: coszens(:)  => null()    ! Cosine SZA for photolysis
     real (kind=kind_phys), pointer :: jo3o1d(:)   => null()    ! instantaneous O3O1D photolysis rate
@@ -2270,7 +2278,6 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: cfrt (:)    => null()    ! Forest Fraction
     real (kind=kind_phys), pointer :: cclu (:)    => null()    ! Clumping Index
     real (kind=kind_phys), pointer :: cpopu(:)    => null()    ! Population density
-!IVAI
 
     ! Auxiliary output arrays for debugging
     real (kind=kind_phys), pointer :: aux2d(:,:)  => null()    !< auxiliary 2d arrays in output (for debugging)
@@ -3432,7 +3439,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
                                  communicator, ntasks, nthreads,    &
                                  tile_num, isc, jsc, nx, ny,  cnx,  &
                                  cny, gnx, gny, ak, bk, hydrostatic)
-    
+
 !--- modules
     use physcons,         only: con_rerth, con_pi
     use mersenne_twister, only: random_setseed, random_number
@@ -3613,7 +3620,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
     logical              :: lrseeds           = .false.      !< flag to use host-provided random seeds
     integer              :: nrstreams         = 2            !< number of random number streams in host-provided random seed array
     logical              :: lextop            = .false.      !< flag for using an extra top layer for radiation
-    real(kind_phys)      :: xr_con            = -999.0       !< Xu-Randall cloud fraction multiplicative constant          
+    real(kind_phys)      :: xr_con            = -999.0       !< Xu-Randall cloud fraction multiplicative constant
     real(kind_phys)      :: xr_exp            = -999.0       !< Xu-Randall cloud fraction exponent constant
     ! RRTMGP
     logical              :: do_RRTMGP           = .false.    !< Use RRTMGP?
@@ -4437,7 +4444,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
           stop
        endif
     endif
-    
+
     ! dtend selection: default is to match all variables:
     dtend_select(1)='*'
     do ipat=2,pat_count
@@ -5575,6 +5582,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
     Model%dtidx = physics_no_tracer
 
     if(Model%ntchm>0) then
+! GFS_v16 n=9 "no2" n=10 "no" n=11 "o3" (n=8,9, 10 in PBL resp.)
+      Model%ntno2 = get_physics_tracer_index('no2', Model)  ! n=11 (index 10 "no2" in PBL scheme) GFS_v17_p8
+      Model%ntno  = get_physics_tracer_index('no', Model)   ! n=12 (index 11 "no"  in PBL scheme) GFS_v17_p8
+      Model%nto3  = get_physics_tracer_index('o3', Model)   ! n=13 (index 12 "o3"  in PBL scheme) GFS_v17_p8
+
       Model%ntdu1 = get_physics_tracer_index('dust1', Model)
       Model%ntdu2 = get_physics_tracer_index('dust2', Model)
       Model%ntdu3 = get_physics_tracer_index('dust3', Model)
@@ -5649,7 +5661,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
            endif
 
            ! More specific chemical tracer names:
+! NB. ntchs is 1st chemical tracer (not so2 tracer)
            call label_dtend_tracer(Model,100+Model%ntchs,'so2','sulfur dioxide concentration','kg kg-1 s-1')
+
            if(Model%ntchm>0) then
               ! Need better descriptions of these.
               call label_dtend_tracer(Model,100+Model%ntchm+Model%ntchs-1,'pp10','pp10 concentration','kg kg-1 s-1')
@@ -5658,14 +5672,14 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
               if(itrac>0) then
                  call label_dtend_tracer(Model,100+itrac,'DMS','DMS concentration','kg kg-1 s-1')
               endif
+
               itrac=get_physics_tracer_index('msa', Model)
               if(itrac>0) then
                  call label_dtend_tracer(Model,100+itrac,'msa','msa concentration','kg kg-1 s-1')
               endif
            endif
         endif
 
-
         call label_dtend_tracer(Model,Model%index_of_temperature,'temp','temperature','K s-1')
         call label_dtend_tracer(Model,Model%index_of_x_wind,'u','x wind','m s-2')
         call label_dtend_tracer(Model,Model%index_of_y_wind,'v','y wind','m s-2')
@@ -5698,6 +5712,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
         call label_dtend_tracer(Model,100+Model%ntia,'ice_aero','number concentration of ice-friendly aerosols','kg-1 s-1')
         call label_dtend_tracer(Model,100+Model%nto,'o_ion','oxygen ion concentration','kg kg-1 s-1')
         call label_dtend_tracer(Model,100+Model%nto2,'o2','oxygen concentration','kg kg-1 s-1')
+! cplaqm tracers CMAQ
+        call label_dtend_tracer(Model,100+Model%ntno2,'no2_cpl','cplaqm NO2   concentration','kg kg-1 s-1')
+        call label_dtend_tracer(Model,100+Model%ntno, 'no_cpl', 'cplaqm NO    concentration','kg kg-1 s-1')
+        call label_dtend_tracer(Model,100+Model%nto3, 'o3_cpl', 'cplaqm ozone concentration','kg kg-1 s-1')
+
         call label_dtend_cause(Model,Model%index_of_process_pbl,'pbl','tendency due to PBL')
         call label_dtend_cause(Model,Model%index_of_process_dcnv,'deepcnv','tendency due to deep convection')
         call label_dtend_cause(Model,Model%index_of_process_scnv,'shalcnv','tendency due to shallow convection')
@@ -5780,6 +5799,12 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
              enddo
           endif
 
+! NB. In PBL scheme chemical tracers indexes are offset by 1
+! (qdiag3d) cplaqm tracers "no2", "no", "o3"
+          call fill_dtidx(Model,dtend_select,100+Model%ntno2,Model%index_of_process_pbl,have_pbl) ! ntno2= 11  (index 10 is "no2" in PBL scheme) GFS_v17_p8
+          call fill_dtidx(Model,dtend_select,100+Model%ntno ,Model%index_of_process_pbl,have_pbl) ! ntno = 12  (index 11 is "no"  in PBL scheme) GFS_v17_p8
+          call fill_dtidx(Model,dtend_select,100+Model%nto3 ,Model%index_of_process_pbl,have_pbl) ! nto3 = 13  (index 12 is "o3"  in PBL scheme) GFS_v17_p8
+
           call fill_dtidx(Model,dtend_select,100+Model%ntoz,Model%index_of_process_pbl,have_pbl)
           call fill_dtidx(Model,dtend_select,100+Model%ntoz,Model%index_of_process_prod_loss,have_oz_phys)
           call fill_dtidx(Model,dtend_select,100+Model%ntoz,Model%index_of_process_ozmix,have_oz_phys)
@@ -6603,7 +6628,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
 !--- BEGIN CODE FROM GLOOPB
 !--- set up random number seed needed for RAS and old SAS and when cal_pre=.true.
 !    Model%imfdeepcnv < 0 when Model%ras = .true.
-    
+
     if (xr_con > 0.0 .and. xr_exp > 0.0) then !values have been read in from namelist, so set them to read values
       Model%xr_con = xr_con
       Model%xr_exp = xr_exp
@@ -6628,9 +6653,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me,         &
           Model%xr_con = 2000.0
           Model%xr_exp = 0.25
         endif
-      endif     
+      endif
     endif
-    
+
     if (Model%imfdeepcnv <= 0 .or. Model%cal_pre ) then
       if (Model%random_clds) then
         seed0 = Model%idate(1) + Model%idate(2) + Model%idate(3) + Model%idate(4)
@@ -7361,6 +7386,9 @@ subroutine control_print(Model)
       print *, ' nqrimef           : ', Model%nqrimef
       print *, ' ntqv              : ', Model%ntqv
       print *, ' ntoz              : ', Model%ntoz
+      print *, ' ntno2             : ', Model%ntno2 ! "no2"  tracer cplaqm/CMAQ
+      print *, ' ntno              : ', Model%ntno  ! "no"   tracer cplaqm/CMAQ
+      print *, ' nto3              : ', Model%nto3  ! "o3"   tracer cplaqm/CMAQ
       print *, ' ntcw              : ', Model%ntcw
       print *, ' ntiw              : ', Model%ntiw
       print *, ' ntrw              : ', Model%ntrw
@@ -7775,7 +7803,7 @@ subroutine cldprop_create (Cldprop, Model)
     Cldprop%cvt = clear_val
     Cldprop%cvb = clear_val
     Cldprop%cnvw = clear_val
-    
+
   end subroutine cldprop_create
 
 
@@ -8297,6 +8325,12 @@ subroutine diag_create (Diag, Model)
     allocate (Diag%dkt(IM,Model%levs))
     allocate (Diag%dku(IM,Model%levs))
 
+    !--- New PBL Diagnostics in 3-layer canopy
+    if (Model%do_canopy .and. Model%cplaqm) then
+      allocate (Diag%dkt_can(IM,Model%levs))
+      allocate (Diag%dku_can(IM,Model%levs))
+    endif
+
     !--  New max hourly diag.
     allocate (Diag%refdmax(IM))
     allocate (Diag%refdmax263k(IM))
@@ -8373,12 +8407,11 @@ subroutine diag_create (Diag, Model)
       Diag%aod = zero
     end if
 
-!IVAI:
     ! Air quality diagnostics
     ! -- initialize diagnostic variables
     if (Model%cplaqm) then
 
-!IVAI: photdiag arrays
+! photdiag arrays
       allocate (Diag%coszens(IM))
       Diag%coszens= zero
 
@@ -8388,7 +8421,7 @@ subroutine diag_create (Diag, Model)
       allocate (Diag%jno2(IM))
       Diag%jno2 = zero
 
-!IVAI: canopy arrays read via aqm_emis_read
+! Canopy arrays read via aqm_emis_read
       if (Model%do_canopy) then
         allocate (Diag%claie(IM))
         Diag%claie = zero
@@ -8407,7 +8440,6 @@ subroutine diag_create (Diag, Model)
       end if! (Model%do_canopy)
 
     end if ! (Model%cplaqm)
-!IVAI
 
     ! Auxiliary arrays in output for debugging
     if (Model%naux2d>0) then
@@ -8675,6 +8707,12 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     Diag%dkt = zero
     Diag%dku = zero
 
+! Extra PBL diagnostics in 3-layer canopy
+    if (Model%do_canopy .and. Model%cplaqm ) then
+      Diag%dkt_can = zero
+      Diag%dku_can = zero
+    endif
+
 ! max hourly diagnostics
     Diag%refl_10cm   = -35.
     Diag%max_hail_diam_sfc = -999.
@@ -8708,22 +8746,22 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     endif
 
   end subroutine diag_phys_zero
-  
+
   function get_physics_tracer_index (name, Model)
     !This function uses the FMS version of get_tracer_index, but changes the missing tracer index to the value used throughout the physics code, rather than the one used in FMS
     use tracer_manager_mod, only: get_tracer_index, NO_TRACER
     use field_manager_mod, only: MODEL_ATMOS
-    
+
     character(len=*),  intent(in) :: name
     type(GFS_control_type), intent(in) :: Model
-    
+
     !--- local variables
     integer :: get_physics_tracer_index
-    
+
     get_physics_tracer_index = get_tracer_index(MODEL_ATMOS, name, verbose = (Model%me == Model%master) .and. Model%debug)
-    
+
     if (get_physics_tracer_index == NO_TRACER) get_physics_tracer_index = physics_no_tracer
-    
+
   end function get_physics_tracer_index
 
 end module GFS_typedefs

ccpp/data/GFS_typedefs.meta

@@ -6824,6 +6824,29 @@
   units = index
   dimensions = ()
   type = integer
+### IVAI
+[ntno2]
+  standard_name = index_for_no2_chemical_tracer_in_tracer_concentration_array
+  long_name = tracer index for coupled AQM/CMAQ NO2 chemical tracer
+  units = index
+  dimensions = ()
+  type = integer
+##  active = (flag_for_air_quality_coupling .and. flag_for_canopy_option)
+[ntno]
+  standard_name = index_for_no_chemical_tracer_in_tracer_concentration_array
+  long_name = tracer index for coupled AQM/CMAQ NO chemical tracer
+  units = index
+  dimensions = ()
+  type = integer
+##  active = (flag_for_air_quality_coupling .and. flag_for_canopy_option)
+[nto3]
+  standard_name = index_for_ozone_chemical_tracer_in_tracer_concentration_array
+  long_name = tracer index for coupled AQM/CMAQ ozone chemical tracer
+  units = index
+  dimensions = ()
+  type = integer
+##  active = (flag_for_air_quality_coupling .and. flag_for_canopy_option)
+### IVAI
 [ntcw]
   standard_name = index_of_cloud_liquid_water_mixing_ratio_in_tracer_concentration_array
   long_name = tracer index for cloud condensate (or liquid water)
@@ -9916,6 +9939,22 @@
   dimensions = (horizontal_dimension,vertical_layer_dimension)
   type = real
   kind = kind_phys
+[dkt_can]
+  standard_name = atmosphere_heat_diffusivity_in_canopy
+  long_name = atmospheric heat diffusivity in canopy
+  units = m2 s-1
+  dimensions = (horizontal_dimension,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = (flag_for_air_quality_coupling .and. flag_for_canopy_option)
+[dku_can]
+  standard_name = atmosphere_momentum_diffusivity_in_canopy
+  long_name = atmospheric momentum diffusivity in canopy
+  units = m2 s-1
+  dimensions = (horizontal_dimension,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  active = (flag_for_air_quality_coupling .and. flag_for_canopy_option)
 [cldfra]
   standard_name = instantaneous_3d_cloud_fraction
   long_name = instantaneous 3D cloud fraction for all MPs

ccpp/suites/suite_FV3_GFS_v16.xml

@@ -57,6 +57,7 @@
       <scheme>GFS_surface_generic_post</scheme>
       <scheme>GFS_PBL_generic_pre</scheme>
       <scheme>satmedmfvdifq</scheme>
+      <scheme>canopy_driver</scheme>
       <scheme>GFS_PBL_generic_post</scheme>
       <scheme>GFS_GWD_generic_pre</scheme>
       <scheme>cires_ugwp</scheme>

ccpp/suites/suite_FV3_GFS_v17_p8.xml

@@ -53,6 +53,7 @@
       <scheme>GFS_surface_generic_post</scheme>
       <scheme>GFS_PBL_generic_pre</scheme>
       <scheme>satmedmfvdifq</scheme>
+      <scheme>canopy_driver</scheme>
       <scheme>GFS_PBL_generic_post</scheme>
       <scheme>GFS_GWD_generic_pre</scheme>
       <scheme>unified_ugwp</scheme>