! source file: /Users/oschlies/UVIC/master/source/mom/termbal.F subroutine utb1 (joff, js, je, is, ie, n) !======================================================================= ! accumulate terms in the momentum equations over the ! volume of the specified regions ! 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 "cregin.h" include "diag.h" include "grdvar.h" include "hmixc.h" include "mw.h" include "scalar.h" include "vmixc.h" include "fdifm.h" !----------------------------------------------------------------------- ! set local constants !----------------------------------------------------------------------- do j=js,je jrow = j + joff fx = csu(jrow)*dyu(jrow) !----------------------------------------------------------------------- ! accumulate terms for all regions within the current jrow !----------------------------------------------------------------------- do k=1,km do i=is,ie nreg = nhreg*(mskvr(k)-1) + mskhr(i,jrow) if (nreg .gt. 0 .and. mskhr(i,jrow) .gt. 0) then boxvol = fx*dxu(i)*dzt(k) !----------------------------------------------------------------------- ! pressure term !----------------------------------------------------------------------- term = -umask(i,k,j)*grad_p(i,k,j,n) call addto (termbm(k,2,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! zonal advection (flux form) of momentum !----------------------------------------------------------------------- term = -umask(i,k,j)*ADV_Ux(i,k,j) call addto (termbm(k,3,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! pure zonal advection of momentum !----------------------------------------------------------------------- ! - U(U)x = U(U)x - (UU)x (when n=1) ! - U(V)x = V(U)x - (UV)x (when n=2) dudx = (adv_veu(i,k,j)-adv_veu(i-1,k,j))*dxur(i) & *csur(jrow) term = umask(i,k,j)*(u(i,k,j,n,tau)*dudx - ADV_Ux(i,k,j)) call addto (termbm(k,14,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! advective metric term !----------------------------------------------------------------------- term = ADV_metric(i,k,j,jrow,n) call addto (termbm(k,13,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! meridional advection (flux form) of momentum !----------------------------------------------------------------------- term = -umask(i,k,j)*ADV_Uy(i,k,j,jrow,n) call addto (termbm(k,4,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! pure meridional advection of momentum !----------------------------------------------------------------------- ! - V(U)y = U(V)y - (VU)y (when n=1) ! - V(V)y = V(V)y - (VV)y (when n=2) dvdy = (adv_vnu(i,k,j)-adv_vnu(i,k,j-1))*dyur(jrow) & *csur(jrow) term = umask(i,k,j)*(u(i,k,j,n,tau)*dvdy & - ADV_Uy(i,k,j,jrow,n)) call addto (termbm(k,15,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! vertical advection (flux form) of momentum !----------------------------------------------------------------------- term = -umask(i,k,j)*ADV_Uz(i,k,j) call addto (termbm(k,5,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! pure vertical advection of momentum !----------------------------------------------------------------------- ! - W(U)z = U(W)z - (WU)z (when n=1) ! - W(V)z = V(W)z - (WV)z (when n=2) dwdz = (adv_vbu(i,k-1,j)-adv_vbu(i,k,j))*dztr(k) term = umask(i,k,j)*(u(i,k,j,n,tau)*dwdz - ADV_Uz(i,k,j)) call addto (termbm(k,16,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! zonal diffusion of momentum !----------------------------------------------------------------------- term = umask(i,k,j)*DIFF_Ux(i,k,j) call addto (termbm(k,6,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! meridional diffusion of momentum !----------------------------------------------------------------------- term = umask(i,k,j)*DIFF_Uy(i,k,j,jrow,n) call addto (termbm(k,7,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! diffusive metric term !----------------------------------------------------------------------- term = umask(i,k,j)*DIFF_metric(i,k,j,jrow,n) call addto (termbm(k,9,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! vertical diffusion of momentum !----------------------------------------------------------------------- term = umask(i,k,j)*DIFF_Uz(i,k,j) call addto (termbm(k,8,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! coriolis term !----------------------------------------------------------------------- term = umask(i,k,j)*CORIOLIS(i,k,j,jrow,n) call addto (termbm(k,10,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! accumulate the source term !----------------------------------------------------------------------- term = umask(i,k,j)*source(i,k,j) call addto (termbm(k,11,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! accumulate u, v, and w !----------------------------------------------------------------------- term = umask(i,k,j)*u(i,k,j,n,tau) call addto (termbm(k,17,n,nreg), term*boxvol) if (n .eq. 2) then term = p5*(adv_vbu(i,k,j)+adv_vbu(i,k-1,j))*umask(i,k,j) call addto (avgw(nreg), term*boxvol) endif !----------------------------------------------------------------------- ! accumulate the surface momentum flux !----------------------------------------------------------------------- if (k .eq. 1) then term = umask(i,k,j)*smf(i,j,n) call addto (smflx(n,nreg), term*fx*dxu(i)) endif endif enddo enddo enddo return end subroutine utb2 (joff, js, je, is, ie, c2dtuv, acor) !======================================================================= ! accumulate external mode parts of d/dt and the implicit coriolis ! term in the momentum equations over the volume in the specified ! regions ! 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 ! c2dtuv = (2*dtuv,dtuv) on (lpfrod,mixing) time steps ! acor = implicit factor ! based on code by: R. C. Pacanowski !======================================================================= include "param.h" include "coord.h" include "cregin.h" include "diag.h" include "grdvar.h" include "mw.h" !----------------------------------------------------------------------- ! local constants !----------------------------------------------------------------------- r2dt = c1/c2dtuv do j=js,je jrow = j + joff fx = csu(jrow)*dyu(jrow) do n=1,2 do k=1,km do i=is,ie nreg = nhreg*(mskvr(k)-1) + mskhr(i,jrow) if (nreg .gt. 0 .and. mskhr(i,jrow) .gt. 0) then boxvol = fx*dxu(i)*dzt(k) !----------------------------------------------------------------------- ! d/dt of velocity (external mode part of tau+1 will be ! added later when the external mode is solved) !----------------------------------------------------------------------- term = umask(i,k,j)*(u(i,k,j,n,taup1) - & u(i,k,j,n,taum1))*r2dt call addto (termbm(k,1,n,nreg), term*boxvol) !----------------------------------------------------------------------- ! implicit coriolis term (external mode part will be added ! later when external mode is solved) !----------------------------------------------------------------------- if (acor .ne. c0) then term = umask(i,k,j)*acor*cori(i,jrow,n)* & (u(i,k,j,3-n,taup1) - u(i,k,j,3-n,taum1)) call addto (termbm(k,10,n,nreg), term*boxvol) endif endif enddo enddo enddo enddo return end subroutine utb3 !======================================================================= ! accumulate external mode parts of d/dt, the implicit coriolis ! term and the surface pressure gradientsover the volume in the ! specified regions. ! based on code by: R. C. Pacanowski !======================================================================= include "param.h" include "coord.h" include "cregin.h" include "grdvar.h" include "levind.h" include "mw.h" include "scalar.h" include "termbal.h" parameter (is=1, ie=1, js=1, je=1) dimension psgrad(is:ie,js:je,2) do jrow=1,jmt-1 fddt = csu(jrow)*dyu(jrow)/c2dtuv fspr = csu(jrow)*dyu(jrow) do i=2,imtm1 atosp = acor*cori(i,jrow,1) f1 = atosp*csu(jrow)*dyu(jrow) kz = kmu(i,jrow) if (kz .ne. 0) then do k=1,kz n = nhreg*(mskvr(k)-1) + mskhr(i,jrow) if (n .gt. 0 .and. mskhr(i,jrow) .gt. 0) then ! construct the surface pressure gradients for pt (i,jrow) if (k .eq. 1) then call calc_psgrad(psgrad, uext, vext, jrow, jrow, i, i) endif boxfac = fddt*dxu(i)*dzt(k) boxspr = fspr*dxu(i)*dzt(k) termbm(k,1,1,n) = termbm(k,1,1,n) + uext*boxfac termbm(k,1,2,n) = termbm(k,1,2,n) + vext*boxfac termbm(k,12,1,n) = termbm(k,12,1,n) - & psgrad(is,js,1)*boxspr termbm(k,12,2,n) = termbm(k,12,2,n) - & psgrad(is,js,2)*boxspr boxacr = f1*dxu(i)*dzt(k) termbm(k,10,1,n) = termbm(k,10,1,n) + vext*boxacr termbm(k,10,2,n) = termbm(k,10,2,n) - uext*boxacr endif enddo endif enddo enddo return end subroutine ttb1 (joff, js, je, is, ie, n) !======================================================================= ! accumulate terms in the tracer equations over the volume in the ! specified regions ! 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 = tracer component ! based on code by: R. C. Pacanowski !======================================================================= include "param.h" include "accel.h" include "coord.h" include "cregin.h" include "diag.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "mw.h" include "scalar.h" include "vmixc.h" include "fdift.h" !----------------------------------------------------------------------- ! limit the longitude indices !----------------------------------------------------------------------- istrt = max(2,is) iend = min(imt-1,ie) do j=js,je jrow = j + joff !----------------------------------------------------------------------- ! set local constants !----------------------------------------------------------------------- fx = cst(jrow)*dyt(jrow) do k=1,km do i=istrt,iend nreg = nhreg*(mskvr(k)-1) + mskhr(i,jrow) if (nreg .gt. 0 .and. mskhr(i,jrow) .gt. 0) then area = fx*dxt(i) boxvol = area*dzt(k) !----------------------------------------------------------------------- ! tracer !----------------------------------------------------------------------- term = tmask(i,k,j)*t(i,k,j,n,tau) call addto (termbt(k,15,n,nreg), term*boxvol) call addto (termbt(k,15,n,0), term*boxvol) !----------------------------------------------------------------------- ! d(tracer)/dt !----------------------------------------------------------------------- r2dt = c1/(c2dtts*dtxcel(k)) term = tmask(i,k,j)*(t(i,k,j,n,taup1) - & t(i,k,j,n,taum1))*r2dt call addto (termbt(k,9,n,nreg), term*boxvol) call addto (termbt(k,9,n,0), term*boxvol) !----------------------------------------------------------------------- ! zonal advection (flux form) of tracer !----------------------------------------------------------------------- term = -tmask(i,k,j)*ADV_Tx(i,k,j) call addto (termbt(k,2,n,nreg), term*boxvol) call addto (termbt(k,2,n,0), term*boxvol) !----------------------------------------------------------------------- ! pure zonal advection of tracer !----------------------------------------------------------------------- ! - U(T)x = T(U)x - (UT)x dudx = (adv_vet(i,k,j)-adv_vet(i-1,k,j))*dxtr(i) & *cstr(jrow) term = tmask(i,k,j)*(t(i,k,j,n,tau)*dudx - ADV_Tx(i,k,j)) call addto (termbt(k,11,n,nreg), term*boxvol) call addto (termbt(k,11,n,0), term*boxvol) !----------------------------------------------------------------------- ! meridional advection (flux form) of tracer !----------------------------------------------------------------------- term = -tmask(i,k,j)*ADV_Ty(i,k,j,jrow,n) call addto (termbt(k,3,n,nreg), term*boxvol) call addto (termbt(k,3,n,0), term*boxvol) !----------------------------------------------------------------------- ! pure meridional advection of tracer !----------------------------------------------------------------------- ! - V(T)y = T(V)y - (VT)y dvdy = (adv_vnt(i,k,j)-adv_vnt(i,k,j-1))*dytr(jrow) & *cstr(jrow) term = tmask(i,k,j)*(t(i,k,j,n,tau)*dvdy & - ADV_Ty(i,k,j,jrow,n)) call addto (termbt(k,12,n,nreg), term*boxvol) call addto (termbt(k,12,n,0), term*boxvol) !----------------------------------------------------------------------- ! vertical advection (flux form) of tracer !----------------------------------------------------------------------- term = -tmask(i,k,j)*ADV_Tz(i,k,j) call addto (termbt(k,4,n,nreg), term*boxvol) call addto (termbt(k,4,n,0), term*boxvol) !----------------------------------------------------------------------- ! pure vertical advection of tracer !----------------------------------------------------------------------- ! - W(T)z = T(W)z - (WT)z dwdz = (adv_vbt(i,k-1,j)-adv_vbt(i,k,j))*dztr(k) term = tmask(i,k,j)*(t(i,k,j,n,tau)*dwdz - ADV_Tz(i,k,j)) call addto (termbt(k,13,n,nreg), term*boxvol) call addto (termbt(k,13,n,0), term*boxvol) !----------------------------------------------------------------------- ! zonal diffusion of tracer !----------------------------------------------------------------------- term = tmask(i,k,j)*DIFF_Tx(i,k,j) call addto (termbt(k,5,n,nreg), term*boxvol) call addto (termbt(k,5,n,0), term*boxvol) !----------------------------------------------------------------------- ! meridional diffusion of tracer !----------------------------------------------------------------------- term = tmask(i,k,j)*DIFF_Ty(i,k,j,jrow,n) call addto (termbt(k,6,n,nreg), term*boxvol) call addto (termbt(k,6,n,0), term*boxvol) !----------------------------------------------------------------------- ! vertical diffusion of tracer !----------------------------------------------------------------------- term = tmask(i,k,j)*DIFF_Tz(i,k,j) & + tmask(i,k,j)*zzi(i,k,j) call addto (termbt(k,7,n,nreg), term*boxvol) call addto (termbt(k,7,n,0), term*boxvol) !----------------------------------------------------------------------- ! tracer source term !----------------------------------------------------------------------- term = tmask(i,k,j)*source(i,k,j) call addto (termbt(k,8,n,nreg), term*boxvol) call addto (termbt(k,8,n,0), term*boxvol) if (k .eq. 1) then !----------------------------------------------------------------------- ! surface tracer !----------------------------------------------------------------------- term = tmask(i,k,j)*t(i,k,j,n,tau) call addto (asst(n,nreg), term*area) call addto (asst(n,0), term*area) !----------------------------------------------------------------------- ! surface tracer flux !----------------------------------------------------------------------- term = tmask(i,k,j)*stf(i,j,n) call addto (stflx(n,nreg), term*area) call addto (stflx(n,0), term*area) endif endif enddo enddo enddo return end subroutine ttb2 (joff, js, je, is, ie, iterm) !======================================================================= ! accumulate d/dt and change in tracer in the tracer equations over ! the volume in the specified regions ! 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 ! iterm = 1 => total change ! iterm = 10 => change due to filtering ! based on code by: R. C. Pacanowski !======================================================================= include "param.h" include "accel.h" include "coord.h" include "cregin.h" include "diag.h" include "grdvar.h" include "mw.h" include "scalar.h" if (iterm .ne. 1 .and. iterm .ne. 10) then write (stdout,*) '=>Error: iterm=',iterm,' in ttb2' stop '=>ttb2' endif do j=js,je jrow = j + joff fx = cst(jrow)*dyt(jrow) do n=1,nt do k=1,km r2dt = c1/(c2dtts*dtxcel(k)) do i=is,ie nreg = nhreg*(mskvr(k)-1) + mskhr(i,jrow) if (nreg .gt. 0 .and. mskhr(i,jrow) .gt. 0) then area = fx*dxt(i) boxvol = area*dzt(k) !----------------------------------------------------------------------- ! d/dt(tracer) !----------------------------------------------------------------------- term = tmask(i,k,j)*(t(i,k,j,n,taup1) - & t(i,k,j,n,taum1))*r2dt call addto (termbt(k,iterm,n,nreg), term*boxvol) call addto (termbt(k,iterm,n,0), term*boxvol) !----------------------------------------------------------------------- ! change in variance of tracer !----------------------------------------------------------------------- if (iterm .eq. 1) then term = tmask(i,k,j)*(t(i,k,j,n,taup1)**2- & t(i,k,j,n,taum1)**2) call addto (termbt(k,14,n,nreg), term*boxvol) call addto (termbt(k,14,n,0), term*boxvol) endif endif enddo enddo enddo enddo return end