C <%Z% %Y% %M% %I% %G%> subroutine VDIFFUSE(idim,km,im,nx, ndet,extmode, & H, U, V, T, S, X , det, B, & Tdeep, Sdeep, Xdeep, detdeep, & ieast,iwest,inorth, isouth, & Area_h, Bu, Kappa, Nu, Ri, Ks, dtime) implicit none c==> Arguments integer idim integer km integer im integer nx,ndet logical extmode real H( idim, km) real U( idim, km) real V( idim, km) real T( idim, km) real S( idim, km) real X( idim, km, 2, nx ) real det( idim, km, 2, ndet ) real B( idim, km) real Tdeep real Sdeep real Xdeep( nx ) real detdeep( ndet ) integer ieast( idim ) integer iwest( idim ) integer inorth( idim ) integer isouth( idim ) real Area_h( idim ) integer Bu( idim ) real hmax1 real dtime c==> Output arrays for diagnostics real Kappa( idim, km) real Nu( idim, km ) real Ri( idim, km ) real Ks( idim, km ) c==> Local Parameters real MINRI ! Richardson number at convection real KAPPA_CONVECT ! Kappa when convecting real NU_CONVECT ! Viscosity when convecting real NU_ML ! Nu at base of mixed layer (if isolated ) real NU_BOTT ! Nu at bottom of ocean (if extmode ) integer NITS ! Number of iterations parameter ( MINRI = 1.E-03 ) real MINN2 parameter ( MINN2 = 1.E-10 ) parameter ( KAPPA_CONVECT = 1.E+03 ) parameter ( NU_CONVECT = 1.E-02 ) parameter ( NU_ML = 1.E-06 ) parameter ( NU_BOTT = 1.E-04 ) ! parameter ( NITS = 4 ) parameter ( NITS = 1 ) c==> Local Variables real E1(idim, km) real A1(idim, km) real N2(idim, km) real Hu(idim, km) logical mix(idim, km) real Area_u_r( idim ) real Es(idim, km) real As(idim, km) real Alphs(idim) real Alpha(idim) real Wk1(idim) real Wk2(idim) ! Insitu conversion real Z(idim,km) real Tbottom(idim) ! Abyssal in situ temperature real Haa(idim) real Riu integer num_new integer nremain real Haas integer i,k,m,ntimes,k2 real dt2df,fact real dh real hmn real hmx real func1 real nu_bottom c==> Functions real FKAPPA real FKAPPA_SALT real FNU real INSITU real POTEMP real BYNCY CFPP$ EXPAND ( FKAPPA ) R CFPP$ EXPAND ( FNU ) R CFPP$ EXPAND ( INSITU ) R CFPP$ EXPAND ( POTEMP ) R CFPP$ EXPAND ( BYNCY ) R c--------------------------------------------------------------------- if ( extmode ) then k2 = km-1 do i=1,idim Kappa(i,km) = 0.0 Ks(i,km) = 0.0 Nu(i,km) = NU_BOTT Ri(i,km) = 0.0 mix(i,km) = .FALSE. enddo else k2 = km endif do k=1,k2 do i=1,idim Ri(i,k) = 0.0 N2(i,k) = 1.0 Nu(i,k) = 0.0 Kappa(i,k) = 0.0 Ks(i,k) = 0.0 mix(i,k)= .FALSE. enddo enddo do i=1,idim Wk1(i) = 0. enddo do i=1,im Wk1(i) = Area_h(i) + Area_h(ieast(i)) enddo do i=1,im Area_u_r(i) = Bu(i) / ( Wk1(i) + Wk1(inorth(i)) + 1.E-36 ) enddo dt2df = 2.*dtime/NITS ! compute Z as the pressure of the center of each layer (not depth) do i=1,im Z(i,1) = 0.5 * H(i,1) enddo do k=2,km do i=1,im Z(i,k) = Z(i,k-1) + 0.5*( H(i,k) + H(i,k-1) ) enddo enddo do i=1,im Tbottom(i)=INSITU(Tdeep,Sdeep,Z(i,km)+H(i,km)*.5) enddo nremain=NITS DO WHILE ( nremain .GT. 0 ) num_new = 0 do k=1,k2 call VHRICH (idim, km, im, idim-im, k, Ri(1,k), N2(1,k), & H,U,V,B, & iwest,isouth ) ! If we mixed on the last iteration, keep the diffusivity ! high, regardless of the present Richardson number. ! Also, update the array which flags high mixing (convection) do i=1,im ri(i,k)= amax1(Ri(i,k), MINRI) N2(i,k)= amax1(N2(i,k), MINN2) Kappa(i,k) = FKAPPA(Ri(i,k),N2(i,k)) Ks(i,k) = Kappa(i,k) Nu(i,k) = FNU(Ri(i,k),N2(i,k)) enddo ! i enddo ! k ! If we have not uncovered any new convection, just go ! ahead and do all the remaining iterations with the ! present values of the diffusivities. This will save ! time. To further save time, we could gather all the ! points for which some_new applies, and re-do them, ! but that will await a later release if ( num_new .EQ. 0 ) then dt2df=dt2df*nremain nremain=0 endif do k=1,km-1 do i=1,im A1(i,k) = Kappa(i,k)*dt2df/(H(i,k) + H(i,k+1)) As(i,k) = Ks(i,k) *dt2df/(H(i,k) + H(i,k+1)) enddo enddo k = km if ( extmode ) then do i=1,im A1(i,km) = 0. As(i,km) = 0. enddo else do i=1,im A1(i,km) = Kappa(i,km)*dt2df/H(i,km) As(i,km) = Ks(i,km) *dt2df/H(i,km) enddo endif ! extmode ! Convert T to Insitu temp do k=1,km do i=1,im T(i,k) = INSITU( T(i,k) , S(i,k) , Z(i,k) ) enddo enddo k=1 do i=1,im Alpha(i) = H(i,k)*A1(i,k) / (H(i,k)+A1(i,k)) Alphs(i) = H(i,k)*As(i,k) / (H(i,k)+As(i,k)) E1(i,k) = A1(i,k) / (H(i,k)+A1(i,k)) Es(i,k) = As(i,k) / (H(i,k)+As(i,k)) T(i,k) = T(i,k)* H(i,k) / (H(i,k)+A1(i,k)) S(i,k) = S(i,k)* H(i,k) / (H(i,k)+As(i,k)) enddo do m=1,nx do i=1,im X(i,k,1,m) = X(i,k,1,m)*H(i,k)/(H(i,k)+A1(i,k)) enddo enddo do m=1,ndet do i=1,im det(i,k,1,m) = det(i,k,1,m)*H(i,k)/(H(i,k)+A1(i,k)) enddo enddo do k=2,km do i=1,im Haa(i) = 1./( H(i,k) + A1(i,k) + Alpha(i) ) Haas = 1./( H(i,k) + As(i,k) + Alphs(i) ) Alpha(i) = A1(i,k)*(H(i,k)+Alpha(i))*Haa(i) Alphs(i) = As(i,k)*(H(i,k)+Alphs(i))*Haas E1(i,k) = A1(i,k)*Haa(i) Es(i,k) = As(i,k)*Haas T(i,k) = (T(i,k)*H(i,k) + A1(i,k-1)*T(i,k-1))*Haa(i) S(i,k) = (S(i,k)*H(i,k) + As(i,k-1)*S(i,k-1))*Haas enddo do m=1,nx do i=1,im X(i,k,1,m) = & (X(i,k,1,m)*H(i,k) + A1(i,k-1)*X(i,k-1,1,m))*Haa(i) enddo enddo do m=1,ndet do i=1,im det(i,k,1,m) = & (det(i,k,1,m)*H(i,k) + A1(i,k-1)*det(i,k-1,1,m))*Haa(i) enddo enddo enddo !k if ( .NOT. extmode ) then do i=1,im T(i,km) = E1(i,km)*Tbottom(i) + T(i,km) S(i,km) = Es(i,km)*Sdeep + S(i,km) enddo do m=1,nx do i=1,im X(i,km,1,m) = E1(i,km)*Xdeep(m) + X(i,km,1,m) enddo enddo do m=1,ndet do i=1,im det(i,km,1,m) = E1(i,km)*detdeep(m) + det(i,km,1,m) enddo enddo endif ! .NOT. extmode do k=km-1,1,-1 do m=1,nx do i=1,im X(i,k,1,m) = E1(i,k)*X(i,k+1,1,m) + X(i,k,1,m) enddo enddo enddo do k=km-1,1,-1 do m=1,ndet do i=1,im det(i,k,1,m) = E1(i,k)*det(i,k+1,1,m) + det(i,k,1,m) enddo enddo enddo do k=km-1,1,-1 do i=1,im T(i,k) = E1(i,k)*T(i,k+1) + T(i,k) S(i,k) = Es(i,k)*S(i,k+1) + S(i,k) enddo enddo do k=1,km do i=1,im T(i,k) = POTEMP( T(i,k), S(i,k) , Z(i,k) ) enddo enddo c==> Average viscosity and layer thicknesses C Set Diffusion coefficients do i = 1,idim Wk1(i)=0. Wk2(i)=0. enddo DO k=1,km do i=1,im Wk1(i) = ( Nu(ieast(i),k) + Nu(i,k) )*Area_h(i) Wk2(i) = ( H(ieast(i),k) + H(i,k) )*Area_h(i) enddo ! i do i=1,im Nu(i,k) = (Wk1(i)+Wk1(inorth(i)))*Area_u_r(i) Hu(i,k) = (Wk2(i)+Wk2(inorth(i)))*Area_u_r(i) + ( 1.-Bu(i) ) enddo ! i ENDDO ! k Finish this loop to compute Hu at all k do k=1,km-1 do i=1,im A1(i,k) = Nu(i,k)*dt2df/( Hu(i,k) + Hu(i,k+1) ) enddo enddo do i=1,im ! Note that Nu(*,km) was set at top if extmode A1(i,km) = Nu(i,km)*dt2df/( Hu(i,km) ) enddo k=1 do i=1,im Alpha(i) = A1(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k)) E1(i,k) = A1(i,k) / (Hu(i,k)+A1(i,k)) U(i,k) = U(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k)) V(i,k) = V(i,k)*Hu(i,k) / (Hu(i,k)+A1(i,k)) enddo do k=2,km do i=1,im fact = Bu(i)/(Hu(i,k) + A1(i,k) + Alpha(i)) Alpha(i) = A1(i,k)*(Hu(i,k)+Alpha(i))*fact E1(i,k) = A1(i,k)*fact U(i,k) = (U(i,k)*Hu(i,k) + A1(i,k-1)*U(i,k-1))*fact V(i,k) = (V(i,k)*Hu(i,k) + A1(i,k-1)*V(i,k-1))*fact enddo enddo do k=km-1,1,-1 do i=1,im U(i,k) = E1(i,k)*U(i,k+1) + U(i,k) V(i,k) = E1(i,k)*V(i,k+1) + V(i,k) enddo enddo do k=1,km do i=1,im B(i,k) = BYNCY( T(i,k), S(i,k), 0. ) enddo end do nremain = nremain-1 END DO ! while end