pkg/generic_advdiff/gad_os7mp_adv_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_r.F,v 1.3 2007/05/11 17:30:06 adcroft Exp $
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,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
      _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
      _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
      _RL Del2MM,Del2M,Del2,Del2P,Del2PP
      _RL Del3MM,Del3M,Del3P,Del3PP
      _RL Del4M,Del4,Del4P
      _RL Del5M,Del5P
      _RL Del6

      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 (wTrans(i,j).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 (wTrans(i,j).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 (wTrans(i,j).ne.0.) THEN
C        2nd order correction [i i-1]
         Fac = 1.
         DelP = (Qip-Qi)*MskI
         Phi = Fac * DelP
C        3rd order correction [i i-1 i-2]
         Fac = Fac * ( cfl + 1. )/3.
         DelM = (Qi-Qim)*MskIm
         Del2 = DelP - DelM
         Phi = Phi - Fac * Del2
C        4th order correction [i+1 i i-1 i-2]
         Fac = Fac * ( cfl - 2. )/4.
         DelPP = (Qipp-Qip)*MskIp*MskI
         Del2P = DelPP - DelP
         Del3P = Del2P - Del2
         Phi = Phi + Fac * Del3p
C        5th order correction [i+1 i i-1 i-2 i-3]
         Fac = Fac * ( cfl - 3. )/5.
         DelMM = (Qim-Qimm)*MskImm*MskIm
         Del2M = DelM - DelMM
         Del3M = Del2 - Del2M
         Del4 = Del3P - Del3M
         Phi = Phi + Fac * Del4
C        6th order correction [i+2 i+1 i i-1 i-2 i-3]
         Fac = Fac * ( cfl + 2. )/6.
         DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
         Del2PP = DelPP - DelP
         Del3PP = Del2PP - Del2P
         Del4P = Del3PP - Del3P
         Del5P = Del4P - Del4
         Phi = Phi + Fac * Del5P
C        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
         Fac = Fac * ( cfl + 2. )/7.
         DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
         Del2MM = DelMM - DelMMM
         Del3MM = Del2M - Del2MM
         Del4M = Del3M - Del3MM
         Del5M = Del4 - Del4M
         Del6 = Del5P - Del5M
         Phi = Phi - Fac * Del6

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

         DelI = ( Qi - Qim ) * MskIm
         rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
     &        *abs(DelI+Eps)/abs(DelIp+Eps)
         rp1h_cfl = rp1h/(cfl+Eps)

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

C        MP limiter
         d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
         d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
         d2m1 = Del2M !( ( 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. _d 0)+min(max(A,B,C,D),0. _d 0)
         A = 4.*d2m1 - d2
         B = 4.*d2 - d2m1
         C = d2m1
         D = d2;
         dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
        !qMD = 0.5*( ( Qi + Qip ) - dp1h )
         qMD = 0.5*( ( 2.*Qi + DelIp ) - dp1h )
         qUL = Qi + (1.-cfl)/(cfl+Eps)*DelI
         qLC = Qi + 0.5*( 1.+dm1h/(DelI+Eps) )*(qUL-Qi)
         PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
         PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
         PhiMin = max(min(0. _d 0,PhiMD),
     &        min(0. _d 0,2.*rp1h_cfl,PhiLC))
         PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),
     &        max(0. _d 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.4	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_r.F,v 1.3 2007/05/11 17:30:06 adcroft Exp $
2	mlosch	1.2	C $Name: $
3	adcroft	1.1
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	adcroft	1.4	_RL wLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
37	adcroft	1.1	_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	adcroft	1.4	_RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
42			_RL Del2MM,Del2M,Del2,Del2P,Del2PP
43			_RL Del3MM,Del3M,Del3P,Del3PP
44			_RL Del4M,Del4,Del4P
45			_RL Del5M,Del5P
46			_RL Del6
47	adcroft	1.1
48			Eps = 1. _d -20
49
50			km4=MAX(1,k-4)
51			km3=MAX(1,k-3)
52			km2=MAX(1,k-2)
53			km1=MAX(1,k-1)
54			kp1=MIN(Nr,k+1)
55			kp2=MIN(Nr,k+2)
56			kp3=MIN(Nr,k+3)
57
58			DO j=1-Oly,sNy+Oly
59			DO i=1-Olx,sNx+Olx
60
61			wLoc = wFld(i,j)
62			cfl = abs(wLocdeltaTlocrecip_drC(k))
63
64	adcroft	1.3	IF (wTrans(i,j).lt.0.) THEN
65	adcroft	1.1	Qippp = Q(i,j,kp2)
66			Qipp = Q(i,j,kp1)
67			Qip = Q(i,j,k)
68			Qi = Q(i,j,km1)
69			Qim = Q(i,j,km2)
70			Qimm = Q(i,j,km3)
71			Qimmm = Q(i,j,km4)
72
73			MskIpp = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
74			MskIp = maskC(i,j,kp1,bi,bj) * float(kp1-k)
75			MskI = maskC(i,j,k,bi,bj) * float(k-km1)
76			MskIm = maskC(i,j,km1,bi,bj) * float(km1-km2)
77			MskImm = maskC(i,j,km2,bi,bj) * float(km2-km3)
78			MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)
79	adcroft	1.3	ELSEIF (wTrans(i,j).gt.0.) THEN
80	adcroft	1.1	Qippp = Q(i,j,km3)
81			Qipp = Q(i,j,km2)
82			Qip = Q(i,j,km1)
83			Qi = Q(i,j,k)
84			Qim = Q(i,j,kp1)
85			Qimm = Q(i,j,kp2)
86			Qimmm = Q(i,j,kp3)
87
88			MskIpp = maskC(i,j,km2,bi,bj) * float(km2-km3)
89			MskIp = maskC(i,j,km1,bi,bj) * float(km1-km2)
90			MskI = maskC(i,j,k,bi,bj) * float(k-km1)
91			MskIm = maskC(i,j,kp1,bi,bj) * float(kp1-k)
92			MskImm = maskC(i,j,kp2,bi,bj) * float(kp2-kp1)
93			MskImmm = maskC(i,j,kp3,bi,bj) * float(kp3-kp2)
94			ENDIF
95
96	adcroft	1.3	IF (wTrans(i,j).ne.0.) THEN
97	adcroft	1.1	C 2nd order correction [i i-1]
98			Fac = 1.
99	adcroft	1.4	DelP = (Qip-Qi)*MskI
100			Phi = Fac * DelP
101	adcroft	1.1	C 3rd order correction [i i-1 i-2]
102			Fac = Fac * ( cfl + 1. )/3.
103	adcroft	1.4	DelM = (Qi-Qim)*MskIm
104			Del2 = DelP - DelM
105			Phi = Phi - Fac * Del2
106	adcroft	1.1	C 4th order correction [i+1 i i-1 i-2]
107			Fac = Fac * ( cfl - 2. )/4.
108	adcroft	1.4	DelPP = (Qipp-Qip)MskIpMskI
109			Del2P = DelPP - DelP
110			Del3P = Del2P - Del2
111			Phi = Phi + Fac * Del3p
112	adcroft	1.1	C 5th order correction [i+1 i i-1 i-2 i-3]
113			Fac = Fac * ( cfl - 3. )/5.
114	adcroft	1.4	DelMM = (Qim-Qimm)MskImmMskIm
115			Del2M = DelM - DelMM
116			Del3M = Del2 - Del2M
117			Del4 = Del3P - Del3M
118			Phi = Phi + Fac * Del4
119	adcroft	1.1	C 6th order correction [i+2 i+1 i i-1 i-2 i-3]
120			Fac = Fac * ( cfl + 2. )/6.
121	adcroft	1.4	DelPPP = (Qippp-Qipp)MskIppMskIp*MskI
122			Del2PP = DelPP - DelP
123			Del3PP = Del2PP - Del2P
124			Del4P = Del3PP - Del3P
125			Del5P = Del4P - Del4
126			Phi = Phi + Fac * Del5P
127	adcroft	1.1	C 7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
128			Fac = Fac * ( cfl + 2. )/7.
129	adcroft	1.4	DelMMM = (Qimm-Qimmm)MskImmmMskImm*MskIm
130			Del2MM = DelMM - DelMMM
131			Del3MM = Del2M - Del2MM
132			Del4M = Del3M - Del3MM
133			Del5M = Del4 - Del4M
134			Del6 = Del5P - Del5M
135			Phi = Phi - Fac * Del6
136	adcroft	1.1
137			DelIp = ( Qip - Qi ) * MskI
138	mlosch	1.2	Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
139			& *abs(Phi+Eps)/abs(DelIp+Eps)
140	adcroft	1.1
141			DelI = ( Qi - Qim ) * MskIm
142	mlosch	1.2	rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
143			& *abs(DelI+Eps)/abs(DelIp+Eps)
144	adcroft	1.3	rp1h_cfl = rp1h/(cfl+Eps)
145	adcroft	1.1
146			C TVD limiter
147	adcroft	1.3	! Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )
148	adcroft	1.1
149			C MP limiter
150	adcroft	1.4	d2 = Del2 !( ( Qip + Qim ) - 2.Qi ) MskI * MskIm
151			d2p1 = Del2P !( ( Qipp + Qi ) - 2.Qip ) MskIp * MskI
152			d2m1 = Del2M !( ( Qi + Qimm ) - 2.Qim ) MskIm * MskImm
153	adcroft	1.1	A = 4.*d2 - d2p1
154			B = 4.*d2p1 - d2
155			C = d2
156			D = d2p1;
157	mlosch	1.2	dp1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
158	adcroft	1.1	A = 4.*d2m1 - d2
159			B = 4.*d2 - d2m1
160			C = d2m1
161			D = d2;
162	mlosch	1.2	dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
163	adcroft	1.4	!qMD = 0.5*( ( Qi + Qip ) - dp1h )
164			qMD = 0.5( ( 2.Qi + DelIp ) - dp1h )
165			qUL = Qi + (1.-cfl)/(cfl+Eps)*DelI
166			qLC = Qi + 0.5( 1.+dm1h/(DelI+Eps) )(qUL-Qi)
167			PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
168	adcroft	1.3	PhiLC = 2.rp1h_cfl(qLC-Qi+Eps)/(qUL-Qi+Eps)
169	mlosch	1.2	PhiMin = max(min(0. _d 0,PhiMD),
170	adcroft	1.3	& min(0. _d 0,2.*rp1h_cfl,PhiLC))
171	mlosch	1.2	PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),
172	adcroft	1.3	& max(0. _d 0,2.*rp1h_cfl,PhiLC))
173	adcroft	1.1	Phi = max(PhiMin,min(Phi,PhiMax))
174
175			Psi = Phi * 0.5 * (1. - cfl)
176			wT(i,j) = wTrans(i,j)( Qi + PsiDelIp )
177			ELSE
178			wT(i,j) = 0.
179			ENDIF
180
181			ENDDO
182			ENDDO
183
184			RETURN
185			END