subroutine ppmix (j)
c
c=======================================================================
c ===
c This routine computes the vertical mixing coefficients VDC and ===
c VVC based on the Pacanowski and Philander (JPO, v11, 11, 1981) ===
c Richardson number mixing scheme. ===
c ===
c ------ ===
c Input: ===
c ------ ===
c ===
c J Current slab number. (integer) ===
c ===
c Common Blocks: ===
c ===
c /CPPMIX/ ===
c ===
c FRICMX maximum visc/diff coeff. (real; cm^2/s) ===
c ===
c /CVMIX/ ===
c ===
c FKPH bkgd tracer vertical viscosity coeff. (real; cm^2/s) ===
c FKPM bkgd momentum vertical viscosity coeff. (real; cm^2/s) ===
c ===
c /VERTSLABS/ ===
c ===
c DZZQZ vertical grid spacing, tracer points. (real array) ===
c DZZVQZ vertical grid spacing, velocity points. (real array) ===
c TDEPTH depth at the center of the tracer box. (real array) ===
c ===
c /WORKSPA/ ===
c ===
c FM masking array, tracer points, row J. (real array) ===
c FMP masking array, tracer points, row J+1. (real array) ===
c GM masking array, velocity points, row J. (real array) ===
c TB tracer variables, row J, previous. (real array) ===
c TBP tracer variables, row J+1, previous. (real array) ===
c UB zonal velocity, row J, previous. (real array) ===
c UBM zonal velocity, row J-1, previous. (real array) ===
c VB meridional velocity, row J, previous. (real array) ===
c VBM meridional velocity, row J-1, previous. (real array) ===
c ===
c ------- ===
c Output: ===
c ------- ===
c ===
c Common Blocks: ===
c ===
c /CVMIX/ ===
c ===
c DRHO d(density)/dz, row J, tracer box bottom. (real array)===
c DRHOP d(density)/dz, row J+1, tracer box bottom. (real array)===
c VDC vertical eddy diffusion coeff. (real array)===
c VVC vertical eddy viscosity coeff. (real array)===
c ===
c /CPPMIX/ ===
c ===
c RIU Richardson #, row J , U,V box bottoms. (real array)===
c RIUM Richardson #, row J-1, U,V box bottoms. (real array)===
c WRIUM mask for U,V grid, row J-1. (real array)===
c ===
c Calls: CONVECT, MIXLAYER, PDENVGRAD, RICH_NO ===
c ===
c=======================================================================
c
c-----------------------------------------------------------------------
c Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <cppmix.h>
#include <cvmix.h>
#include <workspa.h>
#include <vertslabs.h>
c
c-----------------------------------------------------------------------
c Define local data.
c-----------------------------------------------------------------------
c
integer j
#ifndef barotropic
integer i,k,ip1,im1,kp1
FLOAT
& fac
FLOAT
& drhou(imt,km),dudz2(imt,km),dudz2m(imt,km),dudz2t(imt,km),
& rit(imt,km)
c
save dudz2,dudz2m
c
equivalence (riu, drhou)
equivalence (rit, dudz2t)
c
c=======================================================================
c Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c Compute potential density gradient.
c-----------------------------------------------------------------------
c
c Across bottom of tracer boxes.
c
if(j.eq.1) then
call pdenvgrad (tb(1,1,1),tb(1,1,2),tdepth(1,1,jrs),
& dzzqz(1,1,0),drho)
else
do 10 i=1,imt
do 10 k=1,km
drho(i,k)=drhop(i,k)
10 continue
endif
c
call pdenvgrad (tbp(1,1,1),tbp(1,1,2),tdepth(1,1,jrn),
& dzzqz(1,1,1),drhop)
c
c Across bottom of velocity boxes.
c
do 20 k = 1, kmm1
kp1=k+1
do 20 i = 1, imtm1
ip1 = i + 1
drhou(i,k) = p25*(fm(i,kp1)*drho(i,k)+fm(ip1,kp1)*drho(ip1,k)
& +fmp(i,kp1)*drhop(i,k)+fmp(ip1,kp1)*drhop(ip1,k))
20 continue
c
c-----------------------------------------------------------------------
c Compute square of the velocity shear.
c-----------------------------------------------------------------------
c
c Across bottom of velocity boxes.
c
if (j.eq.1) then
do 30 k = 1, kmm1
kp1=k+1
do 30 i = 1, imtm1
dudz2m(i,k) = ( (ubm(i,k)-ubm(i,kp1))**2
& + (vbm(i,k)-vbm(i,kp1))**2 )
& /( dzzvqz(i,kp1,0)*dzzvqz(i,kp1,0) )
30 continue
else
do 40 k = 1, kmm1
do 40 i = 1, imtm1
dudz2m(i,k) = dudz2(i,k)
40 continue
end if
c
do 50 k = 1, kmm1
kp1=k+1
do 50 i = 1, imtm1
dudz2(i,k) = ((ub(i,k)-ub(i,kp1))**2+(vb(i,k)-vb(i,kp1))**2)
& /( dzzvqz(i,kp1,0)*dzzvqz(i,kp1,0) )
50 continue
c
c Across bottom of tracer boxes.
c
do 60 k = 1, kmm1
kp1=k+1
do 60 i = 2, imtm1
im1 = i - 1
dudz2t(i,k) = p25*( gm(i,kp1)*dudz2(i,k)
& + gm(im1,kp1)*dudz2(im1,k)
& + wrium(i,k)*dudz2m(i,k)
& + wrium(im1,k)*dudz2m(im1,k) )
60 continue
c
c-----------------------------------------------------------------------
c Compute Richardson number.
c-----------------------------------------------------------------------
c
c On bottom of velocity boxes.
c
call rich_no (drhou,dudz2,gm,riu)
c
c On bottom of tracer boxes.
c
call rich_no (drho,dudz2t,fm,rit)
c
c-----------------------------------------------------------------------
c Compute vertical viscosity coefficient on U,V and T grid box bottoms
c using Pacanowski and Philander formula.
c-----------------------------------------------------------------------
c
c On bottom of velocity boxes.
c
if(j.gt.1) then
do 70 k=1,kmm1
do 70 i=2,imtm2
if(riu(i,k).le.c0) then
vvc(i,k)=fkpm+fricmx
else
fac=c1/(c1+c5*riu(i,k))
vvc(i,k)=fkpm+fricmx*fac*fac
endif
70 continue
c
c On bottom of tracer boxes.
c
do 80 k=1,kmm1
do 80 i=2,imtm1
if(rit(i,k).le.c0) then
vdc(i,k)=fkph+fricmx
else
fac=c1/(c1+c5*rit(i,k))
vdc(i,k)=fkph+fricmx*fac*fac*fac
endif
80 continue
endif
c
c Transfer Richardson numbers for use with next row.
c
do 90 k=1,kmm1
do 90 i=1,imt
rium(i,k)=riu(i,k)
wrium(i,k)=gm(i,k+1)
90 continue
c
c-----------------------------------------------------------------------
c Determine mixed-layer depth, and set vertical mixing coefficients to
c WVMIX and WDMIX within the mixed layer.
c-----------------------------------------------------------------------
c
call mixlayer (j)
c
c-----------------------------------------------------------------------
c Determine vertical gravitational stability, and if unstable set the
c vertical mixing coefficients to VVCLIM and VDCLIM.
c-----------------------------------------------------------------------
c
call convect (j)
c
c-----------------------------------------------------------------------
c Set viscosity and diffusivity to zero in land.
c-----------------------------------------------------------------------
c
do 100 k=1,kmm1
do 100 i=2,imtm1
vvc(i,k) = gm(i,k+1)*vvc(i,k)
vdc(i,k) = fm(i,k+1)*vdc(i,k)
100 continue
c
c-----------------------------------------------------------------------
c Set viscosity and diffusivity at level KM.
c-----------------------------------------------------------------------
c
do 110 i=2,imtm1
vvc(i,km) = c0
vdc(i,km) = c0
110 continue
#endif
c
return
end