Merge pull request NCAR#142 from mdtoyNOAA/gsl/develop_physics_mods_f…

…rom_RRFS_dev

GWD, LSM and MYNN physics updates from RRFS_dev branch

physics/GFS_surface_composites.F90

@@ -27,8 +27,9 @@ end subroutine GFS_surface_composites_pre_finalize
 !> \section arg_table_GFS_surface_composites_pre_run Argument Table
 !! \htmlinclude GFS_surface_composites_pre_run.html
 !!
-   subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, frac_grid,                                &
-                                 flag_cice, cplflx, cplice, cplwav2atm, landfrac, lakefrac, lakedepth, oceanfrac, frland, &
+   subroutine GFS_surface_composites_pre_run (im, xlat_d, xlon_d, flag_init, lsm_cold_start, lkm, frac_grid,                &
+                                 flag_cice, cplflx, cplice, cplwav2atm, lsm, lsm_ruc,                                     &
+                                 landfrac, lakefrac, lakedepth, oceanfrac, frland,                                        &
                                  dry, icy, lake, use_flake, wet, hice, cice, zorlo, zorll, zorli,                         &
                                  snowd,            snowd_lnd, snowd_ice, tprcp, tprcp_wat,                                &
                                  tprcp_lnd, tprcp_ice, uustar, uustar_wat, uustar_lnd, uustar_ice,                        &
@@ -40,10 +41,11 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
       implicit none
 
       ! Interface variables
-      integer,                             intent(in   ) :: im, lkm, kdt
-      logical,                             intent(in   ) :: flag_init, flag_restart, frac_grid, cplflx, cplice, cplwav2atm
+      integer,                             intent(in   ) :: im, lkm, kdt, lsm, lsm_ruc
+      logical,                             intent(in   ) :: flag_init, lsm_cold_start, frac_grid, cplflx, cplice, cplwav2atm
       logical, dimension(:),              intent(inout)  :: flag_cice
       logical,              dimension(:), intent(inout)  :: dry, icy, lake, use_flake, wet
+      real(kind=kind_phys), dimension(:), intent(in   )  :: xlat_d, xlon_d
       real(kind=kind_phys), dimension(:), intent(in   )  :: landfrac, lakefrac, lakedepth, oceanfrac
       real(kind=kind_phys), dimension(:), intent(inout)  :: cice, hice
       real(kind=kind_phys), dimension(:), intent(  out)  :: frland
@@ -201,12 +203,13 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
             endif
           endif
         enddo
-      endif
+      endif ! frac_grid
 
       do i=1,im
         tprcp_wat(i) = tprcp(i)
         tprcp_lnd(i) = tprcp(i)
         tprcp_ice(i) = tprcp(i)
+
         if (wet(i)) then                   ! Water
           uustar_wat(i) = uustar(i)
             tsfc_wat(i) = tsfco(i)
@@ -219,7 +222,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
         endif
         if (dry(i)) then                   ! Land
           uustar_lnd(i) = uustar(i)
-           weasd_lnd(i) = weasd(i)
+          if(lsm /= lsm_ruc) weasd_lnd(i) = weasd(i)
            tsurf_lnd(i) = tsfcl(i)
         ! DH*
         else
@@ -230,7 +233,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
         endif
         if (icy(i)) then                   ! Ice
           uustar_ice(i) = uustar(i)
-           weasd_ice(i) = weasd(i)
+          if(lsm /= lsm_ruc) weasd_ice(i) = weasd(i)
            tsurf_ice(i) = tisfc(i)
             ep1d_ice(i) = zero
             gflx_ice(i) = zero
@@ -279,7 +282,17 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
           endif
         enddo
       else
+       if(lsm /= lsm_ruc) then ! do not do snow initialization  with RUC lsm
         do i=1,im
+          !-- print ice point
+          !if ( (xlon_d(i) > 298.6) .and.  (xlon_d(i) < 298.7) .and. &
+          !   (xlat_d(i) >  68.6 ) .and.  (xlat_d(i) < 68.7 )) then
+          !  print *,'Composit weasd_ice(i),snowd_ice',kdt,i,xlat_d(i),xlon_d(i),weasd_ice(i),snowd_ice(i)
+          !endif
+          !if ( (xlon_d(i) > 284.35) .and.  (xlon_d(i) < 284.6) .and. &
+          ! (xlat_d(i) >  41.0 ) .and.  (xlat_d(i) < 41.2 )) then
+          !  print *,'Composit2  weasd_lnd(i),snowd_lnd',kdt,i,xlat_d(i),xlon_d(i),weasd_lnd(i),snowd_lnd(i)
+          !endif
           if (icy(i)) then
             if (kdt == 1 .or. (.not. cplflx .or. lakefrac(i) > zero)) then
               snowd_lnd(i) = zero
@@ -290,6 +303,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
             endif
           endif
         enddo
+       endif ! lsm/=lsm_ruc
       endif
 
 !     write(0,*)' minmax of ice snow=',minval(snowd_ice),maxval(snowd_ice)
@@ -644,6 +658,7 @@ subroutine GFS_surface_composites_post_run (
 
         do i=1,im
           if (islmsk(i) == 1) then
+          !-- land
             zorl(i)   = zorll(i)
             cd(i)     = cd_lnd(i)
             cdq(i)    = cdq_lnd(i)
@@ -669,6 +684,7 @@ subroutine GFS_surface_composites_post_run (
             hice(i)   = zero
             cice(i)   = zero
           elseif (islmsk(i) == 0) then
+          !-- water
             zorl(i)   = zorlo(i)
             cd(i)     = cd_wat(i)
             cdq(i)    = cdq_wat(i)
@@ -695,6 +711,7 @@ subroutine GFS_surface_composites_post_run (
             hice(i)   = zero
             cice(i)   = zero
           else ! islmsk(i) == 2
+          !-- ice
             zorl(i)   = zorli(i)
             cd(i)     = cd_ice(i)
             cdq(i)    = cdq_ice(i)

physics/GFS_surface_composites.meta

@@ -14,16 +14,32 @@
   dimensions = ()
   type = integer
   intent = in
+[xlat_d]
+  standard_name = latitude_in_degree
+  long_name = latitude in degree north
+  units = degree_north
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[xlon_d]
+  standard_name = longitude_in_degree
+  long_name = longitude in degree east
+  units = degree_east
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
 [flag_init]
   standard_name = flag_for_first_timestep
   long_name = flag signaling first time step for time integration loop
   units = flag
   dimensions = ()
   type = logical
   intent = in
-[flag_restart]
-  standard_name = flag_for_restart
-  long_name = flag for restart (warmstart) or coldstart
+[lsm_cold_start]
+  standard_name = do_lsm_cold_start
+  long_name = flag to signify LSM is cold-started
   units = flag
   dimensions = ()
   type = logical
@@ -70,6 +86,20 @@
   dimensions = ()
   type = logical
   intent = in
+[lsm]
+  standard_name = control_for_land_surface_scheme
+  long_name = flag for land surface model
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
+[lsm_ruc]
+  standard_name = identifier_for_ruc_land_surface_scheme
+  long_name = flag for RUC land surface model
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
 [landfrac]
   standard_name = land_area_fraction
   long_name = fraction of horizontal grid area occupied by land

physics/drag_suite.F90

@@ -589,17 +589,8 @@ subroutine drag_suite_run(                                           &
    endif
 enddo
 
-do i=1,im
-   if ( dx(i) .ge. dxmax_ss ) then
-      ss_taper(i) = 1.
-   else
-      if ( dx(i) .le. dxmin_ss) then
-         ss_taper(i) = 0.
-      else
-         ss_taper(i) = dxmax_ss * (1. - dxmin_ss/dx(i))/(dxmax_ss-dxmin_ss)
-      endif
-   endif
-enddo
+! Remove ss_tapering
+ss_taper(:) = 1.
 
 ! SPP, if spp_gwd is 0, no perturbations are applied.
 if ( spp_gwd==1 ) then
@@ -986,13 +977,11 @@ subroutine drag_suite_run(                                           &
          enddo
          if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
             if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
-              cleff_ss    = sqrt(dxy(i)**2 + dxyp(i)**2)   ! WRF
+              ! Modify xlinv to represent wave number of "typical" small-scale topography 
 !              cleff_ss    = 3. * max(dx(i),cleff_ss)
 !              cleff_ss    = 10. * max(dxmax_ss,cleff_ss)
-              cleff_ss    = 0.1 * max(dxmax_ss,cleff_ss)  ! WRF
 !               cleff_ss    = 0.1 * 12000.
-              coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
-              xlinv(i) = coefm_ss(i) / cleff_ss
+              xlinv(i) = 0.001*pi   ! 2km horizontal wavelength
               !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts)))
               govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
               !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i))
@@ -1003,8 +992,8 @@ subroutine drag_suite_run(                                           &
                 !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i))
                 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2)
                 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
-                var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) +                       &
-                              MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
+                ! Remove limit on varss_stoch
+                var_temp = varss_stoch(i)
                 ! Note:  This is a semi-implicit treatment of the time differencing
                 var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)  ! this is greater than zero
                 tauwavex0=var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
@@ -1018,8 +1007,8 @@ subroutine drag_suite_run(                                           &
                 !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i))
                 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2)
                 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
-                var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) +                       &
-                              MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
+                ! Remove limit on varss_stoch
+                var_temp = varss_stoch(i)
                 ! Note:  This is a semi-implicit treatment of the time differencing
                 var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)  ! this is greater than zero
                 tauwavey0=var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
@@ -1083,17 +1072,16 @@ subroutine drag_suite_run(                                           &
 
          IF ((xland(i)-1.5) .le. 0.) then
             !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
-            var_temp = MIN(varss_stoch(i),varmax_fd_stoch(i)) +                           &
-                       MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i)))
-            var_temp = MIN(var_temp, 250.)
+            ! Remove limit on varss_stoch
+            var_temp = varss_stoch(i)
+            !var_temp = MIN(var_temp, 250.)
             a1=0.00026615161*var_temp**2
 !           a1=0.00026615161*MIN(varss(i),varmax)**2
 !           a1=0.00026615161*(0.5*varss(i))**2
            ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
             a2=a1*0.005363
-           ! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018
-            H_efold = max(2*varss_stoch(i),hpbl(i))
-            H_efold = min(H_efold,1500.)
+            ! Beljaars H_efold
+            H_efold = 1500.
             DO k=kts,km
                wsp=SQRT(u1(i,k)**2 + v1(i,k)**2)
                ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759

physics/module_sf_mynn.F90

@@ -112,7 +112,7 @@ MODULE module_sf_mynn
                                         !1: check input
                                         !2: everything - heavy I/O
   LOGICAL, PARAMETER :: compute_diag = .false.
-  LOGICAL, PARAMETER :: compute_flux = .false.  !shouldn't need compute 
+  LOGICAL, PARAMETER :: compute_flux = .true.  !shouldn't need compute 
                ! these in FV3. They will be written over anyway.
                ! Computing the fluxes here is leftover from the WRF world.
 

physics/module_sf_ruclsm.F90

@@ -2595,7 +2595,8 @@ SUBROUTINE SOIL (debug_print,                        &
 ! 3feb21 - in RRFS testing (fv3-based), ref*0.5 gives too much direct
 !          evaporation. Therefore , it is replaced with ref*0.7.
         !fc=max(qmin,ref*0.5)
-        fc=max(qmin,ref*0.7)
+        !fc=max(qmin,ref*0.7)
+        fc=ref
         fex_fc=1.
       if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > 0.) then
         soilres = 1.