pkg/generic_advdiff/gad_os7mp_adv_y.F

C $Header: $
C $Name: $

#include "GAD_OPTIONS.h"

      SUBROUTINE GAD_OS7MP_ADV_Y(
     I           bi,bj,k,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
      _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,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 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 = abs(vLoc*deltaTloc*recip_dyC(i,j,bi,bj))

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

       ENDDO
      ENDDO

      RETURN
      END
1	C $Header: $
2	C $Name: $
3
4	#include "GAD_OPTIONS.h"
5
6	SUBROUTINE GAD_OS7MP_ADV_Y(
7	I bi,bj,k,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	_RL deltaTloc
27	_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28	_RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29	_RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30	_RL Q (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31	_RL vT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32	INTEGER myThid
33
34	C == Local variables ==
35	INTEGER i,j
36	_RL cfl,Psi
37	_RL vLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,Msk
38	_RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
39	_RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
40	_RL d2,d2p1,d2m1,A,B,C,D
41	_RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax
42
43	Eps = 1. _d -20
44
45	DO i=1-Olx,sNx+Olx
46	vT(i,1-Oly)=0. _d 0
47	vT(i,2-Oly)=0. _d 0
48	vT(i,3-Oly)=0. _d 0
49	vT(i,4-Oly)=0. _d 0
50	vT(i,sNy+Oly-2)=0. _d 0
51	vT(i,sNy+Oly-1)=0. _d 0
52	vT(i,sNy+Oly)=0. _d 0
53	ENDDO
54	DO j=1-Oly+4,sNy+Oly-3
55	DO i=1-Olx,sNx+Olx
56
57	vLoc = vFld(i,j)
58	cfl = abs(vLocdeltaTlocrecip_dyC(i,j,bi,bj))
59
60	IF (vLoc.gt.0.) THEN
61	Qippp = Q(i,j+2)
62	Qipp = Q(i,j+1)
63	Qip = Q(i,j)
64	Qi = Q(i,j-1)
65	Qim = Q(i,j-2)
66	Qimm = Q(i,j-3)
67	Qimmm = Q(i,j-4)
68
69	MskIpp = maskLocS(i,j+2)
70	MskIp = maskLocS(i,j+1)
71	MskI = maskLocS(i,j)
72	MskIm = maskLocS(i,j-1)
73	MskImm = maskLocS(i,j-2)
74	MskImmm = maskLocS(i,j-3)
75	ELSEIF (vLoc.lt.0.) THEN
76	Qippp = Q(i,j-3)
77	Qipp = Q(i,j-2)
78	Qip = Q(i,j-1)
79	Qi = Q(i,j)
80	Qim = Q(i,j+1)
81	Qimm = Q(i,j+2)
82	Qimmm = Q(i,j+3)
83
84	MskIpp = maskLocS(i,j-2)
85	MskIp = maskLocS(i,j-1)
86	MskI = maskLocS(i,j)
87	MskIm = maskLocS(i,j+1)
88	MskImm = maskLocS(i,j+2)
89	MskImmm = maskLocS(i,j+3)
90	ENDIF
91
92	IF (vLoc.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	Phi = sign(1.,Phi)sign(1.,DelIp)abs(Phi+Eps)/abs(DelIp+Eps)
126
127	DelI = ( Qi - Qim ) * MskIm
128	rp1h =sign(1.,DelI)sign(1.,DelIp)abs(DelI+Eps)/abs(DelIp+Eps)
129
130	C TVD limiter
131	! Phi = max(0., min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )
132
133	C MP limiter
134	d2 = ( ( Qip + Qim ) - 2.Qi ) MskI * MskIm
135	d2p1 = ( ( Qipp + Qi ) - 2.Qip ) MskIp * MskI
136	d2m1 = ( ( Qi + Qimm ) - 2.Qim ) MskIm * MskImm
137	A = 4.*d2 - d2p1
138	B = 4.*d2p1 - d2
139	C = d2
140	D = d2p1;
141	dp1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
142	A = 4.*d2m1 - d2
143	B = 4.*d2 - d2m1
144	C = d2m1
145	D = d2;
146	dm1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)
147	qMD = 0.5*( ( Qi + Qip ) - dp1h )
148	qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )
149	qLC = Qi + 0.5( 1.+dm1h/(Qi-Qim+Eps) )(qUL-Qi)
150	PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
151	PhiLC = 2.rp1h/cfl(qLC-Qi+Eps)/(qUL-Qi+Eps)
152	PhiMin = max(min(0.,PhiMD),min(0.,2.*rp1h/cfl,PhiLC))
153	PhiMax = min(max(2./(1.-cfl),PhiMD),max(0.,2.*rp1h/cfl,PhiLC))
154	Phi = max(PhiMin,min(Phi,PhiMax))
155
156	Psi = Phi * 0.5 * (1. - cfl)
157	vT(i,j) = vTrans(i,j)( Qi + PsiDelIp )
158	ELSE
159	vt(i,j) = 0.
160	ENDIF
161
162	ENDDO
163	ENDDO
164
165	RETURN
166	END