/[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.9 by jmc, Tue Dec 5 22:25:41 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_DST3FL_ADV_R(  CBOP
7       I           bi_arg,bj_arg,k,dTarg,  C !ROUTINE: GAD_DST3FL_ADV_R
8       I           rTrans, wVel,  
9    C !INTERFACE: ==========================================================
10          SUBROUTINE GAD_DST3FL_ADV_R(
11         I           bi,bj,k,dTarg,
12         I           rTrans, wFld,
13       I           tracer,       I           tracer,
14       O           wT,       O           wT,
15       I           myThid )       I           myThid )
16  C     /==========================================================\  
17  C     | SUBROUTINE GAD_DST3_ADV_R                                |  C !DESCRIPTION:
18  C     | o Compute Vertical advective Flux of Tracer using        |  C  Calculates the area integrated vertical flux due to advection of a tracer
19  C     |   3rd Order DST Sceheme with flux limiting               |  C  using 3rd Order DST Scheme with flux limiting
20  C     |==========================================================|  
21    C !USES: ===============================================================
22        IMPLICIT NONE        IMPLICIT NONE
23    
24  C     == GLobal variables ==  C     == GLobal variables ==
25  #include "SIZE.h"  #include "SIZE.h"
26  #include "GRID.h"  #include "GRID.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
27  #include "GAD.h"  #include "GAD.h"
28    
29  C     == Routine arguments ==  C     == Routine arguments ==
30        INTEGER bi_arg,bj_arg,k  C !INPUT PARAMETERS: ===================================================
31    C  bi,bj             :: tile indices
32    C  k                 :: vertical level
33    C  deltaTloc         :: local time-step (s)
34    C  rTrans            :: vertical volume transport
35    C  wFld              :: vertical flow
36    C  tracer            :: tracer field
37    C  myThid            :: thread number
38          INTEGER bi,bj,k
39        _RL dTarg        _RL dTarg
40        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41        _RL wVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
43        INTEGER myThid        INTEGER myThid
44    
45    C !OUTPUT PARAMETERS: ==================================================
46    C  wT                :: vertical advective flux
47          _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48    
49  C     == Local variables ==  C     == Local variables ==
50        INTEGER i,j,kp1,km1,km2,bi,bj  C !LOCAL VARIABLES: ====================================================
51        _RL Rjm,Rj,Rjp,cfl,d0,d1  C  i,j               :: loop indices
52    C  km1               :: =max( k-1 , 1 )
53    C  wLoc              :: velocity, vertical component
54    C  wCFL              :: Courant-Friedrich-Levy number
55          INTEGER i,j,kp1,km1,km2
56          _RL Rjm,Rj,Rjp,wCFL,d0,d1
57        _RL psiP,psiM,thetaP,thetaM        _RL psiP,psiM,thetaP,thetaM
58          _RL wLoc
59        IF (.NOT. multiDimAdvection) THEN        _RL thetaMax
60  C      If using the standard time-stepping/advection schemes (ie. AB-II)        PARAMETER( thetaMax = 1.D+20 )
 C      then the data-structures are all global arrays  
        bi=bi_arg  
        bj=bj_arg  
       ELSE  
 C      otherwise if using the multi-dimensional advection schemes  
 C      then the data-structures are all local arrays except  
 C      for maskC(...) and wVel(...)  
        bi=1  
        bj=1  
       ENDIF  
61    
62        km2=MAX(1,k-2)        km2=MAX(1,k-2)
63        km1=MAX(1,k-1)        km1=MAX(1,k-1)
# Line 56  C      for maskC(...) and wVel(...) Line 65  C      for maskC(...) and wVel(...)
65    
66        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
67         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
68          Rjp=(tracer(i,j,k,bi,bj)-tracer(i,j,kp1,bi,bj))          Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
69       &         *maskC(i,j,kp1,bi_arg,bj_arg)       &         *maskC(i,j,kp1,bi,bj)
70          Rj =(tracer(i,j,km1,bi,bj)-tracer(i,j,k,bi,bj))          Rj =(tracer(i,j,km1)-tracer(i,j,k))
71       &         *maskC(i,j,k,bi_arg,bj_arg)*maskC(i,j,km1,bi_arg,bj_arg)       &         *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
72          Rjm=(tracer(i,j,km2,bi,bj)-tracer(i,j,km1,bi,bj))          Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
73       &         *maskC(i,j,km1,bi_arg,bj_arg)       &         *maskC(i,j,km1,bi,bj)
74    
75          cfl=abs(wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drc(k))          wLoc = wFld(i,j)
76          d0=(2.D0-cfl)*(1.-cfl)*oneSixth          wCFL = ABS(wLoc*dTarg*recip_drC(k))
77          d1=(1.D0-cfl*cfl)*oneSixth          d0=(2. _d 0 -wCFL)*(1. _d 0 -wCFL)*oneSixth
78            d1=(1. _d 0 -wCFL*wCFL)*oneSixth
79    
80    C-      the old version: can produce overflow, division by zero,
81    C       and is wrong for tracer with low concentration:
82    c       thetaP=Rjm/(1.D-20+Rj)
83    c       thetaM=Rjp/(1.D-20+Rj)
84    C-      the right expression, but not bounded:
85  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)  
86  c       thetaM=0.D0  c       thetaM=0.D0
87    c       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
88  c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj  c       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
89    C-      prevent |thetaP,M| to reach too big value:
90            IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
91              thetaP=SIGN(thetaMax,Rjm*Rj)
92            ELSE
93              thetaP=Rjm/Rj
94            ENDIF
95            IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
96              thetaM=SIGN(thetaMax,Rjp*Rj)
97            ELSE
98              thetaM=Rjp/Rj
99            ENDIF
100    
101            psiP=d0+d1*thetaP
102            psiP=MAX(0. _d 0,MIN(MIN(1. _d 0,psiP),
103         &                       thetaP*(1. _d 0 -wCFL)/(wCFL+1. _d -20) ))
104          psiM=d0+d1*thetaM          psiM=d0+d1*thetaM
105          psiM=max(0.D0,min(min(1.D0,psiM),          psiM=MAX(0. _d 0,MIN(MIN(1. _d 0,psiM),
106       &       (1.D0-cfl)/(1.D-20+cfl)*thetaM))       &                       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 ,bi,bj) + 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,bi,bj) - psiP*Rj )       &      *( tracer(i,j,km1) - psiP*Rj )
113    
114         ENDDO         ENDDO
115        ENDDO        ENDDO
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22