! source file: /sfs/fs6/home-geomar/smomw258/UVic_ESCM/2.9/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
!-----------------------------------------------------------------------

      implicit none

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

      real adv_ux, adv_uy, adv_uz, adv_metric, diff_ux, diff_uz
      real diff_uy, diff_metric, coriolis, acor
      real diag1, diag0, fx, uext, uint, boxvol, term, dhdt, gcor

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.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"
      include "fdifm.h"

      real ext(imt,2)

      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
!-----------------------------------------------------------------------

      implicit none

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

      real fx, diag1, diag0, r2dt, c2dtuv, uext, uint, boxvol, term
      real f1, area, grav, rho0r, udxdy, tdxdy

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "emode.h"
      include "diag.h"
      include "grdvar.h"
      include "mw.h"

      integer kmt(imt,jmt), kmu(imt,jmt)

      real ext(imt,2)

      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
!-----------------------------------------------------------------------

      implicit none

      integer i, jrow, m

      real  fx, diag1, diag0, r2dt, c2dtuv, uext, vext, uextn
      real vextn, boxvol, term

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "emode.h"
      include "diag.h"
      include "grdvar.h"
      include "levind.h"
      include "mw.h"

      real  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)

      implicit none

      real sum, term

      sum = sum + term

      return
      end