/[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.71 - (show annotations) (download)
Sun Feb 9 18:46:46 2014 UTC (10 years, 4 months ago) by jmc
Branch: MAIN
Changes since 1.70: +22 -3 lines
fix sideDrag option for thin-walls with Non-Lin Free-Surf
using 2nd hFacZ that is computed from initial (fix domain) hFac

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

  ViewVC Help
Powered by ViewVC 1.1.22