/[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.3 by adcroft, Thu Jan 17 16:48:59 2002 UTC revision 1.5 by jmc, Tue Oct 18 16:03:55 2005 UTC
# Line 33  C     == Routine arguments == Line 33  C     == Routine arguments ==
33        INTEGER myThid        INTEGER myThid
34    
35  C     == Local variables ==  C     == Local variables ==
36    C     wFld     :: velocity, vertical component
37        INTEGER i,j,kp1,km1,km2,bi,bj        INTEGER i,j,kp1,km1,km2,bi,bj
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
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 63  C      for maskC(...) and wVel(...) Line 67  C      for maskC(...) and wVel(...)
67          Rjm=(tracer(i,j,km2,bi,bj)-tracer(i,j,km1,bi,bj))          Rjm=(tracer(i,j,km2,bi,bj)-tracer(i,j,km1,bi,bj))
68       &         *maskC(i,j,km1,bi_arg,bj_arg)       &         *maskC(i,j,km1,bi_arg,bj_arg)
69    
70          cfl=abs(wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drc(k))  c       wFld = wVel(i,j,k,bi_arg,bj_arg)
71          d0=(2.D0-cfl)*(1.-cfl)*oneSixth          wFld = rTrans(i,j)*recip_rA(i,j,bi_arg,bj_arg)
72          d1=(1.D0-cfl*cfl)*oneSixth          cfl=abs(wFld*dTarg*recip_drC(k))
73            d0=(2. _d 0 -cfl)*(1. _d 0 -cfl)*oneSixth
74            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.3  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22