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

Contents of /MITgcm/pkg/mom_vecinv/mom_vecinv.F

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


Revision 1.1 - (show annotations) (download)
Thu Aug 16 17:16:03 2001 UTC (22 years, 9 months ago) by adcroft
Branch: MAIN
Added run-time control of vector-invariant/flux-form momentum eqns.

1 C $Header: $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 SUBROUTINE MOM_VECINV(
7 I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown,
8 I phi_hyd,KappaRU,KappaRV,
9 U fVerU, fVerV,
10 I myCurrentTime, myThid)
11 C /==========================================================\
12 C | S/R MOM_VECINV |
13 C | o Form the right hand-side of the momentum equation. |
14 C |==========================================================|
15 C | Terms are evaluated one layer at a time working from |
16 C | the bottom to the top. The vertically integrated |
17 C | barotropic flow tendency term is evluated by summing the |
18 C | tendencies. |
19 C | Notes: |
20 C | We have not sorted out an entirely satisfactory formula |
21 C | for the diffusion equation bc with lopping. The present |
22 C | form produces a diffusive flux that does not scale with |
23 C | open-area. Need to do something to solidfy this and to |
24 C | deal "properly" with thin walls. |
25 C \==========================================================/
26 IMPLICIT NONE
27
28 C == Global variables ==
29 #include "SIZE.h"
30 #include "DYNVARS.h"
31 #include "EEPARAMS.h"
32 #include "PARAMS.h"
33 #include "GRID.h"
34
35 C == Routine arguments ==
36 C fVerU - Flux of momentum in the vertical
37 C fVerV direction out of the upper face of a cell K
38 C ( flux into the cell above ).
39 C phi_hyd - Hydrostatic pressure
40 C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
41 C results will be set.
42 C kUp, kDown - Index for upper and lower layers.
43 C myThid - Instance number for this innvocation of CALC_MOM_RHS
44 _RL phi_hyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
45 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46 _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
47 _RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
48 _RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
49 INTEGER kUp,kDown
50 INTEGER myThid
51 _RL myCurrentTime
52 INTEGER bi,bj,iMin,iMax,jMin,jMax
53
54 C == Local variables ==
55 _RL aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56 _RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57 _RL vrF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58 _RL uCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59 _RL vCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60 _RL mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 _RL pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 _RL del2u(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63 _RL del2v(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL uDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL vDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 C I,J,K - Loop counters
77 INTEGER i,j,k
78 C rVelMaskOverride - Factor for imposing special surface boundary conditions
79 C ( set according to free-surface condition ).
80 C hFacROpen - Lopped cell factos used tohold fraction of open
81 C hFacRClosed and closed cell wall.
82 _RL rVelMaskOverride
83 C xxxFac - On-off tracer parameters used for switching terms off.
84 _RL uDudxFac
85 _RL AhDudxFac
86 _RL A4DuxxdxFac
87 _RL vDudyFac
88 _RL AhDudyFac
89 _RL A4DuyydyFac
90 _RL rVelDudrFac
91 _RL ArDudrFac
92 _RL fuFac
93 _RL phxFac
94 _RL mtFacU
95 _RL uDvdxFac
96 _RL AhDvdxFac
97 _RL A4DvxxdxFac
98 _RL vDvdyFac
99 _RL AhDvdyFac
100 _RL A4DvyydyFac
101 _RL rVelDvdrFac
102 _RL ArDvdrFac
103 _RL fvFac
104 _RL phyFac
105 _RL vForcFac
106 _RL mtFacV
107 INTEGER km1,kp1
108 _RL wVelBottomOverride
109 LOGICAL bottomDragTerms
110 _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111 _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113 _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114
115 km1=MAX(1,k-1)
116 kp1=MIN(Nr,k+1)
117 rVelMaskOverride=1.
118 IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac
119 wVelBottomOverride=1.
120 IF (k.EQ.Nr) wVelBottomOverride=0.
121
122 C Initialise intermediate terms
123 DO J=1-OLy,sNy+OLy
124 DO I=1-OLx,sNx+OLx
125 aF(i,j) = 0.
126 vF(i,j) = 0.
127 vrF(i,j) = 0.
128 uCf(i,j) = 0.
129 vCf(i,j) = 0.
130 mT(i,j) = 0.
131 pF(i,j) = 0.
132 del2u(i,j) = 0.
133 del2v(i,j) = 0.
134 dStar(i,j) = 0.
135 zStar(i,j) = 0.
136 uDiss(i,j) = 0.
137 vDiss(i,j) = 0.
138 vort3(i,j) = 0.
139 omega3(i,j) = 0.
140 ke(i,j) = 0.
141 ENDDO
142 ENDDO
143
144 C-- Term by term tracer parmeters
145 C o U momentum equation
146 uDudxFac = afFacMom*1.
147 AhDudxFac = vfFacMom*1.
148 A4DuxxdxFac = vfFacMom*1.
149 vDudyFac = afFacMom*1.
150 AhDudyFac = vfFacMom*1.
151 A4DuyydyFac = vfFacMom*1.
152 rVelDudrFac = afFacMom*1.
153 ArDudrFac = vfFacMom*1.
154 mTFacU = mtFacMom*1.
155 fuFac = cfFacMom*1.
156 phxFac = pfFacMom*1.
157 C o V momentum equation
158 uDvdxFac = afFacMom*1.
159 AhDvdxFac = vfFacMom*1.
160 A4DvxxdxFac = vfFacMom*1.
161 vDvdyFac = afFacMom*1.
162 AhDvdyFac = vfFacMom*1.
163 A4DvyydyFac = vfFacMom*1.
164 rVelDvdrFac = afFacMom*1.
165 ArDvdrFac = vfFacMom*1.
166 mTFacV = mtFacMom*1.
167 fvFac = cfFacMom*1.
168 phyFac = pfFacMom*1.
169 vForcFac = foFacMom*1.
170
171 IF ( no_slip_bottom
172 & .OR. bottomDragQuadratic.NE.0.
173 & .OR. bottomDragLinear.NE.0.) THEN
174 bottomDragTerms=.TRUE.
175 ELSE
176 bottomDragTerms=.FALSE.
177 ENDIF
178
179 C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP
180 IF (staggerTimeStep) THEN
181 phxFac = 0.
182 phyFac = 0.
183 ENDIF
184
185 C-- Calculate open water fraction at vorticity points
186 CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
187
188 C---- Calculate common quantities used in both U and V equations
189 C Calculate tracer cell face open areas
190 DO j=1-OLy,sNy+OLy
191 DO i=1-OLx,sNx+OLx
192 xA(i,j) = _dyG(i,j,bi,bj)
193 & *drF(k)*_hFacW(i,j,k,bi,bj)
194 yA(i,j) = _dxG(i,j,bi,bj)
195 & *drF(k)*_hFacS(i,j,k,bi,bj)
196 ENDDO
197 ENDDO
198
199 C Make local copies of horizontal flow field
200 DO j=1-OLy,sNy+OLy
201 DO i=1-OLx,sNx+OLx
202 uFld(i,j) = uVel(i,j,k,bi,bj)
203 vFld(i,j) = vVel(i,j,k,bi,bj)
204 ENDDO
205 ENDDO
206
207 C Calculate velocity field "volume transports" through tracer cell faces.
208 DO j=1-OLy,sNy+OLy
209 DO i=1-OLx,sNx+OLx
210 uTrans(i,j) = uFld(i,j)*xA(i,j)
211 vTrans(i,j) = vFld(i,j)*yA(i,j)
212 ENDDO
213 ENDDO
214
215 CALL MOM_VI_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
216
217 CALL MOM_VI_CALC_HDIV(bi,bj,k,uFld,vFld,hDiv,myThid)
218
219 CALL MOM_VI_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
220
221 CALL MOM_VI_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
222
223 IF (momViscosity) THEN
224 C Calculate del^2 u and del^2 v for bi-harmonic term
225 CALL MOM_VI_DEL2UV(
226 I bi,bj,k,hDiv,vort3,hFacZ,
227 O del2u,del2v,
228 & myThid)
229 CALL MOM_VI_CALC_HDIV(bi,bj,k,del2u,del2v,dStar,myThid)
230 CALL MOM_VI_CALC_RELVORT3(bi,bj,k,del2u,del2v,hFacZ,zStar,myThid)
231 C Calculate dissipation terms for U and V equations
232 CALL MOM_VI_HDISSIP(
233 I bi,bj,k,hDiv,vort3,hFacZ,dStar,zStar,
234 O uDiss,vDiss,
235 & myThid)
236 ENDIF
237
238 C---- Zonal momentum equation starts here
239
240 C-- Vertical flux (fVer is at upper face of "u" cell)
241
242 C Eddy component of vertical flux (interior component only) -> vrF
243 IF (momViscosity.AND..NOT.implicitViscosity)
244 & CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)
245
246 C Combine fluxes
247 DO j=jMin,jMax
248 DO i=iMin,iMax
249 fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)
250 ENDDO
251 ENDDO
252
253 C--- Hydrostatic term ( -1/rhoConst . dphi/dx )
254 IF (momPressureForcing) THEN
255 DO j=1-Olx,sNy+Oly
256 DO i=2-Olx,sNx+Olx
257 pf(i,j) = - _recip_dxC(i,j,bi,bj)
258 & *(phi_hyd(i,j,k)-phi_hyd(i-1,j,k))
259 ENDDO
260 ENDDO
261 ENDIF
262
263 C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
264 DO j=2-Oly,sNy+Oly-1
265 DO i=2-Olx,sNx+Olx-1
266 gU(i,j,k,bi,bj) = uDiss(i,j)
267 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
268 & *recip_rAw(i,j,bi,bj)
269 & *(
270 & +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
271 & )
272 & _PHM( +phxFac * pf(i,j) )
273 ENDDO
274 ENDDO
275
276 C-- No-slip and drag BCs appear as body forces in cell abutting topography
277 IF (momViscosity.AND.no_slip_sides) THEN
278 C- No-slip BCs impose a drag at walls...
279 CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,del2u,hFacZ,vF,myThid)
280 DO j=jMin,jMax
281 DO i=iMin,iMax
282 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
283 ENDDO
284 ENDDO
285 ENDIF
286 C- No-slip BCs impose a drag at bottom
287 IF (momViscosity.AND.bottomDragTerms) THEN
288 CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
289 DO j=jMin,jMax
290 DO i=iMin,iMax
291 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
292 ENDDO
293 ENDDO
294 ENDIF
295
296 C-- Forcing term
297 IF (momForcing)
298 & CALL EXTERNAL_FORCING_U(
299 I iMin,iMax,jMin,jMax,bi,bj,k,
300 I myCurrentTime,myThid)
301
302 C-- Metric terms for curvilinear grid systems
303 c IF (usingSphericalPolarMTerms) THEN
304 C o Spherical polar grid metric terms
305 c CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
306 c DO j=jMin,jMax
307 c DO i=iMin,iMax
308 c gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
309 c ENDDO
310 c ENDDO
311 c ENDIF
312
313 C-- Set du/dt on boundaries to zero
314 DO j=jMin,jMax
315 DO i=iMin,iMax
316 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
317 ENDDO
318 ENDDO
319
320
321 C---- Meridional momentum equation starts here
322
323 C-- Vertical flux (fVer is at upper face of "v" cell)
324
325 C Eddy component of vertical flux (interior component only) -> vrF
326 IF (momViscosity.AND..NOT.implicitViscosity)
327 & CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)
328
329 C Combine fluxes -> fVerV
330 DO j=jMin,jMax
331 DO i=iMin,iMax
332 fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)
333 ENDDO
334 ENDDO
335
336 C--- Hydorstatic term (-1/rhoConst . dphi/dy )
337 IF (momPressureForcing) THEN
338 DO j=jMin,jMax
339 DO i=iMin,iMax
340 pF(i,j) = -_recip_dyC(i,j,bi,bj)
341 & *(phi_hyd(i,j,k)-phi_hyd(i,j-1,k))
342 ENDDO
343 ENDDO
344 ENDIF
345
346 C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
347 DO j=jMin,jMax
348 DO i=iMin,iMax
349 gV(i,j,k,bi,bj) = vDiss(i,j)
350 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
351 & *recip_rAs(i,j,bi,bj)
352 & *(
353 & +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
354 & )
355 & _PHM( +phyFac*pf(i,j) )
356 ENDDO
357 ENDDO
358
359 C-- No-slip and drag BCs appear as body forces in cell abutting topography
360 IF (momViscosity.AND.no_slip_sides) THEN
361 C- No-slip BCs impose a drag at walls...
362 CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,del2v,hFacZ,vF,myThid)
363 DO j=jMin,jMax
364 DO i=iMin,iMax
365 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
366 ENDDO
367 ENDDO
368 ENDIF
369 C- No-slip BCs impose a drag at bottom
370 IF (momViscosity.AND.bottomDragTerms) THEN
371 CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
372 DO j=jMin,jMax
373 DO i=iMin,iMax
374 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
375 ENDDO
376 ENDDO
377 ENDIF
378
379 C-- Forcing term
380 IF (momForcing)
381 & CALL EXTERNAL_FORCING_V(
382 I iMin,iMax,jMin,jMax,bi,bj,k,
383 I myCurrentTime,myThid)
384
385 C-- Metric terms for curvilinear grid systems
386 c IF (usingSphericalPolarMTerms) THEN
387 C o Spherical polar grid metric terms
388 c CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
389 c DO j=jMin,jMax
390 c DO i=iMin,iMax
391 c gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
392 c ENDDO
393 c ENDDO
394 c ENDIF
395
396 C-- Set dv/dt on boundaries to zero
397 DO j=jMin,jMax
398 DO i=iMin,iMax
399 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
400 ENDDO
401 ENDDO
402
403 C-- Horizontal Coriolis terms
404 CALL MOM_VI_CORIOLIS(bi,bj,K,uFld,vFld,omega3,r_hFacZ,
405 & uCf,vCf,myThid)
406 DO j=jMin,jMax
407 DO i=iMin,iMax
408 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
409 & *_maskW(i,j,k,bi,bj)
410 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
411 & *_maskS(i,j,k,bi,bj)
412 ENDDO
413 ENDDO
414 c CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,r_hFacZ,uCf,myThid)
415 CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)
416 c CALL MOM_VI_U_CORIOLIS_C4(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)
417 DO j=jMin,jMax
418 DO i=iMin,iMax
419 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
420 & *_maskW(i,j,k,bi,bj)
421 ENDDO
422 ENDDO
423 c CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,r_hFacZ,vCf,myThid)
424 CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)
425 c CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)
426 DO j=jMin,jMax
427 DO i=iMin,iMax
428 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
429 & *_maskS(i,j,k,bi,bj)
430 ENDDO
431 ENDDO
432
433 IF (momAdvection) THEN
434 C-- Vertical shear terms (Coriolis)
435 CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)
436 DO j=jMin,jMax
437 DO i=iMin,iMax
438 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
439 & *_maskW(i,j,k,bi,bj)
440 ENDDO
441 ENDDO
442 CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)
443 DO j=jMin,jMax
444 DO i=iMin,iMax
445 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
446 & *_maskS(i,j,k,bi,bj)
447 ENDDO
448 ENDDO
449
450 C-- Bernoulli term
451 CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)
452 DO j=jMin,jMax
453 DO i=iMin,iMax
454 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
455 & *_maskW(i,j,k,bi,bj)
456 ENDDO
457 ENDDO
458 CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)
459 DO j=jMin,jMax
460 DO i=iMin,iMax
461 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
462 & *_maskS(i,j,k,bi,bj)
463 ENDDO
464 ENDDO
465 ENDIF
466
467 RETURN
468 END

  ViewVC Help
Powered by ViewVC 1.1.22