GFS.v16.0.4 physics update

gfsphysics/GFS_layer/GFS_physics_driver.F90

@@ -2387,7 +2387,7 @@ subroutine GFS_physics_driver                                 &
                 call satmedmfvdifq(ix, im, levs, nvdiff, ntcw, ntiw, ntke,          &
                        dvdt, dudt, dtdt, dqdt,                                      &
                        Statein%ugrs, Statein%vgrs, Statein%tgrs, Statein%qgrs,      &
-                       Radtend%htrsw, Radtend%htrlw, xmu, garea,                    &
+                       Radtend%htrsw, Radtend%htrlw, xmu, garea, islmsk, snowd3,    &
                        Statein%prsik(1,1), rb, Sfcprop%zorl, Diag%u10m, Diag%v10m,  &
                        Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%tsfc, hflxq, evapq,      &
                        stress, wind, kpbl, Statein%prsi, del, Statein%prsl,         &
@@ -2640,7 +2640,7 @@ subroutine GFS_physics_driver                                 &
                 call satmedmfvdifq(ix, im, levs, nvdiff, ntcw, ntiwx, ntkev,          &
                          dvdt, dudt, dtdt, dvdftra,                                   &
                          Statein%ugrs, Statein%vgrs, Statein%tgrs, vdftra,            &
-                         Radtend%htrsw, Radtend%htrlw, xmu, garea,                    &
+                         Radtend%htrsw, Radtend%htrlw, xmu, garea, islmsk, snowd3,    &
                          Statein%prsik(1,1), rb, Sfcprop%zorl, Diag%u10m, Diag%v10m,  &
                          Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%tsfc, hflxq, evapq,      &
                          stress, wind, kpbl, Statein%prsi, del, Statein%prsl,         &

gfsphysics/physics/radsw_datatb.f

@@ -2552,7 +2552,7 @@ module module_radsw_sflux          !
 !> band index (3rd index in array sfluxref described below)
       integer, dimension(nblow:nbhgh), public :: ibx
 
-      data  layreffr/ 18,30, 6, 3, 3, 8, 2, 6, 1, 2, 0,32,58,49 /
+      data  layreffr/ 18,30, 6, 3, 3, 8, 2, 6, 1, 2, 0,32,42,49 /
       data  ix1     /  1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 3, 0 /
       data  ix2     /  5, 2, 5, 2, 0, 2, 6, 0, 6, 0, 0, 0, 6, 0 /
       data  ibx     /  1, 1, 1, 2, 2, 3, 4, 3, 5, 4, 5, 6, 2, 7 /

gfsphysics/physics/samfdeepcnv.f

@@ -217,7 +217,7 @@ subroutine samfdeepcnv(im,ix,km,delt,itc,ntc,ntk,ntr,delp,
       parameter(clamd=0.03,tkemx=0.65,tkemn=0.05)
       parameter(dtke=tkemx-tkemn)
       parameter(dbeta=0.1)
-      parameter(cthk=200.,dthk=25.)
+      parameter(cthk=150.,dthk=25.)
       parameter(cinpcrmx=180.,cinpcrmn=120.)
 !     parameter(cinacrmx=-120.,cinacrmn=-120.)
       parameter(cinacrmx=-120.,cinacrmn=-80.)
@@ -1191,9 +1191,7 @@ subroutine samfdeepcnv(im,ix,km,delt,itc,ntc,ntk,ntr,delp,
 !
           k = kbcon(i)
           dp = 1000. * del(i,k)
-          xmbmax(i) = dp / (2. * g * dt2)
-!
-!         xmbmax(i) = dp / (g * dt2)
+          xmbmax(i) = dp / (g * dt2)
 !
 !         mbdt(i) = 0.1 * dp / g
 !

gfsphysics/physics/satmedmfvdifq.f

@@ -19,7 +19,8 @@
 !
 !----------------------------------------------------------------------
       subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
-     &     dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,
+!wz  &     dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,
+     &     dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,islimsk,snwdph,
      &     psk,rbsoil,zorl,u10m,v10m,fm,fh,
      &     tsea,heat,evap,stress,spd1,kpbl,
      &     prsi,del,prsl,prslk,phii,phil,delt,
@@ -38,6 +39,10 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
 !----------------------------------------------------------------------
       integer ix, im, km, ntrac, ntcw, ntiw, ntke
       integer kpbl(im), kinver(im)
+!
+!wz
+      integer islimsk(im)
+      real(kind=kind_phys), dimension(im,3), intent(in)    ::  snwdph
 !
       real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s, dspfac,
      &                     bl_upfr, bl_dnfr
@@ -156,6 +161,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
      &                     zlup,    zldn,   bsum,
      &                     tem,     tem1,   tem2,
      &                     ptem,    ptem0,  ptem1,  ptem2
+!wz
+      real(kind=kind_phys)       xkzm_mp, xkzm_hp
 !
       real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck
 !
@@ -173,7 +180,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
       parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1)
       parameter(vk=0.4,rimin=-100.)
       parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3)
-      parameter(rlmn=30.,rlmn1=5.,rlmn2=15.)
+      parameter(rlmn=30.,rlmn1=5.,rlmn2=10.)
       parameter(rlmx=300.,elmx=300.)
       parameter(prmin=0.25,prmax=4.0)
       parameter(pr0=1.0,prtke=1.0,prscu=0.67)
@@ -183,7 +190,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
       parameter(aphi5=5.,aphi16=16.)
       parameter(elmfac=1.0,elefac=1.0,cql=100.)
       parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=5000.)
-      parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15)
+      parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.1)
       parameter(h1=0.33333333)
       parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15)
       parameter(ce0=0.4)
@@ -254,17 +261,37 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
 !  xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.)
 !  xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.)
 !
+!wz
       do i=1,im
+              xkzm_mp = xkzm_m
+              xkzm_hp = xkzm_h
+!
+         if( islimsk(i) == 1 .and. snwdph(i,1) > 10.0 ) then    ! over land 
+            if (rbsoil(i) > 0. .and. rbsoil(i) <= 0.25) then
+              xkzm_mp = xkzm_m * (1.0 - rbsoil(i)/0.25)**2 +
+     &                     0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2)
+              xkzm_hp = xkzm_h * (1.0 - rbsoil(i)/0.25)**2 +
+     &                     0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2)
+            else if (rbsoil(i) > 0.25) then
+              xkzm_mp = 0.1
+              xkzm_hp = 0.1
+            endif
+         endif
+!#
         kx1(i) = 1
         tx1(i) = 1.0 / prsi(i,1)
         tx2(i) = tx1(i)
         if(gdx(i) >= xkgdx) then
-          xkzm_hx(i) = xkzm_h
-          xkzm_mx(i) = xkzm_m
+!wz       xkzm_hx(i) = xkzm_h
+!wz       xkzm_mx(i) = xkzm_m
+          xkzm_hx(i) = xkzm_hp
+          xkzm_mx(i) = xkzm_mp
         else
           tem  = 1. / (xkgdx - 5.)
-          tem1 = (xkzm_h - 0.01) * tem
-          tem2 = (xkzm_m - 0.01) * tem
+!wz       tem1 = (xkzm_h - 0.01) * tem
+!wz       tem2 = (xkzm_m - 0.01) * tem
+          tem1 = (xkzm_hp - 0.01) * tem
+          tem2 = (xkzm_mp - 0.01) * tem
           ptem = gdx(i) - 5.
           xkzm_hx(i) = 0.01 + tem1 * ptem
           xkzm_mx(i) = 0.01 + tem2 * ptem
@@ -699,7 +726,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke,
 !         tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k)
 !         if(tem1 > 1.e-5) then
           tem1 = tvx(i,k+1)-tvx(i,k)
-          if(tem1 > 0.) then
+          if(tem1 > 0. .and. islimsk(i) /= 1 ) then
              xkzo(i,k)  = min(xkzo(i,k), xkzinv)
              xkzmo(i,k) = min(xkzmo(i,k), xkzinv)
              rlmnz(i,k) = min(rlmnz(i,k), rlmn2)

gfsphysics/physics/sfc_diff.f

@@ -139,17 +139,12 @@ subroutine sfc_diff(im,ps,t1,q1,z1, wind,                 !intent(in)
             z0max = max(z0max, 1.0e-6)
 
 !           czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil
-!           czilc = 0.8
-
-!           tem1  = 1.0 - sigmaf(i)
-!           ztmax = z0max*exp( - tem1*tem1
-!    &                       * czilc*ca*sqrt(ustar(i,1)*(0.01/1.5e-05)))
-!
-            czilc = 10.0 ** (- 4. * z0max) ! Trier et al. (2011, WAF)
-            ztmax = z0max * exp( - czilc * ca
-     &            * 258.2 * sqrt(ustar(i,1)*z0max) )
-
+            czilc = 0.8
 
+            tem1  = 1.0 - sigmaf(i)
+            ztmax = z0max*exp( - tem1*tem1
+     &              * czilc*ca*sqrt(ustar(i,1)*(0.01/1.5e-05)))
+
 ! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land
             if (ztpert(i) /= 0.0) then
               ztmax = ztmax * (10.**ztpert(i))
@@ -183,16 +178,11 @@ subroutine sfc_diff(im,ps,t1,q1,z1, wind,                 !intent(in)
 
 !           czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height
 !           dependance of czil
-!           czilc = 0.8
+            czilc = 0.8
 
-!           tem1 = 1.0 - sigmaf(i)
-!           ztmax = z0max*exp( - tem1*tem1
-!    &                       * czilc*ca*sqrt(ustar(i,2)*(0.01/1.5e-05)))
-!
-            czilc = 10.0 ** (- 4. * z0max) ! Trier et al. (2011, WAF)
-            ztmax = z0max * exp( - czilc * ca
-     &            * 258.2 * sqrt(ustar(i,2)*z0max) )
-!
+            tem1 = 1.0 - sigmaf(i)
+            ztmax = z0max*exp( - tem1*tem1
+     &              * czilc*ca*sqrt(ustar(i,2)*(0.01/1.5e-05)))
             ztmax = max(ztmax, 1.0e-6)
 !
             call stability

io/FV3GFS_io.F90

@@ -1176,6 +1176,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
       endif
     endif
 
+#ifdef CCPP
     if (nint(sfc_var3ice(1,1,1)) == -9999) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing tiice')
       do nb = 1, Atm_block%nblks
@@ -1186,6 +1187,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain)
       enddo
     endif
 
+#endif
 !#endif
 
     if(Model%frac_grid) then ! 3-way composite