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

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

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

revision 1.6 by mlosch, Fri Oct 5 10:50:47 2007 UTC revision 1.7 by jmc, Mon Jun 16 13:40:25 2008 UTC
# Line 34  C     == Local variables == Line 34  C     == Local variables ==
34        INTEGER i,j,kp3,kp2,kp1,km1,km2,km3,km4        INTEGER i,j,kp3,kp2,kp1,km1,km2,km3,km4
35        _RL cfl,Psi        _RL cfl,Psi
36        _RL wLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl        _RL wLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
37          _RL recip_DelIp, recip_DelI
38        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
39        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
40        _RL d2,d2p1,d2m1,A,B,C,D        _RL d2,d2p1,d2m1,A,B,C,D
41        _RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax        _RL dp1h,dm1h, PhiMD,PhiLC,PhiMin,PhiMax
42        _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP        _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
43        _RL Del2MM,Del2M,Del2,Del2P,Del2PP        _RL Del2MM,Del2M,Del2,Del2P,Del2PP
44        _RL Del3MM,Del3M,Del3P,Del3PP        _RL Del3MM,Del3M,Del3P,Del3PP
# Line 61  C     == Local variables == Line 62  C     == Local variables ==
62          wLoc = wFld(i,j)          wLoc = wFld(i,j)
63          cfl = abs(wLoc*deltaTloc*recip_drC(k))          cfl = abs(wLoc*deltaTloc*recip_drC(k))
64    
65          IF (wTrans(i,j).lt.0.) THEN          IF (wTrans(i,j).LT.0. _d 0) THEN
66           Qippp = Q(i,j,kp2)           Qippp = Q(i,j,kp2)
67           Qipp  = Q(i,j,kp1)           Qipp  = Q(i,j,kp1)
68           Qip   = Q(i,j,k)           Qip   = Q(i,j,k)
# Line 76  C     == Local variables == Line 77  C     == Local variables ==
77           MskIm   = maskC(i,j,km1,bi,bj) * float(km1-km2)           MskIm   = maskC(i,j,km1,bi,bj) * float(km1-km2)
78           MskImm  = maskC(i,j,km2,bi,bj) * float(km2-km3)           MskImm  = maskC(i,j,km2,bi,bj) * float(km2-km3)
79           MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)           MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)
80          ELSEIF (wTrans(i,j).gt.0.) THEN          ELSEIF (wTrans(i,j).GT.0. _d 0) THEN
81           Qippp = Q(i,j,km3)           Qippp = Q(i,j,km3)
82           Qipp  = Q(i,j,km2)           Qipp  = Q(i,j,km2)
83           Qip   = Q(i,j,km1)           Qip   = Q(i,j,km1)
# Line 108  C     == Local variables == Line 109  C     == Local variables ==
109           MskImmm = 0. _d 0           MskImmm = 0. _d 0
110          ENDIF          ENDIF
111    
112          IF (wTrans(i,j).ne.0.) THEN          IF (wTrans(i,j).NE.0. _d 0) THEN
113  C        2nd order correction [i i-1]  C        2nd order correction [i i-1]
114           Fac = 1.           Fac = 1. _d 0
115           DelP = (Qip-Qi)*MskI           DelP = (Qip-Qi)*MskI
116           Phi = Fac * DelP           Phi = Fac * DelP
117  C        3rd order correction [i i-1 i-2]  C        3rd order correction [i i-1 i-2]
118           Fac = Fac * ( cfl + 1. )/3.           Fac = Fac * ( cfl + 1. _d 0 )/3. _d 0
119           DelM = (Qi-Qim)*MskIm           DelM = (Qi-Qim)*MskIm
120           Del2 = DelP - DelM           Del2 = DelP - DelM
121           Phi = Phi - Fac * Del2           Phi = Phi - Fac * Del2
122  C        4th order correction [i+1 i i-1 i-2]  C        4th order correction [i+1 i i-1 i-2]
123           Fac = Fac * ( cfl - 2. )/4.           Fac = Fac * ( cfl - 2. _d 0 )/4. _d 0
124           DelPP = (Qipp-Qip)*MskIp*MskI           DelPP = (Qipp-Qip)*MskIp*MskI
125           Del2P = DelPP - DelP           Del2P = DelPP - DelP
126           Del3P = Del2P - Del2           Del3P = Del2P - Del2
127           Phi = Phi + Fac * Del3p           Phi = Phi + Fac * Del3p
128  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]
129           Fac = Fac * ( cfl - 3. )/5.           Fac = Fac * ( cfl - 3. _d 0 )/5. _d 0
130           DelMM = (Qim-Qimm)*MskImm*MskIm           DelMM = (Qim-Qimm)*MskImm*MskIm
131           Del2M = DelM - DelMM           Del2M = DelM - DelMM
132           Del3M = Del2 - Del2M           Del3M = Del2 - Del2M
133           Del4 = Del3P - Del3M           Del4 = Del3P - Del3M
134           Phi = Phi + Fac * Del4           Phi = Phi + Fac * Del4
135  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]
136           Fac = Fac * ( cfl + 2. )/6.           Fac = Fac * ( cfl + 2. _d 0 )/6. _d 0
137           DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI           DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
138           Del2PP = DelPP - DelP           Del2PP = DelPP - DelP
139           Del3PP = Del2PP - Del2P           Del3PP = Del2PP - Del2P
# Line 140  C        6th order correction [i+2 i+1 i Line 141  C        6th order correction [i+2 i+1 i
141           Del5P = Del4P - Del4           Del5P = Del4P - Del4
142           Phi = Phi + Fac * Del5P           Phi = Phi + Fac * Del5P
143  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]
144           Fac = Fac * ( cfl + 2. )/7.           Fac = Fac * ( cfl + 2. _d 0 )/7. _d 0
145           DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm           DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
146           Del2MM = DelMM - DelMMM           Del2MM = DelMM - DelMMM
147           Del3MM = Del2M - Del2MM           Del3MM = Del2M - Del2MM
# Line 150  C        7th order correction [i+2 i+1 i Line 151  C        7th order correction [i+2 i+1 i
151           Phi = Phi - Fac * Del6           Phi = Phi - Fac * Del6
152    
153           DelIp = ( Qip - Qi ) * MskI           DelIp = ( Qip - Qi ) * MskI
154           Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)  c        Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
155       &        *abs(Phi+Eps)/abs(DelIp+Eps)  c    &        *abs(Phi+Eps)/abs(DelIp+Eps)
156    C--   simplify and avoid division by zero
157             recip_DelIp = sign(1. _d 0,DelIp)/max(abs(DelIp),Eps)
158             Phi = Phi*recip_DelIp
159    
160           DelI = ( Qi - Qim ) * MskIm           DelI = ( Qi - Qim ) * MskIm
161           rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)  c        rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
162       &        *abs(DelI+Eps)/abs(DelIp+Eps)  c    &        *abs(DelI+Eps)/abs(DelIp+Eps)
163    C--   simplify and avoid division by zero
164             recip_DelI = sign(1. _d 0,DelI)/max(abs(DelI),Eps)
165             rp1h = DelI*recip_DelIp
166           rp1h_cfl = rp1h/(cfl+Eps)           rp1h_cfl = rp1h/(cfl+Eps)
167    
168  C        TVD limiter  C        TVD limiter
169  !        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 ) )
170    
171  C        MP limiter  C        MP limiter
172           d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm           d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
173           d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI           d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
174           d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm           d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
175           A = 4.*d2 - d2p1           A = 4. _d 0*d2 - d2p1
176           B = 4.*d2p1 - d2           B = 4. _d 0*d2p1 - d2
177           C = d2           C = d2
178           D = d2p1           D = d2p1
179           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)
180           A = 4.*d2m1 - d2           A = 4. _d 0*d2m1 - d2
181           B = 4.*d2 - d2m1           B = 4. _d 0*d2 - d2m1
182           C = d2m1           C = d2m1
183           D = d2           D = d2
184           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)
185          !qMD = 0.5*( ( Qi + Qip ) - dp1h )  c        qMD = 0.5*( ( Qi + Qip ) - dp1h )
186           qMD = 0.5*( ( 2.*Qi + DelIp ) - dp1h )  c        qMD = 0.5 _d 0*( ( 2. _d 0*Qi + DelIp ) - dp1h )
187           qUL = Qi + (1.-cfl)/(cfl+Eps)*DelI  c        qUL = Qi + (1. _d 0-cfl)/(cfl+Eps)*DelI
188           qLC = Qi + 0.5*( 1.+dm1h/(DelI+Eps) )*(qUL-Qi)  c        qLC = Qi + 0.5 _d 0*( 1. _d 0+dm1h/(DelI+Eps) )*(qUL-Qi)
189           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)  c        PhiMD = 2. _d 0/(1. _d 0-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
190           PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)  c        PhiLC = 2. _d 0*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
191    C--   simplify and avoid division by zero
192             PhiMD = 1. _d 0/(1. _d 0-cfl)*(DelIp-dp1h)*recip_DelIp
193             PhiLC = rp1h_cfl*( 1. _d 0+dm1h*recip_DelI )
194    C--
195           PhiMin = max(min(0. _d 0,PhiMD),           PhiMin = max(min(0. _d 0,PhiMD),
196       &        min(0. _d 0,2.*rp1h_cfl,PhiLC))       &        min(0. _d 0,2. _d 0*rp1h_cfl,PhiLC))
197           PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),           PhiMax = min(max(2. _d 0/(1. _d 0-cfl),PhiMD),
198       &        max(0. _d 0,2.*rp1h_cfl,PhiLC))       &        max(0. _d 0,2. _d 0*rp1h_cfl,PhiLC))
199           Phi = max(PhiMin,min(Phi,PhiMax))           Phi = max(PhiMin,min(Phi,PhiMax))
200    
201           Psi = Phi * 0.5 * (1. - cfl)           Psi = Phi * 0.5 _d 0 * (1. _d 0 - cfl)
202           wT(i,j) = wTrans(i,j)*( Qi + Psi*DelIp )           wT(i,j) = wTrans(i,j)*( Qi + Psi*DelIp )
203          ELSE          ELSE
204           wT(i,j) = 0.           wT(i,j) = 0. _d 0
205          ENDIF          ENDIF
206    
207         ENDDO         ENDDO
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22