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

Contents of /MITgcm/pkg/mom_vecinv/mom_vecinv.F

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


Revision 1.64 - (show annotations) (download)
Mon Apr 25 02:03:13 2011 UTC (13 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62w, checkpoint62y, checkpoint62x
Changes since 1.63: +11 -1 lines
add output (function of diagFreq) of intermediate var. del2u,del2v,dStar & zStar

1 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vecinv.F,v 1.63 2010/03/16 00:16:50 jmc Exp $
2 C $Name: $
3
4 #include "MOM_VECINV_OPTIONS.h"
5
6 SUBROUTINE MOM_VECINV(
7 I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
8 I KappaRU, KappaRV,
9 U fVerU, fVerV,
10 O guDiss, gvDiss,
11 I myTime, myIter, myThid)
12 C /==========================================================\
13 C | S/R MOM_VECINV |
14 C | o Form the right hand-side of the momentum equation. |
15 C |==========================================================|
16 C | Terms are evaluated one layer at a time working from |
17 C | the bottom to the top. The vertically integrated |
18 C | barotropic flow tendency term is evluated by summing the |
19 C | tendencies. |
20 C | Notes: |
21 C | We have not sorted out an entirely satisfactory formula |
22 C | for the diffusion equation bc with lopping. The present |
23 C | form produces a diffusive flux that does not scale with |
24 C | open-area. Need to do something to solidfy this and to |
25 C | deal "properly" with thin walls. |
26 C \==========================================================/
27 IMPLICIT NONE
28
29 C == Global variables ==
30 #include "SIZE.h"
31 #include "DYNVARS.h"
32 #include "EEPARAMS.h"
33 #include "PARAMS.h"
34 #ifdef ALLOW_MNC
35 #include "MNC_PARAMS.h"
36 #endif
37 #include "GRID.h"
38 #ifdef ALLOW_TIMEAVE
39 #include "TIMEAVE_STATV.h"
40 #endif
41 #ifdef ALLOW_AUTODIFF_TAMC
42 # include "tamc.h"
43 # include "tamc_keys.h"
44 #endif
45
46 C == Routine arguments ==
47 C fVerU :: Flux of momentum in the vertical direction, out of the upper
48 C fVerV :: face of a cell K ( flux into the cell above ).
49 C guDiss :: dissipation tendency (all explicit terms), u component
50 C gvDiss :: dissipation tendency (all explicit terms), v component
51 C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
52 C results will be set.
53 C kUp, kDown - Index for upper and lower layers.
54 C myThid :: my Thread Id number
55 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56 _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
57 _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
58 _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
59 _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60 _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 INTEGER kUp,kDown
62 _RL myTime
63 INTEGER myIter
64 INTEGER myThid
65 INTEGER bi,bj,iMin,iMax,jMin,jMax
66
67 #ifdef ALLOW_MOM_VECINV
68
69 C == Functions ==
70 LOGICAL DIFFERENT_MULTIPLE
71 EXTERNAL DIFFERENT_MULTIPLE
72
73 C == Local variables ==
74 _RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RS hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82 _RL del2u (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 _RL del2v (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84 _RL dStar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85 _RL zStar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86 _RL tension (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL strain (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL KE (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL omega3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL vort3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL hDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 C i,j,k :: Loop counters
97 INTEGER i,j,k
98 C xxxFac - On-off tracer parameters used for switching terms off.
99 _RL ArDudrFac
100 _RL ArDvdrFac
101 _RL sideMaskFac
102 LOGICAL bottomDragTerms
103 LOGICAL writeDiag
104 LOGICAL harmonic,biharmonic,useVariableViscosity
105 #ifdef ALLOW_AUTODIFF_TAMC
106 INTEGER imomkey
107 #endif
108
109 #ifdef ALLOW_MNC
110 INTEGER offsets(9)
111 CHARACTER*(1) pf
112 #endif
113
114 #ifdef ALLOW_AUTODIFF_TAMC
115 C-- only the kDown part of fverU/V is set in this subroutine
116 C-- the kUp is still required
117 C-- In the case of mom_fluxform Kup is set as well
118 C-- (at least in part)
119 fVerU(1,1,kUp) = fVerU(1,1,kUp)
120 fVerV(1,1,kUp) = fVerV(1,1,kUp)
121 #endif
122
123 #ifdef ALLOW_AUTODIFF_TAMC
124 act0 = k - 1
125 max0 = Nr
126 act1 = bi - myBxLo(myThid)
127 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
128 act2 = bj - myByLo(myThid)
129 max2 = myByHi(myThid) - myByLo(myThid) + 1
130 act3 = myThid - 1
131 max3 = nTx*nTy
132 act4 = ikey_dynamics - 1
133 imomkey = (act0 + 1)
134 & + act1*max0
135 & + act2*max0*max1
136 & + act3*max0*max1*max2
137 & + act4*max0*max1*max2*max3
138 #endif /* ALLOW_AUTODIFF_TAMC */
139
140 writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
141
142 #ifdef ALLOW_MNC
143 IF (useMNC .AND. snapshot_mnc .AND. writeDiag) THEN
144 IF ( writeBinaryPrec .EQ. precFloat64 ) THEN
145 pf(1:1) = 'D'
146 ELSE
147 pf(1:1) = 'R'
148 ENDIF
149 IF ((bi .EQ. 1).AND.(bj .EQ. 1).AND.(k .EQ. 1)) THEN
150 CALL MNC_CW_SET_UDIM('mom_vi', -1, myThid)
151 CALL MNC_CW_RL_W_S('D','mom_vi',0,0,'T',myTime,myThid)
152 CALL MNC_CW_SET_UDIM('mom_vi', 0, myThid)
153 CALL MNC_CW_I_W_S('I','mom_vi',0,0,'iter',myIter,myThid)
154 ENDIF
155 DO i = 1,9
156 offsets(i) = 0
157 ENDDO
158 offsets(3) = k
159 c write(*,*) 'offsets = ',(offsets(i),i=1,9)
160 ENDIF
161 #endif /* ALLOW_MNC */
162
163 C-- Initialise intermediate terms
164 DO j=1-OLy,sNy+OLy
165 DO i=1-OLx,sNx+OLx
166 vF(i,j) = 0.
167 vrF(i,j) = 0.
168 uCf(i,j) = 0.
169 vCf(i,j) = 0.
170 del2u(i,j) = 0.
171 del2v(i,j) = 0.
172 dStar(i,j) = 0.
173 zStar(i,j) = 0.
174 guDiss(i,j)= 0.
175 gvDiss(i,j)= 0.
176 vort3(i,j) = 0.
177 omega3(i,j)= 0.
178 KE(i,j) = 0.
179 C- need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
180 hDiv(i,j) = 0.
181 viscAh_Z(i,j) = 0.
182 viscAh_D(i,j) = 0.
183 viscA4_Z(i,j) = 0.
184 viscA4_D(i,j) = 0.
185
186 strain(i,j) = 0. _d 0
187 tension(i,j) = 0. _d 0
188 #ifdef ALLOW_AUTODIFF_TAMC
189 hFacZ(i,j) = 0. _d 0
190 #endif
191 ENDDO
192 ENDDO
193
194 C-- Term by term tracer parmeters
195 C o U momentum equation
196 ArDudrFac = vfFacMom*1.
197 C o V momentum equation
198 ArDvdrFac = vfFacMom*1.
199
200 C note: using standard stencil (no mask) results in under-estimating
201 C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
202 IF ( no_slip_sides ) THEN
203 sideMaskFac = sideDragFactor
204 ELSE
205 sideMaskFac = 0. _d 0
206 ENDIF
207
208 IF ( no_slip_bottom
209 & .OR. bottomDragQuadratic.NE.0.
210 & .OR. bottomDragLinear.NE.0.) THEN
211 bottomDragTerms=.TRUE.
212 ELSE
213 bottomDragTerms=.FALSE.
214 ENDIF
215
216 C-- Calculate open water fraction at vorticity points
217 CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
218
219 C Make local copies of horizontal flow field
220 DO j=1-OLy,sNy+OLy
221 DO i=1-OLx,sNx+OLx
222 uFld(i,j) = uVel(i,j,k,bi,bj)
223 vFld(i,j) = vVel(i,j,k,bi,bj)
224 ENDDO
225 ENDDO
226
227 C note (jmc) : Dissipation and Vort3 advection do not necesary
228 C use the same maskZ (and hFacZ) => needs 2 call(s)
229 c CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFacZ,r_hFacZ,myThid)
230
231 CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
232
233 CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
234
235 IF (momViscosity) THEN
236 C-- For viscous term, compute horizontal divergence, tension & strain
237 C and mask relative vorticity (free-slip case):
238
239 #ifdef ALLOW_AUTODIFF_TAMC
240 CADJ STORE vort3(:,:) =
241 CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
242 #endif
243
244 CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
245
246 CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
247
248 CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
249
250 C- account for no-slip / free-slip BC:
251 DO j=1-Oly,sNy+Oly
252 DO i=1-Olx,sNx+Olx
253 IF ( hFacZ(i,j).EQ.0. ) THEN
254 vort3(i,j) = sideMaskFac*vort3(i,j)
255 strain(i,j) = sideMaskFac*strain(i,j)
256 ENDIF
257 ENDDO
258 ENDDO
259
260 C-- Calculate Viscosities
261 CALL MOM_CALC_VISC(
262 I bi,bj,k,
263 O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
264 O harmonic,biharmonic,useVariableViscosity,
265 I hDiv,vort3,tension,strain,KE,hfacZ,
266 I myThid)
267
268 C Calculate del^2 u and del^2 v for bi-harmonic term
269 IF (biharmonic) THEN
270 CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
271 O del2u,del2v,
272 & myThid)
273 CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
274 CALL MOM_CALC_RELVORT3(bi,bj,k,
275 & del2u,del2v,hFacZ,zStar,myThid)
276 IF ( writeDiag ) THEN
277 CALL WRITE_LOCAL_RL( 'del2u', 'I10', 1, del2u,
278 & bi,bj,k, myIter, myThid )
279 CALL WRITE_LOCAL_RL( 'del2v', 'I10', 1, del2v,
280 & bi,bj,k, myIter, myThid )
281 CALL WRITE_LOCAL_RL( 'dStar', 'I10', 1, dStar,
282 & bi,bj,k, myIter, myThid )
283 CALL WRITE_LOCAL_RL( 'zStar', 'I10', 1, zStar,
284 & bi,bj,k, myIter, myThid )
285 ENDIF
286 ENDIF
287
288 C- Strain diagnostics:
289 IF ( writeDiag ) THEN
290 IF (snapshot_mdsio) THEN
291 CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)
292 ENDIF
293 #ifdef ALLOW_MNC
294 IF (useMNC .AND. snapshot_mnc) THEN
295 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strain,
296 & offsets, myThid)
297 ENDIF
298 #endif /* ALLOW_MNC */
299 ENDIF
300 #ifdef ALLOW_DIAGNOSTICS
301 IF ( useDiagnostics ) THEN
302 CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
303 ENDIF
304 #endif /* ALLOW_DIAGNOSTICS */
305
306 C--- Calculate dissipation terms for U and V equations
307
308 C in terms of tension and strain
309 IF (useStrainTensionVisc) THEN
310 C mask strain as if free-slip since side-drag is computed separately
311 DO j=1-Oly,sNy+Oly
312 DO i=1-Olx,sNx+Olx
313 IF ( hFacZ(i,j).EQ.0. ) strain(i,j) = 0. _d 0
314 ENDDO
315 ENDDO
316 CALL MOM_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
317 I hFacZ,
318 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
319 I harmonic,biharmonic,useVariableViscosity,
320 O guDiss,gvDiss,
321 I myThid)
322 ELSE
323 C in terms of vorticity and divergence
324 CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
325 I hFacZ,dStar,zStar,
326 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
327 I harmonic,biharmonic,useVariableViscosity,
328 O guDiss,gvDiss,
329 & myThid)
330 ENDIF
331 C-- if (momViscosity) end of block.
332 ENDIF
333
334 C- Return to standard hfacZ (min-4) and mask vort3 accordingly:
335 c CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
336
337 C--- Other dissipation terms in Zonal momentum equation
338
339 C-- Vertical flux (fVer is at upper face of "u" cell)
340
341 C Eddy component of vertical flux (interior component only) -> vrF
342 IF (momViscosity.AND..NOT.implicitViscosity) THEN
343 CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,vrF,myThid)
344
345 C Combine fluxes
346 DO j=jMin,jMax
347 DO i=iMin,iMax
348 fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)
349 ENDDO
350 ENDDO
351
352 C-- Tendency is minus divergence of the fluxes
353 DO j=2-Oly,sNy+Oly-1
354 DO i=2-Olx,sNx+Olx-1
355 guDiss(i,j) = guDiss(i,j)
356 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
357 & *recip_rAw(i,j,bi,bj)
358 & *(
359 & fVerU(i,j,kDown) - fVerU(i,j,kUp)
360 & )*rkSign
361 ENDDO
362 ENDDO
363 ENDIF
364
365 C-- No-slip and drag BCs appear as body forces in cell abutting topography
366 IF (momViscosity.AND.no_slip_sides) THEN
367 C- No-slip BCs impose a drag at walls...
368 CALL MOM_U_SIDEDRAG(
369 I bi,bj,k,
370 I uFld, del2u, hFacZ,
371 I viscAh_Z,viscA4_Z,
372 I harmonic,biharmonic,useVariableViscosity,
373 O vF,
374 I myThid)
375 DO j=jMin,jMax
376 DO i=iMin,iMax
377 guDiss(i,j) = guDiss(i,j)+vF(i,j)
378 ENDDO
379 ENDDO
380 ENDIF
381 C- No-slip BCs impose a drag at bottom
382 IF (momViscosity.AND.bottomDragTerms) THEN
383 CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
384 DO j=jMin,jMax
385 DO i=iMin,iMax
386 guDiss(i,j) = guDiss(i,j)+vF(i,j)
387 ENDDO
388 ENDDO
389 ENDIF
390 #ifdef ALLOW_SHELFICE
391 IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
392 CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
393 DO j=jMin,jMax
394 DO i=iMin,iMax
395 guDiss(i,j) = guDiss(i,j) + vF(i,j)
396 ENDDO
397 ENDDO
398 ENDIF
399 #endif /* ALLOW_SHELFICE */
400
401
402 C--- Other dissipation terms in Meridional momentum equation
403
404 C-- Vertical flux (fVer is at upper face of "v" cell)
405
406 C Eddy component of vertical flux (interior component only) -> vrF
407 IF (momViscosity.AND..NOT.implicitViscosity) THEN
408 CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,vrF,myThid)
409
410 C Combine fluxes -> fVerV
411 DO j=jMin,jMax
412 DO i=iMin,iMax
413 fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)
414 ENDDO
415 ENDDO
416
417 C-- Tendency is minus divergence of the fluxes
418 DO j=jMin,jMax
419 DO i=iMin,iMax
420 gvDiss(i,j) = gvDiss(i,j)
421 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
422 & *recip_rAs(i,j,bi,bj)
423 & *(
424 & fVerV(i,j,kDown) - fVerV(i,j,kUp)
425 & )*rkSign
426 ENDDO
427 ENDDO
428 ENDIF
429
430 C-- No-slip and drag BCs appear as body forces in cell abutting topography
431 IF (momViscosity.AND.no_slip_sides) THEN
432 C- No-slip BCs impose a drag at walls...
433 CALL MOM_V_SIDEDRAG(
434 I bi,bj,k,
435 I vFld, del2v, hFacZ,
436 I viscAh_Z,viscA4_Z,
437 I harmonic,biharmonic,useVariableViscosity,
438 O vF,
439 I myThid)
440 DO j=jMin,jMax
441 DO i=iMin,iMax
442 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
443 ENDDO
444 ENDDO
445 ENDIF
446 C- No-slip BCs impose a drag at bottom
447 IF (momViscosity.AND.bottomDragTerms) THEN
448 CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
449 DO j=jMin,jMax
450 DO i=iMin,iMax
451 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
452 ENDDO
453 ENDDO
454 ENDIF
455 #ifdef ALLOW_SHELFICE
456 IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
457 CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
458 DO j=jMin,jMax
459 DO i=iMin,iMax
460 gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
461 ENDDO
462 ENDDO
463 ENDIF
464 #endif /* ALLOW_SHELFICE */
465
466
467 C- Vorticity diagnostics:
468 IF ( writeDiag ) THEN
469 IF (snapshot_mdsio) THEN
470 CALL WRITE_LOCAL_RL('Z3','I10',1,vort3, bi,bj,k,myIter,myThid)
471 ENDIF
472 #ifdef ALLOW_MNC
473 IF (useMNC .AND. snapshot_mnc) THEN
474 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3,
475 & offsets, myThid)
476 ENDIF
477 #endif /* ALLOW_MNC */
478 ENDIF
479 #ifdef ALLOW_DIAGNOSTICS
480 IF ( useDiagnostics ) THEN
481 CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
482 ENDIF
483 #endif /* ALLOW_DIAGNOSTICS */
484
485 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
486
487 C--- Prepare for Advection & Coriolis terms:
488 C- Mask relative vorticity and calculate absolute vorticity
489 DO j=1-Oly,sNy+Oly
490 DO i=1-Olx,sNx+Olx
491 IF ( hFacZ(i,j).EQ.0. ) vort3(i,j) = 0.
492 ENDDO
493 ENDDO
494 IF (useAbsVorticity)
495 & CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
496
497 C-- Horizontal Coriolis terms
498 c IF (useCoriolis .AND. .NOT.useCDscheme
499 c & .AND. .NOT. useAbsVorticity) THEN
500 C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
501 IF ( useCoriolis .AND.
502 & .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
503 & ) THEN
504 IF (useAbsVorticity) THEN
505 CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
506 & uCf,myThid)
507 CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
508 & vCf,myThid)
509 ELSE
510 CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
511 & uCf,vCf,myThid)
512 ENDIF
513 DO j=jMin,jMax
514 DO i=iMin,iMax
515 gU(i,j,k,bi,bj) = uCf(i,j)
516 gV(i,j,k,bi,bj) = vCf(i,j)
517 ENDDO
518 ENDDO
519 IF ( writeDiag ) THEN
520 IF (snapshot_mdsio) THEN
521 CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
522 CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
523 ENDIF
524 #ifdef ALLOW_MNC
525 IF (useMNC .AND. snapshot_mnc) THEN
526 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
527 & offsets, myThid)
528 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
529 & offsets, myThid)
530 ENDIF
531 #endif /* ALLOW_MNC */
532 ENDIF
533 #ifdef ALLOW_DIAGNOSTICS
534 IF ( useDiagnostics ) THEN
535 CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
536 CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
537 ENDIF
538 #endif /* ALLOW_DIAGNOSTICS */
539 ELSE
540 DO j=jMin,jMax
541 DO i=iMin,iMax
542 gU(i,j,k,bi,bj) = 0. _d 0
543 gV(i,j,k,bi,bj) = 0. _d 0
544 ENDDO
545 ENDDO
546 ENDIF
547
548 IF (momAdvection) THEN
549 C-- Horizontal advection of relative (or absolute) vorticity
550 IF ( (highOrderVorticity.OR.upwindVorticity)
551 & .AND.useAbsVorticity ) THEN
552 CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,
553 & uCf,myThid)
554 ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
555 CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,
556 & uCf,myThid)
557 ELSEIF ( useAbsVorticity ) THEN
558 CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
559 & uCf,myThid)
560 ELSE
561 CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,
562 & uCf,myThid)
563 ENDIF
564 DO j=jMin,jMax
565 DO i=iMin,iMax
566 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
567 ENDDO
568 ENDDO
569 IF ( (highOrderVorticity.OR.upwindVorticity)
570 & .AND.useAbsVorticity ) THEN
571 CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ,
572 & vCf,myThid)
573 ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
574 CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ,
575 & vCf,myThid)
576 ELSEIF ( useAbsVorticity ) THEN
577 CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
578 & vCf,myThid)
579 ELSE
580 CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,
581 & vCf,myThid)
582 ENDIF
583 DO j=jMin,jMax
584 DO i=iMin,iMax
585 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
586 ENDDO
587 ENDDO
588
589 IF ( writeDiag ) THEN
590 IF (snapshot_mdsio) THEN
591 CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
592 CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
593 ENDIF
594 #ifdef ALLOW_MNC
595 IF (useMNC .AND. snapshot_mnc) THEN
596 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
597 & offsets, myThid)
598 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
599 & offsets, myThid)
600 ENDIF
601 #endif /* ALLOW_MNC */
602 ENDIF
603
604 #ifdef ALLOW_TIMEAVE
605 IF (taveFreq.GT.0.) THEN
606 CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
607 & Nr, k, bi, bj, myThid)
608 CALL TIMEAVE_CUMUL_1K1T(vZetatave,uCf,deltaTClock,
609 & Nr, k, bi, bj, myThid)
610 ENDIF
611 #endif /* ALLOW_TIMEAVE */
612 #ifdef ALLOW_DIAGNOSTICS
613 IF ( useDiagnostics ) THEN
614 CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
615 CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
616 ENDIF
617 #endif /* ALLOW_DIAGNOSTICS */
618
619 C-- Vertical shear terms (-w*du/dr & -w*dv/dr)
620 IF ( .NOT. momImplVertAdv ) THEN
621 CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)
622 DO j=jMin,jMax
623 DO i=iMin,iMax
624 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
625 ENDDO
626 ENDDO
627 CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)
628 DO j=jMin,jMax
629 DO i=iMin,iMax
630 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
631 ENDDO
632 ENDDO
633 #ifdef ALLOW_DIAGNOSTICS
634 IF ( useDiagnostics ) THEN
635 CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
636 CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
637 ENDIF
638 #endif /* ALLOW_DIAGNOSTICS */
639 ENDIF
640
641 C-- Bernoulli term
642 CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)
643 DO j=jMin,jMax
644 DO i=iMin,iMax
645 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
646 ENDDO
647 ENDDO
648 CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)
649 DO j=jMin,jMax
650 DO i=iMin,iMax
651 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
652 ENDDO
653 ENDDO
654 IF ( writeDiag ) THEN
655 IF (snapshot_mdsio) THEN
656 CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
657 CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
658 ENDIF
659 #ifdef ALLOW_MNC
660 IF (useMNC .AND. snapshot_mnc) THEN
661 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
662 & offsets, myThid)
663 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
664 & offsets, myThid)
665 ENDIF
666 #endif /* ALLOW_MNC */
667 ENDIF
668
669 C-- end if momAdvection
670 ENDIF
671
672 C-- 3.D Coriolis term (horizontal momentum, Eastward component: -fprime*w)
673 IF ( use3dCoriolis ) THEN
674 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
675 DO j=jMin,jMax
676 DO i=iMin,iMax
677 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
678 ENDDO
679 ENDDO
680 IF ( usingCurvilinearGrid ) THEN
681 C- presently, non zero angleSinC array only supported with Curvilinear-Grid
682 CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
683 DO j=jMin,jMax
684 DO i=iMin,iMax
685 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
686 ENDDO
687 ENDDO
688 ENDIF
689 ENDIF
690
691 C-- Non-Hydrostatic (spherical) metric terms
692 IF ( useNHMTerms ) THEN
693 CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
694 DO j=jMin,jMax
695 DO i=iMin,iMax
696 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
697 ENDDO
698 ENDDO
699 CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
700 DO j=jMin,jMax
701 DO i=iMin,iMax
702 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
703 ENDDO
704 ENDDO
705 ENDIF
706
707 C-- Set du/dt & dv/dt on boundaries to zero
708 DO j=jMin,jMax
709 DO i=iMin,iMax
710 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
711 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
712 ENDDO
713 ENDDO
714
715 #ifdef ALLOW_DEBUG
716 IF ( debugLevel .GE. debLevB
717 & .AND. k.EQ.4 .AND. myIter.EQ.nIter0
718 & .AND. nPx.EQ.1 .AND. nPy.EQ.1
719 & .AND. useCubedSphereExchange ) THEN
720 CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
721 & guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
722 ENDIF
723 #endif /* ALLOW_DEBUG */
724
725 IF ( writeDiag ) THEN
726 IF (snapshot_mdsio) THEN
727 CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
728 CALL WRITE_LOCAL_RL('KE','I10',1,KE, bi,bj,k,myIter,myThid)
729 CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv, bi,bj,k,myIter,myThid)
730 CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
731 CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
732 CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
733 ENDIF
734 #ifdef ALLOW_MNC
735 IF (useMNC .AND. snapshot_mnc) THEN
736 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
737 & offsets, myThid)
738 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
739 & offsets, myThid)
740 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
741 & offsets, myThid)
742 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
743 & offsets, myThid)
744 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
745 & offsets, myThid)
746 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
747 & offsets, myThid)
748 ENDIF
749 #endif /* ALLOW_MNC */
750 ENDIF
751
752 #ifdef ALLOW_DIAGNOSTICS
753 IF ( useDiagnostics ) THEN
754 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
755 IF (momViscosity) THEN
756 CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
757 CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
758 CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)
759 CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)
760 ENDIF
761 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
762 & 'Um_Advec',k,1,2,bi,bj,myThid)
763 CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
764 & 'Vm_Advec',k,1,2,bi,bj,myThid)
765 ENDIF
766 #endif /* ALLOW_DIAGNOSTICS */
767
768 #endif /* ALLOW_MOM_VECINV */
769
770 RETURN
771 END

  ViewVC Help
Powered by ViewVC 1.1.22