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