/[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.3 by adcroft, Fri May 11 17:30:06 2007 UTC revision 1.4 by adcroft, Fri May 11 18:24:31 2007 UTC
# Line 33  C     == Routine arguments == Line 33  C     == Routine arguments ==
33  C     == Local variables ==  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,Del,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl,Msk        _RL wLoc,Fac,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl
37        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
38        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
39        _RL d2,d2p1,d2m1,A,B,C,D        _RL d2,d2p1,d2m1,A,B,C,D
40        _RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax        _RL dp1h,dm1h,qMD,qUL,qLC,PhiMD,PhiLC,PhiMin,PhiMax
41          _RL DelM,DelP,DelMM,DelPP,DelMMM,DelPPP
42          _RL Del2MM,Del2M,Del2,Del2P,Del2PP
43          _RL Del3MM,Del3M,Del3P,Del3PP
44          _RL Del4M,Del4,Del4P
45          _RL Del5M,Del5P
46          _RL Del6
47    
48        Eps = 1. _d -20        Eps = 1. _d -20
49    
# Line 88  C     == Local variables == Line 94  C     == Local variables ==
94          ENDIF          ENDIF
95    
96          IF (wTrans(i,j).ne.0.) THEN          IF (wTrans(i,j).ne.0.) THEN
          Phi = 0.  
97  C        2nd order correction [i i-1]  C        2nd order correction [i i-1]
98           Fac = 1.           Fac = 1.
99           Del = Qip-Qi           DelP = (Qip-Qi)*MskI
100           Msk = MskI           Phi = Fac * DelP
          Phi = Msk * Fac * Del  
101  C        3rd order correction [i i-1 i-2]  C        3rd order correction [i i-1 i-2]
102           Fac = Fac * ( cfl + 1. )/3.           Fac = Fac * ( cfl + 1. )/3.
103           Del = Del - ( Qi-Qim )           DelM = (Qi-Qim)*MskIm
104           Msk = Msk * MskIm           Del2 = DelP - DelM
105           Phi = Phi - Msk * Fac * Del           Phi = Phi - Fac * Del2
106  C        4th order correction [i+1 i i-1 i-2]  C        4th order correction [i+1 i i-1 i-2]
107           Fac = Fac * ( cfl - 2. )/4.           Fac = Fac * ( cfl - 2. )/4.
108           Del = ( Qipp-2.*Qip+Qi ) - Del           DelPP = (Qipp-Qip)*MskIp*MskI
109           Msk = Msk * MskIp           Del2P = DelPP - DelP
110           Phi = Phi + Msk * Fac * Del           Del3P = Del2P - Del2
111             Phi = Phi + Fac * Del3p
112  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]
113           Fac = Fac * ( cfl - 3. )/5.           Fac = Fac * ( cfl - 3. )/5.
114           Del = Del - ( Qip-3.*Qi+3.*Qim-Qimm )           DelMM = (Qim-Qimm)*MskImm*MskIm
115           Msk = Msk * MskImm           Del2M = DelM - DelMM
116           Phi = Phi + Msk * Fac * Del           Del3M = Del2 - Del2M
117             Del4 = Del3P - Del3M
118             Phi = Phi + Fac * Del4
119  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]
120           Fac = Fac * ( cfl + 2. )/6.           Fac = Fac * ( cfl + 2. )/6.
121           Del = ( Qippp-4.*Qipp+6.*Qip-4.*Qi+Qim ) - Del           DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
122           Msk = Msk * MskIpp           Del2PP = DelPP - DelP
123           Phi = Phi + Msk * Fac * Del           Del3PP = Del2PP - Del2P
124             Del4P = Del3PP - Del3P
125             Del5P = Del4P - Del4
126             Phi = Phi + Fac * Del5P
127  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]
128           Fac = Fac * ( cfl + 2. )/7.           Fac = Fac * ( cfl + 2. )/7.
129           Del = Del - ( Qipp-5.*Qip+10.*Qi-10.*Qim+5.*Qimm-Qimmm )           DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
130           Msk = Msk * MskImmm           Del2MM = DelMM - DelMMM
131           Phi = Phi - Msk * Fac * Del           Del3MM = Del2M - Del2MM
132             Del4M = Del3M - Del3MM
133             Del5M = Del4 - Del4M
134             Del6 = Del5P - Del5M
135             Phi = Phi - Fac * Del6
136    
137           DelIp = ( Qip - Qi ) * MskI           DelIp = ( Qip - Qi ) * MskI
138           Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)           Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
# Line 133  C        TVD limiter Line 147  C        TVD limiter
147  !        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )  !        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )
148    
149  C        MP limiter  C        MP limiter
150           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm           d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
151           d2p1 = ( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI           d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
152           d2m1 = ( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm           d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
153           A = 4.*d2 - d2p1           A = 4.*d2 - d2p1
154           B = 4.*d2p1 - d2           B = 4.*d2p1 - d2
155           C = d2           C = d2
# Line 146  C        MP limiter Line 160  C        MP limiter
160           C = d2m1           C = d2m1
161           D = d2;           D = d2;
162           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)
163           qMD = 0.5*( ( Qi + Qip ) - dp1h )          !qMD = 0.5*( ( Qi + Qip ) - dp1h )
164           qUL = Qi + (1.-cfl)/(cfl+Eps)*( Qi-Qim )           qMD = 0.5*( ( 2.*Qi + DelIp ) - dp1h )
165           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)           qUL = Qi + (1.-cfl)/(cfl+Eps)*DelI
166           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)           qLC = Qi + 0.5*( 1.+dm1h/(DelI+Eps) )*(qUL-Qi)
167             PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
168           PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)           PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
169           PhiMin = max(min(0. _d 0,PhiMD),           PhiMin = max(min(0. _d 0,PhiMD),
170       &        min(0. _d 0,2.*rp1h_cfl,PhiLC))       &        min(0. _d 0,2.*rp1h_cfl,PhiLC))
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22