From 833f7dc7520a4e57d52b94a1ed5f4df3ee7b1684 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Sat, 18 May 2024 09:56:56 -0400 Subject: [PATCH 01/12] HR4 GWD update for files drag_suite.F90 ugwpv1_gsldrag.F90 unified_ugwp.F90 module_sf_noahmplsm.F90 noahmpdrv.F90 --- physics/GWD/drag_suite.F90 | 1217 +++++++++++++++++ physics/GWD/drag_suite.meta | 15 + physics/GWD/ugwpv1_gsldrag.F90 | 37 +- physics/GWD/ugwpv1_gsldrag.meta | 22 + physics/GWD/unified_ugwp.F90 | 36 +- physics/GWD/unified_ugwp.meta | 46 + .../Land/Noahmp/module_sf_noahmplsm.F90 | 59 +- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 14 +- physics/SFC_Models/Land/Noahmp/noahmpdrv.meta | 23 + 9 files changed, 1446 insertions(+), 23 deletions(-) diff --git a/physics/GWD/drag_suite.F90 b/physics/GWD/drag_suite.F90 index ff68f4216..2dfd0a598 100644 --- a/physics/GWD/drag_suite.F90 +++ b/physics/GWD/drag_suite.F90 @@ -1413,6 +1413,1223 @@ subroutine drag_suite_run( & return end subroutine drag_suite_run !------------------------------------------------------------------- + + subroutine drag_suite_psl( & + & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, & + & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, & + & var,oc1,oa4,ol4, & + & varss,oc1ss,oa4ss,ol4ss, & + & THETA,SIGMA,GAMMA,ELVMAX, & + & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, & + & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, & + & dusfc,dvsfc, & + & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & + & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & + & slmsk,br1,hpbl,vtype, & + & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, & + & lprnt, ipr, rdxzb, dx, gwd_opt, & + & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & + & psl_gwd_dx_factor, & + & dtend, dtidx, index_of_process_orographic_gwd, & + & index_of_temperature, index_of_x_wind, & + & index_of_y_wind, ldiag3d, ldiag_ugwp, ugwp_seq_update, & + & spp_wts_gwd, spp_gwd, errmsg, errflg) + +! ******************************************************************** +! -----> I M P L E M E N T A T I O N V E R S I O N <---------- +! +! ----- This code ----- +!begin WRF code + +! this code handles the time tendencies of u v due to the effect of mountain +! induced gravity wave drag from sub-grid scale orography. this routine +! not only treats the traditional upper-level wave breaking due to mountain +! variance (alpert 1988), but also the enhanced lower-tropospheric wave +! breaking due to mountain convexity and asymmetry (kim and arakawa 1995). +! thus, in addition to the terrain height data in a model grid box, +! additional 10-2d topographic statistics files are needed, including +! orographic standard deviation (var), convexity (oc1), asymmetry (oa4) +! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography +! hong (1999). the current scheme was implmented as in hong et al.(2008) +! +! Originally coded by song-you hong and young-joon kim and implemented by song-you hong +! +! program history log: +! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle +! with blocked height by dividing streamline theory +! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale +! orographic grabity wave drag: +! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the +! topographic form drag of Beljaars et al. (2004, QJRMS) +! Activation of each component is done by specifying the integer-parameters +! (defined below) to 0: inactive or 1: active +! gwd_opt_ls = 0 or 1: large-scale +! gwd_opt_bl = 0 or 1: blocking drag +! gwd_opt_ss = 0 or 1: small-scale gravity wave drag +! gwd_opt_fd = 0 or 1: topographic form drag +! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating +! gsd_diss_ht_opt = 0: dissipation heating off +! gsd_diss_ht_opt = 1: dissipation heating on +! 2020-08-25 Michael Toy changed logic control for drag component selection +! for CCPP. +! Namelist options: +! do_gsl_drag_ls_bl - logical flag for large-scale GWD + blocking +! do_gsl_drag_ss - logical flag for small-scale GWD +! do_gsl_drag_tofd - logical flag for turbulent form drag +! Compile-time options (same as before): +! gwd_opt_ls = 0 or 1: large-scale GWD +! gwd_opt_bl = 0 or 1: blocking drag +! +! References: +! Choi and Hong (2015) J. Geophys. Res. +! Hong et al. (2008), wea. and forecasting +! Kim and Doyle (2005), Q. J. R. Meteor. Soc. +! Kim and Arakawa (1995), j. atmos. sci. +! Alpert et al. (1988), NWP conference. +! Hong (1999), NCEP office note 424. +! Steeneveld et al (2008), JAMC +! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc. +! Beljaars et al. (2004), Q. J. R. Meteor. Soc. +! +! notice : comparible or lower resolution orography files than model resolution +! are desirable in preprocess (wps) to prevent weakening of the drag +!------------------------------------------------------------------------------- +! +! input +! dudt (im,km) non-lin tendency for u wind component +! dvdt (im,km) non-lin tendency for v wind component +! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt +! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt +! t1(im,km) temperature deg k at t0-dt +! q1(im,km) specific humidity at t0-dt +! deltim time step secs +! del(km) positive increment of pressure across layer (pa) +! KPBL(IM) is the index of the top layer of the PBL +! ipr & lprnt for diagnostics +! +! output +! dudt, dvdt wind tendency due to gwdo +! dTdt +! +!------------------------------------------------------------------------------- + +!end wrf code +!----------------------------------------------------------------------C +! USE +! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES) +! +! PURPOSE +! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH- +! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V +! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED +! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING +! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF +! CRITICAL LEVELS +! +! +! ******************************************************************** + USE MACHINE , ONLY : kind_phys + implicit none + + ! Interface variables + integer, intent(in) :: im, km, imx, kdt, ipr, me, master + integer, intent(in) :: gwd_opt + logical, intent(in) :: lprnt + integer, intent(in) :: KPBL(:) + real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:) + real(kind=kind_phys), intent(inout) :: dtend(:,:,:) + logical, intent(in) :: ldiag3d + integer, intent(in) :: dtidx(:,:), index_of_temperature, & + & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind + + integer :: kpblmax + integer, parameter :: ims=1, kms=1, its=1, kts=1 + real(kind=kind_phys), intent(in) :: fv, pi + real(kind=kind_phys) :: rcl, cdmb + real(kind=kind_phys) :: g_inv + + real(kind=kind_phys), intent(inout) :: & + & dudt(:,:),dvdt(:,:), & + & dtdt(:,:) + real(kind=kind_phys), intent(out) :: rdxzb(:) + real(kind=kind_phys), intent(in) :: & + & u1(:,:),v1(:,:), & + & t1(:,:),q1(:,:), & + & PHII(:,:),prsl(:,:), & + & prslk(:,:),PHIL(:,:) + real(kind=kind_phys), intent(in) :: prsi(:,:), & + & del(:,:) + real(kind=kind_phys), intent(in) :: var(:),oc1(:), & + & oa4(:,:),ol4(:,:), & + & dx(:) + real(kind=kind_phys), intent(in) :: varss(:),oc1ss(:), & + & oa4ss(:,:),ol4ss(:,:) + real(kind=kind_phys), intent(in) :: THETA(:),SIGMA(:), & + & GAMMA(:),ELVMAX(:) + +! added for small-scale orographic wave drag + real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx + integer, intent(in) :: vtype(:) + real(kind=kind_phys), intent(in) :: br1(:), & + & hpbl(:), & + & slmsk(:) + real(kind=kind_phys), dimension(im) :: govrth,xland + !real(kind=kind_phys), dimension(im,km) :: dz2 + real(kind=kind_phys) :: tauwavex0,tauwavey0, & + & XNBV,density,tvcon,hpbl2 + integer :: kpbl2,kvar + !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g + real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g + +!SPP + real(kind=kind_phys), dimension(im) :: var_stoch, varss_stoch, & + varmax_ss_stoch, varmax_fd_stoch + real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) + integer, intent(in) :: spp_gwd + + real(kind=kind_phys), dimension(im) :: rstoch + +!Output: + real(kind=kind_phys), intent(out) :: & + & dusfc(:), dvsfc(:) +!Output (optional): + real(kind=kind_phys), intent(out) :: & + & dusfc_ls(:),dvsfc_ls(:), & + & dusfc_bl(:),dvsfc_bl(:), & + & dusfc_ss(:),dvsfc_ss(:), & + & dusfc_fd(:),dvsfc_fd(:) + real(kind=kind_phys), intent(out) :: & + & dtaux2d_ls(:,:),dtauy2d_ls(:,:), & + & dtaux2d_bl(:,:),dtauy2d_bl(:,:), & + & dtaux2d_ss(:,:),dtauy2d_ss(:,:), & + & dtaux2d_fd(:,:),dtauy2d_fd(:,:) + +!Misc arrays + real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d + +!------------------------------------------------------------------------- +! Flags to regulate the activation of specific components of drag suite: +! Each component is tapered off automatically as a function of dx, so best to +! keep them activated (.true.). + logical, intent(in) :: & + do_gsl_drag_ls_bl, & ! large-scale gravity wave drag and blocking + do_gsl_drag_ss, & ! small-scale gravity wave drag (Steeneveld et al. 2008) + do_gsl_drag_tofd ! form drag (Beljaars et al. 2004, QJRMS) +! Flag for diagnostic outputs + logical, intent(in) :: ldiag_ugwp + +! Flag for sequential update of u and v between +! LSGWD + BLOCKING and SSGWD + TOFD calculations + logical, intent(in) :: ugwp_seq_update +! +! Additional flags + integer, parameter :: & + gwd_opt_ls = 1, & ! large-scale gravity wave drag + gwd_opt_bl = 1, & ! blocking drag + gsd_diss_ht_opt = 0 + +! Parameters for bounding the scale-adaptive variability: +! Small-scale GWD + turbulent form drag + real(kind=kind_phys), parameter :: dxmin_ss = 1000., & + & dxmax_ss = 12000. ! min,max range of tapering (m) +! Large-scale GWD + blocking + real(kind=kind_phys), parameter :: dxmin_ls = 3000., & + & dxmax_ls = 13000. ! min,max range of tapering (m) + real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-) +! +! Variables for limiting topographic standard deviation (var) + real(kind=kind_phys), parameter :: varmax_ss = 50., & + varmax_fd = 150., & + beta_ss = 0.1, & + beta_fd = 0.2 + real(kind=kind_phys) :: var_temp, var_temp2 + +! added Beljaars orographic form drag + real(kind=kind_phys), dimension(im,km) :: utendform,vtendform + real(kind=kind_phys) :: a1,a2,wsp + real(kind=kind_phys) :: H_efold + +! multification factor of standard deviation : ! larger drag with larger value +!!! real(kind=kind_phys), parameter :: psl_gwd_dx_factor = 6.0 + real(kind=kind_phys), intent(in) :: psl_gwd_dx_factor + +! critical richardson number for wave breaking : ! larger drag with larger value + real(kind=kind_phys), parameter :: ric = 0.25 + real(kind=kind_phys), parameter :: dw2min = 1. + real(kind=kind_phys), parameter :: rimin = -100. + real(kind=kind_phys), parameter :: bnv2min = 1.0e-5 + real(kind=kind_phys), parameter :: efmin = 0.0 + real(kind=kind_phys), parameter :: efmax = 10.0 + real(kind=kind_phys), parameter :: xl = 4.0e4 + real(kind=kind_phys), parameter :: critac = 1.0e-5 + real(kind=kind_phys), parameter :: gmax = 1. + real(kind=kind_phys), parameter :: veleps = 1.0 + real(kind=kind_phys), parameter :: factop = 0.5 + real(kind=kind_phys), parameter :: frc = 1.0 + real(kind=kind_phys), parameter :: ce = 0.8 + real(kind=kind_phys), parameter :: cg = 1.0 +! real(kind=kind_phys), parameter :: var_min = 100.0 + real(kind=kind_phys), parameter :: var_min = 10.0 + real(kind=kind_phys), parameter :: hmt_min = 50. + real(kind=kind_phys), parameter :: oc_min = 1.0 + real(kind=kind_phys), parameter :: oc_max = 10.0 +! 7.5 mb -- 33 km ... 0.01 kgm-3 reduce gwd drag above cutoff level + real(kind=kind_phys), parameter :: pcutoff = 7.5e2 +! 0.76 mb -- 50 km ...0.001 kgm-3 --- 0.1 mb 65 km 0.0001 kgm-3 + real(kind=kind_phys), parameter :: pcutoff_den = 0.01 ! + + integer,parameter :: kpblmin = 2 + +! +! local variables +! + integer :: i,j,k,lcap,lcapp1,nwd,idir, & + klcap,kp1 +! + real(kind=kind_phys) :: rcs,csg,fdir,cs, & + rcsks,wdir,ti,rdz,tem2,dw2,shr2, & + bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, & + rim,temc,tem1,efact,temv,dtaux,dtauy, & + dtauxb,dtauyb,eng0,eng1 + real(kind=kind_phys) :: denfac +! + logical :: ldrag(im),icrilv(im), & + flag(im) +! + real(kind=kind_phys) :: invgrcs +! + real(kind=kind_phys) :: taub(im),taup(im,km+1), & + xn(im),yn(im), & + ubar(im),vbar(im), & + fr(im),ulow(im), & + rulow(im),bnv(im), & + oa(im),ol(im),oc(im), & + oass(im),olss(im), & + roll(im),dtfac(im,km), & + brvf(im),xlinv(im), & + delks(im),delks1(im), & + bnv2(im,km),usqj(im,km), & + taud_ls(im,km),taud_bl(im,km), & + ro(im,km), & + vtk(im,km),vtj(im,km), & + zlowtop(im),velco(im,km-1), & + coefm(im),coefm_ss(im) + real(kind=kind_phys) :: cleff(im),cleff_ss(im) +! + integer :: kbl(im),klowtop(im) + integer,parameter :: mdir=8 + !integer :: nwdir(mdir) + !data nwdir/6,7,5,8,2,3,1,4/ + integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/) +! +! variables for flow-blocking drag +! + real(kind=kind_phys),parameter :: frmax = 10. + real(kind=kind_phys),parameter :: olmin = 1.e-5 + real(kind=kind_phys),parameter :: odmin = 0.1 + real(kind=kind_phys),parameter :: odmax = 10. + real(kind=kind_phys),parameter :: cdmin = 0.0 + real(kind=kind_phys),parameter :: erad = 6371.315e+3 + integer :: komax(im),kbmax(im),kblk(im) + real(kind=kind_phys) :: hmax(im) + real(kind=kind_phys) :: cd + real(kind=kind_phys) :: zblk,tautem + real(kind=kind_phys) :: pe,ke + real(kind=kind_phys) :: delx,dely + real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4) + real(kind=kind_phys) :: dxy(im),dxyp(im) + real(kind=kind_phys) :: ol4p(4),olp(im),od(im) + real(kind=kind_phys) :: taufb(im,km+1) + + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: udtend, vdtend, Tdtend + + ! Calculate inverse of gravitational acceleration + g_inv = 1./G + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Initialize local variables + var_temp2 = 0. + udtend = -1 + vdtend = -1 + Tdtend = -1 + + if(ldiag3d) then + udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd) + vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd) + Tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd) + endif +! +!---- constants +! + rcl = 1. + rcs = sqrt(rcl) + cs = 1. / sqrt(rcl) + csg = cs * g + lcap = km + lcapp1 = lcap + 1 + fdir = mdir / (2.0*pi) + invgrcs = 1._kind_phys/g*rcs + kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2 + denfac = 1.0 + + do i=1,im + if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3 + xland(i)=1.0 !but land/water = (1/2) in this module + else + xland(i)=2.0 + endif + RDXZB(i) = 0.0 + enddo + +!--- calculate scale-aware tapering factors +do i=1,im + if ( dx(i) .ge. dxmax_ls ) then + ls_taper(i) = 1. + else + if ( dx(i) .le. dxmin_ls) then + ls_taper(i) = 0. + else + ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & + (dxmax_ls-dxmin_ls)) + 1. ) + endif + endif +enddo + +! Remove ss_tapering +ss_taper(:) = 1. + +! SPP, if spp_gwd is 0, no perturbations are applied. +if ( spp_gwd==1 ) then + do i = its,im + var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1) + varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1) + varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1) + varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1) + enddo +else + do i = its,im + var_stoch(i) = var(i) + varss_stoch(i) = varss(i) + varmax_ss_stoch(i) = varmax_ss + varmax_fd_stoch(i) = varmax_fd + enddo +endif + +!--- calculate length of grid for flow-blocking drag +! +do i=1,im + delx = dx(i) + dely = dx(i) + dxy4(i,1) = delx + dxy4(i,2) = dely + dxy4(i,3) = sqrt(delx*delx + dely*dely) + dxy4(i,4) = dxy4(i,3) + dxy4p(i,1) = dxy4(i,2) + dxy4p(i,2) = dxy4(i,1) + dxy4p(i,3) = dxy4(i,4) + dxy4p(i,4) = dxy4(i,3) + cleff(i) = psl_gwd_dx_factor*(delx+dely)*0.5 + cleff_ss(i) = 0.1 * max(dxmax_ss,dxy4(i,3)) +! cleff_ss(i) = cleff(i) ! consider ..... +enddo +! +!-----initialize arrays +! + dtaux = 0.0 + dtauy = 0.0 + do i = its,im + klowtop(i) = 0 + kbl(i) = 0 + enddo +! + do i = its,im + xn(i) = 0.0 + yn(i) = 0.0 + ubar (i) = 0.0 + vbar (i) = 0.0 + roll (i) = 0.0 + taub (i) = 0.0 + oa(i) = 0.0 + ol(i) = 0.0 + oc(i) = 0.0 + oass(i) = 0.0 + olss(i) = 0.0 + ulow (i) = 0.0 + rstoch(i) = 0.0 + ldrag(i) = .false. + icrilv(i) = .false. + enddo + + do k = kts,km + do i = its,im + usqj(i,k) = 0.0 + bnv2(i,k) = 0.0 + vtj(i,k) = 0.0 + vtk(i,k) = 0.0 + taup(i,k) = 0.0 + taud_ls(i,k) = 0.0 + taud_bl(i,k) = 0.0 + dtaux2d(i,k) = 0.0 + dtauy2d(i,k) = 0.0 + dtfac(i,k) = 1.0 + enddo + enddo +! + if ( ldiag_ugwp ) then + do i = its,im + dusfc_ls(i) = 0.0 + dvsfc_ls(i) = 0.0 + dusfc_bl(i) = 0.0 + dvsfc_bl(i) = 0.0 + dusfc_ss(i) = 0.0 + dvsfc_ss(i) = 0.0 + dusfc_fd(i) = 0.0 + dvsfc_fd(i) = 0.0 + enddo + do k = kts,km + do i = its,im + dtaux2d_ls(i,k)= 0.0 + dtauy2d_ls(i,k)= 0.0 + dtaux2d_bl(i,k)= 0.0 + dtauy2d_bl(i,k)= 0.0 + dtaux2d_ss(i,k)= 0.0 + dtauy2d_ss(i,k)= 0.0 + dtaux2d_fd(i,k)= 0.0 + dtauy2d_fd(i,k)= 0.0 + enddo + enddo + endif + + do i = its,im + taup(i,km+1) = 0.0 + xlinv(i) = 1.0/xl + dusfc(i) = 0.0 + dvsfc(i) = 0.0 + enddo +! +! initialize array for flow-blocking drag +! + taufb(1:im,1:km+1) = 0.0 + hmax(1:im) = 0.0 + komax(1:im) = 0 + kbmax(1:im) = 0 + kblk(1:im) = 0 +! + do k = kts,km + do i = its,im + vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k)) + vtk(i,k) = vtj(i,k) / prslk(i,k) + ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3 + enddo + enddo +! +! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2). +! + !zq=0. + do k = kts,km + do i = its,im + !zq(i,k+1) = PHII(i,k+1)*g_inv + !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv + zl(i,k) = PHIL(i,k)*g_inv + enddo + enddo +! +! determine reference level: maximum of 2*var and pbl heights +! + do i = its,im + if(vtype(i)==15) then + zlowtop(i) = 1.0 * var_stoch(i) !!! reduce drag over land ice + else + zlowtop(i) = 2.0 * var_stoch(i) + endif + enddo +! + do i = its,im + flag(i) = .true. + enddo +! + do k = kts+1,km + do i = its,im + if(flag(i).and.zl(i,k).ge.zlowtop(i)) then + klowtop(i) = k+1 + flag(i) = .false. + endif + enddo + enddo +! +! determine the maximum height level +! note taht elvmax and zl are the heights from the model surface whereas +! oro (mean orography) is the height from the sea level +! + do i = its,im + flag(i) = .true. + enddo +! + do k = kts+1,km + do i = its,im + if(flag(i).and.zl(i,k).ge.elvmax(i)) then + komax(i) = k+1 + flag(i) = .false. + endif + enddo + enddo +! +! determine the launching level in determining blocking layer +! + do i = its,im + flag(i) = .true. + enddo +! + do k = kts+1,km + do i = its,im + if(flag(i).and.zl(i,k).ge.elvmax(i)+zlowtop(i)) then + kbmax(i) = k+1 + flag(i) = .false. + endif + enddo + enddo +! +! determing the reference level for gwd and blockding... +! + do i = its,im + hmax(i) = max(elvmax(i),zlowtop(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(min(kbl(i),kpblmax),kpblmin) + enddo +! +! compute low level averages below reference level +! + do i = its,im + delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i))) + delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i))) + enddo + do k = kts,kpblmax + do i = its,im + if (k.lt.kbl(i)) then + rcsks = rcs * del(i,k) * delks(i) + rdelks = del(i,k) * delks(i) + ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean + vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean + roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean + endif + enddo + enddo +! +! figure out low-level horizontal wind direction +! +! nwd 1 2 3 4 5 6 7 8 +! wd w s sw nw e n ne se +! + do i = its,im + wdir = atan2(ubar(i),vbar(i)) + pi + idir = mod(nint(fdir*wdir),mdir) + 1 + nwd = nwdir(idir) + oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1) + ol(i) = max(ol4(i,mod(nwd-1,4)+1),olmin) + oc(i) = min(max(oc1(i),oc_min),oc_max) +! if (var(i).le.var_min) then +! oc(i) = max(oc(i)*var(i)/var_min,oc_min) +! endif + ! Repeat for small-scale gwd + oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1) + olss(i) = ol4ss(i,mod(nwd-1,4)+1) + +! +!----- compute orographic width along (ol) and perpendicular (olp) +!----- the direction of wind +! + ol4p(1) = ol4(i,2) + ol4p(2) = ol4(i,1) + ol4p(3) = ol4(i,4) + ol4p(4) = ol4(i,3) + olp(i) = max(ol4p(mod(nwd-1,4)+1),olmin) +! +!----- compute orographic direction (horizontal orographic aspect ratio) +! + od(i) = olp(i)/ol(i) + od(i) = min(od(i),odmax) + od(i) = max(od(i),odmin) +! +!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions +! + dxy(i) = dxy4(i,MOD(nwd-1,4)+1) + dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1) + enddo +! +! END INITIALIZATION; BEGIN GWD CALCULATIONS: +! +IF ( (do_gsl_drag_ls_bl).and. & + ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + +! +!--- saving richardson number in usqj for migwdi +! + do k = kts,km-1 + ti = 2.0 / (t1(i,k)+t1(i,k+1)) + rdz = 1./(zl(i,k+1) - zl(i,k)) + tem1 = u1(i,k) - u1(i,k+1) + tem2 = v1(i,k) - v1(i,k+1) + dw2 = rcl*(tem1*tem1 + tem2*tem2) + shr2 = max(dw2,dw2min) * rdz * rdz + bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + usqj(i,k) = max(bvf2/shr2,rimin) + bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) + bnv2(i,k) = max( bnv2(i,k), bnv2min ) + enddo +! +!----compute the "low level" or 1/3 wind magnitude (m/s) +! + ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0) + rulow(i) = 1./ulow(i) +! + do k = kts,km-1 + velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) & + + (v1(i,k)+v1(i,k+1)) * vbar(i)) + velco(i,k) = velco(i,k) * rulow(i) + if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then + velco(i,k) = veleps + endif + enddo +! +! no drag when sub-oro is too small.. +! + ldrag(i) = hmax(i).le.hmt_min +! +! no drag when critical level in the base layer +! + ldrag(i) = ldrag(i).or. velco(i,1).le.0. +! +! no drag when velco.lt.0 +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0. + enddo +! +! no drag when bnv2.lt.0 +! + do k = kts,kpblmax + if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0. + enddo +! +!-----the low level weighted average ri is stored in usqj(1,1; im) +!-----the low level weighted average n**2 is stored in bnv2(1,1; im) +!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 +!---- rdelks (del(k)/delks) vert ave factor so we can * instead of / +! + wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i) + bnv2(i,1) = wtkbj * bnv2(i,1) + usqj(i,1) = wtkbj * usqj(i,1) +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) then + rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i) + bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks + usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks + endif + enddo +! + ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0 + ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0 + ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0 + ldrag(i) = ldrag(i) .or. xland(i) .gt. 1.5 +! +! set all ri low level values to the low level value +! + do k = kpblmin,kpblmax + if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1) + enddo +! + if (.not.ldrag(i)) then + bnv(i) = sqrt( bnv2(i,1) ) + fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i) + fr(i) = min(fr(i),frmax) + xn(i) = ubar(i) * rulow(i) + yn(i) = vbar(i) * rulow(i) + endif +! +! compute the base level stress and store it in taub +! calculate enhancement factor, number of mountains & aspect +! ratio const. use simplified relationship between standard +! deviation & critical hgt + + if (.not. ldrag(i)) then + efact = (oa(i) + 2.) ** (ce*fr(i)/frc) + efact = min( max(efact,efmin), efmax ) + coefm(i) = (1. + ol(i)) ** (oa(i)+1.) + xlinv(i) = coefm(i) / cleff(i) + tem = fr(i) * fr(i) * oc(i) + gfobnv = gmax * tem / ((tem + cg)*bnv(i)) + if ( gwd_opt_ls .NE. 0 ) then + taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) & + * ulow(i) * gfobnv * efact + else ! We've gotten what we need for the blocking scheme + taub(i) = 0.0 + end if + else + taub(i) = 0.0 + xn(i) = 0.0 + yn(i) = 0.0 + endif + + endif ! (ls_taper(i).GT.1.E-02) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) + +!========================================================= +! add small-scale wavedrag for stable boundary layer +!========================================================= + XNBV=0. + tauwavex0=0. + tauwavey0=0. + density=1.2 + utendwave=0. + vtendwave=0. +! +IF ( do_gsl_drag_ss ) THEN + + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + ! + ! calculating potential temperature + ! + do k = kts,km + thx(i,k) = t1(i,k)/prslk(i,k) + enddo + ! + do k = kts,km + tvcon = (1.+fv*q1(i,k)) + thvx(i,k) = thx(i,k)*tvcon + enddo + + hpbl2 = hpbl(i)+10. + kpbl2 = kpbl(i) + !kvar = MIN(kpbl, k-level of var) + kvar = 1 + do k=kts+1,MAX(kpbl(i),kts+1) +! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then + IF (zl(i,k)>300.) then + kpbl2 = k + IF (k == kpbl(i)) then + hpbl2 = hpbl(i)+10. + ELSE + hpbl2 = zl(i,k)+10. + ENDIF + exit + ENDIF + 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 + coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.) + xlinv(i) = coefm_ss(i) / cleff_ss(i) + !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)) + XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2) +! + !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then + !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) + ! 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) + tauwavex0=tauwavex0*ss_taper(i) + else + tauwavex0=0. + endif +! + !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then + if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then + !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) + ! 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) + tauwavey0=tauwavey0*ss_taper(i) + else + tauwavey0=0. + endif + + do k=kts,kpbl(i) !MIN(kpbl2+1,km-1) +!original + !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) + !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i) +!new + utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 + vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2 +!mod-to be used in HRRRv3/RAPv4 + !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2 + !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2 + enddo + endif + endif + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendwave(i,k) + dvdt(i,k) = dvdt(i,k) + vtendwave(i,k) + dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim + endif + if ( ldiag_ugwp ) then + do k = kts,km + dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k) + dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k) + dtaux2d_ss(i,k) = utendwave(i,k) + dtauy2d_ss(i,k) = vtendwave(i,k) + enddo + endif + + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im + +ENDIF ! if (do_gsl_drag_ss) + +!================================================================ +! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16): +!================================================================ +IF ( do_gsl_drag_tofd ) THEN + + do i=its,im + + if ( ss_taper(i).GT.1.E-02 ) then + + utendform=0. + vtendform=0. + + IF ((xland(i)-1.5) .le. 0.) then + !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161 + ! 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 + ! 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 + var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & + zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero + ! Note: This is a semi-implicit treatment of the time differencing + ! per Beljaars et al. (2004, QJRMS) + utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp) + vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp) + !IF(zl(i,k) > 4000.) exit + ENDDO + ENDIF + + do k = kts,km + dudt(i,k) = dudt(i,k) + utendform(i,k) + dvdt(i,k) = dvdt(i,k) + vtendform(i,k) + dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k) + dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k) + enddo + if(udtend>0) then + dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim + endif + if(vdtend>0) then + dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim + endif + if ( ldiag_ugwp ) then + do k = kts,km + dtaux2d_fd(i,k) = utendform(i,k) + dtauy2d_fd(i,k) = vtendform(i,k) + dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k) + dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k) + enddo + endif + + endif ! if (ss_taper(i).GT.1.E-02) + + enddo ! i=its,im + +ENDIF ! if (do_gsl_drag_tofd) +!======================================================= +! More for the large-scale gwd component +IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + +! +! now compute vertical structure of the stress. + do k = kts,kpblmax + if (k .le. kbl(i)) taup(i,k) = taub(i) + enddo +! + do k = kpblmin, km-1 ! vertical level k loop! + kp1 = k + 1 +! +! unstablelayer if ri < ric +! unstable layer if upper air vel comp along surf vel <=0 (crit lay) +! at (u-c)=0. crit layer exists and bit vector should be set (.le.) +! + if (k .ge. kbl(i)) then + icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) & + .or. (velco(i,k) .le. 0.0) + brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared + brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency + endif +! + if (k .ge. kbl(i) .and. (.not. ldrag(i))) then + if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then + temv = 1.0 / velco(i,k) + tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* & + velco(i,k)*0.5 + hd = sqrt(taup(i,k) / tem1) + fro = brvf(i) * hd * temv +! +! rim is the minimum-richardson number by shutts (1985) + tem2 = sqrt(usqj(i,k)) + tem = 1. + tem2 * fro + rim = usqj(i,k) * (1.-fro) / (tem * tem) +! +! check stability to employ the 'saturation hypothesis' +! of lindzen (1981) except at tropospheric downstream regions +! + if (rim .le. ric) then ! saturation hypothesis! + temc = 2.0 + 1.0 / tem2 + hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i) + taup(i,kp1) = tem1 * hd * hd + else ! no wavebreaking! + taup(i,kp1) = taup(i,k) + endif + endif + endif + enddo +! + if(lcap.lt.km) then + do klcap = lcapp1,km + taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap) + enddo + endif + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) +!=============================================================== +!COMPUTE BLOCKING COMPONENT +!=============================================================== +IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN + do i = its,im + flag(i) = .true. + enddo + + do i=its,im + + if ( ls_taper(i).GT.1.E-02 ) then + + if (.not.ldrag(i)) then +! +!------- determine the height of flow-blocking layer +! + pe = 0.0 + ke = 0.0 + do k = km, kpblmin, -1 + if(flag(i).and. k.le.kbmax(i)) then + pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* & + del(i,k)/g/ro(i,k) + ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) +! +!---------- apply flow-blocking drag when pe >= ke +! + if(pe.ge.ke.and.zl(i,k).le.hmax(i)) then + kblk(i)= k + zblk = zl(i,k) + RDXZB(i) = real(k,kind=kind_phys) + flag(i) = .false. + endif + endif + enddo + if(.not.flag(i)) then +! +!--------- compute flow-blocking stress +! + cd = max(2.0-1.0/od(i),cdmin) + taufb(i,kts) = 0.5 * roll(i) * coefm(i) / & + max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * & + olp(i) * zblk * ulow(i)**2 + tautem = taufb(i,kts)/float(kblk(i)-kts) + do k = kts+1, kpblmax + if (k .le. kblk(i)) taufb(i,k) = taufb(i,k-1) - tautem + enddo +! +! reset gwd stress below blocking layer +! + do k = kts,kpblmax + if (k .le. kblk(i)) taup(i,k) = taup(i,kblk(i)) + enddo +! if(kblk(i).gt.5) print *,' gwd kbl komax kbmax kblk ',kbl(i),komax(i),kbmax(i),kblk(i) +! if(kblk(i).gt.5) print *,' gwd elvmax zlowtop zblk ',elvmax(i),zlowtop(i),zl(i,kblk(i)) + endif + + endif ! if (.not.ldrag(i)) + + endif ! if ( ls_taper(i).GT.1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) +!=========================================================== +IF ( (do_gsl_drag_ls_bl) .and. & + (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN + + do i=its,im + + if ( ls_taper(i) .GT. 1.E-02 ) then + +! +! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy +! + do k = kts,km + taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k) + taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k) + enddo +! +! limit de-acceleration (momentum deposition ) at top to 1/2 value +! the idea is some stuff must go out the 'top' + do klcap = lcap,km + taud_ls(i,klcap) = taud_ls(i,klcap) * factop + taud_bl(i,klcap) = taud_bl(i,klcap) * factop + enddo +! +! if the gravity wave drag would force a critical line +! in the lower ksmm1 layers during the next deltim timestep, +! then only apply drag until that critical line is reached. +! + do k = kts,kpblmax-1 + if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) then + dtfac(i,k) = min(dtfac(i,k),abs(velco(i,k) & + /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) + endif + enddo +! apply limiter to mesosphere drag, reduce the drag by density factor 10-3 +! prevent wind reversal... +! + do k = kpblmax,km + if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff) then + denfac = min(ro(i,k)/pcutoff_den,1.) + dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k) & + /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k))))) + endif + enddo +! + do k = kts,km + taud_ls(i,k) = taud_ls(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i)) + taud_bl(i,k) = taud_bl(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i)) + dtaux = taud_ls(i,k) * xn(i) + dtauy = taud_ls(i,k) * yn(i) + dtauxb = taud_bl(i,k) * xn(i) + dtauyb = taud_bl(i,k) * yn(i) + + !add blocking and large-scale contributions to tendencies + dudt(i,k) = dtaux + dtauxb + dudt(i,k) + dvdt(i,k) = dtauy + dtauyb + dvdt(i,k) + + if ( gsd_diss_ht_opt .EQ. 1 ) then + ! Calculate dissipation heating + ! Initial kinetic energy (at t0-dt) + eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. ) + ! Kinetic energy after wave-breaking/flow-blocking + eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + & + (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 ) + ! Modify theta tendency + dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim + if ( Tdtend>0 ) then + dtend(i,k,Tdtend) = dtend(i,k,Tdtend) + max((eng0-eng1),0.0)/cp + endif + endif + + dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + & + taud_bl(i,k)*xn(i)*del(i,k) + dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + & + taud_bl(i,k)*yn(i)*del(i,k) + if(udtend>0) then + dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * & + xn(i) + taud_bl(i,k) * xn(i)) * deltim + endif + if(vdtend>0) then + dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * & + yn(i) + taud_bl(i,k) * yn(i)) * deltim + endif + + enddo + + ! Finalize dusfc and dvsfc diagnostics + dusfc(i) = -(invgrcs) * dusfc(i) + dvsfc(i) = -(invgrcs) * dvsfc(i) + + if ( ldiag_ugwp ) then + do k = kts,km + dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i) + dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i) + dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i) + dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i) + dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k) + dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k) + dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k) + dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k) + enddo + endif + + endif ! if ( ls_taper(i) .GT. 1.E-02 ) + + enddo ! do i=its,im + +ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1) + +if ( ldiag_ugwp ) then + ! Finalize dusfc and dvsfc diagnostics + do i = its,im + dusfc_ls(i) = -(invgrcs) * dusfc_ls(i) + dvsfc_ls(i) = -(invgrcs) * dvsfc_ls(i) + dusfc_bl(i) = -(invgrcs) * dusfc_bl(i) + dvsfc_bl(i) = -(invgrcs) * dvsfc_bl(i) + dusfc_ss(i) = -(invgrcs) * dusfc_ss(i) + dvsfc_ss(i) = -(invgrcs) * dvsfc_ss(i) + dusfc_fd(i) = -(invgrcs) * dusfc_fd(i) + dvsfc_fd(i) = -(invgrcs) * dvsfc_fd(i) + enddo +endif +! + return + end subroutine drag_suite_psl ! !> @} diff --git a/physics/GWD/drag_suite.meta b/physics/GWD/drag_suite.meta index 94dddcc93..6981d4761 100644 --- a/physics/GWD/drag_suite.meta +++ b/physics/GWD/drag_suite.meta @@ -438,6 +438,13 @@ type = real kind = kind_phys intent = in +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in [g] standard_name = gravitational_acceleration long_name = gravitational acceleration @@ -573,6 +580,14 @@ dimensions = () type = logical intent = in +[psl_gwd_dx_factor] + standard_name = effective_grid_spacing_of_psl_gwd_suite + long_name = multiplication of grid spacing + units = count + dimensions = () + type = real + kind = kind_phys + intent = in [dtend] standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables diff --git a/physics/GWD/ugwpv1_gsldrag.F90 b/physics/GWD/ugwpv1_gsldrag.F90 index 9303e0221..233377f29 100644 --- a/physics/GWD/ugwpv1_gsldrag.F90 +++ b/physics/GWD/ugwpv1_gsldrag.F90 @@ -44,7 +44,7 @@ module ugwpv1_gsldrag use cires_ugwpv1_solv2, only: cires_ugwpv1_ngw_solv2 use cires_ugwpv1_oro, only: orogw_v1 - use drag_suite, only: drag_suite_run + use drag_suite, only: drag_suite_run, drag_suite_psl implicit none @@ -305,11 +305,13 @@ end subroutine ugwpv1_gsldrag_finalize !! @{ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, & fhzero, kdt, ldiag3d, lssav, flag_for_gwd_generic_tend, do_gsl_drag_ls_bl, & - do_gsl_drag_ss, do_gsl_drag_tofd, do_ugwp_v1, do_ugwp_v1_orog_only, & + do_gsl_drag_ss, do_gsl_drag_tofd, & + do_gwd_opt_psl, psl_gwd_dx_factor, & + do_ugwp_v1, do_ugwp_v1_orog_only, & do_ugwp_v1_w_gsldrag, gwd_opt, do_tofd, ldiag_ugwp, ugwp_seq_update, & cdmbgwd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, & elvmax, clx, oa4, varss,oc1ss,oa4ss,ol4ss, dx, xlat, xlat_d, sinlat, coslat, & - area, rain, br1, hpbl, kpbl, slmsk, & + area, rain, br1, hpbl,vtype, kpbl, slmsk, & ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, del, tau_amf, & dudt_ogw, dvdt_ogw, du_ogwcol, dv_ogwcol, & dudt_obl, dvdt_obl, du_oblcol, dv_oblcol, & @@ -367,7 +369,9 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, real(kind=kind_phys), intent(in) :: dtp, fhzero real(kind=kind_phys), intent(in) :: ak(:), bk(:) integer, intent(in) :: kdt, jdat(:) - +! option for psl gwd + logical, intent(in) :: do_gwd_opt_psl ! option for psl gravity wave drag + real(kind=kind_phys), intent(in) :: psl_gwd_dx_factor ! ! SSO parameters and variables integer, intent(in) :: gwd_opt !gwd_opt and nmtvr are "redundant" controls integer, intent(in) :: nmtvr @@ -396,6 +400,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, real(kind=kind_phys), intent(in), dimension(:,:) :: prsi, phii real(kind=kind_phys), intent(in), dimension(:,:) :: q1 integer, intent(in), dimension(:) :: kpbl + integer, intent(in), dimension(:) :: vtype real(kind=kind_phys), intent(in), dimension(:) :: rain real(kind=kind_phys), intent(in), dimension(:) :: br1, hpbl, slmsk @@ -544,6 +549,27 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, ! dusfcg, dvsfcg ! ! + if (do_gwd_opt_psl) then + call drag_suite_psl(im, levs, Pdvdt, Pdudt, Pdtdt, & + ugrs,vgrs,tgrs,q1, & + kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & + kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, & + ol4ss,theta,sigma,gamma,elvmax, & + dudt_ogw, dvdt_ogw, dudt_obl, dvdt_obl, & + dudt_oss, dvdt_oss, dudt_ofd, dvdt_ofd, & + dusfcg, dvsfcg, & + du_ogwcol, dv_ogwcol, du_oblcol, dv_oblcol, & + du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, & + slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, & + con_fv, con_pi, lonr, & + cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & + psl_gwd_dx_factor, & + dtend, dtidx, index_of_process_orographic_gwd, & + index_of_temperature, index_of_x_wind, & + index_of_y_wind, ldiag3d, ldiag_ugwp, & + ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) + else call drag_suite_run(im, levs, Pdvdt, Pdudt, Pdtdt, & ugrs,vgrs,tgrs,q1, & kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & @@ -554,7 +580,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, dusfcg, dvsfcg, & du_ogwcol, dv_ogwcol, du_oblcol, dv_oblcol, & du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, & - slmsk,br1,hpbl, con_g,con_cp,con_rd,con_rv, & + slmsk,br1,hpbl,con_g,con_cp,con_rd,con_rv, & con_fv, con_pi, lonr, & cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & @@ -562,6 +588,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, index_of_temperature, index_of_x_wind, & index_of_y_wind, ldiag3d, ldiag_ugwp, & ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) + endif ! ! dusfcg = du_ogwcol + du_oblcol + du_osscol + du_ofdcol ! diff --git a/physics/GWD/ugwpv1_gsldrag.meta b/physics/GWD/ugwpv1_gsldrag.meta index 73d7eee1c..76e0c4d47 100644 --- a/physics/GWD/ugwpv1_gsldrag.meta +++ b/physics/GWD/ugwpv1_gsldrag.meta @@ -402,6 +402,21 @@ dimensions = () type = logical intent = in +[do_gwd_opt_psl] + standard_name = do_gsl_drag_suite_with_psl_gwd_option + long_name = flag to activate PSL drag suite - mesoscale GWD and blocking + units = flag + dimensions = () + type = logical + intent = in +[psl_gwd_dx_factor] + standard_name = effective_grid_spacing_of_psl_gwd_suite + long_name = multiplication of grid spacing + units = count + dimensions = () + type = real + kind = kind_phys + intent = in [do_ugwp_v1] standard_name = flag_for_ugwp_version_1 long_name = flag to activate ver 1 CIRES UGWP @@ -641,6 +656,13 @@ type = real kind = kind_phys intent = in +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in [kpbl] standard_name = vertical_index_at_top_of_atmosphere_boundary_layer long_name = vertical index at top atmospheric boundary layer diff --git a/physics/GWD/unified_ugwp.F90 b/physics/GWD/unified_ugwp.F90 index fcdee3b5d..a2f56c0e4 100644 --- a/physics/GWD/unified_ugwp.F90 +++ b/physics/GWD/unified_ugwp.F90 @@ -40,7 +40,7 @@ module unified_ugwp use gwdps, only: gwdps_run use cires_ugwp_triggers use ugwp_driver_v0 - use drag_suite, only: drag_suite_run + use drag_suite, only: drag_suite_run, drag_suite_psl implicit none @@ -249,11 +249,11 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt varss,oc1ss,oa4ss,ol4ss,dx,dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl,dusfc_ss, & dvsfc_ss,dusfc_fd,dvsfc_fd,dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,dtauy2d_bl, & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd,dudt_ngw,dvdt_ngw,dtdt_ngw, & - br1,hpbl,slmsk, do_tofd, ldiag_ugwp, ugwp_seq_update, & + br1,hpbl,vtype,slmsk, do_tofd, ldiag_ugwp, ugwp_seq_update, & cdmbgwd, jdat, xlat, xlat_d, sinlat, coslat, area, & ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, & del, kpbl, dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, & - tau_tofd, tau_mtb, tau_ogw, tau_ngw, & + tau_tofd, tau_mtb, tau_ogw, tau_ngw, zmtb, zlwb, zogw, & dudt_mtb, dudt_tms, du3dt_mtb, du3dt_ogw, du3dt_tms, & dudt, dvdt, dtdt, rdxzb, con_g, con_omega, con_pi, con_cp, con_rd, con_rv, & con_rerth, con_fvirt, rain, ntke, q_tke, dqdt_tke, lprnt, ipr, & @@ -262,6 +262,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt index_of_process_nonorographic_gwd, & lssav, flag_for_gwd_generic_tend, do_ugwp_v0, do_ugwp_v0_orog_only, & do_ugwp_v0_nst_only, do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & + do_gwd_opt_psl, psl_gwd_dx_factor, & gwd_opt, spp_wts_gwd, spp_gwd, errmsg, errflg) implicit none @@ -270,6 +271,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt integer, intent(in) :: me, master, im, levs, ntrac, kdt, lonr, nmtvr integer, intent(in) :: gwd_opt integer, intent(in), dimension(:) :: kpbl + integer, intent(in), dimension(:) :: vtype real(kind=kind_phys), intent(in), dimension(:) :: ak, bk real(kind=kind_phys), intent(in), dimension(:) :: oro, oro_uf, hprime, oc, theta, sigma, gamma real(kind=kind_phys), intent(in), dimension(:) :: varss,oc1ss, dx @@ -309,7 +311,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt & slmsk(:) real(kind=kind_phys), intent(out), dimension(:) :: dusfcg, dvsfcg - real(kind=kind_phys), intent(out), dimension(:) :: rdxzb + real(kind=kind_phys), intent(out), dimension(:) :: zmtb, zlwb, zogw, rdxzb real(kind=kind_phys), intent(out), dimension(:) :: tau_mtb, tau_ogw, tau_tofd, tau_ngw real(kind=kind_phys), intent(out), dimension(:,:) :: gw_dudt, gw_dvdt, gw_dtdt, gw_kdis real(kind=kind_phys), intent(out), dimension(:,:) :: dudt_mtb, dudt_tms @@ -345,6 +347,10 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) integer, intent(in) :: spp_gwd + ! option for psl gwd + logical, intent(in) :: do_gwd_opt_psl ! option for psl gravity wave drag + real(kind=kind_phys), intent(in) :: psl_gwd_dx_factor ! + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -487,7 +493,26 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt ! Note: In case of GSL drag_suite, this includes ss and tofd if ( do_gsl_drag_ls_bl.or.do_gsl_drag_ss.or.do_gsl_drag_tofd ) then - +! +if (do_gwd_opt_psl) then + call drag_suite_psl(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, & + tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & + kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, & + ol4ss,theta,sigma,gamma,elvmax,dtaux2d_ms, & + dtauy2d_ms,dtaux2d_bl,dtauy2d_bl,dtaux2d_ss, & + dtauy2d_ss,dtaux2d_fd,dtauy2d_fd,dusfcg, & + dvsfcg,dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl, & + dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & + slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, & + con_fvirt,con_pi,lonr, & + cdmbgwd,me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & + psl_gwd_dx_factor, & + dtend, dtidx, index_of_process_orographic_gwd, & + index_of_temperature, index_of_x_wind, & + index_of_y_wind, ldiag3d, ldiag_ugwp, & + ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) +else call drag_suite_run(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, & tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, & @@ -504,6 +529,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt index_of_temperature, index_of_x_wind, & index_of_y_wind, ldiag3d, ldiag_ugwp, & ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) +endif ! ! put zeros due to xy GSL-drag style: dtaux2d_bl,dtauy2d_bl,dtaux2d_ss.......dusfc_ms,dvsfc_ms ! diff --git a/physics/GWD/unified_ugwp.meta b/physics/GWD/unified_ugwp.meta index 189f7072c..e0d60459f 100644 --- a/physics/GWD/unified_ugwp.meta +++ b/physics/GWD/unified_ugwp.meta @@ -648,6 +648,13 @@ type = real kind = kind_phys intent = in +[vtype] + standard_name = vegetation_type_classification + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_loop_extent) + type = integer + intent = in [slmsk] standard_name = area_type long_name = landmask: sea/land/ice=0/1/2 @@ -900,6 +907,30 @@ type = real kind = kind_phys intent = out +[zmtb] + standard_name = height_of_mountain_blocking + long_name = height of mountain blocking drag + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[zlwb] + standard_name = height_of_low_level_wave_breaking + long_name = height of low level wave breaking + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[zogw] + standard_name = height_of_launch_level_of_orographic_gravity_wave + long_name = height of launch level of orographic gravity wave + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [dudt_mtb] standard_name = instantaneous_change_in_x_wind_due_to_mountain_blocking_drag long_name = instantaneous change in x wind due to mountain blocking drag @@ -1194,6 +1225,21 @@ dimensions = () type = logical intent = in +[do_gwd_opt_psl] + standard_name = do_gsl_drag_suite_with_psl_gwd_option + long_name = flag to activate PSL drag suite - mesoscale GWD and blocking + units = flag + dimensions = () + type = logical + intent = in +[psl_gwd_dx_factor] + standard_name = effective_grid_spacing_of_psl_gwd_suite + long_name = multiplication of grid spacing + units = count + dimensions = () + type = real + kind = kind_phys + intent = in [gwd_opt] standard_name = control_for_drag_suite_gravity_wave_drag long_name = flag to choose gwd scheme diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 2472f11ba..38f88f487 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -423,8 +423,9 @@ subroutine noahmp_sflx (parameters, & smceq , & ! in : vegetation/soil characteristics sfctmp , sfcprs , psfc , uu , vv , q2, garea1 , & ! in : forcing qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing + varf , psl_gwd_z0m_factor, & ! in : small scale oro std pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -436,7 +437,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -445,9 +446,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP @@ -493,6 +494,8 @@ subroutine noahmp_sflx (parameters, & real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< in exner function + real (kind=kind_phys) , intent(in) :: varf !< small scale oro std + real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! real (kind=kind_phys) , intent(in) :: pblhx !< pbl height integer , intent(in) :: iz0tlnd !< z0t option integer , intent(in) :: itime !< @@ -832,6 +835,7 @@ subroutine noahmp_sflx (parameters, & fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in + varf,psl_gwd_z0m_factor , & !in pblhx ,iz0tlnd, itime ,psi_opt, ep_1, ep_2, epsm1,cp, & z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out @@ -1676,6 +1680,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in + varf,psl_gwd_z0m_factor , & !in pblhx , iz0tlnd, itime,psi_opt,ep_1, ep_2, epsm1, cp, & z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out @@ -1757,6 +1762,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< + real (kind=kind_phys) , intent(in) :: varf + real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -2252,6 +2259,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in + varf , psl_gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tgv ,cmv, ustarx , & !inout #ifdef CCPP @@ -2266,7 +2274,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in !jref:start qc ,qsfc ,psfc , & !in q2v ,chv2 ,chleaf ,chuc , & - rb) !out + rb) !out ! new coupling code @@ -2290,6 +2298,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in emg ,stc ,df ,rsurf ,latheag , & !in gammag ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in + varf ,psl_gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & #ifdef CCPP tgb ,cmb ,chb, ustarx,errmsg ,errflg , & !inout @@ -3703,6 +3712,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in + varf ,psl_gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tg ,cm,ustarx,& !inout #ifdef CCPP @@ -3716,7 +3726,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & csigmaf1, & !out qc ,qsfc ,psfc , & !in q2v ,cah2 ,chleaf ,chuc , & !inout - rb) !out + rb) !out ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3753,6 +3763,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: dt !< time step (s) real (kind=kind_phys), intent(in) :: fsno !< snow fraction + real (kind=kind_phys) , intent(in) :: varf ! + real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -4125,6 +4137,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd ,snowh ,shdfac ,garea1 , & !in + varf,psl_gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -4434,6 +4447,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & emg ,stc ,df ,rsurf ,lathea , & !in gamma ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in + varf,psl_gwd_z0m_factor, & !in pblhx , iz0tlnd , itime ,psi_opt,ep_1,ep_2,epsm1,cp ,& #ifdef CCPP tgb ,cm ,ch,ustarx,errmsg ,errflg , & !inout @@ -4488,6 +4502,8 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(in) :: rhsur !< raltive humidity in surface soil/snow air space (-) real (kind=kind_phys), intent(in) :: fsno !< snow fraction + real (kind=kind_phys), intent(in) :: varf ! + real (kind=kind_phys), intent(in) :: psl_gwd_z0m_factor ! real (kind=kind_phys), intent(in) :: pblhx !< pbl height (m) real (kind=kind_phys), intent(in) :: ep_1 !< real (kind=kind_phys), intent(in) :: ep_2 !< @@ -4731,6 +4747,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd,snowh ,shdfac ,garea1 , & !in + varf,psl_gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5437,6 +5454,7 @@ end subroutine sfcdif2 subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h,zpd ,snowh ,fveg ,garea1 , & !in + varf,psl_gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5466,6 +5484,8 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: snowh !< snow depth [m] real (kind=kind_phys), intent(in ) :: fveg !< fractional vegetation cover real (kind=kind_phys), intent(in ) :: garea1 !< grid area [km2] + real (kind=kind_phys), intent(in ) :: varf ! standard deviation org [m] + real (kind=kind_phys), intent(in ) :: psl_gwd_z0m_factor ! z0m factor in psl tofd real (kind=kind_phys), intent(inout) :: ustarx !< friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm !< momentum stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fh !< sen heat stability correction, weighted by prior iters @@ -5520,7 +5540,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur endif call gfs_stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & - rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) + varf,psl_gwd_z0m_factor,rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) end subroutine sfcdif3 @@ -5530,6 +5550,7 @@ subroutine gfs_stability & ! --- inputs: ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & thsfc_loc, & + varf,psl_gwd_z0m_factor, & ! --- outputs: rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) @@ -5549,6 +5570,8 @@ subroutine gfs_stability & real(kind=kind_phys), intent(in) :: ztmax ! thermal roughness length real(kind=kind_phys), intent(in) :: tvs ! surface virtual temperature real(kind=kind_phys), intent(in) :: grav ! local gravity +real(kind=kind_phys), intent(in) :: varf ! turbulent scale oro std +real(kind=kind_phys), intent(in) :: psl_gwd_z0m_factor ! tofd factor logical, intent(in) :: thsfc_loc ! use local theta reference flag real(kind=kind_phys), intent(out) :: rb ! bulk richardson number [-] @@ -5607,6 +5630,11 @@ subroutine gfs_stability & real(kind=kind_phys) :: zolmax real(kind=kind_phys) xkzo +! tofd + real(kind=kind_phys) :: cf, zf, fri, fphim + real(kind=kind_phys), parameter :: varf_min = 50.0 + real(kind=kind_phys), parameter :: varf_max = 500.0 +!!! real(kind=kind_phys), parameter :: psl_gwd_z0m_factor = 0.003 z1i = one / z1 ! inverse of model height @@ -5727,6 +5755,17 @@ subroutine gfs_stability & endif ! end of if (dtv >= 0 ) then loop ! +! form drag coefficient +! + if ( varf.gt.varf_min .and. psl_gwd_z0m_factor.gt.0.0 ) then + zf = min( min(varf,varf_max)*psl_gwd_z0m_factor,z1 ) + fri = min( max( 1.-rb,0. ), 1.) + fphim = log( ( z1 + zf) / zf ) + cf = ca*ca / (fphim*fphim) * fri + else + cf = 0. + endif +! ! finish the exchange coefficient computation to provide fm and fh ! fm = fm - pm ! phi_m @@ -5738,7 +5777,7 @@ subroutine gfs_stability & tem1 = 0.00001_kind_phys/z1 ! minimum exhange coef (?) cm = max(cm, tem1) ch = max(ch, tem1) - stress = cm * wind * wind ! surface stress = Cm*U*U + stress = (cm + cf)* wind * wind ! cf for tofd 10.15) or (10.19) in Arya ustar = sqrt(stress) ! friction velocity return @@ -5887,7 +5926,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, else z0h_out = z0m_out * 0.01 endif - + elseif (opt_trs == blumel99) then reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index db3ae472f..90cc14a8d 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -137,6 +137,7 @@ subroutine noahmpdrv_run & rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, rhonewsn1,& con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1, & con_fvirt, con_rd, con_hfus, thsfc_loc, cpllnd, cpllnd2atm,& + nmtvr,mntvar,psl_gwd_z0m_factor, & ! --- in/outs: weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & @@ -287,6 +288,7 @@ subroutine noahmpdrv_run & integer , intent(in) :: iopt_stc ! option for snow/soil temperature time scheme (only layer 1) integer , intent(in) :: iopt_trs ! option for thermal roughness scheme integer , intent(in) :: iopt_diag ! option for surface diagnose approach + integer , intent(in) :: nmtvr ! number of subgrid scale oroghic data real(kind=kind_phys), dimension(:) , intent(in) :: xlatin ! latitude real(kind=kind_phys), dimension(:) , intent(in) :: xcoszin ! cosine of zenith angle integer , intent(in) :: iyrlen ! year length [days] @@ -297,6 +299,8 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(in) :: snow_mp ! microphysics snow [mm] real(kind=kind_phys), dimension(:) , intent(in) :: graupel_mp ! microphysics graupel [mm] real(kind=kind_phys), dimension(:) , intent(in) :: ice_mp ! microphysics ice/hail [mm] + real(kind=kind_phys), dimension(:,:) , intent(in) :: mntvar ! subgrid orographic statistics data + real(kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! psl tofd zom factor real(kind=kind_phys), dimension(:) , intent(in) :: rhonewsn1 ! precipitation ice density (kg/m^3) real(kind=kind_phys) , intent(in) :: con_hvap ! latent heat condensation [J/kg] real(kind=kind_phys) , intent(in) :: con_cp ! specific heat air [J/kg/K] @@ -625,7 +629,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: canopy_heat_storage ! out | within-canopy heat [W/m2] real (kind=kind_phys) :: spec_humid_sfc_veg ! out | surface specific humidty over vegetation [kg/kg] real (kind=kind_phys) :: spec_humid_sfc_bare ! out | surface specific humidty over bare soil [kg/kg] - + real (kind=kind_phys) :: ustarx ! inout |surface friction velocity real (kind=kind_phys) :: prslkix ! in exner function real (kind=kind_phys) :: prsik1x ! in exner function @@ -635,6 +639,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: cq2 real (kind=kind_phys) :: qfx real (kind=kind_phys) :: wspd1 ! wind speed with all components + real (kind=kind_phys) :: varf_grid ! std deviation orography real (kind=kind_phys) :: pblhx ! height of pbl integer :: mnice @@ -740,6 +745,7 @@ subroutine noahmpdrv_run & area_grid = garea(i) pblhx = pblh(i) + varf_grid = mntvar(i,15) prslkix = prslki(i) prsik1x = prsik1(i) @@ -971,6 +977,7 @@ subroutine noahmpdrv_run & spec_humidity_forcing ,area_grid ,cloud_water_forcing , & sw_radiation_forcing ,radiation_lw_forcing ,thsfc_loc , & prslkix ,prsik1x ,prslk1x , & + varf_grid,psl_gwd_z0m_factor , & pblhx ,iz0tlnd ,itime , & psi_opt , & precip_convective , & @@ -1062,7 +1069,7 @@ subroutine noahmpdrv_run & chxy (i) = ch_noahmp zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total - + !LAI-scale canopy resistance based on weighted sunlit shaded fraction if(rs_sunlit .le. 0.0 .or. rs_shaded .le. 0.0 .or. & lai_sunlit .eq. 0.0 .or. lai_shaded .eq. 0.0) then @@ -1072,7 +1079,7 @@ subroutine noahmpdrv_run & ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance end if - + smc (i,:) = soil_moisture_vol slc (i,:) = soil_liquid_vol snowxy (i) = float(snow_levels) @@ -1215,6 +1222,7 @@ subroutine noahmpdrv_run & call gfs_stability & (zf(i), zvfun(i), gdx, virtual_temperature, vptemp,wind(i), z0_total, z0h_total, & tvs1, con_g, thsfc_loc, & + varf_grid, psl_gwd_z0m_factor, & rb1(i), fm1(i), fh1(i), fm101(i), fh21(i), cm(i), ch(i), stress1(i), ustar1(i)) rmol1(i) = undefined !not used in GFS sfcdif -> to satsify output diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta index 7fa13f3af..3e24624fa 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta @@ -649,6 +649,29 @@ dimensions = () type = logical intent = in +[nmtvr] + standard_name = number_of_statistical_measures_of_subgrid_orography + long_name = number of statistical measures of subgrid height_above_mean_sea_level + units = count + dimensions = () + type = integer + intent = in +[mntvar] + standard_name = statistical_measures_of_subgrid_orography_collection_array + long_name = array of statistical measures of subgrid height_above_mean_sea_level + units = mixed + dimensions = (horizontal_loop_extent,number_of_statistical_measures_of_subgrid_orography) + type = real + kind = kind_phys + intent = in +[psl_gwd_z0m_factor] + standard_name = momentum_roughness_factor_in_turbulent_form_drag + long_name = multiplication of Z0m in PSL turbulent form drag + units = count + dimensions = () + type = real + kind = kind_phys + intent = in [weasd] standard_name = water_equivalent_accumulated_snow_depth_over_land long_name = water equiv of acc snow depth over land From 3342a4545eaabb3e5997dd78b693c40a961c3d0c Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Sun, 19 May 2024 15:55:12 -0400 Subject: [PATCH 02/12] Change code format in file module_sf_noahmplsm.F90 --- physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 38f88f487..8086831e9 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -425,7 +425,7 @@ subroutine noahmp_sflx (parameters, & qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing varf , psl_gwd_z0m_factor, & ! in : small scale oro std pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -437,7 +437,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -446,9 +446,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP From e142700c7f015d1bd4e52ffb26b41cc8a538fe1b Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Wed, 22 May 2024 09:18:05 -0400 Subject: [PATCH 03/12] update files noahmpdrv.F90 noahmpdrv.meta module_sf_noahmplsm.F90 --- .../Land/Noahmp/module_sf_noahmplsm.F90 | 93 +++++++++---------- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 29 ++---- physics/SFC_Models/Land/Noahmp/noahmpdrv.meta | 19 +--- 3 files changed, 53 insertions(+), 88 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 8086831e9..9e9dfed91 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -423,9 +423,9 @@ subroutine noahmp_sflx (parameters, & smceq , & ! in : vegetation/soil characteristics sfctmp , sfcprs , psfc , uu , vv , q2, garea1 , & ! in : forcing qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing - varf , psl_gwd_z0m_factor, & ! in : small scale oro std + varf , gwd_z0m_factor, & ! in : turbulent orographic form drag pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -437,7 +437,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -446,9 +446,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP @@ -494,8 +494,8 @@ subroutine noahmp_sflx (parameters, & real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< in exner function - real (kind=kind_phys) , intent(in) :: varf !< small scale oro std - real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! + real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] + real (kind=kind_phys) , intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height integer , intent(in) :: iz0tlnd !< z0t option integer , intent(in) :: itime !< @@ -835,7 +835,7 @@ subroutine noahmp_sflx (parameters, & fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in - varf,psl_gwd_z0m_factor , & !in + varf ,gwd_z0m_factor , & !in pblhx ,iz0tlnd, itime ,psi_opt, ep_1, ep_2, epsm1,cp, & z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out @@ -1680,7 +1680,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in - varf,psl_gwd_z0m_factor , & !in + varf,gwd_z0m_factor , & !in pblhx , iz0tlnd, itime,psi_opt,ep_1, ep_2, epsm1, cp, & z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out @@ -1762,8 +1762,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< - real (kind=kind_phys) , intent(in) :: varf - real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor + real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] + real (kind=kind_phys) , intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -2020,8 +2020,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in chuc = 0. chv2 = 0. rb = 0. - laisun = 0. - laisha = 0. cdmnv = 0.0 ezpdv = 0.0 @@ -2259,7 +2257,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in - varf , psl_gwd_z0m_factor, & !in + varf , gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tgv ,cmv, ustarx , & !inout #ifdef CCPP @@ -2273,8 +2271,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in csigmaf1, & !out !jref:start qc ,qsfc ,psfc , & !in - q2v ,chv2 ,chleaf ,chuc , & - rb) !out + q2v ,chv2, chleaf, chuc) !inout ! new coupling code @@ -2298,7 +2295,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in emg ,stc ,df ,rsurf ,latheag , & !in gammag ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in - varf ,psl_gwd_z0m_factor, & !in + varf ,gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & #ifdef CCPP tgb ,cmb ,chb, ustarx,errmsg ,errflg , & !inout @@ -3712,7 +3709,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in - varf ,psl_gwd_z0m_factor, & !in + varf ,gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tg ,cm,ustarx,& !inout #ifdef CCPP @@ -3725,8 +3722,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & t2mv ,psnsun ,psnsha ,canhs , & !out csigmaf1, & !out qc ,qsfc ,psfc , & !in - q2v ,cah2 ,chleaf ,chuc , & !inout - rb) !out + q2v ,cah2 ,chleaf ,chuc ) !inout ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3763,8 +3759,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: dt !< time step (s) real (kind=kind_phys), intent(in) :: fsno !< snow fraction - real (kind=kind_phys) , intent(in) :: varf ! - real (kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! + real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] + real (kind=kind_phys) , intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -3852,7 +3848,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(out) :: chuc !< under canopy exchange coefficient real (kind=kind_phys), intent(out) :: canhs !< canopy heat storage change (w/m2) real (kind=kind_phys), intent(out) :: q2v !< - real (kind=kind_phys), intent(out) :: rb !< bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: cah !< sensible heat conductance, canopy air to zlvl air (m/s) real (kind=kind_phys) :: u10v !< 10 m wind speed in eastward dir (m/s) real (kind=kind_phys) :: v10v !< 10 m wind speed in eastward dir (m/s) @@ -3869,6 +3864,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: z0mo !roughness length for intermediate output only (m) real (kind=kind_phys) :: z0h !roughness length, sensible heat (m) real (kind=kind_phys) :: z0hg !roughness length, sensible heat (m) + real (kind=kind_phys) :: rb !bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: ramc !aerodynamic resistance for momentum (s/m) real (kind=kind_phys) :: rahc !aerodynamic resistance for sensible heat (s/m) real (kind=kind_phys) :: rawc !aerodynamic resistance for water vapor (s/m) @@ -4137,7 +4133,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd ,snowh ,shdfac ,garea1 , & !in - varf,psl_gwd_z0m_factor, & !in + varf,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -4447,7 +4443,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & emg ,stc ,df ,rsurf ,lathea , & !in gamma ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in - varf,psl_gwd_z0m_factor, & !in + varf,gwd_z0m_factor, & !in pblhx , iz0tlnd , itime ,psi_opt,ep_1,ep_2,epsm1,cp ,& #ifdef CCPP tgb ,cm ,ch,ustarx,errmsg ,errflg , & !inout @@ -4502,8 +4498,8 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(in) :: rhsur !< raltive humidity in surface soil/snow air space (-) real (kind=kind_phys), intent(in) :: fsno !< snow fraction - real (kind=kind_phys), intent(in) :: varf ! - real (kind=kind_phys), intent(in) :: psl_gwd_z0m_factor ! + real (kind=kind_phys), intent(in) :: varf !< sub-grid orography standard deviation [m] + real (kind=kind_phys), intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys), intent(in) :: pblhx !< pbl height (m) real (kind=kind_phys), intent(in) :: ep_1 !< real (kind=kind_phys), intent(in) :: ep_2 !< @@ -4747,7 +4743,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd,snowh ,shdfac ,garea1 , & !in - varf,psl_gwd_z0m_factor, & !in + varf,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5454,7 +5450,7 @@ end subroutine sfcdif2 subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h,zpd ,snowh ,fveg ,garea1 , & !in - varf,psl_gwd_z0m_factor, & !in + varf,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5484,8 +5480,8 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: snowh !< snow depth [m] real (kind=kind_phys), intent(in ) :: fveg !< fractional vegetation cover real (kind=kind_phys), intent(in ) :: garea1 !< grid area [km2] - real (kind=kind_phys), intent(in ) :: varf ! standard deviation org [m] - real (kind=kind_phys), intent(in ) :: psl_gwd_z0m_factor ! z0m factor in psl tofd + real (kind=kind_phys), intent(in ) :: varf !< sub-grid orography standard deviation [m] + real (kind=kind_phys), intent(in ) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys), intent(inout) :: ustarx !< friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm !< momentum stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fh !< sen heat stability correction, weighted by prior iters @@ -5540,7 +5536,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur endif call gfs_stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & - varf,psl_gwd_z0m_factor,rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) + varf,gwd_z0m_factor,rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) end subroutine sfcdif3 @@ -5550,7 +5546,7 @@ subroutine gfs_stability & ! --- inputs: ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & thsfc_loc, & - varf,psl_gwd_z0m_factor, & + varf,gwd_z0m_factor, & ! --- outputs: rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) @@ -5570,8 +5566,8 @@ subroutine gfs_stability & real(kind=kind_phys), intent(in) :: ztmax ! thermal roughness length real(kind=kind_phys), intent(in) :: tvs ! surface virtual temperature real(kind=kind_phys), intent(in) :: grav ! local gravity -real(kind=kind_phys), intent(in) :: varf ! turbulent scale oro std -real(kind=kind_phys), intent(in) :: psl_gwd_z0m_factor ! tofd factor +real(kind=kind_phys), intent(in) :: varf ! turbulent orographic standard deviation [m] +real(kind=kind_phys), intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor logical, intent(in) :: thsfc_loc ! use local theta reference flag real(kind=kind_phys), intent(out) :: rb ! bulk richardson number [-] @@ -5630,11 +5626,10 @@ subroutine gfs_stability & real(kind=kind_phys) :: zolmax real(kind=kind_phys) xkzo -! tofd - real(kind=kind_phys) :: cf, zf, fri, fphim - real(kind=kind_phys), parameter :: varf_min = 50.0 - real(kind=kind_phys), parameter :: varf_max = 500.0 -!!! real(kind=kind_phys), parameter :: psl_gwd_z0m_factor = 0.003 +! +real(kind=kind_phys) :: cf, zf, fri, fphim +real(kind=kind_phys), parameter :: varf_min = 50.0 +real(kind=kind_phys), parameter :: varf_max = 500.0 z1i = one / z1 ! inverse of model height @@ -5755,13 +5750,15 @@ subroutine gfs_stability & endif ! end of if (dtv >= 0 ) then loop ! -! form drag coefficient +! compute the effect of orographic form drag : koo et al. (2018,jgr, eqs. 1 and 2) ! - if ( varf.gt.varf_min .and. psl_gwd_z0m_factor.gt.0.0 ) then - zf = min( min(varf,varf_max)*psl_gwd_z0m_factor,z1 ) - fri = min( max( 1.-rb,0. ), 1.) - fphim = log( ( z1 + zf) / zf ) - cf = ca*ca / (fphim*fphim) * fri + if ( varf.gt.varf_min .and. gwd_z0m_factor.gt.0.0_kind_phys ) then + tem1 = min(varf, varf_max) + zf = tem1 * gwd_z0m_factor ! effective roughness length + zf = min(zf, z1) + fri = min( max( 1.0_kind_phys-rb,0.0_kind_phys ), 1.0_kind_phys) ! a function of ther bulk Richardson number + fphim = log( ( z1 + zf) / zf ) ! integrated profile function for momentum + cf = ca*ca / (fphim*fphim) * fri ! exchange coefficient for momentum else cf = 0. endif @@ -5777,7 +5774,7 @@ subroutine gfs_stability & tem1 = 0.00001_kind_phys/z1 ! minimum exhange coef (?) cm = max(cm, tem1) ch = max(ch, tem1) - stress = (cm + cf)* wind * wind ! cf for tofd 10.15) or (10.19) in Arya + stress = (cm + cf)* wind * wind ustar = sqrt(stress) ! friction velocity return diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index 90cc14a8d..93f8719b9 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -137,7 +137,7 @@ subroutine noahmpdrv_run & rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, rhonewsn1,& con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1, & con_fvirt, con_rd, con_hfus, thsfc_loc, cpllnd, cpllnd2atm,& - nmtvr,mntvar,psl_gwd_z0m_factor, & + mntvar,gwd_z0m_factor, & ! --- in/outs: weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & @@ -158,7 +158,7 @@ subroutine noahmpdrv_run & sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, & cmm, chh, evbs, evcw, sbsno, pah, ecan, etran, edir, snowc,& stm, snohf,smcwlt2, smcref2, wet1, t2mmp, q2mp,zvfun, & - ztmax, rca, errmsg, errflg, & + ztmax, errmsg, errflg, & canopy_heat_storage_ccpp, & rainfall_ccpp, & sw_absorbed_total_ccpp, & @@ -288,7 +288,6 @@ subroutine noahmpdrv_run & integer , intent(in) :: iopt_stc ! option for snow/soil temperature time scheme (only layer 1) integer , intent(in) :: iopt_trs ! option for thermal roughness scheme integer , intent(in) :: iopt_diag ! option for surface diagnose approach - integer , intent(in) :: nmtvr ! number of subgrid scale oroghic data real(kind=kind_phys), dimension(:) , intent(in) :: xlatin ! latitude real(kind=kind_phys), dimension(:) , intent(in) :: xcoszin ! cosine of zenith angle integer , intent(in) :: iyrlen ! year length [days] @@ -300,7 +299,7 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(in) :: graupel_mp ! microphysics graupel [mm] real(kind=kind_phys), dimension(:) , intent(in) :: ice_mp ! microphysics ice/hail [mm] real(kind=kind_phys), dimension(:,:) , intent(in) :: mntvar ! subgrid orographic statistics data - real(kind=kind_phys) , intent(in) :: psl_gwd_z0m_factor ! psl tofd zom factor + real(kind=kind_phys) , intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor real(kind=kind_phys), dimension(:) , intent(in) :: rhonewsn1 ! precipitation ice density (kg/m^3) real(kind=kind_phys) , intent(in) :: con_hvap ! latent heat condensation [J/kg] real(kind=kind_phys) , intent(in) :: con_cp ! specific heat air [J/kg/K] @@ -404,8 +403,6 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(out) :: q2mp ! combined q2m from tiles real(kind=kind_phys), dimension(:) , intent(out) :: zvfun ! real(kind=kind_phys), dimension(:) , intent(out) :: ztmax ! thermal roughness length - real(kind=kind_phys), dimension(:) , intent(out) :: rca ! total canopy/stomatal resistance (s/m) - character(len=*) , intent(out) :: errmsg integer , intent(out) :: errflg @@ -639,7 +636,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: cq2 real (kind=kind_phys) :: qfx real (kind=kind_phys) :: wspd1 ! wind speed with all components - real (kind=kind_phys) :: varf_grid ! std deviation orography + real (kind=kind_phys) :: varf_grid ! sub-grid orography standard deviation real (kind=kind_phys) :: pblhx ! height of pbl integer :: mnice @@ -956,10 +953,6 @@ subroutine noahmpdrv_run & ch_vegetated = 0.0 ch_bare_ground = ch_noahmp canopy_heat_storage = 0.0 - lai_sunlit = 0.0 - lai_shaded = 0.0 - rs_sunlit = 0.0 - rs_shaded = 0.0 else ! not glacier @@ -977,7 +970,7 @@ subroutine noahmpdrv_run & spec_humidity_forcing ,area_grid ,cloud_water_forcing , & sw_radiation_forcing ,radiation_lw_forcing ,thsfc_loc , & prslkix ,prsik1x ,prslk1x , & - varf_grid,psl_gwd_z0m_factor , & + varf_grid,gwd_z0m_factor , & pblhx ,iz0tlnd ,itime , & psi_opt , & precip_convective , & @@ -1070,16 +1063,6 @@ subroutine noahmpdrv_run & zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total - !LAI-scale canopy resistance based on weighted sunlit shaded fraction - if(rs_sunlit .le. 0.0 .or. rs_shaded .le. 0.0 .or. & - lai_sunlit .eq. 0.0 .or. lai_shaded .eq. 0.0) then - rca(i) = parameters%rsmax - else !calculate LAI-scale canopy conductance (1/Rs) - rca(i) = ((1.0/(rs_sunlit+leaf_air_resistance)*lai_sunlit) + & - ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) - rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance - end if - smc (i,:) = soil_moisture_vol slc (i,:) = soil_liquid_vol snowxy (i) = float(snow_levels) @@ -1222,7 +1205,7 @@ subroutine noahmpdrv_run & call gfs_stability & (zf(i), zvfun(i), gdx, virtual_temperature, vptemp,wind(i), z0_total, z0h_total, & tvs1, con_g, thsfc_loc, & - varf_grid, psl_gwd_z0m_factor, & + varf_grid, gwd_z0m_factor, & rb1(i), fm1(i), fh1(i), fm101(i), fh21(i), cm(i), ch(i), stress1(i), ustar1(i)) rmol1(i) = undefined !not used in GFS sfcdif -> to satsify output diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta index 3e24624fa..f6b57f406 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta @@ -649,13 +649,6 @@ dimensions = () type = logical intent = in -[nmtvr] - standard_name = number_of_statistical_measures_of_subgrid_orography - long_name = number of statistical measures of subgrid height_above_mean_sea_level - units = count - dimensions = () - type = integer - intent = in [mntvar] standard_name = statistical_measures_of_subgrid_orography_collection_array long_name = array of statistical measures of subgrid height_above_mean_sea_level @@ -664,9 +657,9 @@ type = real kind = kind_phys intent = in -[psl_gwd_z0m_factor] +[gwd_z0m_factor] standard_name = momentum_roughness_factor_in_turbulent_form_drag - long_name = multiplication of Z0m in PSL turbulent form drag + long_name = multiplication of orograhic standard deviation units = count dimensions = () type = real @@ -1383,14 +1376,6 @@ type = real kind = kind_phys intent = out -[rca] - standard_name = aerodynamic_resistance_in_canopy - long_name = canopy resistance - units = s m-1 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 1eb7f7c5dcb9ef0c7a8c9f68fb7be5bfb0e61ff5 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Thu, 23 May 2024 08:58:55 -0400 Subject: [PATCH 04/12] update files module_sf_noahmplsm.F90 noahmpdrv.F90 noahmpdrv.meta module_sf_noahmp_glacier.F90 --- .../Land/Noahmp/module_sf_noahmplsm.F90 | 142 +++++++++--------- 1 file changed, 71 insertions(+), 71 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 9e9dfed91..8ac3b21d7 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -90,7 +90,7 @@ module module_sf_noahmplsm integer :: opt_crs !< options for canopy stomatal resistance ! **1 -> ball-berry - ! 2 -> jarvis + ! 2 -> jarvis integer :: opt_btr !< options for soil moisture factor for stomatal resistance ! **1 -> noah (soil moisture) @@ -103,16 +103,16 @@ module module_sf_noahmplsm ! 3 -> original surface and subsurface runoff (free drainage) ! 4 -> bats surface and subsurface runoff (free drainage) ! 5 -> miguez-macho&fan groundwater scheme (miguez-macho et al. 2007 jgr; fan et al. 2007 jgr) - ! (needs further testing for public use) + ! (needs further testing for public use) integer :: opt_sfc !< options for surface layer drag coeff (ch & cm) ! **1 -> m-o - ! **2 -> original noah (chen97) - ! **3 -> myj consistent; 4->ysu consistent. mb: removed in v3.7 for further testing + ! **2 -> original noah (chen97) + ! **3 -> myj consistent; 4->ysu consistent. mb: removed in v3.7 for further testing integer :: opt_frz !< options for supercooled liquid water (or ice fraction) ! **1 -> no iteration (niu and yang, 2006 jhm) - ! 2 -> koren's iteration + ! 2 -> koren's iteration integer :: opt_inf !< options for frozen soil permeability ! **1 -> linear effects, more permeable (niu and yang, 2006, jhm) @@ -125,13 +125,13 @@ module module_sf_noahmplsm integer :: opt_alb !< options for ground snow surface albedo ! 1 -> bats - ! **2 -> class + ! **2 -> class integer :: opt_snf !< options for partitioning precipitation into rainfall & snowfall ! **1 -> jordan (1991) - ! 2 -> bats: when sfctmp sfctmp < tfrz - ! 4 -> use wrf microphysics output + ! 2 -> bats: when sfctmp sfctmp < tfrz + ! 4 -> use wrf microphysics output integer :: opt_tbot !< options for lower boundary condition of soil temperature ! 1 -> zero heat flux from bottom (zbot and tbot not used) @@ -139,20 +139,20 @@ module module_sf_noahmplsm integer :: opt_stc !< options for snow/soil temperature time scheme (only layer 1) ! **1 -> semi-implicit; flux top boundary condition - ! 2 -> full implicit (original noah); temperature top boundary condition + ! 2 -> full implicit (original noah); temperature top boundary condition ! 3 -> same as 1, but fsno for ts calculation (generally improves snow; v3.7) integer :: opt_rsf !< options for surface resistent to evaporation/sublimation ! **1 -> sakaguchi and zeng, 2009 - ! 2 -> sellers (1992) + ! 2 -> sellers (1992) ! 3 -> adjusted sellers to decrease rsurf for wet soil - ! 4 -> option 1 for non-snow; rsurf = rsurf_snow for snow (set in mptable); ad v3.8 + ! 4 -> option 1 for non-snow; rsurf = rsurf_snow for snow (set in mptable); ad v3.8 integer :: opt_soil !< options for defining soil properties ! **1 -> use input dominant soil texture - ! 2 -> use input soil texture that varies with depth + ! 2 -> use input soil texture that varies with depth ! 3 -> use soil composition (sand, clay, orgm) and pedotransfer functions (opt_pedo) - ! 4 -> use input soil properties (bexp_3d, smcmax_3d, etc.) + ! 4 -> use input soil properties (bexp_3d, smcmax_3d, etc.) integer :: opt_pedo !< options for pedotransfer functions (used when opt_soil = 3) ! **1 -> saxton and rawls (2006) @@ -425,7 +425,7 @@ subroutine noahmp_sflx (parameters, & qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing varf , gwd_z0m_factor, & ! in : turbulent orographic form drag pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -437,7 +437,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -446,9 +446,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP @@ -822,7 +822,7 @@ subroutine noahmp_sflx (parameters, & canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + fwet ,cmc ) !out ! compute energy budget (momentum & energy fluxes and phase changes) @@ -837,7 +837,7 @@ subroutine noahmp_sflx (parameters, & thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in varf ,gwd_z0m_factor , & !in pblhx ,iz0tlnd, itime ,psi_opt, ep_1, ep_2, epsm1,cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -858,7 +858,7 @@ subroutine noahmp_sflx (parameters, & fsrg ,rssun ,rssha ,albd ,albi ,albsnd,albsni, bgap ,wgap, tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb , & !out emissi ,pah ,canhs, & - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out qsfcveg = eah*ep_2/(sfcprs + epsm1*eah) qsfcbare = qsfc @@ -881,7 +881,7 @@ subroutine noahmp_sflx (parameters, & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc,jloc , smceq , & !in bdfall ,fp ,rain ,snow , & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -915,9 +915,9 @@ subroutine noahmp_sflx (parameters, & if (opt_crop == 1 .and. crop_active) then call carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julian , & !in dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in - soldn ,t2m , & !in + soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - lai ,sai ,gdd , & !inout + lai ,sai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out end if @@ -968,7 +968,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & prcpconv,prcpnonc ,prcpshcv,prcpsnow,prcpgrpl,prcphail , & soldn ,cosz ,thair ,qair , & eair ,rhoair ,qprecc ,qprecl ,solad , solai , & - swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) + swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) ! -------------------------------------------------------------------------------------------------- ! re-process atmospheric forcing ! ---------------------------------------------------------------------- @@ -1041,7 +1041,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4) then qprecc = prcpconv + prcpshcv - qprecl = prcpnonc + qprecl = prcpnonc else qprecc = 0.10 * prcp ! should be from the atmospheric model qprecl = 0.90 * prcp ! should be from the atmospheric model @@ -1094,15 +1094,15 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4 .or. opt_snf == 5) then prcp_frozen = prcpsnow + prcpgrpl + prcphail if(prcpnonc > 0. .and. prcp_frozen > 0.) then - fpice = min(1.0,prcp_frozen/prcpnonc) - fpice = max(0.0,fpice) + fpice = min(1.0,prcp_frozen/prcpnonc) + fpice = max(0.0,fpice) if(opt_snf==4) bdfall = bdfall*(prcpsnow/prcp_frozen) + rho_grpl*(prcpgrpl/prcp_frozen) + & rho_hail*(prcphail/prcp_frozen) if(opt_snf==5) bdfall = parameters%prcpiceden - else - fpice = 0.0 + else + fpice = 0.0 endif - + endif rain = prcp * (1.-fpice) @@ -1237,8 +1237,8 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv bdfall ,rain ,snow ,fp , & !in canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out - pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out + fwet ,cmc ) !out ! ------------------------ code history ------------------------------ ! michael barlage: oct 2013 - split canwater to calculate precip movement for @@ -1340,10 +1340,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qintr = 0. qdripr = 0. qthror = rain - if(canliq > 0.) then ! for case of canopy getting buried - qdripr = qdripr + canliq/dt - canliq = 0.0 - end if + if(canliq > 0.) then ! for case of canopy getting buried + qdripr = qdripr + canliq/dt + canliq = 0.0 + end if end if ! heat transported by liquid water @@ -1366,8 +1366,8 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qints = max(qints, 0.) ft = max(0.0,(tv - 270.15) / 1.87e5) fv = sqrt(uu*uu + vv*vv) / 1.56e5 - ! mb: changed below to reflect the rain assumption that all precip gets intercepted - icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt +! mb: changed below to reflect the rain assumption that all precip gets intercepted + icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt qdrips = (fveg * snow - qints) + icedrip qthros = (1.0-fveg) * snow canice= max(0.,canice + (qints - icedrip)*dt) @@ -1375,10 +1375,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qints = 0. qdrips = 0. qthros = snow - if(canice > 0.) then ! for case of canopy getting buried - qdrips = qdrips + canice/dt - canice = 0.0 - end if + if(canice > 0.) then ! for case of canopy getting buried + qdrips = qdrips + canice/dt + canice = 0.0 + end if endif ! print*, "precip_heat canopy through:",3600.0*(fveg * snow - qints) ! print*, "precip_heat canopy drip:",3600.0*max(0.,canice) * (fv+ft) @@ -1408,13 +1408,13 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv if (fveg > 0.0 .and. fveg < 1.0) then pahg = pahg / fveg ! these will be multiplied by fraction later - pahb = pahb / (1.0-fveg) + pahb = pahb / (1.0-fveg) elseif (fveg <= 0.0) then pahb = pahg + pahb ! for case of canopy getting buried pahg = 0.0 - pahv = 0.0 + pahv = 0.0 elseif (fveg >= 1.0) then - pahb = 0.0 + pahb = 0.0 end if pahv = max(pahv,-20.0) ! put some artificial limits here for stability @@ -1680,9 +1680,9 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in - varf,gwd_z0m_factor , & !in + varf ,gwd_z0m_factor , & !in pblhx , iz0tlnd, itime,psi_opt,ep_1, ep_2, epsm1, cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -1702,7 +1702,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2mv ,t2mb ,fsrv , & fsrg ,rssun ,rssha ,albd ,albi,albsnd ,albsni,bgap ,wgap,tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb, emissi,pah,canhs,& - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out !jref:end ! -------------------------------------------------------------------------------------------------- @@ -2216,19 +2216,19 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in if (tv .gt. tfrz) then ! barlage: add distinction between ground and latheav = hvap ! vegetation in v3.6 - frozen_canopy = .false. + frozen_canopy = .false. else latheav = hsub - frozen_canopy = .true. + frozen_canopy = .true. end if gammav = cpair*sfcprs/(ep_2*latheav) if (tg .gt. tfrz) then latheag = hvap - frozen_ground = .false. + frozen_ground = .false. else latheag = hsub - frozen_ground = .true. + frozen_ground = .true. end if gammag = cpair*sfcprs/(ep_2*latheag) @@ -2340,7 +2340,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in ssoil = fveg * ghv + (1.0 - fveg) * ghb fcev = evc fctr = tr - pah = fveg * pahg + (1.0 - fveg) * pahb + pahv + pah = fveg * pahg + (1.0 - fveg) * pahb + pahv tg = fveg * tgv + (1.0 - fveg) * tgb t2m = fveg * t2mv + (1.0 - fveg) * t2mb ts = fveg * tah + (1.0 - fveg) * tgb @@ -2370,7 +2370,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2m = t2mb fcev = 0. fctr = 0. - pah = pahb + pah = pahb ts = tg cm = cmb ch = chb @@ -3541,7 +3541,7 @@ subroutine twostream (parameters,ib ,ic ,vegtyp ,cosz ,vai , & ! kopen = 1.0 else if(opt_rad == 1) then - denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) + denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) hd = parameters%hvt - parameters%hvb bb = 0.5 * hd thetap = atan(bb/parameters%rc * tan(acos(max(0.01,cosz))) ) @@ -4256,11 +4256,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & shc = fveg*rhoair*cpair*cvh * ( tv-tah) evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav - if (tv > tfrz) then + if (tv > tfrz) then evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 - else + else evc = min(canice*latheav/dt,evc) - end if + end if ! canopy heat capacity hcv = fveg*(parameters%cbiom*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice) !j/m2/k @@ -5450,7 +5450,7 @@ end subroutine sfcdif2 subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h,zpd ,snowh ,fveg ,garea1 , & !in - varf,gwd_z0m_factor, & !in + varf ,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -7056,7 +7056,7 @@ subroutine water (parameters,vegtyp ,nsnow ,nsoil ,imelt ,dt ,uu , & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc ,jloc ,smceq , & !in bdfall ,fp ,rain ,snow, & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -7685,19 +7685,19 @@ subroutine combine (parameters,nsnow ,nsoil ,iloc ,jloc , & !in snice(j-1) = snice(j-1) + snice(j) dzsnso(j-1) = dzsnso(j-1) + dzsnso(j) else - if(snice(j) >= 0.) then + if(snice(j) >= 0.) then ponding1 = snliq(j) ! isnow will get set to zero below; ponding1 will get sneqv = snice(j) ! added to ponding from phasechange ponding should be snowh = dzsnso(j) ! zero here because it was calculated for thin snow - else ! snice over-sublimated earlier - ponding1 = snliq(j) + snice(j) - if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil - sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) + else ! snice over-sublimated earlier + ponding1 = snliq(j) + snice(j) + if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil + sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) ponding1 = 0.0 - end if + end if sneqv = 0.0 snowh = 0.0 - end if + end if snliq(j) = 0.0 snice(j) = 0.0 dzsnso(j) = 0.0 @@ -9786,7 +9786,7 @@ subroutine carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julia dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - xlai ,xsai ,gdd , & !inout + xlai ,xsai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out ! ------------------------------------------------------------------------------------------ ! initial crop version created by xing liu @@ -10465,7 +10465,7 @@ end subroutine psn_crop !! subroutine noahmp_options(idveg , iopt_crs , iopt_btr , iopt_run , iopt_sfc , iopt_frz , & iopt_inf, iopt_rad , iopt_alb , iopt_snf , iopt_tbot, iopt_stc , & - iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & + iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & iopt_z0m ) implicit none From 2259ee9b7244ba740ba35ac9c52b28238ba89fc7 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Sun, 26 May 2024 11:45:49 -0400 Subject: [PATCH 05/12] fix the compile problem for HR4-GWD-update --- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index e4623b83a..ca2a0db42 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -158,7 +158,7 @@ subroutine noahmpdrv_run & sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, & cmm, chh, evbs, evcw, sbsno, pah, ecan, etran, edir, snowc,& stm, snohf,smcwlt2, smcref2, wet1, t2mmp, q2mp,zvfun, & - ztmax, errmsg, errflg, & + ztmax, rca, errmsg, errflg, & canopy_heat_storage_ccpp, & rainfall_ccpp, & sw_absorbed_total_ccpp, & @@ -1064,6 +1064,16 @@ subroutine noahmpdrv_run & zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total + !LAI-scale canopy resistance based on weighted sunlit shaded fraction + if(rs_sunlit .le. 0.0 .or. rs_shaded .le. 0.0 .or. & + lai_sunlit .eq. 0.0 .or. lai_shaded .eq. 0.0) then + rca(i) = parameters%rsmax + else !calculate LAI-scale canopy conductance (1/Rs) + rca(i) = ((1.0/(rs_sunlit+leaf_air_resistance)*lai_sunlit) + & + ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) + rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance + end if + smc (i,:) = soil_moisture_vol slc (i,:) = soil_liquid_vol snowxy (i) = float(snow_levels) From 622aa411642175b13c330599e62580e7ac5307b1 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Fri, 31 May 2024 13:45:31 -0400 Subject: [PATCH 06/12] update Noahmp code for HR4 --- physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 | 8 +++++--- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 4 ++++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 8ac3b21d7..a3b466f46 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -2271,7 +2271,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in csigmaf1, & !out !jref:start qc ,qsfc ,psfc , & !in - q2v ,chv2, chleaf, chuc) !inout + q2v ,chv2 ,chleaf ,chuc , & + rb) !out ! new coupling code @@ -3722,7 +3723,8 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & t2mv ,psnsun ,psnsha ,canhs , & !out csigmaf1, & !out qc ,qsfc ,psfc , & !in - q2v ,cah2 ,chleaf ,chuc ) !inout + q2v ,cah2 ,chleaf ,chuc , & !inout + rb) !out ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3848,6 +3850,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(out) :: chuc !< under canopy exchange coefficient real (kind=kind_phys), intent(out) :: canhs !< canopy heat storage change (w/m2) real (kind=kind_phys), intent(out) :: q2v !< + real (kind=kind_phys), intent(out) :: rb !< bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: cah !< sensible heat conductance, canopy air to zlvl air (m/s) real (kind=kind_phys) :: u10v !< 10 m wind speed in eastward dir (m/s) real (kind=kind_phys) :: v10v !< 10 m wind speed in eastward dir (m/s) @@ -3864,7 +3867,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys) :: z0mo !roughness length for intermediate output only (m) real (kind=kind_phys) :: z0h !roughness length, sensible heat (m) real (kind=kind_phys) :: z0hg !roughness length, sensible heat (m) - real (kind=kind_phys) :: rb !bulk leaf boundary layer resistance (s/m) real (kind=kind_phys) :: ramc !aerodynamic resistance for momentum (s/m) real (kind=kind_phys) :: rahc !aerodynamic resistance for sensible heat (s/m) real (kind=kind_phys) :: rawc !aerodynamic resistance for water vapor (s/m) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index ca2a0db42..4af3b93b6 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -954,6 +954,10 @@ subroutine noahmpdrv_run & ch_vegetated = 0.0 ch_bare_ground = ch_noahmp canopy_heat_storage = 0.0 + lai_sunlit = 0.0 + lai_shaded = 0.0 + rs_sunlit = 0.0 + rs_shaded = 0.0 else ! not glacier From 6b783323e2dd8553cc46ef0ffff1adbae3e095ff Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Fri, 31 May 2024 13:46:14 -0400 Subject: [PATCH 07/12] update GWD code for HR4 --- physics/GWD/drag_suite.F90 | 61 +++++++++++++++++++-------------- physics/GWD/drag_suite.meta | 10 +++++- physics/GWD/ugwpv1_gsldrag.F90 | 13 +++---- physics/GWD/ugwpv1_gsldrag.meta | 10 +++++- physics/GWD/unified_ugwp.F90 | 16 +++++---- physics/GWD/unified_ugwp.meta | 10 +++++- 6 files changed, 78 insertions(+), 42 deletions(-) diff --git a/physics/GWD/drag_suite.F90 b/physics/GWD/drag_suite.F90 index 0445cbeb5..91edfbec8 100644 --- a/physics/GWD/drag_suite.F90 +++ b/physics/GWD/drag_suite.F90 @@ -219,8 +219,8 @@ subroutine drag_suite_run( & & dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl, & & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & & slmsk,br1,hpbl, & - & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, & - & lprnt, ipr, rdxzb, dx, gwd_opt, & + & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, & + & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, & & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & & dtend, dtidx, index_of_process_orographic_gwd, & & index_of_temperature, index_of_x_wind, & @@ -327,18 +327,18 @@ subroutine drag_suite_run( & integer, intent(in) :: gwd_opt logical, intent(in) :: lprnt integer, intent(in) :: KPBL(:) - real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:) + real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, & + & cdmbgwd(:), alpha_fd real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:) logical, intent(in) :: ldiag3d - integer, intent(in) :: dtidx(:,:) - integer, intent(in) :: index_of_temperature, & + integer, intent(in) :: dtidx(:,:), index_of_temperature, & & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind integer :: kpblmax integer, parameter :: ims=1, kms=1, its=1, kts=1 real(kind=kind_phys), intent(in) :: fv, pi real(kind=kind_phys) :: rcl, cdmb - real(kind=kind_phys) :: g_inv + real(kind=kind_phys) :: g_inv, rd_inv real(kind=kind_phys), intent(inout) :: & & dudt(:,:),dvdt(:,:), & @@ -444,6 +444,7 @@ subroutine drag_suite_run( & real(kind=kind_phys), dimension(im,km) :: utendform,vtendform real(kind=kind_phys) :: a1,a2,wsp real(kind=kind_phys) :: H_efold + real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3 ! critical richardson number for wave breaking : ! larger drag with larger value real(kind=kind_phys), parameter :: ric = 0.25 @@ -708,11 +709,12 @@ subroutine drag_suite_run( & taufb(1:im,1:km+1) = 0.0 komax(1:im) = 0 ! + rd_inv = 1./rd do k = kts,km do i = its,im vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k)) vtk(i,k) = vtj(i,k) / prslk(i,k) - ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3 + ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3 enddo enddo ! @@ -1363,8 +1365,10 @@ subroutine drag_suite_run( & H_efold = 1500. DO k=kts,km wsp=SQRT(uwnd1(i,k)**2 + vwnd1(i,k)**2) - ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 - var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & + ! Note: In Beljaars et al. (2004): + ! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 + ! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys + var_temp = alpha_fd*coeff_fd*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero ! Note: This is a semi-implicit treatment of the time differencing ! per Beljaars et al. (2004, QJRMS) @@ -1427,8 +1431,8 @@ subroutine drag_suite_psl( & & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, & & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & & slmsk,br1,hpbl,vtype, & - & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, & - & lprnt, ipr, rdxzb, dx, gwd_opt, & + & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, & + & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, & & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, & & psl_gwd_dx_factor, & & dtend, dtidx, index_of_process_orographic_gwd, & @@ -1537,7 +1541,8 @@ subroutine drag_suite_psl( & integer, intent(in) :: gwd_opt logical, intent(in) :: lprnt integer, intent(in) :: KPBL(:) - real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(:) + real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, & + & cdmbgwd(:), alpha_fd real(kind=kind_phys), intent(inout) :: dtend(:,:,:) logical, intent(in) :: ldiag3d integer, intent(in) :: dtidx(:,:), index_of_temperature, & @@ -1547,7 +1552,7 @@ subroutine drag_suite_psl( & integer, parameter :: ims=1, kms=1, its=1, kts=1 real(kind=kind_phys), intent(in) :: fv, pi real(kind=kind_phys) :: rcl, cdmb - real(kind=kind_phys) :: g_inv + real(kind=kind_phys) :: g_inv, g_cp real(kind=kind_phys), intent(inout) :: & & dudt(:,:),dvdt(:,:), & @@ -1649,6 +1654,7 @@ subroutine drag_suite_psl( & real(kind=kind_phys), dimension(im,km) :: utendform,vtendform real(kind=kind_phys) :: a1,a2,wsp real(kind=kind_phys) :: H_efold + real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3 ! multification factor of standard deviation : ! larger drag with larger value !!! real(kind=kind_phys), parameter :: psl_gwd_dx_factor = 6.0 @@ -1789,18 +1795,18 @@ subroutine drag_suite_psl( & enddo !--- calculate scale-aware tapering factors -do i=1,im - if ( dx(i) .ge. dxmax_ls ) then - ls_taper(i) = 1. - else - if ( dx(i) .le. dxmin_ls) then - ls_taper(i) = 0. + do i=1,im + if ( dx(i) .ge. dxmax_ls ) then + ls_taper(i) = 1. else - ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & + if ( dx(i) .le. dxmin_ls) then + ls_taper(i) = 0. + else + ls_taper(i) = 0.5 * ( SIN(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ & (dxmax_ls-dxmin_ls)) + 1. ) + endif endif - endif -enddo + enddo ! Remove ss_tapering ss_taper(:) = 1. @@ -2072,6 +2078,7 @@ subroutine drag_suite_psl( & IF ( (do_gsl_drag_ls_bl).and. & ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then + g_cp = g/cp do i=its,im if ( ls_taper(i).GT.1.E-02 ) then @@ -2086,7 +2093,7 @@ subroutine drag_suite_psl( & tem2 = v1(i,k) - v1(i,k+1) dw2 = rcl*(tem1*tem1 + tem2*tem2) shr2 = max(dw2,dw2min) * rdz * rdz - bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti + bvf2 = g*(g_cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti usqj(i,k) = max(bvf2/shr2,rimin) bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k)) bnv2(i,k) = max( bnv2(i,k), bnv2min ) @@ -2341,8 +2348,10 @@ subroutine drag_suite_psl( & 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 - var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & + ! Note: In Beljaars et al. (2004): + ! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759 + ! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys + var_temp = alpha_fd*coeff_fd*EXP(-(zl(i,k)/H_efold)**1.5)*a2* & zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero ! Note: This is a semi-implicit treatment of the time differencing ! per Beljaars et al. (2004, QJRMS) @@ -2465,7 +2474,7 @@ subroutine drag_suite_psl( & do k = km, kpblmin, -1 if(flag(i).and. k.le.kbmax(i)) then pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* & - del(i,k)/g/ro(i,k) + del(i,k)*g_inv/ro(i,k) ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.) ! !---------- apply flow-blocking drag when pe >= ke diff --git a/physics/GWD/drag_suite.meta b/physics/GWD/drag_suite.meta index 323d89d36..088387cfd 100644 --- a/physics/GWD/drag_suite.meta +++ b/physics/GWD/drag_suite.meta @@ -528,6 +528,14 @@ type = real kind = kind_phys intent = in +[alpha_fd] + standard_name = alpha_coefficient_for_turbulent_orographic_form_drag + long_name = alpha coefficient for Beljaars et al turbulent orographic form drag + units = none + dimensions = () + type = real + kind = kind_phys + intent = in [me] standard_name = mpi_rank long_name = rank of the current MPI task @@ -603,7 +611,7 @@ [psl_gwd_dx_factor] standard_name = effective_grid_spacing_of_psl_gwd_suite long_name = multiplication of grid spacing - units = count + units = 1 dimensions = () type = real kind = kind_phys diff --git a/physics/GWD/ugwpv1_gsldrag.F90 b/physics/GWD/ugwpv1_gsldrag.F90 index d0eba4ae1..fc90955bd 100644 --- a/physics/GWD/ugwpv1_gsldrag.F90 +++ b/physics/GWD/ugwpv1_gsldrag.F90 @@ -309,7 +309,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, do_gwd_opt_psl, psl_gwd_dx_factor, & do_ugwp_v1, do_ugwp_v1_orog_only, & do_ugwp_v1_w_gsldrag, gwd_opt, do_tofd, ldiag_ugwp, ugwp_seq_update, & - cdmbgwd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, & + cdmbgwd, alpha_fd, jdat, nmtvr, hprime, oc, theta, sigma, gamma, & elvmax, clx, oa4, varss,oc1ss,oa4ss,ol4ss, dx, xlat, xlat_d, sinlat, coslat, & area, rain, br1, hpbl,vtype, kpbl, slmsk, & ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, del, tau_amf, & @@ -375,7 +375,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, ! SSO parameters and variables integer, intent(in) :: gwd_opt !gwd_opt and nmtvr are "redundant" controls integer, intent(in) :: nmtvr - real(kind=kind_phys), intent(in) :: cdmbgwd(:) ! for gsl_drag + real(kind=kind_phys), intent(in) :: cdmbgwd(:), alpha_fd ! for gsl_drag real(kind=kind_phys), intent(in), dimension(:) :: hprime, oc, theta, sigma, gamma @@ -404,8 +404,7 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, integer, intent(in), dimension(:) :: vtype real(kind=kind_phys), intent(in), dimension(:) :: rain - real(kind=kind_phys), intent(in), dimension(:) :: br1, slmsk - real(kind=kind_phys), intent(in), dimension(:) :: hpbl + real(kind=kind_phys), intent(in), dimension(:) :: br1, hpbl, slmsk ! ! moved to GFS_phys_time_vary ! real(kind=kind_phys), intent(in), dimension(:) :: ddy_j1tau, ddy_j2tau @@ -563,7 +562,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, & slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, & con_fv, con_pi, lonr, & - cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + cdmbgwd(1:2),alpha_fd,me,master, & + lprnt,ipr,rdxzb,dx,gwd_opt, & do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & psl_gwd_dx_factor, & dtend, dtidx, index_of_process_orographic_gwd, & @@ -583,7 +583,8 @@ subroutine ugwpv1_gsldrag_run(me, master, im, levs, ak, bk, ntrac, lonr, dtp, du_osscol, dv_osscol, du_ofdcol, dv_ofdcol, & slmsk,br1,hpbl,con_g,con_cp,con_rd,con_rv, & con_fv, con_pi, lonr, & - cdmbgwd(1:2),me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + cdmbgwd(1:2),alpha_fd,me,master, & + lprnt,ipr,rdxzb,dx,gwd_opt, & do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & dtend, dtidx, index_of_process_orographic_gwd, & index_of_temperature, index_of_x_wind, & diff --git a/physics/GWD/ugwpv1_gsldrag.meta b/physics/GWD/ugwpv1_gsldrag.meta index 8ed38a222..e1a201d57 100644 --- a/physics/GWD/ugwpv1_gsldrag.meta +++ b/physics/GWD/ugwpv1_gsldrag.meta @@ -412,7 +412,7 @@ [psl_gwd_dx_factor] standard_name = effective_grid_spacing_of_psl_gwd_suite long_name = multiplication of grid spacing - units = count + units = 1 dimensions = () type = real kind = kind_phys @@ -474,6 +474,14 @@ type = real kind = kind_phys intent = in +[alpha_fd] + standard_name = alpha_coefficient_for_turbulent_orographic_form_drag + long_name = alpha coefficient for Beljaars et al turbulent orographic form drag + units = none + dimensions = () + type = real + kind = kind_phys + intent = in [jdat] standard_name = date_and_time_of_forecast_in_united_states_order long_name = current forecast date and time diff --git a/physics/GWD/unified_ugwp.F90 b/physics/GWD/unified_ugwp.F90 index 832f34e47..5e47f2001 100644 --- a/physics/GWD/unified_ugwp.F90 +++ b/physics/GWD/unified_ugwp.F90 @@ -250,7 +250,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt dvsfc_ss,dusfc_fd,dvsfc_fd,dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,dtauy2d_bl, & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd,dudt_ngw,dvdt_ngw,dtdt_ngw, & br1,hpbl,vtype,slmsk, do_tofd, ldiag_ugwp, ugwp_seq_update, & - cdmbgwd, jdat, xlat, xlat_d, sinlat, coslat, area, & + cdmbgwd, alpha_fd, jdat, xlat, xlat_d, sinlat, coslat, area, & ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, & del, kpbl, dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, & tau_tofd, tau_mtb, tau_ogw, tau_ngw, zmtb, zlwb, zogw, & @@ -290,7 +290,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt real(kind=kind_phys), intent(in), dimension(:,:) :: del, ugrs, vgrs, tgrs, prsl, prslk, phil real(kind=kind_phys), intent(in), dimension(:,:) :: prsi, phii real(kind=kind_phys), intent(in), dimension(:,:) :: q1 - real(kind=kind_phys), intent(in) :: dtp, fhzero, cdmbgwd(:) + real(kind=kind_phys), intent(in) :: dtp, fhzero, cdmbgwd(:), alpha_fd integer, intent(in) :: jdat(:) logical, intent(in) :: do_tofd, ldiag_ugwp, ugwp_seq_update @@ -495,7 +495,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt if ( do_gsl_drag_ls_bl.or.do_gsl_drag_ss.or.do_gsl_drag_tofd ) then ! -if (do_gwd_opt_psl) then + if (do_gwd_opt_psl) then call drag_suite_psl(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, & tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, & @@ -506,14 +506,15 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & slmsk,br1,hpbl,vtype,con_g,con_cp,con_rd,con_rv, & con_fvirt,con_pi,lonr, & - cdmbgwd,me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + cdmbgwd,alpha_fd,me,master, & + lprnt,ipr,rdxzb,dx,gwd_opt, & do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & psl_gwd_dx_factor, & dtend, dtidx, index_of_process_orographic_gwd, & index_of_temperature, index_of_x_wind, & index_of_y_wind, ldiag3d, ldiag_ugwp, & ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) -else + else call drag_suite_run(im,levs,dvdt,dudt,dtdt,uwnd1,vwnd1, & tgrs,q1,kpbl,prsi,del,prsl,prslk,phii,phil,dtp, & kdt,hprime,oc,oa4,clx,varss,oc1ss,oa4ss, & @@ -524,13 +525,14 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, & slmsk,br1,hpbl,con_g,con_cp,con_rd,con_rv, & con_fvirt,con_pi,lonr, & - cdmbgwd,me,master,lprnt,ipr,rdxzb,dx,gwd_opt, & + cdmbgwd,alpha_fd,me,master, & + lprnt,ipr,rdxzb,dx,gwd_opt, & do_gsl_drag_ls_bl,do_gsl_drag_ss,do_gsl_drag_tofd, & dtend, dtidx, index_of_process_orographic_gwd, & index_of_temperature, index_of_x_wind, & index_of_y_wind, ldiag3d, ldiag_ugwp, & ugwp_seq_update, spp_wts_gwd, spp_gwd, errmsg, errflg) -endif + endif ! ! put zeros due to xy GSL-drag style: dtaux2d_bl,dtauy2d_bl,dtaux2d_ss.......dusfc_ms,dvsfc_ms ! diff --git a/physics/GWD/unified_ugwp.meta b/physics/GWD/unified_ugwp.meta index c885aa8d8..5dc347650 100644 --- a/physics/GWD/unified_ugwp.meta +++ b/physics/GWD/unified_ugwp.meta @@ -715,6 +715,14 @@ type = real kind = kind_phys intent = in +[alpha_fd] + standard_name = alpha_coefficient_for_turbulent_orographic_form_drag + long_name = alpha coefficient for Beljaars et al turbulent orographic form drag + units = none + dimensions = () + type = real + kind = kind_phys + intent = in [jdat] standard_name = date_and_time_of_forecast_in_united_states_order long_name = current forecast date and time @@ -1262,7 +1270,7 @@ [psl_gwd_dx_factor] standard_name = effective_grid_spacing_of_psl_gwd_suite long_name = multiplication of grid spacing - units = count + units = 1 dimensions = () type = real kind = kind_phys From f94b43df1d4481d88632d80ddce272f967996d23 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 31 May 2024 16:58:35 -0400 Subject: [PATCH 08/12] remove spurious metadata changes in drag_suite.meta --- physics/GWD/drag_suite.meta | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/physics/GWD/drag_suite.meta b/physics/GWD/drag_suite.meta index 088387cfd..29239bdd3 100644 --- a/physics/GWD/drag_suite.meta +++ b/physics/GWD/drag_suite.meta @@ -458,13 +458,6 @@ type = real kind = kind_phys intent = in -[vtype] - standard_name = vegetation_type_classification - long_name = vegetation type for lsm - units = index - dimensions = (horizontal_loop_extent) - type = integer - intent = in [g] standard_name = gravitational_acceleration long_name = gravitational acceleration @@ -608,14 +601,6 @@ dimensions = () type = logical intent = in -[psl_gwd_dx_factor] - standard_name = effective_grid_spacing_of_psl_gwd_suite - long_name = multiplication of grid spacing - units = 1 - dimensions = () - type = real - kind = kind_phys - intent = in [dtend] standard_name = cumulative_change_of_state_variables long_name = diagnostic tendencies for state variables From e59157849677fdaf66aa42b6c84ba6626becb70a Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 31 May 2024 19:59:50 -0400 Subject: [PATCH 09/12] remove some unnecessary variables and add some initialization of outputs --- physics/GWD/drag_suite.F90 | 77 +++++++++---------- physics/GWD/unified_ugwp.F90 | 16 +++- physics/GWD/unified_ugwp.meta | 24 ------ .../Land/Noahmp/module_sf_noahmplsm.F90 | 2 + physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 2 +- 5 files changed, 55 insertions(+), 66 deletions(-) diff --git a/physics/GWD/drag_suite.F90 b/physics/GWD/drag_suite.F90 index 91edfbec8..707cef597 100644 --- a/physics/GWD/drag_suite.F90 +++ b/physics/GWD/drag_suite.F90 @@ -513,7 +513,6 @@ subroutine drag_suite_run( & real(kind=kind_phys),parameter :: olmin = 1.0e-5 real(kind=kind_phys),parameter :: odmin = 0.1 real(kind=kind_phys),parameter :: odmax = 10. - real(kind=kind_phys),parameter :: erad = 6371.315e+3 integer :: komax(im) integer :: kblk real(kind=kind_phys) :: cd @@ -1552,7 +1551,7 @@ subroutine drag_suite_psl( & integer, parameter :: ims=1, kms=1, its=1, kts=1 real(kind=kind_phys), intent(in) :: fv, pi real(kind=kind_phys) :: rcl, cdmb - real(kind=kind_phys) :: g_inv, g_cp + real(kind=kind_phys) :: g_inv, g_cp, rd_inv real(kind=kind_phys), intent(inout) :: & & dudt(:,:),dvdt(:,:), & @@ -1736,7 +1735,6 @@ subroutine drag_suite_psl( & real(kind=kind_phys),parameter :: odmin = 0.1 real(kind=kind_phys),parameter :: odmax = 10. real(kind=kind_phys),parameter :: cdmin = 0.0 - real(kind=kind_phys),parameter :: erad = 6371.315e+3 integer :: komax(im),kbmax(im),kblk(im) real(kind=kind_phys) :: hmax(im) real(kind=kind_phys) :: cd @@ -1808,43 +1806,43 @@ subroutine drag_suite_psl( & endif enddo -! Remove ss_tapering -ss_taper(:) = 1. + ! Remove ss_tapering + ss_taper(:) = 1. -! SPP, if spp_gwd is 0, no perturbations are applied. -if ( spp_gwd==1 ) then - do i = its,im - var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1) - varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1) - varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1) - varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1) - enddo -else - do i = its,im - var_stoch(i) = var(i) - varss_stoch(i) = varss(i) - varmax_ss_stoch(i) = varmax_ss - varmax_fd_stoch(i) = varmax_fd - enddo -endif + ! SPP, if spp_gwd is 0, no perturbations are applied. + if ( spp_gwd==1 ) then + do i = its,im + var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1) + varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1) + varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1) + varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1) + enddo + else + do i = its,im + var_stoch(i) = var(i) + varss_stoch(i) = varss(i) + varmax_ss_stoch(i) = varmax_ss + varmax_fd_stoch(i) = varmax_fd + enddo + endif -!--- calculate length of grid for flow-blocking drag -! -do i=1,im - delx = dx(i) - dely = dx(i) - dxy4(i,1) = delx - dxy4(i,2) = dely - dxy4(i,3) = sqrt(delx*delx + dely*dely) - dxy4(i,4) = dxy4(i,3) - dxy4p(i,1) = dxy4(i,2) - dxy4p(i,2) = dxy4(i,1) - dxy4p(i,3) = dxy4(i,4) - dxy4p(i,4) = dxy4(i,3) - cleff(i) = psl_gwd_dx_factor*(delx+dely)*0.5 - cleff_ss(i) = 0.1 * max(dxmax_ss,dxy4(i,3)) -! cleff_ss(i) = cleff(i) ! consider ..... -enddo + !--- calculate length of grid for flow-blocking drag + ! + do i=1,im + delx = dx(i) + dely = dx(i) + dxy4(i,1) = delx + dxy4(i,2) = dely + dxy4(i,3) = sqrt(delx*delx + dely*dely) + dxy4(i,4) = dxy4(i,3) + dxy4p(i,1) = dxy4(i,2) + dxy4p(i,2) = dxy4(i,1) + dxy4p(i,3) = dxy4(i,4) + dxy4p(i,4) = dxy4(i,3) + cleff(i) = psl_gwd_dx_factor*(delx+dely)*0.5 + cleff_ss(i) = 0.1 * max(dxmax_ss,dxy4(i,3)) + ! cleff_ss(i) = cleff(i) ! consider ..... + enddo ! !-----initialize arrays ! @@ -1928,11 +1926,12 @@ subroutine drag_suite_psl( & kbmax(1:im) = 0 kblk(1:im) = 0 ! + rd_inv = 1./rd do k = kts,km do i = its,im vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k)) vtk(i,k) = vtj(i,k) / prslk(i,k) - ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3 + ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3 enddo enddo ! diff --git a/physics/GWD/unified_ugwp.F90 b/physics/GWD/unified_ugwp.F90 index 5e47f2001..adedeeb15 100644 --- a/physics/GWD/unified_ugwp.F90 +++ b/physics/GWD/unified_ugwp.F90 @@ -253,7 +253,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt cdmbgwd, alpha_fd, jdat, xlat, xlat_d, sinlat, coslat, area, & ugrs, vgrs, tgrs, q1, prsi, prsl, prslk, phii, phil, & del, kpbl, dusfcg, dvsfcg, gw_dudt, gw_dvdt, gw_dtdt, gw_kdis, & - tau_tofd, tau_mtb, tau_ogw, tau_ngw, zmtb, zlwb, zogw, & + tau_tofd, tau_mtb, tau_ogw, tau_ngw, & dudt_mtb, dudt_tms, du3dt_mtb, du3dt_ogw, du3dt_tms, & dudt, dvdt, dtdt, rdxzb, con_g, con_omega, con_pi, con_cp, con_rd, con_rv, & con_rerth, con_fvirt, rain, ntke, q_tke, dqdt_tke, lprnt, ipr, & @@ -311,7 +311,7 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt & slmsk(:) real(kind=kind_phys), intent(out), dimension(:) :: dusfcg, dvsfcg - real(kind=kind_phys), intent(out), dimension(:) :: zmtb, zlwb, zogw, rdxzb + real(kind=kind_phys), intent(out), dimension(:) :: rdxzb real(kind=kind_phys), intent(out), dimension(:) :: tau_mtb, tau_ogw, tau_tofd, tau_ngw real(kind=kind_phys), intent(out), dimension(:,:) :: gw_dudt, gw_dvdt, gw_dtdt, gw_kdis real(kind=kind_phys), intent(out), dimension(:,:) :: dudt_mtb, dudt_tms @@ -385,6 +385,18 @@ subroutine unified_ugwp_run(me, master, im, levs, ak,bk, ntrac, dtp, fhzero, kdt errflg = 0 + ! Initialize intent(out) variables in case they are not set below + dusfcg(:) = 0.0 + dvsfcg(:) = 0.0 + rdxzb(:) = 0.0 + tau_ngw(:) = 0.0 + gw_dudt(:,:) = 0.0 + gw_dvdt(:,:) = 0.0 + gw_dtdt(:,:) = 0.0 + gw_kdis(:,:) = 0.0 + dudt_mtb(:,:) = 0.0 + dudt_tms(:,:) = 0.0 + ! 1) ORO stationary GWs ! ------------------ diff --git a/physics/GWD/unified_ugwp.meta b/physics/GWD/unified_ugwp.meta index 5dc347650..814699375 100644 --- a/physics/GWD/unified_ugwp.meta +++ b/physics/GWD/unified_ugwp.meta @@ -938,30 +938,6 @@ type = real kind = kind_phys intent = out -[zmtb] - standard_name = height_of_mountain_blocking - long_name = height of mountain blocking drag - units = m - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out -[zlwb] - standard_name = height_of_low_level_wave_breaking - long_name = height of low level wave breaking - units = m - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out -[zogw] - standard_name = height_of_launch_level_of_orographic_gravity_wave - long_name = height of launch level of orographic gravity wave - units = m - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [dudt_mtb] standard_name = instantaneous_change_in_x_wind_due_to_mountain_blocking_drag long_name = instantaneous change in x wind due to mountain blocking drag diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index a3b466f46..97f11a1ba 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -2020,6 +2020,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in chuc = 0. chv2 = 0. rb = 0. + laisun = 0. + laisha = 0. cdmnv = 0.0 ezpdv = 0.0 diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index 4af3b93b6..0f62fb68e 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -158,7 +158,7 @@ subroutine noahmpdrv_run & sncovr1, qsurf, gflux, drain, evap, hflx, ep, runoff, & cmm, chh, evbs, evcw, sbsno, pah, ecan, etran, edir, snowc,& stm, snohf,smcwlt2, smcref2, wet1, t2mmp, q2mp,zvfun, & - ztmax, rca, errmsg, errflg, & + ztmax, rca, errmsg, errflg, & canopy_heat_storage_ccpp, & rainfall_ccpp, & sw_absorbed_total_ccpp, & From 1108aab4fa2284663f48c1552024526919787675 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 31 May 2024 20:16:51 -0400 Subject: [PATCH 10/12] fix some units in metadata files --- physics/GWD/drag_suite.meta | 2 +- physics/GWD/ugwpv1_gsldrag.meta | 2 +- physics/GWD/unified_ugwp.meta | 2 +- physics/SFC_Models/Land/Noahmp/noahmpdrv.meta | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/GWD/drag_suite.meta b/physics/GWD/drag_suite.meta index 29239bdd3..5413a5482 100644 --- a/physics/GWD/drag_suite.meta +++ b/physics/GWD/drag_suite.meta @@ -524,7 +524,7 @@ [alpha_fd] standard_name = alpha_coefficient_for_turbulent_orographic_form_drag long_name = alpha coefficient for Beljaars et al turbulent orographic form drag - units = none + units = 1 dimensions = () type = real kind = kind_phys diff --git a/physics/GWD/ugwpv1_gsldrag.meta b/physics/GWD/ugwpv1_gsldrag.meta index e1a201d57..6d9f5426b 100644 --- a/physics/GWD/ugwpv1_gsldrag.meta +++ b/physics/GWD/ugwpv1_gsldrag.meta @@ -477,7 +477,7 @@ [alpha_fd] standard_name = alpha_coefficient_for_turbulent_orographic_form_drag long_name = alpha coefficient for Beljaars et al turbulent orographic form drag - units = none + units = 1 dimensions = () type = real kind = kind_phys diff --git a/physics/GWD/unified_ugwp.meta b/physics/GWD/unified_ugwp.meta index 814699375..91f63f03e 100644 --- a/physics/GWD/unified_ugwp.meta +++ b/physics/GWD/unified_ugwp.meta @@ -718,7 +718,7 @@ [alpha_fd] standard_name = alpha_coefficient_for_turbulent_orographic_form_drag long_name = alpha coefficient for Beljaars et al turbulent orographic form drag - units = none + units = 1 dimensions = () type = real kind = kind_phys diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta index 594903122..7376ba4a2 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta @@ -665,7 +665,7 @@ [gwd_z0m_factor] standard_name = momentum_roughness_factor_in_turbulent_form_drag long_name = multiplication of orograhic standard deviation - units = count + units = 1 dimensions = () type = real kind = kind_phys From 00dc921754a699d361ee0f98b4f321ca021563a1 Mon Sep 17 00:00:00 2001 From: Qingfu Liu Date: Sat, 22 Jun 2024 13:11:48 -0400 Subject: [PATCH 11/12] retract the changes of the Noahmp model --- .../Land/Noahmp/module_sf_noahmplsm.F90 | 188 +++++++----------- physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 | 14 +- physics/SFC_Models/Land/Noahmp/noahmpdrv.meta | 16 -- 3 files changed, 78 insertions(+), 140 deletions(-) diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 index 97f11a1ba..2472f11ba 100644 --- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 +++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 @@ -90,7 +90,7 @@ module module_sf_noahmplsm integer :: opt_crs !< options for canopy stomatal resistance ! **1 -> ball-berry - ! 2 -> jarvis + ! 2 -> jarvis integer :: opt_btr !< options for soil moisture factor for stomatal resistance ! **1 -> noah (soil moisture) @@ -103,16 +103,16 @@ module module_sf_noahmplsm ! 3 -> original surface and subsurface runoff (free drainage) ! 4 -> bats surface and subsurface runoff (free drainage) ! 5 -> miguez-macho&fan groundwater scheme (miguez-macho et al. 2007 jgr; fan et al. 2007 jgr) - ! (needs further testing for public use) + ! (needs further testing for public use) integer :: opt_sfc !< options for surface layer drag coeff (ch & cm) ! **1 -> m-o - ! **2 -> original noah (chen97) - ! **3 -> myj consistent; 4->ysu consistent. mb: removed in v3.7 for further testing + ! **2 -> original noah (chen97) + ! **3 -> myj consistent; 4->ysu consistent. mb: removed in v3.7 for further testing integer :: opt_frz !< options for supercooled liquid water (or ice fraction) ! **1 -> no iteration (niu and yang, 2006 jhm) - ! 2 -> koren's iteration + ! 2 -> koren's iteration integer :: opt_inf !< options for frozen soil permeability ! **1 -> linear effects, more permeable (niu and yang, 2006, jhm) @@ -125,13 +125,13 @@ module module_sf_noahmplsm integer :: opt_alb !< options for ground snow surface albedo ! 1 -> bats - ! **2 -> class + ! **2 -> class integer :: opt_snf !< options for partitioning precipitation into rainfall & snowfall ! **1 -> jordan (1991) - ! 2 -> bats: when sfctmp sfctmp < tfrz - ! 4 -> use wrf microphysics output + ! 2 -> bats: when sfctmp sfctmp < tfrz + ! 4 -> use wrf microphysics output integer :: opt_tbot !< options for lower boundary condition of soil temperature ! 1 -> zero heat flux from bottom (zbot and tbot not used) @@ -139,20 +139,20 @@ module module_sf_noahmplsm integer :: opt_stc !< options for snow/soil temperature time scheme (only layer 1) ! **1 -> semi-implicit; flux top boundary condition - ! 2 -> full implicit (original noah); temperature top boundary condition + ! 2 -> full implicit (original noah); temperature top boundary condition ! 3 -> same as 1, but fsno for ts calculation (generally improves snow; v3.7) integer :: opt_rsf !< options for surface resistent to evaporation/sublimation ! **1 -> sakaguchi and zeng, 2009 - ! 2 -> sellers (1992) + ! 2 -> sellers (1992) ! 3 -> adjusted sellers to decrease rsurf for wet soil - ! 4 -> option 1 for non-snow; rsurf = rsurf_snow for snow (set in mptable); ad v3.8 + ! 4 -> option 1 for non-snow; rsurf = rsurf_snow for snow (set in mptable); ad v3.8 integer :: opt_soil !< options for defining soil properties ! **1 -> use input dominant soil texture - ! 2 -> use input soil texture that varies with depth + ! 2 -> use input soil texture that varies with depth ! 3 -> use soil composition (sand, clay, orgm) and pedotransfer functions (opt_pedo) - ! 4 -> use input soil properties (bexp_3d, smcmax_3d, etc.) + ! 4 -> use input soil properties (bexp_3d, smcmax_3d, etc.) integer :: opt_pedo !< options for pedotransfer functions (used when opt_soil = 3) ! **1 -> saxton and rawls (2006) @@ -423,9 +423,8 @@ subroutine noahmp_sflx (parameters, & smceq , & ! in : vegetation/soil characteristics sfctmp , sfcprs , psfc , uu , vv , q2, garea1 , & ! in : forcing qc , soldn , lwdn,thsfc_loc, prslkix,prsik1x,prslk1x,& ! in : forcing - varf , gwd_z0m_factor, & ! in : turbulent orographic form drag pblhx , iz0tlnd , itime ,psi_opt ,& - prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing + prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, & ! in : forcing tbot , co2air , o2air , foln , ficeold , zlvl , & ! in : forcing ep_1 , ep_2 , epsm1 , cp , & ! in : constants albold , sneqvo , & ! in/out : @@ -437,7 +436,7 @@ subroutine noahmp_sflx (parameters, & cm , ch , tauss , & ! in/out : grain , gdd , pgs , & ! in/out smcwtd ,deeprech , rech , ustarx , & ! in/out : - z0wrf , z0hwrf , ts , & ! out : + z0wrf , z0hwrf , ts , & ! out : fsa , fsr , fira , fsh , ssoil , fcev , & ! out : fgev , fctr , ecan , etran , edir , trad , & ! out : tgb , tgv , t2mv , t2mb , q2v , q2b , & ! out : @@ -446,9 +445,9 @@ subroutine noahmp_sflx (parameters, & qsnbot , ponding , ponding1, ponding2, rssun , rssha , & ! out : albd , albi , albsnd , albsni , & ! out : bgap , wgap , chv , chb , emissi , & ! out : - shg , shc , shb , evg , evb , ghv , & ! out : - ghb , irg , irc , irb , tr , evc , & ! out : - chleaf , chuc , chv2 , chb2 , fpice , pahv , & + shg , shc , shb , evg , evb , ghv , & ! out : + ghb , irg , irc , irb , tr , evc , & ! out : + chleaf , chuc , chv2 , chb2 , fpice , pahv , & pahg , pahb , pah , esnow , canhs , laisun , & laisha , rb , qsfcveg , qsfcbare & #ifdef CCPP @@ -494,8 +493,6 @@ subroutine noahmp_sflx (parameters, & real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< in exner function - real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] - real (kind=kind_phys) , intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height integer , intent(in) :: iz0tlnd !< z0t option integer , intent(in) :: itime !< @@ -822,7 +819,7 @@ subroutine noahmp_sflx (parameters, & canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + fwet ,cmc ) !out ! compute energy budget (momentum & energy fluxes and phase changes) @@ -835,9 +832,8 @@ subroutine noahmp_sflx (parameters, & fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc, jloc , & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in - varf ,gwd_z0m_factor , & !in pblhx ,iz0tlnd, itime ,psi_opt, ep_1, ep_2, epsm1,cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -858,7 +854,7 @@ subroutine noahmp_sflx (parameters, & fsrg ,rssun ,rssha ,albd ,albi ,albsnd,albsni, bgap ,wgap, tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb , & !out emissi ,pah ,canhs, & - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out qsfcveg = eah*ep_2/(sfcprs + epsm1*eah) qsfcbare = qsfc @@ -881,7 +877,7 @@ subroutine noahmp_sflx (parameters, & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc,jloc , smceq , & !in bdfall ,fp ,rain ,snow , & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -915,9 +911,9 @@ subroutine noahmp_sflx (parameters, & if (opt_crop == 1 .and. crop_active) then call carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julian , & !in dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in - soldn ,t2m , & !in + soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - lai ,sai ,gdd , & !inout + lai ,sai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out end if @@ -968,7 +964,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & prcpconv,prcpnonc ,prcpshcv,prcpsnow,prcpgrpl,prcphail , & soldn ,cosz ,thair ,qair , & eair ,rhoair ,qprecc ,qprecl ,solad , solai , & - swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) + swdown ,bdfall ,rain ,snow ,fp , fpice ,prcp ) ! -------------------------------------------------------------------------------------------------- ! re-process atmospheric forcing ! ---------------------------------------------------------------------- @@ -1041,7 +1037,7 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4) then qprecc = prcpconv + prcpshcv - qprecl = prcpnonc + qprecl = prcpnonc else qprecc = 0.10 * prcp ! should be from the atmospheric model qprecl = 0.90 * prcp ! should be from the atmospheric model @@ -1094,15 +1090,15 @@ subroutine atm (parameters,ep_2,epsm1,sfcprs ,sfctmp ,q2 , & if(opt_snf == 4 .or. opt_snf == 5) then prcp_frozen = prcpsnow + prcpgrpl + prcphail if(prcpnonc > 0. .and. prcp_frozen > 0.) then - fpice = min(1.0,prcp_frozen/prcpnonc) - fpice = max(0.0,fpice) + fpice = min(1.0,prcp_frozen/prcpnonc) + fpice = max(0.0,fpice) if(opt_snf==4) bdfall = bdfall*(prcpsnow/prcp_frozen) + rho_grpl*(prcpgrpl/prcp_frozen) + & rho_hail*(prcphail/prcp_frozen) if(opt_snf==5) bdfall = parameters%prcpiceden - else - fpice = 0.0 + else + fpice = 0.0 endif - + endif rain = prcp * (1.-fpice) @@ -1237,8 +1233,8 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv bdfall ,rain ,snow ,fp , & !in canliq ,canice ,tv ,sfctmp ,tg , & !in qintr ,qdripr ,qthror ,qints ,qdrips ,qthros , & !out - pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out - fwet ,cmc ) !out + pahv ,pahg ,pahb ,qrain ,qsnow ,snowhin, & !out + fwet ,cmc ) !out ! ------------------------ code history ------------------------------ ! michael barlage: oct 2013 - split canwater to calculate precip movement for @@ -1340,10 +1336,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qintr = 0. qdripr = 0. qthror = rain - if(canliq > 0.) then ! for case of canopy getting buried - qdripr = qdripr + canliq/dt - canliq = 0.0 - end if + if(canliq > 0.) then ! for case of canopy getting buried + qdripr = qdripr + canliq/dt + canliq = 0.0 + end if end if ! heat transported by liquid water @@ -1366,8 +1362,8 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qints = max(qints, 0.) ft = max(0.0,(tv - 270.15) / 1.87e5) fv = sqrt(uu*uu + vv*vv) / 1.56e5 -! mb: changed below to reflect the rain assumption that all precip gets intercepted - icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt + ! mb: changed below to reflect the rain assumption that all precip gets intercepted + icedrip = max(0.,canice) * (fv+ft) !mb: removed /dt qdrips = (fveg * snow - qints) + icedrip qthros = (1.0-fveg) * snow canice= max(0.,canice + (qints - icedrip)*dt) @@ -1375,10 +1371,10 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv qints = 0. qdrips = 0. qthros = snow - if(canice > 0.) then ! for case of canopy getting buried - qdrips = qdrips + canice/dt - canice = 0.0 - end if + if(canice > 0.) then ! for case of canopy getting buried + qdrips = qdrips + canice/dt + canice = 0.0 + end if endif ! print*, "precip_heat canopy through:",3600.0*(fveg * snow - qints) ! print*, "precip_heat canopy drip:",3600.0*max(0.,canice) * (fv+ft) @@ -1408,13 +1404,13 @@ subroutine precip_heat (parameters,iloc ,jloc ,vegtyp ,dt ,uu ,vv if (fveg > 0.0 .and. fveg < 1.0) then pahg = pahg / fveg ! these will be multiplied by fraction later - pahb = pahb / (1.0-fveg) + pahb = pahb / (1.0-fveg) elseif (fveg <= 0.0) then pahb = pahg + pahb ! for case of canopy getting buried pahg = 0.0 - pahv = 0.0 + pahv = 0.0 elseif (fveg >= 1.0) then - pahb = 0.0 + pahb = 0.0 end if pahv = max(pahv,-20.0) ! put some artificial limits here for stability @@ -1680,9 +1676,8 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in fveg ,shdfac, pahv ,pahg ,pahb , & !in qsnow ,dzsnso ,lat ,canliq ,canice ,iloc , jloc, & !in thsfc_loc, prslkix,prsik1x,prslk1x,garea1, & !in - varf ,gwd_z0m_factor , & !in pblhx , iz0tlnd, itime,psi_opt,ep_1, ep_2, epsm1, cp, & - z0wrf ,z0hwrf , & !out + z0wrf ,z0hwrf , & !out imelt ,snicev ,snliqv ,epore ,t2m ,fsno , & !out sav ,sag ,qmelt ,fsa ,fsr ,taux , & !out tauy ,fira ,fsh ,fcev ,fgev ,fctr , & !out @@ -1702,7 +1697,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2mv ,t2mb ,fsrv , & fsrg ,rssun ,rssha ,albd ,albi,albsnd ,albsni,bgap ,wgap,tgv,tgb,& q1 ,q2v ,q2b ,q2e ,chv ,chb, emissi,pah,canhs,& - shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out + shg,shc,shb,evg,evb,ghv,ghb,irg,irc,irb,tr,evc,chleaf,chuc,chv2,chb2 ) !out !jref:end ! -------------------------------------------------------------------------------------------------- @@ -1762,8 +1757,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) , intent(in) :: prslk1x !< in exner function real (kind=kind_phys) , intent(in) :: garea1 !< - real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] - real (kind=kind_phys) , intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -2218,19 +2211,19 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in if (tv .gt. tfrz) then ! barlage: add distinction between ground and latheav = hvap ! vegetation in v3.6 - frozen_canopy = .false. + frozen_canopy = .false. else latheav = hsub - frozen_canopy = .true. + frozen_canopy = .true. end if gammav = cpair*sfcprs/(ep_2*latheav) if (tg .gt. tfrz) then latheag = hvap - frozen_ground = .false. + frozen_ground = .false. else latheag = hsub - frozen_ground = .true. + frozen_ground = .true. end if gammag = cpair*sfcprs/(ep_2*latheag) @@ -2259,7 +2252,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in - varf , gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tgv ,cmv, ustarx , & !inout #ifdef CCPP @@ -2274,7 +2266,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in !jref:start qc ,qsfc ,psfc , & !in q2v ,chv2 ,chleaf ,chuc , & - rb) !out + rb) !out ! new coupling code @@ -2298,7 +2290,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in emg ,stc ,df ,rsurf ,latheag , & !in gammag ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in - varf ,gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & #ifdef CCPP tgb ,cmb ,chb, ustarx,errmsg ,errflg , & !inout @@ -2343,7 +2334,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in ssoil = fveg * ghv + (1.0 - fveg) * ghb fcev = evc fctr = tr - pah = fveg * pahg + (1.0 - fveg) * pahb + pahv + pah = fveg * pahg + (1.0 - fveg) * pahb + pahv tg = fveg * tgv + (1.0 - fveg) * tgb t2m = fveg * t2mv + (1.0 - fveg) * t2mb ts = fveg * tah + (1.0 - fveg) * tgb @@ -2373,7 +2364,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in t2m = t2mb fcev = 0. fctr = 0. - pah = pahb + pah = pahb ts = tg cm = cmb ch = chb @@ -3544,7 +3535,7 @@ subroutine twostream (parameters,ib ,ic ,vegtyp ,cosz ,vai , & ! kopen = 1.0 else if(opt_rad == 1) then - denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) + denfveg = -log(max(1.0-fveg,0.01))/(pai*parameters%rc**2) hd = parameters%hvt - parameters%hvb bb = 0.5 * hd thetap = atan(bb/parameters%rc * tan(acos(max(0.01,cosz))) ) @@ -3712,7 +3703,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & foln ,co2air ,o2air ,btran ,sfcprs , & !in rhsur ,iloc ,jloc ,q2 ,pahv ,pahg , & !in thsfc_loc, prslkix,prsik1x,prslk1x, garea1, & !in - varf ,gwd_z0m_factor, & !in pblhx ,iz0tlnd ,itime ,psi_opt ,ep_1, ep_2, epsm1, cp, & eah ,tah ,tv ,tg ,cm,ustarx,& !inout #ifdef CCPP @@ -3726,7 +3716,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & csigmaf1, & !out qc ,qsfc ,psfc , & !in q2v ,cah2 ,chleaf ,chuc , & !inout - rb) !out + rb) !out ! -------------------------------------------------------------------------------------------------- ! use newton-raphson iteration to solve for vegetation (tv) and @@ -3763,8 +3753,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & real (kind=kind_phys), intent(in) :: dt !< time step (s) real (kind=kind_phys), intent(in) :: fsno !< snow fraction - real (kind=kind_phys) , intent(in) :: varf !< sub-grid orography standard deviation [m] - real (kind=kind_phys) , intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys) , intent(in) :: pblhx !< pbl height real (kind=kind_phys) , intent(in) :: ep_1 !< real (kind=kind_phys) , intent(in) :: ep_2 !< @@ -4137,7 +4125,6 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tah ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd ,snowh ,shdfac ,garea1 , & !in - varf,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -4260,11 +4247,11 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & shc = fveg*rhoair*cpair*cvh * ( tv-tah) evc = fveg*rhoair*cpair*cew * (estv-eah) / gammav ! barlage: change to v in v3.6 tr = fveg*rhoair*cpair*ctw * (estv-eah) / gammav - if (tv > tfrz) then + if (tv > tfrz) then evc = min(canliq*latheav/dt,evc) ! barlage: add if block for canice in v3.6 - else + else evc = min(canice*latheav/dt,evc) - end if + end if ! canopy heat capacity hcv = fveg*(parameters%cbiom*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice) !j/m2/k @@ -4447,7 +4434,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & emg ,stc ,df ,rsurf ,lathea , & !in gamma ,rhsur ,iloc ,jloc ,q2 ,pahb , & !in thsfc_loc, prslkix,prsik1x,prslk1x,vegtyp,fveg,shdfac,garea1, & !in - varf,gwd_z0m_factor, & !in pblhx , iz0tlnd , itime ,psi_opt,ep_1,ep_2,epsm1,cp ,& #ifdef CCPP tgb ,cm ,ch,ustarx,errmsg ,errflg , & !inout @@ -4502,8 +4488,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & real (kind=kind_phys), intent(in) :: rhsur !< raltive humidity in surface soil/snow air space (-) real (kind=kind_phys), intent(in) :: fsno !< snow fraction - real (kind=kind_phys), intent(in) :: varf !< sub-grid orography standard deviation [m] - real (kind=kind_phys), intent(in) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys), intent(in) :: pblhx !< pbl height (m) real (kind=kind_phys), intent(in) :: ep_1 !< real (kind=kind_phys), intent(in) :: ep_2 !< @@ -4747,7 +4731,6 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h, zpd,snowh ,shdfac ,garea1 , & !in - varf,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5454,7 +5437,6 @@ end subroutine sfcdif2 subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in z0h,zpd ,snowh ,fveg ,garea1 , & !in - varf ,gwd_z0m_factor, & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout fv ,cm ,ch ) !out @@ -5484,8 +5466,6 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur real (kind=kind_phys), intent(in ) :: snowh !< snow depth [m] real (kind=kind_phys), intent(in ) :: fveg !< fractional vegetation cover real (kind=kind_phys), intent(in ) :: garea1 !< grid area [km2] - real (kind=kind_phys), intent(in ) :: varf !< sub-grid orography standard deviation [m] - real (kind=kind_phys), intent(in ) :: gwd_z0m_factor !< turbulent orographic form drag roughness length factor real (kind=kind_phys), intent(inout) :: ustarx !< friction velocity [m/s] real (kind=kind_phys), intent(inout) :: fm !< momentum stability correction, weighted by prior iters real (kind=kind_phys), intent(inout) :: fh !< sen heat stability correction, weighted by prior iters @@ -5540,7 +5520,7 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur endif call gfs_stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & - varf,gwd_z0m_factor,rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) + rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) end subroutine sfcdif3 @@ -5550,7 +5530,6 @@ subroutine gfs_stability & ! --- inputs: ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & thsfc_loc, & - varf,gwd_z0m_factor, & ! --- outputs: rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) @@ -5570,8 +5549,6 @@ subroutine gfs_stability & real(kind=kind_phys), intent(in) :: ztmax ! thermal roughness length real(kind=kind_phys), intent(in) :: tvs ! surface virtual temperature real(kind=kind_phys), intent(in) :: grav ! local gravity -real(kind=kind_phys), intent(in) :: varf ! turbulent orographic standard deviation [m] -real(kind=kind_phys), intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor logical, intent(in) :: thsfc_loc ! use local theta reference flag real(kind=kind_phys), intent(out) :: rb ! bulk richardson number [-] @@ -5630,10 +5607,6 @@ subroutine gfs_stability & real(kind=kind_phys) :: zolmax real(kind=kind_phys) xkzo -! -real(kind=kind_phys) :: cf, zf, fri, fphim -real(kind=kind_phys), parameter :: varf_min = 50.0 -real(kind=kind_phys), parameter :: varf_max = 500.0 z1i = one / z1 ! inverse of model height @@ -5754,19 +5727,6 @@ subroutine gfs_stability & endif ! end of if (dtv >= 0 ) then loop ! -! compute the effect of orographic form drag : koo et al. (2018,jgr, eqs. 1 and 2) -! - if ( varf.gt.varf_min .and. gwd_z0m_factor.gt.0.0_kind_phys ) then - tem1 = min(varf, varf_max) - zf = tem1 * gwd_z0m_factor ! effective roughness length - zf = min(zf, z1) - fri = min( max( 1.0_kind_phys-rb,0.0_kind_phys ), 1.0_kind_phys) ! a function of ther bulk Richardson number - fphim = log( ( z1 + zf) / zf ) ! integrated profile function for momentum - cf = ca*ca / (fphim*fphim) * fri ! exchange coefficient for momentum - else - cf = 0. - endif -! ! finish the exchange coefficient computation to provide fm and fh ! fm = fm - pm ! phi_m @@ -5778,7 +5738,7 @@ subroutine gfs_stability & tem1 = 0.00001_kind_phys/z1 ! minimum exhange coef (?) cm = max(cm, tem1) ch = max(ch, tem1) - stress = (cm + cf)* wind * wind + stress = cm * wind * wind ! surface stress = Cm*U*U ustar = sqrt(stress) ! friction velocity return @@ -5927,7 +5887,7 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, else z0h_out = z0m_out * 0.01 endif - + elseif (opt_trs == blumel99) then reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c @@ -7060,7 +7020,7 @@ subroutine water (parameters,vegtyp ,nsnow ,nsoil ,imelt ,dt ,uu , & esai ,sfctmp ,qvap ,qdew ,zsoil ,btrani , & !in ficeold,ponding,tg ,ist ,fveg ,iloc ,jloc ,smceq , & !in bdfall ,fp ,rain ,snow, & !in mb/an: v3.7 - qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb + qsnow ,qrain ,snowhin,latheav,latheag,frozen_canopy,frozen_ground, & !in mb isnow ,canliq ,canice ,tv ,snowh ,sneqv , & !inout snice ,snliq ,stc ,zsnso ,sh2o ,smc , & !inout sice ,zwt ,wa ,wt ,dzsnso ,wslake , & !inout @@ -7689,19 +7649,19 @@ subroutine combine (parameters,nsnow ,nsoil ,iloc ,jloc , & !in snice(j-1) = snice(j-1) + snice(j) dzsnso(j-1) = dzsnso(j-1) + dzsnso(j) else - if(snice(j) >= 0.) then + if(snice(j) >= 0.) then ponding1 = snliq(j) ! isnow will get set to zero below; ponding1 will get sneqv = snice(j) ! added to ponding from phasechange ponding should be snowh = dzsnso(j) ! zero here because it was calculated for thin snow - else ! snice over-sublimated earlier - ponding1 = snliq(j) + snice(j) - if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil - sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) + else ! snice over-sublimated earlier + ponding1 = snliq(j) + snice(j) + if(ponding1 < 0.) then ! if snice and snliq sublimates remove from soil + sice(1) = max(0.0,sice(1)+ponding1/(dzsnso(1)*1000.)) ponding1 = 0.0 - end if + end if sneqv = 0.0 snowh = 0.0 - end if + end if snliq(j) = 0.0 snice(j) = 0.0 dzsnso(j) = 0.0 @@ -9790,7 +9750,7 @@ subroutine carbon_crop (parameters,nsnow ,nsoil ,vegtyp ,dt ,zsoil ,julia dzsnso ,stc ,smc ,tv ,psn ,foln ,btran , & !in soldn ,t2m , & !in lfmass ,rtmass ,stmass ,wood ,stblcp ,fastcp ,grain , & !inout - xlai ,xsai ,gdd , & !inout + xlai ,xsai ,gdd , & !inout gpp ,npp ,nee ,autors ,heters ,totsc ,totlb, pgs ) !out ! ------------------------------------------------------------------------------------------ ! initial crop version created by xing liu @@ -10469,7 +10429,7 @@ end subroutine psn_crop !! subroutine noahmp_options(idveg , iopt_crs , iopt_btr , iopt_run , iopt_sfc , iopt_frz , & iopt_inf, iopt_rad , iopt_alb , iopt_snf , iopt_tbot, iopt_stc , & - iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & + iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & iopt_z0m ) implicit none diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index 0f62fb68e..6234f46d9 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -137,7 +137,6 @@ subroutine noahmpdrv_run & rainn_mp, rainc_mp, snow_mp, graupel_mp, ice_mp, rhonewsn1,& con_hvap, con_cp, con_jcal, rhoh2o, con_eps, con_epsm1, & con_fvirt, con_rd, con_hfus, thsfc_loc, cpllnd, cpllnd2atm,& - mntvar,gwd_z0m_factor, & ! --- in/outs: weasd, snwdph, tskin, tprcp, srflag, smc, stc, slc, & @@ -298,8 +297,6 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(in), optional :: snow_mp ! microphysics snow [mm] real(kind=kind_phys), dimension(:) , intent(in), optional :: graupel_mp ! microphysics graupel [mm] real(kind=kind_phys), dimension(:) , intent(in), optional :: ice_mp ! microphysics ice/hail [mm] - real(kind=kind_phys), dimension(:,:) , intent(in) :: mntvar ! subgrid orographic statistics data - real(kind=kind_phys) , intent(in) :: gwd_z0m_factor ! turbulent orographic form drag roughness length factor real(kind=kind_phys), dimension(:) , intent(in) :: rhonewsn1 ! precipitation ice density (kg/m^3) real(kind=kind_phys) , intent(in) :: con_hvap ! latent heat condensation [J/kg] real(kind=kind_phys) , intent(in) :: con_cp ! specific heat air [J/kg/K] @@ -404,6 +401,7 @@ subroutine noahmpdrv_run & real(kind=kind_phys), dimension(:) , intent(out) :: zvfun ! real(kind=kind_phys), dimension(:) , intent(out) :: ztmax ! thermal roughness length real(kind=kind_phys), dimension(:) , intent(out), optional :: rca ! total canopy/stomatal resistance (s/m) + character(len=*) , intent(out) :: errmsg integer , intent(out) :: errflg @@ -627,7 +625,7 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: canopy_heat_storage ! out | within-canopy heat [W/m2] real (kind=kind_phys) :: spec_humid_sfc_veg ! out | surface specific humidty over vegetation [kg/kg] real (kind=kind_phys) :: spec_humid_sfc_bare ! out | surface specific humidty over bare soil [kg/kg] - + real (kind=kind_phys) :: ustarx ! inout |surface friction velocity real (kind=kind_phys) :: prslkix ! in exner function real (kind=kind_phys) :: prsik1x ! in exner function @@ -637,7 +635,6 @@ subroutine noahmpdrv_run & real (kind=kind_phys) :: cq2 real (kind=kind_phys) :: qfx real (kind=kind_phys) :: wspd1 ! wind speed with all components - real (kind=kind_phys) :: varf_grid ! sub-grid orography standard deviation real (kind=kind_phys) :: pblhx ! height of pbl integer :: mnice @@ -743,7 +740,6 @@ subroutine noahmpdrv_run & area_grid = garea(i) pblhx = pblh(i) - varf_grid = mntvar(i,15) prslkix = prslki(i) prsik1x = prsik1(i) @@ -975,7 +971,6 @@ subroutine noahmpdrv_run & spec_humidity_forcing ,area_grid ,cloud_water_forcing , & sw_radiation_forcing ,radiation_lw_forcing ,thsfc_loc , & prslkix ,prsik1x ,prslk1x , & - varf_grid ,gwd_z0m_factor , & pblhx ,iz0tlnd ,itime , & psi_opt , & precip_convective , & @@ -1067,7 +1062,7 @@ subroutine noahmpdrv_run & chxy (i) = ch_noahmp zorl (i) = z0_total * 100.0 ! convert to cm ztmax (i) = z0h_total - + !LAI-scale canopy resistance based on weighted sunlit shaded fraction if(rs_sunlit .le. 0.0 .or. rs_shaded .le. 0.0 .or. & lai_sunlit .eq. 0.0 .or. lai_shaded .eq. 0.0) then @@ -1077,7 +1072,7 @@ subroutine noahmpdrv_run & ((1.0/(rs_shaded+leaf_air_resistance))*lai_shaded)) rca(i) = max((1.0/rca(i)),parameters%rsmin) !resistance end if - + smc (i,:) = soil_moisture_vol slc (i,:) = soil_liquid_vol snowxy (i) = float(snow_levels) @@ -1220,7 +1215,6 @@ subroutine noahmpdrv_run & call gfs_stability & (zf(i), zvfun(i), gdx, virtual_temperature, vptemp,wind(i), z0_total, z0h_total, & tvs1, con_g, thsfc_loc, & - varf_grid, gwd_z0m_factor, & rb1(i), fm1(i), fh1(i), fm101(i), fh21(i), cm(i), ch(i), stress1(i), ustar1(i)) rmol1(i) = undefined !not used in GFS sfcdif -> to satsify output diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta index 7376ba4a2..753550016 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta @@ -654,22 +654,6 @@ dimensions = () type = logical intent = in -[mntvar] - standard_name = statistical_measures_of_subgrid_orography_collection_array - long_name = array of statistical measures of subgrid height_above_mean_sea_level - units = mixed - dimensions = (horizontal_loop_extent,number_of_statistical_measures_of_subgrid_orography) - type = real - kind = kind_phys - intent = in -[gwd_z0m_factor] - standard_name = momentum_roughness_factor_in_turbulent_form_drag - long_name = multiplication of orograhic standard deviation - units = 1 - dimensions = () - type = real - kind = kind_phys - intent = in [weasd] standard_name = water_equivalent_accumulated_snow_depth_over_land long_name = water equiv of acc snow depth over land From c6dec90f4cd8e9f7ac5772cf51e9df6cd5f59cdd Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 28 Jun 2024 10:09:15 -0400 Subject: [PATCH 12/12] add optional attribute as necessary to arguments in drag_suite_psl subroutine; change upper bound of loop to prevent array bound excursion --- physics/GWD/drag_suite.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/physics/GWD/drag_suite.F90 b/physics/GWD/drag_suite.F90 index 707cef597..5c2bf6c2c 100644 --- a/physics/GWD/drag_suite.F90 +++ b/physics/GWD/drag_suite.F90 @@ -1542,7 +1542,7 @@ subroutine drag_suite_psl( & integer, intent(in) :: KPBL(:) real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, & & cdmbgwd(:), alpha_fd - real(kind=kind_phys), intent(inout) :: dtend(:,:,:) + real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:) logical, intent(in) :: ldiag3d integer, intent(in) :: dtidx(:,:), index_of_temperature, & & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind @@ -1567,7 +1567,7 @@ subroutine drag_suite_psl( & real(kind=kind_phys), intent(in) :: var(:),oc1(:), & & oa4(:,:),ol4(:,:), & & dx(:) - real(kind=kind_phys), intent(in) :: varss(:),oc1ss(:), & + real(kind=kind_phys), intent(in), optional :: varss(:),oc1ss(:), & & oa4ss(:,:),ol4ss(:,:) real(kind=kind_phys), intent(in) :: THETA(:),SIGMA(:), & & GAMMA(:),ELVMAX(:) @@ -1589,7 +1589,7 @@ subroutine drag_suite_psl( & !SPP real(kind=kind_phys), dimension(im) :: var_stoch, varss_stoch, & varmax_ss_stoch, varmax_fd_stoch - real(kind=kind_phys), intent(in) :: spp_wts_gwd(:,:) + real(kind=kind_phys), intent(in), optional :: spp_wts_gwd(:,:) integer, intent(in) :: spp_gwd real(kind=kind_phys), dimension(im) :: rstoch @@ -1598,12 +1598,12 @@ subroutine drag_suite_psl( & real(kind=kind_phys), intent(out) :: & & dusfc(:), dvsfc(:) !Output (optional): - real(kind=kind_phys), intent(out) :: & + real(kind=kind_phys), intent(out), optional :: & & dusfc_ls(:),dvsfc_ls(:), & & dusfc_bl(:),dvsfc_bl(:), & & dusfc_ss(:),dvsfc_ss(:), & & dusfc_fd(:),dvsfc_fd(:) - real(kind=kind_phys), intent(out) :: & + real(kind=kind_phys), intent(out), optional :: & & dtaux2d_ls(:,:),dtauy2d_ls(:,:), & & dtaux2d_bl(:,:),dtauy2d_bl(:,:), & & dtaux2d_ss(:,:),dtauy2d_ss(:,:), & @@ -2551,7 +2551,7 @@ subroutine drag_suite_psl( & ! apply limiter to mesosphere drag, reduce the drag by density factor 10-3 ! prevent wind reversal... ! - do k = kpblmax,km + do k = kpblmax,km-1 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff) then denfac = min(ro(i,k)/pcutoff_den,1.) dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k) &