pkg/generic_advdiff/gad_c4_adv_r.F

31566b6684 Alis*0001 #include "GAD_OPTIONS.h"
                0002
527a84022c Alis*0003 CBOP
                0004 C !ROUTINE: GAD_C4_ADV_R
                0005
                0006 C !INTERFACE: ==========================================================
e9de1d7682 Jean*0007       SUBROUTINE GAD_C4_ADV_R(
                0008      I           bi, bj, k,
31566b6684 Alis*0009      I           rTrans,
                0010      I           tracer,
                0011      O           wT,
                0012      I           myThid )
                0013
527a84022c Alis*0014 C !DESCRIPTION:
                0015 C Calculates the area integrated vertical flux due to advection of a tracer
                0016 C using centered fourth-order interpolation:
                0017 C \begin{equation*}
                0018 C F^r_{adv} = W \overline{ \theta - \frac{1}{6} \delta_{kk} \theta }^k
                0019 C \end{equation*}
                0020 C Near boundaries, the scheme reduces to a second if the flow is away
                0021 C from the boundary and to third order if the flow is towards
                0022 C the boundary.
                0023
                0024 C !USES: ===============================================================
                0025       IMPLICIT NONE
31566b6684 Alis*0026 #include "SIZE.h"
                0027 #include "GRID.h"
                0028 #include "GAD.h"
                0029
527a84022c Alis*0030 C !INPUT PARAMETERS: ===================================================
                0031 C  bi,bj                :: tile indices
                0032 C  k                    :: vertical level
                0033 C  rTrans               :: vertical volume transport
                0034 C  tracer               :: tracer field
                0035 C  myThid               :: thread number
31566b6684 Alis*0036       INTEGER bi,bj,k
                0037       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e9de1d7682 Jean*0038       _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31566b6684 Alis*0039       INTEGER myThid
                0040
527a84022c Alis*0041 C !OUTPUT PARAMETERS: ==================================================
                0042 C  wT                   :: vertical advective flux
                0043       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044
                0045 C !LOCAL VARIABLES: ====================================================
                0046 C  i,j                  :: loop indices
                0047 C  kp1                  :: =min( k+1 , Nr )
                0048 C  km1                  :: =max( k-1 , 1 )
                0049 C  km2                  :: =max( k-2 , 1 )
                0050 C  Rjm,Rj,Rjp           :: differences at i-1,i,i+1
                0051 C  Rjjm,Rjjp            :: second differences at i-1,i
                0052 C  maskP1               :: =1 for k<Nr, =0 for k>=Nr
31566b6684 Alis*0053       INTEGER i,j,kp1,km1,km2
8a37624a25 Jean*0054       _RL Rjm,Rj,Rjp,Rjjm,Rjjp
5d63e10665 Jean*0055       _RL maskPM, maskBound
527a84022c Alis*0056 CEOP
31566b6684 Alis*0057
                0058       km2=MAX(1,k-2)
                0059       km1=MAX(1,k-1)
                0060       kp1=MIN(Nr,k+1)
5d63e10665 Jean*0061       maskPM = 1.
                0062       IF (k.LE.2 .OR. k.GE.Nr) maskPM = 0.
31566b6684 Alis*0063
ab2f187c7d Jean*0064       IF ( k.EQ.1 .OR. k.GT.Nr) THEN
e9de1d7682 Jean*0065        DO j=1-OLy,sNy+OLy
                0066         DO i=1-OLx,sNx+OLx
31566b6684 Alis*0067          wT(i,j) = 0.
                0068         ENDDO
                0069        ENDDO
                0070       ELSE
e9de1d7682 Jean*0071        DO j=1-OLy,sNy+OLy
                0072         DO i=1-OLx,sNx+OLx
5d63e10665 Jean*0073          maskBound = maskPM*maskC(i,j,km2,bi,bj)*maskC(i,j,kp1,bi,bj)
e9de1d7682 Jean*0074          Rjp = (tracer(i,j,kp1)-tracer(i,j, k ))*maskC(i,j,kp1,bi,bj)
                0075          Rj  = (tracer(i,j, k )-tracer(i,j,km1))
                0076          Rjm = (tracer(i,j,km1)-tracer(i,j,km2))*maskC(i,j,km1,bi,bj)
                0077      &
5d63e10665 Jean*0078          Rjjp=(Rjp-Rj)
                0079          Rjjm=(Rj-Rjm)
8a37624a25 Jean*0080          wT(i,j) = maskC(i,j,km1,bi,bj)*(
e9de1d7682 Jean*0081      &     rTrans(i,j)*( (tracer(i,j,k)+tracer(i,j,km1))*0.5 _d 0
                0082      &                  -oneSixth*(Rjjm+Rjjp)*0.5 _d 0 )
8a37624a25 Jean*0083      &    +ABS(rTrans(i,j))*
e9de1d7682 Jean*0084      &                   oneSixth*(Rjjm-Rjjp)*0.5 _d 0
                0085      &                     *(1. _d 0 - maskBound)
                0086      &                                  )
31566b6684 Alis*0087         ENDDO
                0088        ENDDO
                0089       ENDIF
                0090
                0091       RETURN
                0092       END