C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_biharm_r.F,v 1.2 2014/08/18 12:22:46 jmc Exp $ C $Name: $ #include "GAD_OPTIONS.h" CBOP C !ROUTINE: GAD_BIHARM_R C !INTERFACE: ========================================================== SUBROUTINE GAD_BIHARM_R( I bi, bj, k, I diffKr4, tracer, U d4f, I myThid ) C !DESCRIPTION: C Calculates the vertical flux due to bi-harmonic diffusion of a tracer. C \begin{equation*} C F^r_{diff} = \kappa_4 \partial_r ( \partial^2 / \partial_r^2 \theta ) C \end{equation*} C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" C !INPUT PARAMETERS: =================================================== C bi,bj :: tile indices C k :: vertical level C diffKr4 :: vertical bi-harmonic diffusivity C tracer :: tracer field C myThid :: my thread Id number INTEGER bi, bj, k _RL diffKr4(Nr) _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) INTEGER myThid C !OUTPUT PARAMETERS: ================================================== C dfy :: vertical diffusive flux _RL d4f (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== C i,j,n :: loop indices C del2T :: vertical 2nd derivative of tracer C gradR :: vertical 1rst derivative of tracer (*rkSign) INTEGER i, j, n INTEGER kl, km, kp _RL del2T(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:2) _RL gradR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:3) _RL tmpFac CEOP IF ( k.GE.2 ) THEN C Calculate vertical gradient @ interface k-1, k & k+1 (n=1:3) DO n=1,3 km = k+n-3 kl = k+n-2 IF ( km.LT.1 .OR. kl.GT.Nr ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx gradR(i,j,n) = 0. ENDDO ENDDO ELSE tmpFac = recip_drC(kl)*deepFac2F(kl)*rhoFacF(kl) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx gradR(i,j,n) = ( tracer(i,j,kl)-tracer(i,j,km) ) & *tmpFac*maskC(i,j,kl,bi,bj)*maskC(i,j,km,bi,bj) ENDDO ENDDO ENDIF ENDDO C Calculate vertical 2nd derivative @ level k-1 & k (n=1:2) DO n=1,2 kl = k+n-2 kp = k+n-1 tmpFac = recip_drF(kl)*recip_deepFac2C(kl)*recip_rhoFacC(kl) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx del2T(i,j,n) = ( gradR(i,j,n+1)-gradR(i,j,n) ) & *_recip_hFacC(i,j,kl,bi,bj) ENDDO ENDDO ENDDO C Add bi-harmonic flux (gradient of 2nd derivative) tmpFac = rkSign*recip_drC(k)*deepFac2F(k)*rhoFacF(k) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx d4f(i,j) = d4f(i,j) & + diffKr4(k)*( del2T(i,j,2)-del2T(i,j,1) ) & *tmpFac*_rA(i,j,bi,bj) ENDDO ENDDO ENDIF RETURN END