/[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.4 by adcroft, Thu Jun 29 18:51:18 2000 UTC revision 1.45 by jmc, Tue May 3 19:26:03 2011 UTC
# Line 1  Line 1 
1    C $Header$
2    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        SUBROUTINE CALC_GW(      CBOP
9       I        myThid)  C     !ROUTINE: CALC_GW
10  C     /==========================================================\  C     !INTERFACE:
11  C     | S/R CALC_GW                                              |        SUBROUTINE CALC_GW(
12  C     \==========================================================/       I               bi, bj, KappaRU, KappaRV,
13        IMPLICIT NONE       I               myTime, myIter, myThid )
14    C     !DESCRIPTION: \bv
15    C     *==========================================================*
16    C     | S/R CALC_GW
17    C     | o Calculate vertical velocity tendency terms
18    C     |   ( Non-Hydrostatic only )
19    C     *==========================================================*
20    C     | In NH, the vertical momentum tendency must be
21    C     | calculated explicitly and included as a source term
22    C     | for a 3d pressure eqn. Calculate that term here.
23    C     | This routine is not used in HYD calculations.
24    C     *==========================================================*
25    C     \ev
26    
27    C     !USES:
28          IMPLICIT NONE
29  C     == Global variables ==  C     == Global variables ==
30  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
31  #include "EEPARAMS.h"  #include "EEPARAMS.h"
32  #include "PARAMS.h"  #include "PARAMS.h"
33  #include "GRID.h"  #include "GRID.h"
34  #include "CG2D.h"  #include "RESTART.h"
35  #include "GW.h"  #include "SURFACE.h"
36  #include "CG3D.h"  #include "DYNVARS.h"
37    #include "NH_VARS.h"
38    
39    C     !INPUT/OUTPUT PARAMETERS:
40  C     == Routine arguments ==  C     == Routine arguments ==
41  C     myThid - Instance number for this innvocation of CALC_GW  C     bi,bj   :: current tile indices
42    C     KappaRU :: vertical viscosity at U points
43    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
51          INTEGER myIter
52        INTEGER myThid        INTEGER myThid
53    
54  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
55    
56    C     !LOCAL VARIABLES:
57  C     == Local variables ==  C     == Local variables ==
58        INTEGER bi,bj,iMin,iMax,jMin,jMax  C     iMin,iMax
59        _RL      aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     jMin,jMax
60        _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     xA            :: W-Cell face area normal to X
61        _RL      v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     yA            :: W-Cell face area normal to Y
62        _RL      cF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
63        _RL      mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
64        _RL      pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
65        _RL    fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
66        _RL    fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     flx_NS        :: vertical momentum flux, meridional direction
67    C     flx_EW        :: vertical momentum flux, zonal direction
68        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
69        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
70        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
71        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     gwDiss        :: vertical momentum dissipation tendency
72  C     I,J,K - Loop counters  C     gwAdd         :: other tendencies (Coriolis, Metric-terms)
73        INTEGER i,j,k, kP1, kUp  C     del2w         :: laplacian of wVel
74        _RL  wOverride  C     wFld          :: local copy of wVel
75        _RS  hFacROpen  C     i,j,k         :: Loop counters
76        _RS  hFacRClosed        INTEGER iMin,iMax,jMin,jMax
77        _RL ab15,ab05        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL  Half        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        PARAMETER(Half=0.5D0)        _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83  #define I0 1        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84  #define In sNx        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85  #define J0 1        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86  #define Jn sNy        _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88  C     Adams-Bashforth timestepping weights        _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        ab15=1.5+abeps        _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        ab05=-0.5-abeps        _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL    wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        DO bj=myByLo(myThid),myByHi(myThid)        INTEGER i,j,k, km1, kp1
93         DO bi=myBxLo(myThid),myBxHi(myThid)        _RL  mskM1, mskP1
94          DO K=1,Nr        _RL  tmp_WbarZ
95           DO j=1-OLy,sNy+OLy        _RL  uTrans, vTrans, rTrans
96            DO i=1-OLx,sNx+OLx        _RL  viscLoc
97          _RL  halfRL
98          _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
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---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
111    
112    #ifdef ALLOW_DIAGNOSTICS
113          IF ( useDiagnostics ) THEN
114            diagDiss  = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
115            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
123          DO k=1,Nr
124            DO j=1-OLy,sNy+OLy
125             DO i=1-OLx,sNx+OLx
126             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
           ENDDO  
127           ENDDO           ENDDO
128          ENDDO          ENDDO
129          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         ENDDO
135        ENDDO        ENDDO
136          IF (momViscosity) THEN
137  C Catch barotropic mode  C-    Initialize del2w to zero:
138        IF ( Nr .LT. 2 ) RETURN          DO j=1-Oly,sNy+Oly
139             DO i=1-Olx,sNx+Olx
140  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=J0,Jn  
          DO I=I0,In  
           Flx_Dn(I,J,bi,bj)=0.  
141           ENDDO           ENDDO
142          ENDDO          ENDDO
143          ENDIF
144    
145  C Sweep down column  C--   Boundaries condition at top (vertical advection of vertical momentum):
146          DO K=2,Nr        DO j=1-OLy,sNy+OLy
147           Kp1=K+1         DO i=1-OLx,sNx+OLx
148           wOverRide=1.           flxAdvUp(i,j) = 0.
149           if (K.EQ.Nr) then  c        flxDisUp(i,j) = 0.
           Kp1=Nr  
           wOverRide=0.  
          endif  
 C Flux on Southern face  
          DO J=J0,Jn+1  
           DO I=I0,In  
            tmp_VbarZ=Half*(  
      &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)  
      &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))  
            Flx_NS(I,J,bi,bj)=  
      &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))  
      &    -viscAh*_recip_dyC(I,J,bi,bj)*(  
      &                   wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj) )  
           ENDDO  
          ENDDO  
 C Flux on Western face  
          DO J=J0,Jn  
           DO I=I0,In+1  
            tmp_UbarZ=Half*(  
      &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)  
      &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))  
            Flx_EW(I,J,bi,bj)=  
      &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))  
      &    -viscAh*_recip_dxC(I,J,bi,bj)*(  
      &                   wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj) )  
           ENDDO  
          ENDDO  
 C Flux on Lower face  
          DO J=J0,Jn  
           DO I=I0,In  
            Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)  
            tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))  
            Flx_Dn(I,J,bi,bj)=  
      &     tmp_WbarZ*tmp_WbarZ  
      &    -viscAr*recip_drF(K)*( wVel(I,J,K,bi,bj)  
      &                -wOverRide*wVel(I,J,Kp1,bi,bj) )  
           ENDDO  
          ENDDO  
 C        Divergence of fluxes  
          DO J=J0,Jn  
           DO I=I0,In  
            gW(I,J,K,bi,bj) = 0.  
      &      -(  
      &        +_recip_dxF(I,J,bi,bj)*(  
      &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )  
      &        +_recip_dyF(I,J,bi,bj)*(  
      &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )  
      &        +recip_drC(K)         *(  
      &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )  
      &       )  
 caja    * recip_hFacU(I,J,K,bi,bj)  
 caja   NOTE:  This should be included    
 caja   but we need an hFacUW (above U points)  
 caja           and an hFacUS (above V points) too...  
           ENDDO  
          ENDDO  
         ENDDO  
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152    
153        
154        DO bj=myByLo(myThid),myByHi(myThid)  C---  Sweep down column
155         DO bi=myBxLo(myThid),myBxHi(myThid)        DO k=1,Nr
156          DO K=2,Nr          km1 = MAX( k-1, 1 )
157           DO j=J0,Jn          kp1 = MIN( k+1,Nr )
158            DO i=I0,In          mskM1 = 1.
159             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)          mskP1 = 1.
160       &     +deltatMom*( ab15*gW(i,j,k,bi,bj)          IF ( k.EQ. 1 ) mskM1 = 0.
161       &                 +ab05*gWNM1(i,j,k,bi,bj) )          IF ( k.EQ.Nr ) mskP1 = 0.
162             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.          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            ENDDO
178           ENDDO           ENDDO
179          ENDDO           IF (momViscosity) THEN
180         ENDDO           DO j=1-Oly,sNy+Oly
181        ENDDO            DO i=1-Olx,sNx+Olx
182        DO bj=myByLo(myThid),myByHi(myThid)             rThickC_C(i,j) = MAX( zeroRS,
183         DO bi=myBxLo(myThid),myBxHi(myThid)       &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
184          DO K=1,Nr       &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
185           DO j=J0,Jn       &                         )
           DO i=I0,In  
            gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
186            ENDDO            ENDDO
187           ENDDO           ENDDO
188          ENDDO           DO j=1-Oly,sNy+Oly
189         ENDDO            DO i=1-Olx+1,sNx+Olx
190        ENDDO             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
265    
266  C--   This call is aesthetic: it makes the W field  C--   local copy of wVel:
267  C     consistent with the OBs but this has no algorithmic          DO j=1-Oly,sNy+Oly
268  C     impact. This is purely for diagnostic purposes.            DO i=1-Olx,sNx+Olx
269        DO bj=myByLo(myThid),myByHi(myThid)              wFld(i,j) = wVel(i,j,k,bi,bj)
270         DO bi=myBxLo(myThid),myBxHi(myThid)            ENDDO
         DO K=1,Nr  
          Kup=max(1,K-1)  
 c        CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel,  
 c    I                        myThid )  
271          ENDDO          ENDDO
272         ENDDO  
273    C--   horizontal bi-harmonic dissipation
274            IF ( momViscosity .AND. k.GT.1 .AND. viscA4W.NE.0. ) THEN
275    
276    C-    calculate the horizontal Laplacian of vertical flow
277    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
284               flx_EW(1-Olx,j)=0.
285               DO i=1-Olx+1,sNx+Olx
286                flx_EW(i,j) =
287         &               ( wFld(i,j) - wFld(i-1,j) )
288         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
289    #ifdef COSINEMETH_III
290         &              *sqCosFacU(j,bi,bj)
291    #endif
292    #ifdef ALLOW_OBCS
293         &              *maskInW(i,j,bi,bj)
294    #endif
295               ENDDO
296              ENDDO
297    
298    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
305               flx_NS(i,1-Oly)=0.
306              ENDDO
307              DO j=1-Oly+1,sNy+Oly
308               DO i=1-Olx,sNx+Olx
309                flx_NS(i,j) =
310         &               ( wFld(i,j) - wFld(i,j-1) )
311         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
312    #ifdef ISOTROPIC_COS_SCALING
313    #ifdef COSINEMETH_III
314         &              *sqCosFacV(j,bi,bj)
315    #endif
316    #endif
317    #ifdef ALLOW_OBCS
318         &              *maskInS(i,j,bi,bj)
319    #endif
320               ENDDO
321              ENDDO
322    
323    C     del^2 W
324    C     Divergence of horizontal fluxes
325              DO j=1-Oly,sNy+Oly-1
326               DO i=1-Olx,sNx+Olx-1
327                del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
328         &                    +( flx_NS(i,j+1)-flx_NS(i,j) )
329         &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
330         &                    *recip_deepFac2F(k)
331               ENDDO
332              ENDDO
333    C end if biharmonic viscosity
334            ENDIF
335    
336            IF ( momViscosity .AND. k.GT.1 ) THEN
337    C Viscous Flux on Western face
338              DO j=jMin,jMax
339               DO i=iMin,iMax+1
340                 flx_EW(i,j)=
341         &       - (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))
343         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
344         &              *cosFacU(j,bi,bj)
345         &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
346         &              *(del2w(i,j)-del2w(i-1,j))
347         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
348    #ifdef COSINEMETH_III
349         &              *sqCosFacU(j,bi,bj)
350    #else
351         &              *cosFacU(j,bi,bj)
352    #endif
353               ENDDO
354              ENDDO
355    C Viscous Flux on Southern face
356              DO j=jMin,jMax+1
357               DO i=iMin,iMax
358                 flx_NS(i,j)=
359         &       - (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))
361         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
362    #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
366         &              *(del2w(i,j)-del2w(i,j-1))
367         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
368    #ifdef ISOTROPIC_COS_SCALING
369    #ifdef COSINEMETH_III
370         &              *sqCosFacV(j,bi,bj)
371    #else
372         &              *cosFacV(j,bi,bj)
373    #endif
374    #endif
375               ENDDO
376              ENDDO
377    C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
378              DO j=jMin,jMax
379               DO i=iMin,iMax
380    C     Interpolate vert viscosity to center of tracer-cell (level k):
381                 viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
382         &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
383         &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
384         &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
385         &                 )*0.125 _d 0
386                 flx_Dn(i,j) =
387         &          - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
388         &                     -wVel(i,j, k ,bi,bj) )*rkSign
389         &                   *recip_drF(k)*rA(i,j,bi,bj)
390         &                   *deepFac2C(k)*rhoFacC(k)
391               ENDDO
392              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:
415    C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
416              DO j=jMin,jMax
417               DO i=iMin,iMax
418                 gwDiss(i,j) =
419         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
420         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
421         &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
422         &                                          *recip_rhoFacF(k)
423         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
424         &          *recip_deepFac2F(k)
425    C--        prepare for next level (k+1)
426                 flxDisUp(i,j)=flx_Dn(i,j)
427               ENDDO
428              ENDDO
429            ENDIF
430    
431            IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
432    C-     No-slip BCs impose a drag at walls...
433              CALL MOM_W_SIDEDRAG(
434         I               bi,bj,k,
435         I               wVel, del2w,
436         I               rThickC_C, recip_rThickC,
437         I               viscAh_W, viscA4_W,
438         O               gwAdd,
439         I               myThid )
440              DO j=jMin,jMax
441               DO i=iMin,iMax
442                gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
443               ENDDO
444              ENDDO
445            ENDIF
446    
447    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
448    
449            IF ( momAdvection ) THEN
450    
451             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
452    C Advective Flux on Western face
453              DO j=jMin,jMax
454               DO i=iMin,iMax+1
455    C     transport through Western face area:
456                 uTrans = (
457         &          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)
460         &                  *rhoFacC(k)
461         &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
462                 flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
463    c            flx_EW(i,j)=
464    c    &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
465               ENDDO
466              ENDDO
467    C Advective Flux on Southern face
468              DO j=jMin,jMax+1
469               DO i=iMin,iMax
470    C     transport through Southern face area:
471                 vTrans = (
472         &          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)
475         &                  *rhoFacC(k)
476         &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
477                 flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
478    c            flx_NS(i,j)=
479    c    &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
480               ENDDO
481              ENDDO
482             ENDIF
483    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
486               DO i=iMin,iMax
487    C     NH in p-coord.: advect wSpeed [m/s] with rTrans
488                 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:
492                 rTrans = halfRL*
493         &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
494         &               +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
512             ENDIF
513             IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
514    C     Tendency is minus divergence of advective fluxes:
515    C     anelastic: all transports & advect. fluxes are scaled by rhoFac
516              DO j=jMin,jMax
517               DO i=iMin,iMax
518    C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
519    c            IF (k.EQ.2) flxAdvUp(i,j) = 0.
520                 gW(i,j,k,bi,bj) =
521         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
522         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
523         &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
524         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
525         &          *recip_deepFac2F(k)*recip_rhoFacF(k)
526               ENDDO
527              ENDDO
528             ENDIF
529    
530             DO j=jMin,jMax
531               DO i=iMin,iMax
532    C--          prepare for next level (k+1)
533                 flxAdvUp(i,j)=flx_Dn(i,j)
534               ENDDO
535             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
547    
548    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549    
550            IF ( useNHMTerms .AND. k.GT.1 ) THEN
551              CALL MOM_W_METRIC_NH(
552         I               bi,bj,k,
553         I               uVel, vVel,
554         O               gwAdd,
555         I               myThid )
556              DO j=jMin,jMax
557               DO i=iMin,iMax
558                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
559               ENDDO
560              ENDDO
561            ENDIF
562            IF ( use3dCoriolis .AND. k.GT.1 ) THEN
563              CALL MOM_W_CORIOLIS_NH(
564         I               bi,bj,k,
565         I               uVel, vVel,
566         O               gwAdd,
567         I               myThid )
568              DO j=jMin,jMax
569               DO i=iMin,iMax
570                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
571               ENDDO
572              ENDDO
573            ENDIF
574    
575    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:
593            IF ( momViscosity .AND. momDissip_In_AB) THEN
594              DO j=jMin,jMax
595               DO i=iMin,iMax
596                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
597               ENDDO
598              ENDDO
599            ENDIF
600    
601    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.
603    #ifdef ALLOW_ADAMSBASHFORTH_3
604            CALL ADAMS_BASHFORTH3(
605         I                         bi, bj, k,
606         U                         gW, gwNm,
607         I                         nHydStartAB, myIter, myThid )
608    #else /* ALLOW_ADAMSBASHFORTH_3 */
609            CALL ADAMS_BASHFORTH2(
610         I                         bi, bj, k,
611         U                         gW, gwNm1,
612         I                         nHydStartAB, myIter, myThid )
613    #endif /* ALLOW_ADAMSBASHFORTH_3 */
614    
615    C--   Dissipation term outside the Adams-Bashforth:
616            IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
617              DO j=jMin,jMax
618               DO i=iMin,iMax
619                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
620               ENDDO
621              ENDDO
622            ENDIF
623    
624    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
637        END        END
   
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.22