Merge branch 'ufs/dev' into feature/gfdlmpv3

Author: grantfirl

physics/GWD/drag_suite.F90

@@ -1731,7 +1731,7 @@ subroutine drag_suite_psl(                                           &
    real(kind=kind_phys),parameter       :: odmax  = 10.
    real(kind=kind_phys),parameter       :: cdmin  = 0.0
    integer              :: komax(im),kbmax(im),kblk(im)
-   real(kind=kind_phys)                 :: hmax(im)
+   real(kind=kind_phys)                 :: href(im),hmax(im)
    real(kind=kind_phys)                 :: cd
    real(kind=kind_phys)                 :: zblk,tautem
    real(kind=kind_phys)                 :: pe,ke
@@ -1916,6 +1916,7 @@ subroutine drag_suite_psl(                                           &
 !  initialize array for flow-blocking drag
 !
    taufb(1:im,1:km+1) = 0.0
+   href(1:im) = 0.0
    hmax(1:im) = 0.0
    komax(1:im) = 0
    kbmax(1:im) = 0
@@ -2000,11 +2001,13 @@ subroutine drag_suite_psl(                                           &
 !
    do i = its,im
      hmax(i) = max(elvmax(i),zlowtop(i))
+     href(i) = max(hmax(i),hpbl(i))
    enddo
 !
    do i = its,im
 !!!     kbl(i)   = max(kpbl(i), klowtop(i))    ! do not use pbl height for the time being...
-     kbl(i)   = max(komax(i), klowtop(i))    
+     kbl(i)   = max(komax(i), klowtop(i))
+     kbl(i)   = max(kbl(i), kpbl(i))    
      kbl(i)   = max(min(kbl(i),kpblmax),kpblmin)
    enddo
 !
@@ -2109,7 +2112,7 @@ subroutine drag_suite_psl(                                           &
 !
 !  no drag when sub-oro is too small..
 !
-         ldrag(i) = hmax(i).le.hmt_min
+         ldrag(i) = href(i).le.hmt_min
 !
 !  no drag when critical level in the base layer
 !

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90

@@ -702,14 +702,16 @@ end subroutine GFS_phys_time_vary_init
 !>\section gen_GFS_phys_time_vary_timestep_init GFS_phys_time_vary_timestep_init General Algorithm
 !> @{
       subroutine GFS_phys_time_vary_timestep_init (                                                 &
-            me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, nsswr, fhswr, lsswr, fhour, &
+            me, master, cnx, cny, isc, jsc, nrcm, im, levs, kdt, idate, cplflx,                     &
+            nsswr, fhswr, lsswr, fhour,                                                             &
             imfdeepcnv, cal_pre, random_clds, nscyc, ntoz, h2o_phys, iaerclm, iccn, clstp,          &
             jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, iflip,            &
             jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm,               &
             jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, fn_nml,      &
             imap, jmap, prsl, seed0, rann, nthrds, nx, ny, nsst, tile_num, nlunit, lsoil, lsoil_lsm,&
             kice, ialb, isot, ivegsrc, input_nml_file, use_ufo, nst_anl, frac_grid, fhcyc, phour,   &
-            lakefrac, min_seaice, min_lakeice, smc, slc, stc, smois, sh2o, tslb, tiice, tg3, tref,  &
+            oceanfrac, lakefrac, min_seaice, min_lakeice, smc, slc, stc, smois, sh2o, tslb, tiice,  &
+            tg3, tref,                                                                              &
             tsfc, tsfco, tisfc, hice, fice, facsf, facwf, alvsf, alvwf, alnsf, alnwf, zorli, zorll, &
             zorlo, weasd, slope, snoalb, canopy, vfrac, vtype, stype,scolor, shdmin, shdmax, snowd, &
             cv, cvb, cvt, oro, oro_uf, xlat_d, xlon_d, slmsk, landfrac, ozphys, h2ophys,            &
@@ -723,7 +725,7 @@ subroutine GFS_phys_time_vary_timestep_init (
                                                 nsswr, imfdeepcnv, iccn, nscyc, ntoz, iflip
          integer,              intent(in)    :: idate(:)
          real(kind_phys),      intent(in)    :: fhswr, fhour
-         logical,              intent(in)    :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm
+         logical,              intent(in)    :: lsswr, cal_pre, random_clds, h2o_phys, iaerclm, cplflx
          real(kind_phys),      intent(out)   :: clstp
          integer,              intent(in), optional    :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:)
          real(kind_phys),      intent(in), optional    :: ddy_o3(:),  ddy_h(:)
@@ -753,7 +755,7 @@ subroutine GFS_phys_time_vary_timestep_init (
          character(len=*),     intent(in)    :: fn_nml
          logical,              intent(in)    :: use_ufo, nst_anl, frac_grid
          real(kind_phys),      intent(in)    :: fhcyc, phour, lakefrac(:), min_seaice, min_lakeice,  &
-                                                xlat_d(:), xlon_d(:), landfrac(:)
+                                                xlat_d(:), xlon_d(:), landfrac(:),oceanfrac(:)
          real(kind_phys),      intent(inout) :: smc(:,:), slc(:,:), stc(:,:), tiice(:,:), tg3(:),    &
                                       tsfc(:), tsfco(:), tisfc(:), hice(:), fice(:),                 &
                                       facsf(:), facwf(:), alvsf(:), alvwf(:), alnsf(:), alnwf(:),    &
@@ -926,6 +928,7 @@ subroutine GFS_phys_time_vary_timestep_init (
                  tsfco, tisfc, hice, fice, facsf, facwf, alvsf, alvwf, alnsf, alnwf,         &
                  zorli, zorll, zorlo, weasd, slope, snoalb, canopy, vfrac, vtype,            &
                  stype, scolor, shdmin, shdmax, snowd, cv, cvb, cvt, oro, oro_uf,            &
+                 cplflx, oceanfrac,                                                          &
                  xlat_d, xlon_d, slmsk, imap, jmap, errmsg, errflg)
            endif
          endif

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.meta

@@ -1173,6 +1173,13 @@
   dimensions = (4)
   type = integer
   intent = in
+[cplflx]
+  standard_name = flag_for_surface_flux_coupling
+  long_name = flag controlling cplflx collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [nsswr]
   standard_name = number_of_timesteps_between_shortwave_radiation_calls
   long_name = number of timesteps between shortwave radiation calls
@@ -1638,6 +1645,14 @@
   type = real
   kind = kind_phys
   intent = in
+[oceanfrac]
+  standard_name = sea_area_fraction
+  long_name = fraction of horizontal grid area occupied by ocean
+  units = frac
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
 [lakefrac]
   standard_name = lake_area_fraction
   long_name = fraction of horizontal grid area occupied by lake

physics/Interstitials/UFS_SCM_NEPTUNE/gcycle.F90

@@ -22,6 +22,7 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, fn_nml,
       tsfco, tisfc, hice, fice, facsf, facwf, alvsf, alvwf, alnsf, alnwf,          &
       zorli, zorll, zorlo, weasd, slope, snoalb, canopy, vfrac, vtype,             &
       stype, scolor, shdmin, shdmax, snowd, cv, cvb, cvt, oro, oro_uf,             &
+      cplflx, oceanfrac,                                                           &
       xlat_d, xlon_d, slmsk, imap, jmap, errmsg, errflg)
 !
 !
@@ -33,9 +34,9 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, fn_nml,
     integer,              intent(in)    :: idate(:), ialb, isot, ivegsrc
     character(len = 64), intent(in)     :: fn_nml
     character(len=*),     intent(in)    :: input_nml_file(:)
-    logical,              intent(in)    :: use_ufo, nst_anl, frac_grid
+    logical,              intent(in)    :: use_ufo, nst_anl, frac_grid, cplflx
     real(kind=kind_phys), intent(in)    :: fhcyc, phour, landfrac(:), lakefrac(:), &
-                                           min_seaice, min_lakeice,                &
+                                           min_seaice, min_lakeice,oceanfrac(:), &
                                            xlat_d(:), xlon_d(:)
     real(kind=kind_phys), intent(inout), optional ::   &
                                            smois(:,:), &
@@ -103,6 +104,17 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, fn_nml,
         STCFC1 (nx*ny*max(lsoil,lsoil_lsm)), &
         SLCFC1 (nx*ny*max(lsoil,lsoil_lsm))
 
+!
+! declare the variables (arrays) for cplflx, surface type dependent gcycle changes
+!
+    real(kind=kind_io8) ::                   &
+        hice_save    (nx*ny),                &        ! sea or lake ice thickness
+        fice_save    (nx*ny),                &        ! sea or lake ice fraction
+        snowd_save   (nx*ny),                &        ! water equivalent snow depth
+        snoalb_save  (nx*ny),                &        ! maximum snow albedo
+        tisfc_save   (nx*ny),                &        ! surface skin temperature over (sea or lake) ice
+        weasd_save   (nx*ny)                          ! water equiv of acc snow depth over land and (sea or lake)  ice
+
 
     real (kind=kind_io8) :: min_ice(nx*ny)
     integer              :: i_indx(nx*ny), j_indx(nx*ny)
@@ -131,11 +143,25 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, fn_nml,
       sig1t = 0.0_kind_phys
       npts  = nx*ny
 !
+! Some surface variables need to be updated by gcycle with coupled mode, and nsst mode dependent. A few variables are saved 
+! in order to be able to update them over the specific surface types only after call sfccycle
+!
+      if ( cplflx ) then
+        hice_save   = hice 
+        fice_save   = fice
+        snowd_save  = snowd
+        snoalb_save = snoalb
+        tisfc_save  = tisfc
+        weasd_save  = weasd
+      endif
+
       if ( nsst > 0 ) then
         TSFFCS = tref
       else
         TSFFCS = tsfco
-      end if
+      endif
+
+
 ! integer to real/double precision
       slpfcs = real(slope)
       vegfcs = real(vtype)
@@ -251,11 +277,45 @@ subroutine gcycle (me, nthrds, nx, ny, isc, jsc, nsst, tile_num, nlunit, fn_nml,
       close (nlunit)
 #endif
 !
-      if ( nsst > 0 ) then
-        tref = TSFFCS
+! The gcycle resulted change is applied to some variables in the way of the coupled mode dependent, water surface type (ocean or lake)
+! dependent and nsst mode dependent
+!
+      if ( cplflx ) then
+!       In coupled mode, keep these variables the same as is (before sfccycle is called) over ocean
+        do ix=1,npts
+          if ( oceanfrac(ix) > 0.0_kind_phys ) then      
+            hice(ix)   = hice_save(ix) 
+            fice(ix)   = fice_save(ix)
+            snowd(ix)  = snowd_save(ix)
+            snoalb(ix) = snoalb_save(ix)
+            tisfc(ix)  = tisfc_save(ix)
+            weasd(ix)  = weasd_save(ix)
+          endif
+        enddo
+!       In the coupled mode and when NSST is on, update tref over non-ocean 
+        if ( nsst > 0 ) then       
+          do ix=1,npts
+            if ( oceanfrac(ix) == 0.0_kind_phys ) then 
+              tref(ix)  = TSFFCS(ix) 
+            endif
+          enddo
+!       In the coupled mode and when NSST is off, update tsfc and tsfco over non-ocean
+        else             
+          do ix=1,npts
+            if ( oceanfrac(ix) == 0.0_kind_phys ) then 
+              tsfc(ix)  = TSFFCS(ix)
+              tsfco(ix) = TSFFCS(ix)
+            endif
+          enddo
+        endif
+!     The same as before (this modification) in uncoupled mode
       else
-        tsfc  = TSFFCS
-        tsfco = TSFFCS
+        if ( nsst > 0 ) then
+          tref = TSFFCS
+        else
+          tsfc  = TSFFCS
+          tsfco = TSFFCS
+        endif
       endif
 !
 ! real/double precision to integer

physics/SFC_Layer/UFS/sfc_diff.f

@@ -60,7 +60,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
      &                    sigmaf,vegtype,shdmax,ivegsrc,                &  !intent(in)
      &                    z0pert,ztpert,                                &  ! mg, sfc-perts !intent(in)
      &                    flag_iter,redrag,                             &  !intent(in)
-     &                    flag_lakefreeze,                              &  !intent(in)             
+     &                    flag_lakefreeze,lakefrac,fice,                &  !intent(in)             
      &                    u10m,v10m,sfc_z0_type,                        &  !hafs,z0 type !intent(in)
      &                    u1,v1,usfco,vsfco,icplocn2atm,                &  
      &                    wet,dry,icy,                                  &  !intent(in)
@@ -106,6 +106,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
      &                    ps,t1,q1,z1,garea,prsl1,prslki,prsik1,prslk1, &
      &                    wind,sigmaf,shdmax,                           &
      &                    z0pert,ztpert ! mg, sfc-perts
+      real(kind=kind_phys), dimension(:), intent(in)    :: lakefrac
+      real(kind=kind_phys), dimension(:), intent(in)    :: fice
       real(kind=kind_phys), dimension(:), intent(in)    ::              &
      &                    tskin_wat, tskin_lnd, tskin_ice,              &
      &                    tsurf_wat, tsurf_lnd, tsurf_ice
@@ -340,7 +342,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
 !  ---  outputs:
      &      rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i),
      &      cm_ice(i), ch_ice(i), stress_ice(i), ustar_ice(i))
-      endif ! Icy points
+          endif ! Icy points
 
 ! BWG: Everything from here to end of subroutine was after
 !      the stuff now put into "stability"
@@ -364,20 +366,8 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
               windrel=sqrt((u1(i)-usfco(i))**2+(v1(i)-vsfco(i))**2) 
             endif
 
-            if (sfc_z0_type == -1) then    ! using wave model derived momentum roughness
-              tem1 = 0.11 * vis / ustar_wat(i)
-              z0 = tem1 +  0.01_kp * z0rl_wav(i)
-
-              if (redrag) then
-                z0max = max(min(z0, z0s_max),1.0e-7_kp)
-              else
-                z0max = max(min(z0,0.1_kp), 1.0e-7_kp)
-              endif
-              z0rl_wat(i) = 100.0_kp * z0max   ! cm
-            else
-              z0    = 0.01_kp * z0rl_wat(i)
-              z0max = max(zmin, min(z0,z1(i)))
-            endif
+            z0    = 0.01_kp * z0rl_wat(i)
+            z0max = max(zmin, min(z0,z1(i)))
 !
 !**  test xubin's new z0
 
@@ -416,52 +406,51 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
 !
 !  update z0 over ocean
 !
-            if (sfc_z0_type >= 0) then
-              if (sfc_z0_type == 0) then
-!               z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i)
-                tem1 = 0.11 * vis / ustar_wat(i)
-                z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i)
+            if ((sfc_z0_type == -1) .and.
+     &        (lakefrac(i) == 0.0 .and. fice(i) == 0.0) .and.
+     &        (z0rl_wav(i)>1.0e-7_kp .and. z0rl_wav(i)<0.1_kp)) then
+              ! using wave model derived momentum roughness
+              tem1 = 0.11 * vis / ustar_wat(i)
+              z0 = tem1 +  0.01_kp * z0rl_wav(i)
+
+              if (redrag) then
+                z0rl_wat(i) = 100.0_kp * max(min(z0,z0s_max),1.0e-7_kp)
+              else
+                z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.e-7_kp)
+              endif
+
+            elseif ((sfc_z0_type == 0) .or.
+     &        ((sfc_z0_type == -1) .and.
+     &        (z0rl_wav(i)<=1.0e-7_kp .or. z0rl_wav(i)>=0.1_kp))) then
+!             z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i)
+              tem1 = 0.11 * vis / ustar_wat(i)
+              z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i)
 
 
 ! mbek -- toga-coare flux algorithm
-!               z0 = (charnock / grav) * ustar(i)*ustar(i) +  arnu/ustar(i)
+!             z0 = (charnock / grav) * ustar(i)*ustar(i) +  arnu/ustar(i)
 !  new implementation of z0
-!               cc = ustar(i) * z0 / rnu
-!               pp = cc / (1. + cc)
-!               ff = grav * arnu / (charnock * ustar(i) ** 3)
-!               z0 = arnu / (ustar(i) * ff ** pp)
-
-                if (redrag) then
-                  z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),        &
-     &                                                 1.0e-7_kp)
-                else
-                  z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.e-7_kp)
-                endif
-
-              elseif (sfc_z0_type == 6) then   ! wang
-                 call znot_m_v6(wind10m, z0)   ! wind, m/s, z0, m
-                 z0rl_wat(i) = 100.0_kp * z0   ! cm
-              elseif (sfc_z0_type == 7) then   ! wang
-                 call znot_m_v7(wind10m, z0)   ! wind, m/s, z0, m
-                 z0rl_wat(i) = 100.0_kp * z0   ! cm
+!             cc = ustar(i) * z0 / rnu
+!             pp = cc / (1. + cc)
+!             ff = grav * arnu / (charnock * ustar(i) ** 3)
+!             z0 = arnu / (ustar(i) * ff ** pp)
+
+              if (redrag) then
+                z0rl_wat(i) = 100.0_kp * max(min(z0,z0s_max),1.0e-7_kp)
               else
-                 z0rl_wat(i) = 1.0e-4_kp
+                z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.e-7_kp)
               endif
 
-           elseif (z0rl_wav(i) <= 1.0e-7_kp .or.
-     &             z0rl_wav(i) > 1.0_kp) then
-!            z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i)
-             tem1 = 0.11 * vis / ustar_wat(i)
-             z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i)
-
-             if (redrag) then
-               z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),1.0e-7_kp)
-             else
-               z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.0e-7_kp)
-             endif
-
+            elseif (sfc_z0_type == 6) then   ! wang
+               call znot_m_v6(wind10m, z0)   ! wind, m/s, z0, m
+               z0rl_wat(i) = 100.0_kp * z0   ! cm
+            elseif (sfc_z0_type == 7) then   ! wang
+               call znot_m_v7(wind10m, z0)   ! wind, m/s, z0, m
+               z0rl_wat(i) = 100.0_kp * z0   ! cm
+            else
+               z0rl_wat(i) = 1.0e-4_kp
             endif
-
+!
           endif              ! end of if(open ocean)
 !
         endif                ! end of if(flagiter) loop

physics/SFC_Layer/UFS/sfc_diff.meta

@@ -143,6 +143,22 @@
   type = real
   kind = kind_phys
   intent = in
+[lakefrac]
+  standard_name = lake_area_fraction
+  long_name = fraction of horizontal grid area occupied by lake
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[fice]
+  standard_name = sea_ice_area_fraction_of_sea_area_fraction
+  long_name = ice fraction over open water
+  units = frac
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
 [vegtype]
   standard_name = vegetation_type_classification
   long_name = vegetation type at each grid cell