/[MITgcm]/MITgcm/pkg/mom_fluxform/mom_fluxform.F
ViewVC logotype

Annotation of /MITgcm/pkg/mom_fluxform/mom_fluxform.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.35 - (hide annotations) (download)
Wed May 3 23:35:11 2006 UTC (18 years, 1 month ago) by heimbach
Branch: MAIN
Changes since 1.34: +35 -1 lines
Now rstar adjoint.

1 heimbach 1.35 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.34 2006/03/30 19:49:41 jmc Exp $
2 adcroft 1.2 C $Name: $
3 adcroft 1.1
4 adcroft 1.3 CBOI
5     C !TITLE: pkg/mom\_advdiff
6     C !AUTHORS: adcroft@mit.edu
7 adcroft 1.4 C !INTRODUCTION: Flux-form Momentum Equations Package
8 adcroft 1.3 C
9     C Package "mom\_fluxform" provides methods for calculating explicit terms
10     C in the momentum equation cast in flux-form:
11     C \begin{eqnarray*}
12     C G^u & = & -\frac{1}{\rho} \partial_x \phi_h
13     C -\nabla \cdot {\bf v} u
14     C -fv
15     C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^x
16     C + \mbox{metrics}
17     C \\
18     C G^v & = & -\frac{1}{\rho} \partial_y \phi_h
19     C -\nabla \cdot {\bf v} v
20     C +fu
21     C +\frac{1}{\rho} \nabla \cdot {\bf \tau}^y
22     C + \mbox{metrics}
23     C \end{eqnarray*}
24     C where ${\bf v}=(u,v,w)$ and $\tau$, the stress tensor, includes surface
25     C stresses as well as internal viscous stresses.
26     CEOI
27    
28 edhill 1.13 #include "MOM_FLUXFORM_OPTIONS.h"
29 adcroft 1.1
30 adcroft 1.3 CBOP
31     C !ROUTINE: MOM_FLUXFORM
32    
33     C !INTERFACE: ==========================================================
34 adcroft 1.1 SUBROUTINE MOM_FLUXFORM(
35     I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
36 jmc 1.23 I KappaRU, KappaRV,
37 adcroft 1.1 U fVerU, fVerV,
38 jmc 1.23 O guDiss, gvDiss,
39     I myTime, myIter, myThid)
40 adcroft 1.3
41     C !DESCRIPTION:
42     C Calculates all the horizontal accelerations except for the implicit surface
43     C pressure gradient and implciit vertical viscosity.
44 adcroft 1.1
45 adcroft 1.3 C !USES: ===============================================================
46 adcroft 1.1 C == Global variables ==
47 adcroft 1.3 IMPLICIT NONE
48 adcroft 1.1 #include "SIZE.h"
49     #include "DYNVARS.h"
50     #include "FFIELDS.h"
51     #include "EEPARAMS.h"
52     #include "PARAMS.h"
53     #include "GRID.h"
54     #include "SURFACE.h"
55 heimbach 1.35 #ifdef ALLOW_AUTODIFF_TAMC
56     # include "tamc.h"
57     # include "tamc_keys.h"
58     # include "MOM_FLUXFORM.h"
59     #endif
60 adcroft 1.1
61 adcroft 1.3 C !INPUT PARAMETERS: ===================================================
62     C bi,bj :: tile indices
63     C iMin,iMax,jMin,jMAx :: loop ranges
64     C k :: vertical level
65     C kUp :: =1 or 2 for consecutive k
66     C kDown :: =2 or 1 for consecutive k
67     C KappaRU :: vertical viscosity
68     C KappaRV :: vertical viscosity
69     C fVerU :: vertical flux of U, 2 1/2 dim for pipe-lining
70     C fVerV :: vertical flux of V, 2 1/2 dim for pipe-lining
71 jmc 1.23 C guDiss :: dissipation tendency (all explicit terms), u component
72     C gvDiss :: dissipation tendency (all explicit terms), v component
73 jmc 1.8 C myTime :: current time
74 adcroft 1.3 C myIter :: current time-step number
75     C myThid :: thread number
76     INTEGER bi,bj,iMin,iMax,jMin,jMax
77     INTEGER k,kUp,kDown
78 adcroft 1.1 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
79     _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
80     _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
81     _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
82 jmc 1.23 _RL guDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83     _RL gvDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84 jmc 1.8 _RL myTime
85 adcroft 1.2 INTEGER myIter
86 adcroft 1.1 INTEGER myThid
87    
88 adcroft 1.3 C !OUTPUT PARAMETERS: ==================================================
89     C None - updates gU() and gV() in common blocks
90    
91     C !LOCAL VARIABLES: ====================================================
92     C i,j :: loop indices
93     C vF :: viscous flux
94     C v4F :: bi-harmonic viscous flux
95     C cF :: Coriolis acceleration
96     C mT :: Metric terms
97     C fZon :: zonal fluxes
98     C fMer :: meridional fluxes
99 jmc 1.23 C fVrUp,fVrDw :: vertical viscous fluxes at interface k-1 & k
100 adcroft 1.3 INTEGER i,j
101 heimbach 1.35 #ifdef ALLOW_AUTODIFF_TAMC
102     INTEGER imomkey
103     #endif
104 adcroft 1.3 _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106     _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107     _RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110 jmc 1.23 _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111     _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112 jmc 1.33 C afFacMom :: Tracer parameters for turning terms on and off.
113     C vfFacMom
114 adcroft 1.1 C pfFacMom afFacMom - Advective terms
115     C cfFacMom vfFacMom - Eddy viscosity terms
116 jmc 1.33 C mtFacMom pfFacMom - Pressure terms
117 adcroft 1.1 C cfFacMom - Coriolis terms
118     C foFacMom - Forcing
119 jmc 1.33 C mtFacMom - Metric term
120 jmc 1.23 C uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off
121 adcroft 1.1 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122     _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123     _RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124     _RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127     _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128     _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 jmc 1.8 _RL rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130     _RL rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131 adcroft 1.18 _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 baylor 1.25 _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133     _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134     _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135     _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136     _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137     _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138 adcroft 1.18 _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139     _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140 adcroft 1.1 _RL uDudxFac
141     _RL AhDudxFac
142     _RL vDudyFac
143     _RL AhDudyFac
144     _RL rVelDudrFac
145     _RL ArDudrFac
146     _RL fuFac
147     _RL mtFacU
148 jmc 1.33 _RL mtNHFacU
149 adcroft 1.1 _RL uDvdxFac
150     _RL AhDvdxFac
151     _RL vDvdyFac
152     _RL AhDvdyFac
153     _RL rVelDvdrFac
154     _RL ArDvdrFac
155     _RL fvFac
156     _RL mtFacV
157 jmc 1.33 _RL mtNHFacV
158 jmc 1.29 _RL sideMaskFac
159 baylor 1.25 LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity
160 adcroft 1.3 CEOP
161 adcroft 1.1
162 heimbach 1.35 #ifdef ALLOW_AUTODIFF_TAMC
163     act0 = k - 1
164     max0 = Nr
165     act1 = bi - myBxLo(myThid)
166     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
167     act2 = bj - myByLo(myThid)
168     max2 = myByHi(myThid) - myByLo(myThid) + 1
169     act3 = myThid - 1
170     max3 = nTx*nTy
171     act4 = ikey_dynamics - 1
172     imomkey = (act0 + 1)
173     & + act1*max0
174     & + act2*max0*max1
175     & + act3*max0*max1*max2
176     & + act4*max0*max1*max2*max3
177     #endif /* ALLOW_AUTODIFF_TAMC */
178    
179 adcroft 1.1 C Initialise intermediate terms
180 jmc 1.23 DO j=1-OLy,sNy+OLy
181     DO i=1-OLx,sNx+OLx
182 adcroft 1.1 vF(i,j) = 0.
183     v4F(i,j) = 0.
184     cF(i,j) = 0.
185     mT(i,j) = 0.
186     fZon(i,j) = 0.
187     fMer(i,j) = 0.
188 jmc 1.23 fVrUp(i,j)= 0.
189     fVrDw(i,j)= 0.
190     rTransU(i,j)= 0.
191     rTransV(i,j)= 0.
192 adcroft 1.18 strain(i,j) = 0.
193 jmc 1.23 tension(i,j)= 0.
194     guDiss(i,j) = 0.
195     gvDiss(i,j) = 0.
196 heimbach 1.30 #ifdef ALLOW_AUTODIFF_TAMC
197     vort3(i,j) = 0. _d 0
198     strain(i,j) = 0. _d 0
199     tension(i,j) = 0. _d 0
200     #endif
201 adcroft 1.1 ENDDO
202     ENDDO
203    
204     C-- Term by term tracer parmeters
205     C o U momentum equation
206     uDudxFac = afFacMom*1.
207     AhDudxFac = vfFacMom*1.
208     vDudyFac = afFacMom*1.
209     AhDudyFac = vfFacMom*1.
210     rVelDudrFac = afFacMom*1.
211     ArDudrFac = vfFacMom*1.
212 jmc 1.33 mtFacU = mtFacMom*1.
213     mtNHFacU = 1.
214 adcroft 1.1 fuFac = cfFacMom*1.
215     C o V momentum equation
216     uDvdxFac = afFacMom*1.
217     AhDvdxFac = vfFacMom*1.
218     vDvdyFac = afFacMom*1.
219     AhDvdyFac = vfFacMom*1.
220     rVelDvdrFac = afFacMom*1.
221     ArDvdrFac = vfFacMom*1.
222 jmc 1.33 mtFacV = mtFacMom*1.
223     mtNHFacV = 1.
224 adcroft 1.1 fvFac = cfFacMom*1.
225 jmc 1.23
226     IF (implicitViscosity) THEN
227     ArDudrFac = 0.
228     ArDvdrFac = 0.
229     ENDIF
230 adcroft 1.1
231 jmc 1.29 C note: using standard stencil (no mask) results in under-estimating
232     C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
233     IF ( no_slip_sides ) THEN
234     sideMaskFac = sideDragFactor
235     ELSE
236     sideMaskFac = 0. _d 0
237     ENDIF
238    
239 adcroft 1.1 IF ( no_slip_bottom
240     & .OR. bottomDragQuadratic.NE.0.
241     & .OR. bottomDragLinear.NE.0.) THEN
242     bottomDragTerms=.TRUE.
243     ELSE
244     bottomDragTerms=.FALSE.
245     ENDIF
246    
247     C-- Calculate open water fraction at vorticity points
248     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
249    
250     C---- Calculate common quantities used in both U and V equations
251     C Calculate tracer cell face open areas
252     DO j=1-OLy,sNy+OLy
253     DO i=1-OLx,sNx+OLx
254     xA(i,j) = _dyG(i,j,bi,bj)
255     & *drF(k)*_hFacW(i,j,k,bi,bj)
256     yA(i,j) = _dxG(i,j,bi,bj)
257     & *drF(k)*_hFacS(i,j,k,bi,bj)
258     ENDDO
259     ENDDO
260    
261     C Make local copies of horizontal flow field
262     DO j=1-OLy,sNy+OLy
263     DO i=1-OLx,sNx+OLx
264     uFld(i,j) = uVel(i,j,k,bi,bj)
265     vFld(i,j) = vVel(i,j,k,bi,bj)
266     ENDDO
267     ENDDO
268    
269     C Calculate velocity field "volume transports" through tracer cell faces.
270     DO j=1-OLy,sNy+OLy
271     DO i=1-OLx,sNx+OLx
272     uTrans(i,j) = uFld(i,j)*xA(i,j)
273     vTrans(i,j) = vFld(i,j)*yA(i,j)
274     ENDDO
275     ENDDO
276    
277 baylor 1.25 CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)
278 jmc 1.29 IF ( momViscosity) THEN
279     CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
280     CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
281     CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
282     CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
283     DO j=1-Oly,sNy+Oly
284     DO i=1-Olx,sNx+Olx
285     IF ( hFacZ(i,j).EQ.0. ) THEN
286     vort3(i,j) = sideMaskFac*vort3(i,j)
287     strain(i,j) = sideMaskFac*strain(i,j)
288     ENDIF
289     ENDDO
290     ENDDO
291     #ifdef ALLOW_DIAGNOSTICS
292     IF ( useDiagnostics ) THEN
293     CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
294     CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
295     CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
296     CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
297     ENDIF
298     #endif
299     ENDIF
300 adcroft 1.18
301 jmc 1.8 C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
302     IF (momAdvection.AND.k.EQ.1) THEN
303    
304     C- Calculate vertical transports above U & V points (West & South face):
305 heimbach 1.35
306     #ifdef ALLOW_AUTODIFF_TAMC
307     CADJ STORE dwtransc(:,:,bi,bj) =
308     CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
309     CADJ STORE dwtransu(:,:,bi,bj) =
310     CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
311     CADJ STORE dwtransv(:,:,bi,bj) =
312     CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
313     #endif /* ALLOW_AUTODIFF_TAMC */
314 jmc 1.23 CALL MOM_CALC_RTRANS( k, bi, bj,
315     O rTransU, rTransV,
316     I myTime, myIter, myThid)
317 jmc 1.8
318     C- Free surface correction term (flux at k=1)
319 jmc 1.23 CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
320     O fVerU(1-OLx,1-OLy,kUp), myThid )
321 jmc 1.8
322 jmc 1.23 CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
323     O fVerV(1-OLx,1-OLy,kUp), myThid )
324 jmc 1.8
325     C--- endif momAdvection & k=1
326     ENDIF
327    
328    
329     C--- Calculate vertical transports (at k+1) below U & V points :
330     IF (momAdvection) THEN
331 jmc 1.23 CALL MOM_CALC_RTRANS( k+1, bi, bj,
332     O rTransU, rTransV,
333     I myTime, myIter, myThid)
334 jmc 1.8 ENDIF
335    
336 baylor 1.25 IF (momViscosity) THEN
337     CALL MOM_CALC_VISC(
338     I bi,bj,k,
339     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
340     O harmonic,biharmonic,useVariableViscosity,
341 jmc 1.26 I hDiv,vort3,tension,strain,KE,hFacZ,
342 baylor 1.25 I myThid)
343     ENDIF
344 jmc 1.8
345 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
346    
347 adcroft 1.1 C---- Zonal momentum equation starts here
348    
349 jmc 1.23 IF (momAdvection) THEN
350     C--- Calculate mean fluxes (advection) between cells for zonal flow.
351 adcroft 1.1
352     C-- Zonal flux (fZon is at east face of "u" cell)
353 jmc 1.23 C Mean flow component of zonal flux -> fZon
354     CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,fZon,myThid)
355 adcroft 1.1
356     C-- Meridional flux (fMer is at south face of "u" cell)
357 jmc 1.23 C Mean flow component of meridional flux -> fMer
358     CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,fMer,myThid)
359 adcroft 1.1
360     C-- Vertical flux (fVer is at upper face of "u" cell)
361 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVer
362     CALL MOM_U_ADV_WU(
363     I bi,bj,k+1,uVel,wVel,rTransU,
364     O fVerU(1-OLx,1-OLy,kDown), myThid )
365 adcroft 1.1
366     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
367 jmc 1.23 DO j=jMin,jMax
368     DO i=iMin,iMax
369     gU(i,j,k,bi,bj) =
370 adcroft 1.1 #ifdef OLD_UV_GEOM
371 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
372     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
373 adcroft 1.1 #else
374 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
375     & *recip_rAw(i,j,bi,bj)
376 adcroft 1.1 #endif
377 jmc 1.23 & *( ( fZon(i,j ) - fZon(i-1,j) )*uDudxFac
378     & +( fMer(i,j+1) - fMer(i, j) )*vDudyFac
379     & +(fVerU(i,j,kDown) - fVerU(i,j,kUp))*rkSign*rVelDudrFac
380     & )
381     ENDDO
382     ENDDO
383 adcroft 1.1
384 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
385     IF ( useDiagnostics ) THEN
386     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Um ',k,1,2,bi,bj,myThid)
387     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Um ',k,1,2,bi,bj,myThid)
388     CALL DIAGNOSTICS_FILL(fVerU(1-Olx,1-Oly,kUp),
389     & 'ADVrE_Um',k,1,2,bi,bj,myThid)
390     ENDIF
391     #endif
392    
393 jmc 1.8 #ifdef NONLIN_FRSURF
394     C-- account for 3.D divergence of the flow in rStar coordinate:
395 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
396 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
397     DO j=jMin,jMax
398     DO i=iMin,iMax
399     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
400 jmc 1.8 & - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
401     & *uVel(i,j,k,bi,bj)
402 jmc 1.23 ENDDO
403     ENDDO
404     ENDIF
405     IF ( select_rStar.LT.0 ) THEN
406     DO j=jMin,jMax
407     DO i=iMin,iMax
408     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
409     & - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
410     ENDDO
411     ENDDO
412     ENDIF
413 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
414 jmc 1.23 #endif /* NONLIN_FRSURF */
415    
416     ELSE
417     C- if momAdvection / else
418     DO j=1-OLy,sNy+OLy
419     DO i=1-OLx,sNx+OLx
420     gU(i,j,k,bi,bj) = 0. _d 0
421     ENDDO
422 jmc 1.8 ENDDO
423 jmc 1.23
424     C- endif momAdvection.
425 jmc 1.8 ENDIF
426 jmc 1.23
427     IF (momViscosity) THEN
428     C--- Calculate eddy fluxes (dissipation) between cells for zonal flow.
429    
430     C Bi-harmonic term del^2 U -> v4F
431 baylor 1.25 IF (biharmonic)
432 jmc 1.23 & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
433    
434     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
435 baylor 1.25 CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,
436 baylor 1.27 I viscAh_D,viscA4_D,myThid)
437 jmc 1.23
438     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
439 baylor 1.25 CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,
440 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
441 jmc 1.23
442     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
443     IF (.NOT.implicitViscosity) THEN
444     CALL MOM_U_RVISCFLUX(bi,bj, k, uVel,KappaRU,fVrUp,myThid)
445     CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,fVrDw,myThid)
446     ENDIF
447    
448     C-- Tendency is minus divergence of the fluxes
449     DO j=jMin,jMax
450     DO i=iMin,iMax
451     guDiss(i,j) =
452     #ifdef OLD_UV_GEOM
453     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
454     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
455     #else
456     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
457     & *recip_rAw(i,j,bi,bj)
458     #endif
459     & *( ( fZon(i,j ) - fZon(i-1,j) )*AhDudxFac
460     & +( fMer(i,j+1) - fMer(i, j) )*AhDudyFac
461     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDudrFac
462     & )
463     ENDDO
464 jmc 1.8 ENDDO
465    
466 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
467     IF ( useDiagnostics ) THEN
468     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
469     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
470     IF (.NOT.implicitViscosity)
471     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
472     ENDIF
473     #endif
474    
475 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
476 jmc 1.23 IF (no_slip_sides) THEN
477 adcroft 1.1 C- No-slip BCs impose a drag at walls...
478 baylor 1.27 CALL MOM_U_SIDEDRAG(
479     I bi,bj,k,
480     I uFld, v4f, hFacZ,
481     I viscAh_Z,viscA4_Z,
482     I harmonic,biharmonic,useVariableViscosity,
483     O vF,
484     I myThid)
485 jmc 1.23 DO j=jMin,jMax
486     DO i=iMin,iMax
487     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
488     ENDDO
489     ENDDO
490     ENDIF
491 adcroft 1.1 C- No-slip BCs impose a drag at bottom
492 jmc 1.23 IF (bottomDragTerms) THEN
493     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
494     DO j=jMin,jMax
495     DO i=iMin,iMax
496     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
497     ENDDO
498     ENDDO
499     ENDIF
500    
501 mlosch 1.32 #ifdef ALLOW_SHELFICE
502     IF (useShelfIce) THEN
503     CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
504     DO j=jMin,jMax
505     DO i=iMin,iMax
506     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
507     ENDDO
508     ENDDO
509     ENDIF
510     #endif /* ALLOW_SHELFICE */
511    
512 jmc 1.23 C- endif momViscosity
513 adcroft 1.1 ENDIF
514    
515 jmc 1.12 C-- Forcing term (moved to timestep.F)
516     c IF (momForcing)
517     c & CALL EXTERNAL_FORCING_U(
518     c I iMin,iMax,jMin,jMax,bi,bj,k,
519     c I myTime,myThid)
520 adcroft 1.1
521     C-- Metric terms for curvilinear grid systems
522 adcroft 1.5 IF (useNHMTerms) THEN
523 jmc 1.33 C o Non-Hydrostatic (spherical) metric terms
524 adcroft 1.1 CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
525     DO j=jMin,jMax
526     DO i=iMin,iMax
527 jmc 1.33 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j)
528 adcroft 1.1 ENDDO
529     ENDDO
530 adcroft 1.5 ENDIF
531 jmc 1.33 IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
532     C o Spherical polar grid metric terms
533 adcroft 1.1 CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
534     DO j=jMin,jMax
535     DO i=iMin,iMax
536 jmc 1.33 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
537 adcroft 1.1 ENDDO
538     ENDDO
539 afe 1.20 ENDIF
540 jmc 1.33 IF ( usingCylindricalGrid .AND. metricTerms ) THEN
541     C o Cylindrical grid metric terms
542     CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
543     DO j=jMin,jMax
544     DO i=iMin,iMax
545     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
546     ENDDO
547 afe 1.19 ENDDO
548 adcroft 1.1 ENDIF
549    
550 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
551 adcroft 1.1
552     C---- Meridional momentum equation starts here
553    
554 jmc 1.23 IF (momAdvection) THEN
555     C--- Calculate mean fluxes (advection) between cells for meridional flow.
556     C Mean flow component of zonal flux -> fZon
557     CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,fZon,myThid)
558 adcroft 1.1
559     C-- Meridional flux (fMer is at north face of "v" cell)
560 jmc 1.23 C Mean flow component of meridional flux -> fMer
561     CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,fMer,myThid)
562 adcroft 1.1
563     C-- Vertical flux (fVer is at upper face of "v" cell)
564 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVerV
565     CALL MOM_V_ADV_WV(
566     I bi,bj,k+1,vVel,wVel,rTransV,
567     O fVerV(1-OLx,1-OLy,kDown), myThid )
568 adcroft 1.1
569     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
570 jmc 1.23 DO j=jMin,jMax
571     DO i=iMin,iMax
572     gV(i,j,k,bi,bj) =
573 adcroft 1.1 #ifdef OLD_UV_GEOM
574 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
575     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
576 adcroft 1.1 #else
577 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
578     & *recip_rAs(i,j,bi,bj)
579 adcroft 1.1 #endif
580 jmc 1.23 & *( ( fZon(i+1,j) - fZon(i,j ) )*uDvdxFac
581     & +( fMer(i, j) - fMer(i,j-1) )*vDvdyFac
582     & +(fVerV(i,j,kDown) - fVerV(i,j,kUp))*rkSign*rVelDvdrFac
583     & )
584 jmc 1.24 ENDDO
585     ENDDO
586    
587     #ifdef ALLOW_DIAGNOSTICS
588     IF ( useDiagnostics ) THEN
589     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Vm ',k,1,2,bi,bj,myThid)
590     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Vm ',k,1,2,bi,bj,myThid)
591     CALL DIAGNOSTICS_FILL(fVerV(1-Olx,1-Oly,kUp),
592     & 'ADVrE_Vm',k,1,2,bi,bj,myThid)
593     ENDIF
594     #endif
595 adcroft 1.1
596 jmc 1.8 #ifdef NONLIN_FRSURF
597     C-- account for 3.D divergence of the flow in rStar coordinate:
598 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
599 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
600     DO j=jMin,jMax
601     DO i=iMin,iMax
602     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
603 jmc 1.8 & - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
604     & *vVel(i,j,k,bi,bj)
605 jmc 1.23 ENDDO
606     ENDDO
607     ENDIF
608     IF ( select_rStar.LT.0 ) THEN
609     DO j=jMin,jMax
610     DO i=iMin,iMax
611     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
612     & - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
613     ENDDO
614     ENDDO
615     ENDIF
616 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
617 jmc 1.23 #endif /* NONLIN_FRSURF */
618    
619     ELSE
620     C- if momAdvection / else
621     DO j=1-OLy,sNy+OLy
622     DO i=1-OLx,sNx+OLx
623     gV(i,j,k,bi,bj) = 0. _d 0
624     ENDDO
625 jmc 1.8 ENDDO
626 jmc 1.23
627     C- endif momAdvection.
628 jmc 1.8 ENDIF
629 jmc 1.23
630     IF (momViscosity) THEN
631     C--- Calculate eddy fluxes (dissipation) between cells for meridional flow.
632     C Bi-harmonic term del^2 V -> v4F
633 baylor 1.25 IF (biharmonic)
634 jmc 1.23 & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
635    
636     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
637 baylor 1.25 CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,
638 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
639 jmc 1.23
640     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
641 baylor 1.25 CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,
642 baylor 1.27 I viscAh_D,viscA4_D,myThid)
643 jmc 1.23
644     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
645     IF (.NOT.implicitViscosity) THEN
646     CALL MOM_V_RVISCFLUX(bi,bj, k, vVel,KappaRV,fVrUp,myThid)
647     CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,fVrDw,myThid)
648     ENDIF
649    
650     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
651     DO j=jMin,jMax
652     DO i=iMin,iMax
653     gvDiss(i,j) =
654     #ifdef OLD_UV_GEOM
655     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
656     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
657     #else
658     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
659     & *recip_rAs(i,j,bi,bj)
660     #endif
661     & *( ( fZon(i+1,j) - fZon(i,j ) )*AhDvdxFac
662     & +( fMer(i, j) - fMer(i,j-1) )*AhDvdyFac
663     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDvdrFac
664     & )
665     ENDDO
666 jmc 1.8 ENDDO
667    
668 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
669     IF ( useDiagnostics ) THEN
670     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
671     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
672     IF (.NOT.implicitViscosity)
673     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
674     ENDIF
675     #endif
676    
677 adcroft 1.1 C-- No-slip and drag BCs appear as body forces in cell abutting topography
678 mlosch 1.32 IF (no_slip_sides) THEN
679 adcroft 1.1 C- No-slip BCs impose a drag at walls...
680 baylor 1.27 CALL MOM_V_SIDEDRAG(
681     I bi,bj,k,
682     I vFld, v4f, hFacZ,
683     I viscAh_Z,viscA4_Z,
684     I harmonic,biharmonic,useVariableViscosity,
685     O vF,
686     I myThid)
687 jmc 1.23 DO j=jMin,jMax
688     DO i=iMin,iMax
689     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
690     ENDDO
691     ENDDO
692     ENDIF
693 adcroft 1.1 C- No-slip BCs impose a drag at bottom
694 jmc 1.23 IF (bottomDragTerms) THEN
695     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
696     DO j=jMin,jMax
697     DO i=iMin,iMax
698     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
699     ENDDO
700     ENDDO
701     ENDIF
702    
703 mlosch 1.32 #ifdef ALLOW_SHELFICE
704     IF (useShelfIce) THEN
705     CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
706     DO j=jMin,jMax
707     DO i=iMin,iMax
708     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
709     ENDDO
710     ENDDO
711     ENDIF
712     #endif /* ALLOW_SHELFICE */
713    
714 jmc 1.23 C- endif momViscosity
715 adcroft 1.1 ENDIF
716    
717 jmc 1.12 C-- Forcing term (moved to timestep.F)
718     c IF (momForcing)
719     c & CALL EXTERNAL_FORCING_V(
720     c I iMin,iMax,jMin,jMax,bi,bj,k,
721     c I myTime,myThid)
722 adcroft 1.1
723     C-- Metric terms for curvilinear grid systems
724 adcroft 1.5 IF (useNHMTerms) THEN
725 jmc 1.33 C o Non-Hydrostatic (spherical) metric terms
726 adcroft 1.1 CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
727     DO j=jMin,jMax
728     DO i=iMin,iMax
729 jmc 1.33 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j)
730 adcroft 1.1 ENDDO
731     ENDDO
732 adcroft 1.5 ENDIF
733 jmc 1.33 IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
734     C o Spherical polar grid metric terms
735 adcroft 1.1 CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
736     DO j=jMin,jMax
737     DO i=iMin,iMax
738 jmc 1.33 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
739 adcroft 1.1 ENDDO
740     ENDDO
741     ENDIF
742 jmc 1.33 IF ( usingCylindricalGrid .AND. metricTerms ) THEN
743     C o Cylindrical grid metric terms
744     CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
745     DO j=jMin,jMax
746     DO i=iMin,iMax
747     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
748     ENDDO
749     ENDDO
750 afe 1.19 ENDIF
751 adcroft 1.1
752 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
753 adcroft 1.1
754     C-- Coriolis term
755     C Note. As coded here, coriolis will not work with "thin walls"
756 jmc 1.12 c IF (useCDscheme) THEN
757     c CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
758     c ELSE
759     IF (.NOT.useCDscheme) THEN
760     CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
761     DO j=jMin,jMax
762     DO i=iMin,iMax
763     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
764     ENDDO
765     ENDDO
766 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
767     IF ( useDiagnostics )
768     & CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)
769     #endif
770 jmc 1.12 CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
771     DO j=jMin,jMax
772     DO i=iMin,iMax
773     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
774     ENDDO
775     ENDDO
776 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
777     IF ( useDiagnostics )
778     & CALL DIAGNOSTICS_FILL(cf,'Vm_Cori ',k,1,2,bi,bj,myThid)
779     #endif
780 jmc 1.12 ENDIF
781    
782 jmc 1.34 C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
783     IF ( nonHydrostatic.OR.quasiHydrostatic ) THEN
784     CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
785     DO j=jMin,jMax
786     DO i=iMin,iMax
787     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
788     ENDDO
789     ENDDO
790     IF ( usingCurvilinearGrid ) THEN
791     C- presently, non zero angleSinC array only supported with Curvilinear-Grid
792     CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
793     DO j=jMin,jMax
794     DO i=iMin,iMax
795     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
796     ENDDO
797 adcroft 1.6 ENDDO
798 jmc 1.34 ENDIF
799 adcroft 1.6 ENDIF
800 adcroft 1.1
801 jmc 1.23 C-- Set du/dt & dv/dt on boundaries to zero
802     DO j=jMin,jMax
803     DO i=iMin,iMax
804     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
805     guDiss(i,j) = guDiss(i,j) *_maskW(i,j,k,bi,bj)
806     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
807     gvDiss(i,j) = gvDiss(i,j) *_maskS(i,j,k,bi,bj)
808     ENDDO
809     ENDDO
810    
811 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
812     IF ( useDiagnostics ) THEN
813 baylor 1.28 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
814 jmc 1.24 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
815     & 'Um_Advec',k,1,2,bi,bj,myThid)
816     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
817     & 'Vm_Advec',k,1,2,bi,bj,myThid)
818     IF (momViscosity) THEN
819     CALL DIAGNOSTICS_FILL(guDiss,'Um_Diss ',k,1,2,bi,bj,myThid)
820     CALL DIAGNOSTICS_FILL(gvDiss,'Vm_Diss ',k,1,2,bi,bj,myThid)
821     ENDIF
822     ENDIF
823     #endif /* ALLOW_DIAGNOSTICS */
824    
825 adcroft 1.1 RETURN
826     END

  ViewVC Help
Powered by ViewVC 1.1.22