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

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

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

revision 1.1.4.1 by heimbach, Mon Apr 8 20:10:39 2002 UTC revision 1.11 by jmc, Wed Jul 4 20:22:26 2012 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_DST3_ADV_R(  CBOP
7       I           bi_arg,bj_arg,k,dTarg,  C !ROUTINE: GAD_DST3_ADV_R
8       I           rTrans, wVel,  
9    C !INTERFACE: ==========================================================
10          SUBROUTINE GAD_DST3_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                                  |  C  using 3rd-order Direct Space and Time (DST-3) Advection Scheme
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"  #ifdef OLD_DST3_FORMULATION
27  #include "EEPARAMS.h"  #include "EEPARAMS.h"
28  #include "PARAMS.h"  #include "PARAMS.h"
29    #endif
30    #include "GRID.h"
31  #include "GAD.h"  #include "GAD.h"
32    
33  C     == Routine arguments ==  C     == Routine arguments ==
34        INTEGER bi_arg,bj_arg,k  C !INPUT PARAMETERS: ===================================================
35    C  bi,bj             :: tile indices
36    C  k                 :: vertical level
37    C  deltaTloc         :: local time-step (s)
38    C  rTrans            :: vertical volume transport
39    C  wFld              :: vertical flow
40    C  tracer            :: tracer field
41    C  myThid            :: thread number
42          INTEGER bi,bj,k
43        _RL dTarg        _RL dTarg
44        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45        _RL wVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46        _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)  
47        INTEGER myThid        INTEGER myThid
48    
49    C !OUTPUT PARAMETERS: ==================================================
50    C  wT                :: vertical advective flux
51          _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52    
53  C     == Local variables ==  C     == Local variables ==
54        INTEGER i,j,kp1,km1,km2,bi,bj  C !LOCAL VARIABLES: ====================================================
55    C  i,j               :: loop indices
56    C  km1               :: =max( k-1 , 1 )
57    C  wLoc              :: velocity, vertical component
58    C  wCFL              :: Courant-Friedrich-Levy number
59          INTEGER i,j,kp1,km1,km2
60          _RL wLoc
61        _RL Rjm,Rj,Rjp,cfl,d0,d1        _RL Rjm,Rj,Rjp,cfl,d0,d1
62    #ifdef OLD_DST3_FORMULATION
63        _RL psiP,psiM,thetaP,thetaM        _RL psiP,psiM,thetaP,thetaM
64          _RL smallNo
65    
66        IF (.NOT. multiDimAdvection) THEN  c     IF (inAdMode) THEN
67  C      If using the standard time-stepping/advection schemes (ie. AB-II)  c      smallNo = 1.0D-20
68  C      then the data-structures are all global arrays  c     ELSE
69         bi=bi_arg         smallNo = 1.0D-20
70         bj=bj_arg  c     ENDIF
71        ELSE  #endif
 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  
72    
73        km2=MAX(1,k-2)        km2=MAX(1,k-2)
74        km1=MAX(1,k-1)        km1=MAX(1,k-1)
75        kp1=MIN(Nr,k+1)        kp1=MIN(Nr,k+1)
76    
77        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
78         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
79          Rjp=(tracer(i,j,k,bi,bj)-tracer(i,j,kp1,bi,bj))          Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
80       &         *maskC(i,j,kp1,bi_arg,bj_arg)       &         *maskC(i,j,kp1,bi,bj)
81          Rj =(tracer(i,j,km1,bi,bj)-tracer(i,j,k,bi,bj))          Rj =(tracer(i,j,km1)-tracer(i,j,k))
82       &         *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)
83          Rjm=(tracer(i,j,km2,bi,bj)-tracer(i,j,km1,bi,bj))          Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
84       &         *maskC(i,j,km1,bi_arg,bj_arg)       &         *maskC(i,j,km1,bi,bj)
85    
86          cfl=abs(wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drc(k))          wLoc = wFld(i,j)
87    c       wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
88            cfl=ABS(wLoc*dTarg*recip_drC(k))
89          d0=(2.-cfl)*(1.-cfl)*oneSixth          d0=(2.-cfl)*(1.-cfl)*oneSixth
90          d1=(1.-cfl*cfl)*oneSixth          d1=(1.-cfl*cfl)*oneSixth
91  c       thetaP=0.  #ifdef OLD_DST3_FORMULATION
92  c       IF (Rj.NE.0.) thetaP=Rjm/Rj          IF ( ABS(Rj).LT.smallNo .OR.
93          thetaP=Rjm/(1.D-20+Rj)       &       ABS(Rjm).LT.smallNo ) THEN
94          psiP=d0+d1*thetaP           thetaP=0.
95  c       psiP=max(0.,min(min(1.,psiP),(1.-cfl)/(1.D-20+cfl)*thetaP))           psiP=0.
96          thetaM=Rjp/(1.D-20+Rj)          ELSE
97  c       thetaM=0.           thetaP=(Rjm+smallNo)/(smallNo+Rj)
98  c       IF (Rj.NE.0.) thetaM=Rjp/Rj           psiP=d0+d1*thetaP
99          psiM=d0+d1*thetaM          ENDIF
100  c       psiM=max(0.,min(min(1.,psiM),(1.-cfl)/(1.D-20+cfl)*thetaM))          IF ( ABS(Rj).LT.smallNo .OR.
101         &       ABS(Rjp).LT.smallNo ) THEN
102             thetaM=0.
103             psiM=0.
104            ELSE
105             thetaM=(Rjp+smallNo)/(smallNo+Rj)
106             psiM=d0+d1*thetaM
107            ENDIF
108             wT(i,j)=
109         &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
110         &       *( tracer(i,j, k ) + psiM*Rj )
111         &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
112         &       *( tracer(i,j,km1) - psiP*Rj )
113    #else /* OLD_DST3_FORMULATION */
114          wT(i,j)=          wT(i,j)=
115       &   0.5*(rTrans(i,j)+abs(rTrans(i,j)))       &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
116       &      *( Tracer(i,j, k ,bi,bj) + psiM*Rj )       &       *( tracer(i,j, k ) + (d0*Rj+d1*Rjp) )
117       &  +0.5*(rTrans(i,j)-abs(rTrans(i,j)))       &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
118       &      *( Tracer(i,j,km1,bi,bj) - psiP*Rj )       &       *( tracer(i,j,km1) - (d0*Rj+d1*Rjm) )
119    #endif /* OLD_DST3_FORMULATION */
120    
121         ENDDO         ENDDO
122        ENDDO        ENDDO
Legend:
Removed from v.1.1.4.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22