/[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.60 - (hide annotations) (download)
Thu Nov 23 00:45:21 2006 UTC (17 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58t_post, checkpoint58w_post, mitgcm_mapl_00, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59i, checkpoint59h, checkpoint58v_post, checkpoint58x_post, checkpoint59j, checkpoint58u_post, checkpoint58s_post
Changes since 1.59: +3 -2 lines
initialise tension & strain (Pb reported by Martin).

1 jmc 1.60 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vecinv.F,v 1.59 2006/07/18 03:23:30 heimbach 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     C write(*,*) 'offsets = ',(offsets(i),i=1,9)
160 edhill 1.24 ENDIF
161     #endif /* ALLOW_MNC */
162    
163 adcroft 1.1 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.60 c hDiv(i,j) = 0.
180 baylor 1.50 viscAh_Z(i,j) = 0.
181     viscAh_D(i,j) = 0.
182     viscA4_Z(i,j) = 0.
183     viscA4_D(i,j) = 0.
184    
185 heimbach 1.8 strain(i,j) = 0. _d 0
186     tension(i,j) = 0. _d 0
187 jmc 1.60 #ifdef ALLOW_AUTODIFF_TAMC
188 heimbach 1.55 hFacZ(i,j) = 0. _d 0
189 heimbach 1.8 #endif
190 adcroft 1.1 ENDDO
191     ENDDO
192    
193     C-- Term by term tracer parmeters
194     C o U momentum equation
195     ArDudrFac = vfFacMom*1.
196     C o V momentum equation
197     ArDvdrFac = vfFacMom*1.
198    
199 jmc 1.54 C note: using standard stencil (no mask) results in under-estimating
200     C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
201     IF ( no_slip_sides ) THEN
202     sideMaskFac = sideDragFactor
203     ELSE
204     sideMaskFac = 0. _d 0
205     ENDIF
206    
207 adcroft 1.1 IF ( no_slip_bottom
208     & .OR. bottomDragQuadratic.NE.0.
209     & .OR. bottomDragLinear.NE.0.) THEN
210     bottomDragTerms=.TRUE.
211     ELSE
212     bottomDragTerms=.FALSE.
213     ENDIF
214    
215     C-- Calculate open water fraction at vorticity points
216     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
217    
218     C Make local copies of horizontal flow field
219     DO j=1-OLy,sNy+OLy
220     DO i=1-OLx,sNx+OLx
221     uFld(i,j) = uVel(i,j,k,bi,bj)
222     vFld(i,j) = vVel(i,j,k,bi,bj)
223     ENDDO
224     ENDDO
225    
226 jmc 1.7 C note (jmc) : Dissipation and Vort3 advection do not necesary
227     C use the same maskZ (and hFacZ) => needs 2 call(s)
228     c CALL MOM_VI_HFACZ_DISS(bi,bj,k,hFacZ,r_hFacZ,myThid)
229    
230 jmc 1.52 CALL MOM_CALC_KE(bi,bj,k,selectKEscheme,uFld,vFld,KE,myThid)
231 adcroft 1.1
232 jmc 1.54 CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
233 adcroft 1.1
234 jmc 1.54 IF (momViscosity) THEN
235 jmc 1.57 C-- For viscous term, compute horizontal divergence, tension & strain
236 jmc 1.54 C and mask relative vorticity (free-slip case):
237 adcroft 1.1
238 heimbach 1.59 #ifdef ALLOW_AUTODIFF_TAMC
239     CADJ STORE vort3(:,:) =
240     CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
241     #endif
242    
243 jmc 1.54 CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
244 adcroft 1.1
245 jmc 1.52 CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
246    
247     CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
248    
249 jmc 1.54 C- account for no-slip / free-slip BC:
250     DO j=1-Oly,sNy+Oly
251     DO i=1-Olx,sNx+Olx
252     IF ( hFacZ(i,j).EQ.0. ) THEN
253     vort3(i,j) = sideMaskFac*vort3(i,j)
254     strain(i,j) = sideMaskFac*strain(i,j)
255     ENDIF
256     ENDDO
257     ENDDO
258    
259     C-- Calculate Viscosities
260 jmc 1.52 CALL MOM_CALC_VISC(
261 baylor 1.50 I bi,bj,k,
262     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
263     O harmonic,biharmonic,useVariableViscosity,
264     I hDiv,vort3,tension,strain,KE,hfacZ,
265     I myThid)
266    
267 adcroft 1.1 C Calculate del^2 u and del^2 v for bi-harmonic term
268 baylor 1.50 IF (biharmonic) THEN
269 adcroft 1.2 CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
270     O del2u,del2v,
271     & myThid)
272 jmc 1.48 CALL MOM_CALC_HDIV(bi,bj,k,2,del2u,del2v,dStar,myThid)
273     CALL MOM_CALC_RELVORT3(bi,bj,k,
274     & del2u,del2v,hFacZ,zStar,myThid)
275 adcroft 1.2 ENDIF
276 baylor 1.47
277 jmc 1.54 C- Strain diagnostics:
278     IF ( writeDiag ) THEN
279     IF (snapshot_mdsio) THEN
280     CALL WRITE_LOCAL_RL('Ds','I10',1,strain,bi,bj,k,myIter,myThid)
281     ENDIF
282     #ifdef ALLOW_MNC
283     IF (useMNC .AND. snapshot_mnc) THEN
284     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Ds',strain,
285     & offsets, myThid)
286     ENDIF
287     #endif /* ALLOW_MNC */
288     ENDIF
289     #ifdef ALLOW_DIAGNOSTICS
290     IF ( useDiagnostics ) THEN
291     CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
292     ENDIF
293     #endif /* ALLOW_DIAGNOSTICS */
294    
295     C--- Calculate dissipation terms for U and V equations
296 baylor 1.47
297     C in terms of tension and strain
298     IF (useStrainTensionVisc) THEN
299 jmc 1.54 C mask strain as if free-slip since side-drag is computed separately
300     DO j=1-Oly,sNy+Oly
301     DO i=1-Olx,sNx+Olx
302     IF ( hFacZ(i,j).EQ.0. ) strain(i,j) = 0. _d 0
303     ENDDO
304     ENDDO
305 baylor 1.47 CALL MOM_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
306     I hFacZ,
307 baylor 1.50 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
308     I harmonic,biharmonic,useVariableViscosity,
309 baylor 1.47 O guDiss,gvDiss,
310     I myThid)
311     ELSE
312 adcroft 1.2 C in terms of vorticity and divergence
313 baylor 1.47 CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,tension,strain,KE,
314     I hFacZ,dStar,zStar,
315 baylor 1.50 I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
316     I harmonic,biharmonic,useVariableViscosity,
317 jmc 1.31 O guDiss,gvDiss,
318 jmc 1.57 & myThid)
319 adcroft 1.3 ENDIF
320 jmc 1.54 C-- if (momViscosity) end of block.
321 adcroft 1.1 ENDIF
322    
323 jmc 1.7 C- Return to standard hfacZ (min-4) and mask vort3 accordingly:
324     c CALL MOM_VI_MASK_VORT3(bi,bj,k,hFacZ,r_hFacZ,vort3,myThid)
325    
326 jmc 1.54 C--- Other dissipation terms in Zonal momentum equation
327 adcroft 1.1
328     C-- Vertical flux (fVer is at upper face of "u" cell)
329    
330     C Eddy component of vertical flux (interior component only) -> vrF
331 jmc 1.31 IF (momViscosity.AND..NOT.implicitViscosity) THEN
332 jmc 1.44 CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,vrF,myThid)
333 adcroft 1.1
334     C Combine fluxes
335 jmc 1.31 DO j=jMin,jMax
336     DO i=iMin,iMax
337     fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)
338     ENDDO
339 adcroft 1.1 ENDDO
340    
341 jmc 1.31 C-- Tendency is minus divergence of the fluxes
342     DO j=2-Oly,sNy+Oly-1
343     DO i=2-Olx,sNx+Olx-1
344     guDiss(i,j) = guDiss(i,j)
345 adcroft 1.1 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
346     & *recip_rAw(i,j,bi,bj)
347     & *(
348 jmc 1.42 & fVerU(i,j,kDown) - fVerU(i,j,kUp)
349     & )*rkSign
350 jmc 1.31 ENDDO
351 adcroft 1.1 ENDDO
352 jmc 1.31 ENDIF
353 adcroft 1.1
354 jmc 1.57 C-- No-slip and drag BCs appear as body forces in cell abutting topography
355 adcroft 1.1 IF (momViscosity.AND.no_slip_sides) THEN
356     C- No-slip BCs impose a drag at walls...
357 baylor 1.50 CALL MOM_U_SIDEDRAG(
358     I bi,bj,k,
359     I uFld, del2u, hFacZ,
360     I viscAh_Z,viscA4_Z,
361     I harmonic,biharmonic,useVariableViscosity,
362     O vF,
363     I myThid)
364 adcroft 1.1 DO j=jMin,jMax
365     DO i=iMin,iMax
366 jmc 1.31 guDiss(i,j) = guDiss(i,j)+vF(i,j)
367 adcroft 1.1 ENDDO
368     ENDDO
369     ENDIF
370     C- No-slip BCs impose a drag at bottom
371     IF (momViscosity.AND.bottomDragTerms) THEN
372     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
373     DO j=jMin,jMax
374     DO i=iMin,iMax
375 jmc 1.31 guDiss(i,j) = guDiss(i,j)+vF(i,j)
376 adcroft 1.1 ENDDO
377     ENDDO
378     ENDIF
379 mlosch 1.56 #ifdef ALLOW_SHELFICE
380     IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
381     CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
382     DO j=jMin,jMax
383     DO i=iMin,iMax
384     guDiss(i,j) = guDiss(i,j) + vF(i,j)
385     ENDDO
386     ENDDO
387     ENDIF
388     #endif /* ALLOW_SHELFICE */
389    
390 adcroft 1.1
391 jmc 1.54 C--- Other dissipation terms in Meridional momentum equation
392 adcroft 1.1
393     C-- Vertical flux (fVer is at upper face of "v" cell)
394    
395     C Eddy component of vertical flux (interior component only) -> vrF
396 jmc 1.31 IF (momViscosity.AND..NOT.implicitViscosity) THEN
397 jmc 1.44 CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,vrF,myThid)
398 adcroft 1.1
399     C Combine fluxes -> fVerV
400 jmc 1.31 DO j=jMin,jMax
401     DO i=iMin,iMax
402     fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)
403     ENDDO
404 adcroft 1.1 ENDDO
405    
406 jmc 1.31 C-- Tendency is minus divergence of the fluxes
407     DO j=jMin,jMax
408     DO i=iMin,iMax
409     gvDiss(i,j) = gvDiss(i,j)
410 adcroft 1.1 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
411     & *recip_rAs(i,j,bi,bj)
412     & *(
413 jmc 1.42 & fVerV(i,j,kDown) - fVerV(i,j,kUp)
414     & )*rkSign
415 jmc 1.31 ENDDO
416 adcroft 1.1 ENDDO
417 jmc 1.31 ENDIF
418 adcroft 1.1
419 jmc 1.57 C-- No-slip and drag BCs appear as body forces in cell abutting topography
420 adcroft 1.1 IF (momViscosity.AND.no_slip_sides) THEN
421     C- No-slip BCs impose a drag at walls...
422 baylor 1.50 CALL MOM_V_SIDEDRAG(
423     I bi,bj,k,
424     I vFld, del2v, hFacZ,
425     I viscAh_Z,viscA4_Z,
426     I harmonic,biharmonic,useVariableViscosity,
427     O vF,
428     I myThid)
429 adcroft 1.1 DO j=jMin,jMax
430     DO i=iMin,iMax
431 jmc 1.31 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
432 adcroft 1.1 ENDDO
433     ENDDO
434     ENDIF
435     C- No-slip BCs impose a drag at bottom
436     IF (momViscosity.AND.bottomDragTerms) THEN
437     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
438     DO j=jMin,jMax
439     DO i=iMin,iMax
440 jmc 1.31 gvDiss(i,j) = gvDiss(i,j)+vF(i,j)
441 adcroft 1.1 ENDDO
442     ENDDO
443     ENDIF
444 mlosch 1.56 #ifdef ALLOW_SHELFICE
445     IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) THEN
446     CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
447     DO j=jMin,jMax
448     DO i=iMin,iMax
449     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
450     ENDDO
451     ENDDO
452     ENDIF
453     #endif /* ALLOW_SHELFICE */
454    
455 adcroft 1.1
456 jmc 1.54 C- Vorticity diagnostics:
457     IF ( writeDiag ) THEN
458     IF (snapshot_mdsio) THEN
459     CALL WRITE_LOCAL_RL('Z3','I10',1,vort3, bi,bj,k,myIter,myThid)
460     ENDIF
461     #ifdef ALLOW_MNC
462     IF (useMNC .AND. snapshot_mnc) THEN
463     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Z3',vort3,
464     & offsets, myThid)
465     ENDIF
466     #endif /* ALLOW_MNC */
467     ENDIF
468     #ifdef ALLOW_DIAGNOSTICS
469     IF ( useDiagnostics ) THEN
470     CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
471     ENDIF
472     #endif /* ALLOW_DIAGNOSTICS */
473    
474     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
475    
476     C--- Prepare for Advection & Coriolis terms:
477     C- Mask relative vorticity and calculate absolute vorticity
478     DO j=1-Oly,sNy+Oly
479     DO i=1-Olx,sNx+Olx
480     IF ( hFacZ(i,j).EQ.0. ) vort3(i,j) = 0.
481     ENDDO
482     ENDDO
483     IF (useAbsVorticity)
484     & CALL MOM_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
485 adcroft 1.1
486 jmc 1.5 C-- Horizontal Coriolis terms
487 jmc 1.37 c IF (useCoriolis .AND. .NOT.useCDscheme
488     c & .AND. .NOT. useAbsVorticity) THEN
489     C- jmc: change it to keep the Coriolis terms when useAbsVorticity=T & momAdvection=F
490 jmc 1.46 IF ( useCoriolis .AND.
491 jmc 1.37 & .NOT.( useCDscheme .OR. useAbsVorticity.AND.momAdvection )
492     & ) THEN
493     IF (useAbsVorticity) THEN
494     CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
495     & uCf,myThid)
496     CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
497     & vCf,myThid)
498     ELSE
499     CALL MOM_VI_CORIOLIS(bi,bj,k,uFld,vFld,hFacZ,r_hFacZ,
500     & uCf,vCf,myThid)
501     ENDIF
502 jmc 1.5 DO j=jMin,jMax
503     DO i=iMin,iMax
504 jmc 1.43 gU(i,j,k,bi,bj) = uCf(i,j)
505     gV(i,j,k,bi,bj) = vCf(i,j)
506 jmc 1.5 ENDDO
507 adcroft 1.1 ENDDO
508 jmc 1.15 IF ( writeDiag ) THEN
509 edhill 1.24 IF (snapshot_mdsio) THEN
510     CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
511     CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
512     ENDIF
513     #ifdef ALLOW_MNC
514     IF (useMNC .AND. snapshot_mnc) THEN
515 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fV', uCf,
516 edhill 1.25 & offsets, myThid)
517 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'fU', vCf,
518 edhill 1.25 & offsets, myThid)
519 edhill 1.24 ENDIF
520     #endif /* ALLOW_MNC */
521 jmc 1.15 ENDIF
522 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
523     IF ( useDiagnostics ) THEN
524     CALL DIAGNOSTICS_FILL(uCf,'Um_Cori ',k,1,2,bi,bj,myThid)
525     CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid)
526     ENDIF
527     #endif /* ALLOW_DIAGNOSTICS */
528 jmc 1.31 ELSE
529     DO j=jMin,jMax
530     DO i=iMin,iMax
531 jmc 1.43 gU(i,j,k,bi,bj) = 0. _d 0
532     gV(i,j,k,bi,bj) = 0. _d 0
533 jmc 1.31 ENDDO
534     ENDDO
535 jmc 1.5 ENDIF
536 adcroft 1.1
537 jmc 1.5 IF (momAdvection) THEN
538 jmc 1.41 C-- Horizontal advection of relative (or absolute) vorticity
539     IF (highOrderVorticity.AND.useAbsVorticity) THEN
540     CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ,
541 adcroft 1.20 & uCf,myThid)
542 jmc 1.40 ELSEIF (highOrderVorticity) THEN
543 jmc 1.41 CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ,
544     & uCf,myThid)
545     ELSEIF (useAbsVorticity) THEN
546     CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ,
547 jmc 1.40 & uCf,myThid)
548 adcroft 1.20 ELSE
549 jmc 1.41 CALL MOM_VI_U_CORIOLIS(bi,bj,k,vFld,vort3, hFacZ,r_hFacZ,
550 adcroft 1.20 & uCf,myThid)
551     ENDIF
552 jmc 1.5 DO j=jMin,jMax
553     DO i=iMin,iMax
554     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
555     ENDDO
556 adcroft 1.1 ENDDO
557 jmc 1.41 IF (highOrderVorticity.AND.useAbsVorticity) THEN
558     CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ,
559 adcroft 1.20 & vCf,myThid)
560 jmc 1.40 ELSEIF (highOrderVorticity) THEN
561 jmc 1.41 CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ,
562     & vCf,myThid)
563     ELSEIF (useAbsVorticity) THEN
564     CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ,
565 jmc 1.40 & vCf,myThid)
566 adcroft 1.20 ELSE
567 jmc 1.41 CALL MOM_VI_V_CORIOLIS(bi,bj,k,uFld,vort3, hFacZ,r_hFacZ,
568 adcroft 1.20 & vCf,myThid)
569     ENDIF
570 jmc 1.5 DO j=jMin,jMax
571     DO i=iMin,iMax
572     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
573     ENDDO
574 adcroft 1.1 ENDDO
575    
576 jmc 1.15 IF ( writeDiag ) THEN
577 edhill 1.24 IF (snapshot_mdsio) THEN
578     CALL WRITE_LOCAL_RL('zV','I10',1,uCf,bi,bj,k,myIter,myThid)
579     CALL WRITE_LOCAL_RL('zU','I10',1,vCf,bi,bj,k,myIter,myThid)
580     ENDIF
581     #ifdef ALLOW_MNC
582     IF (useMNC .AND. snapshot_mnc) THEN
583 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zV', uCf,
584 edhill 1.25 & offsets, myThid)
585 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'zU', vCf,
586 edhill 1.25 & offsets, myThid)
587 edhill 1.24 ENDIF
588     #endif /* ALLOW_MNC */
589 jmc 1.15 ENDIF
590 edhill 1.24
591 jmc 1.7 #ifdef ALLOW_TIMEAVE
592     IF (taveFreq.GT.0.) THEN
593     CALL TIMEAVE_CUMUL_1K1T(uZetatave,vCf,deltaTClock,
594     & Nr, k, bi, bj, myThid)
595     CALL TIMEAVE_CUMUL_1K1T(vZetatave,uCf,deltaTClock,
596     & Nr, k, bi, bj, myThid)
597     ENDIF
598 dimitri 1.13 #endif /* ALLOW_TIMEAVE */
599 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
600     IF ( useDiagnostics ) THEN
601     CALL DIAGNOSTICS_FILL(uCf,'Um_AdvZ3',k,1,2,bi,bj,myThid)
602     CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvZ3',k,1,2,bi,bj,myThid)
603     ENDIF
604     #endif /* ALLOW_DIAGNOSTICS */
605 jmc 1.7
606 jmc 1.5 C-- Vertical shear terms (-w*du/dr & -w*dv/dr)
607 jmc 1.12 IF ( .NOT. momImplVertAdv ) THEN
608     CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)
609     DO j=jMin,jMax
610     DO i=iMin,iMax
611     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
612     ENDDO
613 jmc 1.5 ENDDO
614 jmc 1.12 CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)
615     DO j=jMin,jMax
616     DO i=iMin,iMax
617     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
618     ENDDO
619 jmc 1.5 ENDDO
620 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
621     IF ( useDiagnostics ) THEN
622     CALL DIAGNOSTICS_FILL(uCf,'Um_AdvRe',k,1,2,bi,bj,myThid)
623     CALL DIAGNOSTICS_FILL(vCf,'Vm_AdvRe',k,1,2,bi,bj,myThid)
624     ENDIF
625     #endif /* ALLOW_DIAGNOSTICS */
626 jmc 1.12 ENDIF
627 adcroft 1.1
628     C-- Bernoulli term
629 jmc 1.5 CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)
630     DO j=jMin,jMax
631     DO i=iMin,iMax
632     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
633     ENDDO
634     ENDDO
635     CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)
636     DO j=jMin,jMax
637     DO i=iMin,iMax
638     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
639     ENDDO
640 adcroft 1.1 ENDDO
641 jmc 1.15 IF ( writeDiag ) THEN
642 edhill 1.24 IF (snapshot_mdsio) THEN
643     CALL WRITE_LOCAL_RL('KEx','I10',1,uCf,bi,bj,k,myIter,myThid)
644     CALL WRITE_LOCAL_RL('KEy','I10',1,vCf,bi,bj,k,myIter,myThid)
645     ENDIF
646     #ifdef ALLOW_MNC
647     IF (useMNC .AND. snapshot_mnc) THEN
648 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEx', uCf,
649 edhill 1.25 & offsets, myThid)
650 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj, 'KEy', vCf,
651 edhill 1.25 & offsets, myThid)
652 jmc 1.54 ENDIF
653 edhill 1.24 #endif /* ALLOW_MNC */
654 jmc 1.15 ENDIF
655    
656 jmc 1.5 C-- end if momAdvection
657     ENDIF
658    
659 jmc 1.57 C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
660 jmc 1.58 IF ( use3dCoriolis ) THEN
661 jmc 1.57 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid)
662     DO j=jMin,jMax
663     DO i=iMin,iMax
664     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
665     ENDDO
666     ENDDO
667     IF ( usingCurvilinearGrid ) THEN
668     C- presently, non zero angleSinC array only supported with Curvilinear-Grid
669     CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,vCf,myThid)
670     DO j=jMin,jMax
671     DO i=iMin,iMax
672     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
673     ENDDO
674     ENDDO
675     ENDIF
676     ENDIF
677    
678     C-- Non-Hydrostatic (spherical) metric terms
679     IF ( useNHMTerms ) THEN
680     CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid)
681     DO j=jMin,jMax
682     DO i=iMin,iMax
683     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j)
684     ENDDO
685     ENDDO
686     CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid)
687     DO j=jMin,jMax
688     DO i=iMin,iMax
689     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j)
690     ENDDO
691     ENDDO
692     ENDIF
693 jmc 1.54
694 jmc 1.5 C-- Set du/dt & dv/dt on boundaries to zero
695 adcroft 1.1 DO j=jMin,jMax
696     DO i=iMin,iMax
697 jmc 1.5 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
698     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
699 adcroft 1.1 ENDDO
700     ENDDO
701 jmc 1.5
702 jmc 1.22 #ifdef ALLOW_DEBUG
703     IF ( debugLevel .GE. debLevB
704     & .AND. k.EQ.4 .AND. myIter.EQ.nIter0
705     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
706     & .AND. useCubedSphereExchange ) THEN
707 jmc 1.23 CALL DEBUG_CS_CORNER_UV( ' uDiss,vDiss from MOM_VECINV',
708 jmc 1.31 & guDiss,gvDiss, k, standardMessageUnit,bi,bj,myThid )
709 jmc 1.22 ENDIF
710     #endif /* ALLOW_DEBUG */
711 adcroft 1.2
712 jmc 1.15 IF ( writeDiag ) THEN
713 edhill 1.24 IF (snapshot_mdsio) THEN
714 jmc 1.54 CALL WRITE_LOCAL_RL('W3','I10',1,omega3, bi,bj,k,myIter,myThid)
715     CALL WRITE_LOCAL_RL('KE','I10',1,KE, bi,bj,k,myIter,myThid)
716     CALL WRITE_LOCAL_RL('D', 'I10',1,hDiv, bi,bj,k,myIter,myThid)
717     CALL WRITE_LOCAL_RL('Dt','I10',1,tension,bi,bj,k,myIter,myThid)
718     CALL WRITE_LOCAL_RL('Du','I10',1,guDiss, bi,bj,k,myIter,myThid)
719     CALL WRITE_LOCAL_RL('Dv','I10',1,gvDiss, bi,bj,k,myIter,myThid)
720 edhill 1.24 ENDIF
721     #ifdef ALLOW_MNC
722     IF (useMNC .AND. snapshot_mnc) THEN
723 jmc 1.54 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'W3',omega3,
724     & offsets, myThid)
725     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'KE',KE,
726     & offsets, myThid)
727     CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'D', hDiv,
728 edhill 1.25 & offsets, myThid)
729 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dt',tension,
730 edhill 1.25 & offsets, myThid)
731 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Du',guDiss,
732 edhill 1.25 & offsets, myThid)
733 edhill 1.53 CALL MNC_CW_RL_W_OFFSET(pf,'mom_vi',bi,bj,'Dv',gvDiss,
734 edhill 1.25 & offsets, myThid)
735 edhill 1.24 ENDIF
736     #endif /* ALLOW_MNC */
737 adcroft 1.1 ENDIF
738 jmc 1.41
739 jmc 1.46 #ifdef ALLOW_DIAGNOSTICS
740     IF ( useDiagnostics ) THEN
741 jmc 1.54 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
742 jmc 1.46 IF (momViscosity) THEN
743 jmc 1.54 CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
744 jmc 1.52 CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
745 jmc 1.54 CALL DIAGNOSTICS_FILL(guDiss, 'Um_Diss ',k,1,2,bi,bj,myThid)
746     CALL DIAGNOSTICS_FILL(gvDiss, 'Vm_Diss ',k,1,2,bi,bj,myThid)
747 jmc 1.46 ENDIF
748 jmc 1.54 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
749     & 'Um_Advec',k,1,2,bi,bj,myThid)
750     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
751     & 'Vm_Advec',k,1,2,bi,bj,myThid)
752 jmc 1.46 ENDIF
753     #endif /* ALLOW_DIAGNOSTICS */
754    
755 edhill 1.11 #endif /* ALLOW_MOM_VECINV */
756 adcroft 1.1
757     RETURN
758     END

  ViewVC Help
Powered by ViewVC 1.1.22