/[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.29 by jmc, Mon Jul 10 15:17:19 2006 UTC revision 1.30 by jmc, Tue Jul 11 12:51:05 2006 UTC
# Line 13  C     !INTERFACE: Line 13  C     !INTERFACE:
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.
# Line 51  C     myThid  :: Thread number for this Line 52  C     myThid  :: Thread number for this
52    
53  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
54  C     == Local variables ==  C     == Local variables ==
55  C     iMin, iMax,  C     iMin,iMax
56  C     jMin, jMax  C     jMin,jMax
57  C     flx_NS       :: Temp. used for fVol meridional terms.  C     xA            :: W-Cell face area normal to X
58  C     flx_EW       :: Temp. used for fVol zonal terms.  C     yA            :: W-Cell face area normal to Y
59  C     flx_Up       :: Temp. used for fVol vertical terms.  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
60  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
61    C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
62    C     flx_NS        :: vertical momentum flux, meridional direction
63    C     flx_EW        :: vertical momentum flux, zonal direction
64    C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
65    C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
66    C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
67    C     gwDiss        :: vertical momentum dissipation tendency
68    C     i,j,k         :: Loop counters
69        INTEGER iMin,iMax,jMin,jMax        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)        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL    fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL    fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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)        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 C     i,j,k - Loop counters  
83        INTEGER i,j,k, kp1        INTEGER i,j,k, kp1
84        _RL  wOverride        _RL  wOverride
85        _RS  hFacWtmp        _RL  tmp_WbarZ
86        _RS  hFacStmp        _RL  uTrans, vTrans, rTrans
       _RS  hFacCtmp  
       _RS  recip_hFacCtmp  
       _RL  slipSideFac  
       _RL  tmp_VbarZ, tmp_UbarZ, tmp_WbarZ  
87        _RL  viscLoc        _RL  viscLoc
88        _RL  Half        _RL  halfRL
89        PARAMETER(Half=0.5D0)        _RS  halfRS, zeroRS
90          PARAMETER( halfRL = 0.5D0 )
91          PARAMETER( halfRS = 0.5 , zeroRS = 0. )
92  CEOP  CEOP
93    
94  C Catch barotropic mode  C Catch barotropic mode
# Line 87  C Catch barotropic mode Line 99  C Catch barotropic mode
99        jMin = 1        jMin = 1
100        jMax = sNy        jMax = sNy
101    
 C     Lateral friction (no-slip, free slip, or half slip):  
       IF ( no_slip_sides ) THEN  
         slipSideFac = -1. _d 0  
       ELSE  
         slipSideFac =  1. _d 0  
       ENDIF  
 CML   half slip was used before ; keep the line for now, but half slip is  
 CML   not used anywhere in the code as far as I can see.  
 C        slipSideFac = 0. _d 0  
 C- need to fix the side-drag implementation; for now, always impose free-slip  
         slipSideFac =  1. _d 0  
   
102  C--   Initialise gW to zero  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
# Line 107  C--   Initialise gW to zero Line 107  C--   Initialise gW to zero
107           ENDDO           ENDDO
108          ENDDO          ENDDO
109        ENDDO        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
115          ENDDO
116    
117  C--   Boundaries condition at top  C--   Boundaries condition at top
118        DO j=jMin,jMax        DO j=1-OLy,sNy+OLy
119         DO i=iMin,iMax         DO i=1-OLx,sNx+OLx
120           flx_Up(i,j)=0.           flxAdvUp(i,j) = 0.
121             flxDisUp(i,j) = 0.
122         ENDDO         ENDDO
123        ENDDO        ENDDO
124    
# Line 123  C---  Sweep down column Line 130  C---  Sweep down column
130            kp1=Nr            kp1=Nr
131            wOverRide=0.            wOverRide=0.
132          ENDIF          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
156            ENDDO
157    
158  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
159          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
160  C-    calculate the horizontal Laplacian of vertical flow  C-    calculate the horizontal Laplacian of vertical flow
161  C     Zonal flux d/dx W  C     Zonal flux d/dx W
162            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
163             fZon(1-Olx,j)=0.             flx_EW(1-Olx,j)=0.
164             DO i=1-Olx+1,sNx+Olx             DO i=1-Olx+1,sNx+Olx
165  C- Problem here: drF(k)*_hFacC & fZon are not at the same Horiz.&Vert. location              flx_EW(i,j) =
166              fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
167       &           *_dyG(i,j,bi,bj)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
      &           *_recip_dxC(i,j,bi,bj)  
      &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))  
168  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
169       &           *sqcosFacU(j,bi,bj)       &              *sqcosFacU(j,bi,bj)
170  #endif  #endif
171             ENDDO             ENDDO
172            ENDDO            ENDDO
173  C     Meridional flux d/dy W  C     Meridional flux d/dy W
174            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
175             fMer(i,1-Oly)=0.             flx_NS(i,1-Oly)=0.
176            ENDDO            ENDDO
177            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
178             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
179  C- Problem here: drF(k)*_hFacC & fMer are not at the same Horiz.&Vert. location              flx_NS(i,j) =
180              fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)       &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
181       &           *_dxG(i,j,bi,bj)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
      &           *_recip_dyC(i,j,bi,bj)  
      &           *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))  
182  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
183  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
184       &           *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
185  #endif  #endif
186  #endif  #endif
187             ENDDO             ENDDO
# Line 163  C     del^2 W Line 191  C     del^2 W
191  C     Difference of zonal fluxes ...  C     Difference of zonal fluxes ...
192            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
193             DO i=1-Olx,sNx+Olx-1             DO i=1-Olx,sNx+Olx-1
194              del2w(i,j)=fZon(i+1,j)-fZon(i,j)              del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
195             ENDDO             ENDDO
196             del2w(sNx+Olx,j)=0.             del2w(sNx+Olx,j)=0.
197            ENDDO            ENDDO
# Line 171  C     Difference of zonal fluxes ... Line 199  C     Difference of zonal fluxes ...
199  C     ... add difference of meridional fluxes and divide by volume  C     ... add difference of meridional fluxes and divide by volume
200            DO j=1-Oly,sNy+Oly-1            DO j=1-Oly,sNy+Oly-1
201             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
202  C     First compute the fraction of open water for the w-control volume              del2w(i,j) = ( del2w(i,j)
203  C     at the southern face       &                    +(flx_NS(i,j+1)-flx_NS(i,j))
204              hFacCtmp=max( _hFacC(i,j,k-1,bi,bj)-Half,0. _d 0 )       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
      &           +   min( _hFacC(i,j,k  ,bi,bj)     ,Half )  
             IF (hFacCtmp .GT. 0.) THEN  
              recip_hFacCtmp = 1./hFacCtmp  
             ELSE  
              recip_hFacCtmp = 0. _d 0  
             ENDIF  
             del2w(i,j)=recip_rA(i,j,bi,bj)  
      &           *recip_drC(k)*recip_hFacCtmp  
      &           *(  
      &           del2w(i,j)  
      &           +( fMer(i,j+1)-fMer(i,j) )  
      &           )  
205             ENDDO             ENDDO
206            ENDDO            ENDDO
207  C-- No-slip BCs impose a drag at walls...  C-- No-slip BCs impose a drag at walls...
# Line 202  C-    Initialize del2w to zero: Line 218  C-    Initialize del2w to zero:
218            ENDDO            ENDDO
219          ENDIF          ENDIF
220    
221  C Flux on Southern face          IF (momViscosity) THEN
222          DO j=jMin,jMax+1  C Viscous Flux on Western face
223           DO i=iMin,iMax            DO j=jMin,jMax
224  C     First compute the fraction of open water for the w-control volume             DO i=iMin,iMax+1
225  C     at the southern face               flx_EW(i,j)=
226             hFacStmp=max(_hFacS(i,j,k-1,bi,bj)-Half,0. _d 0)       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
227       &         +    min(_hFacS(i,j,k  ,bi,bj),Half)       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
228             tmp_VbarZ=Half*(  c    &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
229       &          _hFacS(i,j,k-1,bi,bj)*vVel( i ,j,k-1,bi,bj)       &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
230       &         +_hFacS(i,j, k ,bi,bj)*vVel( i ,j, k ,bi,bj))       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
231             flx_NS(i,j)=       &              *(del2w(i,j)-del2w(i-1,j))
232       &     tmp_VbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))  c    &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
233       &    -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*Half       &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
      &       *_recip_dyC(i,j,bi,bj)  
      &       *(hFacStmp*(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))  
 C- Problem here: No-slip bc CANNOT be written in term of a flux  
      &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)  
      &         *wVel(i,j,k,bi,bj))  
      &    +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*Half  
      &      *_recip_dyC(i,j,bi,bj)*(del2w(i,j)-del2w(i,j-1))  
 #ifdef ISOTROPIC_COS_SCALING  
234  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
235       &    *sqCosFacV(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
236  #else  #else
237       &    *CosFacV(j,bi,bj)       &              *CosFacU(j,bi,bj)
 #endif  
238  #endif  #endif
239  C     The last term is the weighted average of the viscous stress at the open             ENDDO
240  C     fraction of the w control volume and at the closed fraction of the            ENDDO
241  C     the control volume. A more compact but less intelligible version  C Viscous Flux on Southern face
242  C     of the last three lines is:            DO j=jMin,jMax+1
243  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))             DO i=iMin,iMax
244  CML     &       *wVel(i,j,k,bi,bi) + hFacStmp*wVel(i,j-1,k,bi,bj) )               flx_NS(i,j)=
245           ENDDO       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
246          ENDDO       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
247  C Flux on Western face  c    &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
248          DO j=jMin,jMax       &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
249           DO i=iMin,iMax+1       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
250  C     First compute the fraction of open water for the w-control volume       &              *(del2w(i,j)-del2w(i,j-1))
251  C     at the western face  c    &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
252             hFacWtmp=max(_hFacW(i,j,k-1,bi,bj)-Half,0. _d 0)       &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
253       &         +    min(_hFacW(i,j,k  ,bi,bj),Half)  #ifdef ISOTROPIC_COS_SCALING
            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)=  
      &     tmp_UbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))  
      &    -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*Half  
      &      *_recip_dxC(i,j,bi,bj)  
      &      *(hFacWtmp*(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))  
 C- Problem here: No-slip bc CANNOT be written in term of a flux  
      &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)  
      &        *wVel(i,j,k,bi,bj) )  
      &    +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*Half  
      &      *_recip_dxC(i,j,bi,bj)*(del2w(i,j)-del2w(i-1,j))  
254  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
255       &                *sqCosFacU(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
256  #else  #else
257       &                *CosFacU(j,bi,bj)       &              *CosFacV(j,bi,bj)
258  #endif  #endif
259  C     The last term is the weighted average of the viscous stress at the open  #endif
260  C     fraction of the w control volume and at the closed fraction of the             ENDDO
261  C     the control volume. A more compact but less intelligible version            ENDDO
262  C     of the last three lines is:  C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
263  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))            DO j=jMin,jMax
264  CML     &       *wVel(i,j,k,bi,bi) + hFacWtmp*wVel(i-1,j,k,bi,bj) )             DO i=iMin,iMax
265           ENDDO  C     Interpolate vert viscosity to center of tracer-cell (level k):
266          ENDDO               viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
267  C Flux on Lower face       &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
268          DO j=jMin,jMax       &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
269           DO i=iMin,iMax       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
270  C Interpolate vert viscosity to W points       &                 )*0.125 _d 0
271             viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)               flx_Dn(i,j) =
272       &                +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
273       &                +KappaRV(i,j,k)  +KappaRV(i,j+1,k)       &                     -wVel(i,j, k ,bi,bj) )*rkSign
274       &                +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)  c    &                   *recip_drF(k)*rA(i,j,bi,bj)
275       &               )*0.125 _d 0       &                   *recip_drF(k)
276             tmp_WbarZ = Half*( wVel(i,j, k ,bi,bj)             ENDDO
277       &                       +wVel(i,j,kp1,bi,bj)*wOverRide )            ENDDO
278             flx_Dn(i,j)=  C     Tendency is minus divergence of viscous fluxes:
279       &                  tmp_WbarZ*tmp_WbarZ            DO j=jMin,jMax
280       &                 -viscLoc*recip_drF(k)             DO i=iMin,iMax
281       &                         *( wVel(i,j, k ,bi,bj)  c            gwDiss(i,j) =
282       &                           -wVel(i,j,kp1,bi,bj)*wOverRide )  c    &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
283           ENDDO  c    &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
284          ENDDO  c    &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
285  C       Divergence of fluxes  c    &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
286          DO j=jMin,jMax               gwDiss(i,j) =
287           DO i=iMin,iMax       &        -(
288             gW(i,j,k,bi,bj) = 0.       &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
289       &      -(       &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
290       &        +_recip_dxF(i,j,bi,bj)*(       &          +                      ( flxDisUp(i,j)-flx_Dn(i,j) )
291       &              flx_EW(i+1,j)-flx_EW(i,j) )       &         )*recip_rThickC(i,j)
292       &        +_recip_dyF(i,j,bi,bj)*(  c    &         )*recip_drC(k)
      &              flx_NS(i,j+1)-flx_NS(i,j) )  
      &        +recip_drC(k)         *(  
      &              flx_Up(i,j)  -flx_Dn(i,j) )  
      &       )  
 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...  
293  C--        prepare for next level (k+1)  C--        prepare for next level (k+1)
294             flx_Up(i,j)=flx_Dn(i,j)               flxDisUp(i,j)=flx_Dn(i,j)
295           ENDDO             ENDDO
296          ENDDO            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    c    &                )*halfRL*_dyG(i,j,bi,bj)
323         &                )*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    c    &                )*halfRL*_dxG(i,j,bi,bj)
336         &                )*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    c            flx_Dn(i,j) = rTrans*tmp_WbarZ
349                 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    c            gW(i,j,k,bi,bj) =
356    c    &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
357    c    &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
358    c    &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign
359    c    &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
360                 gW(i,j,k,bi,bj) =
361         &        -(
362         &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
363         &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
364         &          +                      ( flxAdvUp(i,j)-flx_Dn(i,j) )
365         &         )*recip_rThickC(i,j)
366    c    &         )*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 ( useCoriolis ) 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-|--+----|  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)  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.  C     and save gW_[n] into gwNm1 for the next time step.
411  c#ifdef ALLOW_ADAMSBASHFORTH_3  c#ifdef ALLOW_ADAMSBASHFORTH_3
# Line 325  c#else /* ALLOW_ADAMSBASHFORTH_3 */ Line 420  c#else /* ALLOW_ADAMSBASHFORTH_3 */
420       I                         myIter, myThid )       I                         myIter, myThid )
421  c#endif /* ALLOW_ADAMSBASHFORTH_3 */  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  C-    end of the k loop
433        ENDDO        ENDDO
434    
Legend:
Removed from v.1.29  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22