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

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

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

revision 1.6 by jmc, Thu Oct 20 17:03:08 2005 UTC revision 1.7 by jmc, Sat Oct 22 20:17:44 2005 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: GAD_U3C4_IMPL_R  C     !ROUTINE: GAD_U3C4_IMPL_R
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE GAD_U3C4_IMPL_R(        SUBROUTINE GAD_U3C4_IMPL_R(
10       I           bi,bj,k, iMin,iMax,jMin,jMax,       I           bi,bj,k, iMin,iMax,jMin,jMax,
11       I           advectionScheme, deltaTarg, rTrans,       I           advectionScheme, deltaTarg, rTrans,
12       O           a5d, b5d, c5d, d5d, e5d,       O           a5d, b5d, c5d, d5d, e5d,
13       I           myThid )       I           myThid )
# Line 15  C     !INTERFACE: Line 15  C     !INTERFACE:
15  C     !DESCRIPTION:  C     !DESCRIPTION:
16  C     Compute matrix element to solve vertical advection implicitly  C     Compute matrix element to solve vertical advection implicitly
17  C      using 3rd order upwind advection scheme,  C      using 3rd order upwind advection scheme,
18    C         or 3rd order Direct Space and Time advection scheme,
19  C         or 4th order Centered advection scheme.  C         or 4th order Centered advection scheme.
20  C     Method:  C     Method:
21  C      contribution of vertical transport at interface k is added  C      contribution of vertical transport at interface k is added
# Line 39  C     jMin,jMax       :: computation dom Line 40  C     jMin,jMax       :: computation dom
40  C     advectionScheme :: advection scheme to use  C     advectionScheme :: advection scheme to use
41  C     deltaTarg       :: time step  C     deltaTarg       :: time step
42  C     rTrans          :: vertical volume transport  C     rTrans          :: vertical volume transport
 C     tFld            :: tracer field  
43  C     a5d             :: 2nd  lower diag of pentadiagonal matrix  C     a5d             :: 2nd  lower diag of pentadiagonal matrix
44  C     b5d             :: 1rst lower diag of pentadiagonal matrix  C     b5d             :: 1rst lower diag of pentadiagonal matrix
45  C     c5d             :: main diag       of pentadiagonal matrix  C     c5d             :: main diag       of pentadiagonal matrix
# Line 64  C     kp1             :: =min( k+1 , Nr Line 64  C     kp1             :: =min( k+1 , Nr
64  C     km2             :: =max( k-2 , 1 )  C     km2             :: =max( k-2 , 1 )
65  C     rCenter         :: centered contribution  C     rCenter         :: centered contribution
66  C     rUpwind         :: upwind   contribution  C     rUpwind         :: upwind   contribution
67    C     rC4km, rC4kp    :: high order contribution
68    C     rHigh           :: high order term factor
69        LOGICAL flagC4        LOGICAL flagC4
70        INTEGER i,j,kp1,km2        INTEGER i,j,kp1,km2
71        _RL rCenter, rUpwind        _RL wCFL, rCenter, rUpwind
72        _RL rC4km, rC4kp, rU1k, rU3km, rU3kp        _RL rC4km, rC4kp, rHigh
73        _RL mskM, mskP, maskM2, maskP1        _RL mskM, mskP, maskM2, maskP1
74          _RL deltaTcfl
75  CEOP  CEOP
76    
77  C--   process interior interface only:  C--   process interior interface only:
# Line 76  C--   process interior interface only: Line 79  C--   process interior interface only:
79    
80        km2=MAX(1,k-2)        km2=MAX(1,k-2)
81        kp1=MIN(Nr,k+1)        kp1=MIN(Nr,k+1)
82        maskP1 = 1. _d 0        maskP1 = 1. _d 0
83        maskM2 = 1. _d 0        maskM2 = 1. _d 0
84        IF ( k.LE.2 ) maskM2 = 0. _d 0        IF ( k.LE.2 ) maskM2 = 0. _d 0
85        IF ( k.GE.Nr) maskP1 = 0. _d 0        IF ( k.GE.Nr) maskP1 = 0. _d 0
86        flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH        flagC4 = advectionScheme.EQ.ENUM_CENTERED_4TH
87       &         .AND. k.GT.2 .AND. k.LT.Nr       &         .AND. k.GT.2 .AND. k.LT.Nr
88    
89  C--    Add centered & upwind contributions  C--    Add centered, upwind and high-order contributions
90           deltaTcfl = deltaTarg(k)
91         DO j=jMin,jMax         DO j=jMin,jMax
92           DO i=iMin,iMax           DO i=iMin,iMax
93             rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign             rCenter= 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
94             mskM   = maskC(i,j,km2,bi,bj)*maskM2             mskM   = maskC(i,j,km2,bi,bj)*maskM2
95             mskP   = maskC(i,j,kp1,bi,bj)*maskP1             mskP   = maskC(i,j,kp1,bi,bj)*maskP1
            rC4km  = oneSixth*rCenter*mskM  
            rC4kp  = oneSixth*rCenter*mskP  
96             IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN             IF ( flagC4 .AND. mskM*mskP.GT.0. _d 0 ) THEN
97              rUpwind= 0. _d 0              rUpwind= 0. _d 0
98              rU3km  = 0. _d 0              rC4km  = oneSixth*rCenter*mskM
99              rU3kp  = 0. _d 0              rC4kp  = oneSixth*rCenter*mskP
100               ELSEIF ( advectionScheme.EQ.ENUM_DST3 ) THEN
101                wCFL = deltaTcfl*ABS(rTrans(i,j))
102         &            *recip_rA(i,j,bi,bj)*recip_drC(k)
103                rHigh  = (1. _d 0 -wCFL*wCFL)*oneSixth
104    c           rUpwind= (2. _d 0*rHigh - wCFL)*ABS(rCenter)
105                rUpwind= (2. _d 0*rHigh )*ABS(rCenter)
106                rC4km  =  rHigh * (rCenter+ABS(rCenter))*mskM
107                rC4kp  =  rHigh * (rCenter-ABS(rCenter))*mskP
108             ELSE             ELSE
109              rU1k   = oneSixth*ABS(rCenter)              rUpwind=  2. _d 0*oneSixth*ABS(rCenter)
110              rUpwind= rU1k+rU1k              rC4km  = oneSixth*(rCenter+ABS(rCenter))*mskM
111              rU3km  = rU1k*mskM              rC4kp  = oneSixth*(rCenter-ABS(rCenter))*mskP
             rU3kp  = rU1k*mskP  
112             ENDIF             ENDIF
113             a5d(i,j,k)   = a5d(i,j,k)             a5d(i,j,k)   = a5d(i,j,k)
114       &                  + (rC4km + rU3km)       &                  + rC4km
115       &                   *deltaTarg(k)       &                   *deltaTarg(k)
116       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
117             b5d(i,j,k)   = b5d(i,j,k)             b5d(i,j,k)   = b5d(i,j,k)
118       &                  - (rCenter + rC4km + rUpwind + rU3km)       &                  - ( (rCenter+rUpwind) + rC4km )
119       &                   *deltaTarg(k)       &                   *deltaTarg(k)
120       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
121             c5d(i,j,k)   = c5d(i,j,k)             c5d(i,j,k)   = c5d(i,j,k)
122       &                  - (rCenter + rC4kp - rUpwind - rU3kp)       &                  - ( (rCenter-rUpwind) + rC4kp )
123       &                   *deltaTarg(k)       &                   *deltaTarg(k)
124       &                    *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                    *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
125             d5d(i,j,k)   = d5d(i,j,k)             d5d(i,j,k)   = d5d(i,j,k)
126       &                  + (rC4kp - rU3kp)       &                  + rC4kp
127       &                   *deltaTarg(k)       &                   *deltaTarg(k)
128       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
129             b5d(i,j,k-1) = b5d(i,j,k-1)             b5d(i,j,k-1) = b5d(i,j,k-1)
130       &                  - (rC4km + rU3km)       &                  - rC4km
131       &                   *deltaTarg(k-1)       &                   *deltaTarg(k-1)
132       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
133             c5d(i,j,k-1) = c5d(i,j,k-1)             c5d(i,j,k-1) = c5d(i,j,k-1)
134       &                  + (rCenter + rC4km + rUpwind + rU3km)       &                  + ( (rCenter+rUpwind) + rC4km )
135       &                   *deltaTarg(k-1)       &                   *deltaTarg(k-1)
136       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
137             d5d(i,j,k-1) = d5d(i,j,k-1)             d5d(i,j,k-1) = d5d(i,j,k-1)
138       &                  + (rCenter + rC4kp - rUpwind - rU3kp)       &                  + ( (rCenter-rUpwind) + rC4kp )
139       &                   *deltaTarg(k-1)       &                   *deltaTarg(k-1)
140       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
141             e5d(i,j,k-1) = e5d(i,j,k-1)             e5d(i,j,k-1) = e5d(i,j,k-1)
142       &                  - (rC4kp - rU3kp)       &                  - rC4kp
143       &                   *deltaTarg(k-1)       &                   *deltaTarg(k-1)
144       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                   *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
145           ENDDO           ENDDO
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22