/[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.28 by jmc, Fri Jul 7 20:09:37 2006 UTC revision 1.43 by jmc, Sat Jan 23 00:04:03 2010 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
 C     !DESCRIPTION: \bv  
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #define CALC_GW_NEW_THICK
7    
8  CBOP  CBOP
9  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
# Line 13  C     !INTERFACE: Line 14  C     !INTERFACE:
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | S/R CALC_GW  C     | S/R CALC_GW
17  C     | o Calculate vert. velocity tendency terms ( NH, QH only )  C     | o Calculate vertical velocity tendency terms
18    C     |   ( Non-Hydrostatic only )
19  C     *==========================================================*  C     *==========================================================*
20  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH, the vertical momentum tendency must be
21  C     | calculated explicitly and included as a source term  C     | calculated explicitly and included as a source term
22  C     | for a 3d pressure eqn. Calculate that term here.  C     | for a 3d pressure eqn. Calculate that term here.
23  C     | This routine is not used in HYD calculations.  C     | This routine is not used in HYD calculations.
# Line 29  C     == Global variables == Line 31  C     == Global variables ==
31  #include "EEPARAMS.h"  #include "EEPARAMS.h"
32  #include "PARAMS.h"  #include "PARAMS.h"
33  #include "GRID.h"  #include "GRID.h"
34    #include "RESTART.h"
35    #include "SURFACE.h"
36  #include "DYNVARS.h"  #include "DYNVARS.h"
37  #include "NH_VARS.h"  #include "NH_VARS.h"
38    
# Line 51  C     myThid  :: Thread number for this Line 55  C     myThid  :: Thread number for this
55    
56  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
57  C     == Local variables ==  C     == Local variables ==
58  C     iMin, iMax,  C     iMin,iMax
59  C     jMin, jMax  C     jMin,jMax
60  C     flx_NS       :: Temp. used for fVol meridional terms.  C     xA            :: W-Cell face area normal to X
61  C     flx_EW       :: Temp. used for fVol zonal terms.  C     yA            :: W-Cell face area normal to Y
62  C     flx_Up       :: Temp. used for fVol vertical terms.  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
63  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
64    C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
65    C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
66    C     flx_NS        :: vertical momentum flux, meridional direction
67    C     flx_EW        :: vertical momentum flux, zonal direction
68    C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
69    C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
70    C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
71    C     gwDiss        :: vertical momentum dissipation tendency
72    C     gwAdd         :: other tendencies (Coriolis, Metric-terms)
73    C     del2w         :: laplacian of wVel
74    C     wFld          :: local copy of wVel
75    C     i,j,k         :: Loop counters
76        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
77          _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RL    fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RL    fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91  C     i,j,k - Loop counters        _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        INTEGER i,j,k, kp1        INTEGER i,j,k, km1, kp1
93        _RL  wOverride        _RL  mskM1, mskP1
94        _RS  hFacWtmp        _RL  tmp_WbarZ
95        _RS  hFacStmp        _RL  uTrans, vTrans, rTrans
       _RS  hFacCtmp  
       _RS  recip_hFacCtmp  
       _RL  slipSideFac  
       _RL  tmp_VbarZ, tmp_UbarZ, tmp_WbarZ  
96        _RL  viscLoc        _RL  viscLoc
97        _RL  Half        _RL  halfRL
98        PARAMETER(Half=0.5D0)        _RS  halfRS, zeroRS
99          PARAMETER( halfRL = 0.5 _d 0 )
100          PARAMETER( halfRS = 0.5 , zeroRS = 0. )
101          PARAMETER( iMin = 1 , iMax = sNx )
102          PARAMETER( jMin = 1 , jMax = sNy )
103  CEOP  CEOP
104    #ifdef ALLOW_DIAGNOSTICS
105          LOGICAL diagDiss, diagAdvec
106          LOGICAL  DIAGNOSTICS_IS_ON
107          EXTERNAL DIAGNOSTICS_IS_ON
108    #endif /* ALLOW_DIAGNOSTICS */
109    
110  C Catch barotropic mode  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       IF ( Nr .LT. 2 ) RETURN  
111    
112        iMin = 1  #ifdef ALLOW_DIAGNOSTICS
113        iMax = sNx        IF ( useDiagnostics ) THEN
114        jMin = 1          diagDiss  = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
115        jMax = sNy          diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid )
   
 C     Lateral friction (no-slip, free slip, or half slip):  
       IF ( no_slip_sides ) THEN  
         slipSideFac = -1. _d 0  
116        ELSE        ELSE
117          slipSideFac =  1. _d 0          diagDiss  = .FALSE.
118            diagAdvec = .FALSE.
119        ENDIF        ENDIF
120  CML   half slip was used before ; keep the line for now, but half slip is  #endif /* ALLOW_DIAGNOSTICS */
 CML   not used anywhere in the code as far as I can see.  
 C        slipSideFac = 0. _d 0  
121    
122  C--   Initialise gW to zero  C--   Initialise gW to zero
123        DO k=1,Nr        DO k=1,Nr
# Line 105  C--   Initialise gW to zero Line 127  C--   Initialise gW to zero
127           ENDDO           ENDDO
128          ENDDO          ENDDO
129        ENDDO        ENDDO
130    C-    Initialise gwDiss to zero
131          DO j=1-OLy,sNy+OLy
132           DO i=1-OLx,sNx+OLx
133             gwDiss(i,j) = 0.
134           ENDDO
135          ENDDO
136          IF (momViscosity) THEN
137    C-    Initialize del2w to zero:
138            DO j=1-Oly,sNy+Oly
139             DO i=1-Olx,sNx+Olx
140               del2w(i,j) = 0. _d 0
141             ENDDO
142            ENDDO
143          ENDIF
144    
145  C--   Boundaries condition at top  C--   Boundaries condition at top (vertical advection of vertical momentum):
146        DO j=jMin,jMax        DO j=1-OLy,sNy+OLy
147         DO i=iMin,iMax         DO i=1-OLx,sNx+OLx
148           flx_Up(i,j)=0.           flxAdvUp(i,j) = 0.
149    c        flxDisUp(i,j) = 0.
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152    
153    
154  C---  Sweep down column  C---  Sweep down column
155        DO k=2,Nr        DO k=1,Nr
156          kp1=k+1          km1 = MAX( k-1, 1 )
157          wOverRide=1.          kp1 = MIN( k+1,Nr )
158          IF (k.EQ.Nr) THEN          mskM1 = 1.
159            kp1=Nr          mskP1 = 1.
160            wOverRide=0.          IF ( k.EQ. 1 ) mskM1 = 0.
161            IF ( k.EQ.Nr ) mskP1 = 0.
162            IF ( k.GT.1 ) THEN
163    C--   Compute grid factor arround a W-point:
164    #ifdef CALC_GW_NEW_THICK
165             DO j=1-Oly,sNy+Oly
166              DO i=1-Olx,sNx+Olx
167               IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
168         &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
169                 recip_rThickC(i,j) = 0.
170               ELSE
171    C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
172                 recip_rThickC(i,j) = 1. _d 0 /
173         &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
174         &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
175         &        )
176               ENDIF
177              ENDDO
178             ENDDO
179             IF (momViscosity) THEN
180             DO j=1-Oly,sNy+Oly
181              DO i=1-Olx,sNx+Olx
182               rThickC_C(i,j) = MAX( zeroRS,
183         &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
184         &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
185         &                         )
186              ENDDO
187             ENDDO
188             DO j=1-Oly,sNy+Oly
189              DO i=1-Olx+1,sNx+Olx
190               rThickC_W(i,j) = MAX( zeroRS,
191         &                           MIN( rSurfW(i,j,bi,bj), rC(k-1) )
192         &                          -MAX(  rLowW(i,j,bi,bj), rC(k)   )
193         &                         )
194    C     W-Cell Western face area:
195               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
196    c    &                              *deepFacF(k)
197              ENDDO
198             ENDDO
199             DO j=1-Oly+1,sNy+Oly
200              DO i=1-Olx,sNx+Olx
201               rThickC_S(i,j) = MAX( zeroRS,
202         &                           MIN( rSurfS(i,j,bi,bj), rC(k-1) )
203         &                          -MAX(  rLowS(i,j,bi,bj), rC(k)   )
204         &                         )
205    C     W-Cell Southern face area:
206               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
207    c    &                              *deepFacF(k)
208    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
209    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
210              ENDDO
211             ENDDO
212             ENDIF
213    #else /* CALC_GW_NEW_THICK */
214             DO j=1-Oly,sNy+Oly
215              DO i=1-Olx,sNx+Olx
216    C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
217               IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
218                 recip_rThickC(i,j) = 0.
219               ELSE
220                 recip_rThickC(i,j) = 1. _d 0 /
221         &        ( drF(k-1)*halfRS
222         &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
223         &        )
224               ENDIF
225    c          IF (momViscosity) THEN
226    #ifdef NONLIN_FRSURF
227               rThickC_C(i,j) =
228         &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
229         &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
230    #else
231               rThickC_C(i,j) =
232         &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
233         &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
234    #endif
235               rThickC_W(i,j) =
236         &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
237         &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
238               rThickC_S(i,j) =
239         &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
240         &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
241    C     W-Cell Western face area:
242               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
243    c    &                              *deepFacF(k)
244    C     W-Cell Southern face area:
245               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
246    c    &                              *deepFacF(k)
247    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
248    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
249    c          ENDIF
250              ENDDO
251             ENDDO
252    #endif /* CALC_GW_NEW_THICK */
253            ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
254             DO j=1-Oly,sNy+Oly
255              DO i=1-Olx,sNx+Olx
256               recip_rThickC(i,j) = recip_drC(k)
257    c          rThickC_C(i,j) = drC(k)
258    c          rThickC_W(i,j) = drC(k)
259    c          rThickC_S(i,j) = drC(k)
260    c          xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
261    c          yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
262              ENDDO
263             ENDDO
264          ENDIF          ENDIF
265    
266    C--   local copy of wVel:
267            DO j=1-Oly,sNy+Oly
268              DO i=1-Olx,sNx+Olx
269                wFld(i,j) = wVel(i,j,k,bi,bj)
270              ENDDO
271            ENDDO
272    
273  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
274          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF ( momViscosity .AND. k.GT.1 .AND. viscA4W.NE.0. ) THEN
275    
276  C-    calculate the horizontal Laplacian of vertical flow  C-    calculate the horizontal Laplacian of vertical flow
277  C     Zonal flux d/dx W  C     Zonal flux d/dx W
278              IF ( useCubedSphereExchange ) THEN
279    C     to compute d/dx(W), fill corners with appropriate values:
280                CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
281         &                                 wFld, bi,bj, myThid )
282              ENDIF
283            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
284             fZon(1-Olx,j)=0.             flx_EW(1-Olx,j)=0.
285             DO i=1-Olx+1,sNx+Olx             DO i=1-Olx+1,sNx+Olx
286  C- Problem here: drF(k)*_hFacC & fZon are not at the same Horiz.&Vert. location              flx_EW(i,j) =
287              fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)       &               ( wFld(i,j) - wFld(i-1,j) )
288       &           *_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))  
289  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
290       &           *sqcosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
291  #endif  #endif
292             ENDDO             ENDDO
293            ENDDO            ENDDO
294    
295  C     Meridional flux d/dy W  C     Meridional flux d/dy W
296              IF ( useCubedSphereExchange ) THEN
297    C     to compute d/dy(W), fill corners with appropriate values:
298                CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
299         &                                 wFld, bi,bj, myThid )
300              ENDIF
301            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
302             fMer(i,1-Oly)=0.             flx_NS(i,1-Oly)=0.
303            ENDDO            ENDDO
304            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
305             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
306  C- Problem here: drF(k)*_hFacC & fMer are not at the same Horiz.&Vert. location              flx_NS(i,j) =
307              fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)       &               ( wFld(i,j) - wFld(i,j-1) )
308       &           *_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))  
309  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
310  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
311       &           *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
312  #endif  #endif
313  #endif  #endif
314             ENDDO             ENDDO
315            ENDDO            ENDDO
316    
317  C     del^2 W  C     del^2 W
318  C     Difference of zonal fluxes ...  C     Divergence of horizontal fluxes
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx-1  
             del2w(i,j)=fZon(i+1,j)-fZon(i,j)  
            ENDDO  
            del2w(sNx+Olx,j)=0.  
           ENDDO  
   
 C     ... add difference of meridional fluxes and divide by volume  
319            DO j=1-Oly,sNy+Oly-1            DO j=1-Oly,sNy+Oly-1
320             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx-1
321  C     First compute the fraction of open water for the w-control volume              del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
322  C     at the southern face       &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
323              hFacCtmp=max( _hFacC(i,j,k-1,bi,bj)-Half,0. _d 0 )       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
324       &           +   min( _hFacC(i,j,k  ,bi,bj)     ,Half )       &                    *recip_deepFac2F(k)
             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) )  
      &           )  
            ENDDO  
           ENDDO  
 C-- No-slip BCs impose a drag at walls...  
 CML ************************************************************  
 CML   No-slip Boundary conditions for bi-harmonic dissipation  
 CML   need to be implemented here!  
 CML ************************************************************  
         ELSE  
 C-    Initialize del2w to zero:  
           DO j=1-Oly,sNy+Oly  
            DO i=1-Olx,sNx+Olx  
             del2w(i,j) = 0. _d 0  
325             ENDDO             ENDDO
326            ENDDO            ENDDO
327    C end if biharmonic viscosity
328          ENDIF          ENDIF
329    
330  C Flux on Southern face          IF ( momViscosity .AND. k.GT.1 ) THEN
331          DO j=jMin,jMax+1  C Viscous Flux on Western face
332           DO i=iMin,iMax            DO j=jMin,jMax
333  C     First compute the fraction of open water for the w-control volume             DO i=iMin,iMax+1
334  C     at the southern face               flx_EW(i,j)=
335             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
336       &         +    min(_hFacS(i,j,k  ,bi,bj),Half)       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
337             tmp_VbarZ=Half*(       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
338       &          _hFacS(i,j,k-1,bi,bj)*vVel( i ,j,k-1,bi,bj)       &              *cosFacU(j,bi,bj)
339       &         +_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
340             flx_NS(i,j)=       &              *(del2w(i,j)-del2w(i-1,j))
341       &     tmp_VbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
      &    -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*Half  
      &       *_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  
342  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
343       &    *sqCosFacV(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
344  #else  #else
345       &    *CosFacV(j,bi,bj)       &              *cosFacU(j,bi,bj)
346  #endif  #endif
347               ENDDO
348              ENDDO
349    C Viscous Flux on Southern face
350              DO j=jMin,jMax+1
351               DO i=iMin,iMax
352                 flx_NS(i,j)=
353         &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
354         &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
355         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
356    #ifdef ISOTROPIC_COS_SCALING
357         &              *cosFacV(j,bi,bj)
358  #endif  #endif
359  C     The last term is the weighted average of the viscous stress at the open       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
360  C     fraction of the w control volume and at the closed fraction of the       &              *(del2w(i,j)-del2w(i,j-1))
361  C     the control volume. A more compact but less intelligible version       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
362  C     of the last three lines is:  #ifdef ISOTROPIC_COS_SCALING
 CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))  
 CML     &       *wVel(i,j,k,bi,bi) + hFacStmp*wVel(i,j-1,k,bi,bj) )  
          ENDDO  
         ENDDO  
 C Flux on Western face  
         DO j=jMin,jMax  
          DO i=iMin,iMax+1  
 C     First compute the fraction of open water for the w-control volume  
 C     at the western face  
            hFacWtmp=max(_hFacW(i,j,k-1,bi,bj)-Half,0. _d 0)  
      &         +    min(_hFacW(i,j,k  ,bi,bj),Half)  
            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))  
363  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
364       &                *sqCosFacU(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
365  #else  #else
366       &                *CosFacU(j,bi,bj)       &              *cosFacV(j,bi,bj)
367  #endif  #endif
368  C     The last term is the weighted average of the viscous stress at the open  #endif
369  C     fraction of the w control volume and at the closed fraction of the             ENDDO
370  C     the control volume. A more compact but less intelligible version            ENDDO
371  C     of the last three lines is:  C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
372  CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))            DO j=jMin,jMax
373  CML     &       *wVel(i,j,k,bi,bi) + hFacWtmp*wVel(i-1,j,k,bi,bj) )             DO i=iMin,iMax
374           ENDDO  C     Interpolate vert viscosity to center of tracer-cell (level k):
375          ENDDO               viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
376  C Flux on Lower face       &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
377          DO j=jMin,jMax       &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
378           DO i=iMin,iMax       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
379  C Interpolate vert viscosity to W points       &                 )*0.125 _d 0
380             viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)               flx_Dn(i,j) =
381       &                +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
382       &                +KappaRV(i,j,k)  +KappaRV(i,j+1,k)       &                     -wVel(i,j, k ,bi,bj) )*rkSign
383       &                +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)       &                   *recip_drF(k)*rA(i,j,bi,bj)
384       &               )*0.125 _d 0       &                   *deepFac2C(k)*rhoFacC(k)
385             tmp_WbarZ = Half*( wVel(i,j, k ,bi,bj)             ENDDO
386       &                       +wVel(i,j,kp1,bi,bj)*wOverRide )            ENDDO
387             flx_Dn(i,j)=            IF ( k.EQ.2 ) THEN
388       &                  tmp_WbarZ*tmp_WbarZ  C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
389       &                 -viscLoc*recip_drF(k)             DO j=jMin,jMax
390       &                         *( wVel(i,j, k ,bi,bj)              DO i=iMin,iMax
391       &                           -wVel(i,j,kp1,bi,bj)*wOverRide )  C     Interpolate horizontally (but not vertically) vert viscosity to center:
392           ENDDO  C     Although background visc. might be defined at k=1, this is not
393          ENDDO  C     generally true when using variable visc. (from vertical mixing scheme).
394  C       Divergence of fluxes  C     Therefore, no vert. interp. and only horizontal interpolation.
395          DO j=jMin,jMax               viscLoc = ( KappaRU(i,j,k) + KappaRU(i+1,j,k)
396           DO i=iMin,iMax       &                  +KappaRV(i,j,k) + KappaRV(i,j+1,k)
397             gW(i,j,k,bi,bj) = 0.       &                 )*0.25 _d 0
398       &      -(               flxDisUp(i,j) =
399       &        +_recip_dxF(i,j,bi,bj)*(       &          - viscLoc*( wVel(i,j, k ,bi,bj)
400       &              flx_EW(i+1,j)-flx_EW(i,j) )       &                     -wVel(i,j,k-1,bi,bj) )*rkSign
401       &        +_recip_dyF(i,j,bi,bj)*(       &                   *recip_drF(k-1)*rA(i,j,bi,bj)
402       &              flx_NS(i,j+1)-flx_NS(i,j) )       &                   *deepFac2C(k-1)*rhoFacC(k-1)
403       &        +recip_drC(k)         *(  C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
404       &              flx_Up(i,j)  -flx_Dn(i,j) )  c            flxDisUp(i,j) = 0.
405       &       )              ENDDO
406  caja    * recip_hFacU(i,j,k,bi,bj)             ENDDO
407  caja   NOTE:  This should be included            ENDIF
408  caja   but we need an hFacUW (above U points)  C     Tendency is minus divergence of viscous fluxes:
409  caja           and an hFacUS (above V points) too...  C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
410              DO j=jMin,jMax
411               DO i=iMin,iMax
412                 gwDiss(i,j) =
413         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
414         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
415         &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
416         &                                          *recip_rhoFacF(k)
417         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
418         &          *recip_deepFac2F(k)
419  C--        prepare for next level (k+1)  C--        prepare for next level (k+1)
420             flx_Up(i,j)=flx_Dn(i,j)               flxDisUp(i,j)=flx_Dn(i,j)
421               ENDDO
422              ENDDO
423            ENDIF
424    
425            IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
426    C-     No-slip BCs impose a drag at walls...
427              CALL MOM_W_SIDEDRAG(
428         I               bi,bj,k,
429         I               wVel, del2w,
430         I               rThickC_C, recip_rThickC,
431         I               viscAh_W, viscA4_W,
432         O               gwAdd,
433         I               myThid )
434              DO j=jMin,jMax
435               DO i=iMin,iMax
436                gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
437               ENDDO
438              ENDDO
439            ENDIF
440    
441    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
442    
443            IF ( momAdvection ) THEN
444    
445             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
446    C Advective Flux on Western face
447              DO j=jMin,jMax
448               DO i=iMin,iMax+1
449    C     transport through Western face area:
450                 uTrans = (
451         &          drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
452         &                  *rhoFacC(km1)*mskM1
453         &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
454         &                  *rhoFacC(k)
455         &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
456                 flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
457    c            flx_EW(i,j)=
458    c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
459               ENDDO
460              ENDDO
461    C Advective Flux on Southern face
462              DO j=jMin,jMax+1
463               DO i=iMin,iMax
464    C     transport through Southern face area:
465                 vTrans = (
466         &          drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
467         &                  *rhoFacC(km1)*mskM1
468         &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
469         &                  *rhoFacC(k)
470         &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
471                 flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
472    c            flx_NS(i,j)=
473    c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
474               ENDDO
475              ENDDO
476             ENDIF
477    C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
478    c        IF (.TRUE.) THEN
479              DO j=jMin,jMax
480               DO i=iMin,iMax
481    C     NH in p-coord.: advect wSpeed [m/s] with rTrans
482                 tmp_WbarZ = halfRL*
483         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
484         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
485    C     transport through Lower face area:
486                 rTrans = halfRL*
487         &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
488         &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
489         &                                   *mskP1
490         &              )*rA(i,j,bi,bj)
491                 flx_Dn(i,j) = rTrans*tmp_WbarZ
492               ENDDO
493              ENDDO
494    c        ENDIF
495             IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
496    C Advective Flux on Upper face of W-Cell (= at surface)
497               DO j=jMin,jMax
498                DO i=iMin,iMax
499                 tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
500                 rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
501         &               *rA(i,j,bi,bj)
502                 flxAdvUp(i,j) = rTrans*tmp_WbarZ
503    c            flxAdvUp(i,j) = 0.
504                ENDDO
505               ENDDO
506             ENDIF
507             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
508    C     Tendency is minus divergence of advective fluxes:
509    C     anelastic: all transports & advect. fluxes are scaled by rhoFac
510              DO j=jMin,jMax
511               DO i=iMin,iMax
512                 gW(i,j,k,bi,bj) =
513         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
514         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
515         &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
516         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
517         &          *recip_deepFac2F(k)*recip_rhoFacF(k)
518               ENDDO
519              ENDDO
520             ENDIF
521    
522             DO j=jMin,jMax
523               DO i=iMin,iMax
524    C--          prepare for next level (k+1)
525                 flxAdvUp(i,j)=flx_Dn(i,j)
526               ENDDO
527           ENDDO           ENDDO
528          ENDDO  
529    c       ELSE
530    C-    if momAdvection / else
531    c         DO j=1-OLy,sNy+OLy
532    c          DO i=1-OLx,sNx+OLx
533    c            gW(i,j,k,bi,bj) = 0. _d 0
534    c          ENDDO
535    c         ENDDO
536    
537    C-    endif momAdvection.
538            ENDIF
539    
540    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
541    
542            IF ( useNHMTerms .AND. k.GT.1 ) THEN
543              CALL MOM_W_METRIC_NH(
544         I               bi,bj,k,
545         I               uVel, vVel,
546         O               gwAdd,
547         I               myThid )
548              DO j=jMin,jMax
549               DO i=iMin,iMax
550                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
551               ENDDO
552              ENDDO
553            ENDIF
554            IF ( use3dCoriolis .AND. k.GT.1 ) THEN
555              CALL MOM_W_CORIOLIS_NH(
556         I               bi,bj,k,
557         I               uVel, vVel,
558         O               gwAdd,
559         I               myThid )
560              DO j=jMin,jMax
561               DO i=iMin,iMax
562                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
563               ENDDO
564              ENDDO
565            ENDIF
566    
567  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
568    
569    #ifdef ALLOW_DIAGNOSTICS
570            IF ( diagDiss )  THEN
571              CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
572         &                           k, 1, 2, bi,bj, myThid )
573    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
574    C        does it only if k=1 (never the case here)
575    c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
576            ENDIF
577            IF ( diagAdvec ) THEN
578              CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
579         &                           k,Nr, 1, bi,bj, myThid )
580    c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
581            ENDIF
582    #endif /* ALLOW_DIAGNOSTICS */
583    
584    C--   Dissipation term inside the Adams-Bashforth:
585            IF ( momViscosity .AND. momDissip_In_AB) THEN
586              DO j=jMin,jMax
587               DO i=iMin,iMax
588                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
589               ENDDO
590              ENDDO
591            ENDIF
592    
593  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
594  C     and save gW_[n] into gwNm1 for the next time step.  C     and save gW_[n] into gwNm1 for the next time step.
595  c#ifdef ALLOW_ADAMSBASHFORTH_3  #ifdef ALLOW_ADAMSBASHFORTH_3
596  c       CALL ADAMS_BASHFORTH3(          CALL ADAMS_BASHFORTH3(
597  c    I                         bi, bj, k,       I                         bi, bj, k,
598  c    U                         gW, gwNm,       U                         gW, gwNm,
599  c    I                         momStartAB, myIter, myThid )       I                         nHydStartAB, myIter, myThid )
600  c#else /* ALLOW_ADAMSBASHFORTH_3 */  #else /* ALLOW_ADAMSBASHFORTH_3 */
601          CALL ADAMS_BASHFORTH2(          CALL ADAMS_BASHFORTH2(
602       I                         bi, bj, k,       I                         bi, bj, k,
603       U                         gW, gwNm1,       U                         gW, gwNm1,
604       I                         myIter, myThid )       I                         nHydStartAB, myIter, myThid )
605  c#endif /* ALLOW_ADAMSBASHFORTH_3 */  #endif /* ALLOW_ADAMSBASHFORTH_3 */
606    
607    C--   Dissipation term outside the Adams-Bashforth:
608            IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
609              DO j=jMin,jMax
610               DO i=iMin,iMax
611                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
612               ENDDO
613              ENDDO
614            ENDIF
615    
616  C-    end of the k loop  C-    end of the k loop
617        ENDDO        ENDDO
618    
619    #ifdef ALLOW_DIAGNOSTICS
620          IF (useDiagnostics) THEN
621            CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
622            CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
623          ENDIF
624    #endif /* ALLOW_DIAGNOSTICS */
625    
626  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
627    
628        RETURN        RETURN
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.22