/[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.30 by jmc, Tue Jul 11 12:51:05 2006 UTC revision 1.45 by jmc, Tue May 3 19:26:03 2011 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 30  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 58  C     xA            :: W-Cell face area Line 61  C     xA            :: W-Cell face area
61  C     yA            :: W-Cell face area normal to Y  C     yA            :: W-Cell face area normal to Y
62  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
63  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge  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)  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
66  C     flx_NS        :: vertical momentum flux, meridional direction  C     flx_NS        :: vertical momentum flux, meridional direction
67  C     flx_EW        :: vertical momentum flux, zonal direction  C     flx_EW        :: vertical momentum flux, zonal direction
# Line 65  C     flxAdvUp      :: vertical mom. adv Line 69  C     flxAdvUp      :: vertical mom. adv
69  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
70  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
71  C     gwDiss        :: vertical momentum dissipation tendency  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  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)        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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)        _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)
# Line 80  C     i,j,k         :: Loop counters Line 88  C     i,j,k         :: Loop counters
88        _RL    gwDiss(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)        _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        INTEGER i,j,k, kp1        _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL  wOverride        INTEGER i,j,k, km1, kp1
93          _RL  mskM1, mskP1
94        _RL  tmp_WbarZ        _RL  tmp_WbarZ
95        _RL  uTrans, vTrans, rTrans        _RL  uTrans, vTrans, rTrans
96        _RL  viscLoc        _RL  viscLoc
97        _RL  halfRL        _RL  halfRL
98        _RS  halfRS, zeroRS        _RS  halfRS, zeroRS
99        PARAMETER( halfRL = 0.5D0 )        PARAMETER( halfRL = 0.5 _d 0 )
100        PARAMETER( halfRS = 0.5 , zeroRS = 0. )        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 )
116          ELSE
117            diagDiss  = .FALSE.
118            diagAdvec = .FALSE.
119          ENDIF
120    #endif /* ALLOW_DIAGNOSTICS */
121    
122  C--   Initialise gW to zero  C--   Initialise gW to zero
123        DO k=1,Nr        DO k=1,Nr
# Line 113  C-    Initialise gwDiss to zero Line 133  C-    Initialise gwDiss to zero
133           gwDiss(i,j) = 0.           gwDiss(i,j) = 0.
134         ENDDO         ENDDO
135        ENDDO        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=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
147         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
148           flxAdvUp(i,j) = 0.           flxAdvUp(i,j) = 0.
149           flxDisUp(i,j) = 0.  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          ENDIF          IF ( k.EQ.Nr ) mskP1 = 0.
162            IF ( k.GT.1 ) THEN
163  C--   Compute grid factor arround a W-point:  C--   Compute grid factor arround a W-point:
164          DO j=1-Oly,sNy+Oly  #ifdef CALC_GW_NEW_THICK
165           DO i=1-Olx,sNx+Olx           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 !  C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
            rThickC_W(i,j) =  
      &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )  
      &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )  
            rThickC_S(i,j) =  
      &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )  
      &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )  
217             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
218               recip_rThickC(i,j) = 0.               recip_rThickC(i,j) = 0.
219             ELSE             ELSE
220               recip_rThickC(i,j) = 1. _d 0 /               recip_rThickC(i,j) = 1. _d 0 /
221       &        ( drF(k-1)*halfRS +       &        ( drF(k-1)*halfRS
222       &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )       &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
223       &        )       &        )
224             ENDIF             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:  C     W-Cell Western face area:
242             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
243    c    &                              *deepFacF(k)
244  C     W-Cell Southern face area:  C     W-Cell Southern face area:
245             yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)             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           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
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          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             flx_EW(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              flx_EW(i,j) =              flx_EW(i,j) =
287       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &               ( wFld(i,j) - wFld(i-1,j) )
288       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
289  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
290       &              *sqcosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
291    #endif
292    #ifdef ALLOW_OBCS
293         &              *maskInW(i,j,bi,bj)
294  #endif  #endif
295             ENDDO             ENDDO
296            ENDDO            ENDDO
297    
298  C     Meridional flux d/dy W  C     Meridional flux d/dy W
299              IF ( useCubedSphereExchange ) THEN
300    C     to compute d/dy(W), fill corners with appropriate values:
301                CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
302         &                                 wFld, bi,bj, myThid )
303              ENDIF
304            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
305             flx_NS(i,1-Oly)=0.             flx_NS(i,1-Oly)=0.
306            ENDDO            ENDDO
307            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
308             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
309              flx_NS(i,j) =              flx_NS(i,j) =
310       &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))       &               ( wFld(i,j) - wFld(i,j-1) )
311       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
312  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
313  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
314       &              *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
315  #endif  #endif
316  #endif  #endif
317    #ifdef ALLOW_OBCS
318         &              *maskInS(i,j,bi,bj)
319    #endif
320             ENDDO             ENDDO
321            ENDDO            ENDDO
322    
323  C     del^2 W  C     del^2 W
324  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)=flx_EW(i+1,j)-flx_EW(i,j)  
            ENDDO  
            del2w(sNx+Olx,j)=0.  
           ENDDO  
   
 C     ... add difference of meridional fluxes and divide by volume  
325            DO j=1-Oly,sNy+Oly-1            DO j=1-Oly,sNy+Oly-1
326             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx-1
327              del2w(i,j) = ( del2w(i,j)              del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
328       &                    +(flx_NS(i,j+1)-flx_NS(i,j))       &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
329       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
330         &                    *recip_deepFac2F(k)
331             ENDDO             ENDDO
332            ENDDO            ENDDO
333  C-- No-slip BCs impose a drag at walls...  C end if biharmonic viscosity
 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  
            ENDDO  
           ENDDO  
334          ENDIF          ENDIF
335    
336          IF (momViscosity) THEN          IF ( momViscosity .AND. k.GT.1 ) THEN
337  C Viscous Flux on Western face  C Viscous Flux on Western face
338            DO j=jMin,jMax            DO j=jMin,jMax
339             DO i=iMin,iMax+1             DO i=iMin,iMax+1
340               flx_EW(i,j)=               flx_EW(i,j)=
341       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
342       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
343  c    &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
344       &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)       &              *cosFacU(j,bi,bj)
345       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
346       &              *(del2w(i,j)-del2w(i-1,j))       &              *(del2w(i,j)-del2w(i-1,j))
347  c    &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
      &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)  
348  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
349       &              *sqCosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
350  #else  #else
351       &              *CosFacU(j,bi,bj)       &              *cosFacU(j,bi,bj)
352  #endif  #endif
353             ENDDO             ENDDO
354            ENDDO            ENDDO
# Line 244  C Viscous Flux on Southern face Line 358  C Viscous Flux on Southern face
358               flx_NS(i,j)=               flx_NS(i,j)=
359       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
360       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
361  c    &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
362       &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)  #ifdef ISOTROPIC_COS_SCALING
363         &              *cosFacV(j,bi,bj)
364    #endif
365       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
366       &              *(del2w(i,j)-del2w(i,j-1))       &              *(del2w(i,j)-del2w(i,j-1))
367  c    &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
      &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)  
368  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
369  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
370       &              *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
371  #else  #else
372       &              *CosFacV(j,bi,bj)       &              *cosFacV(j,bi,bj)
373  #endif  #endif
374  #endif  #endif
375             ENDDO             ENDDO
# Line 269  C     Interpolate vert viscosity to cent Line 384  C     Interpolate vert viscosity to cent
384       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
385       &                 )*0.125 _d 0       &                 )*0.125 _d 0
386               flx_Dn(i,j) =               flx_Dn(i,j) =
387       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
388       &                     -wVel(i,j, k ,bi,bj) )*rkSign       &                     -wVel(i,j, k ,bi,bj) )*rkSign
389  c    &                   *recip_drF(k)*rA(i,j,bi,bj)       &                   *recip_drF(k)*rA(i,j,bi,bj)
390       &                   *recip_drF(k)       &                   *deepFac2C(k)*rhoFacC(k)
391             ENDDO             ENDDO
392            ENDDO            ENDDO
393              IF ( k.EQ.2 ) THEN
394    C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
395               DO j=jMin,jMax
396                DO i=iMin,iMax
397    C     Interpolate horizontally (but not vertically) vert viscosity to center:
398    C     Although background visc. might be defined at k=1, this is not
399    C     generally true when using variable visc. (from vertical mixing scheme).
400    C     Therefore, no vert. interp. and only horizontal interpolation.
401                 viscLoc = ( KappaRU(i,j,k) + KappaRU(i+1,j,k)
402         &                  +KappaRV(i,j,k) + KappaRV(i,j+1,k)
403         &                 )*0.25 _d 0
404                 flxDisUp(i,j) =
405         &          - viscLoc*( wVel(i,j, k ,bi,bj)
406         &                     -wVel(i,j,k-1,bi,bj) )*rkSign
407         &                   *recip_drF(k-1)*rA(i,j,bi,bj)
408         &                   *deepFac2C(k-1)*rhoFacC(k-1)
409    C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
410    c            flxDisUp(i,j) = 0.
411                ENDDO
412               ENDDO
413              ENDIF
414  C     Tendency is minus divergence of viscous fluxes:  C     Tendency is minus divergence of viscous fluxes:
415    C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
416            DO j=jMin,jMax            DO j=jMin,jMax
417             DO i=iMin,iMax             DO i=iMin,iMax
 c            gwDiss(i,j) =  
 c    &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )  
 c    &           + ( flx_NS(i,j+1)-flx_NS(i,j) )  
 c    &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign  
 c    &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)  
418               gwDiss(i,j) =               gwDiss(i,j) =
419       &        -(       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
420       &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
421       &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
422       &          +                      ( flxDisUp(i,j)-flx_Dn(i,j) )       &                                          *recip_rhoFacF(k)
423       &         )*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
424  c    &         )*recip_drC(k)       &          *recip_deepFac2F(k)
425  C--        prepare for next level (k+1)  C--        prepare for next level (k+1)
426               flxDisUp(i,j)=flx_Dn(i,j)               flxDisUp(i,j)=flx_Dn(i,j)
427             ENDDO             ENDDO
428            ENDDO            ENDDO
429          ENDIF          ENDIF
430    
431          IF (no_slip_sides) THEN          IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
432  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
433  c         CALL MOM_W_SIDEDRAG(            CALL MOM_W_SIDEDRAG(
434  c    I        bi,bj,k,       I               bi,bj,k,
435  c    O        gwAdd,       I               wVel, del2w,
436  c    I        myThid)       I               rThickC_C, recip_rThickC,
437  c         DO j=jMin,jMax       I               viscAh_W, viscA4_W,
438  c          DO i=iMin,iMax       O               gwAdd,
439  c           gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)       I               myThid )
440  c          ENDDO            DO j=jMin,jMax
441  c         ENDDO             DO i=iMin,iMax
442                gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
443               ENDDO
444              ENDDO
445          ENDIF          ENDIF
446    
447  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
448    
449          IF ( momAdvection ) THEN          IF ( momAdvection ) THEN
450    
451             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
452  C Advective Flux on Western face  C Advective Flux on Western face
453            DO j=jMin,jMax            DO j=jMin,jMax
454             DO i=iMin,iMax+1             DO i=iMin,iMax+1
455  C     transport through Western face area:  C     transport through Western face area:
456               uTrans = (               uTrans = (
457       &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)       &          drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
458         &                  *rhoFacC(km1)*mskM1
459       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
460  c    &                )*halfRL*_dyG(i,j,bi,bj)       &                  *rhoFacC(k)
461       &                )*halfRL       &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
462               flx_EW(i,j)=               flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
463       &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL  c            flx_EW(i,j)=
464    c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
465             ENDDO             ENDDO
466            ENDDO            ENDDO
467  C Advective Flux on Southern face  C Advective Flux on Southern face
# Line 330  C Advective Flux on Southern face Line 469  C Advective Flux on Southern face
469             DO i=iMin,iMax             DO i=iMin,iMax
470  C     transport through Southern face area:  C     transport through Southern face area:
471               vTrans = (               vTrans = (
472       &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)       &          drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
473         &                  *rhoFacC(km1)*mskM1
474       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
475  c    &                )*halfRL*_dxG(i,j,bi,bj)       &                  *rhoFacC(k)
476       &                )*halfRL       &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
477               flx_NS(i,j)=               flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
478       &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL  c            flx_NS(i,j)=
479    c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
480             ENDDO             ENDDO
481            ENDDO            ENDDO
482             ENDIF
483  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
484    c        IF (.TRUE.) THEN
485            DO j=jMin,jMax            DO j=jMin,jMax
486             DO i=iMin,iMax             DO i=iMin,iMax
487               tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)  C     NH in p-coord.: advect wSpeed [m/s] with rTrans
488       &                           +wVel(i,j,kp1,bi,bj)*wOverRide )               tmp_WbarZ = halfRL*
489         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
490         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
491  C     transport through Lower face area:  C     transport through Lower face area:
492               rTrans = tmp_WbarZ*rA(i,j,bi,bj)               rTrans = halfRL*
493  c            flx_Dn(i,j) = rTrans*tmp_WbarZ       &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
494               flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ       &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
495         &                                   *mskP1
496         &              )*rA(i,j,bi,bj)
497                 flx_Dn(i,j) = rTrans*tmp_WbarZ
498               ENDDO
499              ENDDO
500    c        ENDIF
501             IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
502    C Advective Flux on Upper face of W-Cell (= at surface)
503               DO j=jMin,jMax
504                DO i=iMin,iMax
505                 tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
506                 rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
507         &               *rA(i,j,bi,bj)
508                 flxAdvUp(i,j) = rTrans*tmp_WbarZ
509    c            flxAdvUp(i,j) = 0.
510                ENDDO
511             ENDDO             ENDDO
512            ENDDO           ENDIF
513             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
514  C     Tendency is minus divergence of advective fluxes:  C     Tendency is minus divergence of advective fluxes:
515    C     anelastic: all transports & advect. fluxes are scaled by rhoFac
516            DO j=jMin,jMax            DO j=jMin,jMax
517             DO i=iMin,iMax             DO i=iMin,iMax
518  c            gW(i,j,k,bi,bj) =  C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
519  c    &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )  c            IF (k.EQ.2) flxAdvUp(i,j) = 0.
 c    &           + ( flx_NS(i,j+1)-flx_NS(i,j) )  
 c    &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign  
 c    &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)  
520               gW(i,j,k,bi,bj) =               gW(i,j,k,bi,bj) =
521       &        -(       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
522       &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
523       &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
524       &          +                      ( flxAdvUp(i,j)-flx_Dn(i,j) )       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
525       &         )*recip_rThickC(i,j)       &          *recip_deepFac2F(k)*recip_rhoFacF(k)
526  c    &         )*recip_drC(k)             ENDDO
527              ENDDO
528             ENDIF
529    
530             DO j=jMin,jMax
531               DO i=iMin,iMax
532  C--          prepare for next level (k+1)  C--          prepare for next level (k+1)
533               flxAdvUp(i,j)=flx_Dn(i,j)               flxAdvUp(i,j)=flx_Dn(i,j)
534             ENDDO             ENDDO
535            ENDDO           ENDDO
536    
537    c       ELSE
538    C-    if momAdvection / else
539    c         DO j=1-OLy,sNy+OLy
540    c          DO i=1-OLx,sNx+OLx
541    c            gW(i,j,k,bi,bj) = 0. _d 0
542    c          ENDDO
543    c         ENDDO
544    
545    C-    endif momAdvection.
546          ENDIF          ENDIF
547    
548          IF ( useNHMTerms ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549    
550            IF ( useNHMTerms .AND. k.GT.1 ) THEN
551            CALL MOM_W_METRIC_NH(            CALL MOM_W_METRIC_NH(
552       I               bi,bj,k,       I               bi,bj,k,
553       I               uVel, vVel,       I               uVel, vVel,
# Line 382  C--          prepare for next level (k+1 Line 559  C--          prepare for next level (k+1
559             ENDDO             ENDDO
560            ENDDO            ENDDO
561          ENDIF          ENDIF
562          IF ( useCoriolis ) THEN          IF ( use3dCoriolis .AND. k.GT.1 ) THEN
563            CALL MOM_W_CORIOLIS_NH(            CALL MOM_W_CORIOLIS_NH(
564       I               bi,bj,k,       I               bi,bj,k,
565       I               uVel, vVel,       I               uVel, vVel,
# Line 397  C--          prepare for next level (k+1 Line 574  C--          prepare for next level (k+1
574    
575  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
576    
577    #ifdef ALLOW_DIAGNOSTICS
578            IF ( diagDiss )  THEN
579              CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
580         &                           k, 1, 2, bi,bj, myThid )
581    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
582    C        does it only if k=1 (never the case here)
583    c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
584            ENDIF
585            IF ( diagAdvec ) THEN
586              CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
587         &                           k,Nr, 1, bi,bj, myThid )
588    c         IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
589            ENDIF
590    #endif /* ALLOW_DIAGNOSTICS */
591    
592  C--   Dissipation term inside the Adams-Bashforth:  C--   Dissipation term inside the Adams-Bashforth:
593          IF ( momViscosity .AND. momDissip_In_AB) THEN          IF ( momViscosity .AND. momDissip_In_AB) THEN
594            DO j=jMin,jMax            DO j=jMin,jMax
# Line 408  C--   Dissipation term inside the Adams- Line 600  C--   Dissipation term inside the Adams-
600    
601  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
602  C     and save gW_[n] into gwNm1 for the next time step.  C     and save gW_[n] into gwNm1 for the next time step.
603  c#ifdef ALLOW_ADAMSBASHFORTH_3  #ifdef ALLOW_ADAMSBASHFORTH_3
604  c       CALL ADAMS_BASHFORTH3(          CALL ADAMS_BASHFORTH3(
605  c    I                         bi, bj, k,       I                         bi, bj, k,
606  c    U                         gW, gwNm,       U                         gW, gwNm,
607  c    I                         momStartAB, myIter, myThid )       I                         nHydStartAB, myIter, myThid )
608  c#else /* ALLOW_ADAMSBASHFORTH_3 */  #else /* ALLOW_ADAMSBASHFORTH_3 */
609          CALL ADAMS_BASHFORTH2(          CALL ADAMS_BASHFORTH2(
610       I                         bi, bj, k,       I                         bi, bj, k,
611       U                         gW, gwNm1,       U                         gW, gwNm1,
612       I                         myIter, myThid )       I                         nHydStartAB, myIter, myThid )
613  c#endif /* ALLOW_ADAMSBASHFORTH_3 */  #endif /* ALLOW_ADAMSBASHFORTH_3 */
614    
615  C--   Dissipation term outside the Adams-Bashforth:  C--   Dissipation term outside the Adams-Bashforth:
616          IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN          IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
# Line 432  C--   Dissipation term outside the Adams Line 624  C--   Dissipation term outside the Adams
624  C-    end of the k loop  C-    end of the k loop
625        ENDDO        ENDDO
626    
627    #ifdef ALLOW_DIAGNOSTICS
628          IF (useDiagnostics) THEN
629            CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
630            CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
631          ENDIF
632    #endif /* ALLOW_DIAGNOSTICS */
633    
634  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
635    
636        RETURN        RETURN
Legend:
Removed from v.1.30  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.22