subroutine isopi (error, am, ah)
#if defined O_mom
# if defined O_isopycmix
!=======================================================================
! Initialization for isopycnal mixing scheme
! -Disopycmix gives either the full or small angle isopycnal
! mixing tensor
! -Disopycmix -Dgent_mcwilliams gives the isopycnal mixing tensor
! plus the advective velocity parameterization of Gent_McWilliams
! input:
! error = logical to signal problems
! am = background mixing coeff for momentum
! ah = horizontal mixing coeff for tracers
! slmx = max slope of isopycnals
! ahisop = isopycnal tracer diffusivity(cm**2/sec)
! athkdf = isopycnal thickness diffusivity (cm**2/sec)
! output:
! The above input can be reset via namelist
!=======================================================================
implicit none
character(120) :: fname, new_file_name
integer i, k, j, ip, kr, jq, io, i_delta, j_delta, k_delta, jrow
integer iou, n, ib(10), ic(10)
logical exists, inqvardef, error
real slmx, s_dm, ah, am, ft, delta1, delta2
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "coord.h"
include "grdvar.h"
include "iounit.h"
include "isopyc.h"
include "levind.h"
include "mw.h"
include "scalar.h"
include "switch.h"
include "vmixc.h"
real tmpijk(imtm2,jmtm2,km)
# if !defined O_gent_mcwilliams
real athkdf
# endif
namelist /isopyc/ slmx, ahisop, athkdf, del_dm, s_dm
write (stdout,'(/,20x,a,/,20x,a,/)')
& 'I S O P Y C M I X I N I T I A L I Z A T I O N'
# if defined O_gent_mcwilliams
&,'Isopycnal mixing tensor + GM thickness diffusion'
# else
&,'Isopycnal mixing tensor.'
# endif
!-----------------------------------------------------------------------
! USER INPUT
! initialize variables (all mixing units are cm**2/sec.)
!-----------------------------------------------------------------------
! maximum isopycnal slope
slmx = 1.0/100.0
! isopycnal diffusion coefficient
ahisop = 1.e7
! isopycnal thickness diffusion coefficient
athkdf = 1.0e7
! transition for scaling diffusion coefficient
del_dm = 4.0/1000.0
! half width scaling for diffusion coefficient
s_dm = 1.0/1000.0
!-----------------------------------------------------------------------
! provide for namelist over-ride of above settings + documentation
!-----------------------------------------------------------------------
call getunit (io, 'control.in'
&, 'formatted sequential rewind')
read (io,isopyc,end=100)
100 continue
write (stdout,isopyc)
call relunit (io)
! set reciprocal of maximum isopycnal slope and scaling half width
slmxr = c1/slmx
s_dmr = c1/s_dm
!-----------------------------------------------------------------------
! check for problems
!-----------------------------------------------------------------------
# if defined O_gent_mcwilliams && !defined O_isopycmix
write (stdout,*)
& '==> Error: "isopycmix" must be enabled with "gent_mcwilliams"'
error = .true.
# endif
# if defined O_full_tensor
write (stdout,*)
& '==> Note: full isopycnal tensor is being used. '
# else
write (stdout,*)
& '==> Note: small angle approximation to isopycnal tensor is'
&,' being used. '
# if defined O_dm_taper
write (stdout,*)
& '==> Note: Danabasoglu-McWilliams hyperbolic tangent taper '
&,' is being used to reduce mixing coeff within '
&,' areas of steeply sloping isopycnals '
# else
write (stdout,*)
& '==> Note: Re-scaling of Gerdes et al is being used to '
&,' reduce mixing coeff within areas of steeply '
&,' sloping isopycnals '
# endif
# endif
# if !defined O_constvmix
write (stdout,*)
& '==> Error: "isopycmix" only works with "constvmix" for now '
error = .true.
# endif
!-----------------------------------------------------------------------
! print out tracer diffusion coefficients in isopycnal case
!-----------------------------------------------------------------------
write(stdout,9102) ahisop, athkdf, ah, am
9102 format(/' ahisop = ',e12.6,' along isopyncal tracer mixing '
&, '(cm**2/sec) '/' athkdf = ',e12.6,' isopycnal thickness '
&,'diffusion (cm**2/sec) '/' ah = ',e12.6,' cm**2/sec for '
&,' am = ',e12.6,' cm**2/sec for horizontal mixing of momentum'/)
ft = c1/(4.0*ahisop*dtts)
delta_iso = dxt(1)*cst(jmt/2)*dzt(1)*ft
i_delta = 1
j_delta = 1
k_delta = 1
do jrow=2,jmt-1
do i=2,imt-1
do k=1,km
delta1 = dxt(i)*cst(jrow)*dzt(k)*ft
delta2 = dyt(jrow)*dzt(k)*ft
if (delta_iso .ge. delta1 .or. delta_iso .ge. delta2) then
i_delta = i
j_delta = jrow
k_delta = k
delta_iso = min(delta1,delta2)
endif
enddo
enddo
enddo
if (delta_iso .lt. p5) then
s_minus = (c1 - sqrt(c1 - 4.0*delta_iso**2))/(c2*delta_iso)
s_plus = c1/s_minus
else
s_minus = c0
s_plus = c0
endif
write(stdout,'(a)')
&'------------------------------------------------'
write(stdout,'(a,e14.7)')
&'The isopycnal diffusion grid factor delta_iso =',delta_iso
write(stdout,'(a)')
&'was determined at the grid point'
write(stdout,'(a,i4,a,e14.7)')
&'dxt(',i_delta,') = ',dxt(i_delta)
write(stdout,'(a,i4,a,e14.7)')
&'dyt(',j_delta,') = ',dyt(j_delta)
write(stdout,'(a,i4,a,e14.7)')
&'dzt(',k_delta,') = ',dzt(k_delta)
write(stdout,'(a,i4,a,e14.7)')
&'cst(',j_delta,') = ',cst(j_delta)
write(stdout,'(a,e14.7)')
&'dtts = ',dtts
write(stdout,'(a,e14.7)')
&'ahisop = ',ahisop
# if defined O_full_tensor
if (delta_iso .lt. p5) then
write(stdout,'(/a,a,/2(a,e14.7)/)')
& ' The full isopycnal tensor will be rescaled'
&,' when the slope is in the range:'
&,' s(-) = ',s_minus,', s_(+) = ',s_plus
else
write(stdout,'(/a,a/)') 'Isopycnal slopes in the full tensor'
&, 'will not be rescaled since delta_iso > 0.5'
endif
# else
write(stdout,'(/a,e14.7/)')
& 'Maximum allowable isopycnal slope is specified as slmx = ',slmx
# endif
write(stdout,'(a)')
&'------------------------------------------------'
!-----------------------------------------------------------------------
! store the square root of the tracer timestep acceleration values
! into variable "dtxsqr" for use in isopycnal mixing
!-----------------------------------------------------------------------
do k=1,km
dtxsqr(k) = sqrt(dtxcel(k))
enddo
write (stdout,'(a)') ' Acceleration with depth "dtxsqr(k)"='
write (stdout,'(5(1x,e12.6))') (dtxsqr(k),k=1,km)
!-----------------------------------------------------------------------
! read ocean diffusion factor
!-----------------------------------------------------------------------
fisop(:,:,:) = 1.
fname = new_file_name ("O_diffac.nc")
inquire (file=trim(fname), exist=exists)
if (exists) then
ib(:) = 1
ic(:) = 1
ic(1) = imtm2
ic(2) = jmtm2
ic(3) = km
call openfile (fname, iou)
exists = inqvardef('O_diffac', iou)
if (exists) then
write(*,*) "reading ocean diffusion factor from: ",fname
call getvara ('O_diffac', iou, imtm2*jmtm2*km, ib, ic, tmpijk
&, c1, c0)
fisop(2:imtm1,2:jmtm1,1:km) = tmpijk(1:imtm2,1:jmtm2,1:km)
endif
endif
!-----------------------------------------------------------------------
! initialize arrays
!-----------------------------------------------------------------------
do j=1,jmw
do k=1,km
do i=1,imt
alphai(i,k,j) = c0
betai(i,k,j) = c0
enddo
enddo
enddo
do j=jsmw,jemw
do n=1,2
do k=1,km
do i=1,imt
ddxt(i,k,j,n) = c0
enddo
enddo
enddo
enddo
do j=1,jemw
do n=1,2
do k=1,km
do i=1,imt
ddyt(i,k,j,n) = c0
enddo
enddo
enddo
enddo
do j=1,jmw
do n=1,2
do k=0,km
do i=1,imt
ddzt(i,k,j,n) = c0
enddo
enddo
enddo
enddo
do j=jsmw,jemw
do k=1,km
do i=1,imt
do kr=0,1
do ip=0,1
Ai_ez(i,k,j,ip,kr) = c0
enddo
enddo
do kr=0,1
do ip=0,1
Ai_bx(i,k,j,ip,kr) = c0
enddo
enddo
do kr=0,1
do jq=0,1
Ai_by(i,k,j,jq,kr) = c0
enddo
enddo
K33(i,k,j) = c0
K11(i,k,j) = c0
# if defined O_gent_mcwilliams
adv_vetiso(i,k,j) = c0
# endif
enddo
enddo
enddo
do j=1,jemw
do k=1,km
do i=1,imt
do kr=0,1
do jq=0,1
Ai_nz(i,k,j,jq,kr) = c0
enddo
enddo
K22(i,k,j) = c0
# if defined O_gent_mcwilliams
adv_vntiso(i,k,j) = c0
# endif
enddo
enddo
enddo
# if defined O_gent_mcwilliams
do j=jsmw,jemw
do k=0,km
do i=1,imt
adv_vbtiso(i,k,j) = c0
adv_fbiso(i,k,j) = c0
enddo
enddo
enddo
# endif
return
end
subroutine elements (joff, js, je, is, ie)
!=======================================================================
! Estimate alpha, beta, and normal gradients on faces of T cells
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, is, ie, n, kp1, jrow, joff
real dens, tq, sq, drodt, drods, tprime, sprime
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "grdvar.h"
include "levind.h"
include "mw.h"
include "state.h"
include "isopyc.h"
include "dens.h"
!-----------------------------------------------------------------------
! alpha and beta at centers of T cells
!-----------------------------------------------------------------------
do j=js,je
do k=1,km
do i=is,ie
tprime = t(i,k,j,1,taum1)-to(k)
sprime = t(i,k,j,2,taum1)-so(k)
alphai(i,k,j) = drodt (tprime, sprime, k)
betai(i,k,j) = drods (tprime, sprime, k)
enddo
enddo
call setbcx (alphai(1,1,j), imt, km)
call setbcx (betai(1,1,j), imt, km)
enddo
!-----------------------------------------------------------------------
! gradients at bottom face of T cells
!-----------------------------------------------------------------------
do j=js,je
do n=1,2
do k=1,km
kp1 = min(k+1,km)
do i=is,ie
ddzt(i,k,j,n) = tmask(i,kp1,j)*dzwr(k)*
& (t(i,k,j,n,taum1) - t(i,kp1,j,n,taum1))
enddo
enddo
do i=is,ie
ddzt(i,0,j,n) = c0
enddo
call setbcx (ddzt(1,0,j,n), imt, km+1)
enddo
enddo
!-----------------------------------------------------------------------
! gradients at eastern face of T cells
!-----------------------------------------------------------------------
# if defined O_full_tensor
do j=js,je
# else
do j=max(js-1,2), je-1
# endif
jrow = j + joff
do n=1,2
do k=1,km
do i=is,ie
ddxt(i,k,j,n) = tmask(i,k,j)*tmask(i+1,k,j)*cstr(jrow)*
& dxur(i)*(t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
enddo
enddo
call setbcx (ddxt(1,1,j,n), imt, km)
enddo
enddo
!-----------------------------------------------------------------------
! gradients at northern face of T cells
!-----------------------------------------------------------------------
do j=max(js-1,1), je-1
jrow = j + joff
do n=1,2
do k=1,km
kp1 = min(k+1,km)
do i=is,ie
ddyt(i,k,j,n) = tmask(i,k,j)*tmask(i,k,j+1)*dyur(jrow)*
& (t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1))
enddo
enddo
call setbcx (ddyt(1,1,j,n), imt, km)
enddo
enddo
return
end
subroutine isopyc (joff, js, je, is, ie)
!=======================================================================
! Compute isopycnal diffusion coefficients Ai_ez, Ai_nz, Ai_bx,
! Ai_by, and the Gent/McWilliams advection velocities.
! input:
! joff = offset relating row "j" in the MW to latitude "jrow"
! js = starting row within the MW for calculations
! je = ending row within the MW for calculations
! is = starting index longitude within the MW
! ie = ending index longitude within the MW
! output:
! Ai_ez = diffusion coefficient centered on east face of T cells
! Ai_nz = diffusion coefficient centered on north face of T cells
! Ai_bx = diffusion coefficient centered on bottom face of T cells
! Ai_by = diffusion coefficient centered on bottom face of T cells
# if defined O_gent_mcwilliams
! adv_vetiso = isopycnal advective vel on east face of "T" cell
! adv_vntiso = isopycnal advective vel on north face of "T" cell
! (Note: this includes the cosine factor as in "adv_vnt")
! adv_vbtiso = isopycnal advective vel on bottom face of "T" cell
# endif
!=======================================================================
implicit none
integer istrt, is, iend, ie, joff, js, je
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
!-----------------------------------------------------------------------
! set local constants
!-----------------------------------------------------------------------
istrt = max(2,is)
iend = min(imt-1,ie)
!-----------------------------------------------------------------------
! estimate alpha, beta, and gradients on sides of T cells
!-----------------------------------------------------------------------
call elements (joff, js, je, is, ie)
!-----------------------------------------------------------------------
! compute Ai_ez centered on eastern face of T cells
! 1 <= MW < last: calculate for rows jsmw..jemw
! last MW : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------
call ai_east (joff, max(js-1,2), je-1, is, ie)
!-----------------------------------------------------------------------
! compute Ai_nz centered on the northern face of T cells
! first MW : calculate for rows 1..jemw
! 1 < MW < last: calculate for rows jsmw..jemw
! last MW : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------
call ai_north (joff, max(js-1,1), je-1, is, ie)
!-----------------------------------------------------------------------
! evaluate Ai_bx & Ai_by centered on bottom face of T cells
! 1 <= MW < last: calculate for rows jsmw..jemw
! last MW : calculate for rows jsmw..<=jemw
!-----------------------------------------------------------------------
call ai_bottom (joff, max(js-1,2), je-1, is, ie)
# if defined O_gent_mcwilliams
!-----------------------------------------------------------------------
! compute isopycnal advective velocities for tracers
! first MW : calculate (viso) for rows 1..jemw
! (uiso,wiso) for rows 2..jemw
! 1 < MW < last: calculate (uiso,viso,wiso) for rows jsmw..jemw
! last MW : calculate (uiso,viso,wiso) for rows jsmw..<=jemw
!-----------------------------------------------------------------------
call isopyc_adv (joff, max(js-1,1), je-1, is, ie)
# endif
return
end
subroutine ai_east (joff, js, je, is, ie)
!=======================================================================
! compute "Ai_ez" & "K11" at the center of eastern face of T cells.
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is
real sc, dzt4r, ai0, sumz, sxe, sumy, facty
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "coord.h"
include "grdvar.h"
include "hmixc.h"
include "isopyc.h"
include "mw.h"
!-----------------------------------------------------------------------
! compute "Ai_ez" on east face of T cell. Note re-scaling factor
! which reduces mixing coefficient "Ai" where slope "sxe"
! exceeds the critical slope "sc" for the small slope approx. For
! the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,km
sc = c1/(slmxr*dtxsqr(k))
dzt4r = p5*dzt2r(k)
do i=2,imtm1
Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k))
# if defined O_full_tensor
sumz = c0
do kr=0,1
do ip=0,1
sxe = abs(drodxe(i,k,j,ip)/(drodze(i,k,j,ip,kr)+epsln))
if (delta_iso .lt. p5 .and.
& (s_minus .lt. sxe .and. sxe .lt. s_plus)) then
Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j)
& *delta_iso*(sxe + c1/sxe)
else
Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j)
endif
sumz = sumz + dzw(k-1+kr)*Ai_ez(i,k,j,ip,kr)
& *drodze(i,k,j,ip,kr)**2/(drodxe(i,k,j,ip)**2
& + p5*(drodye(i,k,j,ip,0)**2 + drodye(i,k,j,ip,1)**2)
& + drodze(i,k,j,ip,kr)**2 + epsln)
enddo
enddo
sumy = c0
do jq=0,1
facty = csu(jrow-1+jq)*dyu(jrow-1+jq)
do ip=0,1
sumy = sumy + facty*Ai0*tmask(i,k,j)*tmask(i+1,k,j)
& *drodye(i,k,j,ip,jq)**2/(drodxe(i,k,j,ip)**2
& + drodye(i,k,j,ip,jq)**2 + epsln
& + p5*(drodze(i,k,j,ip,0)**2 + drodze(i,k,j,ip,1)**2))
enddo
enddo
K11(i,k,j) = dzt4r*sumz + p25*cstdytr(jrow)*sumy
# else
sumz = c0
do kr=0,1
do ip=0,1
sxe = abs(drodxe(i,k,j,ip)/(drodze(i,k,j,ip,kr)+epsln))
# if defined O_dm_taper
Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j)
& *p5*(c1-tanh((sxe-del_dm)*s_dmr))
# else
if (sxe .gt. sc) then
Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j)
& *(sc/(sxe + epsln))**2
else
Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j)
endif
# endif
sumz = sumz + dzw(k-1+kr)*Ai_ez(i,k,j,ip,kr)
enddo
enddo
K11(i,k,j) = dzt4r*sumz
# endif
enddo
enddo
call setbcx (Ai_ez(1,1,j,0,0), imt, km)
call setbcx (Ai_ez(1,1,j,1,0), imt, km)
call setbcx (Ai_ez(1,1,j,0,1), imt, km)
call setbcx (Ai_ez(1,1,j,1,1), imt, km)
call setbcx (K11(1,1,j), imt, km)
enddo
return
end
subroutine ai_north (joff, js, je, is, ie)
!=======================================================================
! compute "Ai_nz" & "K22" at the center of northern face of T cells.
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is
real sc, dzt4r, ai0, sumz, syn, sumx
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "coord.h"
include "grdvar.h"
include "hmixc.h"
include "isopyc.h"
include "mw.h"
!-----------------------------------------------------------------------
! compute "Ai_nz" on north face of T cell. Note re-scaling factor
! which reduces mixing coefficient "Ai" where slope "syn"
! exceeds the critical slope "sc" for the small slope approx. For
! the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,km
sc = c1/(slmxr*dtxsqr(k))
dzt4r = p5*dzt2r(k)
do i=2,imtm1
Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k))
# if defined O_full_tensor
sumz = c0
do kr=0,1
do jq=0,1
syn=abs(drodyn(i,k,j,jq)/(drodzn(i,k,j,jq,kr)+epsln))
if (delta_iso .lt. p5 .and.
& (s_minus .lt. syn .and. syn .lt. s_plus)) then
Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
& *delta_iso*(syn + c1/syn)
else
Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
endif
sumz = sumz + dzw(k-1+kr)*Ai_nz(i,k,j,jq,kr)
& *drodzn(i,k,j,jq,kr)**2/(
& + p5*(drodxn(i,k,j,0,jq)**2 + drodxn(i,k,j,1,jq)**2)
& + drodyn(i,k,j,jq)**2 + drodzn(i,k,j,jq,kr)**2 + epsln)
enddo
enddo
sumx = c0
do ip=0,1
do jq=0,1
sumx =sumx+dxu(i-1+ip)*Ai0*tmask(i,k,j)*tmask(i,k,j+1)
& *drodxn(i,k,j,ip,jq)**2/(drodxn(i,k,j,ip,jq)**2
& + drodyn(i,k,j,jq)**2 + epsln
& + p5*(drodzn(i,k,j,jq,0)**2 + drodzn(i,k,j,jq,1)**2))
enddo
enddo
K22(i,k,j) = dzt4r*sumz + p25*dxtr(i)*sumx
# else
sumz = c0
do kr=0,1
do jq=0,1
syn =abs(drodyn(i,k,j,jq)/(drodzn(i,k,j,jq,kr)+epsln))
# if defined O_dm_taper
Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
& *p5*(c1-tanh((syn-del_dm)*s_dmr))
# else
if (syn .gt. sc) then
Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
& *(sc/(syn + epsln))**2
else
Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1)
endif
# endif
sumz = sumz + dzw(k-1+kr)*Ai_nz(i,k,j,jq,kr)
enddo
enddo
K22(i,k,j) = dzt4r*sumz
# endif
enddo
enddo
call setbcx (Ai_nz(1,1,j,0,0), imt, km)
call setbcx (Ai_nz(1,1,j,1,0), imt, km)
call setbcx (Ai_nz(1,1,j,0,1), imt, km)
call setbcx (Ai_nz(1,1,j,1,1), imt, km)
call setbcx (K22(1,1,j), imt, km)
enddo
return
end
subroutine ai_bottom (joff, js, je, is, ie)
!=======================================================================
! compute Ai_bx, Ai_by, and K33 at the center of the bottom face of
! T cells.
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is
real sc, kp1, ai0, sumx, sxb, sumy, facty, syb
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "grdvar.h"
include "hmixc.h"
include "isopyc.h"
include "mw.h"
!-----------------------------------------------------------------------
! compute "Ai_bx", "Ai_by", & K33 on bottom face of T cell. Note
! re-scaling factor to reduce mixing coefficient "Ai" where slopes
! "sxb" and "syb" exceeds the critical slope "sc" for the small slope approx. For
! the full tensor, it is re-scaled if outside the stable range.
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,kmm1
sc = c1/(slmxr*dtxsqr(k))
kp1 = min(k+1,km)
do i=2,imtm1
Ai0 = ahisop*p5*(fisop(i,jrow,k+1) + fisop(i,jrow,k))
# if defined O_full_tensor
! eastward slopes at the base of T cells
sumx = c0
do ip=0,1
do kr=0,1
sxb = abs(drodxb(i,k,j,ip,kr)/(drodzb(i,k,j,kr)+epsln))
if (delta_iso .lt. p5 .and.
& (s_minus .lt. sxb .and. sxb .lt. s_plus)) then
Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j)
& *delta_iso*(sxb + c1/sxb)
else
Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j)
endif
sumx = sumx + dxu(i-1+ip)*Ai_bx(i,k,j,ip,kr)
& *drodxb(i,k,j,ip,kr)**2/(drodxb(i,k,j,ip,kr)**2
& + p5*(drodyb(i,k,j,0,kr)**2 + drodyb(i,k,j,1,kr)**2)
& + drodzb(i,k,j,kr)**2 + epsln)
enddo
enddo
! northward slopes at the base of T cells
sumy = c0
jrow = j + joff
do jq=0,1
facty = csu(jrow-1+jq)*dyu(jrow-1+jq)
do kr=0,1
syb = abs(drodyb(i,k,j,jq,kr)/(drodzb(i,k,j,kr)+epsln))
if (delta_iso .lt. p5 .and.
& (s_minus .lt. syb .and. syb .lt. s_plus)) then
Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j)
& *delta_iso*(syb + c1/syb)
else
Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j)
endif
sumy = sumy + facty*Ai_by(i,k,j,jq,kr)
& *drodyb(i,k,j,jq,kr)**2/(
& p5*(drodxb(i,k,j,0,kr)**2 + drodxb(i,k,j,1,kr)**2)
& + drodyb(i,k,j,jq,kr)**2
& + drodzb(i,k,j,kr)**2 + epsln)
enddo
enddo
K33(i,k,j) = dxt4r(i)*sumx + dyt4r(jrow)*cstr(jrow)*sumy
# else
! eastward slopes at the base of T cells
sumx = c0
do ip=0,1
do kr=0,1
sxb = abs(drodxb(i,k,j,ip,kr)/(drodzb(i,k,j,kr)+epsln))
# if defined O_dm_taper
Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j)
& *p5*(c1-tanh((sxb-del_dm)*s_dmr))
# else
if (sxb .gt. sc) then
Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j)
& *(sc/(sxb + epsln))**2
else
Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j)
endif
# endif
sumx = sumx + dxu(i-1+ip)*Ai_bx(i,k,j,ip,kr)*sxb**2
enddo
enddo
! northward slopes at the base of T cells
sumy = c0
do jq=0,1
facty = csu(jrow-1+jq)*dyu(jrow-1+jq)
do kr=0,1
syb = abs(drodyb(i,k,j,jq,kr)/(drodzb(i,k,j,kr)+epsln))
# if defined O_dm_taper
Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j)
& *p5*(c1-tanh((syb-del_dm)*s_dmr))
# else
if (syb .gt. sc) then
Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j)
& *(sc/(syb + epsln))**2
else
Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j)
endif
# endif
sumy = sumy + facty*Ai_by(i,k,j,jq,kr)*syb**2
enddo
enddo
K33(i,k,j) = dxt4r(i)*sumx + dyt4r(jrow)*cstr(jrow)*sumy
# endif
enddo
enddo
call setbcx (Ai_bx(1,1,j,1,0), imt, km)
call setbcx (Ai_bx(1,1,j,0,0), imt, km)
call setbcx (Ai_bx(1,1,j,1,1), imt, km)
call setbcx (Ai_bx(1,1,j,0,1), imt, km)
call setbcx (Ai_by(1,1,j,1,0), imt, km)
call setbcx (Ai_by(1,1,j,0,0), imt, km)
call setbcx (Ai_by(1,1,j,1,1), imt, km)
call setbcx (Ai_by(1,1,j,0,1), imt, km)
call setbcx (K33(1,1,j), imt, km)
enddo
return
end
subroutine isoflux (joff, js, je, is, ie, n)
!=======================================================================
! isopycnal diffusive tracer fluxes are computed.
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, jrow, joff, km1kr, kpkr
integer n, ie, is
real dzt4r, sumz, flux_x, ai0, sumy, facty, csu_dzt4r, flux_y
real sumx
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "coord.h"
include "grdvar.h"
include "hmixc.h"
include "isopyc.h"
include "levind.h"
include "mw.h"
include "vmixc.h"
!-----------------------------------------------------------------------
! construct total isopycnal tracer flux at east face of "T" cells
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,km
dzt4r = p5*dzt2r(k)
do i=2,imtm1
sumz = c0
do kr=0,1
km1kr = max(k-1+kr,1)
kpkr = min(k+kr,km)
do ip=0,1
sumz = sumz - Ai_ez(i,k,j,ip,kr)
& *(t(i+ip,km1kr,j,n,taum1)-t(i+ip,kpkr,j,n,taum1))
& *drodxe(i,k,j,ip)
# if defined O_full_tensor
& *drodze(i,k,j,ip,kr)/(drodxe(i,k,j,ip)**2
& + p5*(drodye(i,k,j,ip,0)**2 + drodye(i,k,j,ip,1)**2)
& + drodze(i,k,j,ip,kr)**2 + epsln)
# else
& /(drodze(i,k,j,ip,kr) + epsln)
# endif
enddo
enddo
flux_x = dzt4r*sumz
# if defined O_full_tensor
Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k))
sumy = c0
do jq=0,1
facty = csu(jrow-1+jq)
do ip=0,1
sumy = sumy - facty*Ai0*tmask(i,k,j)*tmask(i+1,k,j)
& *(t(i+ip,k,j+jq,n,taum1)-t(i+ip,k,j-1+jq,n,taum1))
& *drodye(i,k,j,ip,jq)*drodxe(i,k,j,ip)/(drodxe(i,k,j,ip)**2
& + drodye(i,k,j,ip,jq)**2 + epsln
& + p5*(drodze(i,k,j,ip,0)**2 + drodze(i,k,j,ip,1)**2))
enddo
enddo
flux_x = flux_x + p25*cstdytr(jrow)*sumy
# endif
diff_fe(i,k,j) = diff_fe(i,k,j) + K11(i,k,j)*cstdxur(i,j)*
& (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
& + flux_x
enddo
enddo
call setbcx (diff_fe(1,1,j), imt, km)
enddo
!-----------------------------------------------------------------------
! construct total isopycnal tracer flux at north face of "T" cells
!-----------------------------------------------------------------------
do j=js-1,je
jrow = j + joff
do k=1,km
csu_dzt4r = csu(jrow)*p5*dzt2r(k)
do i=2,imtm1
sumz = c0
do kr=0,1
km1kr = max(k-1+kr,1)
kpkr = min(k+kr,km)
do jq=0,1
sumz = sumz - Ai_nz(i,k,j,jq,kr)
& *(t(i,km1kr,j+jq,n,taum1)-t(i,kpkr,j+jq,n,taum1))
& *drodyn(i,k,j,jq)
# if defined O_full_tensor
& *drodzn(i,k,j,jq,kr)/(
& + p5*(drodxn(i,k,j,0,jq)**2 + drodxn(i,k,j,1,jq)**2)
& + drodyn(i,k,j,jq)**2 + drodzn(i,k,j,jq,kr)**2 + epsln)
# else
& /(drodzn(i,k,j,jq,kr) + epsln)
# endif
enddo
enddo
flux_y = csu_dzt4r*sumz
# if defined O_full_tensor
Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k))
sumx = c0
do ip=0,1
do jq=0,1
sumx = sumx - Ai0*tmask(i,k,j)*tmask(i,k,j+1)
& *(t(i+ip,k,j+jq,n,taum1)-t(i-1+ip,k,j+jq,n,taum1))
& *cstr(jrow+jq)
& *drodxn(i,k,j,ip,jq)*drodyn(i,k,j,jq)/(
& drodxn(i,k,j,ip,jq)**2 + drodyn(i,k,j,jq)**2 + epsln
& + p5*(drodzn(i,k,j,jq,0)**2 + drodzn(i,k,j,jq,1)**2))
enddo
enddo
flux_y = flux_y + p25*csu(jrow)*dxtr(i)*sumx
# endif
diff_fn(i,k,j) = diff_fn(i,k,j)+K22(i,k,j)*csu_dyur(jrow)*
& (t(i,k,j+1,n,taum1)-t(i,k,j,n,taum1))
& + flux_y
enddo
enddo
call setbcx (diff_fn(1,1,j), imt, km)
enddo
!-----------------------------------------------------------------------
! compute the vertical tracer flux "diff_fbiso" containing the K31
! and K32 components which are to be solved explicitly. The K33
! component will be treated implicitly. Note that there are some
! cancellations of dxu(i-1+ip) and dyu(jrow-1+jq)
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,kmm1
do i=2,imtm1
sumx = c0
do ip=0,1
do kr=0,1
sumx = sumx
& - Ai_bx(i,k,j,ip,kr)*cstr(jrow)*
& (t(i+ip,k+kr,j,n,taum1) - t(i-1+ip,k+kr,j,n,taum1))
& *drodxb(i,k,j,ip,kr)
# if defined O_full_tensor
& *drodzb(i,k,j,kr)/(drodxb(i,k,j,ip,kr)**2
& + p5*(drodyb(i,k,j,0,kr)**2 + drodyb(i,k,j,1,kr)**2)
& + drodzb(i,k,j,kr)**2 + epsln)
# else
& /(drodzb(i,k,j,kr) + epsln)
# endif
enddo
enddo
sumy = c0
do jq=0,1
do kr=0,1
sumy = sumy
& - Ai_by(i,k,j,jq,kr)*csu(jrow-1+jq)*
& (t(i,k+kr,j+jq,n,taum1)-t(i,k+kr,j-1+jq,n,taum1))
& *drodyb(i,k,j,jq,kr)
# if defined O_full_tensor
& *drodzb(i,k,j,kr)/(
& p5*(drodxb(i,k,j,0,kr)**2 + drodxb(i,k,j,0,kr)**2)
& + drodyb(i,k,j,jq,kr)**2
& + drodzb(i,k,j,kr)**2 + epsln)
# else
& /(drodzb(i,k,j,kr) + epsln)
# endif
enddo
enddo
diff_fbiso(i,k,j) = dxt4r(i)*sumx
& + dyt4r(jrow)*cstr(jrow)*sumy
enddo
enddo
do i=2,imtm1
diff_fbiso(i,0,j) = c0
diff_fbiso(i,km,j) = c0
enddo
call setbcx (diff_fbiso(1,0,j), imt, km+1)
enddo
# if defined O_gent_mcwilliams
!-----------------------------------------------------------------------
! compute advective tracer flux at the center of the bottom face of
! the "T" cells
!-----------------------------------------------------------------------
do j=js,je
do k=1,kmm1
do i=2,imt-1
adv_fbiso(i,k,j) = adv_vbtiso(i,k,j)*
& (t(i,k,j,n,taum1) + t(i,k+1,j,n,taum1))
enddo
enddo
enddo
! now consider the top and bottom boundaries
do j=js,je
do i=2,imt-1
adv_fbiso(i,0,j) = c0
adv_fbiso(i,km,j) = c0
enddo
enddo
# endif
return
end
subroutine isopyc_adv (joff, js, je, is, ie)
!=======================================================================
! compute isopycnal transport velocities.
!=======================================================================
implicit none
integer i, k, j, ip, kr, jq, js, je, jrow, joff, km1, kp1
integer jstrt, ie, is
# if defined O_gent_mcwilliams
real sc, ath0, at, bt, stn, ab, bb, sbn, absstn, abssbn
real ath_t, ath_b, ste, sbe, absste, abssbe
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "accel.h"
include "coord.h"
include "grdvar.h"
include "isopyc.h"
include "levind.h"
include "mw.h"
include "switch.h"
# if defined O_time_averages
include "timeavgs.h"
# endif
real top_bc(km), bot_bc(km)
do k=1,km
top_bc(k) = c1
bot_bc(k) = c1
enddo
top_bc(1) = c0
bot_bc(km) = c0
!-----------------------------------------------------------------------
! compute the meridional component of the isopycnal mixing velocity
! at the center of the northern face of the "T" cells.
!-----------------------------------------------------------------------
do j=js,je
jrow = j + joff
do k=1,km
sc = c1/(slmxr*dtxsqr(k))
km1 = max(k-1,1)
kp1 = min(k+1,km)
do i=1,imt
Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k))
at = (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,km1,j)
& + alphai(i,km1,j+1))
bt = (betai(i,k,j) + betai(i,k,j+1) + betai(i,km1,j)
& + betai(i,km1,j+1))
stn = -(at*(ddyt(i,k,j,1) + ddyt(i,km1,j,1))
& + bt*(ddyt(i,k,j,2) + ddyt(i,km1,j,2))) /
& (at*(ddzt(i,km1,j,1) + ddzt(i,km1,j+1,1))
& + bt*(ddzt(i,km1,j,2) + ddzt(i,km1,j+1,2))+epsln)
ab = (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,kp1,j)
& + alphai(i,kp1,j+1))
bb = (betai(i,k,j) + betai(i,k,j+1) + betai(i,kp1,j)
& + betai(i,kp1,j+1))
sbn = -(ab*(ddyt(i,k,j,1) + ddyt(i,kp1,j,1))
& + bb*(ddyt(i,k,j,2) + ddyt(i,kp1,j,2))) /
& (ab*(ddzt(i,k,j,1) + ddzt(i,k,j+1,1))
& + bb*(ddzt(i,k,j,2) + ddzt(i,k,j+1,2))+epsln)
absstn = abs(stn)
abssbn = abs(sbn)
# if defined O_dm_taper
ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
& *p5*(c1-tanh((absstn-del_dm)*s_dmr))
ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
& *p5*(c1-tanh((abssbn-del_dm)*s_dmr))
# else
if (absstn .gt. sc) then
ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
& *(sc/(absstn + epsln))**2
else
ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1)
endif
if (abssbn .gt. sc) then
ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
& *(sc/(abssbn + epsln))**2
else
ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1)
endif
# endif
adv_vntiso(i,k,j) = -(ath_t*stn*top_bc(k) -
& ath_b*sbn*bot_bc(k))*dztr(k)*csu(jrow)
enddo
enddo
enddo
!-----------------------------------------------------------------------
! compute the zonal component of the isopycnal mixing velocity
! at the center of the eastern face of the "T" grid cell.
!-----------------------------------------------------------------------
jstrt = max(js,jsmw)
do j=jstrt,je
jrow = j + joff
do k=1,km
sc = c1/(slmxr*dtxsqr(k))
km1 = max(k-1,1)
kp1 = min(k+1,km)
do i=1,imtm1
Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k))
at = (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,km1,j)
& + alphai(i+1,km1,j))
bt = (betai(i,k,j) + betai(i+1,k,j) + betai(i,km1,j)
& + betai(i+1,km1,j))
ste =-(at*(ddxt(i,k,j,1) + ddxt(i,km1,j,1))
& + bt*(ddxt(i,k,j,2) + ddxt(i,km1,j,2)))
& / (at*(ddzt(i,km1,j,1) + ddzt(i+1,km1,j,1))
& + bt*(ddzt(i,km1,j,2) + ddzt(i+1,km1,j,2))+epsln)
ab = (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,kp1,j)
& + alphai(i+1,kp1,j))
bb = (betai(i,k,j) + betai(i+1,k,j) + betai(i,kp1,j)
& + betai(i+1,kp1,j))
sbe =-(ab*(ddxt(i,k,j,1) + ddxt(i,kp1,j,1))
& + bb*(ddxt(i,k,j,2) + ddxt(i,kp1,j,2))) /
& (ab*(ddzt(i,k,j,1) + ddzt(i+1,k,j,1))
& + bb*(ddzt(i,k,j,2) + ddzt(i+1,k,j,2))+epsln)
absste = abs(ste)
abssbe = abs(sbe)
# if defined O_dm_taper
ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
& *p5*(c1-tanh((absste-del_dm)*s_dmr))
ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
& *p5*(c1-tanh((abssbe-del_dm)*s_dmr))
# else
if (absste .gt. sc) then
ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
& *(sc/(absste + epsln))**2
else
ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j)
endif
if (abssbe .gt. sc) then
ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
& *(sc/(abssbe + epsln))**2
else
ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j)
endif
# endif
adv_vetiso(i,k,j) = -(ath_t*ste*top_bc(k) -
& ath_b*sbe*bot_bc(k))*dztr(k)
enddo
enddo
enddo
! set the boundary conditions
do j=jstrt,je
call setbcx (adv_vetiso(1,1,j), imt, km)
enddo
!-----------------------------------------------------------------------
! compute the vertical component of the isopycnal mixing velocity
! at the center of the bottom face of the "T" cells, using the
! continuity equation for the isopycnal mixing velocities
!-----------------------------------------------------------------------
do j=jstrt,je
do i=1,imt
adv_vbtiso(i,0,j) = c0
enddo
enddo
do j=jstrt,je
jrow = j + joff
do k=1,kmm1
do i=2,imt
adv_vbtiso(i,k,j) = dzt(k)*cstr(jrow)*(
& (adv_vetiso(i,k,j) - adv_vetiso(i-1,k,j))*dxtr(i) +
& (adv_vntiso(i,k,j) - adv_vntiso(i,k,j-1))*dytr(jrow))
enddo
enddo
enddo
do j=jstrt,je
do k=1,kmm1
do i=2,imt
adv_vbtiso(i,k,j) = adv_vbtiso(i,k,j) + adv_vbtiso(i,k-1,j)
enddo
enddo
enddo
do j=jstrt,je
jrow = j + joff
do i=2,imt
adv_vbtiso(i,kmt(i,jrow),j) = c0
enddo
enddo
! set the boundary conditions
do j=jstrt,je
call setbcx (adv_vbtiso(1,0,j), imt, km+1)
enddo
# if defined O_time_averages
!-----------------------------------------------------------------------
! accumulate time average gm velocities
!-----------------------------------------------------------------------
if (timavgperts .and. .not. euler2) then
do j=jstrt,je
jrow = j + joff
do k=1,kmm1
do i=2,imt
ta_vetiso(i,k,jrow) = ta_vetiso(i,k,jrow) +
& adv_vetiso(i,k,j)
ta_vntiso(i,k,jrow) = ta_vntiso(i,k,jrow) +
& adv_vntiso(i,k,j)
ta_vbtiso(i,k,jrow) = ta_vbtiso(i,k,jrow) +
& adv_vbtiso(i,k,j)
enddo
enddo
enddo
endif
# endif
# endif
# endif
#endif
return
end