Physics update in cumulus convection, PBL & surface layer for UFS_P7

physics/GFS_PBL_generic.F90

@@ -317,7 +317,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
         dqsfc_cpl, dusfci_cpl, dvsfci_cpl, dtsfci_cpl, dqsfci_cpl, dusfc_diag, dvsfc_diag, dtsfc_diag, dqsfc_diag,             &
         dusfci_diag, dvsfci_diag, dtsfci_diag, dqsfci_diag, dt3dt, du3dt_PBL, du3dt_OGWD, dv3dt_PBL, dv3dt_OGWD, dq3dt,        &
         dq3dt_ozone, rd, cp, fvirt, hvap, t1, q1, prsl, hflx, ushfsfci, oceanfrac, kdt, dusfc_cice, dvsfc_cice,                &
-        dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, hffac, hefac, &
+        dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, hffac,                      &
         ugrs, vgrs, tgrs, qgrs, save_u, save_v, save_t, save_q, errmsg, errflg)
 
       use machine,                only : kind_phys
@@ -366,7 +366,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       real(kind=kind_phys), dimension(:), intent(out) ::  ushfsfci
 
       ! From canopy heat storage - reduction factors in latent/sensible heat flux due to surface roughness
-      real(kind=kind_phys), dimension(:), intent(in) :: hffac, hefac
+      real(kind=kind_phys), dimension(:), intent(in) :: hffac
 
       character(len=*), intent(out) :: errmsg
       integer, intent(out) :: errflg
@@ -543,7 +543,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
               else                           !use PBL fluxes when CICE fluxes is unavailable
                 dusfci_cpl(i) = dusfc1(i)
                 dvsfci_cpl(i) = dvsfc1(i)
-                dtsfci_cpl(i) = dtsfc1(i)
+                dtsfci_cpl(i) = dtsfc1(i)*hffac(i)
                 dqsfci_cpl(i) = dqsfc1(i)
               end if
             elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point
@@ -562,7 +562,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
               dusfci_cpl(i) = dusfc1(i)
               dvsfci_cpl(i) = dvsfc1(i)
               dtsfci_cpl(i) = dtsfc1(i)*hffac(i)
-              dqsfci_cpl(i) = dqsfc1(i)*hefac(i)
+              dqsfci_cpl(i) = dqsfc1(i)
             endif
 !
             dusfc_cpl (i) = dusfc_cpl(i) + dusfci_cpl(i) * dtf
@@ -606,7 +606,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
           dusfci_diag(i) = dusfc1(i)
           dvsfci_diag(i) = dvsfc1(i)
           dtsfci_diag(i) = dtsfc1(i)*hffac(i)
-          dqsfci_diag(i) = dqsfc1(i)*hefac(i)
+          dqsfci_diag(i) = dqsfc1(i)
           dtsfc_diag (i) = dtsfc_diag(i) + dtsfci_diag(i) * dtf
           dqsfc_diag (i) = dqsfc_diag(i) + dqsfci_diag(i) * dtf
         enddo

physics/GFS_PBL_generic.meta

@@ -1317,15 +1317,6 @@
   kind = kind_phys
   intent = in
   optional = F
-[hefac]
-  standard_name = surface_upward_latent_heat_flux_reduction_factor
-  long_name = surface upward latent heat flux reduction factor from canopy heat storage
-  units = none
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
-  intent = in
-  optional = F
 [ugrs]
   standard_name = x_wind
   long_name = zonal wind

physics/GFS_debug.F90

@@ -1188,7 +1188,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_ice            ', Interstitial%ep1d_ice                )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_land           ', Interstitial%ep1d_land               )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ep1d_water          ', Interstitial%ep1d_water              )
-                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evapq               ', Interstitial%evapq                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_ice            ', Interstitial%evap_ice                )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_land           ', Interstitial%evap_land               )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%evap_water          ', Interstitial%evap_water              )
@@ -1231,7 +1230,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gflx_water          ', Interstitial%gflx_water              )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gwdcu               ', Interstitial%gwdcu                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%gwdcv               ', Interstitial%gwdcv                   )
-                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hefac               ', Interstitial%hefac                   )
+                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%zvfun               ', Interstitial%zvfun                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hffac               ', Interstitial%hffac                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hflxq               ', Interstitial%hflxq                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%hflx_ice            ', Interstitial%hflx_ice                )

physics/GFS_surface_composites.F90

@@ -419,20 +419,22 @@ end subroutine GFS_surface_composites_post_finalize
 !!
    subroutine GFS_surface_composites_post_run (                                                                                   &
       im, kice, km, rd, rvrdm1, cplflx, cplwav2atm, frac_grid, flag_cice, thsfc_loc, islmsk, dry, wet, icy, wind, t1, q1, prsl1,  &
-      landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli,                                                                   &
+      landfrac, lakefrac, oceanfrac, zorl, zorlo, zorll, zorli, garea,                                                            &
       cd, cd_wat, cd_lnd, cd_ice, cdq, cdq_wat, cdq_lnd, cdq_ice, rb, rb_wat, rb_lnd, rb_ice, stress, stress_wat, stress_lnd,     &
       stress_ice, ffmm, ffmm_wat, ffmm_lnd, ffmm_ice, ffhh, ffhh_wat, ffhh_lnd, ffhh_ice, uustar, uustar_wat, uustar_lnd,         &
       uustar_ice, fm10, fm10_wat, fm10_lnd, fm10_ice, fh2, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice,            &
       cmm, cmm_wat, cmm_lnd, cmm_ice, chh, chh_wat, chh_lnd, chh_ice, gflx, gflx_wat, gflx_lnd, gflx_ice, ep1d, ep1d_wat,         &
       ep1d_lnd, ep1d_ice, weasd, weasd_lnd, weasd_ice, snowd, snowd_lnd, snowd_ice, tprcp, tprcp_wat,                             &
       tprcp_lnd, tprcp_ice, evap, evap_wat, evap_lnd, evap_ice, hflx, hflx_wat, hflx_lnd, hflx_ice, qss, qss_wat, qss_lnd,        &
-      qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice, stc,                 &
+      qss_ice, tsfc, tsfco, tsfcl, tsfc_wat, tsfc_lnd, tsfc_ice, tisfc, tice, hice, cice, min_seaice, tiice,                      &
+      sigmaf, zvfun, lheatstrg, h0facu, h0facs, hflxq, hffac, stc,                                                                &
       grav, prsik1, prslk1, prslki, z1, ztmax_wat, ztmax_lnd, ztmax_ice, errmsg, errflg)
 
       implicit none
 
       integer,                              intent(in) :: im, kice, km
       logical,                              intent(in) :: cplflx, frac_grid, cplwav2atm
+      logical,                              intent(in) :: lheatstrg
       logical, dimension(:),                intent(in) :: flag_cice, dry, wet, icy
       integer, dimension(:),                intent(in) :: islmsk
       real(kind=kind_phys), dimension(:),   intent(in) :: wind, t1, q1, prsl1, landfrac, lakefrac, oceanfrac,                   &
@@ -441,13 +443,15 @@ subroutine GFS_surface_composites_post_run (
         fm10_wat, fm10_lnd, fm10_ice, fh2_wat, fh2_lnd, fh2_ice, tsurf_wat, tsurf_lnd, tsurf_ice, cmm_wat, cmm_lnd, cmm_ice,    &
         chh_wat, chh_lnd, chh_ice, gflx_wat, gflx_lnd, gflx_ice, ep1d_wat, ep1d_lnd, ep1d_ice, weasd_lnd, weasd_ice,            &
         snowd_lnd, snowd_ice, tprcp_wat, tprcp_lnd, tprcp_ice, evap_wat, evap_lnd, evap_ice, hflx_wat, hflx_lnd,                &
-        hflx_ice, qss_wat, qss_lnd, qss_ice, tsfc_wat, tsfc_lnd, tsfc_ice, zorlo, zorll, zorli
+        hflx_ice, qss_wat, qss_lnd, qss_ice, tsfc_wat, tsfc_lnd, tsfc_ice, zorlo, zorll, zorli, garea
 
       real(kind=kind_phys), dimension(:),   intent(inout) :: zorl, cd, cdq, rb, stress, ffmm, ffhh, uustar, fm10,               &
         fh2, cmm, chh, gflx, ep1d, weasd, snowd, tprcp, evap, hflx, qss, tsfc, tsfco, tsfcl, tisfc
 
       real(kind=kind_phys), dimension(:),   intent(in   ) :: tice ! interstitial sea ice temperature
       real(kind=kind_phys), dimension(:),   intent(inout) :: hice, cice
+      real(kind=kind_phys), dimension(:),   intent(inout) :: sigmaf, zvfun, hflxq, hffac                  
+      real(kind=kind_phys),                 intent(in   ) :: h0facu, h0facs
       real(kind=kind_phys),                 intent(in   ) :: min_seaice
       real(kind=kind_phys),                 intent(in   ) :: rd, rvrdm1
 
@@ -468,6 +472,10 @@ subroutine GFS_surface_composites_post_run (
       real(kind=kind_phys) :: txl, txi, txo, wfrac, q0, rho
       ! For calling "stability"
       real(kind=kind_phys) :: tsurf, virtfac, tv1, thv1, tvs, z0max, ztmax
+!
+      real(kind=kind_phys) :: tem1, tem2, gdx
+      real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0
+!
 
       ! Initialize CCPP error handling variables
       errmsg = ''
@@ -490,6 +498,8 @@ subroutine GFS_surface_composites_post_run (
           weasd(i)  = txl*weasd_lnd(i)  + txi*weasd_ice(i)
           snowd(i)  = txl*snowd_lnd(i)  + txi*snowd_ice(i)
          !tprcp(i)  = txl*tprcp_lnd(i)  + txi*tprcp_ice(i)  + txo*tprcp_wat(i)
+!
+          sigmaf(i) = txl*sigmaf(i)
 
           if (.not. flag_cice(i)) then
             if (islmsk(i) == 2) then
@@ -573,8 +583,32 @@ subroutine GFS_surface_composites_post_run (
             stress(i) = stress_ice(i)
             uustar(i) = uustar_ice(i)
           else ! Mix of multiple surface types (land, water, and/or ice)
-            call stability(z1(i), snowd(i), thv1, wind(i), z0max, ztmax, tvs, grav, & ! inputs
-                           tv1, thsfc_loc,                                          & ! inputs
+!
+! re-compute zvfun with composite surface roughness & green vegetation fraction
+!
+            tem1 = (z0max - z0lo) / (z0up - z0lo)
+            tem1 = min(max(tem1, zero), one)
+            tem2 = max(sigmaf(i), 0.1)
+            zvfun(i) = sqrt(tem1 * tem2)
+            gdx = sqrt(garea(i))
+!
+! re-compute variables for canopy heat storage parameterization with the updated zvfun
+!      in the fractional grid
+!
+            hflxq(i) = hflx(i)
+            hffac(i) = 1.0
+            if (lheatstrg) then
+              if(hflx(i) > 0.) then
+                hffac(i) = h0facu * zvfun(i)
+              else
+                hffac(i) = h0facs * zvfun(i)
+              endif
+              hffac(i) = 1. + hffac(i)
+              hflxq(i) = hflx(i) / hffac(i)
+            endif
+!
+            call stability(z1(i), zvfun(i), gdx, tv1, thv1, wind(i),                & ! inputs
+                           z0max, ztmax, tvs, grav, thsfc_loc,                      & ! inputs          
                            rb(i), ffmm(i), ffhh(i), fm10(i), fh2(i), cd(i), cdq(i), & ! outputs
                            stress(i), uustar(i))
           endif ! Checking to see if point is one or multiple surface types

physics/GFS_surface_composites.meta

@@ -1899,6 +1899,77 @@
   kind = kind_phys
   intent = in
   optional = F
+[garea]
+  standard_name = cell_area
+  long_name = area of the grid cell
+  units = m2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[zvfun]
+  standard_name = function_of_surface_roughness_length_and_green_vegetation_fraction
+  long_name = function of surface roughness length and green vegetation fraction
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
+[lheatstrg]
+  standard_name = flag_for_canopy_heat_storage
+  long_name = flag for canopy heat storage parameterization
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+  optional = F
+[h0facu]
+  standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_unstable_surface_layer
+  long_name = canopy heat storage factor for sensible heat flux in unstable surface layer
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[h0facs]
+  standard_name = canopy_heat_storage_factor_for_sensible_heat_flux_in_stable_surface_layer
+  long_name = canopy heat storage factor for sensible heat flux in stable surface layer
+  units = none
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+  optional = F
+[hflxq]
+  standard_name = kinematic_surface_upward_sensible_heat_flux_reduced_by_surface_roughness_and_vegetation
+  long_name = kinematic surface upward sensible heat flux reduced by surface roughness and vegetation
+  units = K m s-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
+[hffac]
+  standard_name = surface_upward_sensible_heat_flux_reduction_factor
+  long_name = surface upward sensible heat flux reduction factor from canopy heat storage
+  units = none
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
+[sigmaf]
+  standard_name = bounded_vegetation_area_fraction
+  long_name = areal fractional cover of green vegetation bounded on the bottom
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+  optional = F
 [ztmax_wat]
   standard_name = bounded_surface_roughness_length_for_heat_over_water
   long_name = bounded surface roughness length for heat over water

physics/GFS_surface_generic.F90

@@ -204,20 +204,21 @@ end subroutine GFS_surface_generic_post_finalize
 !> \section arg_table_GFS_surface_generic_post_run Argument Table
 !! \htmlinclude GFS_surface_generic_post_run.html
 !!
-      subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy, wet, dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,&
+      subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, dry, icy, wet,                                    &
+        dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,                                                                            &
         adjsfcdlw, adjsfcdsw, adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfculw, adjsfculw_wat, adjnirbmu, adjnirdfu,           &
         adjvisbmu, adjvisdfu,t2m, q2m, u10m, v10m, tsfc, tsfc_wat, pgr, xcosz, evbs, evcw, trans, sbsno, snowc, snohf,              &
         epi, gfluxi, t1, q1, u1, v1, dlwsfci_cpl, dswsfci_cpl, dlwsfc_cpl, dswsfc_cpl, dnirbmi_cpl, dnirdfi_cpl, dvisbmi_cpl,       &
         dvisdfi_cpl, dnirbm_cpl, dnirdf_cpl, dvisbm_cpl, dvisdf_cpl, nlwsfci_cpl, nlwsfc_cpl, t2mi_cpl, q2mi_cpl, u10mi_cpl,        &
         v10mi_cpl, tsfci_cpl, psurfi_cpl, nnirbmi_cpl, nnirdfi_cpl, nvisbmi_cpl, nvisdfi_cpl, nswsfci_cpl, nswsfc_cpl, nnirbm_cpl,  &
         nnirdf_cpl, nvisbm_cpl, nvisdf_cpl, gflux, evbsa, evcwa, transa, sbsnoa, snowca, snohfa, ep,                                &
-        runoff, srunoff, runof, drain, lheatstrg, z0fac, e0fac, zorl, hflx, evap, hflxq, evapq, hffac, hefac, errmsg, errflg)
+        runoff, srunoff, runof, drain, lheatstrg, h0facu, h0facs, zvfun, hflx, evap, hflxq, hffac, errmsg, errflg)
 
         implicit none
 
         integer,                                intent(in) :: im
         logical,                                intent(in) :: cplflx, cplchm, cplwav, lssav
-        logical, dimension(:),                  intent(in) :: icy, wet
+        logical, dimension(:),                  intent(in) :: dry, icy, wet
         real(kind=kind_phys),                   intent(in) :: dtf
 
         real(kind=kind_phys), dimension(:),  intent(in)  :: ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1, adjsfcdlw, adjsfcdsw,  &
@@ -235,11 +236,11 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy,
 
         ! For canopy heat storage
         logical, intent(in) :: lheatstrg
-        real(kind=kind_phys), intent(in) :: z0fac, e0fac
-        real(kind=kind_phys), dimension(:), intent(in)  :: zorl
+        real(kind=kind_phys), intent(in) :: h0facu, h0facs
+        real(kind=kind_phys), dimension(:), intent(in)  :: zvfun
         real(kind=kind_phys), dimension(:), intent(in)  :: hflx,  evap
-        real(kind=kind_phys), dimension(:), intent(out) :: hflxq, evapq
-        real(kind=kind_phys), dimension(:), intent(out) :: hffac, hefac
+        real(kind=kind_phys), dimension(:), intent(out) :: hflxq
+        real(kind=kind_phys), dimension(:), intent(out) :: hffac
 
         ! CCPP error handling variables
         character(len=*), intent(out) :: errmsg
@@ -248,13 +249,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy,
         ! Local variables
         real(kind=kind_phys), parameter :: albdf = 0.06_kind_phys
 
-        ! Parameters for canopy heat storage parametrization
-        real(kind=kind_phys), parameter :: z0min=0.2, z0max=1.0
-        real(kind=kind_phys), parameter :: u10min=2.5, u10max=7.5
-
         integer :: i
         real(kind=kind_phys) :: xcosz_loc, ocalnirdf_cpl, ocalnirbm_cpl, ocalvisdf_cpl, ocalvisbm_cpl
-        real(kind=kind_phys) :: tem, tem1, tem2
 
         ! Initialize CCPP error handling variables
         errmsg = ''
@@ -359,32 +355,28 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplchm, cplwav, lssav, icy,
           enddo
         endif
 
-!  --- ...  Boundary Layer and Free atmospheic turbulence parameterization
 !
-!  in order to achieve heat storage within canopy layer, in the canopy heat
-!    storage parameterization the kinematic sensible and latent heat fluxes
-!    (hflx & evap) as surface boundary forcings to the pbl scheme are
-!    reduced as a function of surface roughness
+!  in order to achieve heat storage within canopy layer, in the canopy
+!    heat torage parameterization the kinematic sensible heat flux
+!    (hflx) as surface boundary forcing to the pbl scheme is
+!    reduced in a factor of hffac given as a function of surface roughness &
+!    green vegetation fraction (zvfun) 
 !
         do i=1,im
           hflxq(i) = hflx(i)
-          evapq(i) = evap(i)
-          hffac(i) = one
-          hefac(i) = one
+          hffac(i) = 1.0
         enddo
         if (lheatstrg) then
           do i=1,im
-            tem = 0.01_kind_phys * zorl(i)     ! change unit from cm to m
-            tem1 = (tem - z0min) / (z0max - z0min)
-            hffac(i) = z0fac * min(max(tem1, zero), one)
-            tem  = sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i))
-            tem1 = (tem - u10min) / (u10max - u10min)
-            tem2 = one - min(max(tem1, zero), one)
-            hffac(i) = tem2 * hffac(i)
-            hefac(i) = one + e0fac * hffac(i)
-            hffac(i) = one + hffac(i)
-            hflxq(i) = hflx(i) / hffac(i)
-            evapq(i) = evap(i) / hefac(i)
+            if (dry(i)) then
+              if(hflx(i) > 0.) then
+                hffac(i) = h0facu * zvfun(i)
+              else
+                hffac(i) = h0facs * zvfun(i)
+              endif
+              hffac(i) = 1. + hffac(i)
+              hflxq(i) = hflx(i) / hffac(i)
+            endif
           enddo
         endif