From f896694d06ae6b2f93e3285290ca89fe29117c30 Mon Sep 17 00:00:00 2001 From: Jun Wang Date: Thu, 11 Jun 2020 02:24:48 +0000 Subject: [PATCH 1/4] merge GFSv16 physics update with rad bug fix, reverting changes in sfc_diff, tuning in samfdeepcnv.f, update veritcal mixing in satmedmfvdifq.F --- physics/radsw_datatb.f | 2 +- physics/samfdeepcnv.f | 16 +++------------- physics/satmedmfvdifq.F | 37 ++++++++++++++++++++++++++++--------- physics/satmedmfvdifq.meta | 17 +++++++++++++++++ physics/sfc_diff.f | 24 ++++++++---------------- 5 files changed, 57 insertions(+), 39 deletions(-) diff --git a/physics/radsw_datatb.f b/physics/radsw_datatb.f index 3cc9e2d82..6d88f1989 100644 --- a/physics/radsw_datatb.f +++ b/physics/radsw_datatb.f @@ -2551,7 +2551,7 @@ module module_radsw_sflux ! !> band index (3rd index in array sfluxref described below) integer, dimension(nblow:nbhgh), public :: ibx - data layreffr/ 18,30, 6, 3, 3, 8, 2, 6, 1, 2, 0,32,58,49 / + data layreffr/ 18,30, 6, 3, 3, 8, 2, 6, 1, 2, 0,32,42,49 / data ix1 / 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 3, 0 / data ix2 / 5, 2, 5, 2, 0, 2, 6, 0, 6, 0, 0, 0, 6, 0 / data ibx / 1, 1, 1, 2, 2, 3, 4, 3, 5, 4, 5, 6, 2, 7 / diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 9ec9ba7f3..677a2cee1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -222,7 +222,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & parameter(clamd=0.03,tkemx=0.65,tkemn=0.05) parameter(dtke=tkemx-tkemn) parameter(dbeta=0.1) - parameter(cthk=200.,dthk=25.) + parameter(cthk=150.,dthk=25.) parameter(cinpcrmx=180.,cinpcrmn=120.) ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) @@ -1239,23 +1239,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c specify upper limit of mass flux at cloud base c !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. - if(hwrf_samfdeep) then - do i = 1, im + do i = 1, im if(cnvflg(i)) then k = kbcon(i) dp = 1000. * del(i,k) xmbmax(i) = dp / (grav * dt2) endif - enddo - else - do i = 1, im - if(cnvflg(i)) then - k = kbcon(i) - dp = 1000. * del(i,k) - xmbmax(i) = dp / (2. * grav * dt2) - endif - enddo - endif + enddo c c compute cloud moisture property and precipitation c diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index d465b7c5e..0f4aa5103 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -59,8 +59,8 @@ end subroutine satmedmfvdifq_finalize !! @{ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & - & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea, & - & psk,rbsoil,zorl,u10m,v10m,fm,fh, & + & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea,islimsk, & + & snwdph_lnd,psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, & @@ -75,6 +75,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & !---------------------------------------------------------------------- integer, intent(in) :: ix, im, km, ntrac, ntcw, ntiw, ntke integer, intent(in) :: kinver(im) + integer, intent(in) :: islimsk(im) integer, intent(out) :: kpbl(im) ! real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, & @@ -88,6 +89,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & t1(ix,km), q1(ix,km,ntrac), & & swh(ix,km), hlw(ix,km), & & xmu(im), garea(im), & + & snwdph_lnd(im), & & psk(ix), rbsoil(im), & & zorl(im), tsea(im), & & u10m(im), v10m(im), & @@ -201,6 +203,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & zlup, zldn, bsum, & tem, tem1, tem2, & ptem, ptem0, ptem1, ptem2 +! + real(kind=kind_phys) xkzm_mp, xkzm_hp ! real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck ! @@ -212,7 +216,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) parameter(vk=0.4,rimin=-100.) parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) - parameter(rlmn=30.,rlmn1=5.,rlmn2=15.) + parameter(rlmn=30.,rlmn1=5.,rlmn2=10.) parameter(rlmx=300.,elmx=300.) parameter(prmin=0.25,prmax=4.0) parameter(pr0=1.0,prtke=1.0,prscu=0.67) @@ -222,7 +226,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & parameter(aphi5=5.,aphi16=16.) parameter(elmfac=1.0,elefac=1.0,cql=100.) parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=5000.) - parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15) + parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.1) parameter(h1=0.33333333) parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15) parameter(ce0=0.4) @@ -317,16 +321,31 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & !! \n xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) do i=1,im + xkzm_mp = xkzm_m + xkzm_hp = xkzm_h +! + if( islimsk(i) == 1 .and. snwdph_lnd(i) > 10.0 ) then ! over land + if (rbsoil(i) > 0. .and. rbsoil(i) <= 0.25) then + xkzm_mp = xkzm_m * (1.0 - rbsoil(i)/0.25)**2 + + & 0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2) + xkzm_hp = xkzm_h * (1.0 - rbsoil(i)/0.25)**2 + + & 0.1 * (1.0 - (1.0-rbsoil(i)/0.25)**2) + else if (rbsoil(i) > 0.25) then + xkzm_mp = 0.1 + xkzm_hp = 0.1 + endif + endif +! kx1(i) = 1 tx1(i) = 1.0 / prsi(i,1) tx2(i) = tx1(i) if(gdx(i) >= xkgdx) then - xkzm_hx(i) = xkzm_h - xkzm_mx(i) = xkzm_m + xkzm_hx(i) = xkzm_hp + xkzm_mx(i) = xkzm_mp else tem = 1. / (xkgdx - 5.) - tem1 = (xkzm_h - 0.01) * tem - tem2 = (xkzm_m - 0.01) * tem + tem1 = (xkzm_hp - 0.01) * tem + tem2 = (xkzm_mp - 0.01) * tem ptem = gdx(i) - 5. xkzm_hx(i) = 0.01 + tem1 * ptem xkzm_mx(i) = 0.01 + tem2 * ptem @@ -833,7 +852,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k) ! if(tem1 > 1.e-5) then tem1 = tvx(i,k+1)-tvx(i,k) - if(tem1 > 0.) then + if(tem1 > 0. .and. islimsk(i) /= 1) then xkzo(i,k) = min(xkzo(i,k), xkzinv) xkzmo(i,k) = min(xkzmo(i,k), xkzinv) rlmnz(i,k) = min(rlmnz(i,k), rlmn2) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 26667a627..f8f0c1918 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -284,6 +284,23 @@ kind = kind_phys intent = in optional = F +[islimsk] + standard_name = sea_land_ice_mask + long_name = sea/land/ice mask (=0/1/2) + units = flag + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[snwdph_lnd] + standard_name = surface_snow_thickness_water_equivalent_over_land + long_name = water equivalent snow depth over land + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F [psk] standard_name = dimensionless_exner_function_at_lowest_model_interface long_name = dimensionless Exner function at the surface interface diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 2a723e70c..e55ec90d7 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -220,15 +220,11 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) z0max = max(z0max, 1.0e-6) ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil -! czilc = 0.8 + czilc = 0.8 -! tem1 = 1.0 - sigmaf(i) -! ztmax = z0max*exp( - tem1*tem1 -! & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) -! - czilc = 10.0 ** (- 4. * z0max) ! Trier et al. (2011, WAF) - ztmax = z0max * exp( - czilc * ca - & * 258.2 * sqrt(ustar_lnd(i)*z0max) ) + tem1 = 1.0 - sigmaf(i) + ztmax = z0max*exp( - tem1*tem1 + & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) ! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land @@ -265,15 +261,11 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height ! dependance of czil -! czilc = 0.8 - -! tem1 = 1.0 - sigmaf(i) -! ztmax = z0max*exp( - tem1*tem1 -! & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) - czilc = 10.0 ** (- 4. * z0max) ! Trier et al. (2011, WAF) - ztmax = z0max * exp( - czilc * ca - & * 258.2 * sqrt(ustar_ice(i)*z0max) ) + czilc = 0.8 + tem1 = 1.0 - sigmaf(i) + ztmax = z0max*exp( - tem1*tem1 + & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) ztmax = max(ztmax, 1.0e-6) ! call stability From 56bca411c1882b0c53fe86cef2332a795eed7f72 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Thu, 11 Jun 2020 12:05:59 +0000 Subject: [PATCH 2/4] update sfc_diff.f --- physics/sfc_diff.f | 44 +++++++++++++++++++++----------------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index e55ec90d7..c75d1a36d 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -307,7 +307,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) call znot_t_v6(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m) else if (sfc_z0_type == 7) then call znot_t_v7(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m) - else if (sfc_z0_type > 0) then + else if (sfc_z0_type /= 0) then write(0,*)'no option for sfc_z0_type=',sfc_z0_type stop endif @@ -322,35 +322,33 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! ! update z0 over ocean ! - if (sfc_z0_type >= 0) then - if (sfc_z0_type == 0) then - z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) + if (sfc_z0_type == 0) then + z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) ! mbek -- toga-coare flux algorithm -! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i) +! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i) ! new implementation of z0 -! cc = ustar(i) * z0 / rnu -! pp = cc / (1. + cc) -! ff = grav * arnu / (charnock * ustar(i) ** 3) -! z0 = arnu / (ustar(i) * ff ** pp) - - if (redrag) then - z0rl_wat(i) = 100.0 * max(min(z0, z0s_max), 1.e-7) - else - z0rl_wat(i) = 100.0 * max(min(z0,.1), 1.e-7) - endif - - elseif (sfc_z0_type == 6) then ! wang - call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m - z0rl_wat(i) = 100.0 * z0 ! cm - elseif (sfc_z0_type == 7) then ! wang - call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m - z0rl_wat(i) = 100.0 * z0 ! cm +! cc = ustar(i) * z0 / rnu +! pp = cc / (1. + cc) +! ff = grav * arnu / (charnock * ustar(i) ** 3) +! z0 = arnu / (ustar(i) * ff ** pp) + + if (redrag) then + z0rl_wat(i) = 100.0 * max(min(z0, z0s_max), 1.e-7) else - z0rl_wat(i) = 1.0e-4 + z0rl_wat(i) = 100.0 * max(min(z0,.1), 1.e-7) endif + elseif (sfc_z0_type == 6) then ! wang + call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m + z0rl_wat(i) = 100.0 * z0 ! cm + elseif (sfc_z0_type == 7) then ! wang + call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m + z0rl_wat(i) = 100.0 * z0 ! cm + else + z0rl_wat(i) = 1.0e-4 endif + endif ! end of if(open ocean) ! endif ! end of if(flagiter) loop From 810426e0df91a5883a1471ce099e941dd3f0e4f8 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Thu, 11 Jun 2020 13:27:04 +0000 Subject: [PATCH 3/4] remove whitespaces in sfc_diff.f --- physics/sfc_diff.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index c75d1a36d..7af159a1d 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -224,7 +224,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) tem1 = 1.0 - sigmaf(i) ztmax = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) + & * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) ! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land @@ -265,7 +265,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) tem1 = 1.0 - sigmaf(i) ztmax = z0max*exp( - tem1*tem1 - & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) + & * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05))) ztmax = max(ztmax, 1.0e-6) ! call stability From f84468b1a3a732631aeee01a741335d1ed49dff2 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Thu, 11 Jun 2020 15:47:43 +0000 Subject: [PATCH 4/4] keep z0 unchanged in coupled mode --- physics/sfc_diff.f | 45 +++++++++++++++++++++++---------------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/physics/sfc_diff.f b/physics/sfc_diff.f index 7af159a1d..fd35d5964 100644 --- a/physics/sfc_diff.f +++ b/physics/sfc_diff.f @@ -307,7 +307,7 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) call znot_t_v6(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m) else if (sfc_z0_type == 7) then call znot_t_v7(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m) - else if (sfc_z0_type /= 0) then + else if (sfc_z0_type > 0) then write(0,*)'no option for sfc_z0_type=',sfc_z0_type stop endif @@ -322,33 +322,34 @@ subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav, & !intent(in) ! ! update z0 over ocean ! - if (sfc_z0_type == 0) then - z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) + if (sfc_z0_type >= 0) then + if (sfc_z0_type == 0) then + z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i) ! mbek -- toga-coare flux algorithm -! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i) +! z0 = (charnock / grav) * ustar(i)*ustar(i) + arnu/ustar(i) ! new implementation of z0 -! cc = ustar(i) * z0 / rnu -! pp = cc / (1. + cc) -! ff = grav * arnu / (charnock * ustar(i) ** 3) -! z0 = arnu / (ustar(i) * ff ** pp) - - if (redrag) then - z0rl_wat(i) = 100.0 * max(min(z0, z0s_max), 1.e-7) +! cc = ustar(i) * z0 / rnu +! pp = cc / (1. + cc) +! ff = grav * arnu / (charnock * ustar(i) ** 3) +! z0 = arnu / (ustar(i) * ff ** pp) + + if (redrag) then + z0rl_wat(i) = 100.0 * max(min(z0, z0s_max), 1.e-7) + else + z0rl_wat(i) = 100.0 * max(min(z0,.1), 1.e-7) + endif + + elseif (sfc_z0_type == 6) then ! wang + call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m + z0rl_wat(i) = 100.0 * z0 ! cm + elseif (sfc_z0_type == 7) then ! wang + call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m + z0rl_wat(i) = 100.0 * z0 ! cm else - z0rl_wat(i) = 100.0 * max(min(z0,.1), 1.e-7) + z0rl_wat(i) = 1.0e-4 endif - - elseif (sfc_z0_type == 6) then ! wang - call znot_m_v6(wind10m, z0) ! wind, m/s, z0, m - z0rl_wat(i) = 100.0 * z0 ! cm - elseif (sfc_z0_type == 7) then ! wang - call znot_m_v7(wind10m, z0) ! wind, m/s, z0, m - z0rl_wat(i) = 100.0 * z0 ! cm - else - z0rl_wat(i) = 1.0e-4 endif - endif ! end of if(open ocean) ! endif ! end of if(flagiter) loop