8 subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, &
9 xlatd, sinlat, coslat, sparea, &
10 cdmbgwd, hprime, oc, oa4, clx4, theta, sigmad, &
11 gammad, elvmaxd, sgh30, kpbl, &
12 u1 ,v1, t1, q1, prsi,del,prsl,prslk, zmeti, zmet, &
13 pdvdt, pdudt, pdtdt, pkdis, dusfc, dvsfc,rdxzb , &
14 zobl, zlwb, zogw, tau_ogw, dudt_ogw, dvdt_ogw, &
15 dudt_obl, dvdt_obl,dudt_ofd, dvdt_ofd, &
16 du_ogwcol, dv_ogwcol, du_oblcol, dv_oblcol, &
17 du_ofdcol, dv_ofdcol, errmsg,errflg )
45 use machine ,
only : kind_phys
46 use ugwp_common,
only : dw2min, velmin, grav, omega1, rd, cpd, rv, pi, arad, fv
47 use ugwp_common,
only : rcpdt, grav2, rgrav, rcpd, rcpd2
48 use ugwp_common,
only : rad_to_deg, deg_to_rad, pi2, pih, rdi, gor, grcp, gocp, gr2, bnv2min
51 hpmax, hpmin, sigfaci => sigfac, &
52 dpmin, minwnd, hminmt, hncrit, &
53 rlolev, gmax, veleps, factop, &
54 frc, ce, ceofrc, frmax, cg, &
56 cdmb, cleff, fcrit_v1, &
57 n_tofd, ze_tofd, ztop_tofd
66 real(kind=kind_phys),
parameter :: sigfac = 3
67 real(kind=kind_phys),
parameter :: sigfacs = 0.25
69 real(kind=kind_phys),
parameter :: dbmax = 1./3600./12.
70 character(len=8) :: strsolver=
'pss-1986'
73 real(kind=kind_phys) :: gammin = 0.00999999
74 real(kind=kind_phys),
parameter :: nhilmax = 15.
75 real(kind=kind_phys),
parameter :: sso_min = 3000.
77 real(kind=kind_phys),
parameter :: nfr = 2.+1.
78 real(kind=kind_phys),
parameter :: afr = 1.
79 real(kind=kind_phys),
parameter :: frnorm =afr+1.0
80 real(kind=kind_phys),
parameter :: max_frf =2.0
82 logical,
parameter :: do_adjoro = .false.
85 integer,
intent(in) :: im, km, imx, kdt
86 integer,
intent(in) :: me, master
87 logical,
intent(in) :: do_tofd
89 integer,
intent(in) :: kpbl(im)
90 real(kind=kind_phys),
intent(in) :: dtp
91 real(kind=kind_phys),
intent(in) :: cdmbgwd(2)
93 real(kind=kind_phys),
intent(in) :: hprime(im), oc(im), oa4(im,4), &
94 clx4(im,4), theta(im), &
95 sigmad(im), gammad(im), elvmaxd(im)
97 real(kind=kind_phys),
intent(in) :: sgh30(im)
99 real(kind=kind_phys),
intent(in),
dimension(im,km) :: &
100 u1, v1, t1, q1,del, prsl, prslk, zmet
102 real(kind=kind_phys),
intent(in),
dimension(im,km+1):: prsi, zmeti
104 real(kind=kind_phys),
intent(in) :: xlatd(im),sinlat(im), coslat(im)
105 real(kind=kind_phys),
intent(in) :: sparea(im)
110 real(kind=kind_phys),
dimension(im,km),
intent(out) :: &
111 pdvdt, pdudt, pkdis, pdtdt
113 real(kind=kind_phys),
dimension(im,km),
intent(out) :: &
114 dudt_ogw,dvdt_ogw, dudt_obl,dvdt_obl, dudt_ofd,dvdt_ofd
116 real(kind=kind_phys),
dimension(im),
intent(out) :: dusfc, dvsfc, &
117 du_ogwcol,dv_ogwcol, du_oblcol,dv_oblcol, du_ofdcol,dv_ofdcol
119 real(kind=kind_phys),
dimension(im),
intent(out) :: rdxzb
120 real(kind=kind_phys),
dimension(im),
intent(out) :: zobl, zogw, zlwb, tau_ogw
122 character(len=*),
intent(out) :: errmsg
123 integer,
intent(out) :: errflg
129 real(kind=kind_phys),
dimension(im) :: oa, clx
130 real(kind=kind_phys),
dimension(im) :: sigma, gamma, elvmax
132 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
133 real(kind=kind_phys) :: belpmin, dsmin, dsmax
135 real(kind=kind_phys) :: arhills(im), mkd05_hills(im)
136 real(kind=kind_phys) :: taub_kd05(im)
140 real(kind=kind_phys),
dimension(im,km) :: ri_n, bnv2, ro
141 real(kind=kind_phys),
dimension(im,km) :: vtk, vtj, velco
145 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
146 real(kind=kind_phys),
dimension(im) :: wk, pe, ek, up
148 real(kind=kind_phys),
dimension(im,km) :: db, ang, uds
150 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
151 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
152 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
160 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
161 real(kind=kind_phys),
dimension(km) :: utofd1, vtofd1
162 real(kind=kind_phys),
dimension(km) :: epstofd1, krf_tofd1
163 real(kind=kind_phys),
dimension(km) :: up1, vp1, zpm
167 real(kind=kind_phys) :: xlingfs
168 logical :: icrilv(im)
170 real(kind=kind_phys),
dimension(im) :: xn, yn, ubar, vbar, ulow, &
171 roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
173 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
174 real(kind=kind_phys) :: taub(im), taulin(im), tausat(im), ahdxres(im)
175 real(kind=kind_phys) :: heff, hsat, hdis
177 integer,
dimension(im) :: kref, idxzb, ipt, khtop, iwk, izlow
181 real(kind=kind_phys) :: bnv, fr, ri_gw, brvf, fr2
182 real(kind=kind_phys) :: tem, tem1, tem2, temc, temv
183 real(kind=kind_phys) :: ti, rdz, dw2, shr2, bvf2
184 real(kind=kind_phys) :: rdelks, efact, coefm, gfobnv
185 real(kind=kind_phys) :: scork, rscor, hd, fro, sira
186 real(kind=kind_phys) :: dtaux, dtauy, zmetp, zmetk
188 real(kind=kind_phys) :: windik, wdir
189 real(kind=kind_phys) :: sigmin, dxres,sigres,hdxres, cdmb4, mtbridge
191 real(kind=kind_phys) :: kxridge, inv_b2eff, zw1, zw2
192 real(kind=kind_phys) :: belps, aelps, nhills, selps
197 real(kind=kind_phys) :: cleff_max
198 real(kind=kind_phys) :: nonh_fact
199 real(kind=kind_phys) :: fcrit2
200 real(kind=kind_phys) :: fr_func, frnd
205 integer :: kmm1, kmm2, lcap, lcapp1
207 integer :: kmps, idir, nwd, klcap, kp1, kmpbl, kmll
208 integer :: k_mtb, k_zlow, ktrial, klevm1
235 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin )
then
259 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
260 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
263 cleff_max = pi2/max(dsmin/5.,sso_min)
264 cleff_max = pi2/dsmin
268 gammin = min(sso_min/hdxres, 1.)
269 gammin = max(0.1, gammin)
271 sigmin = hpmin/hdxres
303 dxres = sqrt(sparea(i))
305 if (gamma(i) > 1.0) gamma(i) = 1.0
307 gamma(i) = max(gammin, gamma(i))
311 sigres = max(sigmin, sigma(i))
313 aelps = min( hprime(i)/sigres, dxres)
314 belps = min(aelps/abs(gamma(i)), dxres)
315 gamma(i) = aelps/belps
321 if (hprime(i) > hdxres*sigres) sigres= hprime(i)/dxres
322 aelps = min( hprime(i)/sigres, hdxres)
324 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i), dxres)
329 if( aelps < sso_min )
then
334 if (belps < sso_min )
then
336 belps = aelps*gamma(i)
338 gamma(i) = min(aelps/belps, 1.0)
341 sigma(i) = hprime(i)/aelps
342 gamma(i) = min(aelps/belps, 1.0)
347 selps = belps*belps*gamma(i)*pi
348 nhills = min(nhilmax, sparea(i)/selps)
349 arhills(j) = max(nhills, 1.0)
374 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
375 lcap = km ; lcapp1 = lcap + 1
380 elvmax(j) = min( sigfac * hprime(j), hncrit)
401 ztoph = sigfac * hprime(j)
402 zlowh = sigfacs* hprime(j)
410 if (( ztoph <= zmetp) .and. (ztoph >= zmetk) ) khtop(i) = max(khtop(i), k+1 )
411 if (zlowh <= zmetp .and. zlowh >= zmetk) izlow(i) = max(izlow(i),k)
418 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
419 vtk(i,k) = vtj(i,k) / prslk(j,k)
420 ro(i,k) = rdi * prsl(j,k) / vtj(i,k)
430 rdz = 1. / (zmet(j,k+1) - zmet(j,k))
431 tem1 = u1(j,k) - u1(j,k+1)
432 tem2 = v1(j,k) - v1(j,k+1)
433 dw2 = tem1*tem1 + tem2*tem2
434 shr2 = max(dw2,dw2min) * rdz * rdz
435 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))/ (vtk(i,k+1)+vtk(i,k))
437 bnv2(i,k+1) = max( bvf2, bnv2min )
438 ri_n(i,k+1) = bnv2(i,k)/shr2
447 bnv2(i,k) = bnv2(i,k+1)
455 if (k_zlow == khtop(i)) k_zlow = 1
456 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,khtop(i)))
467 do k = k_zlow, khtop(i)-1
468 rdelks = del(j,k) * delks(i)
469 ubar(i) = ubar(i) + rdelks * u1(j,k)
470 vbar(i) = vbar(i) + rdelks * v1(j,k)
471 roll(i) = roll(i) + rdelks * ro(i,k)
472 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
483 do k = khtop(i), 1, -1
484 phiang = atan2(v1(j,k),u1(j,k))
485 phiang = theta(j)*rad_to_deg - phiang
486 if ( phiang > pih ) phiang = phiang - pi
487 if ( phiang < -pih ) phiang = phiang + pi
489 uds(i,k) = max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
491 if (idxzb(i) == 0 )
then
492 dz_blk = zmeti(j,k+1) - zmeti(j,k)
493 pe(i) = pe(i) + bnv2(i,k) *( elvmax(j) - zmet(j,k) ) * dz_blk
494 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
495 ek(i) = 0.5 * up(i) * up(i)
496 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
501 if(ph_blk >= fcrit_v1 )
then
504 rdxzb(j) = real(k, kind=kind_phys)
512 if ( idxzb(i) > 0 )
then
516 gam2 = gamma(j)*gamma(j)
517 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
518 cgam = 0.48*gamma(j) + 0.30*gam2
519 do k = idxzb(i)-1, 1, -1
523 zlen = sqrt( (zobl(j)-zmet(j,k) )/(zmet(j,k ) + hprime(j)) )
526 sinang2 = 1.0 - cosang2
527 rdem = cosang2 + gam2 * sinang2
528 rnom = cosang2*gam2 + sinang2
534 rdem = max(rdem, 1.e-6)
536 zr = max( 2. - r, 0. )
537 sigres = max(sigmin, sigma(j))
538 mtbridge = zr * sigres*zlen / hprime(j)
540 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam * sinang2)
545 db(i,k)= dbtmp * uds(i,k)
547 db(i,k)= min(db(i,k), dbmax)
572 tem = (prsi(j,1) - prsi(j,k))
573 if (tem < dpmin) iwk(i) = k
588 k_mtb = max(1, idxzb(i))
590 kref(i) = max(iwk(i), kpbl(j)+1 )
591 kref(i) = max(kref(i), khtop(i) )
593 if (kref(i) <= k_mtb) kref(i) = k_mtb + 1
594 kbps = max(kbps, kref(i))
595 kmps = min(kmps, kref(i))
597 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
609 k_mtb = max(1, idxzb(i))
613 rdelks = del(j,k) * delks(i)
614 ubar(i) = ubar(i) + rdelks * u1(j,k)
615 vbar(i) = vbar(i) + rdelks * v1(j,k)
616 roll(i) = roll(i) + rdelks * ro(i,k)
617 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
626 wdir = atan2(ubar(i),vbar(i)) + pi
627 idir = mod(nint(fdir*wdir),mdir) + 1
629 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
630 clx(i) = clx4(j,mod(nwd-1,4)+1)
642 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
643 xn(i) = ubar(i) / ulow(i)
644 yn(i) = vbar(i) / ulow(i)
650 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)+ (v1(j,k)+v1(j,k+1))*yn(i))
655 velco(i,km) = velco(i,kmm1)
670 taub(:) = 0. ; taulin(:)= 0. ;taub_kd05 =0.
671 fcrit2 =fcrit_v1*fcrit_v1
679 bnv = sqrt( bnv2bar(i) )
680 heff = min(hprime(j),hpmax)
681 if( zobl(j) > 0.) heff = max(sigfac*heff-zobl(j), 0.)/sigfac
685 heff = min(heff, hsat)
693 efact = (oa(i) + 2.) ** (ceofrc*fr)
694 efact = min( max(efact,efmin), efmax )
695 gfobnv = efact* gmax /(zw2 + cg)
700 mkd05_hills(i) = (1. + clx(i)) ** (oa(i)+1.)
701 xlinv(i) = min(cleff * mkd05_hills(i), cleff_max)
702 taub_kd05(i) = roll(i)*xlinv(i) *(gfobnv *zw2)* (zw1 * ulow(i)* ulow(i))
708 inv_b2eff = pi*sigres/heff
709 kxridge = pi /ahdxres(i)
710 xlingfs = max(inv_b2eff, kxridge)
711 nonh_fact = 1. - xlinv(i)*zw1 * xlinv(i)*zw1
713 if (nonh_fact <= 0.) cycle
715 taulin(i) = xlinv(i)*roll(i)*bnv*ulow(i)*heff*heff * nonh_fact
716 tausat(i) = xlinv(i)*roll(i)* zw1*ulow(i)*ulow(i) *fcrit2 * nonh_fact
721 if(fr > fcrit_v1 )
then
723 fr_func = frnorm* frnd*frnd/(afr + frnd ** nfr)
724 taub(i) = tausat(i) *max(fr_func, max_frf)
728 xlinv(i) = .5*xlinv(i)
729 k = max(1, kref(i)-1)
730 tem = max(velco(i,k)*velco(i,k), dw2min)
731 scor(i) = bnv2(i,k) / tem
732 zogw(j) = zmeti(j, kref(i) )
740 taup(i, 1:kref(i) ) = taub(i)
747 if (strsolver ==
'pss-1986')
then
760 if (k >= kref(i))
then
761 icrilv(i) = icrilv(i) .or. ( ri_n(i,k) < ric).or. (velco(i,k) <= 0. )
766 if (k >= kref(i))
then
767 if (.not.icrilv(i) .and. taup(i,k) > 0.0 )
then
768 zw1 = max(velco(i,k), velmin)
774 if (oa(i) > 0. .and. kp1 < kref(i))
then
775 scork = bnv2(i,k) * temv * temv
776 rscor = min(1.0, scork / scor(i))
782 brvf = sqrt(bnv2(i,k))
783 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k)) *brvf* zw1
784 hd = sqrt(taup(i,k) / tem1)
785 fro = brvf * hd * temv
789 tem2 = sqrt(ri_n(i,k))
790 tem = 1. + tem2 * fro
791 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
797 if (ri_gw <= ric .and.(oa(i) <= 0. .or. kp1 >= kref(i) ))
then
798 temc = 2.0 + 1.0 / tem2
799 hd = zw1 * (2.*sqrt(temc)-temc) / brvf
800 taup(i,kp1) = tem1 * hd * hd
802 taup(i,kp1) = taup(i,k) * rscor
805 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
813 taup(1:npt,km+1) = taup(1:npt,km)
819 zw1 = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
824 if (abs(zw1) > max_axyz ) zw1 = sign(max_axyz, zw1)
840 if (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev .and. taud(i,k) /= 0.)
then
841 tem = dtp * taud(i,k)
842 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
852 call oro_spectral_solver(im, km, npt, ipt, kref, kdt, me, master, &
853 dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl, &
854 del, sigma, hprime, gamma, theta, sinlat, xlatd, taup, taud, pkdis)
868 zpbl = zmet( j, kpbl(j) )
869 sigflt = min(sgh30(j), 0.3*hprime(j))
877 call ugwp_tofd1d(km, cpd, dtp, sigflt, zsurf, zpbl, &
878 up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
881 dudt_ofd(j,k) = utofd1(k)
882 dvdt_ofd(j,k) = vtofd1(k)
886 pdvdt(j,k) = pdvdt(j,k) + utofd1(k)
887 pdudt(j,k) = pdudt(j,k) + vtofd1(k)
888 pdtdt(j,k) = pdtdt(j,k) + epstofd1(k)
890 du_ofdcol(j) = sum( utofd1(1:km)* del(j,1:km))
891 dv_ofdcol(j) = sum( vtofd1(1:km)* del(j,1:km))
893 dusfc(j) = dusfc(j) + du_ofdcol(j)
894 dvsfc(j) = dvsfc(j) + dv_ofdcol(j)
905 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
906 if ( k < idxzb(i) .and. idxzb(i) /= 0 )
then
910 dbim = db(i,k) / (1.+db(i,k)*dtp)
911 dudt_obl(j,k) = -dbim * u1(j,k)
912 dvdt_obl(j,k) = -dbim * v1(j,k)
914 pdudt(j,k) = dudt_obl(j,k) +pdudt(j,k)
915 pdvdt(j,k) = dvdt_obl(j,k) +pdvdt(j,k)
916 du_oblcol(j) = du_oblcol(j) + dudt_obl(j,k)* del(j,k)
917 dv_oblcol(j) = dv_oblcol(j) + dvdt_obl(j,k)* del(j,k)
918 dusfc(j) = dusfc(j) + du_oblcol(j)
919 dvsfc(j) = dvsfc(j) + dv_oblcol(j)
925 taud(i,k) = taud(i,k) * dtfac(i)
926 dtaux = taud(i,k) * xn(i)
927 dtauy = taud(i,k) * yn(i)
929 dudt_ogw(j,k) = dtaux
930 dvdt_ogw(j,k) = dtauy
932 pdvdt(j,k) = dtauy +pdvdt(j,k)
933 pdudt(j,k) = dtaux +pdudt(j,k)
936 du_ogwcol(j) = du_ogwcol(j) + dtaux * del(j,k)
937 dv_ogwcol(j) = dv_ogwcol(j) + dtauy * del(j,k)
939 dusfc(j) = dusfc(j) + du_ogwcol(j)
940 dvsfc(j) = dvsfc(j) + dv_ogwcol(j)
946 unew = u1(j,k) + pdudt(j,k)*dtp
947 vnew = v1(j,k) + pdvdt(j,k)*dtp
948 eng1 = 0.5*(unew*unew + vnew*vnew)
949 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt + pdtdt(j,k)
956 dusfc(j) = -rgrav * dusfc(j)
957 dvsfc(j) = -rgrav * dvsfc(j)
958 du_ogwcol(j) = -rgrav *du_ogwcol(j)
959 dv_ogwcol(j) = -rgrav *dv_ogwcol(j)
960 du_oblcol(j) = -rgrav *du_oblcol(j)
961 dv_oblcol(j) = -rgrav *dv_oblcol(j)
962 du_ofdcol(j) = -rgrav * du_ofdcol(j)
963 dv_ofdcol(j) = -rgrav * du_ofdcol(j)
970 if (kdt <= 2 .and. me == 0)
then
971 print *,
'vgw-oro done gwdps_v0 in ugwp-v0 step-proc ', kdt, me
973 print *, maxval(pdudt)*86400., minval(pdudt)*86400,
'vgw_axoro'
974 print *, maxval(pdvdt)*86400., minval(pdvdt)*86400,
'vgw_ayoro'
976 print *, maxval(pdtdt)*86400., minval(pdtdt)*86400,
'vgw_epsoro'
978 print *, maxval(zogw),
' z_ogw ', maxval(tau_ogw),
' tau_ogw '
982 if (maxval(abs(pdudt))*86400. > 100.)
then
984 print *, maxval(u1), minval(u1),
' u1 gwdps-v1 '
985 print *, maxval(v1), minval(v1),
' v1 gwdps-v1 '
986 print *, maxval(t1), minval(t1),
' t1 gwdps-v1 '
987 print *, maxval(q1), minval(q1),
' q1 gwdps-v1 '
988 print *, maxval(del), minval(del),
' del gwdps-v1 '
989 print *, maxval(zmet),minval(zmet),
'zmet'
990 print *, maxval(zmeti),minval(zmeti),
'zmeti'
991 print *, maxval(prsi), minval(prsi),
' prsi '
992 print *, maxval(prsl), minval(prsl),
' prsl '
993 print *, maxval(ro), minval(ro),
' ro-dens '
994 print *, maxval(bnv2(1:npt,:)), minval(bnv2(1:npt,:)),
' bnv2 '
995 print *, maxval(kpbl), minval(kpbl),
' kpbl '
996 print *, maxval(sgh30), maxval(hprime), maxval(elvmax),
'oro-d'
1000 print *,zogw(j)/hprime(j), zobl(j)/hprime(j), &
1001 zmet(j,1)*1.e-3, nint(hprime(j)/sigma(j))
1006 errmsg =
'ERROR(orogw_v1): '