! source file: /Users/jfidler/work/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