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