/[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.26 - (hide annotations) (download)
Fri Sep 23 15:19:38 2005 UTC (18 years, 8 months ago) by jmc
Branch: MAIN
Changes since 1.25: +2 -2 lines 
add missing argument in CALL MOM_CALC_VISC
1 jmc 1.26 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.25 2005/09/16 20:57:09 baylor 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 baylor 1.25 LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity
149 adcroft 1.3 CEOP
150 adcroft 1.1
151     C Initialise intermediate terms
152 jmc 1.23 DO j=1-OLy,sNy+OLy
153     DO i=1-OLx,sNx+OLx
154 adcroft 1.1 vF(i,j) = 0.
155     v4F(i,j) = 0.
156     cF(i,j) = 0.
157     mT(i,j) = 0.
158     fZon(i,j) = 0.
159     fMer(i,j) = 0.
160 jmc 1.23 fVrUp(i,j)= 0.
161     fVrDw(i,j)= 0.
162     rTransU(i,j)= 0.
163     rTransV(i,j)= 0.
164 adcroft 1.18 strain(i,j) = 0.
165 jmc 1.23 tension(i,j)= 0.
166     guDiss(i,j) = 0.
167     gvDiss(i,j) = 0.
168 adcroft 1.1 ENDDO
169     ENDDO
170    
171     C-- Term by term tracer parmeters
172     C o U momentum equation
173     uDudxFac = afFacMom*1.
174     AhDudxFac = vfFacMom*1.
175     vDudyFac = afFacMom*1.
176     AhDudyFac = vfFacMom*1.
177     rVelDudrFac = afFacMom*1.
178     ArDudrFac = vfFacMom*1.
179     mTFacU = mtFacMom*1.
180     fuFac = cfFacMom*1.
181     C o V momentum equation
182     uDvdxFac = afFacMom*1.
183     AhDvdxFac = vfFacMom*1.
184     vDvdyFac = afFacMom*1.
185     AhDvdyFac = vfFacMom*1.
186     rVelDvdrFac = afFacMom*1.
187     ArDvdrFac = vfFacMom*1.
188     mTFacV = mtFacMom*1.
189     fvFac = cfFacMom*1.
190 jmc 1.23
191     IF (implicitViscosity) THEN
192     ArDudrFac = 0.
193     ArDvdrFac = 0.
194     ENDIF
195 adcroft 1.1
196     IF ( no_slip_bottom
197     & .OR. bottomDragQuadratic.NE.0.
198     & .OR. bottomDragLinear.NE.0.) THEN
199     bottomDragTerms=.TRUE.
200     ELSE
201     bottomDragTerms=.FALSE.
202     ENDIF
203    
204     C-- Calculate open water fraction at vorticity points
205     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
206    
207     C---- Calculate common quantities used in both U and V equations
208     C Calculate tracer cell face open areas
209     DO j=1-OLy,sNy+OLy
210     DO i=1-OLx,sNx+OLx
211     xA(i,j) = _dyG(i,j,bi,bj)
212     & *drF(k)*_hFacW(i,j,k,bi,bj)
213     yA(i,j) = _dxG(i,j,bi,bj)
214     & *drF(k)*_hFacS(i,j,k,bi,bj)
215     ENDDO
216     ENDDO
217    
218     C Make local copies of horizontal flow field
219     DO j=1-OLy,sNy+OLy
220     DO i=1-OLx,sNx+OLx
221     uFld(i,j) = uVel(i,j,k,bi,bj)
222     vFld(i,j) = vVel(i,j,k,bi,bj)
223     ENDDO
224     ENDDO
225    
226     C Calculate velocity field "volume transports" through tracer cell faces.
227     DO j=1-OLy,sNy+OLy
228     DO i=1-OLx,sNx+OLx
229     uTrans(i,j) = uFld(i,j)*xA(i,j)
230     vTrans(i,j) = vFld(i,j)*yA(i,j)
231     ENDDO
232     ENDDO
233    
234 baylor 1.25 CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)
235     CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
236     CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
237     CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
238     CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
239 adcroft 1.18
240 jmc 1.8 C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
241     IF (momAdvection.AND.k.EQ.1) THEN
242    
243     C- Calculate vertical transports above U & V points (West & South face):
244 jmc 1.23 CALL MOM_CALC_RTRANS( k, bi, bj,
245     O rTransU, rTransV,
246     I myTime, myIter, myThid)
247 jmc 1.8
248     C- Free surface correction term (flux at k=1)
249 jmc 1.23 CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
250     O fVerU(1-OLx,1-OLy,kUp), myThid )
251 jmc 1.8
252 jmc 1.23 CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
253     O fVerV(1-OLx,1-OLy,kUp), myThid )
254 jmc 1.8
255     C--- endif momAdvection & k=1
256     ENDIF
257    
258    
259     C--- Calculate vertical transports (at k+1) below U & V points :
260     IF (momAdvection) THEN
261 jmc 1.23 CALL MOM_CALC_RTRANS( k+1, bi, bj,
262     O rTransU, rTransV,
263     I myTime, myIter, myThid)
264 jmc 1.8 ENDIF
265    
266 baylor 1.25 IF (momViscosity) THEN
267     CALL MOM_CALC_VISC(
268     I bi,bj,k,
269     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
270     O harmonic,biharmonic,useVariableViscosity,
271 jmc 1.26 I hDiv,vort3,tension,strain,KE,hFacZ,
272 baylor 1.25 I myThid)
273     ENDIF
274 jmc 1.8
275 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
276    
277 adcroft 1.1 C---- Zonal momentum equation starts here
278    
279 jmc 1.23 IF (momAdvection) THEN
280     C--- Calculate mean fluxes (advection) between cells for zonal flow.
281 adcroft 1.1
282     C-- Zonal flux (fZon is at east face of "u" cell)
283 jmc 1.23 C Mean flow component of zonal flux -> fZon
284     CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,fZon,myThid)
285 adcroft 1.1
286     C-- Meridional flux (fMer is at south face of "u" cell)
287 jmc 1.23 C Mean flow component of meridional flux -> fMer
288     CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,fMer,myThid)
289 adcroft 1.1
290     C-- Vertical flux (fVer is at upper face of "u" cell)
291 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVer
292     CALL MOM_U_ADV_WU(
293     I bi,bj,k+1,uVel,wVel,rTransU,
294     O fVerU(1-OLx,1-OLy,kDown), myThid )
295 adcroft 1.1
296     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
297 jmc 1.23 DO j=jMin,jMax
298     DO i=iMin,iMax
299     gU(i,j,k,bi,bj) =
300 adcroft 1.1 #ifdef OLD_UV_GEOM
301 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
302     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
303 adcroft 1.1 #else
304 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
305     & *recip_rAw(i,j,bi,bj)
306 adcroft 1.1 #endif
307 jmc 1.23 & *( ( fZon(i,j ) - fZon(i-1,j) )*uDudxFac
308     & +( fMer(i,j+1) - fMer(i, j) )*vDudyFac
309     & +(fVerU(i,j,kDown) - fVerU(i,j,kUp))*rkSign*rVelDudrFac
310     & )
311     ENDDO
312     ENDDO
313 adcroft 1.1
314 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
315     IF ( useDiagnostics ) THEN
316     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Um ',k,1,2,bi,bj,myThid)
317     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Um ',k,1,2,bi,bj,myThid)
318     CALL DIAGNOSTICS_FILL(fVerU(1-Olx,1-Oly,kUp),
319     & 'ADVrE_Um',k,1,2,bi,bj,myThid)
320     ENDIF
321     #endif
322    
323 jmc 1.8 #ifdef NONLIN_FRSURF
324     C-- account for 3.D divergence of the flow in rStar coordinate:
325 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
326     DO j=jMin,jMax
327     DO i=iMin,iMax
328     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
329 jmc 1.8 & - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
330     & *uVel(i,j,k,bi,bj)
331 jmc 1.23 ENDDO
332     ENDDO
333     ENDIF
334     IF ( select_rStar.LT.0 ) THEN
335     DO j=jMin,jMax
336     DO i=iMin,iMax
337     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
338     & - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
339     ENDDO
340     ENDDO
341     ENDIF
342     #endif /* NONLIN_FRSURF */
343    
344     ELSE
345     C- if momAdvection / else
346     DO j=1-OLy,sNy+OLy
347     DO i=1-OLx,sNx+OLx
348     gU(i,j,k,bi,bj) = 0. _d 0
349     ENDDO
350 jmc 1.8 ENDDO
351 jmc 1.23
352     C- endif momAdvection.
353 jmc 1.8 ENDIF
354 jmc 1.23
355     IF (momViscosity) THEN
356     C--- Calculate eddy fluxes (dissipation) between cells for zonal flow.
357    
358     C Bi-harmonic term del^2 U -> v4F
359 baylor 1.25 IF (biharmonic)
360 jmc 1.23 & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
361    
362     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
363 baylor 1.25 CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,
364     I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
365 jmc 1.23
366     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
367 baylor 1.25 CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,
368     I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
369 jmc 1.23
370     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
371     IF (.NOT.implicitViscosity) THEN
372     CALL MOM_U_RVISCFLUX(bi,bj, k, uVel,KappaRU,fVrUp,myThid)
373     CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,fVrDw,myThid)
374     ENDIF
375    
376     C-- Tendency is minus divergence of the fluxes
377     DO j=jMin,jMax
378     DO i=iMin,iMax
379     guDiss(i,j) =
380     #ifdef OLD_UV_GEOM
381     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
382     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
383     #else
384     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
385     & *recip_rAw(i,j,bi,bj)
386     #endif
387     & *( ( fZon(i,j ) - fZon(i-1,j) )*AhDudxFac
388     & +( fMer(i,j+1) - fMer(i, j) )*AhDudyFac
389     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDudrFac
390     & )
391     ENDDO
392 jmc 1.8 ENDDO
393    
394 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
395     IF ( useDiagnostics ) THEN
396     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
397     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
398     IF (.NOT.implicitViscosity)
399     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
400     ENDIF
401     #endif
402    
403 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
404 jmc 1.23 IF (no_slip_sides) THEN
405 adcroft 1.1 C- No-slip BCs impose a drag at walls...
406 jmc 1.23 CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
407     DO j=jMin,jMax
408     DO i=iMin,iMax
409     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
410     ENDDO
411     ENDDO
412     ENDIF
413 adcroft 1.1 C- No-slip BCs impose a drag at bottom
414 jmc 1.23 IF (bottomDragTerms) THEN
415     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
416     DO j=jMin,jMax
417     DO i=iMin,iMax
418     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
419     ENDDO
420     ENDDO
421     ENDIF
422    
423     C- endif momViscosity
424 adcroft 1.1 ENDIF
425    
426 jmc 1.12 C-- Forcing term (moved to timestep.F)
427     c IF (momForcing)
428     c & CALL EXTERNAL_FORCING_U(
429     c I iMin,iMax,jMin,jMax,bi,bj,k,
430     c I myTime,myThid)
431 adcroft 1.1
432     C-- Metric terms for curvilinear grid systems
433 adcroft 1.5 IF (useNHMTerms) THEN
434     C o Non-hydrosatic metric terms
435 adcroft 1.1 CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
436     DO j=jMin,jMax
437     DO i=iMin,iMax
438     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
439     ENDDO
440     ENDDO
441 adcroft 1.5 ENDIF
442     IF (usingSphericalPolarMTerms) THEN
443 adcroft 1.1 CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
444     DO j=jMin,jMax
445     DO i=iMin,iMax
446     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
447     ENDDO
448     ENDDO
449 afe 1.20 ENDIF
450 afe 1.19 IF (usingCylindricalGrid) THEN
451     CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
452     DO j=jMin,jMax
453     DO i=iMin,iMax
454     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
455     ENDDO
456     ENDDO
457 adcroft 1.1 ENDIF
458    
459 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
460 adcroft 1.1
461     C---- Meridional momentum equation starts here
462    
463 jmc 1.23 IF (momAdvection) THEN
464     C--- Calculate mean fluxes (advection) between cells for meridional flow.
465     C Mean flow component of zonal flux -> fZon
466     CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,fZon,myThid)
467 adcroft 1.1
468     C-- Meridional flux (fMer is at north face of "v" cell)
469 jmc 1.23 C Mean flow component of meridional flux -> fMer
470     CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,fMer,myThid)
471 adcroft 1.1
472     C-- Vertical flux (fVer is at upper face of "v" cell)
473 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVerV
474     CALL MOM_V_ADV_WV(
475     I bi,bj,k+1,vVel,wVel,rTransV,
476     O fVerV(1-OLx,1-OLy,kDown), myThid )
477 adcroft 1.1
478     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
479 jmc 1.23 DO j=jMin,jMax
480     DO i=iMin,iMax
481     gV(i,j,k,bi,bj) =
482 adcroft 1.1 #ifdef OLD_UV_GEOM
483 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
484     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
485 adcroft 1.1 #else
486 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
487     & *recip_rAs(i,j,bi,bj)
488 adcroft 1.1 #endif
489 jmc 1.23 & *( ( fZon(i+1,j) - fZon(i,j ) )*uDvdxFac
490     & +( fMer(i, j) - fMer(i,j-1) )*vDvdyFac
491     & +(fVerV(i,j,kDown) - fVerV(i,j,kUp))*rkSign*rVelDvdrFac
492     & )
493 jmc 1.24 ENDDO
494     ENDDO
495    
496     #ifdef ALLOW_DIAGNOSTICS
497     IF ( useDiagnostics ) THEN
498     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Vm ',k,1,2,bi,bj,myThid)
499     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Vm ',k,1,2,bi,bj,myThid)
500     CALL DIAGNOSTICS_FILL(fVerV(1-Olx,1-Oly,kUp),
501     & 'ADVrE_Vm',k,1,2,bi,bj,myThid)
502     ENDIF
503     #endif
504 adcroft 1.1
505 jmc 1.8 #ifdef NONLIN_FRSURF
506     C-- account for 3.D divergence of the flow in rStar coordinate:
507 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
508     DO j=jMin,jMax
509     DO i=iMin,iMax
510     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
511 jmc 1.8 & - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
512     & *vVel(i,j,k,bi,bj)
513 jmc 1.23 ENDDO
514     ENDDO
515     ENDIF
516     IF ( select_rStar.LT.0 ) THEN
517     DO j=jMin,jMax
518     DO i=iMin,iMax
519     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
520     & - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
521     ENDDO
522     ENDDO
523     ENDIF
524     #endif /* NONLIN_FRSURF */
525    
526     ELSE
527     C- if momAdvection / else
528     DO j=1-OLy,sNy+OLy
529     DO i=1-OLx,sNx+OLx
530     gV(i,j,k,bi,bj) = 0. _d 0
531     ENDDO
532 jmc 1.8 ENDDO
533 jmc 1.23
534     C- endif momAdvection.
535 jmc 1.8 ENDIF
536 jmc 1.23
537     IF (momViscosity) THEN
538     C--- Calculate eddy fluxes (dissipation) between cells for meridional flow.
539     C Bi-harmonic term del^2 V -> v4F
540 baylor 1.25 IF (biharmonic)
541 jmc 1.23 & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
542    
543     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
544 baylor 1.25 CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,
545     I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
546 jmc 1.23
547     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
548 baylor 1.25 CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,
549     I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
550 jmc 1.23
551     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
552     IF (.NOT.implicitViscosity) THEN
553     CALL MOM_V_RVISCFLUX(bi,bj, k, vVel,KappaRV,fVrUp,myThid)
554     CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,fVrDw,myThid)
555     ENDIF
556    
557     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
558     DO j=jMin,jMax
559     DO i=iMin,iMax
560     gvDiss(i,j) =
561     #ifdef OLD_UV_GEOM
562     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
563     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
564     #else
565     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
566     & *recip_rAs(i,j,bi,bj)
567     #endif
568     & *( ( fZon(i+1,j) - fZon(i,j ) )*AhDvdxFac
569     & +( fMer(i, j) - fMer(i,j-1) )*AhDvdyFac
570     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDvdrFac
571     & )
572     ENDDO
573 jmc 1.8 ENDDO
574    
575 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
576     IF ( useDiagnostics ) THEN
577     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
578     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
579     IF (.NOT.implicitViscosity)
580     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
581     ENDIF
582     #endif
583    
584 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
585 jmc 1.23 IF (no_slip_sides) THEN
586 adcroft 1.1 C- No-slip BCs impose a drag at walls...
587 jmc 1.23 CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid)
588     DO j=jMin,jMax
589     DO i=iMin,iMax
590     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
591     ENDDO
592     ENDDO
593     ENDIF
594 adcroft 1.1 C- No-slip BCs impose a drag at bottom
595 jmc 1.23 IF (bottomDragTerms) THEN
596     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
597     DO j=jMin,jMax
598     DO i=iMin,iMax
599     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
600     ENDDO
601     ENDDO
602     ENDIF
603    
604     C- endif momViscosity
605 adcroft 1.1 ENDIF
606    
607 jmc 1.12 C-- Forcing term (moved to timestep.F)
608     c IF (momForcing)
609     c & CALL EXTERNAL_FORCING_V(
610     c I iMin,iMax,jMin,jMax,bi,bj,k,
611     c I myTime,myThid)
612 adcroft 1.1
613     C-- Metric terms for curvilinear grid systems
614 adcroft 1.5 IF (useNHMTerms) THEN
615 adcroft 1.1 C o Spherical polar grid metric terms
616     CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
617     DO j=jMin,jMax
618     DO i=iMin,iMax
619     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
620     ENDDO
621     ENDDO
622 adcroft 1.5 ENDIF
623     IF (usingSphericalPolarMTerms) THEN
624 adcroft 1.1 CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
625     DO j=jMin,jMax
626     DO i=iMin,iMax
627     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
628     ENDDO
629     ENDDO
630     ENDIF
631 afe 1.19 IF (usingCylindricalGrid) THEN
632     CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
633     DO j=jMin,jMax
634     DO i=iMin,iMax
635     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
636     ENDDO
637     ENDDO
638     ENDIF
639 adcroft 1.1
640 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
641 adcroft 1.1
642     C-- Coriolis term
643     C Note. As coded here, coriolis will not work with "thin walls"
644 jmc 1.12 c IF (useCDscheme) THEN
645     c CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
646     c ELSE
647     IF (.NOT.useCDscheme) THEN
648     CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
649     DO j=jMin,jMax
650     DO i=iMin,iMax
651     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
652     ENDDO
653     ENDDO
654 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
655     IF ( useDiagnostics )
656     & CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)
657     #endif
658 jmc 1.12 CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
659     DO j=jMin,jMax
660     DO i=iMin,iMax
661     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
662     ENDDO
663     ENDDO
664 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
665     IF ( useDiagnostics )
666     & CALL DIAGNOSTICS_FILL(cf,'Vm_Cori ',k,1,2,bi,bj,myThid)
667     #endif
668 jmc 1.12 ENDIF
669    
670 adcroft 1.7 IF (nonHydrostatic.OR.quasiHydrostatic) THEN
671 adcroft 1.6 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
672     DO j=jMin,jMax
673     DO i=iMin,iMax
674     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
675     ENDDO
676     ENDDO
677     ENDIF
678 adcroft 1.1
679 jmc 1.23 C-- Set du/dt & dv/dt on boundaries to zero
680     DO j=jMin,jMax
681     DO i=iMin,iMax
682     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
683     guDiss(i,j) = guDiss(i,j) *_maskW(i,j,k,bi,bj)
684     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
685     gvDiss(i,j) = gvDiss(i,j) *_maskS(i,j,k,bi,bj)
686     ENDDO
687     ENDDO
688    
689 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
690     IF ( useDiagnostics ) THEN
691     IF (bottomDragTerms)
692     & CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
693     CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
694     & 'Um_Advec',k,1,2,bi,bj,myThid)
695     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
696     & 'Vm_Advec',k,1,2,bi,bj,myThid)
697     IF (momViscosity) THEN
698     CALL DIAGNOSTICS_FILL(guDiss,'Um_Diss ',k,1,2,bi,bj,myThid)
699     CALL DIAGNOSTICS_FILL(gvDiss,'Vm_Diss ',k,1,2,bi,bj,myThid)
700     ENDIF
701     ENDIF
702     #endif /* ALLOW_DIAGNOSTICS */
703    
704 adcroft 1.1 RETURN
705     END
  ViewVC Help
Powered by ViewVC 1.1.22