/[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.66 - (hide annotations) (download)
Sun Mar 18 22:24:01 2012 UTC (12 years, 2 months ago) by jmc
Branch: MAIN
Changes since 1.65: +49 -43 lines
separate fVer?(:,:,kUp) & fVer?(:,:,kDown) in argument list of MOM_FLUXFORM
 & MOM_VECINV subroutines (to help TAF).

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

  ViewVC Help
Powered by ViewVC 1.1.22