C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/jscott/code_rafmod/gad_diff_x.F,v 1.2 2009/09/03 20:40:01 jscott Exp $ C $Name: $ #include "GAD_OPTIONS.h" CBOP C !ROUTINE: GAD_DIFF_X C !INTERFACE: ========================================================== SUBROUTINE GAD_DIFF_X( I bi,bj,k, I xA, diffKh, I tracer, O dfx, I myThid ) C !DESCRIPTION: C Calculates the area integrated zonal flux due to down-gradient diffusion C of a tracer: C \begin{equation*} C F^x_{diff} = - A^x \kappa_h \frac{1}{\Delta x_c} \delta_i \theta C \end{equation*} C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "GRID.h" C !INPUT PARAMETERS: =================================================== C bi,bj :: tile indices C k :: vertical level C xA :: area of face at U points C diffKh :: horizontal diffusivity C tracer :: tracer field C myThid :: thread number INTEGER bi,bj,k _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL diffKh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid C !OUTPUT PARAMETERS: ================================================== C dfx :: zonal diffusive flux _RL dfx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices INTEGER i,j CEOP C print *,'in gad_diff_x' DO j=1-Oly,sNy+Oly dfx(1-Olx,j)=0. DO i=1-Olx+1,sNx+Olx dfx(i,j) = -diffKh(i,j,k)*xA(i,j) & *_recip_dxC(i,j,bi,bj) & *(Tracer(i,j)-Tracer(i-1,j)) & *CosFacU(j,bi,bj) ENDDO ENDDO RETURN END c---------------------------------------------------------- SUBROUTINE CALC_MLD(bi,bj,diffkh3d_x,diffkh3d_y,myThid) implicit none #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" INTEGER bi,bj _RL diffkh3d_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL diffkh3d_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER myThid C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices INTEGER i,j !,k c jrs next _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL Khplus,gradij,gradim1,gradjm1,calc_grad C _RL dep_frac real *8 khconst,khmin,khmax parameter ( khconst=5.d15 ) !2.5d15 parameter ( khmin=1000.0 ) parameter ( khmax=40000.0 ) !15000.0 cjrs try constants 6e14, 3e15, 5e15, 10e15 c jrs next C print *,'in calc_mld for Khplus ' Cold CALL FIND_RHO( C & bi,bj,1-Olx,sNx+OLx,1-OLy,sNy+OLy,1,1, C & theta,salt,rhoK,myThid) CALL FIND_RHO_2D( & 1-Olx,sNx+OLx,1-OLy,sNy+OLy,1, & theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj), & rhoK, 1, bi, bj, myThid) DO j=1-Oly,sNy+Oly DO i=1-Olx+1,sNx+Olx if (hfacc(i,j,1,bi,bj).ne.0) then gradij= calc_grad(i,j,bi,bj,rhoK) endif if ( (hfacc(i,j,1,bi,bj).ne.0).and. & (hfacc(i-1,j,1,bi,bj).ne.0) ) then gradim1= 0.! calc_grad(i-1,j,bi,bj,rhoK) if ((gradij.ne.0).and.(gradim1.ne.0)) then Khplus= khconst*(gradij+gradim1)/2 elseif (gradij.ne.0) then Khplus= khconst*gradij elseif (gradim1.ne.0) then Khplus= khconst*gradim1 else Khplus=0.d0 endif if (Khplus.lt.khmin) Khplus=khmin if (Khplus.gt.khmax) Khplus=khmax if (yC(i,j,bi,bj) .GE. 80.) Khplus=khmin diffkh3d_x(i,j,1)=Khplus C print *,'Khplus_x: ',i,j,Khplus else diffkh3d_x(i,j,1)=0.d0 endif if ( (hfacc(i,j,1,bi,bj).ne.0).and. & (hfacc(i,j-1,1,bi,bj).ne.0) ) then gradjm1=0.!calc_grad(i,j-1,bi,bj,rhoK) if ((gradij.ne.0).and.(gradjm1.ne.0)) then Khplus= khconst*(gradij+gradjm1)/2 elseif (gradij.ne.0) then Khplus= khconst*gradij elseif (gradjm1.ne.0) then Khplus= khconst*gradjm1 else Khplus=0.d0 endif if (Khplus.lt.khmin) Khplus=khmin if (Khplus.gt.khmax) Khplus=khmax if (yC(i,j,bi,bj).GE. 80.) Khplus=khmin diffkh3d_y(i,j,1)=Khplus c diffkh3d_y(i,j,2)=0.5*Khplus c diffkh3d_y(i,j,3)=0.2*Khplus C print *,'Khplus_y: ',i,j,Khplus else diffkh3d_y(i,j,1)=0.d0 endif c DO k=2,7 c if (-hMixLayer(i,j,bi,bj).ge.rF(k)) then c diffkh3d_x(i,j,k)=0.d0 c diffkh3d_y(i,j,k)=0.d0 c elseif (-hMixLayer(i,j,bi,bj).ge.rF(k+1)) then c dep_frac=(hMixLayer(i,j,bi,bj)+rF(k))/drF(k)* c & (240.+rC(k))/240. c diffkh3d_x(i,j,k)= dep_frac*diffkh3d_x(i,j,1) c diffkh3d_y(i,j,k)= dep_frac*diffkh3d_y(i,j,1) c else c dep_frac = (240.+rC(k))/240. c diffkh3d_x(i,j,k)= dep_frac*diffkh3d_x(i,j,1) c diffkh3d_y(i,j,k)= dep_frac*diffkh3d_y(i,j,1) c endif c ENDDO ENDDO ENDDO return end c-------------------------------------------------- real*8 FUNCTION CALC_GRAD(i,j,bi,bj,rhoK) implicit none #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" INTEGER i,j,bi,bj _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy) c Local variables integer xno,yno _RL xgrad,ygrad c print *,'entering calc_gradient...for ',i,j if (hfacC(i,j,1,bi,bj).eq.0) then calc_grad=0.d0 else xgrad=0.d0 xno=0 if (hfacC(i+1,j,1,bi,bj).ne.0) then xgrad= abs((rhoK(i+1,j)-rhoK(i,j))*_recip_dxC(i+1,j,bi,bj)) xno=1 c print *,'gradient with i+1: ',xgrad endif if (hfacC(i-1,j,1,bi,bj).ne.0) then xgrad= xgrad + & abs((rhoK(i,j)-rhoK(i-1,j))*_recip_dxC(i,j,bi,bj)) xno=xno+1 c print *,'gradient with i-1: ', c & abs((rhoK(i,j)-rhoK(i-1,j))*_recip_dxC(i,j,bi,bj)) endif ygrad=0.d0 yno=0 if (hfacC(i,j+1,1,bi,bj).ne.0) then ygrad= abs((rhoK(i,j+1)-rhoK(i,j))*_recip_dyC(i,j+1,bi,bj)) yno=1 endif if (hfacC(i,j-1,1,bi,bj).ne.0) then ygrad= ygrad + & abs((rhoK(i,j)-rhoK(i,j-1))*_recip_dyC(i,j,bi,bj)) yno=yno+1 endif if ((xno.ne.0).and.(yno.ne.0)) then calc_grad=(xgrad/xno)*(xgrad/xno) + & (ygrad/yno)*(ygrad/yno) elseif (xno.ne.0) then calc_grad= (xgrad/xno)*(xgrad/xno) else calc_grad= (ygrad/yno)*(ygrad/yno) endif endif return end