!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