From f7a9d915269de609bf2597e00f262d752bdf0e94 Mon Sep 17 00:00:00 2001 From: Michael Toy Date: Thu, 24 Mar 2022 15:56:20 +0000 Subject: [PATCH] GWD, LSM and MYNN physics updates from RRFS_dev branch --- physics/GFS_surface_composites.F90 | 31 +++++++++++++++++------ physics/GFS_surface_composites.meta | 36 ++++++++++++++++++++++++--- physics/drag_suite.F90 | 38 ++++++++++------------------- physics/module_sf_mynn.F90 | 2 +- physics/module_sf_ruclsm.F90 | 3 ++- 5 files changed, 73 insertions(+), 37 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 510b3f427..f44df5890 100644 --- a/physics/GFS_surface_composites.F90 +++ b/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) diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 89048e487..40f0c940c 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -14,6 +14,22 @@ 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 @@ -21,9 +37,9 @@ 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 diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90 index 7fea98b13..fe9095210 100644 --- a/physics/drag_suite.F90 +++ b/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 diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 5f227750a..65e83c93d 100644 --- a/physics/module_sf_mynn.F90 +++ b/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. diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index b39610bc8..01e9c1100 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/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.