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

Annotation of /MITgcm/pkg/mom_fluxform/mom_fluxform.F

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


Revision 1.32 - (hide annotations) (download)
Tue Feb 7 11:46:18 2006 UTC (18 years, 3 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58b_post, checkpoint58a_post
Changes since 1.31: +24 -2 lines
o add hooks for friction at water-shelfice interface

1 mlosch 1.32 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.31 2005/12/08 15:44:34 heimbach Exp $
2 adcroft 1.2 C $Name: $
3 adcroft 1.1
4 adcroft 1.3 CBOI
5     C !TITLE: pkg/mom\_advdiff
6     C !AUTHORS: adcroft@mit.edu
7 adcroft 1.4 C !INTRODUCTION: Flux-form Momentum Equations Package
8 adcroft 1.3 C
9     C Package "mom\_fluxform" provides methods for calculating explicit terms
10     C in the momentum equation cast in flux-form:
11     C \begin{eqnarray*}
12     C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
13     C -\nabla \cdot {\bf v} u
14     C -fv
15     C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
16     C + \mbox{metrics}
17     C \\
18     C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
19     C -\nabla \cdot {\bf v} v
20     C +fu
21     C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
22     C + \mbox{metrics}
23     C \end{eqnarray*}
24     C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
25     C stresses as well as internal viscous stresses.
26     CEOI
27    
28 edhill 1.13 #include "MOM_FLUXFORM_OPTIONS.h"
29 adcroft 1.1
30 adcroft 1.3 CBOP
31     C !ROUTINE: MOM_FLUXFORM
32    
33     C !INTERFACE: ==========================================================
34 adcroft 1.1 SUBROUTINE MOM_FLUXFORM(
35     I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
36 jmc 1.23 I KappaRU, KappaRV,
37 adcroft 1.1 U fVerU, fVerV,
38 jmc 1.23 O guDiss, gvDiss,
39     I myTime, myIter, myThid)
40 adcroft 1.3
41     C !DESCRIPTION:
42     C Calculates all the horizontal accelerations except for the implicit surface
43     C pressure gradient and implciit vertical viscosity.
44 adcroft 1.1
45 adcroft 1.3 C !USES: ===============================================================
46 adcroft 1.1 C == Global variables ==
47 adcroft 1.3 IMPLICIT NONE
48 adcroft 1.1 #include "SIZE.h"
49     #include "DYNVARS.h"
50     #include "FFIELDS.h"
51     #include "EEPARAMS.h"
52     #include "PARAMS.h"
53     #include "GRID.h"
54     #include "SURFACE.h"
55    
56 adcroft 1.3 C !INPUT PARAMETERS: ===================================================
57     C bi,bj :: tile indices
58     C iMin,iMax,jMin,jMAx :: loop ranges
59     C k :: vertical level
60     C kUp :: =1 or 2 for consecutive k
61     C kDown :: =2 or 1 for consecutive k
62     C KappaRU :: vertical viscosity
63     C KappaRV :: vertical viscosity
64     C fVerU :: vertical flux of U, 2 1/2 dim for pipe-lining
65     C fVerV :: vertical flux of V, 2 1/2 dim for pipe-lining
66 jmc 1.23 C guDiss :: dissipation tendency (all explicit terms), u component
67     C gvDiss :: dissipation tendency (all explicit terms), v component
68 jmc 1.8 C myTime :: current time
69 adcroft 1.3 C myIter :: current time-step number
70     C myThid :: thread number
71     INTEGER bi,bj,iMin,iMax,jMin,jMax
72     INTEGER k,kUp,kDown
73 adcroft 1.1 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74     _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75     _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
76     _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
77 jmc 1.23 _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78     _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 jmc 1.8 _RL myTime
80 adcroft 1.2 INTEGER myIter
81 adcroft 1.1 INTEGER myThid
82    
83 adcroft 1.3 C !OUTPUT PARAMETERS: ==================================================
84     C None - updates gU() and gV() in common blocks
85    
86     C !LOCAL VARIABLES: ====================================================
87     C i,j :: loop indices
88     C vF :: viscous flux
89     C v4F :: bi-harmonic viscous flux
90     C cF :: Coriolis acceleration
91     C mT :: Metric terms
92     C fZon :: zonal fluxes
93     C fMer :: meridional fluxes
94 jmc 1.23 C fVrUp,fVrDw :: vertical viscous fluxes at interface k-1 & k
95 adcroft 1.3 INTEGER i,j
96     _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97     _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101     _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 jmc 1.23 _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103     _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 adcroft 1.1 C afFacMom - Tracer parameters for turning terms
105     C vfFacMom on and off.
106     C pfFacMom afFacMom - Advective terms
107     C cfFacMom vfFacMom - Eddy viscosity terms
108     C mTFacMom pfFacMom - Pressure terms
109     C cfFacMom - Coriolis terms
110     C foFacMom - Forcing
111     C mTFacMom - Metric term
112 jmc 1.23 C uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off
113 adcroft 1.1 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114     _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115     _RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116     _RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119     _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120     _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 jmc 1.8 _RL rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122     _RL rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 adcroft 1.18 _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 baylor 1.25 _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125     _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126     _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127     _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128     _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129     _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130 adcroft 1.18 _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131     _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 adcroft 1.1 _RL uDudxFac
133     _RL AhDudxFac
134     _RL vDudyFac
135     _RL AhDudyFac
136     _RL rVelDudrFac
137     _RL ArDudrFac
138     _RL fuFac
139     _RL mtFacU
140     _RL uDvdxFac
141     _RL AhDvdxFac
142     _RL vDvdyFac
143     _RL AhDvdyFac
144     _RL rVelDvdrFac
145     _RL ArDvdrFac
146     _RL fvFac
147     _RL mtFacV
148 jmc 1.29 _RL sideMaskFac
149 baylor 1.25 LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity
150 adcroft 1.3 CEOP
151 adcroft 1.1
152     C Initialise intermediate terms
153 jmc 1.23 DO j=1-OLy,sNy+OLy
154     DO i=1-OLx,sNx+OLx
155 adcroft 1.1 vF(i,j) = 0.
156     v4F(i,j) = 0.
157     cF(i,j) = 0.
158     mT(i,j) = 0.
159     fZon(i,j) = 0.
160     fMer(i,j) = 0.
161 jmc 1.23 fVrUp(i,j)= 0.
162     fVrDw(i,j)= 0.
163     rTransU(i,j)= 0.
164     rTransV(i,j)= 0.
165 adcroft 1.18 strain(i,j) = 0.
166 jmc 1.23 tension(i,j)= 0.
167     guDiss(i,j) = 0.
168     gvDiss(i,j) = 0.
169 heimbach 1.30 #ifdef ALLOW_AUTODIFF_TAMC
170     vort3(i,j) = 0. _d 0
171     strain(i,j) = 0. _d 0
172     tension(i,j) = 0. _d 0
173     #endif
174 adcroft 1.1 ENDDO
175     ENDDO
176    
177     C-- Term by term tracer parmeters
178     C o U momentum equation
179     uDudxFac = afFacMom*1.
180     AhDudxFac = vfFacMom*1.
181     vDudyFac = afFacMom*1.
182     AhDudyFac = vfFacMom*1.
183     rVelDudrFac = afFacMom*1.
184     ArDudrFac = vfFacMom*1.
185     mTFacU = mtFacMom*1.
186     fuFac = cfFacMom*1.
187     C o V momentum equation
188     uDvdxFac = afFacMom*1.
189     AhDvdxFac = vfFacMom*1.
190     vDvdyFac = afFacMom*1.
191     AhDvdyFac = vfFacMom*1.
192     rVelDvdrFac = afFacMom*1.
193     ArDvdrFac = vfFacMom*1.
194     mTFacV = mtFacMom*1.
195     fvFac = cfFacMom*1.
196 jmc 1.23
197     IF (implicitViscosity) THEN
198     ArDudrFac = 0.
199     ArDvdrFac = 0.
200     ENDIF
201 adcroft 1.1
202 jmc 1.29 C note: using standard stencil (no mask) results in under-estimating
203     C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
204     IF ( no_slip_sides ) THEN
205     sideMaskFac = sideDragFactor
206     ELSE
207     sideMaskFac = 0. _d 0
208     ENDIF
209    
210 adcroft 1.1 IF ( no_slip_bottom
211     & .OR. bottomDragQuadratic.NE.0.
212     & .OR. bottomDragLinear.NE.0.) THEN
213     bottomDragTerms=.TRUE.
214     ELSE
215     bottomDragTerms=.FALSE.
216     ENDIF
217    
218     C-- Calculate open water fraction at vorticity points
219     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
220    
221     C---- Calculate common quantities used in both U and V equations
222     C Calculate tracer cell face open areas
223     DO j=1-OLy,sNy+OLy
224     DO i=1-OLx,sNx+OLx
225     xA(i,j) = _dyG(i,j,bi,bj)
226     & *drF(k)*_hFacW(i,j,k,bi,bj)
227     yA(i,j) = _dxG(i,j,bi,bj)
228     & *drF(k)*_hFacS(i,j,k,bi,bj)
229     ENDDO
230     ENDDO
231    
232     C Make local copies of horizontal flow field
233     DO j=1-OLy,sNy+OLy
234     DO i=1-OLx,sNx+OLx
235     uFld(i,j) = uVel(i,j,k,bi,bj)
236     vFld(i,j) = vVel(i,j,k,bi,bj)
237     ENDDO
238     ENDDO
239    
240     C Calculate velocity field "volume transports" through tracer cell faces.
241     DO j=1-OLy,sNy+OLy
242     DO i=1-OLx,sNx+OLx
243     uTrans(i,j) = uFld(i,j)*xA(i,j)
244     vTrans(i,j) = vFld(i,j)*yA(i,j)
245     ENDDO
246     ENDDO
247    
248 baylor 1.25 CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)
249 jmc 1.29 IF ( momViscosity) THEN
250     CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
251     CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
252     CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
253     CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
254     DO j=1-Oly,sNy+Oly
255     DO i=1-Olx,sNx+Olx
256     IF ( hFacZ(i,j).EQ.0. ) THEN
257     vort3(i,j) = sideMaskFac*vort3(i,j)
258     strain(i,j) = sideMaskFac*strain(i,j)
259     ENDIF
260     ENDDO
261     ENDDO
262     #ifdef ALLOW_DIAGNOSTICS
263     IF ( useDiagnostics ) THEN
264     CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
265     CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
266     CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
267     CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
268     ENDIF
269     #endif
270     ENDIF
271 adcroft 1.18
272 jmc 1.8 C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
273     IF (momAdvection.AND.k.EQ.1) THEN
274    
275     C- Calculate vertical transports above U & V points (West & South face):
276 jmc 1.23 CALL MOM_CALC_RTRANS( k, bi, bj,
277     O rTransU, rTransV,
278     I myTime, myIter, myThid)
279 jmc 1.8
280     C- Free surface correction term (flux at k=1)
281 jmc 1.23 CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
282     O fVerU(1-OLx,1-OLy,kUp), myThid )
283 jmc 1.8
284 jmc 1.23 CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
285     O fVerV(1-OLx,1-OLy,kUp), myThid )
286 jmc 1.8
287     C--- endif momAdvection & k=1
288     ENDIF
289    
290    
291     C--- Calculate vertical transports (at k+1) below U & V points :
292     IF (momAdvection) THEN
293 jmc 1.23 CALL MOM_CALC_RTRANS( k+1, bi, bj,
294     O rTransU, rTransV,
295     I myTime, myIter, myThid)
296 jmc 1.8 ENDIF
297    
298 baylor 1.25 IF (momViscosity) THEN
299     CALL MOM_CALC_VISC(
300     I bi,bj,k,
301     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
302     O harmonic,biharmonic,useVariableViscosity,
303 jmc 1.26 I hDiv,vort3,tension,strain,KE,hFacZ,
304 baylor 1.25 I myThid)
305     ENDIF
306 jmc 1.8
307 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
308    
309 adcroft 1.1 C---- Zonal momentum equation starts here
310    
311 jmc 1.23 IF (momAdvection) THEN
312     C--- Calculate mean fluxes (advection) between cells for zonal flow.
313 adcroft 1.1
314     C-- Zonal flux (fZon is at east face of "u" cell)
315 jmc 1.23 C Mean flow component of zonal flux -> fZon
316     CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,fZon,myThid)
317 adcroft 1.1
318     C-- Meridional flux (fMer is at south face of "u" cell)
319 jmc 1.23 C Mean flow component of meridional flux -> fMer
320     CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,fMer,myThid)
321 adcroft 1.1
322     C-- Vertical flux (fVer is at upper face of "u" cell)
323 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVer
324     CALL MOM_U_ADV_WU(
325     I bi,bj,k+1,uVel,wVel,rTransU,
326     O fVerU(1-OLx,1-OLy,kDown), myThid )
327 adcroft 1.1
328     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
329 jmc 1.23 DO j=jMin,jMax
330     DO i=iMin,iMax
331     gU(i,j,k,bi,bj) =
332 adcroft 1.1 #ifdef OLD_UV_GEOM
333 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
334     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
335 adcroft 1.1 #else
336 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
337     & *recip_rAw(i,j,bi,bj)
338 adcroft 1.1 #endif
339 jmc 1.23 & *( ( fZon(i,j ) - fZon(i-1,j) )*uDudxFac
340     & +( fMer(i,j+1) - fMer(i, j) )*vDudyFac
341     & +(fVerU(i,j,kDown) - fVerU(i,j,kUp))*rkSign*rVelDudrFac
342     & )
343     ENDDO
344     ENDDO
345 adcroft 1.1
346 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
347     IF ( useDiagnostics ) THEN
348     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Um ',k,1,2,bi,bj,myThid)
349     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Um ',k,1,2,bi,bj,myThid)
350     CALL DIAGNOSTICS_FILL(fVerU(1-Olx,1-Oly,kUp),
351     & 'ADVrE_Um',k,1,2,bi,bj,myThid)
352     ENDIF
353     #endif
354    
355 jmc 1.8 #ifdef NONLIN_FRSURF
356     C-- account for 3.D divergence of the flow in rStar coordinate:
357 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
358 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
359     DO j=jMin,jMax
360     DO i=iMin,iMax
361     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
362 jmc 1.8 & - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
363     & *uVel(i,j,k,bi,bj)
364 jmc 1.23 ENDDO
365     ENDDO
366     ENDIF
367     IF ( select_rStar.LT.0 ) THEN
368     DO j=jMin,jMax
369     DO i=iMin,iMax
370     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
371     & - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
372     ENDDO
373     ENDDO
374     ENDIF
375 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
376 jmc 1.23 #endif /* NONLIN_FRSURF */
377    
378     ELSE
379     C- if momAdvection / else
380     DO j=1-OLy,sNy+OLy
381     DO i=1-OLx,sNx+OLx
382     gU(i,j,k,bi,bj) = 0. _d 0
383     ENDDO
384 jmc 1.8 ENDDO
385 jmc 1.23
386     C- endif momAdvection.
387 jmc 1.8 ENDIF
388 jmc 1.23
389     IF (momViscosity) THEN
390     C--- Calculate eddy fluxes (dissipation) between cells for zonal flow.
391    
392     C Bi-harmonic term del^2 U -> v4F
393 baylor 1.25 IF (biharmonic)
394 jmc 1.23 & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
395    
396     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
397 baylor 1.25 CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,
398 baylor 1.27 I viscAh_D,viscA4_D,myThid)
399 jmc 1.23
400     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
401 baylor 1.25 CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,
402 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
403 jmc 1.23
404     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
405     IF (.NOT.implicitViscosity) THEN
406     CALL MOM_U_RVISCFLUX(bi,bj, k, uVel,KappaRU,fVrUp,myThid)
407     CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,fVrDw,myThid)
408     ENDIF
409    
410     C-- Tendency is minus divergence of the fluxes
411     DO j=jMin,jMax
412     DO i=iMin,iMax
413     guDiss(i,j) =
414     #ifdef OLD_UV_GEOM
415     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
416     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
417     #else
418     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
419     & *recip_rAw(i,j,bi,bj)
420     #endif
421     & *( ( fZon(i,j ) - fZon(i-1,j) )*AhDudxFac
422     & +( fMer(i,j+1) - fMer(i, j) )*AhDudyFac
423     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDudrFac
424     & )
425     ENDDO
426 jmc 1.8 ENDDO
427    
428 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
429     IF ( useDiagnostics ) THEN
430     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
431     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
432     IF (.NOT.implicitViscosity)
433     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
434     ENDIF
435     #endif
436    
437 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
438 jmc 1.23 IF (no_slip_sides) THEN
439 adcroft 1.1 C- No-slip BCs impose a drag at walls...
440 baylor 1.27 CALL MOM_U_SIDEDRAG(
441     I bi,bj,k,
442     I uFld, v4f, hFacZ,
443     I viscAh_Z,viscA4_Z,
444     I harmonic,biharmonic,useVariableViscosity,
445     O vF,
446     I myThid)
447 jmc 1.23 DO j=jMin,jMax
448     DO i=iMin,iMax
449     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
450     ENDDO
451     ENDDO
452     ENDIF
453 adcroft 1.1 C- No-slip BCs impose a drag at bottom
454 jmc 1.23 IF (bottomDragTerms) THEN
455     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
456     DO j=jMin,jMax
457     DO i=iMin,iMax
458     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
459     ENDDO
460     ENDDO
461     ENDIF
462    
463 mlosch 1.32 #ifdef ALLOW_SHELFICE
464     IF (useShelfIce) THEN
465     CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
466     DO j=jMin,jMax
467     DO i=iMin,iMax
468     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
469     ENDDO
470     ENDDO
471     ENDIF
472     #endif /* ALLOW_SHELFICE */
473    
474 jmc 1.23 C- endif momViscosity
475 adcroft 1.1 ENDIF
476    
477 jmc 1.12 C-- Forcing term (moved to timestep.F)
478     c IF (momForcing)
479     c & CALL EXTERNAL_FORCING_U(
480     c I iMin,iMax,jMin,jMax,bi,bj,k,
481     c I myTime,myThid)
482 adcroft 1.1
483     C-- Metric terms for curvilinear grid systems
484 adcroft 1.5 IF (useNHMTerms) THEN
485     C o Non-hydrosatic metric terms
486 adcroft 1.1 CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
487     DO j=jMin,jMax
488     DO i=iMin,iMax
489     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
490     ENDDO
491     ENDDO
492 adcroft 1.5 ENDIF
493     IF (usingSphericalPolarMTerms) THEN
494 adcroft 1.1 CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
495     DO j=jMin,jMax
496     DO i=iMin,iMax
497     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
498     ENDDO
499     ENDDO
500 afe 1.20 ENDIF
501 afe 1.19 IF (usingCylindricalGrid) THEN
502     CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
503     DO j=jMin,jMax
504     DO i=iMin,iMax
505     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
506     ENDDO
507     ENDDO
508 adcroft 1.1 ENDIF
509    
510 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
511 adcroft 1.1
512     C---- Meridional momentum equation starts here
513    
514 jmc 1.23 IF (momAdvection) THEN
515     C--- Calculate mean fluxes (advection) between cells for meridional flow.
516     C Mean flow component of zonal flux -> fZon
517     CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,fZon,myThid)
518 adcroft 1.1
519     C-- Meridional flux (fMer is at north face of "v" cell)
520 jmc 1.23 C Mean flow component of meridional flux -> fMer
521     CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,fMer,myThid)
522 adcroft 1.1
523     C-- Vertical flux (fVer is at upper face of "v" cell)
524 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVerV
525     CALL MOM_V_ADV_WV(
526     I bi,bj,k+1,vVel,wVel,rTransV,
527     O fVerV(1-OLx,1-OLy,kDown), myThid )
528 adcroft 1.1
529     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
530 jmc 1.23 DO j=jMin,jMax
531     DO i=iMin,iMax
532     gV(i,j,k,bi,bj) =
533 adcroft 1.1 #ifdef OLD_UV_GEOM
534 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
535     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
536 adcroft 1.1 #else
537 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
538     & *recip_rAs(i,j,bi,bj)
539 adcroft 1.1 #endif
540 jmc 1.23 & *( ( fZon(i+1,j) - fZon(i,j ) )*uDvdxFac
541     & +( fMer(i, j) - fMer(i,j-1) )*vDvdyFac
542     & +(fVerV(i,j,kDown) - fVerV(i,j,kUp))*rkSign*rVelDvdrFac
543     & )
544 jmc 1.24 ENDDO
545     ENDDO
546    
547     #ifdef ALLOW_DIAGNOSTICS
548     IF ( useDiagnostics ) THEN
549     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Vm ',k,1,2,bi,bj,myThid)
550     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Vm ',k,1,2,bi,bj,myThid)
551     CALL DIAGNOSTICS_FILL(fVerV(1-Olx,1-Oly,kUp),
552     & 'ADVrE_Vm',k,1,2,bi,bj,myThid)
553     ENDIF
554     #endif
555 adcroft 1.1
556 jmc 1.8 #ifdef NONLIN_FRSURF
557     C-- account for 3.D divergence of the flow in rStar coordinate:
558 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
559 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
560     DO j=jMin,jMax
561     DO i=iMin,iMax
562     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
563 jmc 1.8 & - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
564     & *vVel(i,j,k,bi,bj)
565 jmc 1.23 ENDDO
566     ENDDO
567     ENDIF
568     IF ( select_rStar.LT.0 ) THEN
569     DO j=jMin,jMax
570     DO i=iMin,iMax
571     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
572     & - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
573     ENDDO
574     ENDDO
575     ENDIF
576 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
577 jmc 1.23 #endif /* NONLIN_FRSURF */
578    
579     ELSE
580     C- if momAdvection / else
581     DO j=1-OLy,sNy+OLy
582     DO i=1-OLx,sNx+OLx
583     gV(i,j,k,bi,bj) = 0. _d 0
584     ENDDO
585 jmc 1.8 ENDDO
586 jmc 1.23
587     C- endif momAdvection.
588 jmc 1.8 ENDIF
589 jmc 1.23
590     IF (momViscosity) THEN
591     C--- Calculate eddy fluxes (dissipation) between cells for meridional flow.
592     C Bi-harmonic term del^2 V -> v4F
593 baylor 1.25 IF (biharmonic)
594 jmc 1.23 & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
595    
596     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
597 baylor 1.25 CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,
598 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
599 jmc 1.23
600     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
601 baylor 1.25 CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,
602 baylor 1.27 I viscAh_D,viscA4_D,myThid)
603 jmc 1.23
604     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
605     IF (.NOT.implicitViscosity) THEN
606     CALL MOM_V_RVISCFLUX(bi,bj, k, vVel,KappaRV,fVrUp,myThid)
607     CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,fVrDw,myThid)
608     ENDIF
609    
610     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
611     DO j=jMin,jMax
612     DO i=iMin,iMax
613     gvDiss(i,j) =
614     #ifdef OLD_UV_GEOM
615     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
616     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
617     #else
618     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
619     & *recip_rAs(i,j,bi,bj)
620     #endif
621     & *( ( fZon(i+1,j) - fZon(i,j ) )*AhDvdxFac
622     & +( fMer(i, j) - fMer(i,j-1) )*AhDvdyFac
623     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDvdrFac
624     & )
625     ENDDO
626 jmc 1.8 ENDDO
627    
628 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
629     IF ( useDiagnostics ) THEN
630     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
631     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
632     IF (.NOT.implicitViscosity)
633     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
634     ENDIF
635     #endif
636    
637 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
638 mlosch 1.32 IF (no_slip_sides) THEN
639 adcroft 1.1 C- No-slip BCs impose a drag at walls...
640 baylor 1.27 CALL MOM_V_SIDEDRAG(
641     I bi,bj,k,
642     I vFld, v4f, hFacZ,
643     I viscAh_Z,viscA4_Z,
644     I harmonic,biharmonic,useVariableViscosity,
645     O vF,
646     I myThid)
647 jmc 1.23 DO j=jMin,jMax
648     DO i=iMin,iMax
649     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
650     ENDDO
651     ENDDO
652     ENDIF
653 adcroft 1.1 C- No-slip BCs impose a drag at bottom
654 jmc 1.23 IF (bottomDragTerms) THEN
655     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
656     DO j=jMin,jMax
657     DO i=iMin,iMax
658     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
659     ENDDO
660     ENDDO
661     ENDIF
662    
663 mlosch 1.32 #ifdef ALLOW_SHELFICE
664     IF (useShelfIce) THEN
665     CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
666     DO j=jMin,jMax
667     DO i=iMin,iMax
668     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
669     ENDDO
670     ENDDO
671     ENDIF
672     #endif /* ALLOW_SHELFICE */
673    
674 jmc 1.23 C- endif momViscosity
675 adcroft 1.1 ENDIF
676    
677 jmc 1.12 C-- Forcing term (moved to timestep.F)
678     c IF (momForcing)
679     c & CALL EXTERNAL_FORCING_V(
680     c I iMin,iMax,jMin,jMax,bi,bj,k,
681     c I myTime,myThid)
682 adcroft 1.1
683     C-- Metric terms for curvilinear grid systems
684 adcroft 1.5 IF (useNHMTerms) THEN
685 adcroft 1.1 C o Spherical polar grid metric terms
686     CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
687     DO j=jMin,jMax
688     DO i=iMin,iMax
689     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
690     ENDDO
691     ENDDO
692 adcroft 1.5 ENDIF
693     IF (usingSphericalPolarMTerms) THEN
694 adcroft 1.1 CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
695     DO j=jMin,jMax
696     DO i=iMin,iMax
697     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
698     ENDDO
699     ENDDO
700     ENDIF
701 afe 1.19 IF (usingCylindricalGrid) THEN
702     CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
703     DO j=jMin,jMax
704     DO i=iMin,iMax
705     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
706     ENDDO
707     ENDDO
708     ENDIF
709 adcroft 1.1
710 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
711 adcroft 1.1
712     C-- Coriolis term
713     C Note. As coded here, coriolis will not work with "thin walls"
714 jmc 1.12 c IF (useCDscheme) THEN
715     c CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
716     c ELSE
717     IF (.NOT.useCDscheme) THEN
718     CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
719     DO j=jMin,jMax
720     DO i=iMin,iMax
721     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
722     ENDDO
723     ENDDO
724 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
725     IF ( useDiagnostics )
726     & CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)
727     #endif
728 jmc 1.12 CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
729     DO j=jMin,jMax
730     DO i=iMin,iMax
731     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
732     ENDDO
733     ENDDO
734 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
735     IF ( useDiagnostics )
736     & CALL DIAGNOSTICS_FILL(cf,'Vm_Cori ',k,1,2,bi,bj,myThid)
737     #endif
738 jmc 1.12 ENDIF
739    
740 adcroft 1.7 IF (nonHydrostatic.OR.quasiHydrostatic) THEN
741 adcroft 1.6 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
742     DO j=jMin,jMax
743     DO i=iMin,iMax
744     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
745     ENDDO
746     ENDDO
747     ENDIF
748 adcroft 1.1
749 jmc 1.23 C-- Set du/dt & dv/dt on boundaries to zero
750     DO j=jMin,jMax
751     DO i=iMin,iMax
752     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
753     guDiss(i,j) = guDiss(i,j) *_maskW(i,j,k,bi,bj)
754     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
755     gvDiss(i,j) = gvDiss(i,j) *_maskS(i,j,k,bi,bj)
756     ENDDO
757     ENDDO
758    
759 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
760     IF ( useDiagnostics ) THEN
761 baylor 1.28 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
762 jmc 1.24 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
763     & 'Um_Advec',k,1,2,bi,bj,myThid)
764     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
765     & 'Vm_Advec',k,1,2,bi,bj,myThid)
766     IF (momViscosity) THEN
767     CALL DIAGNOSTICS_FILL(guDiss,'Um_Diss ',k,1,2,bi,bj,myThid)
768     CALL DIAGNOSTICS_FILL(gvDiss,'Vm_Diss ',k,1,2,bi,bj,myThid)
769     ENDIF
770     ENDIF
771     #endif /* ALLOW_DIAGNOSTICS */
772    
773 adcroft 1.1 RETURN
774     END

  ViewVC Help
Powered by ViewVC 1.1.22