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 #include #include #include #include #include #include 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