! source file: /Users/oschlies/UVIC/master/source/mom/energy.F
      subroutine ge1 (joff, js, je, is, ie, n)

!-----------------------------------------------------------------------
!     compute global energetics by taking u dot the momentum equations
!     and integrating over the ocean volume

!     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
!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

      include "param.h"
      include "coord.h"
      include "emode.h"
      include "diag.h"
      include "grdvar.h"
      include "hmixc.h"
      include "levind.h"
      include "mw.h"
      include "vmixc.h"
      dimension ext(imt,2)
      include "fdifm.h"

      do j=js,je
        jrow = j + joff

!-----------------------------------------------------------------------
!       compute the external mode velocities
!-----------------------------------------------------------------------

        do i=is,ie
          diag1    = psi(i+1,jrow+1,1) - psi(i  ,jrow,1)
          diag0    = psi(i  ,jrow+1,1) - psi(i+1,jrow,1)
          ext(i,1) = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow)
          ext(i,2) =  (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow)
        enddo

!-----------------------------------------------------------------------
!       compute work done by each term on the internal and external
!       modes in the momentum equations
!-----------------------------------------------------------------------

        fx = csu(jrow)*dyu(jrow)

        do k=1,km
          do i=is,ie
            uext   = ext(i,n)
            uint   = u(i,k,j,n,tau) - uext
            boxvol = fx*dxu(i)*dzt(k)

!-----------------------------------------------------------------------
!           pressure term
!-----------------------------------------------------------------------

            term = -umask(i,k,j)*grad_p(i,k,j,n)
            call addto (engint(k,6,jrow), uint*term*boxvol)
            call addto (engext(6,jrow),   uext*term*boxvol)

!-----------------------------------------------------------------------
!           coriolis term does no work: u*fv - v*fu = 0
!           implicit coriolis work will be reflected in the imbalance
!           between horizontal pressure forces and buoyancy
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!           zonal and meridional advection of momentum + metric term
!-----------------------------------------------------------------------

            term =-umask(i,k,j)*(ADV_Ux(i,k,j) + ADV_Uy(i,k,j,jrow,n)
     &           + ADV_metric(i,k,j,jrow,n))
            call addto (engint(k,2,jrow), uint*term*boxvol)
            call addto (engext(2,jrow),   uext*term*boxvol)

!-----------------------------------------------------------------------
!           vertical advection of momentum
!-----------------------------------------------------------------------

            term = -umask(i,k,j)*ADV_Uz(i,k,j)
            call addto (engint(k,3,jrow), uint*term*boxvol)
            call addto (engext(3,jrow),   uext*term*boxvol)

!-----------------------------------------------------------------------
!           zonal and meridional diffusion of momentum + metric term
!-----------------------------------------------------------------------

            term = umask(i,k,j)*(DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j,jrow,n)
     &           + DIFF_metric(i,k,j,jrow,n))
            call addto (engint(k,4,jrow), uint*term*boxvol)
            call addto (engext(4,jrow),   uext*term*boxvol)

!-----------------------------------------------------------------------
!           vertical diffusion of momentum
!-----------------------------------------------------------------------

            term = umask(i,k,j)*DIFF_Uz(i,k,j)
            call addto (engint(k,5,jrow), uint*term*boxvol)
            call addto (engext(5,jrow),   uext*term*boxvol)
          enddo
        enddo

!-----------------------------------------------------------------------
!       work done by wind stress
!-----------------------------------------------------------------------

        k = 1
        do i=is,ie
          uext = ext(i,n)
          uint = u(i,k,j,n,tau) - uext
          term = umask(i,k,j)*smf(i,j,n)
          call addto (engint(k,7,jrow), uint*term*fx*dxu(i))
          call addto (engext(7,jrow),   uext*term*fx*dxu(i))
        enddo

!-----------------------------------------------------------------------
!       work done by bottom drag
!-----------------------------------------------------------------------

        do i=is,ie
          uext = ext(i,n)
          k = kmu(i,jrow)
          if (k .ne. 0) then
            uint = u(i,k,j,n,tau) - uext
            term = -umask(i,k,j)*bmf(i,j,n)
            call addto (engint(k,8,jrow), uint*term*fx*dxu(i))
            call addto (engext(8,jrow),   uext*term*fx*dxu(i))
          endif
        enddo
      enddo

      return
      end

      subroutine ge2 (joff, js, je, is, ie, kmt, kmu, c2dtuv, grav
     &,               rho0r)

!-----------------------------------------------------------------------
!     compute global energetics by taking u dot the momentum equations
!     and integrating over the entire ocean volume

!     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
!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

      include "param.h"
      include "coord.h"
      include "emode.h"
      include "diag.h"
      include "grdvar.h"
      include "mw.h"
      dimension ext(imt,2)
      dimension kmt(imt,jmt), kmu(imt,jmt)

      do j=js,je
        jrow  = j + joff

!-----------------------------------------------------------------------
!       set local constants
!-----------------------------------------------------------------------

        fx = csu(jrow)*dyu(jrow)

!-----------------------------------------------------------------------
!       compute the external mode velocities
!-----------------------------------------------------------------------

        do i=is,ie
          diag1    = psi(i+1,jrow+1,1) - psi(i  ,jrow,1)
          diag0    = psi(i  ,jrow+1,1) - psi(i+1,jrow,1)
          ext(i,1) = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow)
          ext(i,2) =  (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow)
        enddo

!-----------------------------------------------------------------------
!       compute work done by each term on the internal and external
!       modes in the momentum equations
!-----------------------------------------------------------------------

        r2dt = c1/c2dtuv
        do n=1,2
          do k=1,km
            do i=is,ie
              uext = ext(i,n)
              uint   = u(i,k,j,n,tau) - uext
              boxvol = fx*dxu(i)*dzt(k)

!-----------------------------------------------------------------------
!             total change in kinetic energy.
!-----------------------------------------------------------------------

              term = umask(i,k,j)*(u(i,k,j,n,taup1)-
     &                             u(i,k,j,n,taum1))*r2dt
              call addto (engint(k,1,jrow), uint*term*boxvol)
               call addto (engext(1,jrow),   uext*term*boxvol)
            enddo
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     compute the work done by buoyancy integrated over the entire
!     ocean volume.
!-----------------------------------------------------------------------

      do j=js,je
        jrow  = j + joff
        f1 = cst(jrow)*dyt(jrow)
        do i=is,ie
          kz = kmt(i,jrow)
          if (kz .ne. 0) then
            area = f1*dxt(i)
            fx   = area*grav*rho0r*p5
            term =  - c2*fx*dzw(0)*adv_vbt(i,0,j)*rho(i,1,j)
            call addto (buoy(1,jrow), term)
            do k=2,kz
              term =-fx*dzw(k-1)*adv_vbt(i,k-1,j)*
     &               (rho(i,k-1,j) + rho(i,k,j))
              call addto (buoy(k,jrow), term)
            enddo
          endif
        enddo
      enddo

!-----------------------------------------------------------------------
!     find maximum error in continuity for "t" cells and "u" cells
!-----------------------------------------------------------------------

      do j=js,je
        jrow  = j + joff
        do k=1,km
          do i=is,ie
            term =
     &        ((adv_vet(i,k,j) - adv_vet(i-1,k,j))*cstr(jrow)*dxtr(i)
     &        +(adv_vnt(i,k,j) - adv_vnt(i,k,j-1))*cstr(jrow)*dytr(jrow)
     &        +(adv_vbt(i,k-1,j) - adv_vbt(i,k,j))*dztr(k))*tmask(i,k,j)
            if (abs(term) .gt. abs(tcerr(jrow))) then
              tcerr(jrow) = term
              itcerr(jrow) = i
              jtcerr(jrow) = jrow
              ktcerr(jrow) = k
            endif

            term =
     &        ((adv_veu(i,k,j) - adv_veu(i-1,k,j))*csur(jrow)*dxur(i)
     &        +(adv_vnu(i,k,j) - adv_vnu(i,k,j-1))*csur(jrow)*dyur(jrow)
     &        +(adv_vbu(i,k-1,j) - adv_vbu(i,k,j))*dztr(k))*umask(i,k,j)
            if (abs(term) .gt. abs(ucerr(jrow))) then
              ucerr(jrow) = term
              iucerr(jrow) = i
              jucerr(jrow) = jrow
              kucerr(jrow) = k
            endif
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     find max error in "adv_vbt" at bottom and max "adv_vbu" at bottom
!-----------------------------------------------------------------------

      do j=js,je
        jrow  = j + joff
        do i=is,ie
          k = kmt(i,jrow)
          if (k.ne.0 .and.
     &        (abs(adv_vbt(i,k,j)) .gt. abs(wtbot(jrow)))) then
            wtbot(jrow)  = adv_vbt(i,k,j)
            iwtbot(jrow) = i
            jwtbot(jrow) = jrow
            kwtbot(jrow) = k
          endif
        enddo
      enddo

      do j=js,je
        jrow  = j + joff
        do i=is,ie
          k = kmu(i,jrow)
          if (k.ne.0 .and.
     &        (abs(adv_vbu(i,k,j)) .gt. abs(wubot(jrow))))then
            wubot(jrow)  = adv_vbu(i,k,j)
            iwubot(jrow) = i
            jwubot(jrow) = jrow
            kwubot(jrow) = k
          endif
        enddo
      enddo

!-----------------------------------------------------------------------
!     integrate "adv_vbt" and "adv_vbu" for all lat and lon
!-----------------------------------------------------------------------

      do j=js,je
        jrow  = j + joff
        do k=1,km
          do i=is,ie
            udxdy = umask(i,k,j)*dxu(i)*csu(jrow)*dyu(jrow)
            tdxdy = tmask(i,k,j)*dxt(i)*cst(jrow)*dyt(jrow)
            wtlev(k,jrow)  = wtlev(k,jrow) + adv_vbt(i,k,j)*tdxdy
            wulev(k,jrow)  = wulev(k,jrow) + adv_vbu(i,k,j)*udxdy
          enddo
        enddo
      enddo

      return
      end

      subroutine ge3 (c2dtuv)

!-----------------------------------------------------------------------
!     compute global energetics by taking u dot the momentum equations
!     and integrating over the entire ocean volume

!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

      include "param.h"
      include "coord.h"
      include "emode.h"
      include "diag.h"
      include "grdvar.h"
      include "levind.h"
      include "mw.h"
      dimension ext(imt,2,2)

      do jrow=2,jmt-1
        fx = csu(jrow)*dyu(jrow)

!-----------------------------------------------------------------------
!       compute the external mode velocities after external mode has
!       been updated. since the external mode has been updated,
!       m = 1 is "tau+1"
!       m = 2 is "tau"
!-----------------------------------------------------------------------

        do m=1,2
          do i=2,imtm1
            diag1      = psi(i+1,jrow+1,m) - psi(i  ,jrow,m)
            diag0      = psi(i  ,jrow+1,m) - psi(i+1,jrow,m)
            ext(i,1,m) = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow)
            ext(i,2,m) =  (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow)
          enddo
        enddo

!-----------------------------------------------------------------------
!       compute external mode part of work done by "tau+1"
!       component of d/dt. internal mode part is zero. since the
!       external mode has been updated:
!       ubarm1 is "tau"
!       ubar   is "tau+1"
!       note: if using 5 point numerics when solving for the stream
!             function, the total integral is not conserved  for the
!             external mode which shows itself as the "ficticious" term
!             in the global energy integrals.
!-----------------------------------------------------------------------

        r2dt = c1/c2dtuv
        do i=2,imtm1
          uext   = ext(i,1,2)
          vext   = ext(i,2,2)
          uextn  = ext(i,1,1)
          vextn  = ext(i,2,1)
          boxvol = fx*dxu(i)*h(i,jrow)
          term = (uext*uextn + vext*vextn)*r2dt
           call addto (engext(1,jrow),   term*boxvol)
        enddo

      enddo
      return
      end

      subroutine addto (sum, term)
      sum = sum + term
      return
      end