/[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.1 by adcroft, Sat Jan 20 21:20:11 2007 UTC revision 1.4 by adcroft, Fri May 11 17:30:06 2007 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,Del,DelIp,DelI,Phi,Eps,rp1h,rp1h_cfl,Msk
39        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm        _RL Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm
40        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm        _RL MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
41        _RL d2,d2p1,d2m1,A,B,C,D        _RL d2,d2p1,d2m1,A,B,C,D
# Line 53  C     == Local variables == Line 54  C     == Local variables ==
54         DO i=1-Olx+4,sNx+Olx-3         DO i=1-Olx+4,sNx+Olx-3
55    
56          uLoc = uFld(i,j)          uLoc = uFld(i,j)
57          cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))          cfl = uLoc
58            IF ( calcCFL ) cfl = abs(uLoc*deltaTloc*recip_dxC(i,j,bi,bj))
59    
60          IF (uLoc.gt.0.) THEN          IF (uTrans(i,j).gt.0.) THEN
61           Qippp = Q(i+2,j)           Qippp = Q(i+2,j)
62           Qipp  = Q(i+1,j)           Qipp  = Q(i+1,j)
63           Qip   = Q(i,j)           Qip   = Q(i,j)
# Line 70  C     == Local variables == Line 72  C     == Local variables ==
72           MskIm   = maskLocW(i-1,j)           MskIm   = maskLocW(i-1,j)
73           MskImm  = maskLocW(i-2,j)           MskImm  = maskLocW(i-2,j)
74           MskImmm = maskLocW(i-3,j)           MskImmm = maskLocW(i-3,j)
75          ELSEIF (uLoc.lt.0.) THEN          ELSEIF (uTrans(i,j).lt.0.) THEN
76           Qippp = Q(i-3,j)           Qippp = Q(i-3,j)
77           Qipp  = Q(i-2,j)           Qipp  = Q(i-2,j)
78           Qip   = Q(i-1,j)           Qip   = Q(i-1,j)
# Line 87  C     == Local variables == Line 89  C     == Local variables ==
89           MskImmm = maskLocW(i+3,j)           MskImmm = maskLocW(i+3,j)
90          ENDIF          ENDIF
91    
92          IF (uLoc.ne.0.) THEN          IF (uTrans(i,j).ne.0.) THEN
93  C        2nd order correction [i i-1]  C        2nd order correction [i i-1]
94           Fac = 1.           Fac = 1.
95           Del = Qip-Qi           Del = Qip-Qi
# Line 120  C        7th order correction [i+2 i+1 i Line 122  C        7th order correction [i+2 i+1 i
122           Phi = Phi - Msk * Fac * Del           Phi = Phi - Msk * Fac * Del
123    
124           DelIp = ( Qip - Qi ) * MskI           DelIp = ( Qip - Qi ) * MskI
125           Phi = sign(1.,Phi)*sign(1.,DelIp)*abs(Phi+Eps)/abs(DelIp+Eps)           Phi = sign(1. _d 0,Phi)*sign(1. _d 0,DelIp)
126         &        *abs(Phi+Eps)/abs(DelIp+Eps)
127    
128           DelI = ( Qi - Qim ) * MskIm           DelI = ( Qi - Qim ) * MskIm
129           rp1h =sign(1.,DelI)*sign(1.,DelIp)*abs(DelI+Eps)/abs(DelIp+Eps)           rp1h =sign(1. _d 0,DelI)*sign(1. _d 0,DelIp)
130         &        *abs(DelI+Eps)/abs(DelIp+Eps)
131             rp1h_cfl = rp1h/(cfl+Eps)
132    
133  C        TVD limiter  C        TVD limiter
134  !        Phi = max(0., min( 2./(1-cfl), Phi, 2.*rp1h/cfl ) )  !        Phi = max(0. _d 0, min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )
135    
136  C        MP limiter  C        MP limiter
137           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm           d2   = ( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
# Line 136  C        MP limiter Line 141  C        MP limiter
141           B = 4.*d2p1 - d2           B = 4.*d2p1 - d2
142           C = d2           C = d2
143           D = d2p1;           D = d2p1;
144           dp1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)           dp1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
145           A = 4.*d2m1 - d2           A = 4.*d2m1 - d2
146           B = 4.*d2 - d2m1           B = 4.*d2 - d2m1
147           C = d2m1           C = d2m1
148           D = d2;           D = d2;
149           dm1h = max(min(A,B,C,D),0.)+min(max(A,B,C,D),0.)           dm1h = max(min(A,B,C,D),0. _d 0)+min(max(A,B,C,D),0. _d 0)
150           qMD = 0.5*( ( Qi + Qip ) - dp1h )           qMD = 0.5*( ( Qi + Qip ) - dp1h )
151           qUL = Qi + (1.-cfl)/cfl*( Qi-Qim )           qUL = Qi + (1.-cfl)/(cfl+Eps)*( Qi-Qim )
152           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)           qLC = Qi + 0.5*( 1.+dm1h/(Qi-Qim+Eps) )*(qUL-Qi)
153           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)           PhiMD = 2./(1.-cfl)*(qMD-Qi+Eps)/(Qip-Qi+Eps)
154           PhiLC = 2.*rp1h/cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)           PhiLC = 2.*rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
155           PhiMin = max(min(0.,PhiMD),min(0.,2.*rp1h/cfl,PhiLC))           PhiMin = max(min(0. _d 0,PhiMD),
156           PhiMax = min(max(2./(1.-cfl),PhiMD),max(0.,2.*rp1h/cfl,PhiLC))       &        min(0. _d 0,2.*rp1h_cfl,PhiLC))
157             PhiMax = min(max(2. _d 0/(1.-cfl),PhiMD),
158         &        max(0. _d 0,2.*rp1h_cfl,PhiLC))
159           Phi = max(PhiMin,min(Phi,PhiMax))           Phi = max(PhiMin,min(Phi,PhiMax))
160    
161           Psi = Phi * 0.5 * (1. - cfl)           Psi = Phi * 0.5 * (1. - cfl)
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22