/[MITgcm]/MITgcm/pkg/mom_fluxform/mom_fluxform.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_fluxform/mom_fluxform.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.31 by heimbach, Thu Dec 8 15:44:34 2005 UTC revision 1.52 by jmc, Wed Dec 10 22:21:17 2014 UTC
# Line 26  C stresses as well as internal viscous s Line 26  C stresses as well as internal viscous s
26  CEOI  CEOI
27    
28  #include "MOM_FLUXFORM_OPTIONS.h"  #include "MOM_FLUXFORM_OPTIONS.h"
29    #ifdef ALLOW_AUTODIFF
30    # include "AUTODIFF_OPTIONS.h"
31    #endif
32    #ifdef ALLOW_MOM_COMMON
33    # include "MOM_COMMON_OPTIONS.h"
34    #endif
35    
36  CBOP  CBOP
37  C !ROUTINE: MOM_FLUXFORM  C !ROUTINE: MOM_FLUXFORM
38    
39  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
40        SUBROUTINE MOM_FLUXFORM(        SUBROUTINE MOM_FLUXFORM(
41       I        bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,       I        bi,bj,k,iMin,iMax,jMin,jMax,
42       I        KappaRU, KappaRV,       I        KappaRU, KappaRV,
43       U        fVerU, fVerV,       U        fVerUkm, fVerVkm,
44         O        fVerUkp, fVerVkp,
45       O        guDiss, gvDiss,       O        guDiss, gvDiss,
46       I        myTime, myIter, myThid)       I        myTime, myIter, myThid )
47    
48  C !DESCRIPTION:  C !DESCRIPTION:
49  C Calculates all the horizontal accelerations except for the implicit surface  C Calculates all the horizontal accelerations except for the implicit surface
50  C pressure gradient and implciit vertical viscosity.  C pressure gradient and implicit vertical viscosity.
51    
52  C !USES: ===============================================================  C !USES: ===============================================================
53  C     == Global variables ==  C     == Global variables ==
54        IMPLICIT NONE        IMPLICIT NONE
55  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
56  #include "EEPARAMS.h"  #include "EEPARAMS.h"
57  #include "PARAMS.h"  #include "PARAMS.h"
58  #include "GRID.h"  #include "GRID.h"
59    #include "DYNVARS.h"
60    #include "FFIELDS.h"
61  #include "SURFACE.h"  #include "SURFACE.h"
62    #ifdef ALLOW_MOM_COMMON
63    # include "MOM_VISC.h"
64    #endif
65    #ifdef ALLOW_AUTODIFF
66    # include "tamc.h"
67    # include "tamc_keys.h"
68    # include "MOM_FLUXFORM.h"
69    #endif
70    
71  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
72  C  bi,bj                :: tile indices  C  bi,bj                :: current tile indices
73  C  iMin,iMax,jMin,jMAx  :: loop ranges  C  k                    :: current vertical level
74  C  k                    :: vertical level  C  iMin,iMax,jMin,jMax  :: loop ranges
 C  kUp                  :: =1 or 2 for consecutive k  
 C  kDown                :: =2 or 1 for consecutive k  
75  C  KappaRU              :: vertical viscosity  C  KappaRU              :: vertical viscosity
76  C  KappaRV              :: vertical viscosity  C  KappaRV              :: vertical viscosity
77  C  fVerU                :: vertical flux of U, 2 1/2 dim for pipe-lining  C  fVerUkm              :: vertical advective flux of U, interface above (k-1/2)
78  C  fVerV                :: vertical flux of V, 2 1/2 dim for pipe-lining  C  fVerVkm              :: vertical advective flux of V, interface above (k-1/2)
79    C  fVerUkp              :: vertical advective flux of U, interface below (k+1/2)
80    C  fVerVkp              :: vertical advective flux of V, interface below (k+1/2)
81  C  guDiss               :: dissipation tendency (all explicit terms), u component  C  guDiss               :: dissipation tendency (all explicit terms), u component
82  C  gvDiss               :: dissipation tendency (all explicit terms), v component  C  gvDiss               :: dissipation tendency (all explicit terms), v component
83  C  myTime               :: current time  C  myTime               :: current time
84  C  myIter               :: current time-step number  C  myIter               :: current time-step number
85  C  myThid               :: thread number  C  myThid               :: my Thread Id number
86        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,k
87        INTEGER k,kUp,kDown        INTEGER iMin,iMax,jMin,jMax
88        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerUkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerVkm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL fVerUkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          _RL fVerVkp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL     myTime        _RL     myTime
# Line 91  C  cF                   :: Coriolis acce Line 108  C  cF                   :: Coriolis acce
108  C  mT                   :: Metric terms  C  mT                   :: Metric terms
109  C  fZon                 :: zonal fluxes  C  fZon                 :: zonal fluxes
110  C  fMer                 :: meridional fluxes  C  fMer                 :: meridional fluxes
111  C  fVrUp,fVrDw          :: vertical viscous fluxes at interface k-1 & k  C  fVrUp,fVrDw          :: vertical viscous fluxes at interface k & k+1
112        INTEGER i,j        INTEGER i,j
113    #ifdef ALLOW_AUTODIFF_TAMC
114          INTEGER imomkey
115    #endif
116        _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 101  C  fVrUp,fVrDw          :: vertical visc Line 121  C  fVrUp,fVrDw          :: vertical visc
121        _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124  C     afFacMom      - Tracer parameters for turning terms  C     afFacMom     :: Tracer parameters for turning terms on and off.
125  C     vfFacMom        on and off.  C     vfFacMom
126  C     pfFacMom        afFacMom - Advective terms  C     pfFacMom        afFacMom - Advective terms
127  C     cfFacMom        vfFacMom - Eddy viscosity terms  C     cfFacMom        vfFacMom - Eddy viscosity terms
128  C     mTFacMom        pfFacMom - Pressure terms  C     mtFacMom        pfFacMom - Pressure terms
129  C                     cfFacMom - Coriolis terms  C                     cfFacMom - Coriolis terms
130  C                     foFacMom - Forcing  C                     foFacMom - Forcing
131  C                     mTFacMom - Metric term  C                     mtFacMom - Metric term
132  C     uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off  C     uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off
133        _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RS   h0FacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS  r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RS      xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS      xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RS      yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS      yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 137  C     uDudxFac, AhDudxFac, etc ... indiv Line 158  C     uDudxFac, AhDudxFac, etc ... indiv
158        _RL  ArDudrFac        _RL  ArDudrFac
159        _RL  fuFac        _RL  fuFac
160        _RL  mtFacU        _RL  mtFacU
161          _RL  mtNHFacU
162        _RL  uDvdxFac        _RL  uDvdxFac
163        _RL  AhDvdxFac        _RL  AhDvdxFac
164        _RL  vDvdyFac        _RL  vDvdyFac
# Line 145  C     uDudxFac, AhDudxFac, etc ... indiv Line 167  C     uDudxFac, AhDudxFac, etc ... indiv
167        _RL  ArDvdrFac        _RL  ArDvdrFac
168        _RL  fvFac        _RL  fvFac
169        _RL  mtFacV        _RL  mtFacV
170          _RL  mtNHFacV
171        _RL  sideMaskFac        _RL  sideMaskFac
172        LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity        LOGICAL bottomDragTerms
173  CEOP  CEOP
174    #ifdef MOM_BOUNDARY_CONSERVE
175          COMMON / MOM_FLUXFORM_LOCAL / uBnd, vBnd
176          _RL  uBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
177          _RL  vBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
178    #endif /* MOM_BOUNDARY_CONSERVE */
179    
180    #ifdef ALLOW_AUTODIFF_TAMC
181              act0 = k - 1
182              max0 = Nr
183              act1 = bi - myBxLo(myThid)
184              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
185              act2 = bj - myByLo(myThid)
186              max2 = myByHi(myThid) - myByLo(myThid) + 1
187              act3 = myThid - 1
188              max3 = nTx*nTy
189              act4 = ikey_dynamics - 1
190              imomkey = (act0 + 1)
191         &                    + act1*max0
192         &                    + act2*max0*max1
193         &                    + act3*max0*max1*max2
194         &                    + act4*max0*max1*max2*max3
195    #endif /* ALLOW_AUTODIFF_TAMC */
196    
197  C     Initialise intermediate terms  C     Initialise intermediate terms
198        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 162  C     Initialise intermediate terms Line 207  C     Initialise intermediate terms
207          fVrDw(i,j)= 0.          fVrDw(i,j)= 0.
208          rTransU(i,j)= 0.          rTransU(i,j)= 0.
209          rTransV(i,j)= 0.          rTransV(i,j)= 0.
210    c       KE(i,j)     = 0.
211            hDiv(i,j)   = 0.
212            vort3(i,j)  = 0.
213          strain(i,j) = 0.          strain(i,j) = 0.
214          tension(i,j)= 0.          tension(i,j)= 0.
215          guDiss(i,j) = 0.          guDiss(i,j) = 0.
216          gvDiss(i,j) = 0.          gvDiss(i,j) = 0.
 #ifdef ALLOW_AUTODIFF_TAMC  
         vort3(i,j)   = 0. _d 0  
         strain(i,j)  = 0. _d 0  
         tension(i,j) = 0. _d 0  
 #endif  
217         ENDDO         ENDDO
218        ENDDO        ENDDO
219    
# Line 182  C     o U momentum equation Line 225  C     o U momentum equation
225        AhDudyFac    = vfFacMom*1.        AhDudyFac    = vfFacMom*1.
226        rVelDudrFac  = afFacMom*1.        rVelDudrFac  = afFacMom*1.
227        ArDudrFac    = vfFacMom*1.        ArDudrFac    = vfFacMom*1.
228        mTFacU       = mtFacMom*1.        mtFacU       = mtFacMom*1.
229          mtNHFacU     = 1.
230        fuFac        = cfFacMom*1.        fuFac        = cfFacMom*1.
231  C     o V momentum equation  C     o V momentum equation
232        uDvdxFac     = afFacMom*1.        uDvdxFac     = afFacMom*1.
# Line 191  C     o V momentum equation Line 235  C     o V momentum equation
235        AhDvdyFac    = vfFacMom*1.        AhDvdyFac    = vfFacMom*1.
236        rVelDvdrFac  = afFacMom*1.        rVelDvdrFac  = afFacMom*1.
237        ArDvdrFac    = vfFacMom*1.        ArDvdrFac    = vfFacMom*1.
238        mTFacV       = mtFacMom*1.        mtFacV       = mtFacMom*1.
239          mtNHFacV     = 1.
240        fvFac        = cfFacMom*1.        fvFac        = cfFacMom*1.
241    
242        IF (implicitViscosity) THEN        IF (implicitViscosity) THEN
# Line 216  C       vorticity at a no-slip boundary Line 261  C       vorticity at a no-slip boundary
261        ENDIF        ENDIF
262    
263  C--   Calculate open water fraction at vorticity points  C--   Calculate open water fraction at vorticity points
264        CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)        CALL MOM_CALC_HFACZ( bi,bj,k,hFacZ,r_hFacZ,myThid )
265    
266  C---- Calculate common quantities used in both U and V equations  C---- Calculate common quantities used in both U and V equations
267  C     Calculate tracer cell face open areas  C     Calculate tracer cell face open areas
268        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
269         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
270          xA(i,j) = _dyG(i,j,bi,bj)          xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
271       &   *drF(k)*_hFacW(i,j,k,bi,bj)       &          *drF(k)*_hFacW(i,j,k,bi,bj)
272          yA(i,j) = _dxG(i,j,bi,bj)          yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
273       &   *drF(k)*_hFacS(i,j,k,bi,bj)       &          *drF(k)*_hFacS(i,j,k,bi,bj)
274            h0FacZ(i,j) = hFacZ(i,j)
275         ENDDO         ENDDO
276        ENDDO        ENDDO
277    #ifdef NONLIN_FRSURF
278          IF ( momViscosity .AND. no_slip_sides
279         &                  .AND. nonlinFreeSurf.GT.0 ) THEN
280            DO j=2-OLy,sNy+OLy
281             DO i=2-OLx,sNx+OLx
282              h0FacZ(i,j) = MIN(
283         &       MIN( h0FacW(i,j,k,bi,bj), h0FacW(i,j-1,k,bi,bj) ),
284         &       MIN( h0FacS(i,j,k,bi,bj), h0FacS(i-1,j,k,bi,bj) ) )
285             ENDDO
286            ENDDO
287           ENDIF
288    #endif /* NONLIN_FRSURF */
289    
290  C     Make local copies of horizontal flow field  C     Make local copies of horizontal flow field
291        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 238  C     Make local copies of horizontal fl Line 296  C     Make local copies of horizontal fl
296        ENDDO        ENDDO
297    
298  C     Calculate velocity field "volume transports" through tracer cell faces.  C     Calculate velocity field "volume transports" through tracer cell faces.
299    C     anelastic: transports are scaled by rhoFacC (~ mass transport)
300        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
301         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
302          uTrans(i,j) = uFld(i,j)*xA(i,j)          uTrans(i,j) = uFld(i,j)*xA(i,j)*rhoFacC(k)
303          vTrans(i,j) = vFld(i,j)*yA(i,j)          vTrans(i,j) = vFld(i,j)*yA(i,j)*rhoFacC(k)
304         ENDDO         ENDDO
305        ENDDO        ENDDO
306    
307        CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)        CALL MOM_CALC_KE( bi,bj,k,2,uFld,vFld,KE,myThid )
308        IF ( momViscosity) THEN        IF ( useVariableVisc ) THEN
309          CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)          CALL MOM_CALC_HDIV( bi,bj,k,2,uFld,vFld,hDiv,myThid )
310          CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)          CALL MOM_CALC_RELVORT3( bi,bj,k,uFld,vFld,hFacZ,vort3,myThid )
311          CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)          CALL MOM_CALC_TENSION( bi,bj,k,uFld,vFld,tension,myThid )
312          CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)          CALL MOM_CALC_STRAIN( bi,bj,k,uFld,vFld,hFacZ,strain,myThid )
313          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
314           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
315             IF ( hFacZ(i,j).EQ.0. ) THEN             IF ( hFacZ(i,j).EQ.0. ) THEN
316               vort3(i,j)  = sideMaskFac*vort3(i,j)               vort3(i,j)  = sideMaskFac*vort3(i,j)
317               strain(i,j) = sideMaskFac*strain(i,j)               strain(i,j) = sideMaskFac*strain(i,j)
# Line 269  C     Calculate velocity field "volume t Line 328  C     Calculate velocity field "volume t
328  #endif  #endif
329        ENDIF        ENDIF
330    
331  C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)  C---  First call (k=1): compute vertical adv. flux fVerUkm & fVerVkm
332        IF (momAdvection.AND.k.EQ.1) THEN        IF (momAdvection.AND.k.EQ.1) THEN
333    
334    #ifdef MOM_BOUNDARY_CONSERVE
335            CALL MOM_UV_BOUNDARY( bi, bj, k,
336         I                        uVel, vVel,
337         O                        uBnd(1-OLx,1-OLy,k,bi,bj),
338         O                        vBnd(1-OLx,1-OLy,k,bi,bj),
339         I                        myTime, myIter, myThid )
340    #endif /* MOM_BOUNDARY_CONSERVE */
341    
342  C-    Calculate vertical transports above U & V points (West & South face):  C-    Calculate vertical transports above U & V points (West & South face):
343    
344    #ifdef ALLOW_AUTODIFF_TAMC
345    # ifdef NONLIN_FRSURF
346    #  ifndef DISABLE_RSTAR_CODE
347    CADJ STORE dwtransc(:,:,bi,bj) =
348    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
349    CADJ STORE dwtransu(:,:,bi,bj) =
350    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
351    CADJ STORE dwtransv(:,:,bi,bj) =
352    CADJ &     comlev1_bibj_k, key = imomkey, byte = isbyte
353    #  endif
354    # endif /* NONLIN_FRSURF */
355    #endif /* ALLOW_AUTODIFF_TAMC */
356          CALL MOM_CALC_RTRANS( k, bi, bj,          CALL MOM_CALC_RTRANS( k, bi, bj,
357       O                        rTransU, rTransV,       O                        rTransU, rTransV,
358       I                        myTime, myIter, myThid)       I                        myTime, myIter, myThid )
359    
360  C-    Free surface correction term (flux at k=1)  C-    Free surface correction term (flux at k=1)
361          CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,          CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
362       O                     fVerU(1-OLx,1-OLy,kUp), myThid )       O                     fVerUkm, myThid )
363    
364          CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,          CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
365       O                     fVerV(1-OLx,1-OLy,kUp), myThid )       O                     fVerVkm, myThid )
366    
367  C---  endif momAdvection & k=1  C---  endif momAdvection & k=1
368        ENDIF        ENDIF
369    
   
370  C---  Calculate vertical transports (at k+1) below U & V points :  C---  Calculate vertical transports (at k+1) below U & V points :
371        IF (momAdvection) THEN        IF (momAdvection) THEN
372          CALL MOM_CALC_RTRANS( k+1, bi, bj,          CALL MOM_CALC_RTRANS( k+1, bi, bj,
373       O                        rTransU, rTransV,       O                        rTransU, rTransV,
374       I                        myTime, myIter, myThid)       I                        myTime, myIter, myThid )
375          ENDIF
376    
377    #ifdef MOM_BOUNDARY_CONSERVE
378          IF ( momAdvection .AND. k.LT.Nr ) THEN
379            CALL MOM_UV_BOUNDARY( bi, bj, k+1,
380         I                        uVel, vVel,
381         O                        uBnd(1-OLx,1-OLy,k+1,bi,bj),
382         O                        vBnd(1-OLx,1-OLy,k+1,bi,bj),
383         I                        myTime, myIter, myThid )
384        ENDIF        ENDIF
385    #endif /* MOM_BOUNDARY_CONSERVE */
386    
387        IF (momViscosity) THEN        IF (momViscosity) THEN
388         CALL MOM_CALC_VISC(         DO j=1-OLy,sNy+OLy
389       I        bi,bj,k,          DO i=1-OLx,sNx+OLx
390       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,           viscAh_D(i,j) = viscAhD
391       O        harmonic,biharmonic,useVariableViscosity,           viscAh_Z(i,j) = viscAhZ
392       I        hDiv,vort3,tension,strain,KE,hFacZ,           viscA4_D(i,j) = viscA4D
393       I        myThid)           viscA4_Z(i,j) = viscA4Z
394            ENDDO
395           ENDDO
396           IF ( useVariableVisc ) THEN
397            CALL MOM_CALC_VISC( bi, bj, k,
398         O           viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
399         I           hDiv, vort3, tension, strain, KE, hFacZ,
400         I           myThid )
401           ENDIF
402        ENDIF        ENDIF
403    
404  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 311  C---- Zonal momentum equation starts her Line 408  C---- Zonal momentum equation starts her
408        IF (momAdvection) THEN        IF (momAdvection) THEN
409  C---  Calculate mean fluxes (advection)   between cells for zonal flow.  C---  Calculate mean fluxes (advection)   between cells for zonal flow.
410    
411    #ifdef MOM_BOUNDARY_CONSERVE
412            CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
413         O                     fZon,myThid )
414            CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uBnd(1-OLx,1-OLy,k,bi,bj),
415         O                     fMer,myThid )
416            CALL MOM_U_ADV_WU(
417         I                     bi,bj,k+1,uBnd,wVel,rTransU,
418         O                     fVerUkp, myThid )
419    #else /* MOM_BOUNDARY_CONSERVE */
420  C--   Zonal flux (fZon is at east face of "u" cell)  C--   Zonal flux (fZon is at east face of "u" cell)
421  C     Mean flow component of zonal flux -> fZon  C     Mean flow component of zonal flux -> fZon
422          CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,fZon,myThid)          CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uFld,fZon,myThid )
423    
424  C--   Meridional flux (fMer is at south face of "u" cell)  C--   Meridional flux (fMer is at south face of "u" cell)
425  C     Mean flow component of meridional flux -> fMer  C     Mean flow component of meridional flux -> fMer
426          CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,fMer,myThid)          CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uFld,fMer,myThid )
427    
428  C--   Vertical flux (fVer is at upper face of "u" cell)  C--   Vertical flux (fVer is at upper face of "u" cell)
429  C     Mean flow component of vertical flux (at k+1) -> fVer  C     Mean flow component of vertical flux (at k+1) -> fVer
430          CALL MOM_U_ADV_WU(          CALL MOM_U_ADV_WU(
431       I                     bi,bj,k+1,uVel,wVel,rTransU,       I                     bi,bj,k+1,uVel,wVel,rTransU,
432       O                     fVerU(1-OLx,1-OLy,kDown), myThid )       O                     fVerUkp, myThid )
433    #endif /* MOM_BOUNDARY_CONSERVE */
434    
435  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
436          DO j=jMin,jMax          DO j=jMin,jMax
# Line 334  C--   Tendency is minus divergence of th Line 441  C--   Tendency is minus divergence of th
441       &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )       &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
442  #else  #else
443       &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)       &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
444       &     *recip_rAw(i,j,bi,bj)       &     *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
445  #endif  #endif
446       &    *( ( fZon(i,j  )     - fZon(i-1,j) )*uDudxFac       &     *( ( fZon(i,j  )  - fZon(i-1,j)  )*uDudxFac
447       &      +( fMer(i,j+1)     - fMer(i,  j) )*vDudyFac       &       +( fMer(i,j+1)  - fMer(i,  j)  )*vDudyFac
448       &      +(fVerU(i,j,kDown) - fVerU(i,j,kUp))*rkSign*rVelDudrFac       &       +( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign*rVelDudrFac
449       &     )       &     )
450           ENDDO           ENDDO
451          ENDDO          ENDDO
452    
453  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
454          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
455            CALL DIAGNOSTICS_FILL(fZon,'ADVx_Um ',k,1,2,bi,bj,myThid)            CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Um ',k,1,2,bi,bj,myThid)
456            CALL DIAGNOSTICS_FILL(fMer,'ADVy_Um ',k,1,2,bi,bj,myThid)            CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Um ',k,1,2,bi,bj,myThid)
457            CALL DIAGNOSTICS_FILL(fVerU(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL(fVerUkm,'ADVrE_Um',k,1,2,bi,bj,myThid)
      &                               'ADVrE_Um',k,1,2,bi,bj,myThid)  
458          ENDIF          ENDIF
459  #endif  #endif
460    
# Line 359  C-- account for 3.D divergence of the fl Line 465  C-- account for 3.D divergence of the fl
465           DO j=jMin,jMax           DO j=jMin,jMax
466            DO i=iMin,iMax            DO i=iMin,iMax
467             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
468       &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf       &     - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
469       &       *uVel(i,j,k,bi,bj)       &       *uVel(i,j,k,bi,bj)
470            ENDDO            ENDDO
471           ENDDO           ENDDO
# Line 375  C-- account for 3.D divergence of the fl Line 481  C-- account for 3.D divergence of the fl
481  # endif /* DISABLE_RSTAR_CODE */  # endif /* DISABLE_RSTAR_CODE */
482  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
483    
484    #ifdef ALLOW_ADDFLUID
485            IF ( selectAddFluid.GE.1 ) THEN
486             DO j=jMin,jMax
487              DO i=iMin,iMax
488               gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
489         &     + uVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
490         &      *( addMass(i-1,j,k,bi,bj) + addMass(i,j,k,bi,bj) )
491         &      *_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
492         &      * recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
493              ENDDO
494             ENDDO
495            ENDIF
496    #endif /* ALLOW_ADDFLUID */
497    
498        ELSE        ELSE
499  C-    if momAdvection / else  C-    if momAdvection / else
500          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 390  C-    endif momAdvection. Line 510  C-    endif momAdvection.
510  C---  Calculate eddy fluxes (dissipation) between cells for zonal flow.  C---  Calculate eddy fluxes (dissipation) between cells for zonal flow.
511    
512  C     Bi-harmonic term del^2 U -> v4F  C     Bi-harmonic term del^2 U -> v4F
513          IF (biharmonic)          IF ( useBiharmonicVisc )
514       &  CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)       &  CALL MOM_U_DEL2U( bi, bj, k, uFld, hFacZ, h0FacZ,
515         O                    v4f, myThid )
516    
517  C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon  C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
518          CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,          CALL MOM_U_XVISCFLUX( bi,bj,k,uFld,v4F,fZon,
519       I    viscAh_D,viscA4_D,myThid)       I                        viscAh_D,viscA4_D,myThid )
520    
521  C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer  C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
522          CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,          CALL MOM_U_YVISCFLUX( bi,bj,k,uFld,v4F,hFacZ,fMer,
523       I    viscAh_Z,viscA4_Z,myThid)       I                        viscAh_Z,viscA4_Z,myThid )
524    
525  C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw  C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
526         IF (.NOT.implicitViscosity) THEN         IF (.NOT.implicitViscosity) THEN
527          CALL MOM_U_RVISCFLUX(bi,bj, k, uVel,KappaRU,fVrUp,myThid)          CALL MOM_U_RVISCFLUX( bi,bj, k, uVel,KappaRU,fVrUp,myThid )
528          CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,fVrDw,myThid)          CALL MOM_U_RVISCFLUX( bi,bj,k+1,uVel,KappaRU,fVrDw,myThid )
529         ENDIF         ENDIF
530    
531  C--   Tendency is minus divergence of the fluxes  C--   Tendency is minus divergence of the fluxes
532    C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
533          DO j=jMin,jMax          DO j=jMin,jMax
534           DO i=iMin,iMax           DO i=iMin,iMax
535            guDiss(i,j) =            guDiss(i,j) =
# Line 416  C--   Tendency is minus divergence of th Line 538  C--   Tendency is minus divergence of th
538       &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )       &      ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
539  #else  #else
540       &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)       &     -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
541       &     *recip_rAw(i,j,bi,bj)       &     *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
542  #endif  #endif
543       &    *( ( fZon(i,j  ) - fZon(i-1,j) )*AhDudxFac       &     *( ( fZon(i,j  ) - fZon(i-1,j) )*AhDudxFac
544       &      +( fMer(i,j+1) - fMer(i,  j) )*AhDudyFac       &       +( fMer(i,j+1) - fMer(i,  j) )*AhDudyFac
545       &      +( fVrDw(i,j)  - fVrUp(i,j) )*rkSign*ArDudrFac       &       +( fVrDw(i,j)  - fVrUp(i,j)  )*rkSign*ArDudrFac
546         &                                     *recip_rhoFacC(k)
547       &     )       &     )
548           ENDDO           ENDDO
549          ENDDO          ENDDO
# Line 434  C--   Tendency is minus divergence of th Line 557  C--   Tendency is minus divergence of th
557          ENDIF          ENDIF
558  #endif  #endif
559    
560  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
561          IF (no_slip_sides) THEN          IF (no_slip_sides) THEN
562  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
563           CALL MOM_U_SIDEDRAG(           CALL MOM_U_SIDEDRAG( bi, bj, k,
564       I        bi,bj,k,       I        uFld, v4f, h0FacZ,
565       I        uFld, v4f, hFacZ,       I        viscAh_Z, viscA4_Z,
566       I        viscAh_Z,viscA4_Z,       I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
      I        harmonic,biharmonic,useVariableViscosity,  
567       O        vF,       O        vF,
568       I        myThid)       I        myThid )
569           DO j=jMin,jMax           DO j=jMin,jMax
570            DO i=iMin,iMax            DO i=iMin,iMax
571             gUdiss(i,j) = gUdiss(i,j) + vF(i,j)             gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
# Line 452  C-     No-slip BCs impose a drag at wall Line 574  C-     No-slip BCs impose a drag at wall
574          ENDIF          ENDIF
575  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
576          IF (bottomDragTerms) THEN          IF (bottomDragTerms) THEN
577           CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)           CALL MOM_U_BOTTOMDRAG( bi,bj,k,uFld,KE,KappaRU,vF,myThid )
578             DO j=jMin,jMax
579              DO i=iMin,iMax
580               gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
581              ENDDO
582             ENDDO
583            ENDIF
584    
585    #ifdef ALLOW_SHELFICE
586            IF (useShelfIce) THEN
587             CALL SHELFICE_U_DRAG( bi,bj,k,uFld,KE,KappaRU,vF,myThid )
588           DO j=jMin,jMax           DO j=jMin,jMax
589            DO i=iMin,iMax            DO i=iMin,iMax
590             gUdiss(i,j) = gUdiss(i,j) + vF(i,j)             gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
591            ENDDO            ENDDO
592           ENDDO           ENDDO
593          ENDIF          ENDIF
594    #endif /* ALLOW_SHELFICE */
595    
596  C-    endif momViscosity  C-    endif momViscosity
597        ENDIF        ENDIF
# Line 471  c    I     myTime,myThid) Line 604  c    I     myTime,myThid)
604    
605  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
606        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
607  C      o Non-hydrosatic metric terms  C      o Non-Hydrostatic (spherical) metric terms
608         CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)         CALL MOM_U_METRIC_NH( bi,bj,k,uFld,wVel,mT,myThid )
609         DO j=jMin,jMax         DO j=jMin,jMax
610          DO i=iMin,iMax          DO i=iMin,iMax
611           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j)
612          ENDDO          ENDDO
613         ENDDO         ENDDO
614        ENDIF        ENDIF
615        IF (usingSphericalPolarMTerms) THEN        IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
616         CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)  C      o Spherical polar grid metric terms
617           CALL MOM_U_METRIC_SPHERE( bi,bj,k,uFld,vFld,mT,myThid )
618         DO j=jMin,jMax         DO j=jMin,jMax
619          DO i=iMin,iMax          DO i=iMin,iMax
620           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
621          ENDDO          ENDDO
622         ENDDO         ENDDO
623        ENDIF        ENDIF
624        IF (usingCylindricalGrid) THEN        IF ( usingCylindricalGrid .AND. metricTerms ) THEN
625           CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)  C      o Cylindrical grid metric terms
626           DO j=jMin,jMax         CALL MOM_U_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid )
627            DO i=iMin,iMax         DO j=jMin,jMax
628               gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)          DO i=iMin,iMax
629            ENDDO           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
630            ENDDO
631         ENDDO         ENDDO
632        ENDIF        ENDIF
633    
# Line 501  C---+----1----+----2----+----3----+----4 Line 636  C---+----1----+----2----+----3----+----4
636  C---- Meridional momentum equation starts here  C---- Meridional momentum equation starts here
637    
638        IF (momAdvection) THEN        IF (momAdvection) THEN
639    
640    #ifdef MOM_BOUNDARY_CONSERVE
641            CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
642         O                     fZon,myThid )
643            CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vBnd(1-OLx,1-OLy,k,bi,bj),
644         O                     fMer,myThid )
645            CALL MOM_V_ADV_WV( bi,bj,k+1,vBnd,wVel,rTransV,
646         O                     fVerVkp, myThid )
647    #else /* MOM_BOUNDARY_CONSERVE */
648  C---  Calculate mean fluxes (advection)   between cells for meridional flow.  C---  Calculate mean fluxes (advection)   between cells for meridional flow.
649  C     Mean flow component of zonal flux -> fZon  C     Mean flow component of zonal flux -> fZon
650          CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,fZon,myThid)          CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vFld,fZon,myThid )
651    
652  C--   Meridional flux (fMer is at north face of "v" cell)  C--   Meridional flux (fMer is at north face of "v" cell)
653  C     Mean flow component of meridional flux -> fMer  C     Mean flow component of meridional flux -> fMer
654          CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,fMer,myThid)          CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vFld,fMer,myThid )
655    
656  C--   Vertical flux (fVer is at upper face of "v" cell)  C--   Vertical flux (fVer is at upper face of "v" cell)
657  C     Mean flow component of vertical flux (at k+1) -> fVerV  C     Mean flow component of vertical flux (at k+1) -> fVerV
658          CALL MOM_V_ADV_WV(          CALL MOM_V_ADV_WV( bi,bj,k+1,vVel,wVel,rTransV,
659       I                     bi,bj,k+1,vVel,wVel,rTransV,       O                     fVerVkp, myThid )
660       O                     fVerV(1-OLx,1-OLy,kDown), myThid )  #endif /* MOM_BOUNDARY_CONSERVE */
661    
662  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
663          DO j=jMin,jMax          DO j=jMin,jMax
# Line 524  C--   Tendency is minus divergence of th Line 668  C--   Tendency is minus divergence of th
668       &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )       &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
669  #else  #else
670       &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)       &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
671       &      *recip_rAs(i,j,bi,bj)       &     *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
672  #endif  #endif
673       &    *( ( fZon(i+1,j)     - fZon(i,j  ) )*uDvdxFac       &     *( ( fZon(i+1,j)  - fZon(i,j  )  )*uDvdxFac
674       &      +( fMer(i,  j)     - fMer(i,j-1) )*vDvdyFac       &       +( fMer(i,  j)  - fMer(i,j-1)  )*vDvdyFac
675       &      +(fVerV(i,j,kDown) - fVerV(i,j,kUp))*rkSign*rVelDvdrFac       &       +( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign*rVelDvdrFac
676       &     )       &     )
677           ENDDO           ENDDO
678          ENDDO          ENDDO
679    
680  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
681          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
682            CALL DIAGNOSTICS_FILL(fZon,'ADVx_Vm ',k,1,2,bi,bj,myThid)            CALL DIAGNOSTICS_FILL( fZon,  'ADVx_Vm ',k,1,2,bi,bj,myThid)
683            CALL DIAGNOSTICS_FILL(fMer,'ADVy_Vm ',k,1,2,bi,bj,myThid)            CALL DIAGNOSTICS_FILL( fMer,  'ADVy_Vm ',k,1,2,bi,bj,myThid)
684            CALL DIAGNOSTICS_FILL(fVerV(1-Olx,1-Oly,kUp),            CALL DIAGNOSTICS_FILL(fVerVkm,'ADVrE_Vm',k,1,2,bi,bj,myThid)
      &                               'ADVrE_Vm',k,1,2,bi,bj,myThid)  
685          ENDIF          ENDIF
686  #endif  #endif
687    
# Line 549  C-- account for 3.D divergence of the fl Line 692  C-- account for 3.D divergence of the fl
692           DO j=jMin,jMax           DO j=jMin,jMax
693            DO i=iMin,iMax            DO i=iMin,iMax
694             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)             gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
695       &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf       &     - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTFreeSurf
696       &       *vVel(i,j,k,bi,bj)       &       *vVel(i,j,k,bi,bj)
697            ENDDO            ENDDO
698           ENDDO           ENDDO
# Line 565  C-- account for 3.D divergence of the fl Line 708  C-- account for 3.D divergence of the fl
708  # endif /* DISABLE_RSTAR_CODE */  # endif /* DISABLE_RSTAR_CODE */
709  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
710    
711    #ifdef ALLOW_ADDFLUID
712            IF ( selectAddFluid.GE.1 ) THEN
713             DO j=jMin,jMax
714              DO i=iMin,iMax
715               gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
716         &     + vVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
717         &      *( addMass(i,j-1,k,bi,bj) + addMass(i,j,k,bi,bj) )
718         &      *_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)*recip_rhoFacC(k)
719         &      * recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
720              ENDDO
721             ENDDO
722            ENDIF
723    #endif /* ALLOW_ADDFLUID */
724    
725        ELSE        ELSE
726  C-    if momAdvection / else  C-    if momAdvection / else
727          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 579  C-    endif momAdvection. Line 736  C-    endif momAdvection.
736        IF (momViscosity) THEN        IF (momViscosity) THEN
737  C---  Calculate eddy fluxes (dissipation) between cells for meridional flow.  C---  Calculate eddy fluxes (dissipation) between cells for meridional flow.
738  C     Bi-harmonic term del^2 V -> v4F  C     Bi-harmonic term del^2 V -> v4F
739          IF (biharmonic)          IF ( useBiharmonicVisc )
740       &  CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)       &  CALL MOM_V_DEL2V( bi, bj, k, vFld, hFacZ, h0FacZ,
741         O                    v4f, myThid )
742    
743  C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon  C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
744          CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,          CALL MOM_V_XVISCFLUX( bi,bj,k,vFld,v4f,hFacZ,fZon,
745       I    viscAh_Z,viscA4_Z,myThid)       I                        viscAh_Z,viscA4_Z,myThid )
746    
747  C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer  C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
748          CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,          CALL MOM_V_YVISCFLUX( bi,bj,k,vFld,v4f,fMer,
749       I    viscAh_D,viscA4_D,myThid)       I                        viscAh_D,viscA4_D,myThid )
750    
751  C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw  C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
752         IF (.NOT.implicitViscosity) THEN         IF (.NOT.implicitViscosity) THEN
753          CALL MOM_V_RVISCFLUX(bi,bj, k, vVel,KappaRV,fVrUp,myThid)          CALL MOM_V_RVISCFLUX( bi,bj, k, vVel,KappaRV,fVrUp,myThid )
754          CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,fVrDw,myThid)          CALL MOM_V_RVISCFLUX( bi,bj,k+1,vVel,KappaRV,fVrDw,myThid )
755         ENDIF         ENDIF
756    
757  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term  C--   Tendency is minus divergence of the fluxes + coriolis + pressure term
758    C     anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
759          DO j=jMin,jMax          DO j=jMin,jMax
760           DO i=iMin,iMax           DO i=iMin,iMax
761            gvDiss(i,j) =            gvDiss(i,j) =
# Line 605  C--   Tendency is minus divergence of th Line 764  C--   Tendency is minus divergence of th
764       &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )       &      ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
765  #else  #else
766       &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)       &     -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
767       &      *recip_rAs(i,j,bi,bj)       &      *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
768  #endif  #endif
769       &    *( ( fZon(i+1,j)  - fZon(i,j  ) )*AhDvdxFac       &     *( ( fZon(i+1,j)  - fZon(i,j  ) )*AhDvdxFac
770       &      +( fMer(i,  j)  - fMer(i,j-1) )*AhDvdyFac       &       +( fMer(i,  j)  - fMer(i,j-1) )*AhDvdyFac
771       &      +( fVrDw(i,j)   - fVrUp(i,j) )*rkSign*ArDvdrFac       &       +( fVrDw(i,j)   - fVrUp(i,j) )*rkSign*ArDvdrFac
772         &                                     *recip_rhoFacC(k)
773       &     )       &     )
774           ENDDO           ENDDO
775          ENDDO          ENDDO
# Line 623  C--   Tendency is minus divergence of th Line 783  C--   Tendency is minus divergence of th
783          ENDIF          ENDIF
784  #endif  #endif
785    
786  C-- No-slip and drag BCs appear as body forces in cell abutting topography  C-- No-slip and drag BCs appear as body forces in cell abutting topography
787        IF (no_slip_sides) THEN          IF (no_slip_sides) THEN
788  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
789           CALL MOM_V_SIDEDRAG(           CALL MOM_V_SIDEDRAG( bi, bj, k,
790       I        bi,bj,k,       I        vFld, v4f, h0FacZ,
791       I        vFld, v4f, hFacZ,       I        viscAh_Z, viscA4_Z,
792       I        viscAh_Z,viscA4_Z,       I        useHarmonicVisc, useBiharmonicVisc, useVariableVisc,
      I        harmonic,biharmonic,useVariableViscosity,  
793       O        vF,       O        vF,
794       I        myThid)       I        myThid )
795           DO j=jMin,jMax           DO j=jMin,jMax
796            DO i=iMin,iMax            DO i=iMin,iMax
797             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
# Line 641  C-     No-slip BCs impose a drag at wall Line 800  C-     No-slip BCs impose a drag at wall
800          ENDIF          ENDIF
801  C-    No-slip BCs impose a drag at bottom  C-    No-slip BCs impose a drag at bottom
802          IF (bottomDragTerms) THEN          IF (bottomDragTerms) THEN
803           CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)           CALL MOM_V_BOTTOMDRAG( bi,bj,k,vFld,KE,KappaRV,vF,myThid )
804             DO j=jMin,jMax
805              DO i=iMin,iMax
806               gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
807              ENDDO
808             ENDDO
809            ENDIF
810    
811    #ifdef ALLOW_SHELFICE
812            IF (useShelfIce) THEN
813             CALL SHELFICE_V_DRAG( bi,bj,k,vFld,KE,KappaRV,vF,myThid )
814           DO j=jMin,jMax           DO j=jMin,jMax
815            DO i=iMin,iMax            DO i=iMin,iMax
816             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)             gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
817            ENDDO            ENDDO
818           ENDDO           ENDDO
819          ENDIF          ENDIF
820    #endif /* ALLOW_SHELFICE */
821    
822  C-    endif momViscosity  C-    endif momViscosity
823        ENDIF        ENDIF
# Line 660  c    I     myTime,myThid) Line 830  c    I     myTime,myThid)
830    
831  C--   Metric terms for curvilinear grid systems  C--   Metric terms for curvilinear grid systems
832        IF (useNHMTerms) THEN        IF (useNHMTerms) THEN
833  C      o Spherical polar grid metric terms  C      o Non-Hydrostatic (spherical) metric terms
834         CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)         CALL MOM_V_METRIC_NH( bi,bj,k,vFld,wVel,mT,myThid )
835         DO j=jMin,jMax         DO j=jMin,jMax
836          DO i=iMin,iMax          DO i=iMin,iMax
837           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j)
838          ENDDO          ENDDO
839         ENDDO         ENDDO
840        ENDIF        ENDIF
841        IF (usingSphericalPolarMTerms) THEN        IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
842         CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)  C      o Spherical polar grid metric terms
843           CALL MOM_V_METRIC_SPHERE( bi,bj,k,uFld,mT,myThid )
844         DO j=jMin,jMax         DO j=jMin,jMax
845          DO i=iMin,iMax          DO i=iMin,iMax
846           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
847          ENDDO          ENDDO
848         ENDDO         ENDDO
849        ENDIF        ENDIF
850        IF (usingCylindricalGrid) THEN        IF ( usingCylindricalGrid .AND. metricTerms ) THEN
851           CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)  C      o Cylindrical grid metric terms
852           DO j=jMin,jMax         CALL MOM_V_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid )
853              DO i=iMin,iMax         DO j=jMin,jMax
854                 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)          DO i=iMin,iMax
855              ENDDO           gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
856           ENDDO          ENDDO
857           ENDDO
858        ENDIF        ENDIF
859    
860  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
861    
862  C--   Coriolis term  C--   Coriolis term (call to CD_CODE_SCHEME has been moved to timestep.F)
 C     Note. As coded here, coriolis will not work with "thin walls"  
 c     IF (useCDscheme) THEN  
 c       CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)  
 c     ELSE  
863        IF (.NOT.useCDscheme) THEN        IF (.NOT.useCDscheme) THEN
864          CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)          CALL MOM_U_CORIOLIS( bi,bj,k,vFld,cf,myThid )
865          DO j=jMin,jMax          DO j=jMin,jMax
866           DO i=iMin,iMax           DO i=iMin,iMax
867            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
# Line 703  c     ELSE Line 871  c     ELSE
871          IF ( useDiagnostics )          IF ( useDiagnostics )
872       &    CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)       &    CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)
873  #endif  #endif
874          CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)          CALL MOM_V_CORIOLIS( bi,bj,k,uFld,cf,myThid )
875          DO j=jMin,jMax          DO j=jMin,jMax
876           DO i=iMin,iMax           DO i=iMin,iMax
877            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
# Line 715  c     ELSE Line 883  c     ELSE
883  #endif  #endif
884        ENDIF        ENDIF
885    
886        IF (nonHydrostatic.OR.quasiHydrostatic) THEN  C--   3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
887         CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)        IF ( use3dCoriolis ) THEN
888         DO j=jMin,jMax          CALL MOM_U_CORIOLIS_NH( bi,bj,k,wVel,cf,myThid )
889          DO i=iMin,iMax          DO j=jMin,jMax
890           gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)           DO i=iMin,iMax
891              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
892             ENDDO
893          ENDDO          ENDDO
894         ENDDO         IF ( usingCurvilinearGrid ) THEN
895    C-     presently, non zero angleSinC array only supported with Curvilinear-Grid
896            CALL MOM_V_CORIOLIS_NH( bi,bj,k,wVel,cf,myThid )
897            DO j=jMin,jMax
898             DO i=iMin,iMax
899              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
900             ENDDO
901            ENDDO
902           ENDIF
903        ENDIF        ENDIF
904    
905  C--   Set du/dt & dv/dt on boundaries to zero  C--   Set du/dt & dv/dt on boundaries to zero
# Line 737  C--   Set du/dt & dv/dt on boundaries to Line 915  C--   Set du/dt & dv/dt on boundaries to
915  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
916        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
917          CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)          CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)
918          CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),          CALL DIAGNOSTICS_FILL(gU(1-OLx,1-OLy,k,bi,bj),
919       &                               'Um_Advec',k,1,2,bi,bj,myThid)       &                               'Um_Advec',k,1,2,bi,bj,myThid)
920          CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),          CALL DIAGNOSTICS_FILL(gV(1-OLx,1-OLy,k,bi,bj),
921       &                               'Vm_Advec',k,1,2,bi,bj,myThid)       &                               'Vm_Advec',k,1,2,bi,bj,myThid)
        IF (momViscosity) THEN  
         CALL DIAGNOSTICS_FILL(guDiss,'Um_Diss ',k,1,2,bi,bj,myThid)  
         CALL DIAGNOSTICS_FILL(gvDiss,'Vm_Diss ',k,1,2,bi,bj,myThid)  
        ENDIF  
922        ENDIF        ENDIF
923  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
924    
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22