subroutine vmix_aux
c
c=======================================================================
c ===
c This routine handles the auxilary code needed for processing ===
c vertical mixing coefficients. ===
c ===
c Entries: ===
c ===
c MIXLAYER calculates mixed-layer depth, and sets the vertical ===
c mixing coefficients within it to WVMIX and WDMIX. ===
c PDENVGRAD determines the vertical potential density gradient. ===
c CONVECT determines gravitational stability, and if unstable ===
c sets the vertical mixing coefficients to VVCLIM and ===
c VDCLIM. ===
c ===
c=======================================================================
c
c-----------------------------------------------------------------------
c Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
#include <fullwd.h>
#include <workspa.h>
#include <onedim.h>
#include <cvbc.h>
#include <rhomean.h>
#include <vertslabs.h>
#include <cvmix.h>
c
c-----------------------------------------------------------------------
c Define local data.
c-----------------------------------------------------------------------
c
logical first
integer i,j,k
#ifndef barotropic
integer ip1,kp1
FLOAT
* atg,del_P,del_th,drhov,P1,PR,q,sum,th,z1
FLOAT
* rhoa(imt,km),tdown(imt,km),zdown(imt,km)
#endif
#if (defined ext_tide & defined mixtide) | defined posmxtid
FLOAT
* tide_depth,tide_val
#endif
FLOAT
* crho,denom,dept,depv,expans,expant,fac,f0,f0min,mldt,numer,
* s1,sbfu,srfu,stress,t1,ustar
FLOAT
* depth(imt,km),drhox(imt,km),dzz(imt,kmp1),sx(imt,km),
* tx(imt,km)
c
#ifndef barotropic
parameter (PR=c0)
#endif
save crho,f0min,first
data first /.true./
c
c=======================================================================
c Begin executable code. ==============================================
c=======================================================================
c
entry mixlayer(j)
c
c=======================================================================
c ===
c This entry determines the mixed-layer depth, depending on the ===
c scheme option MLDOPT, and then sets the vertical mixing ===
c coefficients within the mixed layer to WVMIX and WDMIX. ===
c ===
c On Input: ===
c ===
c J current slab number (integer). ===
c ===
c On Output: ===
c ===
c VDC vertical diffusion coefficient (real array). ===
c VVC vertical viscosity coefficient (real array). ===
c ===
c=======================================================================
c
if(first) then
f0min=c2*omega*sin(latmf0/radian)
crho=rho0/c1000
first=.false.
endif
c
c-----------------------------------------------------------------------
c Transfer mixed-layer depth computed previously to the south of the
c present row.
c-----------------------------------------------------------------------
c
if((j.gt.1).and.(mldopt.gt.0)) then
do 10 i=1,imt
mldum(i)=mldu(i)
10 continue
endif
c
c-----------------------------------------------------------------------
c Set the wind mixing coefficients to the prescribed depth.
c-----------------------------------------------------------------------
c
if(mldopt.eq.0) then
do 20 i=2,imtm1
vvc(i,1) = max( wvmix, vvc(i,1) )
vdc(i,1) = max( wdmix, vdc(i,1) )
#ifndef barotropic
do 20 k=2,km
depv=vdepth(i,k,jrn)+p5*dzvqz(i,k,jrn)
dept=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs)
if(depv.lt.mldmax) vvc(i,k)=max(wvmix,vvc(i,k))
if(dept.lt.mldmin) vdc(i,k)=max(wdmix,vdc(i,k))
#endif
20 continue
c
c-----------------------------------------------------------------------
c Compute Ekman depth from wind stress forcing.
c-----------------------------------------------------------------------
c
elseif(mldopt.eq.1) then
do 30 i=1,imt
stress=max(sqrt(smf(i,1)*smf(i,1)+smf(i,2)*smf(i,2)),taumin)
f0=max(abs(c2*omega*sine(i,j)),f0min)
ustar=sqrt(stress/crho)
mldu(i)=ekfac*ustar/f0
mldu(i)=max(mldu(i),mldmin)
mldu(i)=min(mldu(i),mldmax)
30 continue
c
c Mixed layer depth is defined at U,V points; interpolate to T points.
c
if(j.gt.1) then
do 40 i=1,imt
vvc(i,1) = max( wvmix, vvc(i,1) )
#ifndef barotropic
do 40 k=2,km
depv=vdepth(i,k,jrn)+p5*dzvqz(i,k,jrn)
if(depv.lt.mldu(i)) vvc(i,k)=max(wvmix,vvc(i,k))
#endif
40 continue
do 60 i=2,imt
vdc(i,1) = max( wdmix, vdc(i,1) )
mldt=p25*(mldum(i-1)+mldum(i)+mldu(i-1)+mldu(i))
#ifndef barotropic
do 50 k=2,km
dept=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs)
if(dept.lt.mldt) vdc(i,k)=max(wdmix,vdc(i,k))
50 continue
#endif
60 continue
endif
c
c-----------------------------------------------------------------------
c Use Niiler and Kraus (1977) formula to determine mixed-layer depth
c from surface momentum, net heat, shortwave radiation, and water
c fluxes. Keep Niiler and Kraus convention of downward flux to be
c negative.
c
c SBF(IMT,1) rate at which buoyancy is removed from water column (B0).
c SBF(IMT,2) buoyancy flux due to penetrating solar radiation (J0).
c-----------------------------------------------------------------------
c
elseif(mldopt.eq.2) then
c
c Compute buoyancy fluxes SBF at row J for T points.
c
if(j.eq.1) then
do 70 i=1,imt
if(fm(i,1).eq.c0) then
expant=tcoef
expans=scoef
else
t1=tb(i,1,1)
s1=tb(i,1,2)+smean
expant=-((((3.268166e-8*t1-4.480332e-6)*t1+3.005055e-4)*t1
* -1.819058e-2)*t1+6.793952e-2
* +s1*(((2.155e-8*t1-2.47401e-6)*t1+1.52876e-4)*t1
* -4.0899e-3)
* +(s1**1.5)*(-3.3092e-6*t1+1.0227e-4))/rho0
expans=((((5.3875e-9*t1-8.2467e-7)*t1+7.6438e-5)*t1
* -4.0899e-3)*t1+0.824493
* +(s1**0.5)*((-2.4819e-6*t1+1.53405e-4)*t1
* -8.58699e-3)+9.6628e-4*s1)/rho0
endif
fac=grav*expant
sbf(i,2)=-srf(i)*fac
sbf(i,1)=-stf(i,1)*fac-sbf(i,2)+stf(i,2)*grav*expans
70 continue
endif
c
c Compute buoyancy fluxes SBF at row J+1 for T points.
c
do 80 i=1,imt
if(fmp(i,1).eq.c0) then
expant=tcoef
expans=scoef
else
t1=tbp(i,1,1)
s1=tbp(i,1,2)+smean
expant=-((((3.268166e-8*t1-4.480332e-6)*t1+3.005055e-4)*t1
* -1.819058e-2)*t1+6.793952e-2
* +s1*(((2.155e-8*t1-2.47401e-6)*t1+1.52876e-4)*t1
* -4.0899e-3)
* +(s1**1.5)*(-3.3092e-6*t1+1.0227e-4))/rho0
expans=((((5.3875e-9*t1-8.2467e-7)*t1+7.6438e-5)*t1
* -4.0899e-3)*t1+0.824493
* +(s1**0.5)*((-2.4819e-6*t1+1.53405e-4)*t1
* -8.58699e-3)+9.6628e-4*s1)/rho0
endif
fac=grav*expant
sbfp(i,2)=-srfp(i)*fac
sbfp(i,1)=-stfp(i,1)*fac-sbfp(i,2)+stfp(i,2)*grav*expans
80 continue
c
c Interpolate SBF to U,V points. Calculate mixed-layer depth
#ifndef nkfix
c at U,V points; the denominator given by Niiler and Kraus is modified
c to avoid dividing by zero.
# else
c at U,V points.
#endif
c
do 90 i=1,imtm1
sbfu=p25*(sbf(i,1)+sbf(i+1,1)+sbfp(i,1)+sbfp(i+1,1))
srfu=p25*(sbf(i,2)+sbf(i+1,2)+sbfp(i,2)+sbfp(i+1,2))
stress=max(sqrt(smf(i,1)*smf(i,1)+smf(i,2)*smf(i,2)),taumin)
ustar=sqrt(stress/crho)
#ifndef nkfix
numer=c2*(mcoef*(ustar**3)-srfu/atth2o)
denom=-srfu-p5*((c1+ncoef)*sbfu-(c1-ncoef)*abs(sbfu))
denom=c2*(mcoef*wsdfac*ustar*ustar+denom)
mldu(i)=numer/denom
mldu(i)=max(mldu(i),mldmin)
mldu(i)=min(mldu(i),mldmax)
# else
numer=mcoef*(ustar**3)-srfu/atth2o
denom=-srfu-p5*((c1+ncoef)*sbfu-(c1-ncoef)*abs(sbfu))
c
if ((numer.ne.c0).or.(denom.ne.c0)) then
c
c If fraction is in range, use Niiler-Kraus directly.
c
if ( ((mldmin*abs(denom)).le.abs(numer)) .and.
& (abs(numer).le.(mldmax*abs(denom))) ) then
mldu(i)=numer/denom
c
c Else if denominator is nonzero, can use correct limit.
c
elseif ((mldmin*abs(denom)).lt.(mldmax*abs(denom))) then
if (abs(numer).lt.(mldmax*abs(denom))) then
mldu(i)=mldmin
else
mldu(i)=mldmax
end if
c
c When all else fails, retreat to default value.
c
else
mldu(i)=mldval
end if
else
mldu(i)=mldval
end if
#endif
90 continue
#ifndef cyclic
mldu(imt)=mldu(imtm1)
# else
mldu(imt)=mldu(2)
mldu(1 )=mldu(imtm1)
#endif
c
c Interpolate mixed-layer depth to T points. Then, if depth is less
c than mixed-layer depth, set diffusivity coefficients to large values
c on U,V and tracer grid box bottoms.
c
if(j.gt.1) then
do 110 i=2,imtm1
vvc(i,1) = max ( wvmix, vvc(i,1) )
vdc(i,1) = max ( wdmix, vdc(i,1) )
mldt=p25*(mldum(i-1)+mldum(i)+mldu(i-1)+mldu(i))
#ifndef barotropic
do 100 k=2,kmm1
depv=vdepth(i,k,jrn)+p5*dzvqz(i,k,jrn)
dept=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs)
if(depv.le.mldu(i)) vvc(i,k)=max(wvmix,vvc(i,k))
if(dept.le.mldt) vdc(i,k)=max(wdmix,vdc(i,k))
100 continue
#endif
110 continue
endif
c
c Transfer quantities computed to the north of the present row to
c be defined to the south in the computation of the next row.
c
do 120 i=1,imt
sbf(i,1)=sbfp(i,1)
sbf(i,2)=sbfp(i,2)
120 continue
endif
#if (defined ext_tide & defined mixtide) | defined posmxtid
c
c Set enhanced mixing due to tides.
c
if(j.gt.1) then
do 140 i=1,imt
dept=tide_depth(vdepth(i,km,jrn)+p5*dzvqz(i,km,jrn),i,j)
do 130 k=1,km
depv=vdepth(i,k,jrn)+p5*dzvqz(i,k,jrn)
vvc(i,k)=max(tide_val(vvclim,dept,depv,i,k), vvc(i,k))
130 continue
140 continue
do 160 i=2,imt
depv=tide_depth(tdepth(i,km,jrs)+p5*dzqz(i,km,jrs),i,j)
do 150 k=1,km
dept=tdepth(i,k,jrs)+p5*dzqz(i,k,jrs)
vdc(i,k)=max(tide_val(vdclim,depv,dept,i,k), vdc(i,k))
150 continue
160 continue
endif
#endif
return
c
entry pdenvgrad(tx,sx,depth,dzz,drhox)
c
c=======================================================================
c ===
c This entry calculates the vertical gradient in potential density, ===
c for use in checking graviational stability, and calculating the ===
c Richardson Number. ===
c ===
c On Input: ===
c ===
c TX temperature (real array; degrees Celsius). ===
c SX salinity anomaly from SMEAN (real array; PSU). ===
c DEPTH depths (real array; cm). ===
c DZZ vertical spacing between tracer levels. (real array) ===
c ===
c On Output: ===
c ===
c DRHOX vertical gradient in potential density (real array; ===
c nondimensional). ===
c ===
c=======================================================================
c
#ifndef barotropic
do 200 i=1,imt
do 200 k=1,km
t1=tx(i,k)
s1=sx(i,k)+smean-35.0
z1=depth(i,k)*cm2m
c theta1
del_P = PR - z1
th = t1
P1 = z1
atg=(((-2.1687E-16*th+1.8676E-14)*th-4.6206E-13)*P1
* +((2.7759E-12*th-1.1351E-10)*s1+((-5.4481E-14*th
* +8.733E-12)*th-6.7795E-10)*th+1.8741E-8))*P1
* +(-4.2393E-8*th+1.8932E-6)*s1
* +((6.6228E-10*th-6.836E-8)*th+8.5258E-6)*th+3.5803E-5
del_th = del_P*atg
th = t1 + 0.5*del_th
q = del_th
c theta2
P1 = z1 + 0.5*del_P
atg=(((-2.1687E-16*th+1.8676E-14)*th-4.6206E-13)*P1
* +((2.7759E-12*th-1.1351E-10)*s1+((-5.4481E-14*th
* +8.733E-12)*th-6.7795E-10)*th+1.8741E-8))*P1
* +(-4.2393E-8*th+1.8932E-6)*s1
* +((6.6228E-10*th-6.836E-8)*th+8.5258E-6)*th+3.5803E-5
del_th = del_P*atg
th = th + (1 - 1/sqrt(2.0))*(del_th - q)
q = (2-sqrt(2.0))*del_th + (-2+3/sqrt(2.0))*q
c theta3
atg=(((-2.1687E-16*th+1.8676E-14)*th-4.6206E-13)*P1
* +((2.7759E-12*th-1.1351E-10)*s1+((-5.4481E-14*th
* +8.733E-12)*th-6.7795E-10)*th+1.8741E-8))*P1
* +(-4.2393E-8*th+1.8932E-6)*s1
* +((6.6228E-10*th-6.836E-8)*th+8.5258E-6)*th+3.5803E-5
del_th = del_P*atg
th = th + (1 + 1/sqrt(2.0))*(del_th - q)
q = (2 + sqrt(2.0))*del_th + (-2-3/sqrt(2.0))*q
c theta4
P1 = z1+del_P
atg=(((-2.1687E-16*th+1.8676E-14)*th-4.6206E-13)*P1
* +((2.7759E-12*th-1.1351E-10)*s1+((-5.4481E-14*th
* +8.733E-12)*th-6.7795E-10)*th+1.8741E-8))*P1
* +(-4.2393E-8*th+1.8932E-6)*s1
* +((6.6228E-10*th-6.836E-8)*th+8.5258E-6)*th+3.5803E-5
del_th = del_P*atg
tdown(i,k)= th + (del_th - 2.0*q)/6.0
zdown(i,k)= PR
200 continue
c
c Compute potential density.
c
call state(tdown,sx,zdown,rhoa)
c
c Compute vertical gradient in potential density.
c
do 230 k=1,kmm1
do 230 i=1,imt
drhox(i,k)=(rhoa(i,k)-rhoa(i,k+1))/dzz(i,k+1)
230 continue
#endif
return
c
entry convect(j)
c
c=======================================================================
c ===
c This entry determines vertical gravitational stability, and if ===
c unstable sets the vertical mixing coefficients to VDCLIM and ===
c VVCLIM. ===
c ===
c On Input: ===
c ===
c J current slab number (integer). ===
c DRHO vertical gradient in potential density (real array). ===
c ===
c On Output: ===
c ===
c VDC vertical diffusion coefficient (real array). ===
c VVC vertical viscosity coefficient (real array). ===
c ===
c=======================================================================
c
#ifndef barotropic
if(j.gt.1) then
c
c Check diffusion coefficients.
c
do 300 i=1,imt
do 300 k=1,kmm1
if(drho(i,k).ge.c0) vdc(i,k)=max(vdclim,vdc(i,k))
300 continue
c
c Check viscosity coefficients after interpolating potential density
c gradient to U,V points.
c
do 310 i=1,imtm1
do 310 k=1,kmm1
kp1=k+1
ip1=i+1
sum=fm(i,kp1)+fm(ip1,kp1)+fmp(i,kp1)+fmp(ip1,kp1)
if (sum.gt.c0) then
drhov=(fm(i,kp1)*drho(i,k)+fm(ip1,kp1)*drho(ip1,k)+
* fmp(i,kp1)*drhop(i,k)+fmp(ip1,kp1)*drhop(ip1,k))/sum
if(drhov.ge.c0) vvc(i,k)=max(vvclim,vvc(i,k))
endif
310 continue
endif
#endif
return
end