!WRF:model_layer:physics ! ! ! ! module g_module_bl_gwdo contains ! Generated by TAPENADE (INRIA, Tropics team) ! Tapenade 3.10 (r5498) - 20 Jan 2015 09:48 ! ! Differentiation of gwdo in forward (tangent) mode: ! variations of useful results: rublten dusfcg dvsfcg dtauy3d ! rvblten dtaux3d ! with respect to varying inputs: rublten p3di u3d dusfcg z dvsfcg ! dtauy3d rvblten t3d qv3d pi3d v3d dtaux3d ! RW status of diff variables: rublten:in-out p3di:in u3d:in ! dusfcg:in-out z:in dvsfcg:in-out dtauy3d:in-out ! rvblten:in-out t3d:in qv3d:in pi3d:in v3d:in dtaux3d:in-out !------------------------------------------------------------------------------- SUBROUTINE G_GWDO(u3d, u3dd, v3d, v3dd, t3d, t3dd, qv3d, qv3dd, p3d, & & p3di, p3did, pi3d, pi3dd, z, zd, rublten, rubltend, rvblten, rvbltend& & , dtaux3d, dtaux3dd, dtauy3d, dtauy3dd, dusfcg, dusfcgd, dvsfcg, & & dvsfcgd, var2d, oc12d, oa2d1, oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3& & , ol2d4, znu, znw, p_top, cp, g, rd, rv, ep1, pi, dt, dx, & & kpbl2d, itimestep, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, & & kms, kme, its, ite, jts, jte, kts, kte) IMPLICIT NONE ! !------------------------------------------------------------------------------- ! !-- u3d 3d u-velocity interpolated to theta points (m/s) !-- v3d 3d v-velocity interpolated to theta points (m/s) !-- t3d temperature (k) !-- qv3d 3d water vapor mixing ratio (kg/kg) !-- p3d 3d pressure (pa) !-- p3di 3d pressure (pa) at interface level !-- pi3d 3d exner function (dimensionless) !-- rublten u tendency due to pbl parameterization (m/s/s) !-- rvblten v tendency due to pbl parameterization (m/s/s) !-- znu eta values (sigma values) !-- cp heat capacity at constant pressure for dry air (j/kg/k) !-- g acceleration due to gravity (m/s^2) !-- rd gas constant for dry air (j/kg/k) !-- z height above sea level (m) !-- rv gas constant for water vapor (j/kg/k) !-- dt time step (s) !-- dx model grid interval (m) !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless) !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !-- its start index for i in tile !-- ite end index for i in tile !-- jts start index for j in tile !-- jte end index for j in tile !-- kts start index for k in tile !-- kte end index for k in tile ! !------------------------------------------------------------------------------- INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, & & jme, kms, kme, its, ite, jts, jte, kts, kte INTEGER, INTENT(IN) :: itimestep ! REAL, INTENT(IN) :: dt, dx, cp, g, rd, rv, ep1, pi ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3d, p3d, & & pi3d, t3d, z REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3dd, pi3dd& & , t3dd, zd REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3di REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3did ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rublten, & & rvblten REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rubltend& & , rvbltend REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3d, & & dtauy3d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3dd& & , dtauy3dd ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3d, v3d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3dd, v3dd ! INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl2d REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcg, dvsfcg REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcgd, dvsfcgd ! REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: var2d, oc12d, oa2d1, & & oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3, ol2d4 ! REAL, DIMENSION(kms:kme), OPTIONAL, INTENT(IN) :: znu, znw ! REAL, OPTIONAL, INTENT(IN) :: p_top ! !local ! REAL, DIMENSION(its:ite, kts:kte) :: delprsi, pdh REAL, DIMENSION(its:ite, kts:kte) :: delprsid, pdhd REAL, DIMENSION(its:ite, kts:kte + 1) :: pdhi REAL, DIMENSION(its:ite, kts:kte+1) :: pdhid REAL, DIMENSION(its:ite, 4) :: oa4, ol4 INTEGER :: i, j, k, kdt, kpblmax ! DO k=kts,kte IF (znu(k) .GT. 0.6) kpblmax = k + 1 END DO delprsid = 0.0 pdhid = 0.0 pdhd = 0.0 ! DO j=jts,jte DO k=kts,kte+1 DO i=its,ite IF (k .LE. kte) THEN pdhd(i, k) = 0.0 pdh(i, k) = p3d(i, k, j) END IF pdhid(i, k) = p3did(i, k, j) pdhi(i, k) = p3di(i, k, j) END DO END DO ! DO k=kts,kte DO i=its,ite delprsid(i, k) = pdhid(i, k) - pdhid(i, k+1) delprsi(i, k) = pdhi(i, k) - pdhi(i, k+1) END DO END DO DO i=its,ite oa4(i, 1) = oa2d1(i, j) oa4(i, 2) = oa2d2(i, j) oa4(i, 3) = oa2d3(i, j) oa4(i, 4) = oa2d4(i, j) ol4(i, 1) = ol2d1(i, j) ol4(i, 2) = ol2d2(i, j) ol4(i, 3) = ol2d3(i, j) ol4(i, 4) = ol2d4(i, j) END DO CALL GWDO2D_D(dudt=rublten(ims, kms, j), dudtd=rubltend(ims, kms, j)& & , dvdt=rvblten(ims, kms, j), dvdtd=rvbltend(ims, kms, j), & & dtaux2d=dtaux3d(ims, kms, j), dtaux2dd=dtaux3dd(ims, kms, j)& & , dtauy2d=dtauy3d(ims, kms, j), dtauy2dd=dtauy3dd(ims, kms, & & j), u1=u3d(ims, kms, j), u1d=u3dd(ims, kms, j), v1=v3d(ims, & & kms, j), v1d=v3dd(ims, kms, j), t1=t3d(ims, kms, j), t1d=& & t3dd(ims, kms, j), q1=qv3d(ims, kms, j), q1d=qv3dd(ims, kms& & , j), del=delprsi(its, kts), deld=delprsid(its, kts), prsi=& & pdhi(its, kts), prsid=pdhid(its, kts), prsl=pdh(its, kts), & & prsld=pdhd(its, kts), prslk=pi3d(ims, kms, j), prslkd=pi3dd(& & ims, kms, j), zl=z(ims, kms, j), zld=zd(ims, kms, j), rcl=& & 1.0, kpblmax=kpblmax, var=var2d(ims, j), oc1=oc12d(ims, j), & & oa4=oa4, ol4=ol4, dusfc=dusfcg(ims, j), dusfcd=dusfcgd(ims, & & j), dvsfc=dvsfcg(ims, j), dvsfcd=dvsfcgd(ims, j), g=g, cp=cp& & , rd=rd, rv=rv, fv=ep1, pi=pi, dxmeter=dx, deltim=dt, kpbl=& & kpbl2d(ims, j), kdt=itimestep, lat=j, ids=ids, ide=ide, jds=& & jds, jde=jde, kds=kds, kde=kde, ims=ims, ime=ime, jms=jms, & & jme=jme, kms=kms, kme=kme, its=its, ite=ite, jts=jts, jte=& & jte, kts=kts, kte=kte) END DO END SUBROUTINE G_GWDO ! Differentiation of gwdo2d in forward (tangent) mode: ! variations of useful results: dvsfc dvdt dtauy2d dusfc dudt ! dtaux2d ! with respect to varying inputs: v1 dvsfc dvdt dtauy2d prsi ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- SUBROUTINE GWDO2D_D(dudt, dudtd, dvdt, dvdtd, dtaux2d, dtaux2dd, dtauy2d& & , dtauy2dd, u1, u1d, v1, v1d, t1, t1d, q1, q1d, del, deld, prsi, prsid& & , prsl, prsld, prslk, prslkd, zl, zld, rcl, kpblmax, var, oc1, oa4, & & ol4, dusfc, dusfcd, dvsfc, dvsfcd, g, cp, rd, rv, fv, pi, dxmeter, & & deltim, kpbl, kdt, lat, ids, ide, jds, jde, kds, kde, ims, ime, jms, & & jme, kms, kme, its, ite, jts, jte, kts, kte) IMPLICIT NONE !------------------------------------------------------------------------------- INTEGER :: kdt, lat, latd, lond, kpblmax, ids, ide, jds, jde, kds, kde& & , ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte ! REAL :: g, rd, rv, fv, cp, pi, dxmeter, deltim, rcl REAL :: dudt(ims:ime, kms:kme), dvdt(ims:ime, kms:kme), dtaux2d(ims:& & ime, kms:kme), dtauy2d(ims:ime, kms:kme), u1(ims:ime, kms:kme), v1(ims& & :ime, kms:kme), t1(ims:ime, kms:kme), q1(ims:ime, kms:kme), zl(ims:ime& & , kms:kme), prsl(its:ite, kts:kte), prslk(ims:ime, kms:kme) REAL :: dudtd(ims:ime, kms:kme), dvdtd(ims:ime, kms:kme), dtaux2dd(ims& & :ime, kms:kme), dtauy2dd(ims:ime, kms:kme), u1d(ims:ime, kms:kme), v1d& & (ims:ime, kms:kme), t1d(ims:ime, kms:kme), q1d(ims:ime, kms:kme), zld(& & ims:ime, kms:kme), prsld(its:ite, kts:kte), prslkd(ims:ime, kms:kme) REAL :: prsi(its:ite, kts:kte+1), del(its:ite, kts:kte) REAL :: prsid(its:ite, kts:kte+1), deld(its:ite, kts:kte) REAL :: oa4(its:ite, 4), ol4(its:ite, 4) ! INTEGER :: kpbl(ims:ime) REAL :: var(ims:ime), oc1(ims:ime), dusfc(ims:ime), dvsfc(ims:ime) REAL :: dusfcd(ims:ime), dvsfcd(ims:ime) ! ! critical richardson number for wave breaking : ! larger drag with larger value ! REAL, PARAMETER :: ric=0.25 ! REAL, PARAMETER :: dw2min=1. REAL, PARAMETER :: rimin=-100. REAL, PARAMETER :: bnv2min=1.0e-5 REAL, PARAMETER :: efmin=0.0 REAL, PARAMETER :: efmax=10.0 REAL, PARAMETER :: xl=4.0e4 REAL, PARAMETER :: critac=1.0e-5 REAL, PARAMETER :: gmax=1. REAL, PARAMETER :: veleps=1.0 REAL, PARAMETER :: factop=0.5 REAL, PARAMETER :: frc=1.0 REAL, PARAMETER :: ce=0.8 REAL, PARAMETER :: cg=0.5 INTEGER, PARAMETER :: kpblmin=2 ! ! local variables ! INTEGER :: i, k, lcap, lcapp1, nwd, idir, klcap, kp1, ikount, kk ! REAL :: rcs, rclcs, csg, fdir, cleff, cs, rcsks, wdir, ti, rdz, temp, & & tem2, dw2, shr2, bvf2, rdelks, wtkbj, tem, gfobnv, hd, fro, rim, temc& & , tem1, efact, temv, dtaux, dtauy REAL :: rcsksd, tid, rdzd, tem2d, dw2d, shr2d, bvf2d, rdelksd, wtkbjd& & , temd, gfobnvd, hdd, temcd, tem1d, efactd, dtauxd, dtauyd ! LOGICAL :: ldrag(its:ite), icrilv(its:ite), flag(its:ite), kloop1(its:& & ite) ! REAL :: taub(its:ite), taup(its:ite, kts:kte+1), xn(its:ite), yn(its:& & ite), ubar(its:ite), vbar(its:ite), fr(its:ite), ulow(its:ite), rulow(& & its:ite), bnv(its:ite), oa(its:ite), ol(its:ite), roll(its:ite), dtfac& & (its:ite), brvf(its:ite), xlinv(its:ite), delks(its:ite), delks1(its:& & ite), bnv2(its:ite, kts:kte), usqj(its:ite, kts:kte), taud(its:ite, & & kts:kte), ro(its:ite, kts:kte), vtk(its:ite, kts:kte), vtj(its:ite, & & kts:kte), zlowtop(its:ite), velco(its:ite, kts:kte-1), coefm(its:ite) REAL :: taubd(its:ite), taupd(its:ite, kts:kte+1), xnd(its:ite), ynd(& & its:ite), ubard(its:ite), vbard(its:ite), frd(its:ite), ulowd(its:ite)& & , rulowd(its:ite), bnvd(its:ite), rolld(its:ite), dtfacd(its:ite), & & brvfd(its:ite), delksd(its:ite), delks1d(its:ite), bnv2d(its:ite, kts:& & kte), usqjd(its:ite, kts:kte), taudd(its:ite, kts:kte), rod(its:ite, & & kts:kte), vtkd(its:ite, kts:kte), vtjd(its:ite, kts:kte), velcod(its:& & ite, kts:kte-1) ! INTEGER :: kbl(its:ite), klowtop(its:ite) ! LOGICAL :: iope INTEGER, PARAMETER :: mdir=8 INTEGER :: nwdir(mdir) ! ! variables for flow-blocking drag ! REAL, PARAMETER :: frmax=10. REAL, PARAMETER :: olmin=1.0e-5 REAL, PARAMETER :: odmin=0.1 REAL, PARAMETER :: odmax=10. REAL, PARAMETER :: erad=6371.315e+3 INTEGER :: komax(its:ite) INTEGER :: kblk REAL :: cd REAL :: zblk, tautem REAL :: zblkd, tautemd REAL :: pe, ke REAL :: delx, dely, dxy4(4), dxy4p(4) REAL :: dxy(its:ite), dxyp(its:ite) REAL :: ol4p(4), olp(its:ite), od(its:ite) REAL :: taufb(its:ite, kts:kte+1) REAL :: taufbd(its:ite, kts:kte+1) REAL :: result1 REAL :: result1d REAL :: arg1 REAL :: arg1d REAL :: pwx1 REAL :: pwy1 REAL :: pwy1d REAL :: arg10 REAL :: pwy10 REAL :: max2d REAL :: x4 REAL :: x3 REAL :: x2 REAL :: x2d INTEGER :: x1 REAL :: x4d REAL :: max3 REAL :: max2 REAL :: max1 REAL :: x3d REAL :: y1 REAL :: y1d DATA nwdir /6, 7, 5, 8, 2, 3, 1, 4/ ! !---- constants ! rcs = SQRT(rcl) result1 = SQRT(rcl) cs = 1./result1 csg = cs*g lcap = kte lcapp1 = lcap + 1 fdir = mdir/(2.0*pi) ! !--- calculate length of grid for flow-blocking drag ! delx = dxmeter dely = dxmeter dxy4(1) = delx dxy4(2) = dely arg1 = delx*delx + dely*dely dxy4(3) = SQRT(arg1) dxy4(4) = dxy4(3) dxy4p(1) = dxy4(2) dxy4p(2) = dxy4(1) dxy4p(3) = dxy4(4) dxy4p(4) = dxy4(3) ! ! !-----initialize arrays ! dtaux = 0.0 dtauy = 0.0 DO i=its,ite klowtop(i) = 0 kbl(i) = 0 END DO ! DO i=its,ite xn(i) = 0.0 yn(i) = 0.0 ubar(i) = 0.0 vbar(i) = 0.0 roll(i) = 0.0 taub(i) = 0.0 taup(i, 1) = 0.0 oa(i) = 0.0 ol(i) = 0.0 ulow(i) = 0.0 dtfac(i) = 1.0 ldrag(i) = .false. icrilv(i) = .false. flag(i) = .true. END DO ! DO k=kts,kte DO i=its,ite usqj(i, k) = 0.0 bnv2(i, k) = 0.0 vtj(i, k) = 0.0 vtk(i, k) = 0.0 taup(i, k) = 0.0 taud(i, k) = 0.0 dtaux2dd(i, k) = 0.0 dtaux2d(i, k) = 0.0 dtauy2dd(i, k) = 0.0 dtauy2d(i, k) = 0.0 END DO END DO ! DO i=its,ite taup(i, kte+1) = 0.0 xlinv(i) = 1.0/xl END DO ! ! initialize array for flow-blocking drag ! taufb(its:ite, kts:kte+1) = 0.0 komax(its:ite) = 0 vtjd = 0.0 vtkd = 0.0 rod = 0.0 ! DO k=kts,kte DO i=its,ite vtjd(i, k) = t1d(i, k)*(1.+fv*q1(i, k)) + t1(i, k)*fv*q1d(i, k) vtj(i, k) = t1(i, k)*(1.+fv*q1(i, k)) vtkd(i, k) = (vtjd(i, k)*prslk(i, k)-vtj(i, k)*prslkd(i, k))/prslk& & (i, k)**2 vtk(i, k) = vtj(i, k)/prslk(i, k) ! density kg/m**3 rod(i, k) = (prsld(i, k)*vtj(i, k)/rd-prsl(i, k)*vtjd(i, k)/rd)/& & vtj(i, k)**2 ro(i, k) = 1./rd*prsl(i, k)/vtj(i, k) END DO END DO ! ! determine reference level: maximum of 2*var and pbl heights ! DO i=its,ite zlowtop(i) = 2.*var(i) END DO ! DO i=its,ite kloop1(i) = .true. END DO ! DO k=kts+1,kte DO i=its,ite IF (kloop1(i) .AND. zl(i, k) - zl(i, 1) .GE. zlowtop(i)) THEN klowtop(i) = k + 1 kloop1(i) = .false. END IF END DO END DO ! DO i=its,ite IF (kpbl(i) .LT. klowtop(i)) THEN kbl(i) = klowtop(i) ELSE kbl(i) = kpbl(i) END IF IF (kbl(i) .GT. kpblmax) THEN x1 = kpblmax ELSE x1 = kbl(i) END IF IF (x1 .LT. kpblmin) THEN kbl(i) = kpblmin ELSE kbl(i) = x1 END IF END DO ! ! determine the level of maximum orographic height ! komax(:) = kbl(:) delksd = 0.0 delks1d = 0.0 ! DO i=its,ite delksd(i) = -((prsid(i, 1)-prsid(i, kbl(i)))/(prsi(i, 1)-prsi(i, kbl& & (i)))**2) delks(i) = 1.0/(prsi(i, 1)-prsi(i, kbl(i))) delks1d(i) = -((prsld(i, 1)-prsld(i, kbl(i)))/(prsl(i, 1)-prsl(i, & & kbl(i)))**2) delks1(i) = 1.0/(prsl(i, 1)-prsl(i, kbl(i))) END DO rolld = 0.0 vbard = 0.0 ubard = 0.0 ! ! compute low level averages within pbl ! DO k=kts,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN rcsksd = rcs*(deld(i, k)*delks(i)+del(i, k)*delksd(i)) rcsks = rcs*del(i, k)*delks(i) rdelksd = deld(i, k)*delks(i) + del(i, k)*delksd(i) rdelks = del(i, k)*delks(i) ! pbl u mean ubard(i) = ubard(i) + rcsksd*u1(i, k) + rcsks*u1d(i, k) ubar(i) = ubar(i) + rcsks*u1(i, k) ! pbl v mean vbard(i) = vbard(i) + rcsksd*v1(i, k) + rcsks*v1d(i, k) vbar(i) = vbar(i) + rcsks*v1(i, k) ! ro mean rolld(i) = rolld(i) + rdelksd*ro(i, k) + rdelks*rod(i, k) roll(i) = roll(i) + rdelks*ro(i, k) END IF END DO END DO ! ! figure out low-level horizontal wind direction ! ! nwd 1 2 3 4 5 6 7 8 ! wd w s sw nw e n ne se ! DO i=its,ite wdir = ATAN2(ubar(i), vbar(i)) + pi idir = MOD(NINT(fdir*wdir), mdir) + 1 nwd = nwdir(idir) oa(i) = (1-2*INT((nwd-1)/4))*oa4(i, MOD(nwd-1, 4)+1) ol(i) = ol4(i, MOD(nwd-1, 4)+1) ! !----- compute orographic width along (ol) and perpendicular (olp) !----- the direction of wind ! ol4p(1) = ol4(i, 2) ol4p(2) = ol4(i, 1) ol4p(3) = ol4(i, 4) ol4p(4) = ol4(i, 3) olp(i) = ol4p(MOD(nwd-1, 4)+1) IF (ol(i) .LT. olmin) THEN max1 = olmin ELSE max1 = ol(i) END IF ! !----- compute orographic direction (horizontal orographic aspect ratio) ! od(i) = olp(i)/max1 IF (od(i) .GT. odmax) THEN od(i) = odmax ELSE od(i) = od(i) END IF IF (od(i) .LT. odmin) THEN od(i) = odmin ELSE od(i) = od(i) END IF ! !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions ! dxy(i) = dxy4(MOD(nwd-1, 4)+1) dxyp(i) = dxy4p(MOD(nwd-1, 4)+1) END DO bnv2d = 0.0 usqjd = 0.0 ! !--- saving richardson number in usqj for migwdi ! DO k=kts,kte-1 DO i=its,ite tid = -(2.0*(t1d(i, k)+t1d(i, k+1))/(t1(i, k)+t1(i, k+1))**2) ti = 2.0/(t1(i, k)+t1(i, k+1)) rdzd = -((zld(i, k+1)-zld(i, k))/(zl(i, k+1)-zl(i, k))**2) rdz = 1./(zl(i, k+1)-zl(i, k)) tem1d = u1d(i, k) - u1d(i, k+1) tem1 = u1(i, k) - u1(i, k+1) tem2d = v1d(i, k) - v1d(i, k+1) tem2 = v1(i, k) - v1(i, k+1) dw2d = rcl*(tem1d*tem1+tem1*tem1d+tem2d*tem2+tem2*tem2d) dw2 = rcl*(tem1*tem1+tem2*tem2) IF (dw2 .LT. dw2min) THEN max2 = dw2min max2d = 0.0 ELSE max2d = dw2d max2 = dw2 END IF shr2d = (max2d*rdz+max2*rdzd)*rdz + max2*rdz*rdzd shr2 = max2*rdz*rdz bvf2d = g*((rdzd*(vtj(i, k+1)-vtj(i, k))+rdz*(vtjd(i, k+1)-vtjd(i& & , k)))*ti+(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*tid) bvf2 = g*(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*ti IF (bvf2/shr2 .LT. rimin) THEN usqjd(i, k) = 0.0 usqj(i, k) = rimin ELSE usqjd(i, k) = (bvf2d*shr2-bvf2*shr2d)/shr2**2 usqj(i, k) = bvf2/shr2 END IF bnv2d(i, k) = (2.0*g*(rdzd*(vtk(i, k+1)-vtk(i, k))+rdz*(vtkd(i, k+& & 1)-vtkd(i, k)))*(vtk(i, k+1)+vtk(i, k))-2.0*g*rdz*(vtk(i, k+1)-& & vtk(i, k))*(vtkd(i, k+1)+vtkd(i, k)))/(vtk(i, k+1)+vtk(i, k))**2 bnv2(i, k) = 2.0*g*rdz*(vtk(i, k+1)-vtk(i, k))/(vtk(i, k+1)+vtk(i& & , k)) IF (bnv2(i, k) .LT. bnv2min) THEN bnv2d(i, k) = 0.0 bnv2(i, k) = bnv2min ELSE bnv2(i, k) = bnv2(i, k) END IF END DO END DO rulowd = 0.0 ulowd = 0.0 ! !----compute the "low level" or 1/3 wind magnitude (m/s) ! DO i=its,ite arg1d = ubard(i)*ubar(i) + ubar(i)*ubard(i) + vbard(i)*vbar(i) + & & vbar(i)*vbard(i) arg1 = ubar(i)*ubar(i) + vbar(i)*vbar(i) IF (arg1 .EQ. 0.0) THEN x2d = 0.0 ELSE x2d = arg1d/(2.0*SQRT(arg1)) END IF x2 = SQRT(arg1) IF (x2 .LT. 1.0) THEN ulowd(i) = 0.0 ulow(i) = 1.0 ELSE ulowd(i) = x2d ulow(i) = x2 END IF rulowd(i) = -(ulowd(i)/ulow(i)**2) rulow(i) = 1./ulow(i) END DO velcod = 0.0 ! DO k=kts,kte-1 DO i=its,ite velcod(i, k) = 0.5*rcs*((u1d(i, k)+u1d(i, k+1))*ubar(i)+(u1(i, k)+& & u1(i, k+1))*ubard(i)+(v1d(i, k)+v1d(i, k+1))*vbar(i)+(v1(i, k)+& & v1(i, k+1))*vbard(i)) velco(i, k) = 0.5*rcs*((u1(i, k)+u1(i, k+1))*ubar(i)+(v1(i, k)+v1(& & i, k+1))*vbar(i)) velcod(i, k) = velcod(i, k)*rulow(i) + velco(i, k)*rulowd(i) velco(i, k) = velco(i, k)*rulow(i) IF (velco(i, k) .LT. veleps .AND. velco(i, k) .GT. 0.) THEN velcod(i, k) = 0.0 velco(i, k) = veleps END IF END DO END DO ! ! no drag when critical level in the base layer ! DO i=its,ite ldrag(i) = velco(i, 1) .LE. 0. END DO ! ! no drag when velco.lt.0 ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. velco(i, k) .LE. 0. END DO END DO ! ! no drag when bnv2.lt.0 ! DO k=kts,kpblmax DO i=its,ite IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. bnv2(i, k) .LT. 0. END DO END DO ! !-----the low level weighted average ri is stored in usqj(1,1; im) !-----the low level weighted average n**2 is stored in bnv2(1,1; im) !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of / ! DO i=its,ite wtkbjd = (prsld(i, 1)-prsld(i, 2))*delks1(i) + (prsl(i, 1)-prsl(i, 2& & ))*delks1d(i) wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i) bnv2d(i, 1) = wtkbjd*bnv2(i, 1) + wtkbj*bnv2d(i, 1) bnv2(i, 1) = wtkbj*bnv2(i, 1) usqjd(i, 1) = wtkbjd*usqj(i, 1) + wtkbj*usqjd(i, 1) usqj(i, 1) = wtkbj*usqj(i, 1) END DO ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN rdelksd = (prsld(i, k)-prsld(i, k+1))*delks1(i) + (prsl(i, k)-& & prsl(i, k+1))*delks1d(i) rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i) bnv2d(i, 1) = bnv2d(i, 1) + bnv2d(i, k)*rdelks + bnv2(i, k)*& & rdelksd bnv2(i, 1) = bnv2(i, 1) + bnv2(i, k)*rdelks usqjd(i, 1) = usqjd(i, 1) + usqjd(i, k)*rdelks + usqj(i, k)*& & rdelksd usqj(i, 1) = usqj(i, 1) + usqj(i, k)*rdelks END IF END DO END DO ! DO i=its,ite ldrag(i) = ldrag(i) .OR. bnv2(i, 1) .LE. 0.0 ldrag(i) = ldrag(i) .OR. ulow(i) .EQ. 1.0 ldrag(i) = ldrag(i) .OR. var(i) .LE. 0.0 END DO ! ! set all ri low level values to the low level value ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN usqjd(i, k) = usqjd(i, 1) usqj(i, k) = usqj(i, 1) END IF END DO END DO bnvd = 0.0 xnd = 0.0 ynd = 0.0 frd = 0.0 ! DO i=its,ite IF (.NOT.ldrag(i)) THEN IF (bnv2(i, 1) .EQ. 0.0) THEN bnvd(i) = 0.0 ELSE bnvd(i) = bnv2d(i, 1)/(2.0*SQRT(bnv2(i, 1))) END IF bnv(i) = SQRT(bnv2(i, 1)) frd(i) = 2.*var(i)*od(i)*(bnvd(i)*rulow(i)+bnv(i)*rulowd(i)) fr(i) = bnv(i)*rulow(i)*2.*var(i)*od(i) IF (fr(i) .GT. frmax) THEN frd(i) = 0.0 fr(i) = frmax ELSE fr(i) = fr(i) END IF xnd(i) = ubard(i)*rulow(i) + ubar(i)*rulowd(i) xn(i) = ubar(i)*rulow(i) ynd(i) = vbard(i)*rulow(i) + vbar(i)*rulowd(i) yn(i) = vbar(i)*rulow(i) END IF END DO taubd = 0.0 ! ! compute the base level stress and store it in taub ! calculate enhancement factor, number of mountains & aspect ! ratio const. use simplified relationship between standard ! deviation & critical hgt ! DO i=its,ite IF (.NOT.ldrag(i)) THEN pwx1 = oa(i) + 2. pwy1d = ce*frd(i)/frc pwy1 = ce*fr(i)/frc IF (pwx1 .GT. 0.0) THEN efactd = LOG(pwx1)*pwx1**pwy1*pwy1d ELSE efactd = 0.0 END IF efact = pwx1**pwy1 IF (efact .LT. efmin) THEN x3 = efmin x3d = 0.0 ELSE x3d = efactd x3 = efact END IF IF (x3 .GT. efmax) THEN efact = efmax efactd = 0.0 ELSE efactd = x3d efact = x3 END IF !!!!!!! cleff (effective grid length) is highly tunable parameter !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag arg10 = dxy(i)**2. + dxyp(i)**2. cleff = SQRT(arg10) IF (dxmeter .LT. cleff) THEN max3 = cleff ELSE max3 = dxmeter END IF cleff = 3.*max3 pwx1 = 1. + ol(i) pwy10 = oa(i) + 1. coefm(i) = pwx1**pwy10 xlinv(i) = coefm(i)/cleff temd = oc1(i)*(frd(i)*fr(i)+fr(i)*frd(i)) tem = fr(i)*fr(i)*oc1(i) gfobnvd = (gmax*temd*(tem+cg)*bnv(i)-gmax*tem*(temd*bnv(i)+(tem+cg& & )*bnvd(i)))/((tem+cg)*bnv(i))**2 gfobnv = gmax*tem/((tem+cg)*bnv(i)) taubd(i) = xlinv(i)*(((rolld(i)*ulow(i)+roll(i)*ulowd(i))*gfobnv*& & efact+roll(i)*ulow(i)*(gfobnvd*efact+gfobnv*efactd))*ulow(i)**2+& & roll(i)*ulow(i)*gfobnv*efact*(ulowd(i)*ulow(i)+ulow(i)*ulowd(i))& & ) taub(i) = xlinv(i)*roll(i)*ulow(i)*ulow(i)*ulow(i)*gfobnv*efact ELSE taubd(i) = 0.0 taub(i) = 0.0 xnd(i) = 0.0 xn(i) = 0.0 ynd(i) = 0.0 yn(i) = 0.0 END IF END DO taupd = 0.0 ! ! now compute vertical structure of the stress. ! DO k=kts,kpblmax DO i=its,ite IF (k .LE. kbl(i)) THEN taupd(i, k) = taubd(i) taup(i, k) = taub(i) END IF END DO END DO brvfd = 0.0 ! ! vertical level k loop! DO k=kpblmin,kte-1 kp1 = k + 1 DO i=its,ite ! ! unstablelayer if ri < ric ! unstable layer if upper air vel comp along surf vel <=0 (crit lay) ! at (u-c)=0. crit layer exists and bit vector should be set (.le.) ! IF (k .GE. kbl(i)) THEN icrilv(i) = (icrilv(i) .OR. usqj(i, k) .LT. ric) .OR. velco(i, k& & ) .LE. 0.0 IF (bnv2(i, k) .LT. bnv2min) THEN brvfd(i) = 0.0 brvf(i) = bnv2min ELSE brvfd(i) = bnv2d(i, k) brvf(i) = bnv2(i, k) END IF ! brunt-vaisala frequency IF (brvf(i) .EQ. 0.0) THEN brvfd(i) = 0.0 ELSE brvfd(i) = brvfd(i)/(2.0*SQRT(brvf(i))) END IF brvf(i) = SQRT(brvf(i)) END IF END DO ! DO i=its,ite IF (k .GE. kbl(i) .AND. (.NOT.ldrag(i))) THEN IF (.NOT.icrilv(i) .AND. taup(i, k) .GT. 0.0) THEN temv = 1.0/velco(i, k) tem1d = coefm(i)*0.5*((rod(i, kp1)+rod(i, k))*brvf(i)*velco(i& & , k)+(ro(i, kp1)+ro(i, k))*(brvfd(i)*velco(i, k)+brvf(i)*& & velcod(i, k)))/dxy(i) tem1 = coefm(i)/dxy(i)*(ro(i, kp1)+ro(i, k))*brvf(i)*velco(i, & & k)*0.5 hd = SQRT(taup(i, k)/tem1) fro = brvf(i)*hd*temv ! ! rim is the minimum-richardson number by shutts (1985) ! IF (usqj(i, k) .EQ. 0.0) THEN tem2d = 0.0 ELSE tem2d = usqjd(i, k)/(2.0*SQRT(usqj(i, k))) END IF tem2 = SQRT(usqj(i, k)) tem = 1. + tem2*fro rim = usqj(i, k)*(1.-fro)/(tem*tem) ! ! check stability to employ the 'saturation hypothesis' ! of lindzen (1981) except at tropospheric downstream regions ! IF (rim .LE. ric) THEN ! saturation hypothesis! IF (oa(i) .LE. 0. .OR. kp1 .GE. kpblmin) THEN temcd = -(tem2d/tem2**2) temc = 2.0 + 1.0/tem2 IF (temc .EQ. 0.0) THEN result1d = 0.0 ELSE result1d = temcd/(2.0*SQRT(temc)) END IF result1 = SQRT(temc) hdd = ((velcod(i, k)*(2.*result1-temc)+velco(i, k)*(2.*& & result1d-temcd))*brvf(i)-velco(i, k)*(2.*result1-temc)*& & brvfd(i))/brvf(i)**2 hd = velco(i, k)*(2.*result1-temc)/brvf(i) taupd(i, kp1) = (tem1d*hd+tem1*hdd)*hd + tem1*hd*hdd taup(i, kp1) = tem1*hd*hd END IF ELSE ! no wavebreaking! taupd(i, kp1) = taupd(i, k) taup(i, kp1) = taup(i, k) END IF END IF END IF END DO END DO ! IF (lcap .LT. kte) THEN DO klcap=lcapp1,kte DO i=its,ite taupd(i, klcap) = (prsid(i, klcap)*prsi(i, lcap)-prsi(i, klcap)*& & prsid(i, lcap))*taup(i, lcap)/prsi(i, lcap)**2 + prsi(i, klcap& & )*taupd(i, lcap)/prsi(i, lcap) taup(i, klcap) = prsi(i, klcap)/prsi(i, lcap)*taup(i, lcap) END DO END DO zblkd = 0.0 taufbd = 0.0 ELSE zblkd = 0.0 taufbd = 0.0 END IF DO i=its,ite IF (.NOT.ldrag(i)) THEN ! !------- determine the height of flow-blocking layer ! kblk = 0 pe = 0.0 DO k=kte,kpblmin,-1 IF (kblk .EQ. 0 .AND. k .LE. komax(i)) THEN pe = pe + bnv2(i, k)*(zl(i, komax(i))-zl(i, k))*del(i, k)/g/ro& & (i, k) ke = 0.5*((rcs*u1(i, k))**2.+(rcs*v1(i, k))**2.) ! !---------- apply flow-blocking drag when pe >= ke ! IF (pe .GE. ke) THEN kblk = k IF (kblk .GT. kbl(i)) THEN kblk = kbl(i) ELSE kblk = kblk END IF zblkd = zld(i, kblk) - zld(i, kts) zblk = zl(i, kblk) - zl(i, kts) END IF END IF END DO IF (kblk .NE. 0) THEN ! !--------- compute flow-blocking stress ! IF (2.0 - 1.0/od(i) .LT. 0.0) THEN cd = 0.0 ELSE cd = 2.0 - 1.0/od(i) END IF taufbd(i, kts) = cd*dxyp(i)*olp(i)*(0.5*coefm(i)*rolld(i)*zblk*& & ulow(i)**2/dxy(i)**2+0.5*roll(i)*coefm(i)*(zblkd*ulow(i)**2+& & zblk*2*ulow(i)*ulowd(i))/dxy(i)**2) taufb(i, kts) = 0.5*roll(i)*coefm(i)/dxy(i)**2*cd*dxyp(i)*olp(i)& & *zblk*ulow(i)**2 tautemd = taufbd(i, kts)/FLOAT(kblk-kts) tautem = taufb(i, kts)/FLOAT(kblk-kts) DO k=kts+1,kblk taufbd(i, k) = taufbd(i, k-1) - tautemd taufb(i, k) = taufb(i, k-1) - tautem END DO ! !----------sum orographic GW stress and flow-blocking stress ! taupd(i, :) = taupd(i, :) + taufbd(i, :) taup(i, :) = taup(i, :) + taufb(i, :) END IF END IF END DO taudd = 0.0 ! ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy ! DO k=kts,kte DO i=its,ite taudd(i, k) = (csg*(taupd(i, k+1)-taupd(i, k))*del(i, k)-(taup(i, & & k+1)-taup(i, k))*csg*deld(i, k))/del(i, k)**2 taud(i, k) = 1.*(taup(i, k+1)-taup(i, k))*csg/del(i, k) END DO END DO ! ! limit de-acceleration (momentum deposition ) at top to 1/2 value ! the idea is some stuff must go out the 'top' ! DO klcap=lcap,kte DO i=its,ite taudd(i, klcap) = factop*taudd(i, klcap) taud(i, klcap) = taud(i, klcap)*factop END DO END DO dtfacd = 0.0 ! ! if the gravity wave drag would force a critical line ! in the lower ksmm1 layers during the next deltim timestep, ! then only apply drag until that critical line is reached. ! DO k=kts,kpblmax-1 DO i=its,ite IF (k .LE. kbl(i)) THEN IF (taud(i, k) .NE. 0.) THEN x4d = (velcod(i, k)*deltim*rcs*taud(i, k)-velco(i, k)*deltim*& & rcs*taudd(i, k))/(deltim*rcs*taud(i, k))**2 x4 = velco(i, k)/(deltim*rcs*taud(i, k)) IF (x4 .GE. 0.) THEN y1d = x4d y1 = x4 ELSE y1d = -x4d y1 = -x4 END IF IF (dtfac(i) .GT. y1) THEN dtfacd(i) = y1d dtfac(i) = y1 ELSE dtfac(i) = dtfac(i) END IF END IF END IF END DO END DO ! DO i=its,ite dusfcd(i) = 0.0 dusfc(i) = 0. dvsfcd(i) = 0.0 dvsfc(i) = 0. END DO ! DO k=kts,kte DO i=its,ite taudd(i, k) = taudd(i, k)*dtfac(i) + taud(i, k)*dtfacd(i) taud(i, k) = taud(i, k)*dtfac(i) dtauxd = taudd(i, k)*xn(i) + taud(i, k)*xnd(i) dtaux = taud(i, k)*xn(i) dtauyd = taudd(i, k)*yn(i) + taud(i, k)*ynd(i) dtauy = taud(i, k)*yn(i) dtaux2dd(i, k) = dtauxd dtaux2d(i, k) = dtaux dtauy2dd(i, k) = dtauyd dtauy2d(i, k) = dtauy dudtd(i, k) = dtauxd + dudtd(i, k) dudt(i, k) = dtaux + dudt(i, k) dvdtd(i, k) = dtauyd + dvdtd(i, k) dvdt(i, k) = dtauy + dvdt(i, k) dusfcd(i) = dusfcd(i) + dtauxd*del(i, k) + dtaux*deld(i, k) dusfc(i) = dusfc(i) + dtaux*del(i, k) dvsfcd(i) = dvsfcd(i) + dtauyd*del(i, k) + dtauy*deld(i, k) dvsfc(i) = dvsfc(i) + dtauy*del(i, k) END DO END DO ! DO i=its,ite dusfcd(i) = -(rcs*dusfcd(i)/g) dusfc(i) = (-(1./g*rcs))*dusfc(i) dvsfcd(i) = -(rcs*dvsfcd(i)/g) dvsfc(i) = (-(1./g*rcs))*dvsfc(i) END DO ! RETURN END SUBROUTINE GWDO2D_D end module g_module_bl_gwdo