pkg/generic_advdiff/gad_os7mp_adv_x.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_x.F,v 1.2 2007/01/21 17:25:31 mlosch Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

      SUBROUTINE GAD_OS7MP_ADV_X(
     I           bi,bj,k, calcCFL, deltaTloc,
     I           uTrans, uFld,
     I           maskLocW, Q,
     O           uT,
     I           myThid )
C     /==========================================================\
C     | SUBROUTINE GAD_OS7MP_ADV_X                               |
C     | o Compute Zonal 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
      LOGICAL calcCFL
      _RL deltaTloc
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL Q     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

C     == Local variables ==
      INTEGER i,j
      _RL cfl,Psi
      _RL uLoc,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

      DO j=1-Oly,sNy+Oly
       uT(1-Olx,j)=0. _d 0
       uT(2-Olx,j)=0. _d 0
       uT(3-Olx,j)=0. _d 0
       uT(4-Olx,j)=0. _d 0
       uT(sNx+Olx-2,j)=0. _d 0
       uT(sNx+Olx-1,j)=0. _d 0
       uT(sNx+Olx,j)=0. _d 0
       DO i=1-Olx+4,sNx+Olx-3

        uLoc = uFld(i,j)
        cfl = uLoc
        IF ( calcCFL ) cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))

        IF (uLoc.gt.0.) THEN
         Qippp = Q(i+2,j)
         Qipp  = Q(i+1,j)
         Qip   = Q(i,j)
         Qi    = Q(i-1,j)
         Qim   = Q(i-2,j)
         Qimm  = Q(i-3,j)
         Qimmm = Q(i-4,j)

         MskIpp  = maskLocW(i+2,j)
         MskIp   = maskLocW(i+1,j)
         MskI    = maskLocW(i,j)
         MskIm   = maskLocW(i-1,j)
         MskImm  = maskLocW(i-2,j)
         MskImmm = maskLocW(i-3,j)
        ELSEIF (uLoc.lt.0.) THEN
         Qippp = Q(i-3,j)
         Qipp  = Q(i-2,j)
         Qip   = Q(i-1,j)
         Qi    = Q(i,j)
         Qim   = Q(i+1,j)
         Qimm  = Q(i+2,j)
         Qimmm = Q(i+3,j)

         MskIpp  = maskLocW(i-2,j)
         MskIp   = maskLocW(i-1,j)
         MskI    = maskLocW(i,j)
         MskIm   = maskLocW(i+1,j)
         MskImm  = maskLocW(i+2,j)
         MskImmm = maskLocW(i+3,j)
        ENDIF

        IF (uLoc.ne.0.) THEN
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. _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)

C        TVD limiter
!        Phi = max(0. _d 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. _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 )
         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. _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)
         uT(i,j) = uTrans(i,j)*( Qi + Psi*DelIp )
        ELSE
         uT(i,j) = 0.
        ENDIF

       ENDDO
      ENDDO

      RETURN
      END
1	jmc	1.3	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_x.F,v 1.2 2007/01/21 17:25:31 mlosch Exp $
2	mlosch	1.2	C $Name: $
3	adcroft	1.1
4			#include "GAD_OPTIONS.h"
5
6			SUBROUTINE GAD_OS7MP_ADV_X(
7	jmc	1.3	I bi,bj,k, calcCFL, deltaTloc,
8	adcroft	1.1	I uTrans, uFld,
9			I maskLocW, Q,
10			O uT,
11			I myThid )
12			C /==========================================================\
13			C \| SUBROUTINE GAD_OS7MP_ADV_X \|
14			C \| o Compute Zonal 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	jmc	1.3	LOGICAL calcCFL
27	adcroft	1.1	_RL deltaTloc
28			_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29			_RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30			_RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31			_RL Q (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32			_RL uT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33			INTEGER myThid
34
35			C == Local variables ==
36			INTEGER i,j
37			_RL cfl,Psi
38			_RL uLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,Msk
39			_RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
40			_RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
41			_RL d2,d2p1,d2m1,A,B,C,D
42			_RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax
43
44			Eps = 1. _d -20
45
46			DO j=1-Oly,sNy+Oly
47			uT(1-Olx,j)=0. _d 0
48			uT(2-Olx,j)=0. _d 0
49			uT(3-Olx,j)=0. _d 0
50			uT(4-Olx,j)=0. _d 0
51			uT(sNx+Olx-2,j)=0. _d 0
52			uT(sNx+Olx-1,j)=0. _d 0
53			uT(sNx+Olx,j)=0. _d 0
54			DO i=1-Olx+4,sNx+Olx-3
55
56			uLoc = uFld(i,j)
57	jmc	1.3	cfl = uLoc
58			IF ( calcCFL ) cfl = abs(uLocdeltaTlocrecip_dxC(i,j,bi,bj))
59	adcroft	1.1
60			IF (uLoc.gt.0.) THEN
61			Qippp = Q(i+2,j)
62			Qipp = Q(i+1,j)
63			Qip = Q(i,j)
64			Qi = Q(i-1,j)
65			Qim = Q(i-2,j)
66			Qimm = Q(i-3,j)
67			Qimmm = Q(i-4,j)
68
69			MskIpp = maskLocW(i+2,j)
70			MskIp = maskLocW(i+1,j)
71			MskI = maskLocW(i,j)
72			MskIm = maskLocW(i-1,j)
73			MskImm = maskLocW(i-2,j)
74			MskImmm = maskLocW(i-3,j)
75			ELSEIF (uLoc.lt.0.) THEN
76			Qippp = Q(i-3,j)
77			Qipp = Q(i-2,j)
78			Qip = Q(i-1,j)
79			Qi = Q(i,j)
80			Qim = Q(i+1,j)
81			Qimm = Q(i+2,j)
82			Qimmm = Q(i+3,j)
83
84			MskIpp = maskLocW(i-2,j)
85			MskIp = maskLocW(i-1,j)
86			MskI = maskLocW(i,j)
87			MskIm = maskLocW(i+1,j)
88			MskImm = maskLocW(i+2,j)
89			MskImmm = maskLocW(i+3,j)
90			ENDIF
91
92			IF (uLoc.ne.0.) THEN
93			C 2nd order correction [i i-1]
94			Fac = 1.
95			Del = Qip-Qi
96			Msk = MskI
97			Phi = Msk * Fac * Del
98			C 3rd order correction [i i-1 i-2]
99			Fac = Fac * ( cfl + 1. )/3.
100			Del = Del - ( Qi-Qim )
101			Msk = Msk * MskIm
102			Phi = Phi - Msk * Fac * Del
103			C 4th order correction [i+1 i i-1 i-2]
104			Fac = Fac * ( cfl - 2. )/4.
105			Del = ( Qipp-2.*Qip+Qi ) - Del
106			Msk = Msk * MskIp
107			Phi = Phi + Msk * Fac * Del
108			C 5th order correction [i+1 i i-1 i-2 i-3]
109			Fac = Fac * ( cfl - 3. )/5.
110			Del = Del - ( Qip-3.Qi+3.Qim-Qimm )
111			Msk = Msk * MskImm
112			Phi = Phi + Msk * Fac * Del
113			C 6th order correction [i+2 i+1 i i-1 i-2 i-3]
114			Fac = Fac * ( cfl + 2. )/6.
115			Del = ( Qippp-4.Qipp+6.Qip-4.*Qi+Qim ) - Del
116			Msk = Msk * MskIpp
117			Phi = Phi + Msk * Fac * Del
118			C 7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
119			Fac = Fac * ( cfl + 2. )/7.
120			Del = Del - ( Qipp-5.Qip+10.Qi-10.Qim+5.Qimm-Qimmm )
121			Msk = Msk * MskImmm
122			Phi = Phi - Msk * Fac * Del
123
124			DelIp = ( Qip - Qi ) * MskI
125	mlosch	1.2	Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
126			& *abs(Phi+Eps)/abs(DelIp+Eps)
127	adcroft	1.1
128			DelI = ( Qi - Qim ) * MskIm
129	mlosch	1.2	rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
130			& *abs(DelI+Eps)/abs(DelIp+Eps)
131	adcroft	1.1
132			C TVD limiter
133	mlosch	1.2	! Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )
134	adcroft	1.1
135			C MP limiter
136			d2 = ( ( Qip + Qim ) - 2.Qi ) MskI * MskIm
137			d2p1 = ( ( Qipp + Qi ) - 2.Qip ) MskIp * MskI
138			d2m1 = ( ( Qi + Qimm ) - 2.Qim ) MskIm * MskImm
139			A = 4.*d2 - d2p1
140			B = 4.*d2p1 - d2
141			C = d2
142			D = d2p1;
143	mlosch	1.2	dp1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
144	adcroft	1.1	A = 4.*d2m1 - d2
145			B = 4.*d2 - d2m1
146			C = d2m1
147			D = d2;
148	mlosch	1.2	dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
149	adcroft	1.1	qMD = 0.5*( ( Qi + Qip ) - dp1h )
150			qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )
151			qLC = Qi + 0.5( 1.+dm1h/(Qi-Qim+Eps) )(qUL-Qi)
152			PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
153			PhiLC = 2.rp1h/cfl(qLC-Qi+Eps)/(qUL-Qi+Eps)
154	mlosch	1.2	PhiMin = max(min(0. _d 0,PhiMD),
155			& min(0. _d 0,2.*rp1h/cfl,PhiLC))
156			PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),
157			& max(0. _d 0,2.*rp1h/cfl,PhiLC))
158	adcroft	1.1	Phi = max(PhiMin,min(Phi,PhiMax))
159
160			Psi = Phi * 0.5 * (1. - cfl)
161			uT(i,j) = uTrans(i,j)( Qi + PsiDelIp )
162			ELSE
163			uT(i,j) = 0.
164			ENDIF
165
166			ENDDO
167			ENDDO
168
169			RETURN
170			END