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

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

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

revision 1.4 by jmc, Wed Mar 6 01:29:36 2002 UTC revision 1.5 by jmc, Tue Oct 18 16:03:55 2005 UTC
# Line 38  C     wFld     :: velocity, vertical com Line 38  C     wFld     :: velocity, vertical com
38        _RL Rjm,Rj,Rjp,cfl,d0,d1        _RL Rjm,Rj,Rjp,cfl,d0,d1
39        _RL psiP,psiM,thetaP,thetaM        _RL psiP,psiM,thetaP,thetaM
40        _RL wFld        _RL wFld
41          _RL thetaMax
42          PARAMETER( thetaMax = 1.D+20 )
43    
44        IF (.NOT. multiDimAdvection) THEN        IF (.NOT. multiDimAdvection) THEN
45  C      If using the standard time-stepping/advection schemes (ie. AB-II)  C      If using the standard time-stepping/advection schemes (ie. AB-II)
# Line 68  C      for maskC(...) and wVel(...) Line 70  C      for maskC(...) and wVel(...)
70  c       wFld = wVel(i,j,k,bi_arg,bj_arg)  c       wFld = wVel(i,j,k,bi_arg,bj_arg)
71          wFld = rTrans(i,j)*recip_rA(i,j,bi_arg,bj_arg)          wFld = rTrans(i,j)*recip_rA(i,j,bi_arg,bj_arg)
72          cfl=abs(wFld*dTarg*recip_drC(k))          cfl=abs(wFld*dTarg*recip_drC(k))
73          d0=(2.D0-cfl)*(1.-cfl)*oneSixth          d0=(2. _d 0 -cfl)*(1. _d 0 -cfl)*oneSixth
74          d1=(1.D0-cfl*cfl)*oneSixth          d1=(1. _d 0 -cfl*cfl)*oneSixth
75    
76    C-      the old version: can produce overflow, division by zero,
77    C       and is wrong for tracer with low concentration:
78    c       thetaP=Rjm/(1.D-20+Rj)
79    c       thetaM=Rjp/(1.D-20+Rj)
80    C-      the right expression, but not bounded:
81  c       thetaP=0.D0  c       thetaP=0.D0
 c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj  
         thetaP=Rjm/(1.D-20+Rj)  
         psiP=d0+d1*thetaP  
         psiP=max(0.D0,min(min(1.D0,psiP),  
      &       (1.D0-cfl)/(1.D-20+cfl)*thetaP))  
         thetaM=Rjp/(1.D-20+Rj)  
82  c       thetaM=0.D0  c       thetaM=0.D0
83    c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
84  c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj  c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
85    C-      prevent |thetaP,M| to reach too big value:
86            IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
87              thetaP=SIGN(thetaMax,Rjm*Rj)
88            ELSE
89              thetaP=Rjm/Rj
90            ENDIF
91            IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
92              thetaM=SIGN(thetaMax,Rjp*Rj)
93            ELSE
94              thetaM=Rjp/Rj
95            ENDIF
96    
97            psiP=d0+d1*thetaP
98            psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
99         &                       thetaP*(1. _d 0 -cfl)/(cfl+1. _d -20) ))
100          psiM=d0+d1*thetaM          psiM=d0+d1*thetaM
101          psiM=max(0.D0,min(min(1.D0,psiM),          psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
102       &       (1.D0-cfl)/(1.D-20+cfl)*thetaM))       &                       thetaM*(1. _d 0 -cfl)/(cfl+1. _d -20) ))
103    
104          wT(i,j)=          wT(i,j)=
105       &   0.5*(rTrans(i,j)+abs(rTrans(i,j)))       &   0.5*(rTrans(i,j)+abs(rTrans(i,j)))
106       &      *( Tracer(i,j, k ,bi,bj) + psiM*Rj )       &      *( Tracer(i,j, k ,bi,bj) + psiM*Rj )
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22