      subroutine clinic (joff, js, je, is, ie)

#if defined O_mom
!=======================================================================
!     compute internal mode velocity components for rows js through je
!     in the MW.

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!=======================================================================

      implicit none

      integer istrt, iend
      integer i, k, j, jrow, n, js, je, limit, joff, kb, is, ie

      real adv_ux, adv_uy, adv_uz, adv_metric, diff_ux, diff_uz, fx
      real diff_uy, diff_metric, coriolis, aprime, grav_rho0r, fxa
      real fxb, t1, t2, ambi_csur, del2, ambi_cst_dytr, detmr, unep

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
# if defined O_neptune
      include "cnep.h"
# endif
      include "coord.h"
      include "csbc.h"
      include "grdvar.h"
      include "hmixc.h"
      include "emode.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "vmixc.h"
      include "fdifm.h"

      parameter (istrt=2, iend=imt-1)

      real tempik(imt,km,jsmw:jmw)
      real baru(imt,jsmw:jemw,2)

!-----------------------------------------------------------------------
!     bail out if starting row exceeds ending row
!-----------------------------------------------------------------------

      if (js .gt. je) return

!-----------------------------------------------------------------------
!     limit the longitude indices based on those from the argument list
!     Note: these are currently bypassed. istrt and iend are set as
!           parameters to optimize performance
!-----------------------------------------------------------------------

!      istrt = max(2,is)
!      iend  = min(imt-1,ie)

!-----------------------------------------------------------------------
!     build coefficients to minimize advection and diffusion computation
!-----------------------------------------------------------------------

# if defined O_biharmonic
      limit = min(je+1+joff,jmt) - joff
      do j=js,limit
# else
      do j=js,je
# endif
        jrow = j + joff
# if defined O_anisotropic_viscosity
        do k=1,km
           do i=istrt-1,iend
              csudxur(i,j)  = csur(jrow)*dxur(i)
              csudxu2r(i,j) = csur(jrow)*dxur(i)*p5
              am_csudxtr(i,k,j) = visc_ceu(i,k,j)*csur(jrow)*dxtr(i+1)
           enddo
        enddo
# else
        do i=istrt-1,iend
          csudxur(i,j)  = csur(jrow)*dxur(i)
          csudxu2r(i,j) = csur(jrow)*dxur(i)*p5
#  if defined O_consthmix
          am_csudxtr(i,j)  = am*csur(jrow)*dxtr(i+1)
#  endif
        enddo
# endif
      enddo

!-----------------------------------------------------------------------
!     construct the hydrostatic pressure gradients: 1 = dp/dx; 2 = dp/dy
!-----------------------------------------------------------------------

!     compute horizontal pressure gradient at the first level

# if defined O_pressure_gradient_average

!     construct density as in Brown and Campana (1978) (see manual)

      call state (t(1,1,1,1,taum1), t(1,1,1,2,taum1), rhotaum1(1,1,jsmw)
     &,           js, je+1, istrt-1, iend+1)
      call state (t(1,1,1,1,taup1), t(1,1,1,2,taup1), rhotaup1(1,1,jsmw)
     &,           js, je+1, istrt-1, iend+1)
      aprime = 0.25
      do j=js,je+1
        do k=1,km
          do i=1,imt
            rhotilde(i,k,j) = aprime*(rhotaum1(i,k,j) + rhotaup1(i,k,j))
     &       + (1.0-2.0*aprime)*rho(i,k,j)
          enddo
        enddo
      enddo
# endif

      grav_rho0r = grav*rho0r
      do j=js,je
        jrow = j + joff
        fxa  = grav_rho0r*dzw(0)*csur(jrow)
        fxb  = grav_rho0r*dzw(0)*dyu2r(jrow)
        do i=istrt-1,iend
# if defined O_pressure_gradient_average
          t1              = rhotilde(i+1,1,j+1) - rhotilde(i  ,1,j)
          t2              = rhotilde(i  ,1,j+1) - rhotilde(i+1,1,j)
# else
          t1              = rho(i+1,1,j+1) - rho(i  ,1,j)
          t2              = rho(i  ,1,j+1) - rho(i+1,1,j)
# endif
          grad_p(i,1,j,1) = (t1-t2)*fxa*dxu2r(i)
          grad_p(i,1,j,2) = (t1+t2)*fxb
        enddo
      enddo

!     compute the change in pressure gradient between levels

      do j=js,je+1
        do k=2,km
          do i=istrt-1,iend+1
# if defined O_pressure_gradient_average
            tempik(i,k,j) = rhotilde(i,k-1,j) + rhotilde(i,k,j)
# else
            tempik(i,k,j) = rho(i,k-1,j) + rho(i,k,j)
# endif
          enddo
        enddo
      enddo

      do j=js,je
        jrow = j + joff
        fxa = grav_rho0r*csur(jrow)*p5
        fxb = grav_rho0r*dyu4r(jrow)
        do k=2,km
          do i=istrt-1,iend
            t1              = tempik(i+1,k,j+1) - tempik(i  ,k,j)
            t2              = tempik(i  ,k,j+1) - tempik(i+1,k,j)
            grad_p(i,k,j,1) = fxa*(t1-t2)*dzw(k-1)*dxu2r(i)
            grad_p(i,k,j,2) = fxb*(t1+t2)*dzw(k-1)
          enddo
        enddo
      enddo

!     integrate downward from the first level

      do j=js,je
        do k=1,kmm1
          do i=istrt-1,iend
            grad_p(i,k+1,j,1) = grad_p(i,k,j,1) + grad_p(i,k+1,j,1)
            grad_p(i,k+1,j,2) = grad_p(i,k,j,2) + grad_p(i,k+1,j,2)
          enddo
        enddo
      enddo

      do j=js,je
        call setbcx (grad_p(1,1,j,1), imt, km)
        call setbcx (grad_p(1,1,j,2), imt, km)
      enddo

!-----------------------------------------------------------------------
!     solve for one component of velocity at a time
!     n = 1 => zonal component
!     n = 2 => meridional component
!-----------------------------------------------------------------------

      do n=1,2

# if !defined O_linearized_advection

!-----------------------------------------------------------------------
!       calculate 2*advective flux (for speed) across east face of
!       "u" cells.
!-----------------------------------------------------------------------

        do j=js,je
          do k=1,km
            do i=istrt-1,iend
              adv_fe(i,k,j) = adv_veu(i,k,j)*(u(i,  k,j,n,tau) +
     &                                        u(i+1,k,j,n,tau))
            enddo
          enddo
        enddo

!-----------------------------------------------------------------------
!       2*advective flux across northern face of "u" cells is built
!       into ADV_Uy. (It's done this way for performance issues)
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!       diffusive flux across east face of "u" cell
!       diffusive flux across north face of "u" cell
!-----------------------------------------------------------------------

#  if defined O_consthmix && !defined O_biharmonic

!       build diffusive flux on eastern face of "u" cells

        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=istrt-1,iend
#   if defined O_anisotropic_viscosity
              diff_fe(i,k,j) = am_csudxtr(i,k,j)*
#   else
              diff_fe(i,k,j) = am_csudxtr(i,j)*
#   endif
     &                        (u(i+1,k,j,n,taum1) - u(i,k,j,n,taum1)
#   if defined O_neptune
     &                        - unep(i+1,jrow,n)*umask(i+1,k,j) +
     &                          unep(i,jrow,n)*umask(i,k,j)
#   endif
     &                        )
            enddo
          enddo
        enddo

!       diffusive flux on northern face of "u" cells is built
!       into DIFF_Uy
#  endif
# endif
# if defined O_consthmix && defined O_biharmonic

!-----------------------------------------------------------------------
!       diffusive flux across eastern face of "u" cell
!       diffusive flux across northern face of "u" cell
!-----------------------------------------------------------------------

        do j=js,je
          jrow = j + joff
          ambi_csur = visc_ceu*csur(jrow)
          do k=1,km
            do i=istrt-1,iend
              diff_fe(i,k,j) = ambi_csur*dxtr(i+1)*
     &                        (del2(i+1,k,j,n) - del2(i,k,j,n))
            enddo
          enddo
        enddo
        do j=js-1,je
          jrow = j + joff
          ambi_cst_dytr = visc_cnu*cst(jrow+1)*dytr(jrow+1)
          do k=1,km
            do i=istrt,iend
              diff_fn(i,k,j) = ambi_cst_dytr*
     &                        (del2(i,k,j+1,n) - del2(i,k,j,n))
            enddo
          enddo
        enddo
# endif
# if defined O_smagnlmix && !defined O_consthmix

!       calculate diffusive flux on eastern and northern faces of
!      "u" cells

        call smagnlm (joff, js, je, istrt, iend, n)
# endif

# if !defined O_linearized_advection

!-----------------------------------------------------------------------
!       calculate 2*advective flux (for speed) on bottom face of
!       "u" cell. also diffusive flux on bottom face of "u" cell
!-----------------------------------------------------------------------

        do j=js,je
          do k=1,kmm1
            do i=istrt,iend
              adv_fb(i,k,j)  = adv_vbu(i,k,j)*(u(i,k,  j,n,tau) +
     &                                        u(i,k+1,j,n,tau))
              diff_fb(i,k,j) = visc_cbu(i,k,j)*dzwr(k)*
     &                         (u(i,k,j,n,taum1) - u(i,k+1,j,n,taum1))
            enddo
          enddo
        enddo
# endif

!-----------------------------------------------------------------------
!       set surface and bottom vert b.c. on "u" cells for mixing
!       and advection.
!       note: the b.c. at adv_fb(i,k=bottom,j) is set by the above code.
!             However, it is not set when k=km so it is set below.
!             adv_fb(i,km,j) is always zero (to within roundoff).
!-----------------------------------------------------------------------

        do j=js,je
          jrow = j + joff
          do i=istrt,iend
            kb              = kmu(i,jrow)
            diff_fb(i,0,j)  = smf(i,j,n)
            diff_fb(i,kb,j) = bmf(i,j,n)
            adv_fb(i,0,j)   = adv_vbu(i,0,j)*(u(i,1,j,n,tau) +
     &                                        u(i,1,j,n,tau))
            adv_fb(i,km,j)  = adv_vbu(i,km,j)*u(i,km,j,n,tau)
          enddo
        enddo

# if defined O_source_term || defined O_npzd || defined O_carbon_14

!-----------------------------------------------------------------------
!       set source term for "u" cell
!-----------------------------------------------------------------------

        do j=js,je
          do k=1,km
            do i=istrt,iend
              source(i,k,j) = c0
            enddo
          enddo
        enddo
# endif

!-----------------------------------------------------------------------
!       solve for the internal mode part of du/dt at center of
!       "u" cells by neglecting the surface pressure gradients. use
!       statement functions to represent each component of the
!       calculation.
!-----------------------------------------------------------------------

        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              u(i,k,j,n,taup1) =
     &          (DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j,jrow,n) + DIFF_Uz(i,k,j)
     &          + DIFF_metric(i,k,j,jrow,n)
# if !defined O_linearized_advection
     &          - ADV_Ux(i,k,j) - ADV_Uy(i,k,j,jrow,n) - ADV_Uz(i,k,j)
     &          + ADV_metric(i,k,j,jrow,n)
# endif
     &          - grad_p(i,k,j,n) + CORIOLIS(i,k,j,jrow,n)
# if defined O_source_term || defined O_npzd || defined O_carbon_14
     &          + source(i,k,j)
# endif
     &          )*umask(i,k,j)
            enddo
          enddo
        enddo
# if defined O_implicitvmix

!-----------------------------------------------------------------------
!       add in du/dt component due to implicit vertical diffusion
!-----------------------------------------------------------------------

        call ivdifu (joff, js, je, istrt, iend, n)
# endif

!-----------------------------------------------------------------------
!       construct diagnostics associated with velocity component "n"
!-----------------------------------------------------------------------

        call diagc1 (joff, js, je, istrt, iend, n)

!-----------------------------------------------------------------------
!       construct the vertical average of du/dt for forcing
!       the barotropic equation
!-----------------------------------------------------------------------

        do j=js,je
          jrow = j + joff
          do i=istrt,iend
            zu(i,jrow,n) = c0
          enddo
        enddo
        do j=js,je
          jrow = j + joff
          do k=1,km
            fx = dzt(k)
            do i=istrt,iend
              zu(i,jrow,n) = zu(i,jrow,n) + u(i,k,j,n,taup1)*fx
            enddo
          enddo
        enddo

        do j=js,je
          jrow = j + joff
          do i=istrt,iend
            zu(i,jrow,n) = zu(i,jrow,n)*hr(i,jrow)
          enddo
        enddo
# if defined O_symmetry
        do j=js,je
          jrow = j + joff
          if (jrow .eq. jmtm1 .and. n .eq. 2) then
            do i=istrt,iend
              zu(i,jrow,2)   =  c0
              zu(i,jrow+1,2) = -zu(i,jrow-1,2)
              zu(i,jrow+1,1) =  zu(i,jrow-1,1)
            enddo
          endif
        enddo
# endif

!-----------------------------------------------------------------------
!       end of velocity component "n" loop
!-----------------------------------------------------------------------

      enddo

!-----------------------------------------------------------------------
!     compute "tau+1" velocities accounting for implicit part of the
!     coriolis term if treated implicitly. velocities are in error by an
!     arbitrary constant related to neglecting the unknown surface
!     pressure gradients
!-----------------------------------------------------------------------

# if defined O_damp_inertial_oscillation
      do j=js,je
        jrow  = j + joff
        do k=1,km
          do i=istrt,iend
            fx    = c2dtuv*acor*cori(i,jrow,1)
            detmr = c1/(c1 + fx*fx)
            t1 = (u(i,k,j,1,taup1) + fx*u(i,k,j,2,taup1))*detmr
            t2 = (u(i,k,j,2,taup1) - fx*u(i,k,j,1,taup1))*detmr
            u(i,k,j,1,taup1) = u(i,k,j,1,taum1) + c2dtuv*t1
            u(i,k,j,2,taup1) = u(i,k,j,2,taum1) + c2dtuv*t2
          enddo
        enddo
      enddo
# else
      do n=1,2
        do j=js,je
          do k=1,km
            do i=istrt,iend
              u(i,k,j,n,taup1) = u(i,k,j,n,taum1)
     &                            + c2dtuv*u(i,k,j,n,taup1)
            enddo
          enddo
        enddo
      enddo
# endif

!-----------------------------------------------------------------------
!     subtract incorrect vertical means (related to ignoring horizontal
!     gradients of the surface pressure) to get pure internal modes.
!-----------------------------------------------------------------------

      do n=1,2
        do j=js,je
          do i=istrt,iend
            baru(i,j,n) = c0
          enddo
        enddo
        do j=js,je
          do k=1,km
            do i=istrt,iend
              baru(i,j,n) = baru(i,j,n) + u(i,k,j,n,taup1)*dzt(k)
            enddo
          enddo
        enddo
        do j=js,je
          jrow  = j + joff
          do i=istrt,iend
            baru(i,j,n) = baru(i,j,n)*hr(i,jrow)
          enddo
        enddo
        do j=js,je
          do k=1,km
            do i=istrt,iend
              u(i,k,j,n,taup1) = u(i,k,j,n,taup1)
     &                          - umask(i,k,j)*baru(i,j,n)
            enddo
          enddo
          call setbcx (u(1,1,j,n,taup1), imt, km)
        enddo
      enddo

!-----------------------------------------------------------------------
!     construct diagnostics involving internal mode velocity at "tau+1"
!-----------------------------------------------------------------------

      call diagc2 (joff, js, je, is, ie)

# if defined O_fourfil || defined O_firfil

!-----------------------------------------------------------------------
!     filter velocity components at high latitudes
!-----------------------------------------------------------------------

      if (istrt .eq. 2 .and. iend .eq. imt-1) then
        call filuv (joff, js, je)
      else
        write (stdout,'(a)')
     &  'Error: filtering requires is=2 and ie=imt-1 in clinic'
        stop '=>clinic'
      endif
# endif
      do j=js,je
        call setbcx (u(1,1,j,1,taup1), imt, km)
        call setbcx (u(1,1,j,2,taup1), imt, km)
      enddo

# if defined O_ice_evp

!-----------------------------------------------------------------------
!     if needed, construct the Ice S.B.C.(surface boundary conditions)
!     averaged over this segment
!-----------------------------------------------------------------------

      call isbcu (joff, js, je, istrt, iend, igu, igv)
      call asbcu (joff, js, je, istrt, iend, isu, isv)
# endif

      return
      end

      subroutine diagc1 (joff, js, je, is, ie, n)

!-----------------------------------------------------------------------
!     construct diagnostics which don`t require internal mode velocity
!     at "tau+1" for each velocity component "n"

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!       n    = (1,2) = (u,v) velocity component
!-----------------------------------------------------------------------

      implicit none

      integer n, j, js, je, jrow, joff, k, i, is, ie

      real dudx, ce, dudy, cn, dudz, cb, fx, weight

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "diag.h"
      include "diaga.h"
      include "grdvar.h"
      include "hmixc.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"

      real temp(imt,km)

# if defined O_save_mixing_coeff

!-----------------------------------------------------------------------
!     diagnostic: estimate mixing coefficients on east, north, and
!                 bottom face of U cells from the flux
!-----------------------------------------------------------------------

      if (cmixts .and. n .eq. 1 .and. eots) then
        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=2,imt-1
              dudx = (u(i+1,k,j,1,taum1)-u(i,k,j,1,taum1))
     &                *csur(jrow)*dxtr(i+1) + epsln
              ce(i,k,j,1) = diff_fe(i,k,j)/dudx
#  if !defined O_consthmix || defined O_biharmonic
              dudy = (u(i,k,j+1,1,taum1)-u(i,k,j,1,taum1))
     &                *dytr(jrow+1) + epsln
              cn(i,k,j,1) = cstr(jrow+1)*diff_fn(i,k,j)/dudy
#  else
              cn(i,k,j,1) = cstr(jrow+1)*am
#  endif
            enddo
          enddo
        enddo
        do j=js,je
          jrow = j + joff
          do k=1,km-1
            do i=2,imt-1
              dudz = (u(i,k,j,1,taum1)-u(i,k+1,j,1,taum1))
     &                *dzwr(k) + epsln
              cb(i,k,j,1) = diff_fb(i,k,j)/dudz
            enddo
          enddo
          do i=2,imt-1
            cb(i,km,j,1) = 0.0
          enddo
        enddo

        do j=js,je
          call setbcx (ce(1,1,j,1), imt, km)
          call setbcx (cn(1,1,j,1), imt, km)
          call setbcx (cb(1,1,j,1), imt, km)
        enddo
      endif
# endif

# if defined O_time_step_monitor

!-----------------------------------------------------------------------
!     diagnostic: accumulate global kinetic energy on "tau" velocity
!-----------------------------------------------------------------------

      if (tsiperts .and. eots) then
        do j=js,je
          jrow = j + joff
          fx = rho0*p5*csu(jrow)*dyu(jrow)
#  if defined O_symmetry
          if (jrow .eq. jmtm1) fx = fx*p5
#  endif
          do k=1,km
            do i=is,ie
              weight    = fx*dzt(k)*dxu(i)
              temp(i,k) = u(i,k,j,n,tau)**2*weight
            enddo
            do i=is,ie
              ektot(k,jrow) = ektot(k,jrow) + temp(i,k)
            enddo
          enddo
        enddo
      endif
# endif

# if defined O_energy_analysis

!-----------------------------------------------------------------------
!     diagnostic: integrate work done by the r.h.s. terms in the
!                  momentum equations.
!-----------------------------------------------------------------------

      if (glents .and. eots) call ge1 (joff, js, je, is, ie, n)
# endif

# if defined O_term_balances

!-----------------------------------------------------------------------
!     diagnostic: integrate r.h.s. terms in the momentum equations
!                 over specified regional volumes
!-----------------------------------------------------------------------

      if (trmbts .and. eots) call utb1 (joff, js, je, is, ie, n)
# endif
# if defined O_xbts

!-----------------------------------------------------------------------
!     diagnostic: accumulate r.h.s terms in the momentum equation
!-----------------------------------------------------------------------

      if (xbtperts .and. eots) call uxbt1 (joff, js, je, n)
# endif
      return
      end

      subroutine diagc2 (joff, js, je, is, ie)

!-----------------------------------------------------------------------
!     construct diagnostics requiring internal mode velocity at "tau+1"
!     and those not dependent on velocity component fluxes.

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!-----------------------------------------------------------------------

      implicit none

      integer joff, js, je, is, ie

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "levind.h"
      include "scalar.h"
      include "switch.h"

# if defined O_energy_analysis

!-----------------------------------------------------------------------
!     diagnostic: integrate work done by du/dt in the momentum equations
!                 the external mode part of "u" at "tau+1" will be
!                 accounted for after the external mode is solved.
!                 also, integrate the work done by buoyancy.
!-----------------------------------------------------------------------

      if (glents .and. eots) then
        call ge2 (joff, js, je, is, ie, kmt, kmu, c2dtuv, grav, rho0r)
      endif
# endif

# if defined O_term_balances

!-----------------------------------------------------------------------
!     diagnostic: add du/dt and implicit coriolis terms to the integrals
!                 over specified volumes. the external mode parts will
!                 be accounted for after the external mode is solved.
!-----------------------------------------------------------------------

      if (trmbts .and. eots) then
        call utb2 (joff, js, je, is, ie, c2dtuv, acor)
      endif
# endif

# if defined O_xbts

!-----------------------------------------------------------------------
!     diagnostic: accumulate du/dt and implicit coriolis terms from the
!                 momentum equations
!-----------------------------------------------------------------------

      if (xbtperts .and. eots) call uxbt2 (joff, js, je, c2dtuv, acor)
# endif
      return
      end

      subroutine asbcu (joff, js, je, is, ie, iu, iv)

!-----------------------------------------------------------------------
!     construct the Atmos S.B.C.(surface boundary conditions)

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!       iu   = index for u component
!       iv   = index for v component

!     reference: Pacanowski, R.C., Effect of Equatorial Currents
!                on Surface Stress (JPO, Vol 17, No. 6, June 1987)
!-----------------------------------------------------------------------

      implicit none

      integer iu, iv, j, js, je, jrow, joff, i, is, ie

      real rts

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"

!     initialize S.B.C. at the beginning of each ocean segment
!     (do not alter values in land)

      if (eots .and. osegs .and. iu .ne. 0 .and. iv .ne. 0) then
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0) then
              sbc(i,jrow,iu) = c0
              sbc(i,jrow,iv) = c0
            endif
          enddo
        enddo
      endif

!     accumulate surface currents for the Atmos S.B.C. every time step

      if (eots .and. iu .ne. 0 .and. iv .ne. 0) then
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            sbc(i,jrow,iu) = sbc(i,jrow,iu) + p25*(
     &                       u(i,1,j,1,tau) + u(i-1,1,j,1,tau)
     &                     + u(i,1,j-1,1,tau) + u(i-1,1,j-1,1,tau))
            sbc(i,jrow,iv) = sbc(i,jrow,iv) + p25*(
     &                       u(i,1,j,2,tau) + u(i-1,1,j,2,tau)
     &                     + u(i,1,j-1,2,tau) + u(i-1,1,j-1,2,tau))
          enddo
        enddo
      endif

!     average the surface currents for the Atmos S.B.C. at the end
!     of each ocean segment. (do not alter values in land)

      if (eots .and. osege .and. iu .ne. 0 .and. iv .ne. 0) then
        rts = c1/ntspos
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0) then
              sbc(i,jrow,iu) = rts*sbc(i,jrow,iu)
              sbc(i,jrow,iv) = rts*sbc(i,jrow,iv)
            endif
          enddo
        enddo
      endif
#endif

      return
      end

#if defined O_mom && defined O_ice_evp
      subroutine isbcu (joff, js, je, is, ie, iu, iv)

!-----------------------------------------------------------------------
!     construct the Ice S.B.C.(surface boundary conditions)

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!       iu   = index for u component
!       iv   = index for v component

!-----------------------------------------------------------------------

      implicit none

      integer iu, iv, j, js, je, jrow, joff, i, is, ie

      real rts

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"

!     initialize S.B.C. at the beginning of each ocean segment
!     (do not alter values in land)

      if (eots .and. osegs .and. iu .ne. 0 .and. iv .ne. 0) then
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0) then
              sbc(i,jrow,iu) = c0
              sbc(i,jrow,iv) = c0
            endif
          enddo
        enddo
      endif

!     accumulate geostrophic currents (level 2) for the Ice S.B.C.
!     every time step

      if (eots .and. iu .ne. 0 .and. iv .ne. 0) then
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            sbc(i,jrow,iu) = sbc(i,jrow,iu) + u(i,2,j,1,tau)
            sbc(i,jrow,iv) = sbc(i,jrow,iv) + u(i,2,j,2,tau)
          enddo
        enddo
      endif

!     average the currents for the Ice S.B.C. at the end of each ocean
!     segment. (do not alter values in land)

      if (eots .and. osege .and. iu .ne. 0 .and. iv .ne. 0) then
        rts = c1/ntspos
        do j=js,je
          jrow  = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0) then
              sbc(i,jrow,iu) = rts*sbc(i,jrow,iu)
              sbc(i,jrow,iv) = rts*sbc(i,jrow,iv)
            endif
          enddo
        enddo
      endif

      return
      end
#endif

#if defined O_implicitvmix
      subroutine ivdifu (joff, js, je, is, ie, n)

!-----------------------------------------------------------------------
!     solve vertical diffusion of velocity implicitly

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       n     = velocity component
!       twodt = (2*dtuv, dtuv) on (leapfrog, mixing) time steps
!-----------------------------------------------------------------------

      implicit none

      integer j, js, je, k, i, is, ie, n, joff

      real c1, up1, r2dtuv, vvca

      include "size.h"
      include "param.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "vmixc.h"

      real twodt(km)

!     set some constants

      c1 = 1.0

!     construct the "tau+1" velocity without implicit vertical diffusion

      do j=js,je
        do k=1,km
          do i=is,ie
            u(i,k,j,n,taup1) = u(i,k,j,n,taum1)+c2dtuv*u(i,k,j,n,taup1)
          enddo
        enddo
      enddo
# if defined O_xbts || defined O_energy_analysis || defined O_term_balances

!     store terms to compute implicit vertical diffusion on
!     diagnostic time steps

      if ((xbtperts .or. glents .or. trmbts) .and. eots) then
        do j=js,je
          do k=1,km
            do i=is,ie
              zzi(i,k,j) = u(i,k,j,n,taup1)
            enddo
          enddo
        enddo
      endif
# endif

!     add in the implicit vertical diffusion

      do k=1,km
        twodt(k) = c2dtuv
      enddo
      call invtri (u(1,1,1,n,taup1), smf(1,1,n), bmf(1,1,n)
     &, visc_cbu(1,1,jsmw), twodt, kmu, umask(1,1,1), is, ie
     &, joff, js, je)

      r2dtuv = c1/c2dtuv
# if defined O_xbts || defined O_energy_analysis || defined O_term_balances

!     compute residual implicit vertical diffusion for diagnostics

      if ((xbtperts .or. glents .or. trmbts) .and. eots) then
        do j=js,je
          do k=1,km
            do i=is,ie
              zzi(i,k,j) = r2dtuv*(u(i,k,j,n,taup1) - zzi(i,k,j))
            enddo
          enddo
        enddo
      endif
# endif

!     convert back to time change of velocity

      do j=js,je
        do k=1,km
          do i=is,ie
            u(i,k,j,n,taup1) =r2dtuv*(u(i,k,j,n,taup1)-u(i,k,j,n,taum1))
          enddo
        enddo
      enddo

      return
      end
#endif
