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

Annotation of /MITgcm/pkg/mom_vecinv/mom_vecinv.F

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


Revision 1.62 - (hide annotations) (download)
Tue Apr 1 01:27:33 2008 UTC (16 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62a, checkpoint60, checkpoint61, checkpoint62, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint62b, checkpoint61f, checkpoint61n, checkpoint61q, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.61: +9 -7 lines
clarify highOrderVorticity & upwindVorticity (now exclusive);
mom_vi_u/v_coriolis_c4.F now also deal with upwindVorticity ;

1 jmc 1.62 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vecinv.F,v 1.61 2007/11/27 20:20:26 jmc Exp $
2 adcroft 1.2 C $Name: $
3 adcroft 1.1
4 adcroft 1.21 #include "MOM_VECINV_OPTIONS.h"
5 adcroft 1.1
6 jmc 1.57 SUBROUTINE MOM_VECINV(
7 adcroft 1.1 I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
8 jmc 1.43 I KappaRU, KappaRV,
9 adcroft 1.1 U fVerU, fVerV,
10 jmc 1.31 O guDiss, gvDiss,
11 jmc 1.15 I myTime, myIter, myThid)
12 adcroft 1.1 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 edhill 1.27 #ifdef ALLOW_MNC
35     #include "MNC_PARAMS.h"
36     #endif
37 adcroft 1.1 #include "GRID.h"
38 jmc 1.7 #ifdef ALLOW_TIMEAVE
39     #include "TIMEAVE_STATV.h"
40     #endif
41 heimbach 1.59 #ifdef ALLOW_AUTODIFF_TAMC
42     # include "tamc.h"
43     # include "tamc_keys.h"
44     #endif
45 adcroft 1.1
46     C == Routine arguments ==
47 jmc 1.31 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 adcroft 1.1 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 jmc 1.54 C myThid :: my Thread Id number
55 adcroft 1.1 _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 jmc 1.31 _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60     _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 adcroft 1.1 INTEGER kUp,kDown
62 jmc 1.15 _RL myTime
63 adcroft 1.2 INTEGER myIter
64 adcroft 1.1 INTEGER myThid
65     INTEGER bi,bj,iMin,iMax,jMin,jMax
66    
67 edhill 1.11 #ifdef ALLOW_MOM_VECINV
68 jmc 1.7
69 adcroft 1.2 C == Functions ==
70 jmc 1.38 LOGICAL DIFFERENT_MULTIPLE
71     EXTERNAL DIFFERENT_MULTIPLE
72 adcroft 1.2
73 adcroft 1.1 C == Local variables ==
74     _RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 jmc 1.54 _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 adcroft 1.1 INTEGER i,j,k
98     C xxxFac - On-off tracer parameters used for switching terms off.
99     _RL ArDudrFac
100     _RL ArDvdrFac
101 jmc 1.54 _RL sideMaskFac
102 adcroft 1.1 LOGICAL bottomDragTerms
103 jmc 1.15 LOGICAL writeDiag
104 baylor 1.50 LOGICAL harmonic,biharmonic,useVariableViscosity
105 heimbach 1.59 #ifdef ALLOW_AUTODIFF_TAMC
106     INTEGER imomkey
107     #endif
108 adcroft 1.1
109 edhill 1.25 #ifdef ALLOW_MNC
110     INTEGER offsets(9)
111 edhill 1.53 CHARACTER*(1) pf
112 edhill 1.25 #endif
113    
114 heimbach 1.9 #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 heimbach 1.59 #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 jmc 1.38 writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock)
141 adcroft 1.1
142 edhill 1.24 #ifdef ALLOW_MNC
143     IF (useMNC .AND. snapshot_mnc .AND. writeDiag) THEN
144 edhill 1.53 IF ( writeBinaryPrec .EQ. precFloat64 ) THEN
145     pf(1:1) = 'D'
146     ELSE
147     pf(1:1) = 'R'
148     ENDIF
149 edhill 1.25 IF ((bi .EQ. 1).AND.(bj .EQ. 1).AND.(k .EQ. 1)) THEN
150     CALL MNC_CW_SET_UDIM('mom_vi', -1, myThid)
151 edhill 1.39 CALL MNC_CW_RL_W_S('D','mom_vi',0,0,'T',myTime,myThid)
152 edhill 1.25 CALL MNC_CW_SET_UDIM('mom_vi', 0, myThid)
153 edhill 1.39 CALL MNC_CW_I_W_S('I','mom_vi',0,0,'iter',myIter,myThid)
154 edhill 1.25 ENDIF
155     DO i = 1,9
156     offsets(i) = 0
157     ENDDO
158     offsets(3) = k
159 jmc 1.61 c write(*,*) 'offsets = ',(offsets(i),i=1,9)
160 edhill 1.24 ENDIF
161     #endif /* ALLOW_MNC */
162    
163 jmc 1.61 C-- Initialise intermediate terms
164     DO j=1-OLy,sNy+OLy
165     DO i=1-OLx,sNx+OLx
166 jmc 1.31 vF(i,j) = 0.
167     vrF(i,j) = 0.
168 adcroft 1.1 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 jmc 1.31 guDiss(i,j)= 0.
175     gvDiss(i,j)= 0.
176 adcroft 1.1 vort3(i,j) = 0.
177 jmc 1.31 omega3(i,j)= 0.
178 jmc 1.54 KE(i,j) = 0.
179 jmc 1.61 C- need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL)
180     hDiv(i,j) = 0.
181 baylor 1.50 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 heimbach 1.8 strain(i,j) = 0. _d 0
187     tension(i,j) = 0. _d 0
188 jmc 1.60 #ifdef ALLOW_AUTODIFF_TAMC
189 heimbach 1.55 hFacZ(i,j) = 0. _d 0
190 heimbach 1.8 #endif
191 adcroft 1.1 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 jmc 1.54 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 adcroft 1.1 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 jmc 1.7 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 jmc 1.52 CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
232 adcroft 1.1
233 jmc 1.54 CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
234 adcroft 1.1
235 jmc 1.54 IF (momViscosity) THEN
236 jmc 1.57 C-- For viscous term, compute horizontal divergence, tension & strain
237 jmc 1.54 C and mask relative vorticity (free-slip case):
238 adcroft 1.1
239 heimbach 1.59 #ifdef ALLOW_AUTODIFF_TAMC
240 jmc 1.61 CADJ STORE vort3(:,:) =
241 heimbach 1.59 CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
242     #endif
243    
244 jmc 1.54 CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
245 adcroft 1.1
246 jmc 1.52 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 jmc 1.54 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 jmc 1.52 CALL MOM_CALC_VISC(
262 baylor 1.50 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 adcroft 1.1 C Calculate del^2 u and del^2 v for bi-harmonic term
269 baylor 1.50 IF (biharmonic) THEN
270 adcroft 1.2 CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
271     O del2u,del2v,
272     & myThid)
273 jmc 1.48 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 adcroft 1.2 ENDIF
277 baylor 1.47
278 jmc 1.54 C- Strain diagnostics:
279     IF ( writeDiag ) THEN
280     IF (snapshot_mdsio) THEN
281     CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)
282     ENDIF
283     #ifdef ALLOW_MNC
284     IF (useMNC .AND. snapshot_mnc) THEN
285     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strain,
286     & offsets, myThid)
287     ENDIF
288     #endif /* ALLOW_MNC */
289     ENDIF
290     #ifdef ALLOW_DIAGNOSTICS
291     IF ( useDiagnostics ) THEN
292     CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
293     ENDIF
294     #endif /* ALLOW_DIAGNOSTICS */
295    
296     C--- Calculate dissipation terms for U and V equations
297 baylor 1.47
298     C in terms of tension and strain
299     IF (useStrainTensionVisc) THEN
300 jmc 1.54 C mask strain as if free-slip since side-drag is computed separately
301     DO j=1-Oly,sNy+Oly
302     DO i=1-Olx,sNx+Olx
303     IF ( hFacZ(i,j).EQ.0. ) strain(i,j) = 0. _d 0
304     ENDDO
305     ENDDO
306 baylor 1.47 CALL MOM_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
307     I hFacZ,
308 baylor 1.50 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
309     I harmonic,biharmonic,useVariableViscosity,
310 baylor 1.47 O guDiss,gvDiss,
311     I myThid)
312     ELSE
313 adcroft 1.2 C in terms of vorticity and divergence
314 baylor 1.47 CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
315     I hFacZ,dStar,zStar,
316 baylor 1.50 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
317     I harmonic,biharmonic,useVariableViscosity,
318 jmc 1.31 O guDiss,gvDiss,
319 jmc 1.57 & myThid)
320 adcroft 1.3 ENDIF
321 jmc 1.54 C-- if (momViscosity) end of block.
322 adcroft 1.1 ENDIF
323    
324 jmc 1.7 C- Return to standard hfacZ (min-4) and mask vort3 accordingly:
325     c CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
326    
327 jmc 1.54 C--- Other dissipation terms in Zonal momentum equation
328 adcroft 1.1
329     C-- Vertical flux (fVer is at upper face of "u" cell)
330    
331     C Eddy component of vertical flux (interior component only) -> vrF
332 jmc 1.31 IF (momViscosity.AND..NOT.implicitViscosity) THEN
333 jmc 1.44 CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,vrF,myThid)
334 adcroft 1.1
335     C Combine fluxes
336 jmc 1.31 DO j=jMin,jMax
337     DO i=iMin,iMax
338     fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)
339     ENDDO
340 adcroft 1.1 ENDDO
341    
342 jmc 1.31 C-- Tendency is minus divergence of the fluxes
343     DO j=2-Oly,sNy+Oly-1
344     DO i=2-Olx,sNx+Olx-1
345     guDiss(i,j) = guDiss(i,j)
346 adcroft 1.1 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
347     & *recip_rAw(i,j,bi,bj)
348     & *(
349 jmc 1.42 & fVerU(i,j,kDown) - fVerU(i,j,kUp)
350     & )*rkSign
351 jmc 1.31 ENDDO
352 adcroft 1.1 ENDDO
353 jmc 1.31 ENDIF
354 adcroft 1.1
355 jmc 1.57 C-- No-slip and drag BCs appear as body forces in cell abutting topography
356 adcroft 1.1 IF (momViscosity.AND.no_slip_sides) THEN
357     C- No-slip BCs impose a drag at walls...
358 baylor 1.50 CALL MOM_U_SIDEDRAG(
359     I bi,bj,k,
360     I uFld, del2u, hFacZ,
361     I viscAh_Z,viscA4_Z,
362     I harmonic,biharmonic,useVariableViscosity,
363     O vF,
364     I myThid)
365 adcroft 1.1 DO j=jMin,jMax
366     DO i=iMin,iMax
367 jmc 1.31 guDiss(i,j) = guDiss(i,j)+vF(i,j)
368 adcroft 1.1 ENDDO
369     ENDDO
370     ENDIF
371     C- No-slip BCs impose a drag at bottom
372     IF (momViscosity.AND.bottomDragTerms) THEN
373     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
374     DO j=jMin,jMax
375     DO i=iMin,iMax
376 jmc 1.31 guDiss(i,j) = guDiss(i,j)+vF(i,j)
377 adcroft 1.1 ENDDO
378     ENDDO
379     ENDIF
380 mlosch 1.56 #ifdef ALLOW_SHELFICE
381     IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
382     CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
383     DO j=jMin,jMax
384     DO i=iMin,iMax
385     guDiss(i,j) = guDiss(i,j) + vF(i,j)
386     ENDDO
387     ENDDO
388     ENDIF
389     #endif /* ALLOW_SHELFICE */
390    
391 adcroft 1.1
392 jmc 1.54 C--- Other dissipation terms in Meridional momentum equation
393 adcroft 1.1
394     C-- Vertical flux (fVer is at upper face of "v" cell)
395    
396     C Eddy component of vertical flux (interior component only) -> vrF
397 jmc 1.31 IF (momViscosity.AND..NOT.implicitViscosity) THEN
398 jmc 1.44 CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,vrF,myThid)
399 adcroft 1.1
400     C Combine fluxes -> fVerV
401 jmc 1.31 DO j=jMin,jMax
402     DO i=iMin,iMax
403     fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)
404     ENDDO
405 adcroft 1.1 ENDDO
406    
407 jmc 1.31 C-- Tendency is minus divergence of the fluxes
408     DO j=jMin,jMax
409     DO i=iMin,iMax
410     gvDiss(i,j) = gvDiss(i,j)
411 adcroft 1.1 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
412     & *recip_rAs(i,j,bi,bj)
413     & *(
414 jmc 1.42 & fVerV(i,j,kDown) - fVerV(i,j,kUp)
415     & )*rkSign
416 jmc 1.31 ENDDO
417 adcroft 1.1 ENDDO
418 jmc 1.31 ENDIF
419 adcroft 1.1
420 jmc 1.57 C-- No-slip and drag BCs appear as body forces in cell abutting topography
421 adcroft 1.1 IF (momViscosity.AND.no_slip_sides) THEN
422     C- No-slip BCs impose a drag at walls...
423 baylor 1.50 CALL MOM_V_SIDEDRAG(
424     I bi,bj,k,
425     I vFld, del2v, hFacZ,
426     I viscAh_Z,viscA4_Z,
427     I harmonic,biharmonic,useVariableViscosity,
428     O vF,
429     I myThid)
430 adcroft 1.1 DO j=jMin,jMax
431     DO i=iMin,iMax
432 jmc 1.31 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
433 adcroft 1.1 ENDDO
434     ENDDO
435     ENDIF
436     C- No-slip BCs impose a drag at bottom
437     IF (momViscosity.AND.bottomDragTerms) THEN
438     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
439     DO j=jMin,jMax
440     DO i=iMin,iMax
441 jmc 1.31 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
442 adcroft 1.1 ENDDO
443     ENDDO
444     ENDIF
445 mlosch 1.56 #ifdef ALLOW_SHELFICE
446     IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
447     CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
448     DO j=jMin,jMax
449     DO i=iMin,iMax
450     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
451     ENDDO
452     ENDDO
453     ENDIF
454     #endif /* ALLOW_SHELFICE */
455    
456 adcroft 1.1
457 jmc 1.54 C- Vorticity diagnostics:
458     IF ( writeDiag ) THEN
459     IF (snapshot_mdsio) THEN
460     CALL WRITE_LOCAL_RL('Z3','I10',1,vort3, bi,bj,k,myIter,myThid)
461     ENDIF
462     #ifdef ALLOW_MNC
463     IF (useMNC .AND. snapshot_mnc) THEN
464     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3,
465     & offsets, myThid)
466     ENDIF
467     #endif /* ALLOW_MNC */
468     ENDIF
469     #ifdef ALLOW_DIAGNOSTICS
470     IF ( useDiagnostics ) THEN
471     CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
472     ENDIF
473     #endif /* ALLOW_DIAGNOSTICS */
474    
475     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
476    
477     C--- Prepare for Advection & Coriolis terms:
478     C- Mask relative vorticity and calculate absolute vorticity
479     DO j=1-Oly,sNy+Oly
480     DO i=1-Olx,sNx+Olx
481     IF ( hFacZ(i,j).EQ.0. ) vort3(i,j) = 0.
482     ENDDO
483     ENDDO
484     IF (useAbsVorticity)
485     & CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
486 adcroft 1.1
487 jmc 1.5 C-- Horizontal Coriolis terms
488 jmc 1.37 c IF (useCoriolis .AND. .NOT.useCDscheme
489     c & .AND. .NOT. useAbsVorticity) THEN
490     C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
491 jmc 1.46 IF ( useCoriolis .AND.
492 jmc 1.37 & .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
493     & ) THEN
494     IF (useAbsVorticity) THEN
495     CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
496     & uCf,myThid)
497     CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
498     & vCf,myThid)
499     ELSE
500     CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
501     & uCf,vCf,myThid)
502     ENDIF
503 jmc 1.5 DO j=jMin,jMax
504     DO i=iMin,iMax
505 jmc 1.43 gU(i,j,k,bi,bj) = uCf(i,j)
506     gV(i,j,k,bi,bj) = vCf(i,j)
507 jmc 1.5 ENDDO
508 adcroft 1.1 ENDDO
509 jmc 1.15 IF ( writeDiag ) THEN
510 edhill 1.24 IF (snapshot_mdsio) THEN
511     CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
512     CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
513     ENDIF
514     #ifdef ALLOW_MNC
515     IF (useMNC .AND. snapshot_mnc) THEN
516 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
517 edhill 1.25 & offsets, myThid)
518 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
519 edhill 1.25 & offsets, myThid)
520 edhill 1.24 ENDIF
521     #endif /* ALLOW_MNC */
522 jmc 1.15 ENDIF
523 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
524     IF ( useDiagnostics ) THEN
525     CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
526     CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
527     ENDIF
528     #endif /* ALLOW_DIAGNOSTICS */
529 jmc 1.31 ELSE
530     DO j=jMin,jMax
531     DO i=iMin,iMax
532 jmc 1.43 gU(i,j,k,bi,bj) = 0. _d 0
533     gV(i,j,k,bi,bj) = 0. _d 0
534 jmc 1.31 ENDDO
535     ENDDO
536 jmc 1.5 ENDIF
537 adcroft 1.1
538 jmc 1.5 IF (momAdvection) THEN
539 jmc 1.41 C-- Horizontal advection of relative (or absolute) vorticity
540 jmc 1.62 IF ( (highOrderVorticity.OR.upwindVorticity)
541     & .AND.useAbsVorticity ) THEN
542 jmc 1.41 CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,
543 adcroft 1.20 & uCf,myThid)
544 jmc 1.62 ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
545 jmc 1.41 CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,
546     & uCf,myThid)
547 jmc 1.62 ELSEIF ( useAbsVorticity ) THEN
548 jmc 1.41 CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
549 jmc 1.40 & uCf,myThid)
550 adcroft 1.20 ELSE
551 jmc 1.41 CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,
552 adcroft 1.20 & uCf,myThid)
553     ENDIF
554 jmc 1.5 DO j=jMin,jMax
555     DO i=iMin,iMax
556     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
557     ENDDO
558 adcroft 1.1 ENDDO
559 jmc 1.62 IF ( (highOrderVorticity.OR.upwindVorticity)
560     & .AND.useAbsVorticity ) THEN
561 jmc 1.41 CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ,
562 adcroft 1.20 & vCf,myThid)
563 jmc 1.62 ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN
564 jmc 1.41 CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ,
565     & vCf,myThid)
566 jmc 1.62 ELSEIF ( useAbsVorticity ) THEN
567 jmc 1.41 CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
568 jmc 1.40 & vCf,myThid)
569 adcroft 1.20 ELSE
570 jmc 1.41 CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,
571 adcroft 1.20 & vCf,myThid)
572     ENDIF
573 jmc 1.5 DO j=jMin,jMax
574     DO i=iMin,iMax
575     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
576     ENDDO
577 adcroft 1.1 ENDDO
578    
579 jmc 1.15 IF ( writeDiag ) THEN
580 edhill 1.24 IF (snapshot_mdsio) THEN
581     CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
582     CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
583     ENDIF
584     #ifdef ALLOW_MNC
585     IF (useMNC .AND. snapshot_mnc) THEN
586 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
587 edhill 1.25 & offsets, myThid)
588 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
589 edhill 1.25 & offsets, myThid)
590 edhill 1.24 ENDIF
591     #endif /* ALLOW_MNC */
592 jmc 1.15 ENDIF
593 edhill 1.24
594 jmc 1.7 #ifdef ALLOW_TIMEAVE
595     IF (taveFreq.GT.0.) THEN
596     CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
597     & Nr, k, bi, bj, myThid)
598     CALL TIMEAVE_CUMUL_1K1T(vZetatave,uCf,deltaTClock,
599     & Nr, k, bi, bj, myThid)
600     ENDIF
601 dimitri 1.13 #endif /* ALLOW_TIMEAVE */
602 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
603     IF ( useDiagnostics ) THEN
604     CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
605     CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
606     ENDIF
607     #endif /* ALLOW_DIAGNOSTICS */
608 jmc 1.7
609 jmc 1.5 C-- Vertical shear terms (-w*du/dr & -w*dv/dr)
610 jmc 1.12 IF ( .NOT. momImplVertAdv ) THEN
611     CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)
612     DO j=jMin,jMax
613     DO i=iMin,iMax
614     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
615     ENDDO
616 jmc 1.5 ENDDO
617 jmc 1.12 CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)
618     DO j=jMin,jMax
619     DO i=iMin,iMax
620     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
621     ENDDO
622 jmc 1.5 ENDDO
623 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
624     IF ( useDiagnostics ) THEN
625     CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
626     CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
627     ENDIF
628     #endif /* ALLOW_DIAGNOSTICS */
629 jmc 1.12 ENDIF
630 adcroft 1.1
631     C-- Bernoulli term
632 jmc 1.5 CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)
633     DO j=jMin,jMax
634     DO i=iMin,iMax
635     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
636     ENDDO
637     ENDDO
638     CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)
639     DO j=jMin,jMax
640     DO i=iMin,iMax
641     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
642     ENDDO
643 adcroft 1.1 ENDDO
644 jmc 1.15 IF ( writeDiag ) THEN
645 edhill 1.24 IF (snapshot_mdsio) THEN
646     CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
647     CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
648     ENDIF
649     #ifdef ALLOW_MNC
650     IF (useMNC .AND. snapshot_mnc) THEN
651 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
652 edhill 1.25 & offsets, myThid)
653 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
654 edhill 1.25 & offsets, myThid)
655 jmc 1.54 ENDIF
656 edhill 1.24 #endif /* ALLOW_MNC */
657 jmc 1.15 ENDIF
658    
659 jmc 1.5 C-- end if momAdvection
660     ENDIF
661    
662 jmc 1.57 C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
663 jmc 1.58 IF ( use3dCoriolis ) THEN
664 jmc 1.57 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
665     DO j=jMin,jMax
666     DO i=iMin,iMax
667     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
668     ENDDO
669     ENDDO
670     IF ( usingCurvilinearGrid ) THEN
671     C- presently, non zero angleSinC array only supported with Curvilinear-Grid
672     CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
673     DO j=jMin,jMax
674     DO i=iMin,iMax
675     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
676     ENDDO
677     ENDDO
678     ENDIF
679     ENDIF
680    
681     C-- Non-Hydrostatic (spherical) metric terms
682     IF ( useNHMTerms ) THEN
683     CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
684     DO j=jMin,jMax
685     DO i=iMin,iMax
686     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
687     ENDDO
688     ENDDO
689     CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
690     DO j=jMin,jMax
691     DO i=iMin,iMax
692     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
693     ENDDO
694     ENDDO
695     ENDIF
696 jmc 1.54
697 jmc 1.5 C-- Set du/dt & dv/dt on boundaries to zero
698 adcroft 1.1 DO j=jMin,jMax
699     DO i=iMin,iMax
700 jmc 1.5 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
701     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
702 adcroft 1.1 ENDDO
703     ENDDO
704 jmc 1.5
705 jmc 1.22 #ifdef ALLOW_DEBUG
706     IF ( debugLevel .GE. debLevB
707     & .AND. k.EQ.4 .AND. myIter.EQ.nIter0
708     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
709     & .AND. useCubedSphereExchange ) THEN
710 jmc 1.23 CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
711 jmc 1.31 & guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
712 jmc 1.22 ENDIF
713     #endif /* ALLOW_DEBUG */
714 adcroft 1.2
715 jmc 1.15 IF ( writeDiag ) THEN
716 edhill 1.24 IF (snapshot_mdsio) THEN
717 jmc 1.54 CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
718     CALL WRITE_LOCAL_RL('KE','I10',1,KE, bi,bj,k,myIter,myThid)
719     CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv, bi,bj,k,myIter,myThid)
720     CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
721     CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
722     CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
723 edhill 1.24 ENDIF
724     #ifdef ALLOW_MNC
725     IF (useMNC .AND. snapshot_mnc) THEN
726 jmc 1.54 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
727     & offsets, myThid)
728     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
729     & offsets, myThid)
730     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
731 edhill 1.25 & offsets, myThid)
732 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
733 edhill 1.25 & offsets, myThid)
734 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
735 edhill 1.25 & offsets, myThid)
736 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
737 edhill 1.25 & offsets, myThid)
738 edhill 1.24 ENDIF
739     #endif /* ALLOW_MNC */
740 adcroft 1.1 ENDIF
741 jmc 1.41
742 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
743     IF ( useDiagnostics ) THEN
744 jmc 1.54 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
745 jmc 1.46 IF (momViscosity) THEN
746 jmc 1.54 CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
747 jmc 1.52 CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
748 jmc 1.54 CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)
749     CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)
750 jmc 1.46 ENDIF
751 jmc 1.54 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
752     & 'Um_Advec',k,1,2,bi,bj,myThid)
753     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
754     & 'Vm_Advec',k,1,2,bi,bj,myThid)
755 jmc 1.46 ENDIF
756     #endif /* ALLOW_DIAGNOSTICS */
757    
758 edhill 1.11 #endif /* ALLOW_MOM_VECINV */
759 adcroft 1.1
760     RETURN
761     END

  ViewVC Help
Powered by ViewVC 1.1.22