/[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.8 by jmc, Mon Jun 19 14:40:43 2006 UTC revision 1.9 by jmc, Tue Dec 5 22:25:41 2006 UTC
# Line 53  C  km1               :: =max( k-1 , 1 ) Line 53  C  km1               :: =max( k-1 , 1 )
53  C  wLoc              :: velocity, vertical component  C  wLoc              :: velocity, vertical component
54  C  wCFL              :: Courant-Friedrich-Levy number  C  wCFL              :: Courant-Friedrich-Levy number
55        INTEGER i,j,kp1,km1,km2        INTEGER i,j,kp1,km1,km2
56        _RL Rjm,Rj,Rjp,cfl,d0,d1        _RL Rjm,Rj,Rjp,wCFL,d0,d1
57        _RL psiP,psiM,thetaP,thetaM        _RL psiP,psiM,thetaP,thetaM
58        _RL wLoc        _RL wLoc
59        _RL thetaMax        _RL thetaMax
# Line 73  C  wCFL              :: Courant-Friedric Line 73  C  wCFL              :: Courant-Friedric
73       &         *maskC(i,j,km1,bi,bj)       &         *maskC(i,j,km1,bi,bj)
74    
75          wLoc = wFld(i,j)          wLoc = wFld(i,j)
76  c       wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)          wCFL = ABS(wLoc*dTarg*recip_drC(k))
77          cfl=abs(wLoc*dTarg*recip_drC(k))          d0=(2. _d 0 -wCFL)*(1. _d 0 -wCFL)*oneSixth
78          d0=(2. _d 0 -cfl)*(1. _d 0 -cfl)*oneSixth          d1=(1. _d 0 -wCFL*wCFL)*oneSixth
         d1=(1. _d 0 -cfl*cfl)*oneSixth  
79    
80  C-      the old version: can produce overflow, division by zero,  C-      the old version: can produce overflow, division by zero,
81  C       and is wrong for tracer with low concentration:  C       and is wrong for tracer with low concentration:
# Line 101  C-      prevent |thetaP,M| to reach too Line 100  C-      prevent |thetaP,M| to reach too
100    
101          psiP=d0+d1*thetaP          psiP=d0+d1*thetaP
102          psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),          psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
103       &                       thetaP*(1. _d 0 -cfl)/(cfl+1. _d -20) ))       &                       thetaP*(1. _d 0 -wCFL)/(wCFL+1. _d -20) ))
104          psiM=d0+d1*thetaM          psiM=d0+d1*thetaM
105          psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),          psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
106       &                       thetaM*(1. _d 0 -cfl)/(cfl+1. _d -20) ))       &                       thetaM*(1. _d 0 -wCFL)/(wCFL+1. _d -20) ))
107    
108          wT(i,j)=          wT(i,j)=
109       &   0.5*(rTrans(i,j)+abs(rTrans(i,j)))       &   0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
110       &      *( tracer(i,j, k ) + psiM*Rj )       &      *( tracer(i,j, k ) + psiM*Rj )
111       &  +0.5*(rTrans(i,j)-abs(rTrans(i,j)))       &  +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
112       &      *( tracer(i,j,km1) - psiP*Rj )       &      *( tracer(i,j,km1) - psiP*Rj )
113    
114         ENDDO         ENDDO
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22