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