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