/[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.2 by mlosch, Sun Jan 21 17:25:31 2007 UTC revision 1.3 by adcroft, Fri May 11 17:30:06 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,Msk        _RL wLoc,Fac,Del,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl,Msk
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
# Line 55  C     == Local variables == Line 55  C     == Local variables ==
55          wLoc = wFld(i,j)          wLoc = wFld(i,j)
56          cfl = abs(wLoc*deltaTloc*recip_drC(k))          cfl = abs(wLoc*deltaTloc*recip_drC(k))
57    
58          IF (wLoc.lt.0.) THEN          IF (wTrans(i,j).lt.0.) THEN
59           Qippp = Q(i,j,kp2)           Qippp = Q(i,j,kp2)
60           Qipp  = Q(i,j,kp1)           Qipp  = Q(i,j,kp1)
61           Qip   = Q(i,j,k)           Qip   = Q(i,j,k)
# Line 70  C     == Local variables == Line 70  C     == Local variables ==
70           MskIm   = maskC(i,j,km1,bi,bj) * float(km1-km2)           MskIm   = maskC(i,j,km1,bi,bj) * float(km1-km2)
71           MskImm  = maskC(i,j,km2,bi,bj) * float(km2-km3)           MskImm  = maskC(i,j,km2,bi,bj) * float(km2-km3)
72           MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)           MskImmm = maskC(i,j,km3,bi,bj) * float(km3-km4)
73          ELSEIF (wLoc.gt.0.) THEN          ELSEIF (wTrans(i,j).gt.0.) THEN
74           Qippp = Q(i,j,km3)           Qippp = Q(i,j,km3)
75           Qipp  = Q(i,j,km2)           Qipp  = Q(i,j,km2)
76           Qip   = Q(i,j,km1)           Qip   = Q(i,j,km1)
# Line 87  C     == Local variables == Line 87  C     == Local variables ==
87           MskImmm = maskC(i,j,kp3,bi,bj) * float(kp3-kp2)           MskImmm = maskC(i,j,kp3,bi,bj) * float(kp3-kp2)
88          ENDIF          ENDIF
89    
90          IF (wLoc.ne.0.) THEN          IF (wTrans(i,j).ne.0.) THEN
91           Phi = 0.           Phi = 0.
92  C        2nd order correction [i i-1]  C        2nd order correction [i i-1]
93           Fac = 1.           Fac = 1.
# Line 127  C        7th order correction [i+2 i+1 i Line 127  C        7th order correction [i+2 i+1 i
127           DelI = ( Qi - Qim ) * MskIm           DelI = ( Qi - Qim ) * MskIm
128           rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)           rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
129       &        *abs(DelI+Eps)/abs(DelIp+Eps)       &        *abs(DelI+Eps)/abs(DelIp+Eps)
130             rp1h_cfl = rp1h/(cfl+Eps)
131    
132  C        TVD limiter  C        TVD limiter
133  !        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 ) )
134    
135  C        MP limiter  C        MP limiter
136           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
# Line 146  C        MP limiter Line 147  C        MP limiter
147           D = d2;           D = d2;
148           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)
149           qMD = 0.5*( ( Qi + Qip ) - dp1h )           qMD = 0.5*( ( Qi + Qip ) - dp1h )
150           qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )           qUL = Qi + (1.-cfl)/(cfl+Eps)*( Qi-Qim )
151           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)
152           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
153           PhiLC = 2.*rp1h/cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)           PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
154           PhiMin = max(min(0. _d 0,PhiMD),           PhiMin = max(min(0. _d 0,PhiMD),
155       &        min(0. _d 0,2.*rp1h/cfl,PhiLC))       &        min(0. _d 0,2.*rp1h_cfl,PhiLC))
156           PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),           PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),
157       &        max(0. _d 0,2.*rp1h/cfl,PhiLC))       &        max(0. _d 0,2.*rp1h_cfl,PhiLC))
158           Phi = max(PhiMin,min(Phi,PhiMax))           Phi = max(PhiMin,min(Phi,PhiMax))
159    
160           Psi = Phi * 0.5 * (1. - cfl)           Psi = Phi * 0.5 * (1. - cfl)
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22