32 SUBROUTINE gwdps_v0(IM, km, imx, do_tofd,
33 & Pdvdt, Pdudt, Pdtdt, Pkdis, U1,V1,T1,Q1,KPBL,
34 & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DTP,KDT,
35 & sgh30, HPRIME,OC,OA4,CLX4,THETA,vSIGMA,vGAMMA,ELVMAXD,
36 & DUSFC, DVSFC, xlatd, sinlat, coslat, sparea,
37 & cdmbgwd, me, master, rdxzb, con_g, con_omega,
38 & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
39 & dudt_mtb, dudt_ogw, dudt_tms)
48 USE machine ,
ONLY : kind_phys
50 &, pi, rad_to_deg, deg_to_rad, pi2
51 &, rdi, gor, grcp, gocp, fv, gr2
52 &, bnv2min, dw2min, velmin, arad
55 &, hpmax, hpmin, sigfaci => sigfac
56 &, dpmin, minwnd, hminmt, hncrit
57 &, rlolev, gmax, veleps, factop
58 &, frc, ce, ceofrc, frmax, cg
60 &, cdmb, cleff, fcrit_gfs, fcrit_mtb
61 &, n_tofd, ze_tofd, ztop_tofd
67 integer,
parameter :: kp = kind_phys
68 character(len=8) :: strsolver=
'PSS-1986'
69 integer,
intent(in) :: im, km, imx, kdt
70 integer,
intent(in) :: me, master
71 logical,
intent(in) :: do_tofd
72 real(kind=kind_phys),
parameter :: sigfac = 3, sigfacs = 0.5
73 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
74 integer,
intent(in) :: KPBL(IM)
75 real(kind=kind_phys),
intent(in) :: dtp
76 real(kind=kind_phys),
intent(in) :: cdmbgwd(2)
78 real(kind=kind_phys),
intent(in),
dimension(im,km) ::
80 & del, prsl, prslk, phil
81 real(kind=kind_phys),
intent(in),
dimension(im,km+1):: prsi, phii
82 real(kind=kind_phys),
intent(in) :: xlatd(im),sinlat(im),
84 real(kind=kind_phys),
intent(in) :: sparea(im)
86 real(kind=kind_phys),
intent(in) :: oc(im), oa4(im,4), clx4(im,4)
87 real(kind=kind_phys),
intent(in) :: hprime(im), sgh30(im)
88 real(kind=kind_phys),
intent(in) :: elvmaxd(im), theta(im)
89 real(kind=kind_phys),
intent(in) :: vsigma(im), vgamma(im)
90 real(kind=kind_phys) :: sigma(im), gamma(im)
92 real(kind=kind_phys),
intent(in) :: con_g, con_omega
95 real(kind=kind_phys),
dimension(im,km),
intent(out) ::
96 & pdvdt, pdudt, pkdis, pdtdt
98 &, dudt_mtb, dudt_ogw, dudt_tms
100 real(kind=kind_phys),
dimension(im) :: rdxzb, zmtb, zogw
101 &, tau_ogw, tau_mtb, tau_tofd
115 real(kind=kind_phys) :: gammin = 0.00999999_kp
116 real(kind=kind_phys),
parameter :: nhilmax = 25.0_kp
117 real(kind=kind_phys),
parameter :: sso_min = 3000.0_kp
118 logical,
parameter :: do_adjoro = .true.
120 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
121 real(kind=kind_phys) :: belpmin, dsmin, dsmax
123 real(kind=kind_phys) :: xlingfs
128 real(kind=kind_phys),
dimension(im,km) :: ri_n, bnv2, ro
131 real(kind=kind_phys),
dimension(im) :: oa, clx , elvmax, wk
134 real(kind=kind_phys),
dimension(im,km) :: db, ang, uds
136 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
137 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
138 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
144 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
145 real(kind=kind_phys),
dimension(km) :: utofd1, vtofd1
146 &, epstofd1, krf_tofd1
148 real(kind=kind_phys),
dimension(im, km) :: axtms, aytms
154 real(kind=kind_phys),
dimension(im) :: xn, yn, ubar, vbar, ulow,
155 & roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
157 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
158 real(kind=kind_phys) :: taub(im), taulin(im), heff, hsat, hdis
160 integer,
dimension(im) :: kref, idxzb, ipt, kreflm,
165 real(kind=kind_phys) :: bnv, fr, ri_gw
166 &, brvf, tem, tem1, tem2, temc, temv
167 &, ti, rdz, dw2, shr2, bvf2
168 &, rdelks, efact, coefm, gfobnv
169 &, scork, rscor, hd, fro, sira
170 &, dtaux, dtauy, pkp1log, pklog
171 &, grav2, rcpdt, windik, wdir
172 &, sigmin, dxres,sigres,hdxres
174 &, kxridge, inv_b2eff, zw1, zw2
175 &, belps, aelps, nhills, selps
177 integer :: kmm1, kmm2, lcap, lcapp1
178 &, npt, kbps, kbpsp1,kbpsm1
179 &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
180 &, k_mtb, k_zlow, ktrial, klevm1, i, j, k
182 rcpdt = 1.0 / (cpd*dtp)
187 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
188 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
190 dxres = pi2*arad/float(imx)
191 hdxres = 0.5_kp*dxres
195 gammin = min(sso_min/dxres, 1.)
198 sigmin = 2.*hpmin/dxres
202 kxridge = float(imx)/arad * cdmbgwd(2)
216 sigma(i) = max(vsigma(i), sigmin)
217 gamma(i) = max(vgamma(i), gammin)
236 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin )
then
242 sigres = max(sigmin, sigma(i))
244 dxres = sqrt(sparea(i))
245 if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
246 aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
247 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i),.5*dxres)
251 if( aelps < sso_min .and. do_adjoro)
then
256 if (belps < sso_min )
then
258 belps = aelps*gamma(i)
260 gamma(i) = min(aelps/belps, 1.0)
262 sigma(i) = 2.*hprime(i)/aelps
263 gamma(i) = min(aelps/belps, 1.0)
266 selps = belps*belps*gamma(i)*4.
267 nhills = min(nhilmax, sparea(i)/selps)
297 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
298 lcap = km ; lcapp1 = lcap + 1
302 elvmax(j) = min(elvmaxd(j)*0. + sigfac * hprime(j), hncrit)
309 ztoph = sigfac * hprime(j)
310 zlowh = sigfacs* hprime(j)
311 pkp1log = phil(j,k+1) * rgrav
312 pklog = phil(j,k) * rgrav
315 if (( ztoph <= pkp1log) .and. (ztoph >= pklog) )
316 & iwklm(i) = max(iwklm(i), k+1 )
318 if (zlowh <= pkp1log .and. zlowh >= pklog)
319 & izlow(i) = max(izlow(i),k)
326 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
327 vtk(i,k) = vtj(i,k) / prslk(j,k)
328 ro(i,k) = rdi * prsl(j,k) / vtj(i,k)
338 rdz = grav / (phil(j,k+1) - phil(j,k))
339 tem1 = u1(j,k) - u1(j,k+1)
340 tem2 = v1(j,k) - v1(j,k+1)
341 dw2 = tem1*tem1 + tem2*tem2
342 shr2 = max(dw2,dw2min) * rdz * rdz
347 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))
348 & / (vtk(i,k+1)+vtk(i,k))
349 bnv2(i,k+1) = max( bvf2, bnv2min )
350 ri_n(i,k+1) = bnv2(i,k)/shr2
358 bnv2(i,k) = bnv2(i,k+1)
366 if (k_zlow == iwklm(i)) k_zlow = 1
367 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,iwklm(i)))
379 if (k_zlow == iwklm(i)) k_zlow = 1
380 DO k = k_zlow, iwklm(i)-1
382 rdelks = del(j,k) * delks(i)
383 ubar(i) = ubar(i) + rdelks * u1(j,k)
384 vbar(i) = vbar(i) + rdelks * v1(j,k)
385 roll(i) = roll(i) + rdelks * ro(i,k)
387 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
398 DO k = iwklm(i), 1, -1
399 phiang = atan2(v1(j,k),u1(j,k))*rad_to_deg
400 ang(i,k) = ( theta(j) - phiang )
401 if ( ang(i,k) > 90. ) ang(i,k) = ang(i,k) - 180.
402 if ( ang(i,k) < -90. ) ang(i,k) = ang(i,k) + 180.
403 ang(i,k) = ang(i,k) * deg_to_rad
405 & max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
407 IF (idxzb(i) == 0 )
then
408 dz_blk = ( phii(j,k+1) - phii(j,k) ) *rgrav
409 pe(i) = pe(i) + bnv2(i,k) *
410 & ( elvmax(j) - phil(j,k)*rgrav ) * dz_blk
412 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
413 ek(i) = 0.5 * up(i) * up(i)
415 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
419 IF ( ph_blk >= fcrit_gfs )
THEN
421 zmtb(j) = phil(j, k)*rgrav
422 rdxzb(j) = real(k, kind=kind_phys)
433 bnv = sqrt( bnv2bar(i) )
434 heff = 2.*min(hprime(j),hpmax)
435 zw2 = ubar(i)*ubar(i)+vbar(i)*vbar(i)
436 ulow(i) = sqrt(max(zw2,dw2min))
437 fr = heff*bnv/ulow(i)
438 zw1 = max(heff*(1. -fcrit_gfs/fr), 0.0)
439 zw2 = phil(j,2)*rgrav
440 if (fr > fcrit_gfs .and. zw1 > zw2 )
then
442 pkp1log = phil(j,k+1) * rgrav
443 pklog = phil(j,k) * rgrav
444 if (zw1 <= pkp1log .and. zw1 >= pklog)
exit
447 zmtb(j) = phil(j, k)*rgrav
462 IF ( idxzb(i) > 0 )
then
464 gam2 = gamma(j)*gamma(j)
465 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
466 cgam = 0.48*gamma(j) + 0.30*gam2
467 DO k = idxzb(i)-1, 1, -1
469 zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) /
470 & ( phil(j,k ) + grav * hprime(j) ) )
474 sinang2 = 1.0 - cosang2
479 rdem = cosang2 + gam2 * sinang2
480 rnom = cosang2*gam2 + sinang2
487 rdem = max(rdem, 1.e-6)
489 zr = max( 2. - r, 0. )
491 sigres = max(sigmin, sigma(j))
492 if (hprime(j)/sigres > dxres) sigres = hprime(j)/dxres
493 mtbridge = zr * sigres*zlen / hprime(j)
498 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam* sinang2)
499 db(i,k)= dbtmp * uds(i,k)
525 tem = (prsi(j,1) - prsi(j,k))
526 if (tem < dpmin) iwk(i) = k
547 k_mtb = max(1, idxzb(i))
549 kref(i) = max(iwk(i), kpbl(j)+1 )
550 kref(i) = max(kref(i), iwklm(i) )
552 if (kref(i) <= idxzb(i)) kref(i) = idxzb(i) + 1
553 kbps = max(kbps, kref(i))
554 kmps = min(kmps, kref(i))
556 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
568 k_mtb = max(1, idxzb(i))
570 IF (k < kref(i))
THEN
572 rdelks = del(j,k) * delks(i)
573 ubar(i) = ubar(i) + rdelks * u1(j,k)
574 vbar(i) = vbar(i) + rdelks * v1(j,k)
575 roll(i) = roll(i) + rdelks * ro(i,k)
576 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
584 wdir = atan2(ubar(i),vbar(i)) + pi
585 idir = mod(nint(fdir*wdir),mdir) + 1
587 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
588 clx(i) = clx4(j,mod(nwd-1,4)+1)
594 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
595 xn(i) = ubar(i) / ulow(i)
596 yn(i) = vbar(i) / ulow(i)
602 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)
603 & + (v1(j,k)+v1(j,k+1))*yn(i))
618 taub(:) = 0. ; taulin(:)= 0.
621 bnv = sqrt( bnv2bar(i) )
622 heff = min(hprime(j),hpmax)
624 if( zmtb(j) > 0.) heff = max(sigfac*heff-zmtb(j), 0.)/sigfac
627 hsat = fcrit_gfs*ulow(i)/bnv
628 heff = min(heff, hsat)
630 fr = min(bnv * heff /ulow(i), frmax)
632 efact = (oa(i) + 2.) ** (ceofrc*fr)
633 efact = min( max(efact,efmin), efmax )
635 coefm = (1. + clx(i)) ** (oa(i)+1.)
637 xlinv(i) = coefm * cleff
638 xlingfs = coefm * cleff
640 tem = fr * fr * oc(j)
641 gfobnv = gmax * tem / ((tem + cg)*bnv)
645 sigres = max(sigmin, sigma(j))
646 if (heff/sigres > hdxres) sigres = heff/hdxres
647 inv_b2eff = 0.5*sigres/heff
648 kxridge = 1.0 / sqrt(sparea(j))
650 taulin(i) = 0.5*roll(i)*xlinv(i)*bnv*ulow(i)*
653 if ( fr > fcrit_gfs )
then
654 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
655 & * ulow(i) * gfobnv * efact
659 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
660 & * ulow(i) * gfobnv * efact
667 k = max(1, kref(i)-1)
668 tem = max(velco(i,k)*velco(i,k), dw2min)
669 scor(i) = bnv2(i,k) / tem
673 zogw(j) = phii(j, kref(i)) *rgrav
680 IF (k <= kref(i)) taup(i,k) = taub(i)
684 if (strsolver ==
'PSS-1986')
then
698 IF (k >= kref(i))
THEN
699 icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) < ric)
700 & .OR. (velco(i,k) <= 0.0)
705 IF (k >= kref(i))
THEN
706 IF (.NOT.icrilv(i) .AND. taup(i,k) > 0.0 )
THEN
707 temv = 1.0 / max(velco(i,k), velmin)
709 IF (oa(i) > 0. .AND. kp1 < kref(i))
THEN
710 scork = bnv2(i,k) * temv * temv
711 rscor = min(1.0, scork / scor(i))
717 brvf = sqrt(bnv2(i,k))
720 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
721 & * max(velco(i,k), velmin)
722 hd = sqrt(taup(i,k) / tem1)
723 fro = brvf * hd * temv
728 tem2 = sqrt(ri_n(i,k))
729 tem = 1. + tem2 * fro
730 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
735 IF (ri_gw <= ric .AND.
736 & (oa(i) <= 0. .OR. kp1 >= kref(i) ))
THEN
737 temc = 2.0 + 1.0 / tem2
738 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
739 taup(i,kp1) = tem1 * hd * hd
741 taup(i,kp1) = taup(i,k) * rscor
743 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
751 taup(1:npt,km+1) = taup(1:npt,km)
757 taud(i,k) = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
778 IF (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev)
THEN
780 IF(taud(i,k) /= 0.)
THEN
781 tem = dtp * taud(i,k)
782 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
801 call oro_wam_2017(im, km, npt, ipt, kref, kdt, me, master,
802 & dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl,
803 & del, sigma, hprime, gamma, theta,
804 & sinlat, xlatd, taup, taud, pkdis)
814 axtms(:,:) = 0.0 ; aytms(:,:) = 0.0
819 zpbl =rgrav*phil( j, kpbl(j) )
821 sigflt = min(sgh30(j), 0.3*hprime(j))
823 zsurf = phii(j,1)*rgrav
825 zpm(k) = phil(j,k)*rgrav
830 call ugwpv0_tofd1d(km, sigflt, elvmaxd(j), zsurf, zpbl,
831 & up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
834 axtms(j,k) = utofd1(k)
835 aytms(j,k) = vtofd1(k)
839 pdvdt(j,k) = pdvdt(j,k) + aytms(j,k)
840 pdudt(j,k) = pdudt(j,k) + axtms(j,k)
843 tau_tofd(j) = sum( utofd1(1:km)* del(j,1:km))
860 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
862 if ( k < idxzb(i) .AND. idxzb(i) /= 0 )
then
866 dbim = db(i,k) / (1.+db(i,k)*dtp)
867 pdvdt(j,k) = - dbim * v1(j,k) +pdvdt(j,k)
868 pdudt(j,k) = - dbim * u1(j,k) +pdudt(j,k)
869 eng1 = eng0*(1.0-dbim*dtp)*(1.-dbim*dtp)
871 dusfc(j) = dusfc(j) - dbim * u1(j,k) * del(j,k)
872 dvsfc(j) = dvsfc(j) - dbim * v1(j,k) * del(j,k)
874 dudt_mtb(j,k) = -dbim * u1(j,k)
875 tau_mtb(j) = tau_mtb(j) + dudt_mtb(j,k)* del(j,k)
881 taud(i,k) = taud(i,k) * dtfac(i)
882 dtaux = taud(i,k) * xn(i)
883 dtauy = taud(i,k) * yn(i)
885 pdvdt(j,k) = dtauy +pdvdt(j,k)
886 pdudt(j,k) = dtaux +pdudt(j,k)
888 unew = u1(j,k) + dtaux*dtp
889 vnew = v1(j,k) + dtauy*dtp
890 eng1 = 0.5*(unew*unew + vnew*vnew)
892 dusfc(j) = dusfc(j) + dtaux * del(j,k)
893 dvsfc(j) = dvsfc(j) + dtauy * del(j,k)
895 dudt_ogw(j,k) = dtaux
896 tau_ogw(j) = tau_ogw(j) +dtaux*del(j,k)
901 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt
907 dusfc(j) = -rgrav * dusfc(j)
908 dvsfc(j) = -rgrav * dvsfc(j)
909 tau_mtb(j) = -rgrav * tau_mtb(j)
910 tau_ogw(j) = -rgrav * tau_ogw(j)
911 tau_tofd(j) = -rgrav * tau_tofd(j)
927 subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
928 & tm1 , um1, vm1, qm1,
929 & prsl, prsi, philg, xlatd, sinlat, coslat,
930 & pdudt, pdvdt, pdtdt, dked, tau_ngw,
931 & mpi_id, master, kdt)
941 use machine,
only : kind_phys
944 &, omega2, rcpd2, pi, pi2, fv
945 &, rad_to_deg, deg_to_rad
946 &, rdi, gor, grcp, gocp
947 &, bnv2min, dw2min, velmin, gr2
950 &, v_kxw, v_kxw2, tamp_mpa, zfluxglob
951 &, maxdudt, gw_eff, dked_min
952 &, nslope, ilaunch, zmsi
953 &, zci, zdci, zci4, zci3, zci2
954 &, zaz_fct, zcosang, zsinang
955 &, nwav, nazd, zcimin, zcimax
960 integer,
parameter :: kp = kind_phys
961 integer,
intent(in) :: klev
962 integer,
intent(in) :: klon
964 real,
intent(in) :: dtime
965 real,
intent(in) :: vm1(klon,klev)
966 real,
intent(in) :: um1(klon,klev)
967 real,
intent(in) :: qm1(klon,klev)
968 real,
intent(in) :: tm1(klon,klev)
970 real,
intent(in) :: prsl(klon,klev)
971 real,
intent(in) :: philg(klon,klev)
972 real,
intent(in) :: prsi(klon,klev+1)
973 real,
intent(in) :: xlatd(klon)
974 real,
intent(in) :: sinlat(klon)
975 real,
intent(in) :: coslat(klon)
976 real,
intent(in) :: tau_ngw(klon)
978 integer,
intent(in) :: mpi_id, master, kdt
983 real,
intent(out) :: pdudt(klon,klev)
984 real,
intent(out) :: pdvdt(klon,klev)
985 real,
intent(out) :: pdtdt(klon,klev)
986 real,
intent(out) :: dked(klon,klev)
987 real,
parameter :: minvel = 0.5_kp
988 real,
parameter :: epsln = 1.0e-12_kp
989 real,
parameter :: zero = 0.0_kp, one = 1.0_kp, half = 0.5_kp
993 real :: taux(klon,klev+1)
994 real :: tauy(klon,klev+1)
995 real :: phil(klon,klev)
1002 real :: zbvfhm1(klon,ilaunch:klev)
1003 real :: zbn2(klon,ilaunch:klev)
1004 real :: zrhohm1(klon,ilaunch:klev)
1005 real :: zuhm1(klon,ilaunch:klev)
1006 real :: zvhm1(klon,ilaunch:klev)
1007 real :: v_zmet(klon,ilaunch:klev)
1008 real :: vueff(klon,ilaunch:klev)
1013 real :: zul(klon,nazd)
1014 real :: zci_min(klon,nazd)
1016 real :: zact(klon, nwav, nazd)
1019 real :: zpu(klon,klev, nazd)
1021 real :: zfct(klon,klev)
1022 real :: zfnorm(klon)
1024 real :: zfluxlaun(klon)
1025 real :: zui(klon, klev,nazd)
1027 real :: zdfdz_v(klon,klev, nazd)
1028 real :: zflux(klon, nwav, nazd)
1030 real :: zflux_z (klon, nwav,klev)
1032 real :: vm_zflx_mode, vc_zflx_mode
1033 real :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2
1036 real :: zu, zcin, zcpeak, zcin4, zbvfl4
1037 real :: zcin2, zbvfl2, zcin3, zbvfl3, zcinc
1038 real :: zatmp, zfluxs, zdep, zfluxsq, zulm, zdft, ze1, ze2
1041 real :: zdelp,zrgpts
1042 real :: zthstd,zrhostd,zbvfstd
1043 real :: tvc1, tvm1, tem1, tem2, tem3
1044 real :: zhook_handle
1045 real :: delpi(klon,ilaunch:klev)
1049 real,
parameter :: rcpdl = cpd/grav
1050 &, grav2cpd = grav/rcpdl
1053 real :: expdis, fdis
1055 real :: v_kzi, v_kzw, v_cdp, v_wdp, sc, tx1
1057 integer :: j, k, inc, jk, jl, iazi
1067 phil(j,k) = philg(j,k) * rgrav
1075 zpu(jl,jk,iazi) = zero
1086 zci_min(jl,iazi) = zcimin
1091 do jk=max(ilaunch,2),klev
1093 tvc1 = tm1(jl,jk) * (one +fv*qm1(jl,jk))
1094 tvm1 = tm1(jl,jk-1) * (one +fv*qm1(jl,jk-1))
1096 zthm1 = 2.0_kp / (tvc1+tvm1)
1097 zuhm1(jl,jk) = half *(um1(jl,jk-1)+um1(jl,jk))
1098 zvhm1(jl,jk) = half *(vm1(jl,jk-1)+vm1(jl,jk))
1100 zrhohm1(jl,jk) = prsi(jl,jk)*rdi*zthm1
1101 zdelp = phil(jl,jk)-phil(jl,jk-1)
1102 v_zmet(jl,jk) = zdelp + zdelp
1103 delpi(jl,jk) = grav / (prsi(jl,jk-1) - prsi(jl,jk))
1105 & 2.e-5_kp*exp( (phil(jl,jk)+phil(jl,jk-1))*rhp2)+dked_min
1108 zbn2(jl,jk) = grav2cpd*zthm1
1109 & * (one+rcpdl*(tm1(jl,jk)-tm1(jl,jk-1))/zdelp)
1110 zbn2(jl,jk) = max(min(zbn2(jl,jk), gssec), bv2min)
1111 zbvfhm1(jl,jk) = sqrt(zbn2(jl,jk))
1115 if (ilaunch == 1)
then
1119 zuhm1(jl,jk) = um1(jl,jk)
1120 zvhm1(jl,jk) = vm1(jl,jk)
1121 zbvfhm1(jl,1) = zbvfhm1(jl,2)
1122 v_zmet(jl,1) = v_zmet(jl,2)
1123 vueff(jl,1) = dked_min
1124 zbn2(jl,1) = zbn2(jl,2)
1128 tx1 = omega2 * sinlat(jl) / v_kxw
1129 c2f2(jl) = tx1 * tx1
1130 zbvfl(jl) = zbvfhm1(jl,ilaunch)
1137 zul(jl,iazi) = zcosang(iazi) * zuhm1(jl,ilaunch)
1138 & + zsinang(iazi) * zvhm1(jl,ilaunch)
1142 do jk=ilaunch, klev-1
1145 zu = zcosang(iazi)*zuhm1(jl,jk)
1146 & + zsinang(iazi)*zvhm1(jl,jk)
1147 zui(jl,jk,iazi) = zu - zul(jl,iazi)
1154 do jk=ilaunch, klev-1
1156 zfct(jl,jk) = zrhohm1(jl,jk) / zbvfhm1(jl,jk)
1164 if(nslope == 1)
then
1171 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1172 zbvfl4 = zbvfl4 * zbvfl4
1173 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl4*zcin
1177 elseif(nslope == 2)
then
1183 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1184 zbvfl4 = zbvfl4 * zbvfl4
1185 zcpeak = zbvfl(jl) * zmsi
1186 zflux(jl,inc,1) = zfct(jl,ilaunch)*
1187 & zbvfl4*zcin*zcpeak/(zbvfl4*zcpeak+zcin4*zcin)
1190 elseif(nslope == -1)
then
1196 zbvfl2 = zbvfl(jl)*zbvfl(jl)
1197 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl2*zcin
1201 elseif(nslope == 0)
then
1208 zbvfl3 = zbvfl(jl)**3
1209 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl3*zcin
1222 zpu(jl,ilaunch,1) = zpu(jl,ilaunch,1) + zflux(jl,inc,1)*zcinc
1231 zfluxlaun(jl) = tau_ngw(jl)
1232 zfnorm(jl) = zfluxlaun(jl) / zpu(jl,ilaunch,1)
1237 zpu(jl,ilaunch,iazi) = zfluxlaun(jl)
1243 do jk=ilaunch, klev-1
1245 zfct(jl,jk) = zfnorm(jl)*zfct(jl,jk)
1252 zflux(jl,inc,1) = zfnorm(jl)*zflux(jl,inc,1)
1263 zflux(jl,inc,iazi) = zflux(jl,inc,1)
1276 do jk=ilaunch, klev-1
1280 zci_min(jl,iazi) = max(zci_min(jl,iazi),zui(jl,jk,iazi))
1290 zact(jl,inc,iazi) = minvel
1291 & + sign(minvel,zci(inc)-zci_min(jl,iazi))
1326 if (abs(zcin) > epsln)
then
1337 v_cdp = abs(zcin-zui(jl,jk,iazi))
1339 wdop2 = v_wdp* v_wdp
1340 cdf2 = v_cdp*v_cdp - c2f2(jl)
1341 if (cdf2 > zero)
then
1342 kzw2 = (zbn2(jl,jk)-wdop2)/cdf2 - v_kxw2
1346 if ( kzw2 > zero .and. cdf2 > zero)
then
1354 v_cdp = sqrt( cdf2 )
1355 v_wdp = v_kxw * v_cdp
1356 v_kzi = abs(v_kzw*v_kzw*vueff(jl,jk)/v_wdp*v_kzw)
1357 expdis = exp(-v_kzi*v_zmet(jl,jk))
1367 fdis = expdis * zflux(jl,inc,iazi)
1372 zfluxs = zfct(jl,jk)*v_cdp*v_cdp*zcinc
1377 zdep = zact(jl,inc,iazi)* (fdis-zfluxs)
1378 if(zdep > zero )
then
1380 zflux(jl,inc,iazi) = zfluxs
1381 zflux_z(jl,inc,jk) = zfluxs
1384 zflux(jl,inc,iazi) = fdis
1385 zflux_z(jl,inc,jk) = fdis
1392 zdfdz_v(:,jk,iazi) = zero
1397 vc_zflx_mode = zact(jl,inc,iazi)*zflux(jl,inc,iazi)
1398 zpu(jl,jk,iazi) = zpu(jl,jk,iazi) + vc_zflx_mode*zcinc
1405 if (jk > ilaunch)
then
1408 zdelp = delpi(jl,jk) * abs(zcin-zui(jl,jk,iazi)) *zcinc
1409 vm_zflx_mode = zact(jl,inc,iazi)* zflux_z(jl,inc,jk-1)
1411 if (vc_zflx_mode > vm_zflx_mode)
1412 & vc_zflx_mode = vm_zflx_mode
1413 zdfdz_v( jl,jk,iazi) = zdfdz_v( jl,jk,iazi) +
1414 & (vm_zflx_mode-vc_zflx_mode)*zdelp
1439 tem1 = zaz_fct*zcosang(iazi)
1440 tem2 = zaz_fct*zsinang(iazi)
1441 do jk=ilaunch, klev-1
1443 taux(jl,jk) = taux(jl,jk) + tem1 * zpu(jl,jk,iazi)
1444 tauy(jl,jk) = tauy(jl,jk) + tem2 * zpu(jl,jk,iazi)
1445 pdtdt(jl,jk) = pdtdt(jl,jk) + tem3 * zdfdz_v(jl,jk,iazi)
1458 zdelp = delpi(jl,jk)
1459 ze1 = (taux(jl,jk)-taux(jl,jk-1))*zdelp
1460 ze2 = (tauy(jl,jk)-tauy(jl,jk-1))*zdelp
1461 if (abs(ze1) >= maxdudt )
then
1462 ze1 = sign(maxdudt, ze1)
1464 if (abs(ze2) >= maxdudt )
then
1465 ze2 = sign(maxdudt, ze2)
1472 pdtdt(jl,jk) = (ze1*um1(jl,jk) + ze2*vm1(jl,jk)) * cpdi
1474 dked(jl,jk) = max(dked_min, pdtdt(jl,jk)/zbn2(jl,jk))
1483 pdudt(jl,jk) = gw_eff * pdudt(jl,jk)
1484 pdvdt(jl,jk) = gw_eff * pdvdt(jl,jk)
1485 pdtdt(jl,jk) = gw_eff * pdtdt(jl,jk)
1486 dked(jl,jk) = gw_eff * dked(jl,jk)