pkg/generic_advdiff/gad_os7mp_adv_r.F

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