/[MITgcm]/MITgcm/model/src/calc_gw.F
ViewVC logotype

Diff of /MITgcm/model/src/calc_gw.F

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

revision 1.10 by edhill, Thu Oct 9 04:19:18 2003 UTC revision 1.32 by jmc, Thu Jul 13 03:01:21 2006 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
3  C $Name$  C $Name$
4    
 #include "PACKAGES_CONFIG.h"  
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_GW(            SUBROUTINE CALC_GW(
11       I        myThid)       I               bi, bj, KappaRU, KappaRV,
12         I               myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | S/R CALC_GW                                                C     | S/R CALC_GW
16  C     | o Calculate vert. velocity tendency terms ( NH, QH only )  C     | o Calculate vertical velocity tendency terms
17    C     |   ( Non-Hydrostatic only )
18  C     *==========================================================*  C     *==========================================================*
19  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH, the vertical momentum tendency must be
20  C     | calculated explicitly and included as a source term  C     | calculated explicitly and included as a source term
21  C     | for a 3d pressure eqn. Calculate that term here.  C     | for a 3d pressure eqn. Calculate that term here.
22  C     | This routine is not used in HYD calculations.  C     | This routine is not used in HYD calculations.
23  C     *==========================================================*  C     *==========================================================*
24  C     \ev  C     \ev
25    
26  C     !USES:  C     !USES:
27        IMPLICIT NONE        IMPLICIT NONE
28  C     == Global variables ==  C     == Global variables ==
29  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
33  #include "GW.h"  #include "DYNVARS.h"
34  #include "CG3D.h"  #include "NH_VARS.h"
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38  C     myThid - Instance number for this innvocation of CALC_GW  C     bi,bj   :: current tile indices
39    C     KappaRU :: vertical viscosity at U points
40    C     KappaRV :: vertical viscosity at V points
41    C     myTime  :: Current time in simulation
42    C     myIter  :: Current iteration number in simulation
43    C     myThid  :: Thread number for this instance of the routine.
44          INTEGER bi,bj
45          _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46          _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
47          _RL     myTime
48          INTEGER myIter
49        INTEGER myThid        INTEGER myThid
50    
51  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
52    
53  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
54  C     == Local variables ==  C     == Local variables ==
55  C     bi, bj,      :: Loop counters  C     iMin,iMax
56  C     iMin, iMax,  C     jMin,jMax
57  C     jMin, jMax  C     xA            :: W-Cell face area normal to X
58  C     flx_NS       :: Temp. used for fVol meridional terms.  C     yA            :: W-Cell face area normal to Y
59  C     flx_EW       :: Temp. used for fVol zonal terms.  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
60  C     flx_Up       :: Temp. used for fVol vertical terms.  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
61  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
62        INTEGER bi,bj,iMin,iMax,jMin,jMax  C     flx_NS        :: vertical momentum flux, meridional direction
63        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_EW        :: vertical momentum flux, zonal direction
64        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
65        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
66        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
67  C     I,J,K - Loop counters  C     gwDiss        :: vertical momentum dissipation tendency
68        INTEGER i,j,k, kP1, kUp  C     i,j,k         :: Loop counters
69          INTEGER iMin,iMax,jMin,jMax
70          _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71          _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72          _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73          _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74          _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75          _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76          _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          INTEGER i,j,k, kp1
84        _RL  wOverride        _RL  wOverride
85        _RS  hFacROpen        _RL  tmp_WbarZ
86        _RS  hFacRClosed        _RL  uTrans, vTrans, rTrans
87        _RL ab15,ab05        _RL  viscLoc
88        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL  halfRL
89          _RS  halfRS, zeroRS
90        _RL  Half        PARAMETER( halfRL = 0.5D0 )
91        PARAMETER(Half=0.5D0)        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
   
 #define I0 1  
 #define In sNx  
 #define J0 1  
 #define Jn sNy  
92  CEOP  CEOP
93    
94  ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  C Catch barotropic mode
95          IF ( Nr .LT. 2 ) RETURN
96    
97  C     Adams-Bashforth timestepping weights        iMin = 1
98        ab15=1.5+abeps        iMax = sNx
99        ab05=-0.5-abeps        jMin = 1
100          jMax = sNy
101        DO bj=myByLo(myThid),myByHi(myThid)  
102         DO bi=myBxLo(myThid),myBxHi(myThid)  C--   Initialise gW to zero
103          DO K=1,Nr        DO k=1,Nr
104           DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
105            DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
            gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
106             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
           ENDDO  
107           ENDDO           ENDDO
108          ENDDO          ENDDO
109          ENDDO
110    C-    Initialise gwDiss to zero
111          DO j=1-OLy,sNy+OLy
112           DO i=1-OLx,sNx+OLx
113             gwDiss(i,j) = 0.
114         ENDDO         ENDDO
115        ENDDO        ENDDO
116    
117  C Catch barotropic mode  C--   Boundaries condition at top
118        IF ( Nr .LT. 2 ) RETURN        DO j=1-OLy,sNy+OLy
119           DO i=1-OLx,sNx+OLx
120             flxAdvUp(i,j) = 0.
121             flxDisUp(i,j) = 0.
122           ENDDO
123          ENDDO
124    
125  C For each tile  C---  Sweep down column
126        DO bj=myByLo(myThid),myByHi(myThid)        DO k=2,Nr
127         DO bi=myBxLo(myThid),myBxHi(myThid)          kp1=k+1
128            wOverRide=1.
129  C Boundaries condition at top          IF (k.EQ.Nr) THEN
130          DO J=J0,Jn            kp1=Nr
131           DO I=I0,In            wOverRide=0.
132            Flx_Dn(I,J,bi,bj)=0.          ENDIF
133    C--   Compute grid factor arround a W-point:
134            DO j=1-Oly,sNy+Oly
135             DO i=1-Olx,sNx+Olx
136    C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
137               rThickC_W(i,j) =
138         &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
139         &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
140               rThickC_S(i,j) =
141         &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
142         &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
143               IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
144                 recip_rThickC(i,j) = 0.
145               ELSE
146                 recip_rThickC(i,j) = 1. _d 0 /
147         &        ( drF(k-1)*halfRS +
148         &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
149         &        )
150               ENDIF
151    C     W-Cell Western face area:
152               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
153    C     W-Cell Southern face area:
154               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
155           ENDDO           ENDDO
156          ENDDO          ENDDO
157    
158  C Sweep down column  C--   horizontal bi-harmonic dissipation
159          DO K=2,Nr          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
160           Kp1=K+1  C-    calculate the horizontal Laplacian of vertical flow
161           wOverRide=1.  C     Zonal flux d/dx W
162           if (K.EQ.Nr) then            DO j=1-Oly,sNy+Oly
163            Kp1=Nr             flx_EW(1-Olx,j)=0.
164            wOverRide=0.             DO i=1-Olx+1,sNx+Olx
165           endif              flx_EW(i,j) =
166  C Flux on Southern face       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
167           DO J=J0,Jn+1       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
168            DO I=I0,In  #ifdef COSINEMETH_III
169             tmp_VbarZ=Half*(       &              *sqcosFacU(j,bi,bj)
170       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)  #endif
171       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))             ENDDO
            Flx_NS(I,J,bi,bj)=  
      &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))  
      &    -viscAh*_recip_dyC(I,J,bi,bj)*(  
      &                   wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj) )  
172            ENDDO            ENDDO
173           ENDDO  C     Meridional flux d/dy W
174  C Flux on Western face            DO i=1-Olx,sNx+Olx
175           DO J=J0,Jn             flx_NS(i,1-Oly)=0.
           DO I=I0,In+1  
            tmp_UbarZ=Half*(  
      &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)  
      &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))  
            Flx_EW(I,J,bi,bj)=  
      &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))  
      &    -viscAh*_recip_dxC(I,J,bi,bj)*(  
      &                   wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj) )  
176            ENDDO            ENDDO
177           ENDDO            DO j=1-Oly+1,sNy+Oly
178  C Flux on Lower face             DO i=1-Olx,sNx+Olx
179           DO J=J0,Jn              flx_NS(i,j) =
180            DO I=I0,In       &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
181             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
182             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))  #ifdef ISOTROPIC_COS_SCALING
183             Flx_Dn(I,J,bi,bj)=  #ifdef COSINEMETH_III
184       &     tmp_WbarZ*tmp_WbarZ       &              *sqCosFacV(j,bi,bj)
185       &    -viscAr*recip_drF(K)*( wVel(I,J,K,bi,bj)  #endif
186       &                -wOverRide*wVel(I,J,Kp1,bi,bj) )  #endif
187               ENDDO
188            ENDDO            ENDDO
189           ENDDO  
190  C        Divergence of fluxes  C     del^2 W
191           DO J=J0,Jn  C     Difference of zonal fluxes ...
192            DO I=I0,In            DO j=1-Oly,sNy+Oly
193             gW(I,J,K,bi,bj) = 0.             DO i=1-Olx,sNx+Olx-1
194       &      -(              del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
195       &        +_recip_dxF(I,J,bi,bj)*(             ENDDO
196       &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )             del2w(sNx+Olx,j)=0.
      &        +_recip_dyF(I,J,bi,bj)*(  
      &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )  
      &        +recip_drC(K)         *(  
      &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )  
      &       )  
 caja    * recip_hFacU(I,J,K,bi,bj)  
 caja   NOTE:  This should be included    
 caja   but we need an hFacUW (above U points)  
 caja           and an hFacUS (above V points) too...  
197            ENDDO            ENDDO
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
198    
199        C     ... add difference of meridional fluxes and divide by volume
200        DO bj=myByLo(myThid),myByHi(myThid)            DO j=1-Oly,sNy+Oly-1
201         DO bi=myBxLo(myThid),myBxHi(myThid)             DO i=1-Olx,sNx+Olx
202          DO K=2,Nr              del2w(i,j) = ( del2w(i,j)
203           DO j=J0,Jn       &                    +(flx_NS(i,j+1)-flx_NS(i,j))
204            DO i=I0,In       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
205             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)             ENDDO
      &     +deltatMom*( ab15*gW(i,j,k,bi,bj)  
      &                 +ab05*gWNM1(i,j,k,bi,bj) )  
            IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.  
206            ENDDO            ENDDO
207           ENDDO  C-- No-slip BCs impose a drag at walls...
208          ENDDO  CML ************************************************************
209         ENDDO  CML   No-slip Boundary conditions for bi-harmonic dissipation
210        ENDDO  CML   need to be implemented here!
211    CML ************************************************************
212            ELSE
213    C-    Initialize del2w to zero:
214              DO j=1-Oly,sNy+Oly
215               DO i=1-Olx,sNx+Olx
216                del2w(i,j) = 0. _d 0
217               ENDDO
218              ENDDO
219            ENDIF
220    
221  #ifdef ALLOW_OBCS          IF (momViscosity) THEN
222        IF (useOBCS) THEN  C Viscous Flux on Western face
223  C--   This call is aesthetic: it makes the W field            DO j=jMin,jMax
224  C     consistent with the OBs but this has no algorithmic             DO i=iMin,iMax+1
225  C     impact. This is purely for diagnostic purposes.               flx_EW(i,j)=
226         DO bj=myByLo(myThid),myByHi(myThid)       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
227          DO bi=myBxLo(myThid),myBxHi(myThid)       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
228           DO K=1,Nr       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
229            CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )  cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
230           ENDDO       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
231          ENDDO       &              *(del2w(i,j)-del2w(i-1,j))
232         ENDDO       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
233        ENDIF  cOld &              *_recip_dxC(i,j,bi,bj)*drC(k)
234  #endif /* ALLOW_OBCS */  #ifdef COSINEMETH_III
235         &              *sqCosFacU(j,bi,bj)
236    #else
237         &              *CosFacU(j,bi,bj)
238    #endif
239               ENDDO
240              ENDDO
241    C Viscous Flux on Southern face
242              DO j=jMin,jMax+1
243               DO i=iMin,iMax
244                 flx_NS(i,j)=
245         &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
246         &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
247         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
248    cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
249         &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
250         &              *(del2w(i,j)-del2w(i,j-1))
251         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
252    cOld &              *_recip_dyC(i,j,bi,bj)*drC(k)
253    #ifdef ISOTROPIC_COS_SCALING
254    #ifdef COSINEMETH_III
255         &              *sqCosFacV(j,bi,bj)
256    #else
257         &              *CosFacV(j,bi,bj)
258    #endif
259    #endif
260               ENDDO
261              ENDDO
262    C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
263              DO j=jMin,jMax
264               DO i=iMin,iMax
265    C     Interpolate vert viscosity to center of tracer-cell (level k):
266                 viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
267         &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
268         &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
269         &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
270         &                 )*0.125 _d 0
271                 flx_Dn(i,j) =
272         &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
273         &                     -wVel(i,j, k ,bi,bj) )*rkSign
274         &                   *recip_drF(k)*rA(i,j,bi,bj)
275    cOld &                   *recip_drF(k)
276               ENDDO
277              ENDDO
278    C     Tendency is minus divergence of viscous fluxes:
279              DO j=jMin,jMax
280               DO i=iMin,iMax
281                 gwDiss(i,j) =
282         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
283         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
284         &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
285         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
286    cOld         gwDiss(i,j) =
287    cOld &        -(
288    cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
289    cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
290    cOld &          +                      ( flxDisUp(i,j)-flx_Dn(i,j) )
291    c    &         )*recip_rThickC(i,j)
292    cOld &         )*recip_drC(k)
293    C--        prepare for next level (k+1)
294                 flxDisUp(i,j)=flx_Dn(i,j)
295               ENDDO
296              ENDDO
297            ENDIF
298    
299            IF (no_slip_sides) THEN
300    C-     No-slip BCs impose a drag at walls...
301    c         CALL MOM_W_SIDEDRAG(
302    c    I        bi,bj,k,
303    c    O        gwAdd,
304    c    I        myThid)
305    c         DO j=jMin,jMax
306    c          DO i=iMin,iMax
307    c           gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
308    c          ENDDO
309    c         ENDDO
310            ENDIF
311    
312    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313    
314            IF ( momAdvection ) THEN
315    C Advective Flux on Western face
316              DO j=jMin,jMax
317               DO i=iMin,iMax+1
318    C     transport through Western face area:
319                 uTrans = (
320         &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)
321         &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
322         &                )*halfRL*_dyG(i,j,bi,bj)
323    cOld &                )*halfRL
324                 flx_EW(i,j)=
325         &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
326               ENDDO
327              ENDDO
328    C Advective Flux on Southern face
329              DO j=jMin,jMax+1
330               DO i=iMin,iMax
331    C     transport through Southern face area:
332                 vTrans = (
333         &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)
334         &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
335         &                )*halfRL*_dxG(i,j,bi,bj)
336    cOld &                )*halfRL
337                 flx_NS(i,j)=
338         &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
339               ENDDO
340              ENDDO
341    C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
342              DO j=jMin,jMax
343               DO i=iMin,iMax
344                 tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)
345         &                           +wVel(i,j,kp1,bi,bj)*wOverRide )
346    C     transport through Lower face area:
347                 rTrans = tmp_WbarZ*rA(i,j,bi,bj)
348                 flx_Dn(i,j) = rTrans*tmp_WbarZ
349    cOld         flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ
350               ENDDO
351              ENDDO
352    C     Tendency is minus divergence of advective fluxes:
353              DO j=jMin,jMax
354               DO i=iMin,iMax
355                 gW(i,j,k,bi,bj) =
356         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
357         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
358         &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign
359         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
360    cOld         gW(i,j,k,bi,bj) =
361    cOld &        -(
362    cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
363    cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
364    cOld &          +                      ( flxAdvUp(i,j)-flx_Dn(i,j) )
365    c    &         )*recip_rThickC(i,j)
366    cOld &         )*recip_drC(k)
367    C--          prepare for next level (k+1)
368                 flxAdvUp(i,j)=flx_Dn(i,j)
369               ENDDO
370              ENDDO
371            ENDIF
372    
373            IF ( useNHMTerms ) THEN
374              CALL MOM_W_METRIC_NH(
375         I               bi,bj,k,
376         I               uVel, vVel,
377         O               gwAdd,
378         I               myThid )
379              DO j=jMin,jMax
380               DO i=iMin,iMax
381                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
382               ENDDO
383              ENDDO
384            ENDIF
385            IF ( use3dCoriolis ) THEN
386              CALL MOM_W_CORIOLIS_NH(
387         I               bi,bj,k,
388         I               uVel, vVel,
389         O               gwAdd,
390         I               myThid )
391              DO j=jMin,jMax
392               DO i=iMin,iMax
393                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
394               ENDDO
395              ENDDO
396            ENDIF
397    
398    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
399    
400    C--   Dissipation term inside the Adams-Bashforth:
401            IF ( momViscosity .AND. momDissip_In_AB) THEN
402              DO j=jMin,jMax
403               DO i=iMin,iMax
404                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
405               ENDDO
406              ENDDO
407            ENDIF
408    
409    C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
410    C     and save gW_[n] into gwNm1 for the next time step.
411    c#ifdef ALLOW_ADAMSBASHFORTH_3
412    c       CALL ADAMS_BASHFORTH3(
413    c    I                         bi, bj, k,
414    c    U                         gW, gwNm,
415    c    I                         momStartAB, myIter, myThid )
416    c#else /* ALLOW_ADAMSBASHFORTH_3 */
417            CALL ADAMS_BASHFORTH2(
418         I                         bi, bj, k,
419         U                         gW, gwNm1,
420         I                         myIter, myThid )
421    c#endif /* ALLOW_ADAMSBASHFORTH_3 */
422    
423    C--   Dissipation term outside the Adams-Bashforth:
424            IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
425              DO j=jMin,jMax
426               DO i=iMin,iMax
427                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
428               ENDDO
429              ENDDO
430            ENDIF
431    
432    C-    end of the k loop
433          ENDDO
434    
435  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
436    
437        RETURN        RETURN
438        END        END
   
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22