/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_os7mp_adv_x.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_os7mp_adv_x.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by mlosch, Sun Jan 21 17:25:31 2007 UTC revision 1.9 by jmc, Mon Jun 16 13:40:25 2008 UTC
# Line 4  C $Name$ Line 4  C $Name$
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_OS7MP_ADV_X(        SUBROUTINE GAD_OS7MP_ADV_X(
7       I           bi,bj,k,deltaTloc,       I           bi,bj,k, calcCFL, deltaTloc,
8       I           uTrans, uFld,       I           uTrans, uFld,
9       I           maskLocW, Q,       I           maskLocW, Q,
10       O           uT,       O           uT,
# Line 23  C     == GLobal variables == Line 23  C     == GLobal variables ==
23    
24  C     == Routine arguments ==  C     == Routine arguments ==
25        INTEGER bi,bj,k        INTEGER bi,bj,k
26          LOGICAL calcCFL
27        _RL deltaTloc        _RL deltaTloc
28        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29        _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 34  C     == Routine arguments == Line 35  C     == Routine arguments ==
35  C     == Local variables ==  C     == Local variables ==
36        INTEGER i,j        INTEGER i,j
37        _RL cfl,Psi        _RL cfl,Psi
38        _RL uLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,Msk        _RL uLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
39          _RL recip_DelIp, recip_DelI
40        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
41        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
42        _RL d2,d2p1,d2m1,A,B,C,D        _RL d2,d2p1,d2m1,A,B,C,D
43        _RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax        _RL dp1h,dm1h, PhiMD,PhiLC,PhiMin,PhiMax
44          _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
45          _RL Del2MM,Del2M,Del2,Del2P,Del2PP
46          _RL Del3MM,Del3M,Del3P,Del3PP
47          _RL Del4M,Del4,Del4P
48          _RL Del5M,Del5P
49          _RL Del6
50    
51        Eps = 1. _d -20        Eps = 1. _d -20
52    
# Line 50  C     == Local variables == Line 58  C     == Local variables ==
58         uT(sNx+Olx-2,j)=0. _d 0         uT(sNx+Olx-2,j)=0. _d 0
59         uT(sNx+Olx-1,j)=0. _d 0         uT(sNx+Olx-1,j)=0. _d 0
60         uT(sNx+Olx,j)=0. _d 0         uT(sNx+Olx,j)=0. _d 0
61          ENDDO
62          DO j=1-Oly,sNy+Oly
63         DO i=1-Olx+4,sNx+Olx-3         DO i=1-Olx+4,sNx+Olx-3
64    
65          uLoc = uFld(i,j)          uLoc = uFld(i,j)
66          cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))          cfl = uLoc
67            IF ( calcCFL ) cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))
68    
69          IF (uLoc.gt.0.) THEN          IF (uTrans(i,j).GT.0. _d 0) THEN
70           Qippp = Q(i+2,j)           Qippp = Q(i+2,j)
71           Qipp  = Q(i+1,j)           Qipp  = Q(i+1,j)
72           Qip   = Q(i,j)           Qip   = Q(i,j)
# Line 70  C     == Local variables == Line 81  C     == Local variables ==
81           MskIm   = maskLocW(i-1,j)           MskIm   = maskLocW(i-1,j)
82           MskImm  = maskLocW(i-2,j)           MskImm  = maskLocW(i-2,j)
83           MskImmm = maskLocW(i-3,j)           MskImmm = maskLocW(i-3,j)
84          ELSEIF (uLoc.lt.0.) THEN          ELSEIF (uTrans(i,j).LT.0. _d 0) THEN
85           Qippp = Q(i-3,j)           Qippp = Q(i-3,j)
86           Qipp  = Q(i-2,j)           Qipp  = Q(i-2,j)
87           Qip   = Q(i-1,j)           Qip   = Q(i-1,j)
# Line 85  C     == Local variables == Line 96  C     == Local variables ==
96           MskIm   = maskLocW(i+1,j)           MskIm   = maskLocW(i+1,j)
97           MskImm  = maskLocW(i+2,j)           MskImm  = maskLocW(i+2,j)
98           MskImmm = maskLocW(i+3,j)           MskImmm = maskLocW(i+3,j)
99            ELSE
100             Qippp = 0. _d 0
101             Qipp  = 0. _d 0
102             Qip   = 0. _d 0
103             Qi    = 0. _d 0
104             Qim   = 0. _d 0
105             Qimm  = 0. _d 0
106             Qimmm = 0. _d 0
107    
108             MskIpp  = 0. _d 0
109             MskIp   = 0. _d 0
110             MskI    = 0. _d 0
111             MskIm   = 0. _d 0
112             MskImm  = 0. _d 0
113             MskImmm = 0. _d 0
114          ENDIF          ENDIF
115    
116          IF (uLoc.ne.0.) THEN          IF (uTrans(i,j).NE.0. _d 0) THEN
117  C        2nd order correction [i i-1]  C        2nd order correction [i i-1]
118           Fac = 1.           Fac = 1. _d 0
119           Del = Qip-Qi           DelP = (Qip-Qi)*MskI
120           Msk = MskI           Phi = Fac * DelP
          Phi = Msk * Fac * Del  
121  C        3rd order correction [i i-1 i-2]  C        3rd order correction [i i-1 i-2]
122           Fac = Fac * ( cfl + 1. )/3.           Fac = Fac * ( cfl + 1. _d 0 )/3. _d 0
123           Del = Del - ( Qi-Qim )           DelM = (Qi-Qim)*MskIm
124           Msk = Msk * MskIm           Del2 = DelP - DelM
125           Phi = Phi - Msk * Fac * Del           Phi = Phi - Fac * Del2
126  C        4th order correction [i+1 i i-1 i-2]  C        4th order correction [i+1 i i-1 i-2]
127           Fac = Fac * ( cfl - 2. )/4.           Fac = Fac * ( cfl - 2. _d 0 )/4. _d 0
128           Del = ( Qipp-2.*Qip+Qi ) - Del           DelPP = (Qipp-Qip)*MskIp*MskI
129           Msk = Msk * MskIp           Del2P = DelPP - DelP
130           Phi = Phi + Msk * Fac * Del           Del3P = Del2P - Del2
131             Phi = Phi + Fac * Del3p
132  C        5th order correction [i+1 i i-1 i-2 i-3]  C        5th order correction [i+1 i i-1 i-2 i-3]
133           Fac = Fac * ( cfl - 3. )/5.           Fac = Fac * ( cfl - 3. _d 0 )/5. _d 0
134           Del = Del - ( Qip-3.*Qi+3.*Qim-Qimm )           DelMM = (Qim-Qimm)*MskImm*MskIm
135           Msk = Msk * MskImm           Del2M = DelM - DelMM
136           Phi = Phi + Msk * Fac * Del           Del3M = Del2 - Del2M
137             Del4 = Del3P - Del3M
138             Phi = Phi + Fac * Del4
139  C        6th order correction [i+2 i+1 i i-1 i-2 i-3]  C        6th order correction [i+2 i+1 i i-1 i-2 i-3]
140           Fac = Fac * ( cfl + 2. )/6.           Fac = Fac * ( cfl + 2. _d 0 )/6. _d 0
141           Del = ( Qippp-4.*Qipp+6.*Qip-4.*Qi+Qim ) - Del           DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
142           Msk = Msk * MskIpp           Del2PP = DelPP - DelP
143           Phi = Phi + Msk * Fac * Del           Del3PP = Del2PP - Del2P
144             Del4P = Del3PP - Del3P
145             Del5P = Del4P - Del4
146             Phi = Phi + Fac * Del5P
147  C        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]  C        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
148           Fac = Fac * ( cfl + 2. )/7.           Fac = Fac * ( cfl + 2. _d 0 )/7. _d 0
149           Del = Del - ( Qipp-5.*Qip+10.*Qi-10.*Qim+5.*Qimm-Qimmm )           DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
150           Msk = Msk * MskImmm           Del2MM = DelMM - DelMMM
151           Phi = Phi - Msk * Fac * Del           Del3MM = Del2M - Del2MM
152             Del4M = Del3M - Del3MM
153             Del5M = Del4 - Del4M
154             Del6 = Del5P - Del5M
155             Phi = Phi - Fac * Del6
156    
157           DelIp = ( Qip - Qi ) * MskI           DelIp = ( Qip - Qi ) * MskI
158           Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)  c        Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
159       &        *abs(Phi+Eps)/abs(DelIp+Eps)  c    &        *abs(Phi+Eps)/abs(DelIp+Eps)
160    C--   simplify and avoid division by zero
161             recip_DelIp = sign(1. _d 0,DelIp)/max(abs(DelIp),Eps)
162             Phi = Phi*recip_DelIp
163    
164           DelI = ( Qi - Qim ) * MskIm           DelI = ( Qi - Qim ) * MskIm
165           rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)  c        rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
166       &        *abs(DelI+Eps)/abs(DelIp+Eps)  c    &        *abs(DelI+Eps)/abs(DelIp+Eps)
167    C--   simplify and avoid division by zero
168             recip_DelI = sign(1. _d 0,DelI)/max(abs(DelI),Eps)
169             rp1h = DelI*recip_DelIp
170             rp1h_cfl = rp1h/(cfl+Eps)
171    
172  C        TVD limiter  C        TVD limiter
173  !        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )  c        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )
174    
175  C        MP limiter  C        MP limiter
176           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm           d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
177           d2p1 = ( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI           d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
178           d2m1 = ( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm           d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
179           A = 4.*d2 - d2p1           A = 4. _d 0*d2 - d2p1
180           B = 4.*d2p1 - d2           B = 4. _d 0*d2p1 - d2
181           C = d2           C = d2
182           D = d2p1;           D = d2p1
183           dp1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)           dp1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
184           A = 4.*d2m1 - d2           A = 4. _d 0*d2m1 - d2
185           B = 4.*d2 - d2m1           B = 4. _d 0*d2 - d2m1
186           C = d2m1           C = d2m1
187           D = d2;           D = d2
188           dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)           dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
189           qMD = 0.5*( ( Qi + Qip ) - dp1h )  c        qMD = 0.5*( ( Qi + Qip ) - dp1h )
190           qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )  c        qMD = 0.5 _d 0*( ( 2. _d 0*Qi + DelIp ) - dp1h )
191           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)  c        qUL = Qi + (1. _d 0-cfl)/(cfl+Eps)*DelI
192           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)  c        qLC = Qi + 0.5 _d 0*( 1. _d 0+dm1h/(DelI+Eps) )*(qUL-Qi)
193           PhiLC = 2.*rp1h/cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)  c        PhiMD = 2. _d 0/(1. _d 0-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
194           PhiMin = max(min(0. _d 0,PhiMD),  c        PhiLC = 2. _d 0*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
195       &        min(0. _d 0,2.*rp1h/cfl,PhiLC))  C--   simplify and avoid division by zero
196           PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),           PhiMD = 1. _d 0/(1. _d 0-cfl)*(DelIp-dp1h)*recip_DelIp
197       &        max(0. _d 0,2.*rp1h/cfl,PhiLC))           PhiLC = rp1h_cfl*( 1. _d 0+dm1h*recip_DelI )
198    C--
199             PhiMin = max( min(0. _d 0,PhiMD),
200         &                 min(0. _d 0,2. _d 0*rp1h_cfl,PhiLC) )
201             PhiMax = min( max(2. _d 0/(1. _d 0-cfl),PhiMD),
202         &                 max(0. _d 0,2. _d 0*rp1h_cfl,PhiLC) )
203           Phi = max(PhiMin,min(Phi,PhiMax))           Phi = max(PhiMin,min(Phi,PhiMax))
204    
205           Psi = Phi * 0.5 * (1. - cfl)           Psi = Phi * 0.5 _d 0 * (1. _d 0 - cfl)
206           uT(i,j) = uTrans(i,j)*( Qi + Psi*DelIp )           uT(i,j) = uTrans(i,j)*( Qi + Psi*DelIp )
207          ELSE          ELSE
208           uT(i,j) = 0.           uT(i,j) = 0. _d 0
209          ENDIF          ENDIF
210    
211         ENDDO         ENDDO
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22