pkg/generic_advdiff/gad_os7mp_adv_y.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_y.F,v 1.4 2007/05/11 17:30:06 adcroft Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

      SUBROUTINE GAD_OS7MP_ADV_Y(
     I           bi,bj,k, calcCFL, deltaTloc,
     I           vTrans, vFld,
     I           maskLocS, Q,
     O           vT,
     I           myThid )
C     /==========================================================\
C     | SUBROUTINE GAD_OS7MP_ADV_Y                               |
C     | o Compute Meridional 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 vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL Q     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

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

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

        vLoc = vFld(i,j)
        cfl = vLoc
        IF ( calcCFL ) cfl = abs(vLoc*deltaTloc*recip_dyC(i,j,bi,bj))

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

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

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

        IF (vTrans(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)
         vT(i,j) = vTrans(i,j)*( Qi + Psi*DelIp )
        ELSE
         vT(i,j) = 0.
        ENDIF

       ENDDO
      ENDDO

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