diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 90ead827..b91a22c0 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1202,7 +1202,9 @@ subroutine GFS_physics_driver & Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & - tsurf, flag_iter) !, Model%redrag, Model%z0s_max, & + tsurf, flag_iter, & + Model%alpha_stable, Model%alpha_unstable, Model%tune_ocean_surface_layer) +! , Model%redrag, Model%z0s_max, & !Model%do_z0_moon, Model%do_z0_hwrf15, & !Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & !Model%wind_th_hwrf) @@ -1225,7 +1227,8 @@ subroutine GFS_physics_driver & tsurf, flag_iter, Model%redrag, Model%z0s_max, & Model%do_z0_moon, Model%do_z0_hwrf15, & Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & - Model%wind_th_hwrf) + Model%wind_th_hwrf, & + Model%alpha_stable, Model%alpha_unstable, Model%tune_ocean_surface_layer) else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& @@ -1709,7 +1712,8 @@ subroutine GFS_physics_driver & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl,& gamt, gamq, dkt, kinver, Model%xkzm_m, Model%xkzm_h, & Model%xkzm_s, lprnt, ipr, & - Model%xkzminv, Model%moninq_fac) + Model%xkzminv, Model%moninq_fac, Model%alpha_stable, & + Model%alpha_unstable, Model%tune_ocean_surface_layer) ! if (lprnt) write(0,*)' dtdtm=',(dtdt(ipr,k),k=1,15) ! if (lprnt) write(0,*)' dqdtm=',(dqdt(ipr,k,1),k=1,15) @@ -1731,7 +1735,9 @@ subroutine GFS_physics_driver & Model%xkzm_ml, Model%xkzm_hl, Model%xkzm_mi, Model%xkzm_hi, & Model%xkzm_s, Model%xkzminv, Model%do_dk_hb19, & Model%xkzm_lim, Model%xkgdx, & - Model%rlmn, Model%rlmx, Model%cap_k0_land, dkt) + Model%rlmn, Model%rlmx, Model%cap_k0_land, dkt, & + Model%alpha_stable, Model%alpha_unstable, & + Model%tune_ocean_surface_layer) elseif (Model%isatmedmf == 1) then do i=1,im @@ -1758,7 +1764,8 @@ subroutine GFS_physics_driver & Model%xkzm_s, Model%xkzminv, Model%rlmx, Model%zolcru, & Model%cs0, Model%do_dk_hb19, Model%xkgdx, & Model%dspfac, Model%bl_upfr, Model%bl_dnfr, dkt, & - flux_cg, flux_en) !cg as up and en as down + flux_cg, flux_en, Model%alpha_stable, Model%alpha_unstable, & + Model%tune_ocean_surface_layer) !cg as up and en as down endif elseif (Model%ysupbl) then diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index c75ed8b5..561ab78c 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -680,6 +680,9 @@ module GFS_typedefs logical :: do_z0_hwrf17 !< flag for using z0 scheme in 2017 HWRF (kgao) logical :: do_z0_hwrf17_hwonly !< flag for using z0 scheme in 2017 HWRF only under high wind (kgao) real(kind=kind_phys) :: wind_th_hwrf !< wind speed threshold when z0 level off as in HWRF (kgao) + real(kind=kind_phys) :: alpha_stable !< tuning parameter for the dimensionless momentum and scalar gradient function in the surface layer + real(kind=kind_phys) :: alpha_unstable !< tuning parameter for the dimensionless momentum gradient function in the surface layer + logical :: tune_ocean_surface_layer !< if true, use alphas above to modify the surface layer dimensionless gradients over oceans logical :: hybedmf !< flag for hybrid edmf pbl scheme logical :: myj_pbl !< flag for NAM MYJ tke scheme logical :: ysupbl !< flag for ysu pbl scheme (version in WRFV3.8) @@ -2401,6 +2404,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: do_z0_hwrf17 = .false. !< flag for using z0 scheme in 2017 HWRF logical :: do_z0_hwrf17_hwonly = .false. !< flag for using z0 scheme in 2017 HWRF only under high wind real(kind=kind_phys) :: wind_th_hwrf = 33. !< wind speed threshold when z0 level off as in HWRF + real(kind=kind_phys) :: alpha_stable = 5. !< tuning parameter for the dimensionless momentum and scalar gradient function in the surface layer + real(kind=kind_phys) :: alpha_unstable = 16. !< tuning parameter for the dimensionless momentum gradient function in the surface layer + logical :: tune_ocean_surface_layer = .false. !< if true, use alphas above to modify the surface layer dimensionless gradients over oceans logical :: hybedmf = .false. !< flag for hybrid edmf pbl scheme logical :: myj_pbl = .false. !< flag for NAM MYJ tke-based scheme logical :: ysupbl = .false. !< flag for hybrid edmf pbl scheme @@ -2633,7 +2639,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & h2o_phys, pdfcld, shcnvcw, redrag, sfc_gfdl, z0s_max, & sfc_coupled, do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & - do_z0_hwrf17_hwonly, wind_th_hwrf, & + do_z0_hwrf17_hwonly, wind_th_hwrf, alpha_stable, & + alpha_unstable, tune_ocean_surface_layer, & hybedmf, dspheat, lheatstrg, hour_canopy, afac_canopy, & cnvcld, no_pbl, xkzm_lim, xkzm_fac, xkgdx, & rlmn, rlmx, zolcru, cs0, & @@ -2878,6 +2885,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%do_z0_hwrf17 = do_z0_hwrf17 Model%do_z0_hwrf17_hwonly = do_z0_hwrf17_hwonly Model%wind_th_hwrf = wind_th_hwrf + Model%alpha_stable = alpha_stable + Model%alpha_unstable = alpha_unstable + Model%tune_ocean_surface_layer = tune_ocean_surface_layer Model%hybedmf = hybedmf Model%myj_pbl = myj_pbl Model%ysupbl = ysupbl @@ -3585,6 +3595,9 @@ subroutine control_print(Model) print *, ' do_z0_hwrf17 : ', Model%do_z0_hwrf17 print *, ' do_z0_hwrf17_hwonly : ', Model%do_z0_hwrf17_hwonly print *, ' wind_th_hwrf : ', Model%wind_th_hwrf + print *, ' alpha_stable : ', Model%alpha_stable + print *, ' alpha_unstable : ', Model%alpha_unstable + print *, ' tune_ocean_surface_layer : ', Model%tune_ocean_surface_layer print *, ' hybedmf : ', Model%hybedmf print *, ' myj_pbl : ', Model%myj_pbl print *, ' ysupbl : ', Model%ysupbl diff --git a/gsmphys/moninedmf.f b/gsmphys/moninedmf.f index 1c42c6d7..28e010b5 100755 --- a/gsmphys/moninedmf.f +++ b/gsmphys/moninedmf.f @@ -93,7 +93,8 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & & prsi,del,prsl,prslk,phii,phil,delt,dspheat, & & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt, & & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr, & - & xkzminv,moninq_fac) + & xkzminv,moninq_fac, alpha_stable, alpha_unstable, & + & tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -103,11 +104,12 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & ! ! arguments ! - logical lprnt + logical lprnt, tune_ocean_surface_layer integer ipr integer ix, im, km, ntrac, ntcw, kpbl(im), kinver(im) ! - real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s + real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s, alpha_stable, + & alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & & tau(im,km), rtg(im,km,ntrac), & & u1(ix,km), v1(ix,km), & @@ -512,24 +514,12 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & !! \f[ !! w_s = (u_*^3 + 7\epsilon k w_*^3)^{1/3} !! \f] + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) + do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then -! phim(i) = (1.-aphi16*zol1)**(-1./4.) -! phih(i) = (1.-aphi16*zol1)**(-1./2.) - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif wscale(i) = ustar(i)/phim(i) ustmin(i) = ustar(i)/aphi5 wscale(i) = max(wscale(i),ustmin(i)) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index 4c7e32f4..dbb31c8d 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -36,7 +36,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, & kinver,xkzm_mo,xkzm_ho,xkzm_ml,xkzm_hl,xkzm_mi,xkzm_hi, & xkzm_s,xkzinv,do_dk_hb19,xkzm_lim,xkgdx, - & rlmn, rlmx, cap_k0_land, dkt_out) + & rlmn, rlmx, cap_k0_land, dkt_out, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -53,7 +54,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! real(kind=kind_phys) delt, xkzm_s, xkzm_lim, & xkzm_mo, xkzm_ho, xkzm_ml, xkzm_hl, - & xkzm_mi, xkzm_hi + & xkzm_mi, xkzm_hi, + & alpha_stable, alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & tdt(im,km), rtg(im,km,ntrac), & u1(ix,km), v1(ix,km), @@ -77,7 +79,7 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! kgao note - q1 and rtg are local var now ! - logical dspheat, cap_k0_land, do_dk_hb19 + logical dspheat, cap_k0_land, do_dk_hb19, tune_ocean_surface_layer ! flag for tke dissipative heating real(kind=kind_phys),dimension(1:im,1:km),intent(OUT)::dkt_out @@ -592,24 +594,10 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! ! compute similarity parameters ! - do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif -! - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif - enddo + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) ! do i=1,im if(pblflg(i)) then @@ -1564,3 +1552,61 @@ subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) c----------------------------------------------------------------------- return end + + subroutine get_similarity_parameters( + & im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl, + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol + & ) + + use machine, only : kind_phys + + implicit none + ! input + integer, intent(in) :: im + real(kind=kind_phys), dimension(1:im), intent(in) :: + & rbsoil, fm, fh, hpbl, zl + real(kind=kind_phys), intent(in) :: rimin, zfmin, aphi16, aphi5, + & alpha_unstable, alpha_stable, sfcfrac + logical, dimension(1:im), intent(in) :: sfcflg + logical, intent(in) :: tune_ocean_surface_layer + + ! output + real(kind=kind_phys), dimension(1:im), intent(out) :: phim, phih, + & zol + + ! local + integer i + real(kind=kind_phys) :: zol1, tem + + do i=1,im + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) + if(sfcflg(i)) then + zol(i) = min(zol(i),-zfmin) + else + zol(i) = max(zol(i),zfmin) + endif + ! + zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i) + if(sfcflg(i)) then + tem = 1.0 / (1. - aphi16*zol1) + phih(i) = sqrt(tem) + if (tune_ocean_surface_layer) then + tem = 1.0 / (1. - alpha_unstable*zol1) + phim(i) = sqrt(sqrt(tem)) + else + phim(i) = sqrt(phih(i)) + end if + else + if (tune_ocean_surface_layer) then + phim(i) = 1. + alpha_stable*zol1 + else + phim(i) = 1. + aphi5*zol1 + end if + phih(i) = phim(i) + endif + enddo + return + end subroutine + diff --git a/gsmphys/satmedmfvdifq.f b/gsmphys/satmedmfvdifq.f index 4b233c04..4f99b6f2 100644 --- a/gsmphys/satmedmfvdifq.f +++ b/gsmphys/satmedmfvdifq.f @@ -56,7 +56,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, & kinver,xkzm_mo,xkzm_ho,xkzm_ml,xkzm_hl,xkzm_mi,xkzm_hi, & xkzm_s,xkzinv,rlmx,zolcru,cs0, & do_dk_hb19,xkgdx,dspfac,bl_upfr,bl_dnfr,dkt_out, - & flux_up, flux_dn) + & flux_up, flux_dn, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -73,7 +74,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, ! real(kind=kind_phys) delt, xkzm_mo, xkzm_ho, xkzm_s, dspfac, & bl_upfr, bl_dnfr, xkzm_ml, xkzm_hl, - & xkzm_mi, xkzm_hi + & xkzm_mi, xkzm_hi, alpha_stable, + & alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & tdt(im,km), rtg(im,km,ntrac), & u1(ix,km), v1(ix,km), @@ -97,7 +99,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, & rtg_in(im,km,ntrac) ! kgao note - q1 and rtg are local var now ! - logical dspheat, do_dk_hb19 + logical dspheat, do_dk_hb19, tune_ocean_surface_layer ! flag for tke dissipative heating real(kind=kind_phys)::dkt_out(im,km),flux_up(im,km),flux_dn(im,km) ! @@ -619,24 +621,10 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, ! ! compute similarity parameters ! - do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif -! - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif - enddo + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) ! do i=1,im if(pblflg(i)) then diff --git a/gsmphys/sfc_diff_coupled.f b/gsmphys/sfc_diff_coupled.f index dfe08bcb..e8eb2d17 100644 --- a/gsmphys/sfc_diff_coupled.f +++ b/gsmphys/sfc_diff_coupled.f @@ -4,7 +4,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, & stress,fm,fh, & ustar,wind,ddvel,fm10,fh2, & sigmaf,vegtype,shdmax,ivegsrc, - & tsurf,flag_iter) !,redrag, + & tsurf,flag_iter, alpha_stable, alpha_unstable, + & tune_ocean_surface_layer) +! ,redrag, ! & z0s_max) ! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, ! & do_z0_hwrf17_hwonly, wind_th_hwrf) @@ -26,10 +28,12 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, &, prsl1, prslki, stress &, fm, fh, ustar, wind, ddvel &, fm10, fh2, sigmaf, shdmax - &, tsurf, snwdph + &, tsurf, snwdph, alpha_stable + &, alpha_unstable + integer, dimension(im) ::vegtype, islimsk - logical flag_iter(im) + logical flag_iter(im), tune_ocean_surface_layer ! logical redrag ! logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao ! &, do_z0_hwrf17_hwonly ! kgao @@ -143,7 +147,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), stress(i), ustar(i)) @@ -165,7 +171,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, ! --- call similarity ! kgao: use z0/zt to do sfc diagnosis call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 87aeaf7f..373d8f1f 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -22,7 +22,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & tsurf,flag_iter,redrag, & z0s_max, & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, - & do_z0_hwrf17_hwonly, wind_th_hwrf) + & do_z0_hwrf17_hwonly, wind_th_hwrf, + & alpha_stable, alpha_unstable, + & tune_ocean_surface_layer) ! oct 2019 - a clean and updated version by Kun Gao at GFDL (Kun.Gao@noaa.gov) @@ -50,10 +52,10 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, &, fm10, fh2, sigmaf, shdmax &, tsurf, snwdph &, fm_neutral, fm10_neutral ! Sofar added Spring 2023 - real(kind=kind_phys) :: ws1, ws10n ! Sofar added 10/19/23 + real(kind=kind_phys) :: ws1, ws10n, alpha_stable, alpha_unstable ! Sofar added 10/19/23 integer, dimension(im) ::vegtype, islimsk - logical flag_iter(im) + logical flag_iter(im), tune_ocean_surface_layer logical redrag logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao &, do_z0_hwrf17_hwonly ! kgao @@ -171,7 +173,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -202,7 +206,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -302,7 +308,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -508,17 +516,23 @@ end subroutine cal_z0_moon ! ======================================================================= subroutine monin_obukhov_similarity - & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, - & rb, fm, fh, fm10, fh2, + & ( ilsimask, z1, snwdph, thv1, wind, z0max, ztmax, tvs, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, + & rb, fm, fh, fm10, fh2, & fm_neutral, fm10_neutral, !(ADDED by Sofar) & cm, ch, stress, ustar) ! --- input +! ilsimask - land/sea/ice mask ! z1 - lowest model level height ! snwdph - surface snow thickness ! wind - wind speed at lowest model layer ! thv1 - virtual potential temp at lowest model layer ! tvs - surface temp +! alpha_stable - tuning parameter (inverse crit. Richarson number) in dimensinless shear +! and scalar gradient; stable case +! alpha_unstable - tuning parameter in dimensionless wind shear; unstable case +! tune_ocean_surface_layer - if true, use alphas to tune surface layer over oceans ! z0max - surface roughness length for momentum ! ztmax - surface roughness length for heat ! @@ -536,7 +550,11 @@ subroutine monin_obukhov_similarity ! --- inputs: real(kind=kind_phys), intent(in) :: - & z1, snwdph, thv1, wind, z0max, ztmax, tvs + & z1, snwdph, thv1, wind, z0max, ztmax, tvs, alpha_stable, + & alpha_unstable + + integer, intent(in) :: ilsimask + logical, intent(in) :: tune_ocean_surface_layer ! --- outputs: real(kind=kind_phys), intent(out) :: @@ -547,7 +565,7 @@ subroutine monin_obukhov_similarity ! --- locals: real(kind=kind_phys), parameter :: alpha=5., a0=-3.975 - &, a1=12.32, alpha4=4.0*alpha + &, a1=12.32, alpha4_nonocean=4.0*alpha &, b1=-7.755, b2=6.041, alpha2=alpha+alpha, beta=1.0 &, a0p=-7.941, a1p=24.75, b1p=-8.705, b2p=7.899 &, ztmin1=-999.0, ca=.4 @@ -557,10 +575,16 @@ subroutine monin_obukhov_similarity & z1i, & fms, fhs, hl0, hl0inf, hlinf, & hl110, hlt, hltinf, olinf, - & tem1, tem2, ztmax1 + & tem1, tem2, ztmax1, alpha4 z1i = 1.0 / z1 + if (ilsimask == 0 .and. tune_ocean_surface_layer) then + alpha4 = 4. * alpha_stable + else + alpha4 = alpha4_nonocean + end if + tem1 = z0max/z1 if (abs(1.0-tem1) > 1.0e-6) then ztmax1 = - beta*log(tem1)/(alpha2*(1.-tem1)) @@ -642,32 +666,49 @@ subroutine monin_obukhov_similarity ! ! get pm and ph ! + if (hlinf >= -0.5) then hl1 = hlinf - pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) - ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) hl110 = hl1 * 10. * z1i hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = (a0 + a1*hl110) * hl110 / (1.+(b1+b2*hl110)*hl110) + if (ilsimask == 0 .and. tune_ocean_surface_layer) then + call get_psim_unstable(hlinf, z1i, z0max, + & alpha_unstable, pm) + call get_psim_unstable(hl110, 1./10., z0max, + & alpha_unstable, pm10) + else + pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) + pm10 = (a0 + a1*hl110) * hl110 / + & (1.+(b1+b2*hl110)*hl110) + endif + ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) hl12 = (hl1+hl1) * z1i hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = (a0p + a1p*hl12) * hl12 / (1.+(b1p+b2p*hl12)*hl12) else ! hlinf < 0.05 hl1 = -hlinf + hl110 = hl1 * 10. * z1i + hl110 = min(max(hl110, ztmin1), ztmax1) tem1 = 1.0 / sqrt(hl1) - pm = log(hl1) + 2. * sqrt(tem1) - .8776 + if (ilsimask == 0 .and. tune_ocean_surface_layer) then + call get_psim_unstable(hlinf, z1i, z0max, + & alpha_unstable, pm) + call get_psim_unstable(-hl110, 1./10., z0max, + & alpha_unstable, pm10) + else + pm = log(hl1) + 2. * sqrt(tem1) - .8776 + pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 + endif ph = log(hl1) + .5 * tem1 + 1.386 ! pm = log(hl1) + 2.0 * hl1 ** (-.25) - .8776 ! ph = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386 - hl110 = hl1 * 10. * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 ! pm10 = log(hl110) + 2. * hl110 ** (-.25) - .8776 hl12 = (hl1+hl1) * z1i hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = log(hl12) + 0.5 / sqrt(hl12) + 1.386 ! ph2 = log(hl12) + .5 * hl12 ** (-.5) + 1.386 endif + endif ! end of if (dtv >= 0 ) then loop ! @@ -687,3 +728,18 @@ subroutine monin_obukhov_similarity return end subroutine monin_obukhov_similarity + + subroutine get_psim_unstable(zeta, z1i, z0, alpha, psi_m) + use machine , only : kind_phys + implicit none + real(kind=kind_phys),intent(in) :: zeta, z1i, z0, alpha + real(kind=kind_phys),intent(out) :: psi_m + real(kind=kind_phys) :: tem1, zeta0, Rz, R0 + tem1 = zeta * z1i + zeta0 = z0 * tem1 + Rz = (1. - alpha * zeta)**0.25 + R0 = (1. - alpha * zeta0)**0.25 + psi_m = log((1. + Rz)**2. * (1. + Rz**2.) / ((1. + R0)**2. + & *(1. + R0**2.))) + 2. * atan(R0) - 2. * atan(Rz) + return + end subroutine