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