!WRF:model_layer:physics
!
!
!
!
module a_module_bl_gwdo
contains
!
!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.10 (r5498) - 20 Jan 2015 09:48
!
!  Differentiation of gwdo in reverse (adjoint) mode:
!   gradient     of useful results: rublten p3di u3d dusfcg z dvsfcg
!                dtauy3d rvblten t3d qv3d pi3d v3d 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:incr u3d:incr
!                dusfcg:in-out z:incr dvsfcg:in-out dtauy3d:in-out
!                rvblten:in-out t3d:incr qv3d:incr pi3d:incr v3d:incr
!                dtaux3d:in-out
!-------------------------------------------------------------------------------
SUBROUTINE GWDO_B(u3d, u3db, v3d, v3db, t3d, t3db, qv3d, qv3db, p3d, &
& p3di, p3dib, pi3d, pi3db, z, zb, rublten, rubltenb, rvblten, rvbltenb&
& , dtaux3d, dtaux3db, dtauy3d, dtauy3db, dusfcg, dusfcgb, dvsfcg, &
& dvsfcgb, 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) :: qv3db, pi3db, t3db, zb
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3di
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: p3dib
!
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rublten, &
& rvblten
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rubltenb
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3d, &
& dtauy3d
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3db
!
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3d, v3d
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: u3db, v3db
!
  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) :: dusfcgb
!
  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) :: delprsib, pdhb
  REAL, DIMENSION(its:ite, kts:kte + 1) :: pdhi
  REAL, DIMENSION(its:ite, kts:kte+1) :: pdhib
  REAL, DIMENSION(its:ite, 4) :: oa4, ol4
  INTEGER :: i, j, k, kdt, kpblmax
  INTEGER :: branch
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rvbltenb
  REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dvsfcgb
  REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtauy3db
!
  DO k=kts,kte
    IF (znu(k) .GT. 0.6) kpblmax = k + 1
  END DO
!
  DO j=jts,jte
     DO k=kts,kte+1
        DO i=its,ite
          IF (k .LE. kte) THEN
            CALL PUSHREAL8(pdh(i, k))
            pdh(i, k) = p3d(i, k, j)
            CALL PUSHCONTROL1B(0)
          ELSE
            CALL PUSHCONTROL1B(1)
          END IF
          CALL PUSHREAL8(pdhi(i, k))
          pdhi(i, k) = p3di(i, k, j)
        END DO
     END DO
!
    DO k=kts,kte
      DO i=its,ite
        CALL PUSHREAL8(delprsi(i, k))
        delprsi(i, k) = pdhi(i, k) - pdhi(i, k+1)
      END DO
    END DO
    DO i=its,ite
      CALL PUSHREAL8(oa4(i, 1))
      oa4(i, 1) = oa2d1(i, j)
      CALL PUSHREAL8(oa4(i, 2))
      oa4(i, 2) = oa2d2(i, j)
      CALL PUSHREAL8(oa4(i, 3))
      oa4(i, 3) = oa2d3(i, j)
      CALL PUSHREAL8(oa4(i, 4))
      oa4(i, 4) = oa2d4(i, j)
      CALL PUSHREAL8(ol4(i, 1))
      ol4(i, 1) = ol2d1(i, j)
      CALL PUSHREAL8(ol4(i, 2))
      ol4(i, 2) = ol2d2(i, j)
      CALL PUSHREAL8(ol4(i, 3))
      ol4(i, 3) = ol2d3(i, j)
      CALL PUSHREAL8(ol4(i, 4))
      ol4(i, 4) = ol2d4(i, j)
    END DO
  END DO
  delprsib = 0.0
  pdhib = 0.0
  pdhb = 0.0
  DO j=jte,jts,-1
    CALL GWDO2D_B(rublten(ims, kms, j), rubltenb(ims, kms, j), rvblten(&
&           ims, kms, j), rvbltenb(ims, kms, j), dtaux3d(ims, kms, j), &
&           dtaux3db(ims, kms, j), dtauy3d(ims, kms, j), dtauy3db(ims, &
&           kms, j), u3d(ims, kms, j), u3db(ims, kms, j), v3d(ims, kms, &
&           j), v3db(ims, kms, j), t3d(ims, kms, j), t3db(ims, kms, j), &
&           qv3d(ims, kms, j), qv3db(ims, kms, j), delprsi(its, kts), &
&           delprsib(its, kts), pdhi(its, kts), pdhib(its, kts), pdh(its&
&           , kts), pdhb(its, kts), pi3d(ims, kms, j), pi3db(ims, kms, j&
&           ), z(ims, kms, j), zb(ims, kms, j), rcl=1.0, kpblmax=kpblmax&
&           , var=var2d(ims, j), oc1=oc12d(ims, j), oa4=oa4, ol4=ol4, &
&           dusfc=dusfcg(ims, j), dusfcb=dusfcgb(ims, j), dvsfc=dvsfcg(&
&           ims, j), dvsfcb=dvsfcgb(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)
    DO i=ite,its,-1
      CALL POPREAL8(ol4(i, 4))
      CALL POPREAL8(ol4(i, 3))
      CALL POPREAL8(ol4(i, 2))
      CALL POPREAL8(ol4(i, 1))
      CALL POPREAL8(oa4(i, 4))
      CALL POPREAL8(oa4(i, 3))
      CALL POPREAL8(oa4(i, 2))
      CALL POPREAL8(oa4(i, 1))
    END DO
    DO k=kte,kts,-1
      DO i=ite,its,-1
        CALL POPREAL8(delprsi(i, k))
        pdhib(i, k) = pdhib(i, k) + delprsib(i, k)
        pdhib(i, k+1) = pdhib(i, k+1) - delprsib(i, k)
        delprsib(i, k) = 0.0
      END DO
    END DO
    DO k=kte+1,kts,-1
       DO i=ite,its,-1
          CALL POPREAL8(pdhi(i, k))
          p3dib(i, k, j) = p3dib(i, k, j) + pdhib(i, k)
          pdhib(i, k) = 0.0
          CALL POPCONTROL1B(branch)
          IF (branch .EQ. 0) THEN
            CALL POPREAL8(pdh(i, k))
            pdhb(i, k) = 0.0
          END IF
       END DO
    END DO
  END DO
END SUBROUTINE GWDO_B

!  Differentiation of gwdo2d in reverse (adjoint) mode:
!   gradient     of useful results: v1 dvsfc dvdt dtauy2d prsi
!                prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk
!   with respect to varying inputs: v1 dvsfc dvdt dtauy2d prsi
!                prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
SUBROUTINE GWDO2D_B(dudt, dudtb, dvdt, dvdtb, dtaux2d, dtaux2db, dtauy2d&
& , dtauy2db, u1, u1b, v1, v1b, t1, t1b, q1, q1b, del, delb, prsi, prsib&
& , prsl, prslb, prslk, prslkb, zl, zlb, rcl, kpblmax, var, oc1, oa4, &
& ol4, dusfc, dusfcb, dvsfc, dvsfcb, 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 :: dudtb(ims:ime, kms:kme), dvdtb(ims:ime, kms:kme), dtaux2db(ims&
& :ime, kms:kme), dtauy2db(ims:ime, kms:kme), u1b(ims:ime, kms:kme), v1b&
& (ims:ime, kms:kme), t1b(ims:ime, kms:kme), q1b(ims:ime, kms:kme), zlb(&
& ims:ime, kms:kme), prslb(its:ite, kts:kte), prslkb(ims:ime, kms:kme)
  REAL :: prsi(its:ite, kts:kte+1), del(its:ite, kts:kte)
  REAL :: prsib(its:ite, kts:kte+1), delb(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 :: dusfcb(ims:ime), dvsfcb(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 :: rcsksb, tib, rdzb, tem2b, dw2b, shr2b, bvf2b, rdelksb, wtkbjb&
& , temb, gfobnvb, hdb, temcb, tem1b, efactb, dtauxb, dtauyb
!
  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 :: taubb(its:ite), taupb(its:ite, kts:kte+1), xnb(its:ite), ynb(&
& its:ite), ubarb(its:ite), vbarb(its:ite), frb(its:ite), ulowb(its:ite)&
& , rulowb(its:ite), bnvb(its:ite), rollb(its:ite), dtfacb(its:ite), &
& brvfb(its:ite), delksb(its:ite), delks1b(its:ite), bnv2b(its:ite, kts:&
& kte), usqjb(its:ite, kts:kte), taudb(its:ite, kts:kte), rob(its:ite, &
& kts:kte), vtkb(its:ite, kts:kte), vtjb(its:ite, kts:kte), velcob(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 :: zblkb, tautemb
  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 :: taufbb(its:ite, kts:kte+1)
  REAL :: tmp
  REAL :: tmp0
  REAL :: tmp1
  REAL :: tmp2
  REAL :: tmp3
  INTEGER :: branch
  INTEGER :: ad_to
  REAL :: temp3
  REAL :: x3b
  REAL :: temp2
  REAL :: y1b
  REAL :: temp1
  REAL :: temp0
  REAL :: tempb9
  REAL :: tempb8
  REAL :: tempb7
  REAL :: max2b
  REAL :: tempb6
  REAL :: tempb5
  REAL :: tempb4
  REAL :: tempb19
  REAL :: tempb3
  REAL :: tempb18
  REAL :: tempb2
  REAL :: tempb17
  REAL :: tempb1
  REAL :: tempb16
  REAL :: tempb0
  REAL :: tempb15
  REAL :: tempb14
  REAL :: tempb13
  REAL :: tmpb3
  REAL :: tmpb
  REAL :: tempb12
  REAL :: tmpb2
  REAL :: tmpb1
  REAL :: tempb11
  REAL :: tmpb0
  REAL :: tempb10
  REAL :: x4
  REAL :: x3
  REAL :: x2
  INTEGER :: x1
  REAL :: x2b
  REAL :: temp12
  REAL :: temp11
  REAL :: temp10
  REAL :: tempb
  REAL :: x4b
  REAL :: max3
  REAL :: max2
  REAL :: max1
  REAL :: temp9
  REAL :: temp8
  REAL :: tempb20
  REAL :: temp7
  REAL :: temp6
  REAL :: temp5
  REAL :: y1
  REAL :: temp4
  DATA nwdir /6, 7, 5, 8, 2, 3, 1, 4/
!
!---- constants
!
  rcs = SQRT(rcl)
  cs = 1./SQRT(rcl)
  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
  dxy4(3) = SQRT(delx*delx + dely*dely)
  dxy4(4) = dxy4(3)
  dxy4p(1) = dxy4(2)
  dxy4p(2) = dxy4(1)
  dxy4p(3) = dxy4(4)
  dxy4p(4) = dxy4(3)
!
!
!-----initialize arrays
!
  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.
  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
    END DO
  END DO
!
  DO i=its,ite
    taup(i, kte+1) = 0.0
  END DO
!
!  initialize array for flow-blocking drag
!
  taufb(its:ite, kts:kte+1) = 0.0
!
  DO k=kts,kte
    DO i=its,ite
      vtj(i, k) = t1(i, k)*(1.+fv*q1(i, k))
      vtk(i, k) = vtj(i, k)/prslk(i, k)
! density kg/m**3
      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(:)
!
  DO i=its,ite
    delks(i) = 1.0/(prsi(i, 1)-prsi(i, kbl(i)))
    delks1(i) = 1.0/(prsl(i, 1)-prsl(i, kbl(i)))
  END DO
!
!  compute low level averages within pbl
!
  DO k=kts,kpblmax
    DO i=its,ite
      IF (k .LT. kbl(i)) THEN
        rcsks = rcs*del(i, k)*delks(i)
        rdelks = del(i, k)*delks(i)
! pbl u  mean
        ubar(i) = ubar(i) + rcsks*u1(i, k)
! pbl v  mean
        vbar(i) = vbar(i) + rcsks*v1(i, k)
! ro mean
        roll(i) = roll(i) + rdelks*ro(i, k)
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      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
!
!---  saving richardson number in usqj for migwdi
!
  DO k=kts,kte-1
    DO i=its,ite
      ti = 2.0/(t1(i, k)+t1(i, k+1))
      rdz = 1./(zl(i, k+1)-zl(i, k))
      CALL PUSHREAL8(tem1)
      tem1 = u1(i, k) - u1(i, k+1)
      CALL PUSHREAL8(tem2)
      tem2 = v1(i, k) - v1(i, k+1)
      dw2 = rcl*(tem1*tem1+tem2*tem2)
      IF (dw2 .LT. dw2min) THEN
        CALL PUSHREAL8(max2)
        max2 = dw2min
        CALL PUSHCONTROL1B(0)
      ELSE
        CALL PUSHREAL8(max2)
        max2 = dw2
        CALL PUSHCONTROL1B(1)
      END IF
      shr2 = max2*rdz*rdz
      CALL PUSHREAL8(bvf2)
      bvf2 = g*(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*ti
      IF (bvf2/shr2 .LT. rimin) THEN
        usqj(i, k) = rimin
        CALL PUSHCONTROL1B(0)
      ELSE
        usqj(i, k) = bvf2/shr2
        CALL PUSHCONTROL1B(1)
      END IF
      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
        bnv2(i, k) = bnv2min
        CALL PUSHCONTROL1B(0)
      ELSE
        CALL PUSHCONTROL1B(1)
        bnv2(i, k) = bnv2(i, k)
      END IF
    END DO
  END DO
!
!----compute the "low level" or 1/3 wind magnitude (m/s)
!
  DO i=its,ite
    x2 = SQRT(ubar(i)*ubar(i) + vbar(i)*vbar(i))
    IF (x2 .LT. 1.0) THEN
      ulow(i) = 1.0
      CALL PUSHCONTROL1B(0)
    ELSE
      ulow(i) = x2
      CALL PUSHCONTROL1B(1)
    END IF
    rulow(i) = 1./ulow(i)
  END DO
!
  DO k=kts,kte-1
    DO i=its,ite
      velco(i, k) = 0.5*rcs*((u1(i, k)+u1(i, k+1))*ubar(i)+(v1(i, k)+v1(&
&       i, k+1))*vbar(i))
      CALL PUSHREAL8(velco(i, k))
      velco(i, k) = velco(i, k)*rulow(i)
      IF (velco(i, k) .LT. veleps .AND. velco(i, k) .GT. 0.) THEN
        velco(i, k) = veleps
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      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
    wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i)
    CALL PUSHREAL8(bnv2(i, 1))
    bnv2(i, 1) = wtkbj*bnv2(i, 1)
    CALL PUSHREAL8(usqj(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
        rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i)
        tmp = bnv2(i, 1) + bnv2(i, k)*rdelks
        CALL PUSHREAL8(bnv2(i, 1))
        bnv2(i, 1) = tmp
        tmp0 = usqj(i, 1) + usqj(i, k)*rdelks
        CALL PUSHREAL8(usqj(i, 1))
        usqj(i, 1) = tmp0
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      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
        tmp1 = usqj(i, 1)
        CALL PUSHREAL8(usqj(i, k))
        usqj(i, k) = tmp1
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      END IF
    END DO
  END DO
!
  DO i=its,ite
    IF (.NOT.ldrag(i)) THEN
      bnv(i) = SQRT(bnv2(i, 1))
      fr(i) = bnv(i)*rulow(i)*2.*var(i)*od(i)
      IF (fr(i) .GT. frmax) THEN
        fr(i) = frmax
        CALL PUSHCONTROL1B(0)
      ELSE
        CALL PUSHCONTROL1B(1)
        fr(i) = fr(i)
      END IF
      xn(i) = ubar(i)*rulow(i)
      yn(i) = vbar(i)*rulow(i)
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
  END DO
!
!  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
      CALL PUSHREAL8(efact)
      efact = (oa(i)+2.)**(ce*fr(i)/frc)
      IF (efact .LT. efmin) THEN
        x3 = efmin
        CALL PUSHCONTROL1B(0)
      ELSE
        x3 = efact
        CALL PUSHCONTROL1B(1)
      END IF
      IF (x3 .GT. efmax) THEN
        efact = efmax
        CALL PUSHCONTROL1B(0)
      ELSE
        efact = x3
        CALL PUSHCONTROL1B(1)
      END IF
!!!!!!! cleff (effective grid length) is highly tunable parameter
!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
      cleff = SQRT(dxy(i)**2. + dxyp(i)**2.)
      IF (dxmeter .LT. cleff) THEN
        max3 = cleff
      ELSE
        max3 = dxmeter
      END IF
      cleff = 3.*max3
      coefm(i) = (1.+ol(i))**(oa(i)+1.)
      xlinv(i) = coefm(i)/cleff
      tem = fr(i)*fr(i)*oc1(i)
      gfobnv = gmax*tem/((tem+cg)*bnv(i))
      taub(i) = xlinv(i)*roll(i)*ulow(i)*ulow(i)*ulow(i)*gfobnv*efact
      CALL PUSHCONTROL1B(1)
    ELSE
      taub(i) = 0.0
      xn(i) = 0.0
      yn(i) = 0.0
      CALL PUSHCONTROL1B(0)
    END IF
  END DO
!
!   now compute vertical structure of the stress.
!
  DO k=kts,kpblmax
    DO i=its,ite
      IF (k .LE. kbl(i)) THEN
        taup(i, k) = taub(i)
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      END IF
    END DO
  END DO
!
! vertical level k loop!
  DO k=kpblmin,kte-1
    CALL PUSHINTEGER4(kp1)
    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
          CALL PUSHREAL8(brvf(i))
          brvf(i) = bnv2min
          CALL PUSHCONTROL1B(0)
        ELSE
          CALL PUSHREAL8(brvf(i))
          brvf(i) = bnv2(i, k)
          CALL PUSHCONTROL1B(1)
        END IF
! brunt-vaisala frequency
        CALL PUSHREAL8(brvf(i))
        brvf(i) = SQRT(brvf(i))
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      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)
          CALL PUSHREAL8(tem1)
          tem1 = coefm(i)/dxy(i)*(ro(i, kp1)+ro(i, k))*brvf(i)*velco(i, &
&           k)*0.5
          CALL PUSHREAL8(hd)
          hd = SQRT(taup(i, k)/tem1)
          fro = brvf(i)*hd*temv
!
!  rim is the  minimum-richardson number by shutts (1985)
!
          CALL PUSHREAL8(tem2)
          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
              temc = 2.0 + 1.0/tem2
              hd = velco(i, k)*(2.*SQRT(temc)-temc)/brvf(i)
              taup(i, kp1) = tem1*hd*hd
              CALL PUSHCONTROL3B(4)
            ELSE
              CALL PUSHCONTROL3B(3)
            END IF
          ELSE
! no wavebreaking!
            tmp2 = taup(i, k)
            taup(i, kp1) = tmp2
            CALL PUSHCONTROL3B(2)
          END IF
        ELSE
          CALL PUSHCONTROL3B(1)
        END IF
      ELSE
        CALL PUSHCONTROL3B(0)
      END IF
    END DO
  END DO
!
  IF (lcap .LT. kte) THEN
    DO klcap=lcapp1,kte
      DO i=its,ite
        tmp3 = prsi(i, klcap)/prsi(i, lcap)*taup(i, lcap)
        CALL PUSHREAL8(taup(i, klcap))
        taup(i, klcap) = tmp3
      END DO
    END DO
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  DO i=its,ite
    IF (.NOT.ldrag(i)) THEN
!
!------- determine the height of flow-blocking layer
!
      CALL PUSHINTEGER4(kblk)
      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
            CALL PUSHINTEGER4(kblk)
            kblk = k
            IF (kblk .GT. kbl(i)) THEN
              kblk = kbl(i)
            ELSE
              kblk = kblk
            END IF
            CALL PUSHREAL8(zblk)
            zblk = zl(i, kblk) - zl(i, kts)
            CALL PUSHCONTROL2B(2)
          ELSE
            CALL PUSHCONTROL2B(1)
          END IF
        ELSE
          CALL PUSHCONTROL2B(0)
        END IF
      END DO
      IF (kblk .NE. 0) THEN
!
!--------- compute flow-blocking stress
!
        IF (2.0 - 1.0/od(i) .LT. 0.0) THEN
          CALL PUSHREAL8(cd)
          cd = 0.0
          CALL PUSHCONTROL1B(0)
        ELSE
          CALL PUSHREAL8(cd)
          cd = 2.0 - 1.0/od(i)
          CALL PUSHCONTROL1B(1)
        END IF
        taufb(i, kts) = 0.5*roll(i)*coefm(i)/dxy(i)**2*cd*dxyp(i)*olp(i)&
&         *zblk*ulow(i)**2
        tautem = taufb(i, kts)/FLOAT(kblk-kts)
        DO k=kts+1,kblk
          taufb(i, k) = taufb(i, k-1) - tautem
        END DO
        CALL PUSHINTEGER4(k - 1)
!
!----------sum orographic GW stress and flow-blocking stress
!
        CALL PUSHREAL8ARRAY(taup(i, :), kte - kts + 2)
        taup(i, :) = taup(i, :) + taufb(i, :)
        CALL PUSHCONTROL2B(2)
      ELSE
        CALL PUSHCONTROL2B(1)
      END IF
    ELSE
      CALL PUSHCONTROL2B(0)
    END IF
  END DO
!
!  calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
!
  DO k=kts,kte
    DO i=its,ite
      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
      taud(i, klcap) = taud(i, klcap)*factop
    END DO
  END DO
!
!  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
          x4 = velco(i, k)/(deltim*rcs*taud(i, k))
          IF (x4 .GE. 0.) THEN
            y1 = x4
            CALL PUSHCONTROL1B(0)
          ELSE
            y1 = -x4
            CALL PUSHCONTROL1B(1)
          END IF
          IF (dtfac(i) .GT. y1) THEN
            dtfac(i) = y1
            CALL PUSHCONTROL2B(2)
          ELSE
            dtfac(i) = dtfac(i)
            CALL PUSHCONTROL2B(3)
          END IF
        ELSE
          CALL PUSHCONTROL2B(1)
        END IF
      ELSE
        CALL PUSHCONTROL2B(0)
      END IF
    END DO
  END DO
!
  DO k=kts,kte
    DO i=its,ite
      CALL PUSHREAL8(taud(i, k))
      taud(i, k) = taud(i, k)*dtfac(i)
    END DO
  END DO
  DO i=ite,its,-1
    dvsfcb(i) = -(rcs*dvsfcb(i)/g)
    dusfcb(i) = -(rcs*dusfcb(i)/g)
  END DO
  taudb = 0.0
  xnb = 0.0
  ynb = 0.0
  dtfacb = 0.0
  DO k=kte,kts,-1
    DO i=ite,its,-1
      dtaux = taud(i, k)*xn(i)
      dtauy = taud(i, k)*yn(i)
      dtauyb = dvdtb(i, k) + dtauy2db(i, k) + del(i, k)*dvsfcb(i)
      delb(i, k) = delb(i, k) + dtaux*dusfcb(i) + dtauy*dvsfcb(i)
      dtauxb = dudtb(i, k) + dtaux2db(i, k) + del(i, k)*dusfcb(i)
      dtauy2db(i, k) = 0.0
      dtaux2db(i, k) = 0.0
      taudb(i, k) = taudb(i, k) + xn(i)*dtauxb + yn(i)*dtauyb
      ynb(i) = ynb(i) + taud(i, k)*dtauyb
      xnb(i) = xnb(i) + taud(i, k)*dtauxb
      CALL POPREAL8(taud(i, k))
      dtfacb(i) = dtfacb(i) + taud(i, k)*taudb(i, k)
      taudb(i, k) = dtfac(i)*taudb(i, k)
    END DO
  END DO
  DO i=ite,its,-1
    dvsfcb(i) = 0.0
    dusfcb(i) = 0.0
  END DO
  velcob = 0.0
  DO k=kpblmax-1,kts,-1
    DO i=ite,its,-1
      CALL POPCONTROL2B(branch)
      IF (branch .GE. 2) THEN
        IF (branch .EQ. 2) THEN
          y1b = dtfacb(i)
          dtfacb(i) = 0.0
        ELSE
          y1b = 0.0
        END IF
        CALL POPCONTROL1B(branch)
        IF (branch .EQ. 0) THEN
          x4b = y1b
        ELSE
          x4b = -y1b
        END IF
        temp12 = deltim*rcs*taud(i, k)
        velcob(i, k) = velcob(i, k) + x4b/temp12
        taudb(i, k) = taudb(i, k) - velco(i, k)*deltim*rcs*x4b/temp12**2
      END IF
    END DO
  END DO
  DO klcap=kte,lcap,-1
    DO i=ite,its,-1
      taudb(i, klcap) = factop*taudb(i, klcap)
    END DO
  END DO
  taupb = 0.0
  DO k=kte,kts,-1
    DO i=ite,its,-1
      tempb20 = csg*taudb(i, k)/del(i, k)
      taupb(i, k+1) = taupb(i, k+1) + tempb20
      taupb(i, k) = taupb(i, k) - tempb20
      delb(i, k) = delb(i, k) - (taup(i, k+1)-taup(i, k))*tempb20/del(i&
&       , k)
      taudb(i, k) = 0.0
    END DO
  END DO
  rollb = 0.0
  zblkb = 0.0
  taufbb = 0.0
  ulowb = 0.0
  DO i=ite,its,-1
    CALL POPCONTROL2B(branch)
    IF (branch .NE. 0) THEN
      IF (branch .NE. 1) THEN
        CALL POPREAL8ARRAY(taup(i, :), kte - kts + 2)
        taufbb(i, :) = taufbb(i, :) + taupb(i, :)
        tautemb = 0.0
        CALL POPINTEGER4(ad_to)
        DO k=ad_to,kts+1,-1
          taufbb(i, k-1) = taufbb(i, k-1) + taufbb(i, k)
          tautemb = tautemb - taufbb(i, k)
          taufbb(i, k) = 0.0
        END DO
        taufbb(i, kts) = taufbb(i, kts) + tautemb/FLOAT(kblk-kts)
        tempb18 = cd*0.5*olp(i)*coefm(i)*dxyp(i)*taufbb(i, kts)
        tempb19 = ulow(i)**2*tempb18/dxy(i)**2
        rollb(i) = rollb(i) + zblk*tempb19
        zblkb = zblkb + roll(i)*tempb19
        ulowb(i) = ulowb(i) + roll(i)*zblk*2*ulow(i)*tempb18/dxy(i)**2
        taufbb(i, kts) = 0.0
        CALL POPCONTROL1B(branch)
        IF (branch .EQ. 0) THEN
          CALL POPREAL8(cd)
        ELSE
          CALL POPREAL8(cd)
        END IF
      END IF
      DO k=kpblmin,kte,1
        CALL POPCONTROL2B(branch)
        IF (branch .NE. 0) THEN
          IF (branch .NE. 1) THEN
            CALL POPREAL8(zblk)
            zlb(i, kblk) = zlb(i, kblk) + zblkb
            zlb(i, kts) = zlb(i, kts) - zblkb
            CALL POPINTEGER4(kblk)
            zblkb = 0.0
          END IF
        END IF
      END DO
      CALL POPINTEGER4(kblk)
    END IF
  END DO
  CALL POPCONTROL1B(branch)
  IF (branch .NE. 0) THEN
    DO klcap=kte,lcapp1,-1
      DO i=ite,its,-1
        CALL POPREAL8(taup(i, klcap))
        tmpb3 = taupb(i, klcap)
        taupb(i, klcap) = 0.0
        tempb17 = tmpb3/prsi(i, lcap)
        prsib(i, klcap) = prsib(i, klcap) + taup(i, lcap)*tempb17
        taupb(i, lcap) = taupb(i, lcap) + prsi(i, klcap)*tempb17
        prsib(i, lcap) = prsib(i, lcap) - prsi(i, klcap)*taup(i, lcap)*&
&         tempb17/prsi(i, lcap)
      END DO
    END DO
  END IF
  brvfb = 0.0
  bnv2b = 0.0
  rob = 0.0
  usqjb = 0.0
  DO k=kte-1,kpblmin,-1
    DO i=ite,its,-1
      CALL POPCONTROL3B(branch)
      IF (branch .GE. 2) THEN
        IF (branch .EQ. 2) THEN
          tmpb2 = taupb(i, kp1)
          taupb(i, kp1) = 0.0
          taupb(i, k) = taupb(i, k) + tmpb2
          tem1b = 0.0
          tem2b = 0.0
        ELSE IF (branch .EQ. 3) THEN
          tem1b = 0.0
          tem2b = 0.0
        ELSE
          tem1b = hd**2*taupb(i, kp1)
          hdb = tem1*2*hd*taupb(i, kp1)
          taupb(i, kp1) = 0.0
          temc = 2.0 + 1.0/tem2
          temp11 = SQRT(temc)
          tempb16 = (2.*temp11-temc)*hdb/brvf(i)
          temp10 = velco(i, k)/brvf(i)
          velcob(i, k) = velcob(i, k) + tempb16
          brvfb(i) = brvfb(i) - temp10*tempb16
          IF (temc .EQ. 0.0) THEN
            temcb = -(temp10*hdb)
          ELSE
            temcb = (2.*temp10/(2.0*temp11)-temp10)*hdb
          END IF
          tem2b = -(temcb/tem2**2)
        END IF
        CALL POPREAL8(tem2)
        IF (.NOT.usqj(i, k) .EQ. 0.0) usqjb(i, k) = usqjb(i, k) + tem2b/&
&           (2.0*SQRT(usqj(i, k)))
        CALL POPREAL8(hd)
        CALL POPREAL8(tem1)
        temp9 = brvf(i)/dxy(i)
        tempb13 = coefm(i)*0.5*tem1b
        tempb14 = velco(i, k)*temp9*tempb13
        tempb15 = (ro(i, kp1)+ro(i, k))*tempb13
        rob(i, kp1) = rob(i, kp1) + tempb14
        rob(i, k) = rob(i, k) + tempb14
        velcob(i, k) = velcob(i, k) + temp9*tempb15
        brvfb(i) = brvfb(i) + velco(i, k)*tempb15/dxy(i)
      END IF
    END DO
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) THEN
        CALL POPREAL8(brvf(i))
        IF (brvf(i) .EQ. 0.0) THEN
          brvfb(i) = 0.0
        ELSE
          brvfb(i) = brvfb(i)/(2.0*SQRT(brvf(i)))
        END IF
        CALL POPCONTROL1B(branch)
        IF (branch .EQ. 0) THEN
          CALL POPREAL8(brvf(i))
          brvfb(i) = 0.0
        ELSE
          CALL POPREAL8(brvf(i))
          bnv2b(i, k) = bnv2b(i, k) + brvfb(i)
          brvfb(i) = 0.0
        END IF
      END IF
    END DO
    CALL POPINTEGER4(kp1)
  END DO
  taubb = 0.0
  DO k=kpblmax,kts,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) THEN
        taubb(i) = taubb(i) + taupb(i, k)
        taupb(i, k) = 0.0
      END IF
    END DO
  END DO
  bnvb = 0.0
  frb = 0.0
  DO i=ite,its,-1
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      ynb(i) = 0.0
      xnb(i) = 0.0
      taubb(i) = 0.0
    ELSE
      tem = fr(i)*fr(i)*oc1(i)
      gfobnv = gmax*tem/((tem+cg)*bnv(i))
      tempb10 = xlinv(i)*ulow(i)**3*taubb(i)
      ulowb(i) = ulowb(i) + roll(i)*gfobnv*efact*xlinv(i)*3*ulow(i)**2*&
&       taubb(i)
      rollb(i) = rollb(i) + gfobnv*efact*tempb10
      gfobnvb = roll(i)*efact*tempb10
      efactb = roll(i)*gfobnv*tempb10
      taubb(i) = 0.0
      temp8 = (cg+tem)*bnv(i)
      tempb11 = gmax*gfobnvb/temp8
      tempb12 = -(tem*tempb11/temp8)
      temb = bnv(i)*tempb12 + tempb11
      bnvb(i) = bnvb(i) + (cg+tem)*tempb12
      frb(i) = frb(i) + oc1(i)*2*fr(i)*temb
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) THEN
        x3b = 0.0
      ELSE
        x3b = efactb
      END IF
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) THEN
        efactb = 0.0
      ELSE
        efactb = x3b
      END IF
      CALL POPREAL8(efact)
      IF (.NOT.oa(i) + 2. .LE. 0.0) frb(i) = frb(i) + ce*(oa(i)+2.)**(ce&
&         *(fr(i)/frc))*LOG(oa(i)+2.)*efactb/frc
    END IF
  END DO
  vbarb = 0.0
  ubarb = 0.0
  rulowb = 0.0
  DO i=ite,its,-1
    CALL POPCONTROL1B(branch)
    IF (branch .NE. 0) THEN
      vbarb(i) = vbarb(i) + rulow(i)*ynb(i)
      rulowb(i) = rulowb(i) + ubar(i)*xnb(i) + vbar(i)*ynb(i)
      ynb(i) = 0.0
      ubarb(i) = ubarb(i) + rulow(i)*xnb(i)
      xnb(i) = 0.0
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) frb(i) = 0.0
      tempb9 = var(i)*2.*od(i)*frb(i)
      bnvb(i) = bnvb(i) + rulow(i)*tempb9
      rulowb(i) = rulowb(i) + bnv(i)*tempb9
      frb(i) = 0.0
      IF (.NOT.bnv2(i, 1) .EQ. 0.0) bnv2b(i, 1) = bnv2b(i, 1) + bnvb(i)/&
&         (2.0*SQRT(bnv2(i, 1)))
      bnvb(i) = 0.0
    END IF
  END DO
  DO k=kpblmax,kpblmin,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) THEN
        CALL POPREAL8(usqj(i, k))
        tmpb1 = usqjb(i, k)
        usqjb(i, k) = 0.0
        usqjb(i, 1) = usqjb(i, 1) + tmpb1
      END IF
    END DO
  END DO
  delks1b = 0.0
  DO k=kpblmax,kpblmin,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) THEN
        tmpb0 = bnv2b(i, 1)
        rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i)
        CALL POPREAL8(usqj(i, 1))
        tmpb = usqjb(i, 1)
        usqjb(i, 1) = tmpb
        usqjb(i, k) = usqjb(i, k) + rdelks*tmpb
        CALL POPREAL8(bnv2(i, 1))
        rdelksb = bnv2(i, k)*tmpb0 + usqj(i, k)*tmpb
        bnv2b(i, 1) = tmpb0
        bnv2b(i, k) = bnv2b(i, k) + rdelks*tmpb0
        prslb(i, k) = prslb(i, k) + delks1(i)*rdelksb
        prslb(i, k+1) = prslb(i, k+1) - delks1(i)*rdelksb
        delks1b(i) = delks1b(i) + (prsl(i, k)-prsl(i, k+1))*rdelksb
      END IF
    END DO
  END DO
  DO i=ite,its,-1
    wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i)
    CALL POPREAL8(usqj(i, 1))
    CALL POPREAL8(bnv2(i, 1))
    wtkbjb = bnv2(i, 1)*bnv2b(i, 1) + usqj(i, 1)*usqjb(i, 1)
    usqjb(i, 1) = wtkbj*usqjb(i, 1)
    bnv2b(i, 1) = wtkbj*bnv2b(i, 1)
    prslb(i, 1) = prslb(i, 1) + delks1(i)*wtkbjb
    prslb(i, 2) = prslb(i, 2) - delks1(i)*wtkbjb
    delks1b(i) = delks1b(i) + (prsl(i, 1)-prsl(i, 2))*wtkbjb
  END DO
  DO k=kte-1,kts,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) velcob(i, k) = 0.0
      CALL POPREAL8(velco(i, k))
      rulowb(i) = rulowb(i) + velco(i, k)*velcob(i, k)
      velcob(i, k) = rulow(i)*velcob(i, k)
      tempb8 = rcs*0.5*velcob(i, k)
      u1b(i, k) = u1b(i, k) + ubar(i)*tempb8
      u1b(i, k+1) = u1b(i, k+1) + ubar(i)*tempb8
      ubarb(i) = ubarb(i) + (u1(i, k)+u1(i, k+1))*tempb8
      v1b(i, k) = v1b(i, k) + vbar(i)*tempb8
      v1b(i, k+1) = v1b(i, k+1) + vbar(i)*tempb8
      vbarb(i) = vbarb(i) + (v1(i, k)+v1(i, k+1))*tempb8
      velcob(i, k) = 0.0
    END DO
  END DO
  DO i=ite,its,-1
    ulowb(i) = ulowb(i) - rulowb(i)/ulow(i)**2
    rulowb(i) = 0.0
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      ulowb(i) = 0.0
      x2b = 0.0
    ELSE
      x2b = ulowb(i)
      ulowb(i) = 0.0
    END IF
    IF (ubar(i)**2 + vbar(i)**2 .EQ. 0.0) THEN
      tempb7 = 0.0
    ELSE
      tempb7 = x2b/(2.0*SQRT(ubar(i)**2+vbar(i)**2))
    END IF
    ubarb(i) = ubarb(i) + 2*ubar(i)*tempb7
    vbarb(i) = vbarb(i) + 2*vbar(i)*tempb7
  END DO
  vtjb = 0.0
  vtkb = 0.0
  DO k=kte-1,kts,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) bnv2b(i, k) = 0.0
      rdz = 1./(zl(i, k+1)-zl(i, k))
      temp6 = vtk(i, k+1) + vtk(i, k)
      temp7 = vtk(i, k+1) - vtk(i, k)
      tempb5 = g*2.0*bnv2b(i, k)/temp6
      tempb6 = -(rdz*temp7*tempb5/temp6)
      rdzb = temp7*tempb5
      vtkb(i, k+1) = vtkb(i, k+1) + tempb6 + rdz*tempb5
      vtkb(i, k) = vtkb(i, k) + tempb6 - rdz*tempb5
      bnv2b(i, k) = 0.0
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) THEN
        usqjb(i, k) = 0.0
        bvf2b = 0.0
        shr2b = 0.0
      ELSE
        shr2 = max2*rdz*rdz
        bvf2b = usqjb(i, k)/shr2
        shr2b = -(bvf2*usqjb(i, k)/shr2**2)
        usqjb(i, k) = 0.0
      END IF
      ti = 2.0/(t1(i, k)+t1(i, k+1))
      CALL POPREAL8(bvf2)
      temp5 = vtj(i, k+1) - vtj(i, k)
      tempb4 = g*ti*bvf2b
      rdzb = rdzb + max2*2*rdz*shr2b + temp5*tempb4
      vtjb(i, k+1) = vtjb(i, k+1) + rdz*tempb4
      vtjb(i, k) = vtjb(i, k) - rdz*tempb4
      tib = g*(g/cp+rdz*temp5)*bvf2b
      max2b = rdz**2*shr2b
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) THEN
        CALL POPREAL8(max2)
        dw2b = 0.0
      ELSE
        CALL POPREAL8(max2)
        dw2b = max2b
      END IF
      tem1b = rcl*2*tem1*dw2b
      tem2b = rcl*2*tem2*dw2b
      CALL POPREAL8(tem2)
      v1b(i, k) = v1b(i, k) + tem2b
      v1b(i, k+1) = v1b(i, k+1) - tem2b
      CALL POPREAL8(tem1)
      u1b(i, k) = u1b(i, k) + tem1b
      u1b(i, k+1) = u1b(i, k+1) - tem1b
      temp4 = zl(i, k+1) - zl(i, k)
      tempb2 = -(rdzb/temp4**2)
      zlb(i, k+1) = zlb(i, k+1) + tempb2
      zlb(i, k) = zlb(i, k) - tempb2
      temp3 = t1(i, k) + t1(i, k+1)
      tempb3 = -(2.0*tib/temp3**2)
      t1b(i, k) = t1b(i, k) + tempb3
      t1b(i, k+1) = t1b(i, k+1) + tempb3
    END DO
  END DO
  delksb = 0.0
  DO k=kpblmax,kts,-1
    DO i=ite,its,-1
      CALL POPCONTROL1B(branch)
      IF (branch .NE. 0) THEN
        rdelks = del(i, k)*delks(i)
        rdelksb = ro(i, k)*rollb(i)
        rob(i, k) = rob(i, k) + rdelks*rollb(i)
        rcsks = rcs*del(i, k)*delks(i)
        rcsksb = u1(i, k)*ubarb(i) + v1(i, k)*vbarb(i)
        v1b(i, k) = v1b(i, k) + rcsks*vbarb(i)
        u1b(i, k) = u1b(i, k) + rcsks*ubarb(i)
        delb(i, k) = delb(i, k) + rcs*delks(i)*rcsksb + delks(i)*rdelksb
        delksb(i) = delksb(i) + rcs*del(i, k)*rcsksb + del(i, k)*rdelksb
      END IF
    END DO
  END DO
  DO i=ite,its,-1
    temp2 = prsl(i, 1) - prsl(i, kbl(i))
    tempb0 = -(delks1b(i)/temp2**2)
    prslb(i, 1) = prslb(i, 1) + tempb0
    prslb(i, kbl(i)) = prslb(i, kbl(i)) - tempb0
    delks1b(i) = 0.0
    temp1 = prsi(i, 1) - prsi(i, kbl(i))
    tempb1 = -(delksb(i)/temp1**2)
    prsib(i, 1) = prsib(i, 1) + tempb1
    prsib(i, kbl(i)) = prsib(i, kbl(i)) - tempb1
    delksb(i) = 0.0
  END DO
  DO k=kte,kts,-1
    DO i=ite,its,-1
      tempb = vtkb(i, k)/prslk(i, k)
      temp0 = rd*vtj(i, k)
      prslb(i, k) = prslb(i, k) + rob(i, k)/temp0
      vtjb(i, k) = vtjb(i, k) + tempb - prsl(i, k)*rd*rob(i, k)/temp0**2
      rob(i, k) = 0.0
      prslkb(i, k) = prslkb(i, k) - vtj(i, k)*tempb/prslk(i, k)
      vtkb(i, k) = 0.0
      t1b(i, k) = t1b(i, k) + (fv*q1(i, k)+1.)*vtjb(i, k)
      q1b(i, k) = q1b(i, k) + t1(i, k)*fv*vtjb(i, k)
      vtjb(i, k) = 0.0
    END DO
  END DO
  DO k=kte,kts,-1
    DO i=ite,its,-1
      dtauy2db(i, k) = 0.0
      dtaux2db(i, k) = 0.0
    END DO
  END DO
END SUBROUTINE GWDO2D_B

end module a_module_bl_gwdo