pkg/generic_advdiff/gad_os7mp_adv_r.F

C $Header: $
C $Name: $

#include "GAD_OPTIONS.h"

      SUBROUTINE GAD_OS7MP_ADV_R(
     I           bi,bj,k,deltaTloc,
     I           wTrans, wFld,
     I           Q,
     O           wT,
     I           myThid )
C     /==========================================================\
C     | SUBROUTINE GAD_OS7MP_ADV_R                               |
C     | o Compute Vertical advective Flux of tracer Q using      |
C     |   7th Order DST Sceheme with monotone preserving limiter |
C     |==========================================================|
      IMPLICIT NONE

C     == GLobal variables ==
#include "SIZE.h"
#include "GRID.h"
#include "GAD.h"

C     == Routine arguments ==
      INTEGER bi,bj,k
      _RL deltaTloc
      _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL Q     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

C     == Local variables ==
      INTEGER i,j,kp3,kp2,kp1,km1,km2,km3,km4
      _RL cfl,Psi
      _RL wLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,Msk
      _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
      _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
      _RL d2,d2p1,d2m1,A,B,C,D
      _RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax

      Eps = 1. _d -20

      km4=MAX(1,k-4)
      km3=MAX(1,k-3)
      km2=MAX(1,k-2)
      km1=MAX(1,k-1)
      kp1=MIN(Nr,k+1)
      kp2=MIN(Nr,k+2)
      kp3=MIN(Nr,k+3)

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx

        wLoc = wFld(i,j)
        cfl = abs(wLoc*deltaTloc*recip_drC(k))

        IF (wLoc.lt.0.) THEN
         Qippp = Q(i,j,kp2)
         Qipp  = Q(i,j,kp1)
         Qip   = Q(i,j,k)
         Qi    = Q(i,j,km1)
         Qim   = Q(i,j,km2)
         Qimm  = Q(i,j,km3)
         Qimmm = Q(i,j,km4)

         MskIpp  = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
         MskIp   = maskC(i,j,kp1,bi,bj) * float(kp1-k)
         MskI    = maskC(i,j,k,bi,bj)   * float(k-km1)
         MskIm   = maskC(i,j,km1,bi,bj) * float(km1-km2)
         MskImm  = maskC(i,j,km2,bi,bj) * float(km2-km3)
         MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)
        ELSEIF (wLoc.gt.0.) THEN
         Qippp = Q(i,j,km3)
         Qipp  = Q(i,j,km2)
         Qip   = Q(i,j,km1)
         Qi    = Q(i,j,k)
         Qim   = Q(i,j,kp1)
         Qimm  = Q(i,j,kp2)
         Qimmm = Q(i,j,kp3)

         MskIpp  = maskC(i,j,km2,bi,bj) * float(km2-km3)
         MskIp   = maskC(i,j,km1,bi,bj) * float(km1-km2)
         MskI    = maskC(i,j,k,bi,bj)   * float(k-km1)
         MskIm   = maskC(i,j,kp1,bi,bj) * float(kp1-k)
         MskImm  = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
         MskImmm = maskC(i,j,kp3,bi,bj) * float(kp3-kp2)
        ENDIF

        IF (wLoc.ne.0.) THEN
         Phi = 0.
C        2nd order correction [i i-1]
         Fac = 1.
         Del = Qip-Qi
         Msk = MskI
         Phi = Msk * Fac * Del
C        3rd order correction [i i-1 i-2]
         Fac = Fac * ( cfl + 1. )/3.
         Del = Del - ( Qi-Qim )
         Msk = Msk * MskIm
         Phi = Phi - Msk * Fac * Del
C        4th order correction [i+1 i i-1 i-2]
         Fac = Fac * ( cfl - 2. )/4.
         Del = ( Qipp-2.*Qip+Qi ) - Del
         Msk = Msk * MskIp
         Phi = Phi + Msk * Fac * Del
C        5th order correction [i+1 i i-1 i-2 i-3]
         Fac = Fac * ( cfl - 3. )/5.
         Del = Del - ( Qip-3.*Qi+3.*Qim-Qimm )
         Msk = Msk * MskImm
         Phi = Phi + Msk * Fac * Del
C        6th order correction [i+2 i+1 i i-1 i-2 i-3]
         Fac = Fac * ( cfl + 2. )/6.
         Del = ( Qippp-4.*Qipp+6.*Qip-4.*Qi+Qim ) - Del
         Msk = Msk * MskIpp
         Phi = Phi + Msk * Fac * Del
C        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
         Fac = Fac * ( cfl + 2. )/7.
         Del = Del - ( Qipp-5.*Qip+10.*Qi-10.*Qim+5.*Qimm-Qimmm ) 
         Msk = Msk * MskImmm
         Phi = Phi - Msk * Fac * Del

         DelIp = ( Qip - Qi ) * MskI
         Phi = sign(1.,Phi)*sign(1.,DelIp)*abs(Phi+Eps)/abs(DelIp+Eps)

         DelI = ( Qi - Qim ) * MskIm
         rp1h =sign(1.,DelI)*sign(1.,DelIp)*abs(DelI+Eps)/abs(DelIp+Eps)

C        TVD limiter
!        Phi = max(0., min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )

C        MP limiter
         d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
         d2p1 = ( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
         d2m1 = ( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
         A = 4.*d2 - d2p1
         B = 4.*d2p1 - d2
         C = d2
         D = d2p1;
         dp1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
         A = 4.*d2m1 - d2
         B = 4.*d2 - d2m1
         C = d2m1
         D = d2;
         dm1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
         qMD = 0.5*( ( Qi + Qip ) - dp1h )
         qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )
         qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)
         PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
         PhiLC = 2.*rp1h/cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
         PhiMin = max(min(0.,PhiMD),min(0.,2.*rp1h/cfl,PhiLC))
         PhiMax = min(max(2./(1.-cfl),PhiMD),max(0.,2.*rp1h/cfl,PhiLC))
         Phi = max(PhiMin,min(Phi,PhiMax))

         Psi = Phi * 0.5 * (1. - cfl)
         wT(i,j) = wTrans(i,j)*( Qi + Psi*DelIp )
        ELSE
         wT(i,j) = 0.
        ENDIF

       ENDDO
      ENDDO

      RETURN
      END
1	adcroft	1.1	C $Header: $
2			C $Name: $
3
4			#include "GAD_OPTIONS.h"
5
6			SUBROUTINE GAD_OS7MP_ADV_R(
7			I bi,bj,k,deltaTloc,
8			I wTrans, wFld,
9			I Q,
10			O wT,
11			I myThid )
12			C /==========================================================\
13			C \| SUBROUTINE GAD_OS7MP_ADV_R \|
14			C \| o Compute Vertical advective Flux of tracer Q using \|
15			C \| 7th Order DST Sceheme with monotone preserving limiter \|
16			C \|==========================================================\|
17			IMPLICIT NONE
18
19			C == GLobal variables ==
20			#include "SIZE.h"
21			#include "GRID.h"
22			#include "GAD.h"
23
24			C == Routine arguments ==
25			INTEGER bi,bj,k
26			_RL deltaTloc
27			_RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28			_RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29			_RL Q (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
30			_RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31			INTEGER myThid
32
33			C == Local variables ==
34			INTEGER i,j,kp3,kp2,kp1,km1,km2,km3,km4
35			_RL cfl,Psi
36			_RL wLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,Msk
37			_RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
38			_RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
39			_RL d2,d2p1,d2m1,A,B,C,D
40			_RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax
41
42			Eps = 1. _d -20
43
44			km4=MAX(1,k-4)
45			km3=MAX(1,k-3)
46			km2=MAX(1,k-2)
47			km1=MAX(1,k-1)
48			kp1=MIN(Nr,k+1)
49			kp2=MIN(Nr,k+2)
50			kp3=MIN(Nr,k+3)
51
52			DO j=1-Oly,sNy+Oly
53			DO i=1-Olx,sNx+Olx
54
55			wLoc = wFld(i,j)
56			cfl = abs(wLocdeltaTlocrecip_drC(k))
57
58			IF (wLoc.lt.0.) THEN
59			Qippp = Q(i,j,kp2)
60			Qipp = Q(i,j,kp1)
61			Qip = Q(i,j,k)
62			Qi = Q(i,j,km1)
63			Qim = Q(i,j,km2)
64			Qimm = Q(i,j,km3)
65			Qimmm = Q(i,j,km4)
66
67			MskIpp = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
68			MskIp = maskC(i,j,kp1,bi,bj) * float(kp1-k)
69			MskI = maskC(i,j,k,bi,bj) * float(k-km1)
70			MskIm = maskC(i,j,km1,bi,bj) * float(km1-km2)
71			MskImm = maskC(i,j,km2,bi,bj) * float(km2-km3)
72			MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)
73			ELSEIF (wLoc.gt.0.) THEN
74			Qippp = Q(i,j,km3)
75			Qipp = Q(i,j,km2)
76			Qip = Q(i,j,km1)
77			Qi = Q(i,j,k)
78			Qim = Q(i,j,kp1)
79			Qimm = Q(i,j,kp2)
80			Qimmm = Q(i,j,kp3)
81
82			MskIpp = maskC(i,j,km2,bi,bj) * float(km2-km3)
83			MskIp = maskC(i,j,km1,bi,bj) * float(km1-km2)
84			MskI = maskC(i,j,k,bi,bj) * float(k-km1)
85			MskIm = maskC(i,j,kp1,bi,bj) * float(kp1-k)
86			MskImm = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
87			MskImmm = maskC(i,j,kp3,bi,bj) * float(kp3-kp2)
88			ENDIF
89
90			IF (wLoc.ne.0.) THEN
91			Phi = 0.
92			C 2nd order correction [i i-1]
93			Fac = 1.
94			Del = Qip-Qi
95			Msk = MskI
96			Phi = Msk * Fac * Del
97			C 3rd order correction [i i-1 i-2]
98			Fac = Fac * ( cfl + 1. )/3.
99			Del = Del - ( Qi-Qim )
100			Msk = Msk * MskIm
101			Phi = Phi - Msk * Fac * Del
102			C 4th order correction [i+1 i i-1 i-2]
103			Fac = Fac * ( cfl - 2. )/4.
104			Del = ( Qipp-2.*Qip+Qi ) - Del
105			Msk = Msk * MskIp
106			Phi = Phi + Msk * Fac * Del
107			C 5th order correction [i+1 i i-1 i-2 i-3]
108			Fac = Fac * ( cfl - 3. )/5.
109			Del = Del - ( Qip-3.Qi+3.Qim-Qimm )
110			Msk = Msk * MskImm
111			Phi = Phi + Msk * Fac * Del
112			C 6th order correction [i+2 i+1 i i-1 i-2 i-3]
113			Fac = Fac * ( cfl + 2. )/6.
114			Del = ( Qippp-4.Qipp+6.Qip-4.*Qi+Qim ) - Del
115			Msk = Msk * MskIpp
116			Phi = Phi + Msk * Fac * Del
117			C 7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
118			Fac = Fac * ( cfl + 2. )/7.
119			Del = Del - ( Qipp-5.Qip+10.Qi-10.Qim+5.Qimm-Qimmm )
120			Msk = Msk * MskImmm
121			Phi = Phi - Msk * Fac * Del
122
123			DelIp = ( Qip - Qi ) * MskI
124			Phi = sign(1.,Phi)sign(1.,DelIp)abs(Phi+Eps)/abs(DelIp+Eps)
125
126			DelI = ( Qi - Qim ) * MskIm
127			rp1h =sign(1.,DelI)sign(1.,DelIp)abs(DelI+Eps)/abs(DelIp+Eps)
128
129			C TVD limiter
130			! Phi = max(0., min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )
131
132			C MP limiter
133			d2 = ( ( Qip + Qim ) - 2.Qi ) MskI * MskIm
134			d2p1 = ( ( Qipp + Qi ) - 2.Qip ) MskIp * MskI
135			d2m1 = ( ( Qi + Qimm ) - 2.Qim ) MskIm * MskImm
136			A = 4.*d2 - d2p1
137			B = 4.*d2p1 - d2
138			C = d2
139			D = d2p1;
140			dp1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
141			A = 4.*d2m1 - d2
142			B = 4.*d2 - d2m1
143			C = d2m1
144			D = d2;
145			dm1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
146			qMD = 0.5*( ( Qi + Qip ) - dp1h )
147			qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )
148			qLC = Qi + 0.5( 1.+dm1h/(Qi-Qim+Eps) )(qUL-Qi)
149			PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
150			PhiLC = 2.rp1h/cfl(qLC-Qi+Eps)/(qUL-Qi+Eps)
151			PhiMin = max(min(0.,PhiMD),min(0.,2.*rp1h/cfl,PhiLC))
152			PhiMax = min(max(2./(1.-cfl),PhiMD),max(0.,2.*rp1h/cfl,PhiLC))
153			Phi = max(PhiMin,min(Phi,PhiMax))
154
155			Psi = Phi * 0.5 * (1. - cfl)
156			wT(i,j) = wTrans(i,j)( Qi + PsiDelIp )
157			ELSE
158			wT(i,j) = 0.
159			ENDIF
160
161			ENDDO
162			ENDDO
163
164			RETURN
165			END