/[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.2 - (show annotations) (download)
Fri Aug 17 18:40:30 2001 UTC (22 years, 9 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint40pre8
Changes since 1.1: +34 -14 lines
Added method for dumping intermediate local arrays:
 mdsio_writetile - same as mdsio_writefield except works from inside bi,bj loop
 mdsio_writelocal - same as mdsio_writetile except works for local arrays
 write_local_r? - higher-level wrapper for mdsio_writelocal

Controlled by diagFreq. Defaults to zero (ie. no dumps)

Example given at end of mom_vecinv.F that dumps some local arrays.

1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/mom_vecinv/mom_vecinv.F,v 1.1 2001/08/16 17:16:03 adcroft Exp $
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, myIter, 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 _RL myCurrentTime
51 INTEGER myIter
52 INTEGER myThid
53 INTEGER bi,bj,iMin,iMax,jMin,jMax
54
55 C == Functions ==
56 LOGICAL DIFFERENT_MULTIPLE
57 EXTERNAL DIFFERENT_MULTIPLE
58
59 C == Local variables ==
60 _RL aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 _RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 _RL vrF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63 _RL uCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RL vCf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 _RL mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RL pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL del2u(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RL del2v(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 _RL uDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80 _RL vDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 C I,J,K - Loop counters
82 INTEGER i,j,k
83 C rVelMaskOverride - Factor for imposing special surface boundary conditions
84 C ( set according to free-surface condition ).
85 C hFacROpen - Lopped cell factos used tohold fraction of open
86 C hFacRClosed and closed cell wall.
87 _RL rVelMaskOverride
88 C xxxFac - On-off tracer parameters used for switching terms off.
89 _RL uDudxFac
90 _RL AhDudxFac
91 _RL A4DuxxdxFac
92 _RL vDudyFac
93 _RL AhDudyFac
94 _RL A4DuyydyFac
95 _RL rVelDudrFac
96 _RL ArDudrFac
97 _RL fuFac
98 _RL phxFac
99 _RL mtFacU
100 _RL uDvdxFac
101 _RL AhDvdxFac
102 _RL A4DvxxdxFac
103 _RL vDvdyFac
104 _RL AhDvdyFac
105 _RL A4DvyydyFac
106 _RL rVelDvdrFac
107 _RL ArDvdrFac
108 _RL fvFac
109 _RL phyFac
110 _RL vForcFac
111 _RL mtFacV
112 INTEGER km1,kp1
113 _RL wVelBottomOverride
114 LOGICAL bottomDragTerms
115 _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 _RL omega3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117 _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118 _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119
120 km1=MAX(1,k-1)
121 kp1=MIN(Nr,k+1)
122 rVelMaskOverride=1.
123 IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac
124 wVelBottomOverride=1.
125 IF (k.EQ.Nr) wVelBottomOverride=0.
126
127 C Initialise intermediate terms
128 DO J=1-OLy,sNy+OLy
129 DO I=1-OLx,sNx+OLx
130 aF(i,j) = 0.
131 vF(i,j) = 0.
132 vrF(i,j) = 0.
133 uCf(i,j) = 0.
134 vCf(i,j) = 0.
135 mT(i,j) = 0.
136 pF(i,j) = 0.
137 del2u(i,j) = 0.
138 del2v(i,j) = 0.
139 dStar(i,j) = 0.
140 zStar(i,j) = 0.
141 uDiss(i,j) = 0.
142 vDiss(i,j) = 0.
143 vort3(i,j) = 0.
144 omega3(i,j) = 0.
145 ke(i,j) = 0.
146 ENDDO
147 ENDDO
148
149 C-- Term by term tracer parmeters
150 C o U momentum equation
151 uDudxFac = afFacMom*1.
152 AhDudxFac = vfFacMom*1.
153 A4DuxxdxFac = vfFacMom*1.
154 vDudyFac = afFacMom*1.
155 AhDudyFac = vfFacMom*1.
156 A4DuyydyFac = vfFacMom*1.
157 rVelDudrFac = afFacMom*1.
158 ArDudrFac = vfFacMom*1.
159 mTFacU = mtFacMom*1.
160 fuFac = cfFacMom*1.
161 phxFac = pfFacMom*1.
162 C o V momentum equation
163 uDvdxFac = afFacMom*1.
164 AhDvdxFac = vfFacMom*1.
165 A4DvxxdxFac = vfFacMom*1.
166 vDvdyFac = afFacMom*1.
167 AhDvdyFac = vfFacMom*1.
168 A4DvyydyFac = vfFacMom*1.
169 rVelDvdrFac = afFacMom*1.
170 ArDvdrFac = vfFacMom*1.
171 mTFacV = mtFacMom*1.
172 fvFac = cfFacMom*1.
173 phyFac = pfFacMom*1.
174 vForcFac = foFacMom*1.
175
176 IF ( no_slip_bottom
177 & .OR. bottomDragQuadratic.NE.0.
178 & .OR. bottomDragLinear.NE.0.) THEN
179 bottomDragTerms=.TRUE.
180 ELSE
181 bottomDragTerms=.FALSE.
182 ENDIF
183
184 C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP
185 IF (staggerTimeStep) THEN
186 phxFac = 0.
187 phyFac = 0.
188 ENDIF
189
190 C-- Calculate open water fraction at vorticity points
191 CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
192
193 C---- Calculate common quantities used in both U and V equations
194 C Calculate tracer cell face open areas
195 DO j=1-OLy,sNy+OLy
196 DO i=1-OLx,sNx+OLx
197 xA(i,j) = _dyG(i,j,bi,bj)
198 & *drF(k)*_hFacW(i,j,k,bi,bj)
199 yA(i,j) = _dxG(i,j,bi,bj)
200 & *drF(k)*_hFacS(i,j,k,bi,bj)
201 ENDDO
202 ENDDO
203
204 C Make local copies of horizontal flow field
205 DO j=1-OLy,sNy+OLy
206 DO i=1-OLx,sNx+OLx
207 uFld(i,j) = uVel(i,j,k,bi,bj)
208 vFld(i,j) = vVel(i,j,k,bi,bj)
209 ENDDO
210 ENDDO
211
212 C Calculate velocity field "volume transports" through tracer cell faces.
213 DO j=1-OLy,sNy+OLy
214 DO i=1-OLx,sNx+OLx
215 uTrans(i,j) = uFld(i,j)*xA(i,j)
216 vTrans(i,j) = vFld(i,j)*yA(i,j)
217 ENDDO
218 ENDDO
219
220 CALL MOM_VI_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
221
222 CALL MOM_VI_CALC_HDIV(bi,bj,k,uFld,vFld,hDiv,myThid)
223
224 CALL MOM_VI_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
225
226 CALL MOM_VI_CALC_ABSVORT3(bi,bj,k,vort3,omega3,myThid)
227
228 IF (momViscosity) THEN
229 C Calculate del^2 u and del^2 v for bi-harmonic term
230 IF (viscA4.NE.0.) THEN
231 CALL MOM_VI_DEL2UV(bi,bj,k,hDiv,vort3,hFacZ,
232 O del2u,del2v,
233 & myThid)
234 CALL MOM_VI_CALC_HDIV(bi,bj,k,del2u,del2v,dStar,myThid)
235 CALL MOM_VI_CALC_RELVORT3(
236 & bi,bj,k,del2u,del2v,hFacZ,zStar,myThid)
237 ENDIF
238 C Calculate dissipation terms for U and V equations
239 C in terms of vorticity and divergence
240 IF (viscAh.NE.0. .OR. viscA4.NE.0.) THEN
241 CALL MOM_VI_HDISSIP(bi,bj,k,hDiv,vort3,hFacZ,dStar,zStar,
242 O uDiss,vDiss,
243 & myThid)
244 ENDIF
245 ENDIF
246
247 C---- Zonal momentum equation starts here
248
249 C-- Vertical flux (fVer is at upper face of "u" cell)
250
251 C Eddy component of vertical flux (interior component only) -> vrF
252 IF (momViscosity.AND..NOT.implicitViscosity)
253 & CALL MOM_U_RVISCFLUX(bi,bj,k,uVel,KappaRU,vrF,myThid)
254
255 C Combine fluxes
256 DO j=jMin,jMax
257 DO i=iMin,iMax
258 fVerU(i,j,kDown) = ArDudrFac*vrF(i,j)
259 ENDDO
260 ENDDO
261
262 C--- Hydrostatic term ( -1/rhoConst . dphi/dx )
263 IF (momPressureForcing) THEN
264 DO j=1-Olx,sNy+Oly
265 DO i=2-Olx,sNx+Olx
266 pf(i,j) = - _recip_dxC(i,j,bi,bj)
267 & *(phi_hyd(i,j,k)-phi_hyd(i-1,j,k))
268 ENDDO
269 ENDDO
270 ENDIF
271
272 C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
273 DO j=2-Oly,sNy+Oly-1
274 DO i=2-Olx,sNx+Olx-1
275 gU(i,j,k,bi,bj) = uDiss(i,j)
276 & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)
277 & *recip_rAw(i,j,bi,bj)
278 & *(
279 & +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac
280 & )
281 & _PHM( +phxFac * pf(i,j) )
282 ENDDO
283 ENDDO
284
285 C-- No-slip and drag BCs appear as body forces in cell abutting topography
286 IF (momViscosity.AND.no_slip_sides) THEN
287 C- No-slip BCs impose a drag at walls...
288 CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,del2u,hFacZ,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 C- No-slip BCs impose a drag at bottom
296 IF (momViscosity.AND.bottomDragTerms) THEN
297 CALL MOM_U_BOTTOMDRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid)
298 DO j=jMin,jMax
299 DO i=iMin,iMax
300 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+vF(i,j)
301 ENDDO
302 ENDDO
303 ENDIF
304
305 C-- Forcing term
306 IF (momForcing)
307 & CALL EXTERNAL_FORCING_U(
308 I iMin,iMax,jMin,jMax,bi,bj,k,
309 I myCurrentTime,myThid)
310
311 C-- Metric terms for curvilinear grid systems
312 c IF (usingSphericalPolarMTerms) THEN
313 C o Spherical polar grid metric terms
314 c CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid)
315 c DO j=jMin,jMax
316 c DO i=iMin,iMax
317 c gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
318 c ENDDO
319 c ENDDO
320 c ENDIF
321
322 C-- Set du/dt on boundaries to zero
323 DO j=jMin,jMax
324 DO i=iMin,iMax
325 gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj)
326 ENDDO
327 ENDDO
328
329
330 C---- Meridional momentum equation starts here
331
332 C-- Vertical flux (fVer is at upper face of "v" cell)
333
334 C Eddy component of vertical flux (interior component only) -> vrF
335 IF (momViscosity.AND..NOT.implicitViscosity)
336 & CALL MOM_V_RVISCFLUX(bi,bj,k,vVel,KappaRV,vrf,myThid)
337
338 C Combine fluxes -> fVerV
339 DO j=jMin,jMax
340 DO i=iMin,iMax
341 fVerV(i,j,kDown) = ArDvdrFac*vrF(i,j)
342 ENDDO
343 ENDDO
344
345 C--- Hydorstatic term (-1/rhoConst . dphi/dy )
346 IF (momPressureForcing) THEN
347 DO j=jMin,jMax
348 DO i=iMin,iMax
349 pF(i,j) = -_recip_dyC(i,j,bi,bj)
350 & *(phi_hyd(i,j,k)-phi_hyd(i,j-1,k))
351 ENDDO
352 ENDDO
353 ENDIF
354
355 C-- Tendency is minus divergence of the fluxes + coriolis + pressure term
356 DO j=jMin,jMax
357 DO i=iMin,iMax
358 gV(i,j,k,bi,bj) = vDiss(i,j)
359 & -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)
360 & *recip_rAs(i,j,bi,bj)
361 & *(
362 & +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac
363 & )
364 & _PHM( +phyFac*pf(i,j) )
365 ENDDO
366 ENDDO
367
368 C-- No-slip and drag BCs appear as body forces in cell abutting topography
369 IF (momViscosity.AND.no_slip_sides) THEN
370 C- No-slip BCs impose a drag at walls...
371 CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,del2v,hFacZ,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 C- No-slip BCs impose a drag at bottom
379 IF (momViscosity.AND.bottomDragTerms) THEN
380 CALL MOM_V_BOTTOMDRAG(bi,bj,k,vFld,KE,KappaRV,vF,myThid)
381 DO j=jMin,jMax
382 DO i=iMin,iMax
383 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vF(i,j)
384 ENDDO
385 ENDDO
386 ENDIF
387
388 C-- Forcing term
389 IF (momForcing)
390 & CALL EXTERNAL_FORCING_V(
391 I iMin,iMax,jMin,jMax,bi,bj,k,
392 I myCurrentTime,myThid)
393
394 C-- Metric terms for curvilinear grid systems
395 c IF (usingSphericalPolarMTerms) THEN
396 C o Spherical polar grid metric terms
397 c CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid)
398 c DO j=jMin,jMax
399 c DO i=iMin,iMax
400 c gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
401 c ENDDO
402 c ENDDO
403 c ENDIF
404
405 C-- Set dv/dt on boundaries to zero
406 DO j=jMin,jMax
407 DO i=iMin,iMax
408 gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)
409 ENDDO
410 ENDDO
411
412 C-- Horizontal Coriolis terms
413 CALL MOM_VI_CORIOLIS(bi,bj,K,uFld,vFld,omega3,r_hFacZ,
414 & uCf,vCf,myThid)
415 DO j=jMin,jMax
416 DO i=iMin,iMax
417 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
418 & *_maskW(i,j,k,bi,bj)
419 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
420 & *_maskS(i,j,k,bi,bj)
421 ENDDO
422 ENDDO
423 c CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,r_hFacZ,uCf,myThid)
424 CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)
425 c CALL MOM_VI_U_CORIOLIS_C4(bi,bj,K,vFld,vort3,r_hFacZ,uCf,myThid)
426 DO j=jMin,jMax
427 DO i=iMin,iMax
428 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
429 & *_maskW(i,j,k,bi,bj)
430 ENDDO
431 ENDDO
432 c CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,r_hFacZ,vCf,myThid)
433 CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)
434 c CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3,r_hFacZ,vCf,myThid)
435 DO j=jMin,jMax
436 DO i=iMin,iMax
437 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
438 & *_maskS(i,j,k,bi,bj)
439 ENDDO
440 ENDDO
441
442 IF (momAdvection) THEN
443 C-- Vertical shear terms (Coriolis)
444 CALL MOM_VI_U_VERTSHEAR(bi,bj,K,uVel,wVel,uCf,myThid)
445 DO j=jMin,jMax
446 DO i=iMin,iMax
447 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
448 & *_maskW(i,j,k,bi,bj)
449 ENDDO
450 ENDDO
451 CALL MOM_VI_V_VERTSHEAR(bi,bj,K,vVel,wVel,vCf,myThid)
452 DO j=jMin,jMax
453 DO i=iMin,iMax
454 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
455 & *_maskS(i,j,k,bi,bj)
456 ENDDO
457 ENDDO
458
459 C-- Bernoulli term
460 CALL MOM_VI_U_GRAD_KE(bi,bj,K,KE,uCf,myThid)
461 DO j=jMin,jMax
462 DO i=iMin,iMax
463 gU(i,j,k,bi,bj) = (gU(i,j,k,bi,bj)+uCf(i,j))
464 & *_maskW(i,j,k,bi,bj)
465 ENDDO
466 ENDDO
467 CALL MOM_VI_V_GRAD_KE(bi,bj,K,KE,vCf,myThid)
468 DO j=jMin,jMax
469 DO i=iMin,iMax
470 gV(i,j,k,bi,bj) = (gV(i,j,k,bi,bj)+vCf(i,j))
471 & *_maskS(i,j,k,bi,bj)
472 ENDDO
473 ENDDO
474 ENDIF
475
476 IF (
477 & DIFFERENT_MULTIPLE(diagFreq,myCurrentTime,
478 & myCurrentTime-deltaTClock)
479 & ) THEN
480 CALL WRITE_LOCAL_RL('Ph','I10',Nr,phi_hyd,bi,bj,1,myIter,myThid)
481 CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid)
482 CALL WRITE_LOCAL_RL('fU','I10',1,vCf,bi,bj,k,myIter,myThid)
483 CALL WRITE_LOCAL_RL('Du','I10',1,uDiss,bi,bj,k,myIter,myThid)
484 CALL WRITE_LOCAL_RL('Dv','I10',1,vDiss,bi,bj,k,myIter,myThid)
485 ENDIF
486
487 RETURN
488 END

  ViewVC Help
Powered by ViewVC 1.1.22