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

Contents of /MITgcm/pkg/mom_fluxform/mom_fluxform.F

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


Revision 1.36 - (show annotations) (download)
Thu May 4 12:29:07 2006 UTC (18 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58f_post, checkpoint58e_post, checkpoint58k_post, checkpoint58l_post, checkpoint58g_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post
Changes since 1.35: +5 -1 lines
Need to bracket STORE block of rstar-related variables.

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

  ViewVC Help
Powered by ViewVC 1.1.22