/[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.39 - (hide annotations) (download)
Tue Dec 5 05:30:38 2006 UTC (17 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58t_post, checkpoint58w_post, mitgcm_mapl_00, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59i, checkpoint59h, checkpoint58v_post, checkpoint58x_post, checkpoint58u_post, checkpoint58s_post
Changes since 1.38: +32 -27 lines
start to implement deep-atmosphere and/or anelastic formulation

1 jmc 1.39 C $Header: /u/gcmpack/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.38 2006/11/23 00:45:18 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 jmc 1.37 SUBROUTINE MOM_FLUXFORM(
35 adcroft 1.1 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 jmc 1.39 C pressure gradient and implicit 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 jmc 1.37 #ifdef ALLOW_AUTODIFF_TAMC
56 heimbach 1.35 # 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 jmc 1.37 #ifdef ALLOW_AUTODIFF_TAMC
102 heimbach 1.35 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 jmc 1.37 C vfFacMom
114     C pfFacMom afFacMom - Advective terms
115 adcroft 1.1 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 jmc 1.37 imomkey = (act0 + 1)
173 heimbach 1.35 & + 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 jmc 1.38 c KE(i,j) = 0.
193     c hDiv(i,j) = 0.
194     vort3(i,j) = 0.
195 adcroft 1.18 strain(i,j) = 0.
196 jmc 1.23 tension(i,j)= 0.
197     guDiss(i,j) = 0.
198     gvDiss(i,j) = 0.
199 adcroft 1.1 ENDDO
200     ENDDO
201    
202     C-- Term by term tracer parmeters
203     C o U momentum equation
204     uDudxFac = afFacMom*1.
205     AhDudxFac = vfFacMom*1.
206     vDudyFac = afFacMom*1.
207     AhDudyFac = vfFacMom*1.
208     rVelDudrFac = afFacMom*1.
209     ArDudrFac = vfFacMom*1.
210 jmc 1.33 mtFacU = mtFacMom*1.
211     mtNHFacU = 1.
212 adcroft 1.1 fuFac = cfFacMom*1.
213     C o V momentum equation
214     uDvdxFac = afFacMom*1.
215     AhDvdxFac = vfFacMom*1.
216     vDvdyFac = afFacMom*1.
217     AhDvdyFac = vfFacMom*1.
218     rVelDvdrFac = afFacMom*1.
219     ArDvdrFac = vfFacMom*1.
220 jmc 1.33 mtFacV = mtFacMom*1.
221     mtNHFacV = 1.
222 adcroft 1.1 fvFac = cfFacMom*1.
223 jmc 1.23
224     IF (implicitViscosity) THEN
225     ArDudrFac = 0.
226     ArDvdrFac = 0.
227     ENDIF
228 adcroft 1.1
229 jmc 1.29 C note: using standard stencil (no mask) results in under-estimating
230     C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
231     IF ( no_slip_sides ) THEN
232     sideMaskFac = sideDragFactor
233     ELSE
234     sideMaskFac = 0. _d 0
235     ENDIF
236    
237 adcroft 1.1 IF ( no_slip_bottom
238     & .OR. bottomDragQuadratic.NE.0.
239     & .OR. bottomDragLinear.NE.0.) THEN
240     bottomDragTerms=.TRUE.
241     ELSE
242     bottomDragTerms=.FALSE.
243     ENDIF
244    
245     C-- Calculate open water fraction at vorticity points
246     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
247    
248     C---- Calculate common quantities used in both U and V equations
249     C Calculate tracer cell face open areas
250     DO j=1-OLy,sNy+OLy
251     DO i=1-OLx,sNx+OLx
252 jmc 1.39 xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
253     & *drF(k)*_hFacW(i,j,k,bi,bj)
254     yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
255     & *drF(k)*_hFacS(i,j,k,bi,bj)
256 adcroft 1.1 ENDDO
257     ENDDO
258    
259     C Make local copies of horizontal flow field
260     DO j=1-OLy,sNy+OLy
261     DO i=1-OLx,sNx+OLx
262     uFld(i,j) = uVel(i,j,k,bi,bj)
263     vFld(i,j) = vVel(i,j,k,bi,bj)
264     ENDDO
265     ENDDO
266    
267     C Calculate velocity field "volume transports" through tracer cell faces.
268 jmc 1.39 C anelastic: transports are scaled by rhoFacC (~ mass transport)
269 adcroft 1.1 DO j=1-OLy,sNy+OLy
270     DO i=1-OLx,sNx+OLx
271 jmc 1.39 uTrans(i,j) = uFld(i,j)*xA(i,j)*rhoFacC(k)
272     vTrans(i,j) = vFld(i,j)*yA(i,j)*rhoFacC(k)
273 adcroft 1.1 ENDDO
274     ENDDO
275    
276 baylor 1.25 CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)
277 jmc 1.29 IF ( momViscosity) THEN
278     CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
279     CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
280     CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
281     CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
282     DO j=1-Oly,sNy+Oly
283     DO i=1-Olx,sNx+Olx
284     IF ( hFacZ(i,j).EQ.0. ) THEN
285     vort3(i,j) = sideMaskFac*vort3(i,j)
286     strain(i,j) = sideMaskFac*strain(i,j)
287     ENDIF
288     ENDDO
289     ENDDO
290     #ifdef ALLOW_DIAGNOSTICS
291     IF ( useDiagnostics ) THEN
292     CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid)
293     CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid)
294     CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
295     CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid)
296     ENDIF
297     #endif
298     ENDIF
299 adcroft 1.18
300 jmc 1.8 C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
301     IF (momAdvection.AND.k.EQ.1) THEN
302    
303     C- Calculate vertical transports above U & V points (West & South face):
304 heimbach 1.35
305     #ifdef ALLOW_AUTODIFF_TAMC
306 heimbach 1.36 # ifdef NONLIN_FRSURF
307     # ifndef DISABLE_RSTAR_CODE
308 jmc 1.39 CADJ STORE dwtransc(:,:,bi,bj) =
309 heimbach 1.35 CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
310 jmc 1.39 CADJ STORE dwtransu(:,:,bi,bj) =
311 heimbach 1.35 CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
312 jmc 1.39 CADJ STORE dwtransv(:,:,bi,bj) =
313 heimbach 1.35 CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte
314 heimbach 1.36 # endif
315     # endif /* NONLIN_FRSURF */
316 heimbach 1.35 #endif /* ALLOW_AUTODIFF_TAMC */
317 jmc 1.23 CALL MOM_CALC_RTRANS( k, bi, bj,
318     O rTransU, rTransV,
319     I myTime, myIter, myThid)
320 jmc 1.8
321     C- Free surface correction term (flux at k=1)
322 jmc 1.23 CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU,
323     O fVerU(1-OLx,1-OLy,kUp), myThid )
324 jmc 1.8
325 jmc 1.23 CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV,
326     O fVerV(1-OLx,1-OLy,kUp), myThid )
327 jmc 1.8
328     C--- endif momAdvection & k=1
329     ENDIF
330    
331    
332     C--- Calculate vertical transports (at k+1) below U & V points :
333     IF (momAdvection) THEN
334 jmc 1.23 CALL MOM_CALC_RTRANS( k+1, bi, bj,
335     O rTransU, rTransV,
336     I myTime, myIter, myThid)
337 jmc 1.8 ENDIF
338    
339 baylor 1.25 IF (momViscosity) THEN
340     CALL MOM_CALC_VISC(
341     I bi,bj,k,
342     O viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
343     O harmonic,biharmonic,useVariableViscosity,
344 jmc 1.26 I hDiv,vort3,tension,strain,KE,hFacZ,
345 baylor 1.25 I myThid)
346     ENDIF
347 jmc 1.8
348 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
349    
350 adcroft 1.1 C---- Zonal momentum equation starts here
351    
352 jmc 1.23 IF (momAdvection) THEN
353     C--- Calculate mean fluxes (advection) between cells for zonal flow.
354 adcroft 1.1
355     C-- Zonal flux (fZon is at east face of "u" cell)
356 jmc 1.23 C Mean flow component of zonal flux -> fZon
357     CALL MOM_U_ADV_UU(bi,bj,k,uTrans,uFld,fZon,myThid)
358 adcroft 1.1
359     C-- Meridional flux (fMer is at south face of "u" cell)
360 jmc 1.23 C Mean flow component of meridional flux -> fMer
361     CALL MOM_U_ADV_VU(bi,bj,k,vTrans,uFld,fMer,myThid)
362 adcroft 1.1
363     C-- Vertical flux (fVer is at upper face of "u" cell)
364 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVer
365     CALL MOM_U_ADV_WU(
366     I bi,bj,k+1,uVel,wVel,rTransU,
367     O fVerU(1-OLx,1-OLy,kDown), myThid )
368 adcroft 1.1
369     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
370 jmc 1.23 DO j=jMin,jMax
371     DO i=iMin,iMax
372     gU(i,j,k,bi,bj) =
373 adcroft 1.1 #ifdef OLD_UV_GEOM
374 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
375     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
376 adcroft 1.1 #else
377 jmc 1.23 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
378 jmc 1.39 & *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
379 adcroft 1.1 #endif
380 jmc 1.39 & *( ( fZon(i,j ) - fZon(i-1,j) )*uDudxFac
381     & +( fMer(i,j+1) - fMer(i, j) )*vDudyFac
382     & +(fVerU(i,j,kDown) - fVerU(i,j,kUp))*rkSign*rVelDudrFac
383 jmc 1.23 & )
384     ENDDO
385     ENDDO
386 adcroft 1.1
387 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
388     IF ( useDiagnostics ) THEN
389     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Um ',k,1,2,bi,bj,myThid)
390     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Um ',k,1,2,bi,bj,myThid)
391     CALL DIAGNOSTICS_FILL(fVerU(1-Olx,1-Oly,kUp),
392     & 'ADVrE_Um',k,1,2,bi,bj,myThid)
393     ENDIF
394     #endif
395    
396 jmc 1.8 #ifdef NONLIN_FRSURF
397     C-- account for 3.D divergence of the flow in rStar coordinate:
398 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
399 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
400     DO j=jMin,jMax
401     DO i=iMin,iMax
402     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
403 jmc 1.8 & - (rStarExpW(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
404     & *uVel(i,j,k,bi,bj)
405 jmc 1.23 ENDDO
406     ENDDO
407     ENDIF
408     IF ( select_rStar.LT.0 ) THEN
409     DO j=jMin,jMax
410     DO i=iMin,iMax
411     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
412     & - rStarDhWDt(i,j,bi,bj)*uVel(i,j,k,bi,bj)
413     ENDDO
414     ENDDO
415     ENDIF
416 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
417 jmc 1.23 #endif /* NONLIN_FRSURF */
418    
419     ELSE
420     C- if momAdvection / else
421     DO j=1-OLy,sNy+OLy
422     DO i=1-OLx,sNx+OLx
423     gU(i,j,k,bi,bj) = 0. _d 0
424     ENDDO
425 jmc 1.8 ENDDO
426 jmc 1.23
427     C- endif momAdvection.
428 jmc 1.8 ENDIF
429 jmc 1.23
430     IF (momViscosity) THEN
431     C--- Calculate eddy fluxes (dissipation) between cells for zonal flow.
432    
433     C Bi-harmonic term del^2 U -> v4F
434 jmc 1.37 IF (biharmonic)
435 jmc 1.23 & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid)
436    
437     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
438 baylor 1.25 CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,
439 baylor 1.27 I viscAh_D,viscA4_D,myThid)
440 jmc 1.23
441     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
442 baylor 1.25 CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,
443 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
444 jmc 1.23
445     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
446     IF (.NOT.implicitViscosity) THEN
447     CALL MOM_U_RVISCFLUX(bi,bj, k, uVel,KappaRU,fVrUp,myThid)
448     CALL MOM_U_RVISCFLUX(bi,bj,k+1,uVel,KappaRU,fVrDw,myThid)
449     ENDIF
450    
451     C-- Tendency is minus divergence of the fluxes
452 jmc 1.39 C anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
453 jmc 1.23 DO j=jMin,jMax
454     DO i=iMin,iMax
455     guDiss(i,j) =
456     #ifdef OLD_UV_GEOM
457     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/
458     & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) )
459     #else
460     & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
461 jmc 1.39 & *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)
462 jmc 1.23 #endif
463 jmc 1.39 & *( ( fZon(i,j ) - fZon(i-1,j) )*AhDudxFac
464     & +( fMer(i,j+1) - fMer(i, j) )*AhDudyFac
465     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDudrFac
466     & *recip_rhoFacC(k)
467 jmc 1.23 & )
468     ENDDO
469 jmc 1.8 ENDDO
470    
471 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
472     IF ( useDiagnostics ) THEN
473     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Um',k,1,2,bi,bj,myThid)
474     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Um',k,1,2,bi,bj,myThid)
475     IF (.NOT.implicitViscosity)
476     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Um',k,1,2,bi,bj,myThid)
477     ENDIF
478     #endif
479    
480 jmc 1.37 C-- No-slip and drag BCs appear as body forces in cell abutting topography
481 jmc 1.23 IF (no_slip_sides) THEN
482 adcroft 1.1 C- No-slip BCs impose a drag at walls...
483 baylor 1.27 CALL MOM_U_SIDEDRAG(
484     I bi,bj,k,
485     I uFld, v4f, hFacZ,
486     I viscAh_Z,viscA4_Z,
487     I harmonic,biharmonic,useVariableViscosity,
488     O vF,
489     I myThid)
490 jmc 1.23 DO j=jMin,jMax
491     DO i=iMin,iMax
492     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
493     ENDDO
494     ENDDO
495     ENDIF
496 adcroft 1.1 C- No-slip BCs impose a drag at bottom
497 jmc 1.23 IF (bottomDragTerms) THEN
498     CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
499     DO j=jMin,jMax
500     DO i=iMin,iMax
501     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
502     ENDDO
503     ENDDO
504     ENDIF
505    
506 mlosch 1.32 #ifdef ALLOW_SHELFICE
507     IF (useShelfIce) THEN
508     CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
509     DO j=jMin,jMax
510     DO i=iMin,iMax
511     gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
512     ENDDO
513     ENDDO
514     ENDIF
515     #endif /* ALLOW_SHELFICE */
516    
517 jmc 1.23 C- endif momViscosity
518 adcroft 1.1 ENDIF
519    
520 jmc 1.12 C-- Forcing term (moved to timestep.F)
521     c IF (momForcing)
522     c & CALL EXTERNAL_FORCING_U(
523     c I iMin,iMax,jMin,jMax,bi,bj,k,
524     c I myTime,myThid)
525 adcroft 1.1
526     C-- Metric terms for curvilinear grid systems
527 adcroft 1.5 IF (useNHMTerms) THEN
528 jmc 1.33 C o Non-Hydrostatic (spherical) metric terms
529 adcroft 1.1 CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
530     DO j=jMin,jMax
531     DO i=iMin,iMax
532 jmc 1.33 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j)
533 adcroft 1.1 ENDDO
534     ENDDO
535 adcroft 1.5 ENDIF
536 jmc 1.33 IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
537     C o Spherical polar grid metric terms
538 adcroft 1.1 CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid)
539     DO j=jMin,jMax
540     DO i=iMin,iMax
541 jmc 1.33 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
542 adcroft 1.1 ENDDO
543     ENDDO
544 afe 1.20 ENDIF
545 jmc 1.33 IF ( usingCylindricalGrid .AND. metricTerms ) THEN
546     C o Cylindrical grid metric terms
547     CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
548     DO j=jMin,jMax
549     DO i=iMin,iMax
550     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j)
551     ENDDO
552 afe 1.19 ENDDO
553 adcroft 1.1 ENDIF
554    
555 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
556 adcroft 1.1
557     C---- Meridional momentum equation starts here
558    
559 jmc 1.23 IF (momAdvection) THEN
560     C--- Calculate mean fluxes (advection) between cells for meridional flow.
561     C Mean flow component of zonal flux -> fZon
562     CALL MOM_V_ADV_UV(bi,bj,k,uTrans,vFld,fZon,myThid)
563 adcroft 1.1
564     C-- Meridional flux (fMer is at north face of "v" cell)
565 jmc 1.23 C Mean flow component of meridional flux -> fMer
566     CALL MOM_V_ADV_VV(bi,bj,k,vTrans,vFld,fMer,myThid)
567 adcroft 1.1
568     C-- Vertical flux (fVer is at upper face of "v" cell)
569 jmc 1.23 C Mean flow component of vertical flux (at k+1) -> fVerV
570     CALL MOM_V_ADV_WV(
571     I bi,bj,k+1,vVel,wVel,rTransV,
572     O fVerV(1-OLx,1-OLy,kDown), myThid )
573 adcroft 1.1
574     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
575 jmc 1.23 DO j=jMin,jMax
576     DO i=iMin,iMax
577     gV(i,j,k,bi,bj) =
578 adcroft 1.1 #ifdef OLD_UV_GEOM
579 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
580     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
581 adcroft 1.1 #else
582 jmc 1.23 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
583 jmc 1.39 & *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
584 adcroft 1.1 #endif
585 jmc 1.39 & *( ( fZon(i+1,j) - fZon(i,j ) )*uDvdxFac
586     & +( fMer(i, j) - fMer(i,j-1) )*vDvdyFac
587     & +(fVerV(i,j,kDown) - fVerV(i,j,kUp))*rkSign*rVelDvdrFac
588 jmc 1.23 & )
589 jmc 1.24 ENDDO
590     ENDDO
591    
592     #ifdef ALLOW_DIAGNOSTICS
593     IF ( useDiagnostics ) THEN
594     CALL DIAGNOSTICS_FILL(fZon,'ADVx_Vm ',k,1,2,bi,bj,myThid)
595     CALL DIAGNOSTICS_FILL(fMer,'ADVy_Vm ',k,1,2,bi,bj,myThid)
596     CALL DIAGNOSTICS_FILL(fVerV(1-Olx,1-Oly,kUp),
597     & 'ADVrE_Vm',k,1,2,bi,bj,myThid)
598     ENDIF
599     #endif
600 adcroft 1.1
601 jmc 1.8 #ifdef NONLIN_FRSURF
602     C-- account for 3.D divergence of the flow in rStar coordinate:
603 heimbach 1.31 # ifndef DISABLE_RSTAR_CODE
604 jmc 1.23 IF ( select_rStar.GT.0 ) THEN
605     DO j=jMin,jMax
606     DO i=iMin,iMax
607     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
608 jmc 1.8 & - (rStarExpS(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
609     & *vVel(i,j,k,bi,bj)
610 jmc 1.23 ENDDO
611     ENDDO
612     ENDIF
613     IF ( select_rStar.LT.0 ) THEN
614     DO j=jMin,jMax
615     DO i=iMin,iMax
616     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
617     & - rStarDhSDt(i,j,bi,bj)*vVel(i,j,k,bi,bj)
618     ENDDO
619     ENDDO
620     ENDIF
621 heimbach 1.31 # endif /* DISABLE_RSTAR_CODE */
622 jmc 1.23 #endif /* NONLIN_FRSURF */
623    
624     ELSE
625     C- if momAdvection / else
626     DO j=1-OLy,sNy+OLy
627     DO i=1-OLx,sNx+OLx
628     gV(i,j,k,bi,bj) = 0. _d 0
629     ENDDO
630 jmc 1.8 ENDDO
631 jmc 1.23
632     C- endif momAdvection.
633 jmc 1.8 ENDIF
634 jmc 1.23
635     IF (momViscosity) THEN
636     C--- Calculate eddy fluxes (dissipation) between cells for meridional flow.
637     C Bi-harmonic term del^2 V -> v4F
638 jmc 1.37 IF (biharmonic)
639 jmc 1.23 & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid)
640    
641     C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon
642 baylor 1.25 CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,
643 baylor 1.27 I viscAh_Z,viscA4_Z,myThid)
644 jmc 1.23
645     C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
646 baylor 1.25 CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,
647 baylor 1.27 I viscAh_D,viscA4_D,myThid)
648 jmc 1.23
649     C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
650     IF (.NOT.implicitViscosity) THEN
651     CALL MOM_V_RVISCFLUX(bi,bj, k, vVel,KappaRV,fVrUp,myThid)
652     CALL MOM_V_RVISCFLUX(bi,bj,k+1,vVel,KappaRV,fVrDw,myThid)
653     ENDIF
654    
655     C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
656 jmc 1.39 C anelastic: hor.visc.fluxes are not scaled by rhoFac (by vert.visc.flx is)
657 jmc 1.23 DO j=jMin,jMax
658     DO i=iMin,iMax
659     gvDiss(i,j) =
660     #ifdef OLD_UV_GEOM
661     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/
662     & ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) )
663     #else
664     & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
665 jmc 1.39 & *recip_rAs(i,j,bi,bj)*recip_deepFac2C(k)
666 jmc 1.23 #endif
667 jmc 1.39 & *( ( fZon(i+1,j) - fZon(i,j ) )*AhDvdxFac
668     & +( fMer(i, j) - fMer(i,j-1) )*AhDvdyFac
669     & +( fVrDw(i,j) - fVrUp(i,j) )*rkSign*ArDvdrFac
670     & *recip_rhoFacC(k)
671 jmc 1.23 & )
672     ENDDO
673 jmc 1.8 ENDDO
674    
675 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
676     IF ( useDiagnostics ) THEN
677     CALL DIAGNOSTICS_FILL(fZon, 'VISCx_Vm',k,1,2,bi,bj,myThid)
678     CALL DIAGNOSTICS_FILL(fMer, 'VISCy_Vm',k,1,2,bi,bj,myThid)
679     IF (.NOT.implicitViscosity)
680     & CALL DIAGNOSTICS_FILL(fVrUp,'VISrE_Vm',k,1,2,bi,bj,myThid)
681     ENDIF
682     #endif
683    
684 jmc 1.37 C-- No-slip and drag BCs appear as body forces in cell abutting topography
685 mlosch 1.32 IF (no_slip_sides) THEN
686 adcroft 1.1 C- No-slip BCs impose a drag at walls...
687 baylor 1.27 CALL MOM_V_SIDEDRAG(
688     I bi,bj,k,
689     I vFld, v4f, hFacZ,
690     I viscAh_Z,viscA4_Z,
691     I harmonic,biharmonic,useVariableViscosity,
692     O vF,
693     I myThid)
694 jmc 1.23 DO j=jMin,jMax
695     DO i=iMin,iMax
696     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
697     ENDDO
698     ENDDO
699     ENDIF
700 adcroft 1.1 C- No-slip BCs impose a drag at bottom
701 jmc 1.23 IF (bottomDragTerms) THEN
702     CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
703     DO j=jMin,jMax
704     DO i=iMin,iMax
705     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
706     ENDDO
707     ENDDO
708     ENDIF
709    
710 mlosch 1.32 #ifdef ALLOW_SHELFICE
711     IF (useShelfIce) THEN
712     CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid)
713     DO j=jMin,jMax
714     DO i=iMin,iMax
715     gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
716     ENDDO
717     ENDDO
718     ENDIF
719     #endif /* ALLOW_SHELFICE */
720    
721 jmc 1.23 C- endif momViscosity
722 adcroft 1.1 ENDIF
723    
724 jmc 1.12 C-- Forcing term (moved to timestep.F)
725     c IF (momForcing)
726     c & CALL EXTERNAL_FORCING_V(
727     c I iMin,iMax,jMin,jMax,bi,bj,k,
728     c I myTime,myThid)
729 adcroft 1.1
730     C-- Metric terms for curvilinear grid systems
731 adcroft 1.5 IF (useNHMTerms) THEN
732 jmc 1.33 C o Non-Hydrostatic (spherical) metric terms
733 adcroft 1.1 CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
734     DO j=jMin,jMax
735     DO i=iMin,iMax
736 jmc 1.33 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j)
737 adcroft 1.1 ENDDO
738     ENDDO
739 adcroft 1.5 ENDIF
740 jmc 1.33 IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN
741     C o Spherical polar grid metric terms
742 adcroft 1.1 CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid)
743     DO j=jMin,jMax
744     DO i=iMin,iMax
745 jmc 1.33 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
746 adcroft 1.1 ENDDO
747     ENDDO
748     ENDIF
749 jmc 1.33 IF ( usingCylindricalGrid .AND. metricTerms ) THEN
750     C o Cylindrical grid metric terms
751     CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
752     DO j=jMin,jMax
753     DO i=iMin,iMax
754     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j)
755     ENDDO
756     ENDDO
757 afe 1.19 ENDIF
758 adcroft 1.1
759 jmc 1.23 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
760 adcroft 1.1
761     C-- Coriolis term
762     C Note. As coded here, coriolis will not work with "thin walls"
763 jmc 1.12 c IF (useCDscheme) THEN
764     c CALL MOM_CDSCHEME(bi,bj,k,dPhiHydX,dPhiHydY,myThid)
765     c ELSE
766     IF (.NOT.useCDscheme) THEN
767     CALL MOM_U_CORIOLIS(bi,bj,k,vFld,cf,myThid)
768     DO j=jMin,jMax
769     DO i=iMin,iMax
770     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
771     ENDDO
772     ENDDO
773 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
774     IF ( useDiagnostics )
775     & CALL DIAGNOSTICS_FILL(cf,'Um_Cori ',k,1,2,bi,bj,myThid)
776     #endif
777 jmc 1.12 CALL MOM_V_CORIOLIS(bi,bj,k,uFld,cf,myThid)
778     DO j=jMin,jMax
779     DO i=iMin,iMax
780     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
781     ENDDO
782     ENDDO
783 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
784     IF ( useDiagnostics )
785     & CALL DIAGNOSTICS_FILL(cf,'Vm_Cori ',k,1,2,bi,bj,myThid)
786     #endif
787 jmc 1.12 ENDIF
788    
789 jmc 1.34 C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w)
790 jmc 1.37 IF ( use3dCoriolis ) THEN
791 jmc 1.34 CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
792     DO j=jMin,jMax
793     DO i=iMin,iMax
794     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j)
795     ENDDO
796     ENDDO
797     IF ( usingCurvilinearGrid ) THEN
798     C- presently, non zero angleSinC array only supported with Curvilinear-Grid
799     CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid)
800     DO j=jMin,jMax
801     DO i=iMin,iMax
802     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j)
803     ENDDO
804 adcroft 1.6 ENDDO
805 jmc 1.34 ENDIF
806 adcroft 1.6 ENDIF
807 adcroft 1.1
808 jmc 1.23 C-- Set du/dt & dv/dt on boundaries to zero
809     DO j=jMin,jMax
810     DO i=iMin,iMax
811     gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
812     guDiss(i,j) = guDiss(i,j) *_maskW(i,j,k,bi,bj)
813     gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
814     gvDiss(i,j) = gvDiss(i,j) *_maskS(i,j,k,bi,bj)
815     ENDDO
816     ENDDO
817    
818 jmc 1.24 #ifdef ALLOW_DIAGNOSTICS
819     IF ( useDiagnostics ) THEN
820 baylor 1.28 CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid)
821 jmc 1.24 CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
822     & 'Um_Advec',k,1,2,bi,bj,myThid)
823     CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),
824     & 'Vm_Advec',k,1,2,bi,bj,myThid)
825     IF (momViscosity) THEN
826     CALL DIAGNOSTICS_FILL(guDiss,'Um_Diss ',k,1,2,bi,bj,myThid)
827     CALL DIAGNOSTICS_FILL(gvDiss,'Vm_Diss ',k,1,2,bi,bj,myThid)
828     ENDIF
829     ENDIF
830     #endif /* ALLOW_DIAGNOSTICS */
831    
832 adcroft 1.1 RETURN
833     END

  ViewVC Help
Powered by ViewVC 1.1.22