pkg/generic_advdiff/gad_biharm_r.F

d76a0ac2c0 Jean*0001 #include "GAD_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: GAD_BIHARM_R
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE GAD_BIHARM_R(
                0008      I           bi, bj, k,
0e9223926c Jean*0009      I           maskUp, diffKr4, tracer,
d76a0ac2c0 Jean*0010      U           d4f,
                0011      I           myThid )
                0012
                0013 C !DESCRIPTION:
                0014 C Calculates the vertical flux due to bi-harmonic diffusion of a tracer.
                0015 C \begin{equation*}
                0016 C F^r_{diff} = \kappa_4 \partial_r ( \partial^2 / \partial_r^2 \theta )
                0017 C \end{equation*}
                0018
                0019 C !USES: ===============================================================
                0020       IMPLICIT NONE
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025
                0026 C !INPUT PARAMETERS: ===================================================
                0027 C  bi,bj                :: tile indices
                0028 C  k                    :: vertical level
0e9223926c Jean*0029 C  maskUp               :: 2-D array for mask at W points
d76a0ac2c0 Jean*0030 C  diffKr4              :: vertical bi-harmonic diffusivity
                0031 C  tracer               :: tracer field
                0032 C  myThid               :: my thread Id number
                0033       INTEGER bi, bj, k
0e9223926c Jean*0034       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d76a0ac2c0 Jean*0035       _RL diffKr4(Nr)
e9de1d7682 Jean*0036       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
d76a0ac2c0 Jean*0037       INTEGER myThid
                0038
                0039 C !OUTPUT PARAMETERS: ==================================================
                0040 C  dfy                  :: vertical diffusive flux
                0041       _RL d4f    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042
                0043 C !LOCAL VARIABLES: ====================================================
                0044 C  i,j,n                :: loop indices
                0045 C  del2T                :: vertical 2nd  derivative of tracer
                0046 C  gradR                :: vertical 1rst derivative of tracer (*rkSign)
                0047       INTEGER i, j, n
                0048       INTEGER kl, km, kp
                0049       _RL del2T(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:2)
                0050       _RL gradR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:3)
                0051       _RL tmpFac
                0052 CEOP
                0053
                0054       IF ( k.GE.2 ) THEN
                0055
                0056 C     Calculate vertical gradient @ interface k-1, k & k+1 (n=1:3)
                0057        DO n=1,3
                0058         km = k+n-3
                0059         kl = k+n-2
                0060         IF ( km.LT.1 .OR. kl.GT.Nr ) THEN
                0061          DO j=1-OLy,sNy+OLy
                0062           DO i=1-OLx,sNx+OLx
                0063            gradR(i,j,n) = 0.
                0064           ENDDO
                0065          ENDDO
                0066         ELSE
                0067          tmpFac = recip_drC(kl)*deepFac2F(kl)*rhoFacF(kl)
                0068          DO j=1-OLy,sNy+OLy
                0069           DO i=1-OLx,sNx+OLx
e9de1d7682 Jean*0070            gradR(i,j,n) = ( tracer(i,j,kl)-tracer(i,j,km) )
d76a0ac2c0 Jean*0071      &                *tmpFac*maskC(i,j,kl,bi,bj)*maskC(i,j,km,bi,bj)
                0072           ENDDO
                0073          ENDDO
                0074         ENDIF
                0075        ENDDO
                0076
                0077 C     Calculate vertical 2nd derivative @ level k-1 & k (n=1:2)
                0078        DO n=1,2
                0079         kl = k+n-2
                0080         kp = k+n-1
                0081         tmpFac = recip_drF(kl)*recip_deepFac2C(kl)*recip_rhoFacC(kl)
                0082         DO j=1-OLy,sNy+OLy
                0083          DO i=1-OLx,sNx+OLx
                0084            del2T(i,j,n) = ( gradR(i,j,n+1)-gradR(i,j,n) )
                0085      &                  *_recip_hFacC(i,j,kl,bi,bj)
                0086          ENDDO
                0087         ENDDO
                0088        ENDDO
                0089
                0090 C     Add bi-harmonic flux (gradient of 2nd derivative)
                0091        tmpFac = rkSign*recip_drC(k)*deepFac2F(k)*rhoFacF(k)
                0092        DO j=1-OLy,sNy+OLy
                0093         DO i=1-OLx,sNx+OLx
                0094           d4f(i,j) = d4f(i,j)
                0095      &             + diffKr4(k)*( del2T(i,j,2)-del2T(i,j,1) )
0e9223926c Jean*0096      &              *tmpFac*_rA(i,j,bi,bj)*maskUp(i,j)
d76a0ac2c0 Jean*0097         ENDDO
                0098        ENDDO
                0099
                0100       ENDIF
                0101
                0102       RETURN
                0103       END