pkg/generic_advdiff/gad_os7mp_adv_x.F

09cd9ab332 Alis*0001 #include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_OS7MP_ADV_X(
692dd30681 Jean*0004      I           bi,bj,k, calcCFL, deltaTloc,
09cd9ab332 Alis*0005      I           uTrans, uFld,
                0006      I           maskLocW, Q,
                0007      O           uT,
                0008      I           myThid )
                0009 C     /==========================================================\
                0010 C     | SUBROUTINE GAD_OS7MP_ADV_X                               |
                0011 C     | o Compute Zonal 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 uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0026       _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0027       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0028       _RL Q     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0029       _RL uT    (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 uLoc,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 j=1-OLy,sNy+OLy
                0051        uT(1-OLx,j)=0. _d 0
                0052        uT(2-OLx,j)=0. _d 0
                0053        uT(3-OLx,j)=0. _d 0
                0054        uT(4-OLx,j)=0. _d 0
                0055        uT(sNx+OLx-2,j)=0. _d 0
                0056        uT(sNx+OLx-1,j)=0. _d 0
                0057        uT(sNx+OLx,j)=0. _d 0
360ad14abb Mart*0058       ENDDO
b4daa24319 Shre*0059       DO j=1-OLy,sNy+OLy
                0060        DO i=1-OLx+4,sNx+OLx-3
09cd9ab332 Alis*0061
                0062         uLoc = uFld(i,j)
692dd30681 Jean*0063         cfl = uLoc
                0064         IF ( calcCFL ) cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))
09cd9ab332 Alis*0065
46896f866d Jean*0066         IF (uTrans(i,j).GT.0. _d 0) THEN
09cd9ab332 Alis*0067          Qippp = Q(i+2,j)
                0068          Qipp  = Q(i+1,j)
                0069          Qip   = Q(i,j)
                0070          Qi    = Q(i-1,j)
                0071          Qim   = Q(i-2,j)
                0072          Qimm  = Q(i-3,j)
                0073          Qimmm = Q(i-4,j)
                0074
                0075          MskIpp  = maskLocW(i+2,j)
                0076          MskIp   = maskLocW(i+1,j)
                0077          MskI    = maskLocW(i,j)
                0078          MskIm   = maskLocW(i-1,j)
                0079          MskImm  = maskLocW(i-2,j)
                0080          MskImmm = maskLocW(i-3,j)
46896f866d Jean*0081         ELSEIF (uTrans(i,j).LT.0. _d 0) THEN
09cd9ab332 Alis*0082          Qippp = Q(i-3,j)
                0083          Qipp  = Q(i-2,j)
                0084          Qip   = Q(i-1,j)
                0085          Qi    = Q(i,j)
                0086          Qim   = Q(i+1,j)
                0087          Qimm  = Q(i+2,j)
                0088          Qimmm = Q(i+3,j)
                0089
                0090          MskIpp  = maskLocW(i-2,j)
                0091          MskIp   = maskLocW(i-1,j)
                0092          MskI    = maskLocW(i,j)
                0093          MskIm   = maskLocW(i+1,j)
                0094          MskImm  = maskLocW(i+2,j)
                0095          MskImmm = maskLocW(i+3,j)
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 (uTrans(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          uT(i,j) = uTrans(i,j)*( Qi + Psi*DelIp )
                0206         ELSE
46896f866d Jean*0207          uT(i,j) = 0. _d 0
09cd9ab332 Alis*0208         ENDIF
                0209
                0210        ENDDO
                0211       ENDDO
                0212
                0213       RETURN
                0214       END