337 SUBROUTINE gwdps(IM,IX,IY,KM,A,B,C,U1,V1,T1,Q1,KPBL, &
338 & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT, &
339 & HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX, &
340 & DUSFC,DVSFC,G, CP, RD, RV, IMX, &
341 & nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb)
434 USE machine
, ONLY : kind_phys
436 integer im, iy, ix, km, imx, kdt, ipr, me
438 real(kind=kind_phys) deltim, G, CP, RD, RV, cdmbgwd(2)
439 real(kind=kind_phys) A(iy,km), B(iy,km), C(iy,km), &
440 & U1(IX,KM), V1(IX,KM), T1(IX,KM), &
441 & Q1(IX,KM), PRSI(IX,KM+1), DEL(IX,KM), &
442 & PRSL(IX,KM), PRSLK(IX,KM), PHIL(IX,KM), &
443 & PHII(IX,KM+1),RDXZB(IY)
444 real(kind=kind_phys) OC(im), OA4(iy,4), CLX4(iy,4) &
447 real(kind=kind_phys) ELVMAX(im),THETA(im),SIGMA(im),GAMMA(im)
448 real(kind=kind_phys) wk(im)
449 real(kind=kind_phys) bnv2lm(im,km),PE(im),EK(im),ZBK(im),UP(im)
450 real(kind=kind_phys) DB(im,km),ANG(im,km),UDS(im,km)
451 real(kind=kind_phys) ZLEN, DBTMP, R, PHIANG, CDmb, DBIM
452 real(kind=kind_phys) ENG0, ENG1
456 real(kind=kind_phys) pi, dw2min, rimin, ric, bnv2min, efmin
457 &, efmax,hpmax,hpmin, rad_to_deg, deg_to_rad
458 parameter(pi=3.1415926535897931)
459 parameter(rad_to_deg=180.0/pi, deg_to_rad=pi/180.0)
460 parameter(dw2min=1., rimin=-100., ric=0.25, bnv2min=1.0e-5)
462 parameter(efmin=0.0, efmax=10.0, hpmax=2400.0, hpmin=1.0)
464 real(kind=kind_phys) FRC, CE, CEOFRC, frmax, CG, GMAX
465 &, veleps, factop, rlolev, rdi
467 parameter(frc=1.0, ce=0.8, ceofrc=ce/frc, frmax=100., cg=0.5)
468 parameter(gmax=1.0, veleps=1.0, factop=0.5)
470 parameter(rlolev=50000.0)
474 real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac
477 parameter(hncrit=8000.)
479 parameter(sigfac=4.0)
480 parameter(hminmt=50.)
481 parameter(minwnd=0.1)
487 parameter(dpmin=5000.0)
490 real(kind=kind_phys) FDIR
492 parameter(mdir=8, fdir=mdir/(pi+pi))
494 data nwdir/6,7,5,8,2,3,1,4/
501 real(kind=kind_phys) TAUB(im), XN(im), YN(im), UBAR(im) &
502 &, VBAR(IM), ULOW(IM), OA(IM), CLX(IM) &
503 &, ROLL(IM), ULOI(IM), DUSFC(IM), DVSFC(IM) &
504 &, DTFAC(IM), XLINV(IM), DELKS(IM), DELKS1(IM)
506 real(kind=kind_phys) BNV2(im,km), TAUP(im,km+1), ri_n(im,km) &
507 &, TAUD(IM,KM), RO(IM,KM), VTK(IM,KM) &
508 &, VTJ(IM,KM), SCOR(IM), VELCO(IM,KM-1) &
512 Integer kref(im), kint(im), iwk(im), ipt(im)
514 Integer kreflm(im), iwklm(im)
515 Integer idxzb(im), ktrial, klevm1, nmtvr
517 real(kind=kind_phys) gor, gocp, fv, gr2, bnv, fr &
518 &, brvf, cleff, tem, tem1, tem2, temc, temv &
519 &, wdir, ti, rdz, dw2, shr2, bvf2 &
520 &, rdelks, efact, coefm, gfobnv &
521 &, scork, rscor, hd, fro, rim, sira &
522 &, dtaux, dtauy, pkp1log, pklog
523 integer kmm1, kmm2, lcap, lcapp1, kbps, kbpsp1,kbpsm1 &
524 &, kmps, idir, nwd, i, j, k, klcap, kp1, kmpbl, npt, npr &
533 cdmb = 4.0 * 192.0/float(imx)
534 if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
563 IF ( nmtvr .eq. 14)
then 569 IF ( (elvmax(i) .GT. hminmt)
570 & .and. (hprime(i) .GT. hpmin) )
then 573 if (ipr .eq. i) npr = npt
576 IF (npt .eq. 0)
RETURN 608 elvmax(j) = min(elvmax(j) + sigfac * hprime(j), hncrit)
616 pkp1log = phil(j,k+1) / g
617 pklog = phil(j,k) / g
619 if ( ( elvmax(j) .le. pkp1log ) .and.
620 & ( elvmax(j) .ge. pklog ) )
THEN 623 wk(i) = g * elvmax(j) / ( phil(j,k+1) - phil(j,k) )
624 iwklm(i) = max(iwklm(i), k+1 )
630 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
631 vtk(i,k) = vtj(i,k) / prslk(j,k)
632 ro(i,k) = rdi * prsl(j,k) / vtj(i,k)
654 rdz = g / ( phil(j,k+1) - phil(j,k) )
657 bnv2lm(i,k) = (g+g) * rdz * ( vtk(i,k+1)-vtk(i,k) )
658 & / ( vtk(i,k+1)+vtk(i,k) )
659 bnv2lm(i,k) = max( bnv2lm(i,k), bnv2min )
666 delks(i) = 1.0 / (prsi(j,1) - prsi(j,iwklm(i)))
667 delks1(i) = 1.0 / (prsl(j,1) - prsl(j,iwklm(i)))
673 bnv2bar(i) = (prsl(j,1)-prsl(j,2)) * delks1(i) * bnv2lm(i,1)
681 DO ktrial = kmll, 1, -1
683 IF ( ktrial .LT. iwklm(i) .and. kreflm(i) .eq. 0 )
then 698 rdelks = del(j,k) * delks(i)
699 ubar(i) = ubar(i) + rdelks * u1(j,k)
700 vbar(i) = vbar(i) + rdelks * v1(j,k)
701 roll(i) = roll(i) + rdelks * ro(i,k)
702 rdelks = (prsl(j,k)-prsl(j,k+1)) * delks1(i)
703 bnv2bar(i) = bnv2bar(i) + bnv2lm(i,k) * rdelks
715 DO k = iwklm(i), 1, -1
716 phiang = atan2(v1(j,k),u1(j,k))*rad_to_deg
717 ang(i,k) = ( theta(j) - phiang )
718 if ( ang(i,k) .gt. 90. ) ang(i,k) = ang(i,k) - 180.
719 if ( ang(i,k) .lt. -90. ) ang(i,k) = ang(i,k) + 180.
720 ang(i,k) = ang(i,k) * deg_to_rad
729 & max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), minwnd)
731 IF (idxzb(i) .eq. 0 )
then 732 pe(i) = pe(i) + bnv2lm(i,k) *
733 & ( g * elvmax(j) - phil(j,k) ) *
734 & ( phii(j,k+1) - phii(j,k) ) / (g*g)
739 up(i) = uds(i,k) * cos(ang(i,k))
740 ek(i) = 0.5 * up(i) * up(i)
743 IF ( pe(i) .ge. ek(i) )
THEN 745 rdxzb(j) =
real(k,kind=kind_phys)
779 & - sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i))/bnv2bar(i)
794 IF ( idxzb(i) .gt. 0 )
then 795 DO k = idxzb(i), 1, -1
796 IF ( phil(j,idxzb(i)) .gt. phil(j,k) )
then 805 zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) /
806 & ( phil(j,k ) + g * hprime(j) ) )
817 r = (cos(ang(i,k))**2 + gamma(j) * sin(ang(i,k))**2) /
818 & (gamma(j) * cos(ang(i,k))**2 + sin(ang(i,k))**2)
830 dbtmp = 0.25 * cdmb *
831 & max( 2. - 1. / r, 0. ) * sigma(j) *
832 & max(cos(ang(i,k)), gamma(j)*sin(ang(i,k))) *
834 db(i,k) = dbtmp * uds(i,k)
852 ELSEIF ( nmtvr .ne. 14)
then 857 IF ( hprime(i) .GT. hpmin )
then 860 if (ipr .eq. i) npr = npt
863 IF (npt .eq. 0)
RETURN 888 cleff = 0.5e-5 / sqrt(float(imx)/192.0)
894 if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
899 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
900 vtk(i,k) = vtj(i,k) / prslk(j,k)
901 ro(i,k) = rdi * prsl(j,k) / vtj(i,k)
908 ti = 2.0 / (t1(j,k)+t1(j,k+1))
909 tem = ti / (prsl(j,k)-prsl(j,k+1))
910 rdz = g / (phil(j,k+1) - phil(j,k))
911 tem1 = u1(j,k) - u1(j,k+1)
912 tem2 = v1(j,k) - v1(j,k+1)
913 dw2 = tem1*tem1 + tem2*tem2
914 shr2 = max(dw2,dw2min) * rdz * rdz
915 bvf2 = g*(gocp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
916 ri_n(i,k) = max(bvf2/shr2,rimin)
920 bnv2(i,k) = (g+g) * rdz * (vtk(i,k+1)-vtk(i,k))
921 & / (vtk(i,k+1)+vtk(i,k))
922 bnv2(i,k) = max( bnv2(i,k), bnv2min )
948 tem = (prsi(j,1) - prsi(j,k))
949 if (tem .lt. dpmin) iwk(i) = k
959 kref(i) = max(iwk(i), kpbl(j)+1 )
960 delks(i) = 1.0 / (prsi(j,1) - prsi(j,kref(i)))
961 delks1(i) = 1.0 / (prsl(j,1) - prsl(j,kref(i)))
965 kbps = max(kbps, kref(i))
966 kmps = min(kmps, kref(i))
968 bnv2bar(i) = (prsl(j,1)-prsl(j,2)) * delks1(i) * bnv2(i,1)
975 IF (k .LT. kref(i))
THEN 977 rdelks = del(j,k) * delks(i)
978 ubar(i) = ubar(i) + rdelks * u1(j,k)
979 vbar(i) = vbar(i) + rdelks * v1(j,k)
981 roll(i) = roll(i) + rdelks * ro(i,k)
982 rdelks = (prsl(j,k)-prsl(j,k+1)) * delks1(i)
983 bnv2bar(i) = bnv2bar(i) + bnv2(i,k) * rdelks
998 wdir = atan2(ubar(i),vbar(i)) + pi
999 idir = mod(nint(fdir*wdir),mdir) + 1
1001 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
1002 clx(i) = clx4(j,mod(nwd-1,4)+1)
1028 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
1029 uloi(i) = 1.0 / ulow(i)
1035 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*ubar(i)
1036 & + (v1(j,k)+v1(j,k+1))*vbar(i))
1037 velco(i,k) = velco(i,k) * uloi(i)
1070 bnv = sqrt( bnv2bar(i) )
1071 fr = bnv * uloi(i) * min(hprime(j),hpmax)
1073 xn(i) = ubar(i) * uloi(i)
1074 yn(i) = vbar(i) * uloi(i)
1119 efact = (oa(i) + 2.) ** (ceofrc*fr)
1120 efact = min( max(efact,efmin), efmax )
1122 coefm = (1. + clx(i)) ** (oa(i)+1.)
1124 xlinv(i) = coefm * cleff
1126 tem = fr * fr * oc(j)
1127 gfobnv = gmax * tem / ((tem + cg)*bnv)
1129 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
1130 & * ulow(i) * gfobnv * efact
1135 k = max(1, kref(i)-1)
1136 tem = max(velco(i,k)*velco(i,k), 0.1)
1137 scor(i) = bnv2(i,k) / tem
1145 IF (k .LE. kref(i)) taup(i,k) = taub(i)
1159 IF (k .GE. kref(i))
THEN 1160 icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) .LT. ric)
1161 & .OR. (velco(i,k) .LE. 0.0)
1178 IF (k .GE. kref(i))
THEN 1179 IF (.NOT.icrilv(i) .AND. taup(i,k) .GT. 0.0 )
THEN 1180 temv = 1.0 / max(velco(i,k), 0.01)
1182 IF (oa(i).GT.0. .AND. kp1 .lt. kint(i))
THEN 1183 scork = bnv2(i,k) * temv * temv
1184 rscor = min(1.0, scork / scor(i))
1202 brvf = sqrt(bnv2(i,k))
1204 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
1205 & * max(velco(i,k),0.01)
1206 hd = sqrt(taup(i,k) / tem1)
1207 fro = brvf * hd * temv
1220 tem2 = sqrt(ri_n(i,k))
1221 tem = 1. + tem2 * fro
1222 rim = ri_n(i,k) * (1.-fro) / (tem * tem)
1243 IF (rim .LE. ric .AND.
1245 & (oa(i) .LE. 0. .OR. kp1 .ge. kint(i) ))
THEN 1246 temc = 2.0 + 1.0 / tem2
1247 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
1248 taup(i,kp1) = tem1 * hd * hd
1250 taup(i,kp1) = taup(i,k) * rscor
1252 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
1262 IF(lcap .LE. km)
THEN 1263 DO klcap = lcapp1, km+1
1265 sira = prsi(ipt(i),klcap) / prsi(ipt(i),lcap)
1266 taup(i,klcap) = sira * taup(i,lcap)
1275 taud(i,k) = g * (taup(i,k+1) - taup(i,k)) / del(ipt(i),k)
1284 taud(i,klcap) = taud(i,klcap) * factop
1294 IF (k .GT. kref(i) .and. prsi(ipt(i),k) .GE. rlolev)
THEN 1295 IF(taud(i,k).NE.0.)
THEN 1296 tem = deltim * taud(i,k)
1297 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
1315 taud(i,k) = taud(i,k) * dtfac(i)
1316 dtaux = taud(i,k) * xn(i)
1317 dtauy = taud(i,k) * yn(i)
1318 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
1320 if ( k .lt. idxzb(i) .AND. idxzb(i) .ne. 0 )
then 1321 dbim = db(i,k) / (1.+db(i,k)*deltim)
1322 a(j,k) = - dbim * v1(j,k) + a(j,k)
1323 b(j,k) = - dbim * u1(j,k) + b(j,k)
1324 eng1 = eng0*(1.0-dbim*deltim)*(1.0-dbim*deltim)
1328 dusfc(j) = dusfc(j) - dbim * u1(j,k) * del(j,k)
1329 dvsfc(j) = dvsfc(j) - dbim * v1(j,k) * del(j,k)
1332 a(j,k) = dtauy + a(j,k)
1333 b(j,k) = dtaux + b(j,k)
1335 & (u1(j,k)+dtaux*deltim)*(u1(j,k)+dtaux*deltim)
1336 & + (v1(j,k)+dtauy*deltim)*(v1(j,k)+dtauy*deltim))
1337 dusfc(j) = dusfc(j) + dtaux * del(j,k)
1338 dvsfc(j) = dvsfc(j) + dtauy * del(j,k)
1340 c(j,k) = c(j,k) + max(eng0-eng1,0.)/cp/deltim
1352 dusfc(j) = tem * dusfc(j)
1353 dvsfc(j) = tem * dvsfc(j)
subroutine gwdps(IM, IX, IY, KM, A, B, C, U1, V1, T1, Q1, KPBL, PRSI, DEL, PRSL, PRSLK, PHII, PHIL, DELTIM, KDT, HPRIME, OC, OA4, CLX4, THETA, SIGMA, GAMMA, ELVMAX, DUSFC, DVSFC, G, CP, RD, RV, IMX, nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb)