34 SUBROUTINE gwdps_v0(IM, km, imx, do_tofd,
35 & Pdvdt, Pdudt, Pdtdt, Pkdis, U1,V1,T1,Q1,KPBL,
36 & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DTP,KDT,
37 & sgh30, HPRIME,OC,OA4,CLX4,THETA,vSIGMA,vGAMMA,ELVMAXD,
38 & DUSFC, DVSFC, xlatd, sinlat, coslat, sparea,
39 & cdmbgwd, me, master, rdxzb, con_g, con_omega,
40 & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
41 & dudt_mtb, dudt_ogw, dudt_tms)
50 USE machine ,
ONLY : kind_phys
52 &, pi, rad_to_deg, deg_to_rad, pi2
53 &, rdi, gor, grcp, gocp, fv, gr2
54 &, bnv2min, dw2min, velmin, arad
57 &, hpmax, hpmin, sigfaci => sigfac
58 &, dpmin, minwnd, hminmt, hncrit
59 &, rlolev, gmax, veleps, factop
60 &, frc, ce, ceofrc, frmax, cg
62 &, cdmb, cleff, fcrit_gfs, fcrit_mtb
63 &, n_tofd, ze_tofd, ztop_tofd
69 integer,
parameter :: kp = kind_phys
70 character(len=8) :: strsolver=
'PSS-1986'
71 integer,
intent(in) :: im, km, imx, kdt
72 integer,
intent(in) :: me, master
73 logical,
intent(in) :: do_tofd
74 real(kind=kind_phys),
parameter :: sigfac = 3, sigfacs = 0.5
75 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
76 integer,
intent(in) :: KPBL(IM)
77 real(kind=kind_phys),
intent(in) :: dtp
78 real(kind=kind_phys),
intent(in) :: cdmbgwd(2)
80 real(kind=kind_phys),
intent(in),
dimension(im,km) ::
82 & del, prsl, prslk, phil
83 real(kind=kind_phys),
intent(in),
dimension(im,km+1):: prsi, phii
84 real(kind=kind_phys),
intent(in) :: xlatd(im),sinlat(im),
86 real(kind=kind_phys),
intent(in) :: sparea(im)
88 real(kind=kind_phys),
intent(in) :: oc(im), oa4(im,4), clx4(im,4)
89 real(kind=kind_phys),
intent(in) :: hprime(im), sgh30(im)
90 real(kind=kind_phys),
intent(in) :: elvmaxd(im), theta(im)
91 real(kind=kind_phys),
intent(in) :: vsigma(im), vgamma(im)
92 real(kind=kind_phys) :: sigma(im), gamma(im)
94 real(kind=kind_phys),
intent(in) :: con_g, con_omega
97 real(kind=kind_phys),
dimension(im,km),
intent(out) ::
98 & pdvdt, pdudt, pkdis, pdtdt
100 &, dudt_mtb, dudt_ogw, dudt_tms
102 real(kind=kind_phys),
dimension(im) :: rdxzb, zmtb, zogw
103 &, tau_ogw, tau_mtb, tau_tofd
117 real(kind=kind_phys) :: gammin = 0.00999999_kp
118 real(kind=kind_phys),
parameter :: nhilmax = 25.0_kp
119 real(kind=kind_phys),
parameter :: sso_min = 3000.0_kp
120 logical,
parameter :: do_adjoro = .true.
122 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
123 real(kind=kind_phys) :: belpmin, dsmin, dsmax
125 real(kind=kind_phys) :: xlingfs
130 real(kind=kind_phys),
dimension(im,km) :: ri_n, bnv2, ro
133 real(kind=kind_phys),
dimension(im) :: oa, clx , elvmax, wk
136 real(kind=kind_phys),
dimension(im,km) :: db, ang, uds
138 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
139 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
140 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
146 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
147 real(kind=kind_phys),
dimension(km) :: utofd1, vtofd1
148 &, epstofd1, krf_tofd1
150 real(kind=kind_phys),
dimension(im, km) :: axtms, aytms
156 real(kind=kind_phys),
dimension(im) :: xn, yn, ubar, vbar, ulow,
157 & roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
159 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
160 real(kind=kind_phys) :: taub(im), taulin(im), heff, hsat, hdis
162 integer,
dimension(im) :: kref, idxzb, ipt, kreflm,
167 real(kind=kind_phys) :: bnv, fr, ri_gw
168 &, brvf, tem, tem1, tem2, temc, temv
169 &, ti, rdz, dw2, shr2, bvf2
170 &, rdelks, efact, coefm, gfobnv
171 &, scork, rscor, hd, fro, sira
172 &, dtaux, dtauy, pkp1log, pklog
173 &, grav2, rcpdt, windik, wdir
174 &, sigmin, dxres,sigres,hdxres
176 &, kxridge, inv_b2eff, zw1, zw2
177 &, belps, aelps, nhills, selps
179 integer :: kmm1, kmm2, lcap, lcapp1
180 &, npt, kbps, kbpsp1,kbpsm1
181 &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
182 &, k_mtb, k_zlow, ktrial, klevm1, i, j, k
184 rcpdt = 1.0 / (cpd*dtp)
189 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
190 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
192 dxres = pi2*arad/float(imx)
193 hdxres = 0.5_kp*dxres
197 gammin = min(sso_min/dxres, 1.)
200 sigmin = 2.*hpmin/dxres
204 kxridge = float(imx)/arad * cdmbgwd(2)
218 sigma(i) = max(vsigma(i), sigmin)
219 gamma(i) = max(vgamma(i), gammin)
238 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin )
then
244 sigres = max(sigmin, sigma(i))
246 dxres = sqrt(sparea(i))
247 if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
248 aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
249 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i),.5*dxres)
253 if( aelps < sso_min .and. do_adjoro)
then
258 if (belps < sso_min )
then
260 belps = aelps*gamma(i)
262 gamma(i) = min(aelps/belps, 1.0)
264 sigma(i) = 2.*hprime(i)/aelps
265 gamma(i) = min(aelps/belps, 1.0)
268 selps = belps*belps*gamma(i)*4.
269 nhills = min(nhilmax, sparea(i)/selps)
299 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
300 lcap = km ; lcapp1 = lcap + 1
304 elvmax(j) = min(elvmaxd(j)*0. + sigfac * hprime(j), hncrit)
311 ztoph = sigfac * hprime(j)
312 zlowh = sigfacs* hprime(j)
313 pkp1log = phil(j,k+1) * rgrav
314 pklog = phil(j,k) * rgrav
317 if (( ztoph <= pkp1log) .and. (ztoph >= pklog) )
318 & iwklm(i) = max(iwklm(i), k+1 )
320 if (zlowh <= pkp1log .and. zlowh >= pklog)
321 & izlow(i) = max(izlow(i),k)
328 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
329 vtk(i,k) = vtj(i,k) / prslk(j,k)
330 ro(i,k) = rdi * prsl(j,k) / vtj(i,k)
340 rdz = grav / (phil(j,k+1) - phil(j,k))
341 tem1 = u1(j,k) - u1(j,k+1)
342 tem2 = v1(j,k) - v1(j,k+1)
343 dw2 = tem1*tem1 + tem2*tem2
344 shr2 = max(dw2,dw2min) * rdz * rdz
349 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))
350 & / (vtk(i,k+1)+vtk(i,k))
351 bnv2(i,k+1) = max( bvf2, bnv2min )
352 ri_n(i,k+1) = bnv2(i,k)/shr2
360 bnv2(i,k) = bnv2(i,k+1)
368 if (k_zlow == iwklm(i)) k_zlow = 1
369 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,iwklm(i)))
381 if (k_zlow == iwklm(i)) k_zlow = 1
382 DO k = k_zlow, iwklm(i)-1
384 rdelks = del(j,k) * delks(i)
385 ubar(i) = ubar(i) + rdelks * u1(j,k)
386 vbar(i) = vbar(i) + rdelks * v1(j,k)
387 roll(i) = roll(i) + rdelks * ro(i,k)
389 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
400 DO k = iwklm(i), 1, -1
401 phiang = atan2(v1(j,k),u1(j,k))*rad_to_deg
402 ang(i,k) = ( theta(j) - phiang )
403 if ( ang(i,k) > 90. ) ang(i,k) = ang(i,k) - 180.
404 if ( ang(i,k) < -90. ) ang(i,k) = ang(i,k) + 180.
405 ang(i,k) = ang(i,k) * deg_to_rad
407 & max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
409 IF (idxzb(i) == 0 )
then
410 dz_blk = ( phii(j,k+1) - phii(j,k) ) *rgrav
411 pe(i) = pe(i) + bnv2(i,k) *
412 & ( elvmax(j) - phil(j,k)*rgrav ) * dz_blk
414 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
415 ek(i) = 0.5 * up(i) * up(i)
417 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
421 IF ( ph_blk >= fcrit_gfs )
THEN
423 zmtb(j) = phil(j, k)*rgrav
424 rdxzb(j) = real(k, kind=kind_phys)
435 bnv = sqrt( bnv2bar(i) )
436 heff = 2.*min(hprime(j),hpmax)
437 zw2 = ubar(i)*ubar(i)+vbar(i)*vbar(i)
438 ulow(i) = sqrt(max(zw2,dw2min))
439 fr = heff*bnv/ulow(i)
440 zw1 = max(heff*(1. -fcrit_gfs/fr), 0.0)
441 zw2 = phil(j,2)*rgrav
442 if (fr > fcrit_gfs .and. zw1 > zw2 )
then
444 pkp1log = phil(j,k+1) * rgrav
445 pklog = phil(j,k) * rgrav
446 if (zw1 <= pkp1log .and. zw1 >= pklog)
exit
449 zmtb(j) = phil(j, k)*rgrav
464 IF ( idxzb(i) > 0 )
then
466 gam2 = gamma(j)*gamma(j)
467 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
468 cgam = 0.48*gamma(j) + 0.30*gam2
469 DO k = idxzb(i)-1, 1, -1
471 zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) /
472 & ( phil(j,k ) + grav * hprime(j) ) )
476 sinang2 = 1.0 - cosang2
481 rdem = cosang2 + gam2 * sinang2
482 rnom = cosang2*gam2 + sinang2
489 rdem = max(rdem, 1.e-6)
491 zr = max( 2. - r, 0. )
493 sigres = max(sigmin, sigma(j))
494 if (hprime(j)/sigres > dxres) sigres = hprime(j)/dxres
495 mtbridge = zr * sigres*zlen / hprime(j)
500 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam* sinang2)
501 db(i,k)= dbtmp * uds(i,k)
527 tem = (prsi(j,1) - prsi(j,k))
528 if (tem < dpmin) iwk(i) = k
549 k_mtb = max(1, idxzb(i))
551 kref(i) = max(iwk(i), kpbl(j)+1 )
552 kref(i) = max(kref(i), iwklm(i) )
554 if (kref(i) <= idxzb(i)) kref(i) = idxzb(i) + 1
555 kbps = max(kbps, kref(i))
556 kmps = min(kmps, kref(i))
558 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
570 k_mtb = max(1, idxzb(i))
572 IF (k < kref(i))
THEN
574 rdelks = del(j,k) * delks(i)
575 ubar(i) = ubar(i) + rdelks * u1(j,k)
576 vbar(i) = vbar(i) + rdelks * v1(j,k)
577 roll(i) = roll(i) + rdelks * ro(i,k)
578 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
586 wdir = atan2(ubar(i),vbar(i)) + pi
587 idir = mod(nint(fdir*wdir),mdir) + 1
589 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
590 clx(i) = clx4(j,mod(nwd-1,4)+1)
596 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
597 xn(i) = ubar(i) / ulow(i)
598 yn(i) = vbar(i) / ulow(i)
604 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)
605 & + (v1(j,k)+v1(j,k+1))*yn(i))
620 taub(:) = 0. ; taulin(:)= 0.
623 bnv = sqrt( bnv2bar(i) )
624 heff = min(hprime(j),hpmax)
626 if( zmtb(j) > 0.) heff = max(sigfac*heff-zmtb(j), 0.)/sigfac
629 hsat = fcrit_gfs*ulow(i)/bnv
630 heff = min(heff, hsat)
632 fr = min(bnv * heff /ulow(i), frmax)
634 efact = (oa(i) + 2.) ** (ceofrc*fr)
635 efact = min( max(efact,efmin), efmax )
637 coefm = (1. + clx(i)) ** (oa(i)+1.)
639 xlinv(i) = coefm * cleff
640 xlingfs = coefm * cleff
642 tem = fr * fr * oc(j)
643 gfobnv = gmax * tem / ((tem + cg)*bnv)
647 sigres = max(sigmin, sigma(j))
648 if (heff/sigres > hdxres) sigres = heff/hdxres
649 inv_b2eff = 0.5*sigres/heff
650 kxridge = 1.0 / sqrt(sparea(j))
652 taulin(i) = 0.5*roll(i)*xlinv(i)*bnv*ulow(i)*
655 if ( fr > fcrit_gfs )
then
656 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
657 & * ulow(i) * gfobnv * efact
661 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
662 & * ulow(i) * gfobnv * efact
669 k = max(1, kref(i)-1)
670 tem = max(velco(i,k)*velco(i,k), dw2min)
671 scor(i) = bnv2(i,k) / tem
675 zogw(j) = phii(j, kref(i)) *rgrav
682 IF (k <= kref(i)) taup(i,k) = taub(i)
686 if (strsolver ==
'PSS-1986')
then
700 IF (k >= kref(i))
THEN
701 icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) < ric)
702 & .OR. (velco(i,k) <= 0.0)
707 IF (k >= kref(i))
THEN
708 IF (.NOT.icrilv(i) .AND. taup(i,k) > 0.0 )
THEN
709 temv = 1.0 / max(velco(i,k), velmin)
711 IF (oa(i) > 0. .AND. kp1 < kref(i))
THEN
712 scork = bnv2(i,k) * temv * temv
713 rscor = min(1.0, scork / scor(i))
719 brvf = sqrt(bnv2(i,k))
722 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
723 & * max(velco(i,k), velmin)
724 hd = sqrt(taup(i,k) / tem1)
725 fro = brvf * hd * temv
730 tem2 = sqrt(ri_n(i,k))
731 tem = 1. + tem2 * fro
732 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
737 IF (ri_gw <= ric .AND.
738 & (oa(i) <= 0. .OR. kp1 >= kref(i) ))
THEN
739 temc = 2.0 + 1.0 / tem2
740 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
741 taup(i,kp1) = tem1 * hd * hd
743 taup(i,kp1) = taup(i,k) * rscor
745 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
753 taup(1:npt,km+1) = taup(1:npt,km)
759 taud(i,k) = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
780 IF (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev)
THEN
782 IF(taud(i,k) /= 0.)
THEN
783 tem = dtp * taud(i,k)
784 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
803 call oro_wam_2017(im, km, npt, ipt, kref, kdt, me, master,
804 & dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl,
805 & del, sigma, hprime, gamma, theta,
806 & sinlat, xlatd, taup, taud, pkdis)
816 axtms(:,:) = 0.0 ; aytms(:,:) = 0.0
821 zpbl =rgrav*phil( j, kpbl(j) )
823 sigflt = min(sgh30(j), 0.3*hprime(j))
825 zsurf = phii(j,1)*rgrav
827 zpm(k) = phil(j,k)*rgrav
832 call ugwpv0_tofd1d(km, sigflt, elvmaxd(j), zsurf, zpbl,
833 & up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
836 axtms(j,k) = utofd1(k)
837 aytms(j,k) = vtofd1(k)
841 pdvdt(j,k) = pdvdt(j,k) + aytms(j,k)
842 pdudt(j,k) = pdudt(j,k) + axtms(j,k)
845 tau_tofd(j) = sum( utofd1(1:km)* del(j,1:km))
862 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
864 if ( k < idxzb(i) .AND. idxzb(i) /= 0 )
then
868 dbim = db(i,k) / (1.+db(i,k)*dtp)
869 pdvdt(j,k) = - dbim * v1(j,k) +pdvdt(j,k)
870 pdudt(j,k) = - dbim * u1(j,k) +pdudt(j,k)
871 eng1 = eng0*(1.0-dbim*dtp)*(1.-dbim*dtp)
873 dusfc(j) = dusfc(j) - dbim * u1(j,k) * del(j,k)
874 dvsfc(j) = dvsfc(j) - dbim * v1(j,k) * del(j,k)
876 dudt_mtb(j,k) = -dbim * u1(j,k)
877 tau_mtb(j) = tau_mtb(j) + dudt_mtb(j,k)* del(j,k)
883 taud(i,k) = taud(i,k) * dtfac(i)
884 dtaux = taud(i,k) * xn(i)
885 dtauy = taud(i,k) * yn(i)
887 pdvdt(j,k) = dtauy +pdvdt(j,k)
888 pdudt(j,k) = dtaux +pdudt(j,k)
890 unew = u1(j,k) + dtaux*dtp
891 vnew = v1(j,k) + dtauy*dtp
892 eng1 = 0.5*(unew*unew + vnew*vnew)
894 dusfc(j) = dusfc(j) + dtaux * del(j,k)
895 dvsfc(j) = dvsfc(j) + dtauy * del(j,k)
897 dudt_ogw(j,k) = dtaux
898 tau_ogw(j) = tau_ogw(j) +dtaux*del(j,k)
903 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt
909 dusfc(j) = -rgrav * dusfc(j)
910 dvsfc(j) = -rgrav * dvsfc(j)
911 tau_mtb(j) = -rgrav * tau_mtb(j)
912 tau_ogw(j) = -rgrav * tau_ogw(j)
913 tau_tofd(j) = -rgrav * tau_tofd(j)
929 subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
930 & tm1 , um1, vm1, qm1,
931 & prsl, prsi, philg, xlatd, sinlat, coslat,
932 & pdudt, pdvdt, pdtdt, dked, tau_ngw,
933 & mpi_id, master, kdt)
943 use machine,
only : kind_phys
946 &, omega2, rcpd2, pi, pi2, fv
947 &, rad_to_deg, deg_to_rad
948 &, rdi, gor, grcp, gocp
949 &, bnv2min, dw2min, velmin, gr2
952 &, v_kxw, v_kxw2, tamp_mpa, zfluxglob
953 &, maxdudt, gw_eff, dked_min
954 &, nslope, ilaunch, zmsi
955 &, zci, zdci, zci4, zci3, zci2
956 &, zaz_fct, zcosang, zsinang
957 &, nwav, nazd, zcimin, zcimax
962 integer,
parameter :: kp = kind_phys
963 integer,
intent(in) :: klev
964 integer,
intent(in) :: klon
966 real,
intent(in) :: dtime
967 real,
intent(in) :: vm1(klon,klev)
968 real,
intent(in) :: um1(klon,klev)
969 real,
intent(in) :: qm1(klon,klev)
970 real,
intent(in) :: tm1(klon,klev)
972 real,
intent(in) :: prsl(klon,klev)
973 real,
intent(in) :: philg(klon,klev)
974 real,
intent(in) :: prsi(klon,klev+1)
975 real,
intent(in) :: xlatd(klon)
976 real,
intent(in) :: sinlat(klon)
977 real,
intent(in) :: coslat(klon)
978 real,
intent(in) :: tau_ngw(klon)
980 integer,
intent(in) :: mpi_id, master, kdt
985 real,
intent(out) :: pdudt(klon,klev)
986 real,
intent(out) :: pdvdt(klon,klev)
987 real,
intent(out) :: pdtdt(klon,klev)
988 real,
intent(out) :: dked(klon,klev)
989 real,
parameter :: minvel = 0.5_kp
990 real,
parameter :: epsln = 1.0e-12_kp
991 real,
parameter :: zero = 0.0_kp, one = 1.0_kp, half = 0.5_kp
995 real :: taux(klon,klev+1)
996 real :: tauy(klon,klev+1)
997 real :: phil(klon,klev)
1004 real :: zbvfhm1(klon,ilaunch:klev)
1005 real :: zbn2(klon,ilaunch:klev)
1006 real :: zrhohm1(klon,ilaunch:klev)
1007 real :: zuhm1(klon,ilaunch:klev)
1008 real :: zvhm1(klon,ilaunch:klev)
1009 real :: v_zmet(klon,ilaunch:klev)
1010 real :: vueff(klon,ilaunch:klev)
1015 real :: zul(klon,nazd)
1016 real :: zci_min(klon,nazd)
1018 real :: zact(klon, nwav, nazd)
1021 real :: zpu(klon,klev, nazd)
1023 real :: zfct(klon,klev)
1024 real :: zfnorm(klon)
1026 real :: zfluxlaun(klon)
1027 real :: zui(klon, klev,nazd)
1029 real :: zdfdz_v(klon,klev, nazd)
1030 real :: zflux(klon, nwav, nazd)
1032 real :: zflux_z (klon, nwav,klev)
1034 real :: vm_zflx_mode, vc_zflx_mode
1035 real :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2
1038 real :: zu, zcin, zcpeak, zcin4, zbvfl4
1039 real :: zcin2, zbvfl2, zcin3, zbvfl3, zcinc
1040 real :: zatmp, zfluxs, zdep, zfluxsq, zulm, zdft, ze1, ze2
1043 real :: zdelp,zrgpts
1044 real :: zthstd,zrhostd,zbvfstd
1045 real :: tvc1, tvm1, tem1, tem2, tem3
1046 real :: zhook_handle
1047 real :: delpi(klon,ilaunch:klev)
1051 real,
parameter :: rcpdl = cpd/grav
1052 &, grav2cpd = grav/rcpdl
1055 real :: expdis, fdis
1057 real :: v_kzi, v_kzw, v_cdp, v_wdp, sc, tx1
1059 integer :: j, k, inc, jk, jl, iazi
1069 phil(j,k) = philg(j,k) * rgrav
1077 zpu(jl,jk,iazi) = zero
1088 zci_min(jl,iazi) = zcimin
1093 do jk=max(ilaunch,2),klev
1095 tvc1 = tm1(jl,jk) * (one +fv*qm1(jl,jk))
1096 tvm1 = tm1(jl,jk-1) * (one +fv*qm1(jl,jk-1))
1098 zthm1 = 2.0_kp / (tvc1+tvm1)
1099 zuhm1(jl,jk) = half *(um1(jl,jk-1)+um1(jl,jk))
1100 zvhm1(jl,jk) = half *(vm1(jl,jk-1)+vm1(jl,jk))
1102 zrhohm1(jl,jk) = prsi(jl,jk)*rdi*zthm1
1103 zdelp = phil(jl,jk)-phil(jl,jk-1)
1104 v_zmet(jl,jk) = zdelp + zdelp
1105 delpi(jl,jk) = grav / (prsi(jl,jk-1) - prsi(jl,jk))
1107 & 2.e-5_kp*exp( (phil(jl,jk)+phil(jl,jk-1))*rhp2)+dked_min
1110 zbn2(jl,jk) = grav2cpd*zthm1
1111 & * (one+rcpdl*(tm1(jl,jk)-tm1(jl,jk-1))/zdelp)
1112 zbn2(jl,jk) = max(min(zbn2(jl,jk), gssec), bv2min)
1113 zbvfhm1(jl,jk) = sqrt(zbn2(jl,jk))
1117 if (ilaunch == 1)
then
1121 zuhm1(jl,jk) = um1(jl,jk)
1122 zvhm1(jl,jk) = vm1(jl,jk)
1123 zbvfhm1(jl,1) = zbvfhm1(jl,2)
1124 v_zmet(jl,1) = v_zmet(jl,2)
1125 vueff(jl,1) = dked_min
1126 zbn2(jl,1) = zbn2(jl,2)
1130 tx1 = omega2 * sinlat(jl) / v_kxw
1131 c2f2(jl) = tx1 * tx1
1132 zbvfl(jl) = zbvfhm1(jl,ilaunch)
1139 zul(jl,iazi) = zcosang(iazi) * zuhm1(jl,ilaunch)
1140 & + zsinang(iazi) * zvhm1(jl,ilaunch)
1144 do jk=ilaunch, klev-1
1147 zu = zcosang(iazi)*zuhm1(jl,jk)
1148 & + zsinang(iazi)*zvhm1(jl,jk)
1149 zui(jl,jk,iazi) = zu - zul(jl,iazi)
1156 do jk=ilaunch, klev-1
1158 zfct(jl,jk) = zrhohm1(jl,jk) / zbvfhm1(jl,jk)
1166 if(nslope == 1)
then
1173 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1174 zbvfl4 = zbvfl4 * zbvfl4
1175 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl4*zcin
1179 elseif(nslope == 2)
then
1185 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1186 zbvfl4 = zbvfl4 * zbvfl4
1187 zcpeak = zbvfl(jl) * zmsi
1188 zflux(jl,inc,1) = zfct(jl,ilaunch)*
1189 & zbvfl4*zcin*zcpeak/(zbvfl4*zcpeak+zcin4*zcin)
1192 elseif(nslope == -1)
then
1198 zbvfl2 = zbvfl(jl)*zbvfl(jl)
1199 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl2*zcin
1203 elseif(nslope == 0)
then
1210 zbvfl3 = zbvfl(jl)**3
1211 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl3*zcin
1224 zpu(jl,ilaunch,1) = zpu(jl,ilaunch,1) + zflux(jl,inc,1)*zcinc
1233 zfluxlaun(jl) = tau_ngw(jl)
1234 zfnorm(jl) = zfluxlaun(jl) / zpu(jl,ilaunch,1)
1239 zpu(jl,ilaunch,iazi) = zfluxlaun(jl)
1245 do jk=ilaunch, klev-1
1247 zfct(jl,jk) = zfnorm(jl)*zfct(jl,jk)
1254 zflux(jl,inc,1) = zfnorm(jl)*zflux(jl,inc,1)
1265 zflux(jl,inc,iazi) = zflux(jl,inc,1)
1278 do jk=ilaunch, klev-1
1282 zci_min(jl,iazi) = max(zci_min(jl,iazi),zui(jl,jk,iazi))
1292 zact(jl,inc,iazi) = minvel
1293 & + sign(minvel,zci(inc)-zci_min(jl,iazi))
1328 if (abs(zcin) > epsln)
then
1339 v_cdp = abs(zcin-zui(jl,jk,iazi))
1341 wdop2 = v_wdp* v_wdp
1342 cdf2 = v_cdp*v_cdp - c2f2(jl)
1343 if (cdf2 > zero)
then
1344 kzw2 = (zbn2(jl,jk)-wdop2)/cdf2 - v_kxw2
1348 if ( kzw2 > zero .and. cdf2 > zero)
then
1356 v_cdp = sqrt( cdf2 )
1357 v_wdp = v_kxw * v_cdp
1358 v_kzi = abs(v_kzw*v_kzw*vueff(jl,jk)/v_wdp*v_kzw)
1359 expdis = exp(-v_kzi*v_zmet(jl,jk))
1369 fdis = expdis * zflux(jl,inc,iazi)
1374 zfluxs = zfct(jl,jk)*v_cdp*v_cdp*zcinc
1379 zdep = zact(jl,inc,iazi)* (fdis-zfluxs)
1380 if(zdep > zero )
then
1382 zflux(jl,inc,iazi) = zfluxs
1383 zflux_z(jl,inc,jk) = zfluxs
1386 zflux(jl,inc,iazi) = fdis
1387 zflux_z(jl,inc,jk) = fdis
1394 zdfdz_v(:,jk,iazi) = zero
1399 vc_zflx_mode = zact(jl,inc,iazi)*zflux(jl,inc,iazi)
1400 zpu(jl,jk,iazi) = zpu(jl,jk,iazi) + vc_zflx_mode*zcinc
1407 if (jk > ilaunch)
then
1410 zdelp = delpi(jl,jk) * abs(zcin-zui(jl,jk,iazi)) *zcinc
1411 vm_zflx_mode = zact(jl,inc,iazi)* zflux_z(jl,inc,jk-1)
1413 if (vc_zflx_mode > vm_zflx_mode)
1414 & vc_zflx_mode = vm_zflx_mode
1415 zdfdz_v( jl,jk,iazi) = zdfdz_v( jl,jk,iazi) +
1416 & (vm_zflx_mode-vc_zflx_mode)*zdelp
1441 tem1 = zaz_fct*zcosang(iazi)
1442 tem2 = zaz_fct*zsinang(iazi)
1443 do jk=ilaunch, klev-1
1445 taux(jl,jk) = taux(jl,jk) + tem1 * zpu(jl,jk,iazi)
1446 tauy(jl,jk) = tauy(jl,jk) + tem2 * zpu(jl,jk,iazi)
1447 pdtdt(jl,jk) = pdtdt(jl,jk) + tem3 * zdfdz_v(jl,jk,iazi)
1460 zdelp = delpi(jl,jk)
1461 ze1 = (taux(jl,jk)-taux(jl,jk-1))*zdelp
1462 ze2 = (tauy(jl,jk)-tauy(jl,jk-1))*zdelp
1463 if (abs(ze1) >= maxdudt )
then
1464 ze1 = sign(maxdudt, ze1)
1466 if (abs(ze2) >= maxdudt )
then
1467 ze2 = sign(maxdudt, ze2)
1474 pdtdt(jl,jk) = (ze1*um1(jl,jk) + ze2*vm1(jl,jk)) * cpdi
1476 dked(jl,jk) = max(dked_min, pdtdt(jl,jk)/zbn2(jl,jk))
1485 pdudt(jl,jk) = gw_eff * pdudt(jl,jk)
1486 pdvdt(jl,jk) = gw_eff * pdvdt(jl,jk)
1487 pdtdt(jl,jk) = gw_eff * pdtdt(jl,jk)
1488 dked(jl,jk) = gw_eff * dked(jl,jk)