/[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.1 - (hide 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 adcroft 1.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