pkg/generic_advdiff/gad_os7mp_adv_y.F

09cd9ab332 Alis*0001 #include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_OS7MP_ADV_Y(
692dd30681 Jean*0004      I           bi,bj,k, calcCFL, deltaTloc,
09cd9ab332 Alis*0005      I           vTrans, vFld,
                0006      I           maskLocS, Q,
                0007      O           vT,
                0008      I           myThid )
                0009 C     /==========================================================\
                0010 C     | SUBROUTINE GAD_OS7MP_ADV_Y                               |
                0011 C     | o Compute Meridional advective Flux of tracer Q using    |
                0012 C     |   7th Order DST Sceheme with monotone preserving limiter |
                0013 C     |==========================================================|
                0014       IMPLICIT NONE
                0015
                0016 C     == GLobal variables ==
                0017 #include "SIZE.h"
                0018 #include "GRID.h"
                0019 #include "GAD.h"
                0020
                0021 C     == Routine arguments ==
                0022       INTEGER bi,bj,k
692dd30681 Jean*0023       LOGICAL calcCFL
09cd9ab332 Alis*0024       _RL deltaTloc
                0025       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0026       _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0027       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0028       _RL Q     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0029       _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0030       INTEGER myThid
                0031
                0032 C     == Local variables ==
                0033       INTEGER i,j
                0034       _RL cfl,Psi
d9afd45777 Alis*0035       _RL vLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
46896f866d Jean*0036       _RL recip_DelIp, recip_DelI
09cd9ab332 Alis*0037       _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
                0038       _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
                0039       _RL d2,d2p1,d2m1,A,B,C,D
46896f866d Jean*0040       _RL dp1h,dm1h, PhiMD,PhiLC,PhiMin,PhiMax
d9afd45777 Alis*0041       _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
                0042       _RL Del2MM,Del2M,Del2,Del2P,Del2PP
                0043       _RL Del3MM,Del3M,Del3P,Del3PP
                0044       _RL Del4M,Del4,Del4P
                0045       _RL Del5M,Del5P
                0046       _RL Del6
09cd9ab332 Alis*0047
                0048       Eps = 1. _d -20
                0049
b4daa24319 Shre*0050       DO i=1-OLx,sNx+OLx
                0051        vT(i,1-OLy)=0. _d 0
                0052        vT(i,2-OLy)=0. _d 0
                0053        vT(i,3-OLy)=0. _d 0
                0054        vT(i,4-OLy)=0. _d 0
                0055        vT(i,sNy+OLy-2)=0. _d 0
                0056        vT(i,sNy+OLy-1)=0. _d 0
                0057        vT(i,sNy+OLy)=0. _d 0
09cd9ab332 Alis*0058       ENDDO
b4daa24319 Shre*0059       DO j=1-OLy+4,sNy+OLy-3
                0060        DO i=1-OLx,sNx+OLx
09cd9ab332 Alis*0061
                0062         vLoc = vFld(i,j)
692dd30681 Jean*0063         cfl = vLoc
                0064         IF ( calcCFL ) cfl = abs(vLoc*deltaTloc*recip_dyC(i,j,bi,bj))
09cd9ab332 Alis*0065
46896f866d Jean*0066         IF (vTrans(i,j).GT.0. _d 0) THEN
09cd9ab332 Alis*0067          Qippp = Q(i,j+2)
                0068          Qipp  = Q(i,j+1)
                0069          Qip   = Q(i,j)
                0070          Qi    = Q(i,j-1)
                0071          Qim   = Q(i,j-2)
                0072          Qimm  = Q(i,j-3)
                0073          Qimmm = Q(i,j-4)
                0074
                0075          MskIpp  = maskLocS(i,j+2)
                0076          MskIp   = maskLocS(i,j+1)
                0077          MskI    = maskLocS(i,j)
                0078          MskIm   = maskLocS(i,j-1)
                0079          MskImm  = maskLocS(i,j-2)
                0080          MskImmm = maskLocS(i,j-3)
46896f866d Jean*0081         ELSEIF (vTrans(i,j).LT.0. _d 0) THEN
09cd9ab332 Alis*0082          Qippp = Q(i,j-3)
                0083          Qipp  = Q(i,j-2)
                0084          Qip   = Q(i,j-1)
                0085          Qi    = Q(i,j)
                0086          Qim   = Q(i,j+1)
                0087          Qimm  = Q(i,j+2)
                0088          Qimmm = Q(i,j+3)
                0089
                0090          MskIpp  = maskLocS(i,j-2)
                0091          MskIp   = maskLocS(i,j-1)
                0092          MskI    = maskLocS(i,j)
                0093          MskIm   = maskLocS(i,j+1)
                0094          MskImm  = maskLocS(i,j+2)
                0095          MskImmm = maskLocS(i,j+3)
fbb39f7b1d Mart*0096         ELSE
                0097          Qippp = 0. _d 0
                0098          Qipp  = 0. _d 0
                0099          Qip   = 0. _d 0
                0100          Qi    = 0. _d 0
                0101          Qim   = 0. _d 0
                0102          Qimm  = 0. _d 0
                0103          Qimmm = 0. _d 0
                0104
                0105          MskIpp  = 0. _d 0
                0106          MskIp   = 0. _d 0
                0107          MskI    = 0. _d 0
                0108          MskIm   = 0. _d 0
                0109          MskImm  = 0. _d 0
                0110          MskImmm = 0. _d 0
09cd9ab332 Alis*0111         ENDIF
                0112
46896f866d Jean*0113         IF (vTrans(i,j).NE.0. _d 0) THEN
09cd9ab332 Alis*0114 C        2nd order correction [i i-1]
46896f866d Jean*0115          Fac = 1. _d 0
d9afd45777 Alis*0116          DelP = (Qip-Qi)*MskI
                0117          Phi = Fac * DelP
09cd9ab332 Alis*0118 C        3rd order correction [i i-1 i-2]
46896f866d Jean*0119          Fac = Fac * ( cfl + 1. _d 0 )/3. _d 0
d9afd45777 Alis*0120          DelM = (Qi-Qim)*MskIm
                0121          Del2 = DelP - DelM
                0122          Phi = Phi - Fac * Del2
09cd9ab332 Alis*0123 C        4th order correction [i+1 i i-1 i-2]
46896f866d Jean*0124          Fac = Fac * ( cfl - 2. _d 0 )/4. _d 0
d9afd45777 Alis*0125          DelPP = (Qipp-Qip)*MskIp*MskI
                0126          Del2P = DelPP - DelP
                0127          Del3P = Del2P - Del2
                0128          Phi = Phi + Fac * Del3p
09cd9ab332 Alis*0129 C        5th order correction [i+1 i i-1 i-2 i-3]
46896f866d Jean*0130          Fac = Fac * ( cfl - 3. _d 0 )/5. _d 0
d9afd45777 Alis*0131          DelMM = (Qim-Qimm)*MskImm*MskIm
                0132          Del2M = DelM - DelMM
                0133          Del3M = Del2 - Del2M
                0134          Del4 = Del3P - Del3M
                0135          Phi = Phi + Fac * Del4
09cd9ab332 Alis*0136 C        6th order correction [i+2 i+1 i i-1 i-2 i-3]
46896f866d Jean*0137          Fac = Fac * ( cfl + 2. _d 0 )/6. _d 0
d9afd45777 Alis*0138          DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
                0139          Del2PP = DelPP - DelP
                0140          Del3PP = Del2PP - Del2P
                0141          Del4P = Del3PP - Del3P
                0142          Del5P = Del4P - Del4
                0143          Phi = Phi + Fac * Del5P
09cd9ab332 Alis*0144 C        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
46896f866d Jean*0145          Fac = Fac * ( cfl + 2. _d 0 )/7. _d 0
d9afd45777 Alis*0146          DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
                0147          Del2MM = DelMM - DelMMM
                0148          Del3MM = Del2M - Del2MM
                0149          Del4M = Del3M - Del3MM
                0150          Del5M = Del4 - Del4M
                0151          Del6 = Del5P - Del5M
                0152          Phi = Phi - Fac * Del6
09cd9ab332 Alis*0153
                0154          DelIp = ( Qip - Qi ) * MskI
46896f866d Jean*0155 c        Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
                0156 c    &        *abs(Phi+Eps)/abs(DelIp+Eps)
                0157 C--   simplify and avoid division by zero
                0158          recip_DelIp = sign(1. _d 0,DelIp)/max(abs(DelIp),Eps)
                0159          Phi = Phi*recip_DelIp
09cd9ab332 Alis*0160
                0161          DelI = ( Qi - Qim ) * MskIm
46896f866d Jean*0162 c        rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
                0163 c    &        *abs(DelI+Eps)/abs(DelIp+Eps)
                0164 C--   simplify and avoid division by zero
                0165          recip_DelI = sign(1. _d 0,DelI)/max(abs(DelI),Eps)
                0166          rp1h = DelI*recip_DelIp
baa188100f Alis*0167          rp1h_cfl = rp1h/(cfl+Eps)
09cd9ab332 Alis*0168
                0169 C        TVD limiter
46896f866d Jean*0170 c        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )
09cd9ab332 Alis*0171
                0172 C        MP limiter
d9afd45777 Alis*0173          d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
                0174          d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
                0175          d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
46896f866d Jean*0176          A = 4. _d 0*d2 - d2p1
                0177          B = 4. _d 0*d2p1 - d2
09cd9ab332 Alis*0178          C = d2
209a396f41 Jean*0179          D = d2p1
b4daa24319 Shre*0180          dp1h = max(min(min(A,B),min(C,D)),0. _d 0)
                0181      &        + min(max(max(A,B),max(C,D)),0. _d 0)
46896f866d Jean*0182          A = 4. _d 0*d2m1 - d2
                0183          B = 4. _d 0*d2 - d2m1
09cd9ab332 Alis*0184          C = d2m1
209a396f41 Jean*0185          D = d2
b4daa24319 Shre*0186          dm1h = max(min(min(A,B),min(C,D)),0. _d 0)
                0187      &        + min(max(max(A,B),max(C,D)),0. _d 0)
46896f866d Jean*0188 c        qMD = 0.5*( ( Qi + Qip ) - dp1h )
                0189 c        qMD = 0.5 _d 0*( ( 2. _d 0*Qi + DelIp ) - dp1h )
                0190 c        qUL = Qi + (1. _d 0-cfl)/(cfl+Eps)*DelI
                0191 c        qLC = Qi + 0.5 _d 0*( 1. _d 0+dm1h/(DelI+Eps) )*(qUL-Qi)
                0192 c        PhiMD = 2. _d 0/(1. _d 0-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
                0193 c        PhiLC = 2. _d 0*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
                0194 C--   simplify and avoid division by zero
                0195          PhiMD = 1. _d 0/(1. _d 0-cfl)*(DelIp-dp1h)*recip_DelIp
                0196          PhiLC = rp1h_cfl*( 1. _d 0+dm1h*recip_DelI )
                0197 C--
                0198          PhiMin = max( min(0. _d 0,PhiMD),
b4daa24319 Shre*0199      &                 min(min(0. _d 0,2. _d 0*rp1h_cfl),PhiLC) )
46896f866d Jean*0200          PhiMax = min( max(2. _d 0/(1. _d 0-cfl),PhiMD),
b4daa24319 Shre*0201      &                 max(max(0. _d 0,2. _d 0*rp1h_cfl),PhiLC) )
09cd9ab332 Alis*0202          Phi = max(PhiMin,min(Phi,PhiMax))
                0203
46896f866d Jean*0204          Psi = Phi * 0.5 _d 0 * (1. _d 0 - cfl)
09cd9ab332 Alis*0205          vT(i,j) = vTrans(i,j)*( Qi + Psi*DelIp )
                0206         ELSE
46896f866d Jean*0207          vT(i,j) = 0. _d 0
09cd9ab332 Alis*0208         ENDIF
                0209
                0210        ENDDO
                0211       ENDDO
                0212
                0213       RETURN
                0214       END