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