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