1 |
C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_mom_rhs.F,v 1.50 2001/05/29 14:01:36 adcroft Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "CPP_OPTIONS.h" |
5 |
|
6 |
#define COSINEMETH_III |
7 |
#undef ISOTROPIC_COS_SCALING |
8 |
|
9 |
C /==========================================================\ |
10 |
C | S/R CALC_MOM_RHS | |
11 |
C | o Form the right hand-side of the momentum equation. | |
12 |
C |==========================================================| |
13 |
C | Terms are evaluated one layer at a time working from | |
14 |
C | the bottom to the top. The vertically integrated | |
15 |
C | barotropic flow tendency term is evluated by summing the | |
16 |
C | tendencies. | |
17 |
C | Notes: | |
18 |
C | We have not sorted out an entirely satisfactory formula | |
19 |
C | for the diffusion equation bc with lopping. The present | |
20 |
C | form produces a diffusive flux that does not scale with | |
21 |
C | open-area. Need to do something to solidfy this and to | |
22 |
C | deal "properly" with thin walls. | |
23 |
C \==========================================================/ |
24 |
SUBROUTINE CALC_MOM_RHS( |
25 |
I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown, |
26 |
I phiHyd,KappaRU,KappaRV, |
27 |
U fVerU, fVerV, |
28 |
I myCurrentTime, myThid) |
29 |
|
30 |
IMPLICIT NONE |
31 |
|
32 |
C == Global variables == |
33 |
#include "SIZE.h" |
34 |
#include "DYNVARS.h" |
35 |
#include "FFIELDS.h" |
36 |
#include "EEPARAMS.h" |
37 |
#include "PARAMS.h" |
38 |
#include "GRID.h" |
39 |
#include "SURFACE.h" |
40 |
|
41 |
C == Routine arguments == |
42 |
C fZon - Work array for flux of momentum in the east-west |
43 |
C direction at the west face of a cell. |
44 |
C fMer - Work array for flux of momentum in the north-south |
45 |
C direction at the south face of a cell. |
46 |
C fVerU - Flux of momentum in the vertical |
47 |
C fVerV direction out of the upper face of a cell K |
48 |
C ( flux into the cell above ). |
49 |
C phiHyd - Hydrostatic Potential (pressure/rho) |
50 |
C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation |
51 |
C results will be set. |
52 |
C kUp, kDown - Index for upper and lower layers. |
53 |
C myThid - Instance number for this innvocation of CALC_MOM_RHS |
54 |
_RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
55 |
_RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
56 |
_RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
57 |
_RL fVerU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
58 |
_RL fVerV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) |
59 |
INTEGER kUp,kDown |
60 |
INTEGER myThid |
61 |
_RL myCurrentTime |
62 |
INTEGER bi,bj,iMin,iMax,jMin,jMax,jG |
63 |
|
64 |
C == Local variables == |
65 |
C ab15, ab05 - Weights for Adams-Bashforth time stepping scheme. |
66 |
C i,j,k - Loop counters |
67 |
C wMaskOverride - Land sea flag override for top layer. |
68 |
C afFacMom - Tracer parameters for turning terms |
69 |
C vfFacMom on and off. |
70 |
C pfFacMom afFacMom - Advective terms |
71 |
C cfFacMom vfFacMom - Eddy viscosity terms |
72 |
C mTFacMom pfFacMom - Pressure terms |
73 |
C cfFacMom - Coriolis terms |
74 |
C foFacMom - Forcing |
75 |
C mTFacMom - Metric term |
76 |
C vF - Temporary holding viscous term (Laplacian) |
77 |
C v4F - Temporary holding viscous term (Biharmonic) |
78 |
C cF - Temporary holding coriolis term. |
79 |
C mT - Temporary holding metric terms(s). |
80 |
C pF - Temporary holding pressure|potential gradient terms. |
81 |
C uDudxFac, AhDudxFac, etc ... individual term tracer parameters |
82 |
_RL aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
83 |
_RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
84 |
_RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
85 |
_RL cF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
86 |
_RL mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
87 |
_RL pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
88 |
_RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
89 |
_RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
90 |
_RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
91 |
_RS xA(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
92 |
_RS yA(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
93 |
_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
94 |
_RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
95 |
C I,J,K - Loop counters |
96 |
INTEGER i,j,k |
97 |
C rVelMaskOverride - Factor for imposing special surface boundary conditions |
98 |
C ( set according to free-surface condition ). |
99 |
C hFacROpen - Lopped cell factos used tohold fraction of open |
100 |
C hFacRClosed and closed cell wall. |
101 |
_RL rVelMaskOverride |
102 |
_RS hFacZOpen |
103 |
_RS hFacZClosedN |
104 |
_RS hFacZClosedS |
105 |
_RS hFacZClosedE |
106 |
_RS hFacZClosedW |
107 |
C xxxFac - On-off tracer parameters used for switching terms off. |
108 |
_RL uDudxFac |
109 |
_RL AhDudxFac |
110 |
_RL A4DuxxdxFac |
111 |
_RL vDudyFac |
112 |
_RL AhDudyFac |
113 |
_RL A4DuyydyFac |
114 |
_RL rVelDudrFac |
115 |
_RL ArDudrFac |
116 |
_RL fuFac |
117 |
_RL phxFac |
118 |
_RL mtFacU |
119 |
_RL uDvdxFac |
120 |
_RL AhDvdxFac |
121 |
_RL A4DvxxdxFac |
122 |
_RL vDvdyFac |
123 |
_RL AhDvdyFac |
124 |
_RL A4DvyydyFac |
125 |
_RL rVelDvdrFac |
126 |
_RL ArDvdrFac |
127 |
_RL fvFac |
128 |
_RL phyFac |
129 |
_RL vForcFac |
130 |
_RL mtFacV |
131 |
C ab05, ab15 - Adams-Bashforth time-stepping weights. |
132 |
_RL ab05, ab15 |
133 |
INTEGER km1,kp1 |
134 |
_RL maskDown,wVelBottomOverride,rdrckp1 |
135 |
LOGICAL bottomDragTerms |
136 |
_RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
137 |
|
138 |
km1=MAX(1,k-1) |
139 |
kp1=MIN(Nr,k+1) |
140 |
rVelMaskOverride=afFacMom |
141 |
IF ( k .EQ. 1 ) rVelMaskOverride=freeSurfFac*afFacMom |
142 |
wVelBottomOverride=1. |
143 |
IF (k.EQ.Nr) wVelBottomOverride=0. |
144 |
|
145 |
C Initialise intermediate terms |
146 |
DO J=1-OLy,sNy+OLy |
147 |
DO I=1-OLx,sNx+OLx |
148 |
aF(i,j) = 0. |
149 |
vF(i,j) = 0. |
150 |
v4F(i,j) = 0. |
151 |
cF(i,j) = 0. |
152 |
mT(i,j) = 0. |
153 |
pF(i,j) = 0. |
154 |
fZon(i,j) = 0. |
155 |
fMer(i,j) = 0. |
156 |
ENDDO |
157 |
ENDDO |
158 |
DO J=1-OLy,sNy+OLy-1 |
159 |
DO I=1-OLx,sNx+OLx-1 |
160 |
KE(i,j) = 0.25*( |
161 |
& uVel( i , j ,k,bi,bj)*uVel( i , j ,k,bi,bj) |
162 |
& +uVel(i+1, j ,k,bi,bj)*uVel(i+1, j ,k,bi,bj) |
163 |
& +vVel( i , j ,k,bi,bj)*vVel( i , j ,k,bi,bj) |
164 |
& +vVel( i ,j+1,k,bi,bj)*vVel( i ,j+1,k,bi,bj) ) |
165 |
ENDDO |
166 |
ENDDO |
167 |
C-- Term by term tracer parmeters |
168 |
C o U momentum equation |
169 |
uDudxFac = afFacMom*1. |
170 |
AhDudxFac = vfFacMom*1. |
171 |
A4DuxxdxFac = vfFacMom*1. |
172 |
vDudyFac = afFacMom*1. |
173 |
AhDudyFac = vfFacMom*1. |
174 |
A4DuyydyFac = vfFacMom*1. |
175 |
rVelDudrFac = afFacMom*1. |
176 |
ArDudrFac = vfFacMom*1. |
177 |
mTFacU = mtFacMom*1. |
178 |
fuFac = cfFacMom*1. |
179 |
phxFac = pfFacMom*1. |
180 |
C o V momentum equation |
181 |
uDvdxFac = afFacMom*1. |
182 |
AhDvdxFac = vfFacMom*1. |
183 |
A4DvxxdxFac = vfFacMom*1. |
184 |
vDvdyFac = afFacMom*1. |
185 |
AhDvdyFac = vfFacMom*1. |
186 |
A4DvyydyFac = vfFacMom*1. |
187 |
rVelDvdrFac = afFacMom*1. |
188 |
ArDvdrFac = vfFacMom*1. |
189 |
mTFacV = mtFacMom*1. |
190 |
fvFac = cfFacMom*1. |
191 |
phyFac = pfFacMom*1. |
192 |
vForcFac = foFacMom*1. |
193 |
|
194 |
IF (no_slip_bottom) THEN |
195 |
bottomDragTerms=.TRUE. |
196 |
ELSE |
197 |
bottomDragTerms=.FALSE. |
198 |
ENDIF |
199 |
|
200 |
C-- with stagger time stepping, grad Phi_Hyp is directly incoporated in TIMESTEP |
201 |
IF (staggerTimeStep) THEN |
202 |
phxFac = 0. |
203 |
phyFac = 0. |
204 |
ENDIF |
205 |
|
206 |
C-- Adams-Bashforth weighting factors |
207 |
ab15 = 1.5 _d 0 + abEps |
208 |
ab05 = -0.5 _d 0 - abEps |
209 |
|
210 |
C-- Calculate open water fraction at vorticity points |
211 |
DO i=1-Olx,sNx+Olx |
212 |
hFacZ(i,1-Oly)=0. |
213 |
ENDDO |
214 |
DO j=2-Oly,sNy+Oly |
215 |
hFacZ(1-Olx,j)=0. |
216 |
DO i=2-Olx,sNx+Olx |
217 |
hFacZOpen=min(_hFacW(i,j,k,bi,bj), |
218 |
& _hFacW(i,j-1,k,bi,bj)) |
219 |
hFacZOpen=min(_hFacS(i,j,k,bi,bj),hFacZOpen) |
220 |
hFacZOpen=min(_hFacS(i-1,j,k,bi,bj),hFacZOpen) |
221 |
hFacZ(i,j)=hFacZOpen |
222 |
ENDDO |
223 |
ENDDO |
224 |
|
225 |
C---- Calculate common quantities used in both U and V equations |
226 |
C Calculate tracer cell face open areas |
227 |
DO j=1-OLy,sNy+OLy |
228 |
DO i=1-OLx,sNx+OLx |
229 |
xA(i,j) = _dyG(i,j,bi,bj) |
230 |
& *drF(k)*_hFacW(i,j,k,bi,bj) |
231 |
yA(i,j) = _dxG(i,j,bi,bj) |
232 |
& *drF(k)*_hFacS(i,j,k,bi,bj) |
233 |
ENDDO |
234 |
ENDDO |
235 |
C Calculate velocity field "volume transports" through tracer cell faces. |
236 |
DO j=1-OLy,sNy+OLy |
237 |
DO i=1-OLx,sNx+OLx |
238 |
uTrans(i,j) = uVel(i,j,k,bi,bj)*xA(i,j) |
239 |
vTrans(i,j) = vVel(i,j,k,bi,bj)*yA(i,j) |
240 |
ENDDO |
241 |
ENDDO |
242 |
|
243 |
C---- Zonal momentum equation starts here |
244 |
|
245 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
246 |
C-- del^2 U (for use in del^4 U) |
247 |
C Zonal flux d/dx U |
248 |
DO j=1-Oly,sNy+Oly |
249 |
DO i=1-Olx,sNx+Olx-1 |
250 |
fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj) |
251 |
& *_dyF(i,j,bi,bj) |
252 |
& *_recip_dxF(i,j,bi,bj) |
253 |
& *(uVel(i+1,j,k,bi,bj)-uVel(i,j,k,bi,bj)) |
254 |
ENDDO |
255 |
ENDDO |
256 |
C Meridional flux d/dy U |
257 |
DO j=1-Oly+1,sNy+Oly |
258 |
DO i=1-Olx,sNx+Olx |
259 |
fMer(i,j) = drF(k)*hFacZ(i,j) |
260 |
& *_dxV(i,j,bi,bj) |
261 |
& *_recip_dyU(i,j,bi,bj) |
262 |
& *(uVel(i,j,k,bi,bj)-uVel(i,j-1,k,bi,bj)) |
263 |
ENDDO |
264 |
ENDDO |
265 |
C del^2 U |
266 |
DO j=0,sNy+1 |
267 |
DO i=0,sNx+2 |
268 |
v4F(i,j) = recip_drF(k)*_recip_hFacW(i,j,k,bi,bj) |
269 |
& *recip_rAw(i,j,bi,bj) |
270 |
& *( (fZon(i,j ) - fZon(i-1,j)) |
271 |
#ifdef COSINEMETH_III |
272 |
& *sqcosFacU(J,bi,bj) |
273 |
#endif |
274 |
& +fMer(i,j+1) - fMer(i ,j) |
275 |
& )*_maskW(i,j,k,bi,bj) |
276 |
ENDDO |
277 |
ENDDO |
278 |
IF (no_slip_sides) THEN |
279 |
C-- No-slip BCs impose a drag at walls... |
280 |
DO j=0,sNy+1 |
281 |
DO i=0,sNx+2 |
282 |
hFacZClosedS = _hFacW(i,j,k,bi,bj) - hFacZ(i,j) |
283 |
hFacZClosedN = _hFacW(i,j,k,bi,bj) - hFacZ(i,j+1) |
284 |
v4F(i,j) = v4F(i,j) |
285 |
& - _recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
286 |
& *recip_rAw(i,j,bi,bj) |
287 |
& *( hFacZClosedS*dxV(i, j ,bi,bj) |
288 |
& *_recip_dyU(i, j ,bi,bj) |
289 |
& +hFacZClosedN*dxV(i,j+1,bi,bj) |
290 |
& *_recip_dyU(i,j+1,bi,bj) |
291 |
& )*drF(k)*2.*uVel(i,j,k,bi,bj) |
292 |
& *_maskW(i,j,k,bi,bj) |
293 |
ENDDO |
294 |
ENDDO |
295 |
ENDIF |
296 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
297 |
|
298 |
C--- Calculate mean and eddy fluxes between cells for zonal flow. |
299 |
C-- Zonal flux (fZon is at east face of "u" cell) |
300 |
#define _XAatP( i,j,k,bi,bj) ( XA(i,j)+XA(i+1,j) )* 0.5 _d 0 |
301 |
#define _YAatUBY( i,j,k,bi,bj) ( YA(i,j)+YA(i-1,j) )* 0.5 _d 0 |
302 |
#define _RAatU1( i,j,k,bi,bj) ( rA(i,j,bi,bj)*maskC(i,j) |
303 |
#define _RAatU2( i,j,k,bi,bj) + rA(i-1,j,bi,bj) |
304 |
#define _RAatU3( i,j,k,bi,bj) *maskC(i-1,j) ) * 0.5 _d 0 |
305 |
#define _RAatU_Bot(i,j,k,bi,bj) ( rA(i,j,bi,bj) + rA(i-1,j,bi,bj) ) * 0.5 _d 0 |
306 |
C o Mean flow component of zonal flux |
307 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
308 |
DO j=jMin,jMax |
309 |
DO i=iMin,iMax |
310 |
af(i,j) = |
311 |
& (uTrans(i,j)+uTrans(i+1,j) ) |
312 |
& *(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))*0.25 _d 0 |
313 |
ENDDO |
314 |
ENDDO |
315 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
316 |
C o Eddy component of zonal flux |
317 |
C - Laplacian and bi-harmonic terms |
318 |
DO j=jMin,jMax |
319 |
DO i=iMin,iMax |
320 |
vF(i,j) = |
321 |
#ifdef OLD_UV_GEOMETRY |
322 |
& ( XA(i,j)+XA(i+1,j) )* 0.5 _d 0 |
323 |
#else |
324 |
& _dyF(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj) |
325 |
#endif /* OLD_UV_GEOMETRY */ |
326 |
& *( |
327 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
328 |
& -viscAh*(uVel(i+1,j,k,bi,bj)-uVel(i,j,k,bi,bj)) |
329 |
& *cosFacU(J,bi,bj) |
330 |
#else |
331 |
& 0. |
332 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
333 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
334 |
& +viscA4*( v4F(i+1,j) -v4F(i,j) ) |
335 |
#ifdef COSINEMETH_III |
336 |
& *sqcosFacU(J,bi,bj) |
337 |
#else |
338 |
& *cosFacU(J,bi,bj) |
339 |
#endif |
340 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
341 |
& )*_recip_dxF(i,j,bi,bj) |
342 |
|
343 |
ENDDO |
344 |
ENDDO |
345 |
DO j=jMin,jMax |
346 |
DO i=iMin,iMax |
347 |
fZon(i,j) = 0. |
348 |
& _ADM( + uDudxFac * aF (i,j) ) |
349 |
& + AhDudxFac * vF (i,j) |
350 |
ENDDO |
351 |
ENDDO |
352 |
|
353 |
C-- Meridional flux (fMer is at south face of "u" cell) |
354 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
355 |
C o Mean flow component of meridional flux |
356 |
DO j=jMin,jMax |
357 |
DO i=iMin,iMax |
358 |
af(i,j) = |
359 |
& (vTrans(i,j)+vTrans(i-1,j)) |
360 |
& *(uVel(i,j,k,bi,bj)+uVel(i,j-1,k,bi,bj))*0.25 _d 0 |
361 |
#ifdef OLD_ADV_BCS |
362 |
& *_maskW(i,j,k,bi,bj) |
363 |
& *_maskW(i,j-1,k,bi,bj) |
364 |
#endif /* OLD_ADV_BCS */ |
365 |
C Note!!!! The line "*maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)" is |
366 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
367 |
C Note!!!! FV paper used a different bc in which this line is |
368 |
C Note!!!! omitted. |
369 |
ENDDO |
370 |
ENDDO |
371 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
372 |
C Eddy component of meridional flux |
373 |
C o Laplacian and bi-harmonic term |
374 |
DO j=jMin,jMax |
375 |
DO i=iMin,iMax |
376 |
vF(i,j) = _dxV(i,j,bi,bj)*drF(k)*hFacZ(i,j) |
377 |
& *( |
378 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
379 |
#ifdef ISOTROPIC_COS_SCALING |
380 |
& -viscAh*(uVel(i,j,k,bi,bj)-uVel(i,j-1,k,bi,bj)) |
381 |
& *cosFacV(J,bi,bj) |
382 |
#else |
383 |
& -viscAh*(uVel(i,j,k,bi,bj)-uVel(i,j-1,k,bi,bj)) |
384 |
#endif |
385 |
#else |
386 |
& 0. |
387 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
388 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
389 |
#ifdef ISOTROPIC_COS_SCALING |
390 |
& +viscA4*( v4F(i,j) -v4F(i,j-1) ) |
391 |
& *cosFacV(J,bi,bj) |
392 |
#else |
393 |
& +viscA4*( v4F(i,j) -v4F(i,j-1) ) |
394 |
#endif |
395 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
396 |
& ) |
397 |
& *_recip_dyU(i,j,bi,bj) |
398 |
ENDDO |
399 |
ENDDO |
400 |
DO j=jMin,jMax |
401 |
DO i=iMin,iMax |
402 |
fMer(i,j) = 0. |
403 |
& _ADM( + vDudyFac * aF (i,j) ) |
404 |
& + AhDudyFac * vF (i,j) |
405 |
ENDDO |
406 |
ENDDO |
407 |
C-- Vertical flux (fVer is at upper face of "u" cell) |
408 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
409 |
C-- Free surface correction temr |
410 |
IF (k.EQ.1) THEN |
411 |
DO j=jMin,jMax |
412 |
DO i=iMin,iMax |
413 |
fVerU(i,j,kUp)=rVelMaskOverride*( |
414 |
& wVel( i ,j,k,bi,bj)*rA( i ,j,bi,bj) |
415 |
& +wVel(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj) |
416 |
& )*0.5 |
417 |
& *uVel(i,j,k,bi,bj) |
418 |
ENDDO |
419 |
ENDDO |
420 |
ENDIF |
421 |
C Mean flow component of vertical flux |
422 |
DO j=jMin,jMax |
423 |
DO i=iMin,iMax |
424 |
af(i,j) = |
425 |
& wVelBottomOverride*( |
426 |
& wVel( i ,j,kp1,bi,bj)*rA( i ,j,bi,bj) |
427 |
& +wVel(i-1,j,kp1,bi,bj)*rA(i-1,j,bi,bj) |
428 |
& ) |
429 |
& *( uVel(i,j,kp1,bi,bj)+uVel(i,j,k,bi,bj) ) |
430 |
& *0.25 _d 0 |
431 |
#ifdef OLD_ADV_BCS |
432 |
& *_maskW(i,j,kp1,bi,bj) |
433 |
& *_maskW(i,j,k,bi,bj) |
434 |
#endif /* OLD_ADV_BCS */ |
435 |
C Note!!!! The line "*(maskW(i,j,k,bi,bj)*maskW(i,j,kM1,bi,bj))" is |
436 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
437 |
C Note!!!! FV paper used a different bc in which this line is |
438 |
C Note!!!! omitted. |
439 |
ENDDO |
440 |
ENDDO |
441 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
442 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
443 |
#ifdef OLD_UV_GEOMETRY |
444 |
#define _RA_AT_W (RA(i,j,bi,bj)*maskC(i,j)+RA(i-1,j,bi,bj)*maskC(i-1,j))*.5 |
445 |
#else |
446 |
#define _RA_AT_W rAw(i,j,bi,bj) |
447 |
#endif /* OLD_UV_GEOMETRY */ |
448 |
C Eddy component of vertical flux (interior component only) |
449 |
IF (.NOT.implicitViscosity) THEN |
450 |
DO j=jMin,jMax |
451 |
DO i=iMin,iMax |
452 |
vf(i,j) = |
453 |
& -KappaRU(i,j,kp1) |
454 |
& *_RA_AT_W |
455 |
& *( uVel(i,j,k,bi,bj)-uVel(i,j,kp1,bi,bj) |
456 |
& )*rkFac*recip_drC(kp1) |
457 |
& *_maskW(i,j,kp1,bi,bj) |
458 |
ENDDO |
459 |
ENDDO |
460 |
ENDIF |
461 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
462 |
IF (implicitViscosity) THEN |
463 |
DO j=jMin,jMax |
464 |
DO i=iMin,iMax |
465 |
fVerU(i,j,kDown) = 0. |
466 |
& _ADM( + rVelDudrFac * af(i,j) ) |
467 |
ENDDO |
468 |
ENDDO |
469 |
ELSE |
470 |
DO j=jMin,jMax |
471 |
DO i=iMin,iMax |
472 |
fVerU(i,j,kDown) = 0. |
473 |
& _ADM( + rVelDudrFac * af(i,j) ) |
474 |
& _LPM( + ArDudrFac * vf(i,j) ) |
475 |
ENDDO |
476 |
ENDDO |
477 |
ENDIF |
478 |
|
479 |
#ifdef INCLUDE_GRADPH_CODE |
480 |
C--- Hydrostatic term ( -1/rhoConst . dphi/dx ) |
481 |
DO j=jMin,jMax |
482 |
DO i=iMin,iMax |
483 |
pf(i,j) = - _recip_dxC(i,j,bi,bj) |
484 |
& *(phiHyd(i,j,k)-phiHyd(i-1,j,k)) |
485 |
ENDDO |
486 |
ENDDO |
487 |
#endif /* INCLUDE_GRADPH_CODE */ |
488 |
|
489 |
C-- Tendency is minus divergence of the fluxes + coriolis + pressure |
490 |
C-- term |
491 |
DO j=jMin,jMax |
492 |
DO i=iMin,iMax |
493 |
gU(i,j,k,bi,bj) = |
494 |
#ifdef OLD_UV_GEOM |
495 |
& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/ |
496 |
& ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) ) |
497 |
#else |
498 |
& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
499 |
& *recip_rAw(i,j,bi,bj) |
500 |
#endif |
501 |
& *(fZon(i,j ) - fZon(i-1,j) |
502 |
& +fMer(i,j+1) - fMer(i ,j) |
503 |
& +fVerU(i,j,kUp)*rkFac - fVerU(i,j,kDown)*rkFac |
504 |
& ) |
505 |
& _PHM( +phxFac * pf(i,j) ) |
506 |
ENDDO |
507 |
ENDDO |
508 |
|
509 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
510 |
C-- No-slip and drag BCs appear as body forces in cell abutting topography |
511 |
IF (no_slip_sides) THEN |
512 |
C- No-slip BCs impose a drag at walls... |
513 |
DO j=jMin,jMax |
514 |
DO i=iMin,iMax |
515 |
hFacZClosedS = _hFacW(i,j,k,bi,bj) - hFacZ(i,j) |
516 |
hFacZClosedN = _hFacW(i,j,k,bi,bj) - hFacZ(i,j+1) |
517 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) |
518 |
& -_recip_hFacW(i,j,k,bi,bj) |
519 |
& *recip_drF(k)*recip_rAw(i,j,bi,bj) |
520 |
& *( hFacZClosedS*_dxV(i, j ,bi,bj) |
521 |
& *_recip_dyU(i, j ,bi,bj) |
522 |
& +hFacZClosedN*_dxV(i,j+1,bi,bj) |
523 |
& *_recip_dyU(i,j+1,bi,bj) ) |
524 |
& *drF(k)*2.*( viscAh*uVel(i,j,k,bi,bj) |
525 |
#ifdef ISOTROPIC_COS_SCALING |
526 |
& -viscA4*v4F(i,j) ) *cosFacU(J,bi,bj) |
527 |
#else |
528 |
& -viscA4*v4F(i,j) ) |
529 |
#endif |
530 |
ENDDO |
531 |
ENDDO |
532 |
ENDIF |
533 |
C- No-slip BCs impose a drag at bottom |
534 |
IF (bottomDragTerms) THEN |
535 |
rdrckp1=recip_drC(kp1) |
536 |
IF (k.EQ.Nr) rdrckp1=recip_drF(k) |
537 |
DO j=jMin,jMax |
538 |
DO i=iMin,iMax |
539 |
maskDown=_maskW(i,j,kp1,bi,bj) |
540 |
IF (k.EQ.Nr) maskDown=0. _d 0 |
541 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) |
542 |
& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
543 |
& *( |
544 |
& 2.*KappaRU(i,j,kp1)*rkFac*rdrckp1 |
545 |
& + bottomDragLinear |
546 |
& )*(1.-maskDown)*uVel(i,j,k,bi,bj) |
547 |
IF (KE(i,j)+KE(i-1,j) .NE. 0.) |
548 |
& gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) |
549 |
& -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) |
550 |
& *bottomDragQuadratic*sqrt(KE(i,j)+KE(i-1,j)) |
551 |
& *(1.-maskDown)*uVel(i,j,k,bi,bj) |
552 |
ENDDO |
553 |
ENDDO |
554 |
ENDIF |
555 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
556 |
|
557 |
#ifdef INCLUDE_MOMENTUM_FORCING_CODE |
558 |
C-- Forcing term |
559 |
CALL EXTERNAL_FORCING_U( |
560 |
I iMin,iMax,jMin,jMax,bi,bj,k, |
561 |
I myCurrentTime,myThid) |
562 |
#endif /* INCLUDE_MOMENTUM_FORCING_CODE */ |
563 |
|
564 |
#ifdef INCLUDE_MOMENTUM_METRIC_TERM_CODE |
565 |
C-- Metric terms for curvilinear grid systems |
566 |
IF ( usingSphericalPolarMTerms ) THEN |
567 |
C o Spherical polar grid metric terms |
568 |
DO j=jMin,jMax |
569 |
DO i=iMin,iMax |
570 |
mT(i,j) = -1.* uVel(i,j,k,bi,bj)*recip_RSphere |
571 |
& *0.25 _d 0*(wVelBottomOverride* |
572 |
& (wVel(i-1,j,kp1,bi,bj)+wVel(i,j,kp1,bi,bj)) |
573 |
& +wVel(i-1,j, k ,bi,bj)+wVel(i,j, k ,bi,bj) |
574 |
& )*rkFac*recip_horiVertRatio |
575 |
ENDDO |
576 |
ENDDO |
577 |
DO j=jMin,jMax |
578 |
DO i=iMin,iMax |
579 |
mT(i,j) = mT(i,j)+uVel(i,j,k,bi,bj)*recip_RSphere |
580 |
& *0.25 _d 0*( vVel(i,j ,k,bi,bj)+vVel(i-1,j ,k,bi,bj) |
581 |
& +vVel(i,j+1,k,bi,bj)+vVel(i-1,j+1,k,bi,bj) |
582 |
& )*_tanPhiAtU(i,j,bi,bj) |
583 |
ENDDO |
584 |
ENDDO |
585 |
DO j=jMin,jMax |
586 |
DO i=iMin,iMax |
587 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+ |
588 |
& mTFacU*mT(i,j) |
589 |
ENDDO |
590 |
ENDDO |
591 |
ENDIF |
592 |
#endif /* INCLUDE_MOMENTUM_METRIC_TERM_CODE */ |
593 |
|
594 |
C-- Set du/dt on boundaries to zero |
595 |
DO j=jMin,jMax |
596 |
DO i=iMin,iMax |
597 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj) |
598 |
ENDDO |
599 |
ENDDO |
600 |
|
601 |
C---- Meridional momentum equation starts here |
602 |
|
603 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
604 |
C-- del^2 V (for use in del^4 V) |
605 |
C Zonal flux d/dx V |
606 |
DO j=1-Oly,sNy+Oly |
607 |
DO i=1-Olx+1,sNx+Olx |
608 |
fZon(i,j) = drF(k)*hFacZ(i,j) |
609 |
& *_dyU(i,j,bi,bj) |
610 |
& *_recip_dxV(i,j,bi,bj) |
611 |
& *(vVel(i,j,k,bi,bj)-vVel(i-1,j,k,bi,bj)) |
612 |
ENDDO |
613 |
ENDDO |
614 |
C Meridional flux d/dy V |
615 |
DO j=1-Oly,sNy+Oly-1 |
616 |
DO i=1-Olx,sNx+Olx |
617 |
fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj) |
618 |
& *_dxF(i,j,bi,bj) |
619 |
& *_recip_dyF(i,j,bi,bj) |
620 |
& *(vVel(i,j+1,k,bi,bj)-vVel(i,j,k,bi,bj)) |
621 |
ENDDO |
622 |
ENDDO |
623 |
C del^2 U |
624 |
DO j=0,sNy+2 |
625 |
DO i=0,sNx+1 |
626 |
v4F(i,j) = recip_drF(k)*_recip_hFacS(i,j,k,bi,bj) |
627 |
& *recip_rAs(i,j,bi,bj) |
628 |
& *( (fZon(i+1,j) - fZon(i, j )) |
629 |
#ifdef COSINEMETH_III |
630 |
& *sqcosFacV(J,bi,bj) |
631 |
#endif |
632 |
& +fMer( i ,j) - fMer(i,j-1) |
633 |
& )*_maskS(i,j,k,bi,bj) |
634 |
ENDDO |
635 |
ENDDO |
636 |
IF (no_slip_sides) THEN |
637 |
C-- No-slip BCs impose a drag at walls... |
638 |
DO j=0,sNy+2 |
639 |
DO i=0,sNx+1 |
640 |
hFacZClosedW = _hFacS(i,j,k,bi,bj) - hFacZ(i,j) |
641 |
hFacZClosedE = _hFacS(i,j,k,bi,bj) - hFacZ(i+1,j) |
642 |
v4F(i,j) = v4F(i,j) |
643 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k) |
644 |
& *recip_rAs(i,j,bi,bj) |
645 |
& *( hFacZClosedW*dyU( i ,j,bi,bj) |
646 |
& *_recip_dxV( i ,j,bi,bj) |
647 |
& +hFacZClosedE*dyU(i+1,j,bi,bj) |
648 |
& *_recip_dxV(i+1,j,bi,bj) |
649 |
& )*drF(k)*2.*vVel(i,j,k,bi,bj) |
650 |
& *_maskS(i,j,k,bi,bj) |
651 |
#ifdef COSINEMETH_III |
652 |
& *sqcosFacV(J,bi,bj) |
653 |
#endif |
654 |
ENDDO |
655 |
ENDDO |
656 |
ENDIF |
657 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
658 |
|
659 |
C--- Calculate mean and eddy fluxes between cells for meridional |
660 |
C--- velocity. Zonal flux (fZon is at west face of "v" cell) |
661 |
#define _XAatVBX( i,j,k,bi,bj) ( XA(i,j)+XA(i,j-1) )* 0.5 _d 0 |
662 |
#define _YAatP( i,j,k,bi,bj) ( YA(i,j) + YA(i,j+1) ) *0.5 _d 0 |
663 |
#define _RAatV1( i,j,k,bi,bj) ( _rA(i,j,bi,bj)*maskC(i,j) |
664 |
#define _RAatV2( i,j,k,bi,bj) + _rA(i,j-1,bi,bj) |
665 |
#define _RAatV3( i,j,k,bi,bj) *maskC(i,j-1) ) * 0.5 _d 0 |
666 |
#define _RAatV_Bot1(i,j,k,bi,bj) ( _rA(i,j,bi,bj) |
667 |
#define _RAatV_Bot2(i,j,k,bi,bj) +_rA(i,j-1,bi,bj))*0.5 _d 0 |
668 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
669 |
C o Mean flow component of zonal flux |
670 |
DO j=jMin,jMax |
671 |
DO i=iMin,iMax |
672 |
af(i,j) = |
673 |
& (uTrans(i,j)+uTrans(i,j-1) ) |
674 |
& *(vVel(i,j,k,bi,bj)+vVel(i-1,j,k,bi,bj))*0.25 _d 0 |
675 |
#ifdef OLD_ADV_BCS |
676 |
& *_maskS(i,j,k,bi,bj)*_maskS(i-1,j,k,bi,bj) |
677 |
#endif /* OLD_ADV_BCS */ |
678 |
C Note!!!! The line "*maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)" is |
679 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
680 |
C Note!!!! FV paper used a different bc in which this line is |
681 |
C Note!!!! omitted. |
682 |
ENDDO |
683 |
ENDDO |
684 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
685 |
C Eddy component of zonal flux |
686 |
C o Laplacian term |
687 |
DO j=jMin,jMax |
688 |
DO i=iMin,iMax |
689 |
vf(i,j) = _dyU(i,j,bi,bj)*drF(K)*hFacZ(i,j) |
690 |
& *( |
691 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
692 |
& -viscAh*(vVel(i,j,k,bi,bj)-vVel(i-1,j,k,bi,bj)) |
693 |
& *cosFacV(J,bi,bj) |
694 |
#else |
695 |
& 0. |
696 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
697 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
698 |
& +viscA4*(v4F(i,j) -v4F(i-1,j) ) |
699 |
#ifdef COSINEMETH_III |
700 |
& *sqcosFacV(J,bi,bj) |
701 |
#else |
702 |
& *cosFacV(J,bi,bj) |
703 |
#endif |
704 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
705 |
& )*_recip_dxV(i,j,bi,bj) |
706 |
ENDDO |
707 |
ENDDO |
708 |
DO j=jMin,jMax |
709 |
DO i=iMin,iMax |
710 |
fZon(i,j) = 0. |
711 |
& _ADM( + uDvdxFac * aF (i,j) ) |
712 |
& + AhDvdxFac * vF (i,j) |
713 |
ENDDO |
714 |
ENDDO |
715 |
C-- Meridional flux (fMer is at north face of "v" cell) |
716 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
717 |
C o Mean flow component of meridional flux |
718 |
DO j=jMin,jMax |
719 |
DO i=iMin,iMax |
720 |
af(i,j) = |
721 |
& (vTrans(i,j)+vTrans(i,j+1)) |
722 |
& *(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj))*0.25 _d 0 |
723 |
ENDDO |
724 |
ENDDO |
725 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
726 |
C Eddy component of meridional flux |
727 |
C o Laplacian term |
728 |
DO j=jMin,jMax |
729 |
DO i=iMin,iMax |
730 |
vF(i,j) = |
731 |
#ifdef OLD_UV_GEOMETRY |
732 |
& ( YA(i,j) + YA(i,j+1) ) *0.5 _d 0 |
733 |
#else |
734 |
& _dxF(i,j,bi,bj)*drF(k)*_hFacC(i,j,k,bi,bj) |
735 |
#endif /* OLD_UV_GEOMETRY */ |
736 |
& *( 0. |
737 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
738 |
#ifdef ISOTROPIC_COS_SCALING |
739 |
& -viscAh*(vVel(i,j+1,k,bi,bj)-vVel(i,j,k,bi,bj)) |
740 |
& *cosFacU(J,bi,bj) |
741 |
#else |
742 |
& -viscAh*(vVel(i,j+1,k,bi,bj)-vVel(i,j,k,bi,bj)) |
743 |
#endif |
744 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
745 |
#ifdef INCLUDE_BH_MOMENTUM_DIFFUSION_CODE |
746 |
#ifdef ISOTROPIC_COS_SCALING |
747 |
& +viscA4*(v4F(i,j+1) -v4F(i,j) ) |
748 |
& *cosFacU(J,bi,bj) |
749 |
#else |
750 |
& +viscA4*(v4F(i,j+1) -v4F(i,j) ) |
751 |
#endif |
752 |
#endif /* INCLUDE_BH_MOMENTUM_DIFFUSION_CODE */ |
753 |
& )*_recip_dyF(i,j,bi,bj) |
754 |
|
755 |
ENDDO |
756 |
ENDDO |
757 |
DO j=jMin,jMax |
758 |
DO i=iMin,iMax |
759 |
fMer(i,j) = 0. |
760 |
& _ADM( + vDvdyFac * af(i,j) ) |
761 |
& + AhDvdyFac * vf(i,j) |
762 |
ENDDO |
763 |
ENDDO |
764 |
C-- Vertical flux (fVer is at upper face of "v" cell) |
765 |
#ifdef INCLUDE_MOMENTUM_ADVECTION_CODE |
766 |
C-- Free surface correction temr |
767 |
IF (k.EQ.1) THEN |
768 |
DO j=jMin,jMax |
769 |
DO i=iMin,iMax |
770 |
fVerV(i,j,kUp)=rVelMaskOverride*( |
771 |
& wVel(i, j ,k,bi,bj)*rA(i, j ,bi,bj) |
772 |
& +wVel(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj) |
773 |
& )*0.5 |
774 |
& *vVel(i,j,k,bi,bj) |
775 |
ENDDO |
776 |
ENDDO |
777 |
ENDIF |
778 |
C o Mean flow component of vertical flux |
779 |
DO j=jMin,jMax |
780 |
DO i=iMin,iMax |
781 |
af(i,j) = |
782 |
& wVelBottomOverride*( |
783 |
& wVel(i, j ,kp1,bi,bj)*rA(i, j ,bi,bj) |
784 |
& +wVel(i,j-1,kp1,bi,bj)*rA(i,j-1,bi,bj) |
785 |
& ) |
786 |
& *( vVel(i,j,kp1,bi,bj)+vVel(i,j,k,bi,bj) ) |
787 |
& *0.25 _d 0 |
788 |
#ifdef OLD_ADV_BCS |
789 |
& *_maskS(i,j,kp1,bi,bj)*_maskS(i,j,k,bi,bj) |
790 |
#endif /* OLD_ADV_BCS */ |
791 |
C Note!!!! The line "*(maskS(i,j,k,bi,bj)*maskS(i,j,kM1,bi,bj))" is |
792 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
793 |
C Note!!!! FV paper used a different bc in which this line is |
794 |
C Note!!!! omitted. |
795 |
C Note!!!! Is rVelMaskOverridde right for implicit free-surface? |
796 |
ENDDO |
797 |
ENDDO |
798 |
#endif /* INCLUDE_MOMENTUM_ADVECTION_CODE */ |
799 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
800 |
#ifdef OLD_UV_GEOMETRY |
801 |
#define _RA_AT_S (RA(i,j,bi,bj)*maskC(i,j)+RA(i,j-1,bi,bj)*maskC(i,j-1))*.5 |
802 |
#else |
803 |
#define _RA_AT_S rAs(i,j,bi,bj) |
804 |
#endif /* OLD_UV_GEOMETRY */ |
805 |
C Eddy component of vertical flux |
806 |
IF (.NOT.implicitViscosity) THEN |
807 |
DO j=jMin,jMax |
808 |
DO i=iMin,iMax |
809 |
vf(i,j) = |
810 |
& -KappaRV(i,j,kp1) |
811 |
& *_RA_AT_S |
812 |
& *(vVel(i,j,k,bi,bj)-vVel(i,j,kp1,bi,bj)) |
813 |
& *rkFac*recip_drC(kp1) |
814 |
& *_maskS(i,j,kp1,bi,bj) |
815 |
ENDDO |
816 |
ENDDO |
817 |
ENDIF |
818 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
819 |
IF (implicitViscosity) THEN |
820 |
DO j=jMin,jMax |
821 |
DO i=iMin,iMax |
822 |
fVerV(i,j,kDown) = rVelDvdrFac * af(i,j) |
823 |
ENDDO |
824 |
ENDDO |
825 |
ELSE |
826 |
DO j=jMin,jMax |
827 |
DO i=iMin,iMax |
828 |
fVerV(i,j,kDown) = rVelDvdrFac * af(i,j) |
829 |
& + ArDvdrFac * vf(i,j) |
830 |
ENDDO |
831 |
ENDDO |
832 |
ENDIF |
833 |
|
834 |
#ifdef INCLUDE_GRADPH_CODE |
835 |
C--- Hydorstatic term (-1/rho . dp/dy ) |
836 |
DO j=jMin,jMax |
837 |
DO i=iMin,iMax |
838 |
pF(i,j) = -_recip_dyC(i,j,bi,bj) |
839 |
& *(phiHyd(i,j,k)-phiHyd(i,j-1,k)) |
840 |
ENDDO |
841 |
ENDDO |
842 |
#endif /* INCLUDE_GRADPH_CODE */ |
843 |
|
844 |
C-- Tendency is minus divergence of the fluxes + coriolis + pressure |
845 |
C-- term |
846 |
DO j=jMin,jMax |
847 |
DO i=iMin,iMax |
848 |
gV(i,j,k,bi,bj) = |
849 |
#ifdef OLD_UV_GEOM |
850 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k)/ |
851 |
& ( 0.5 _d 0*(_rA(i,j,bi,bj)+_rA(i,j-1,bi,bj)) ) |
852 |
#else |
853 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k) |
854 |
& *recip_rAs(i,j,bi,bj) |
855 |
#endif |
856 |
& *(fZon(i+1,j) - fZon(i,j ) |
857 |
& +fMer(i,j ) - fMer(i,j-1) |
858 |
& +fVerV(i,j,kUp)*rkFac - fVerV(i,j,kDown)*rkFac |
859 |
& ) |
860 |
& _PHM( +phyFac*pf(i,j) ) |
861 |
ENDDO |
862 |
ENDDO |
863 |
|
864 |
#ifdef INCLUDE_LP_MOMENTUM_DIFFUSION_CODE |
865 |
C-- No-slip and drag BCs appear as body forces in cell abutting topography |
866 |
IF (no_slip_sides) THEN |
867 |
C- No-slip BCs impose a drag at walls... |
868 |
DO j=jMin,jMax |
869 |
DO i=iMin,iMax |
870 |
hFacZClosedW = _hFacS(i,j,k,bi,bj) - hFacZ(i,j) |
871 |
hFacZClosedE = _hFacS(i,j,k,bi,bj) - hFacZ(i+1,j) |
872 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) |
873 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k) |
874 |
& /rAs(i,j,bi,bj) |
875 |
& *( hFacZClosedW*_dyU( i ,j,bi,bj) |
876 |
& *_recip_dxV( i ,j,bi,bj) |
877 |
& +hFacZClosedE*_dyU(i+1,j,bi,bj) |
878 |
& *_recip_dxV(i+1,j,bi,bj) ) |
879 |
& *rkFac*drF(k)*2.*( viscAh*vVel(i,j,k,bi,bj)*cosFacV(J,bi,bj) |
880 |
& -viscA4*v4F(i,j) |
881 |
#ifdef COSINEMETH_III |
882 |
& *sqcosFacV(J,bi,bj) |
883 |
#else |
884 |
& *cosFacV(J,bi,bj) |
885 |
#endif |
886 |
& ) |
887 |
ENDDO |
888 |
ENDDO |
889 |
ENDIF |
890 |
C- No-slip BCs impose a drag at bottom |
891 |
IF (bottomDragTerms) THEN |
892 |
rdrckp1=recip_drC(kp1) |
893 |
IF (k.EQ.Nr) rdrckp1=recip_drF(k) |
894 |
DO j=jMin,jMax |
895 |
DO i=iMin,iMax |
896 |
maskDown=_maskS(i,j,kp1,bi,bj) |
897 |
IF (k.EQ.Nr) maskDown=0. |
898 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) |
899 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k) |
900 |
& *( |
901 |
& 2.*KappaRV(i,j,kp1)*rkFac*rdrckp1 |
902 |
& + bottomDragLinear |
903 |
& )*(1.-maskDown)*vVel(i,j,k,bi,bj) |
904 |
IF (KE(i,j)+KE(i,j-1) .NE. 0.) |
905 |
& gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) |
906 |
& -_recip_hFacS(i,j,k,bi,bj)*recip_drF(k) |
907 |
& *bottomDragQuadratic*sqrt(KE(i,j)+KE(i,j-1)) |
908 |
& *(1.-maskDown)*vVel(i,j,k,bi,bj) |
909 |
ENDDO |
910 |
ENDDO |
911 |
ENDIF |
912 |
#endif /* INCLUDE_LP_MOMENTUM_DIFFUSION_CODE */ |
913 |
|
914 |
#ifdef INCLUDE_MOMENTUM_FORCING_CODE |
915 |
C-- Forcing term |
916 |
CALL EXTERNAL_FORCING_V( |
917 |
I iMin,iMax,jMin,jMax,bi,bj,k, |
918 |
I myCurrentTime,myThid) |
919 |
#endif |
920 |
|
921 |
#ifdef INCLUDE_MOMENTUM_METRIC_TERM_CODE |
922 |
C-- Metric terms for curvilinear grid systems |
923 |
IF ( usingSphericalPolarMTerms ) THEN |
924 |
C o Spherical polar grid metric terms |
925 |
DO j=jMin,jMax |
926 |
DO i=iMin,iMax |
927 |
mT(i,j) = -vVel(i,j,k,bi,bj)*recip_RSphere |
928 |
& *0.25 _d 0*(wVelBottomOverride* |
929 |
& (wVel(i,j,kp1,bi,bj)+wVel(i,j-1,kp1,bi,bj)) |
930 |
& + wVel(i,j, k ,bi,bj)+wVel(i,j-1, k ,bi,bj) |
931 |
& )*rkFac*recip_horiVertRatio |
932 |
ENDDO |
933 |
ENDDO |
934 |
DO j=jMin,jMax |
935 |
DO i=iMin,iMax |
936 |
mT(i,j) = mT(i,j)-recip_RSphere |
937 |
& *0.25 _d 0*( uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj) |
938 |
& +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj) |
939 |
& ) |
940 |
& *0.25 _d 0*( uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj) |
941 |
& +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj) |
942 |
& ) |
943 |
& *_tanPhiAtV(i,j,bi,bj) |
944 |
ENDDO |
945 |
ENDDO |
946 |
DO j=jMin,jMax |
947 |
DO i=iMin,iMax |
948 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+ |
949 |
& mTFacV*mT(i,j) |
950 |
ENDDO |
951 |
ENDDO |
952 |
ENDIF |
953 |
#endif |
954 |
|
955 |
C-- Set dv/dt on boundaries to zero |
956 |
DO j=jMin,jMax |
957 |
DO i=iMin,iMax |
958 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj) |
959 |
ENDDO |
960 |
ENDDO |
961 |
|
962 |
C-- Coriolis term |
963 |
C Note. As coded here, coriolis will not work with "thin walls" |
964 |
#ifdef INCLUDE_CD_CODE |
965 |
C Pressure extrapolated forward in time |
966 |
DO j=jMin,jMax |
967 |
DO i=iMin,iMax |
968 |
pf(i,j) = |
969 |
& ab15*( etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) ) |
970 |
& +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) ) |
971 |
ENDDO |
972 |
ENDDO |
973 |
IF (staggerTimeStep) THEN |
974 |
DO j=jMin,jMax |
975 |
DO i=iMin,iMax |
976 |
pf(i,j) = pf(i,j)+phiHyd(i,j,k) |
977 |
ENDDO |
978 |
ENDDO |
979 |
ENDIF |
980 |
#endif /* INCLUDE_CD_CODE */ |
981 |
C-- Zonal velocity coriolis term |
982 |
C Note. As coded here, coriolis will not work with "thin walls" |
983 |
#ifdef INCLUDE_CD_CODE |
984 |
C-- Coriolis with CD scheme allowed |
985 |
C grady(p) + gV |
986 |
DO j=jMin,jMax |
987 |
DO i=iMin,iMax |
988 |
af(i,j) = -_maskS(i,j,k,bi,bj) |
989 |
& *_recip_dyC(i,j,bi,bj) |
990 |
& *(pf(i,j)-pf(i,j-1)) |
991 |
& +gV(i,j,k,bi,bj) |
992 |
ENDDO |
993 |
ENDDO |
994 |
C Average to Vd point and add coriolis |
995 |
DO j=jMin,jMax |
996 |
DO i=iMin,iMax |
997 |
vf(i,j) = |
998 |
& 0.25 _d 0*( af(i ,j)+af(i ,j+1) |
999 |
& +af(i-1,j)+af(i-1,j+1) |
1000 |
& )*_maskW(i,j,k,bi,bj) |
1001 |
& -0.5 _d 0*(_fCori(i,j,bi,bj)+ |
1002 |
& _fCori(i-1,j,bi,bj)) |
1003 |
& *uVel(i,j,k,bi,bj) |
1004 |
ENDDO |
1005 |
ENDDO |
1006 |
C Step forward Vd |
1007 |
DO j=jMin,jMax |
1008 |
DO i=iMin,iMax |
1009 |
vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) + |
1010 |
& deltaTmom*vf(i,j) |
1011 |
ENDDO |
1012 |
ENDDO |
1013 |
C Relax D grid V to C grid V |
1014 |
DO j=jMin,jMax |
1015 |
DO i=iMin,iMax |
1016 |
vVelD(i,j,k,bi,bj) = rCD*vVelD(i,j,k,bi,bj) |
1017 |
& +(1. _d 0 - rCD)*( |
1018 |
& ab15*0.25 _d 0*( |
1019 |
& vVel(i ,j ,k,bi,bj)+vVel(i ,j+1,k,bi,bj) |
1020 |
& +vVel(i-1,j ,k,bi,bj)+vVel(i-1,j+1,k,bi,bj) |
1021 |
& )*_maskW(i,j,k,bi,bj) |
1022 |
& + |
1023 |
& ab05*0.25 _d 0*( |
1024 |
& vNM1(i ,j ,k,bi,bj)+vNM1(i ,j+1,k,bi,bj) |
1025 |
& +vNM1(i-1,j ,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj) |
1026 |
& )*_maskW(i,j,k,bi,bj) |
1027 |
& ) |
1028 |
ENDDO |
1029 |
ENDDO |
1030 |
C Calculate coriolis force on U |
1031 |
DO j=jMin,jMax |
1032 |
DO i=iMin,iMax |
1033 |
guCD(i,j,k,bi,bj) = |
1034 |
& 0.5 _d 0 *( _fCori(i ,j,bi,bj) + |
1035 |
& _fCori(i-1,j,bi,bj) ) |
1036 |
& *vVelD(i,j,k,bi,bj)*fuFac |
1037 |
ENDDO |
1038 |
ENDDO |
1039 |
#endif /* INCLUDE_CD_CODE */ |
1040 |
#ifndef INCLUDE_CD_CODE |
1041 |
C-- No CD scheme support |
1042 |
DO j=jMin,jMax |
1043 |
DO i=iMin,iMax |
1044 |
cf(i,j) = |
1045 |
& 0.5 _d 0 *( _fCori(i ,j,bi,bj) + |
1046 |
& _fCori(i-1,j,bi,bj) ) |
1047 |
& *0.25 _d 0 *( |
1048 |
& vVel(i ,j,k,bi,bj)+vVel(i ,j+1,k,bi,bj) |
1049 |
& +vVel(i-1,j,k,bi,bj)+vVel(i-1,j+1,k,bi,bj) |
1050 |
& ) |
1051 |
ENDDO |
1052 |
ENDDO |
1053 |
DO j=jMin,jMax |
1054 |
DO i=iMin,iMax |
1055 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) |
1056 |
& +fuFac*cf(i,j) |
1057 |
ENDDO |
1058 |
ENDDO |
1059 |
#endif /* !INCLUDE_CD_CODE */ |
1060 |
|
1061 |
|
1062 |
C-- Meridional velocity coriolis term |
1063 |
#ifdef INCLUDE_CD_CODE |
1064 |
C gradx(p)+gU |
1065 |
DO j=jMin,jMax |
1066 |
DO i=iMin,iMax |
1067 |
af(i,j) = -_maskW(i,j,k,bi,bj) |
1068 |
& *_recip_dxC(i,j,bi,bj)* |
1069 |
& (pf(i,j)-pf(i-1,j)) |
1070 |
& +gU(i,j,k,bi,bj) |
1071 |
ENDDO |
1072 |
ENDDO |
1073 |
C Average to Ud point and add coriolis |
1074 |
DO j=jMin,jMax |
1075 |
DO i=iMin,iMax |
1076 |
vf(i,j) = |
1077 |
& 0.25 _d 0*( af(i ,j)+af(i ,j-1) |
1078 |
& +af(i+1,j)+af(i+1,j-1) |
1079 |
& )*_maskS(i,j,k,bi,bj) |
1080 |
& +0.5 _d 0*( _fCori(i,j,bi,bj) |
1081 |
& +_fCori(i,j-1,bi,bj)) |
1082 |
& *vVel(i,j,k,bi,bj) |
1083 |
ENDDO |
1084 |
ENDDO |
1085 |
C Step forward Ud |
1086 |
DO j=jMin,jMax |
1087 |
DO i=iMin,iMax |
1088 |
uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) + |
1089 |
& deltaTmom*vf(i,j)*_maskS(i,j,k,bi,bj) |
1090 |
ENDDO |
1091 |
ENDDO |
1092 |
C Relax D grid U to C grid U |
1093 |
DO j=jMin,jMax |
1094 |
DO i=iMin,iMax |
1095 |
uVelD(i,j,k,bi,bj) = rCD*uVelD(i,j,k,bi,bj) |
1096 |
& +(1. _d 0 - rCD)*( |
1097 |
& ab15*0.25 _d 0*( |
1098 |
& uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj) |
1099 |
& +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj) |
1100 |
& )*_maskS(i,j,k,bi,bj) |
1101 |
& + |
1102 |
& ab05*0.25 _d 0*( |
1103 |
& uNM1(i,j ,k,bi,bj)+uNM1(i+1,j ,k,bi,bj) |
1104 |
& +uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj) |
1105 |
& )*_maskS(i,j,k,bi,bj) |
1106 |
& ) |
1107 |
ENDDO |
1108 |
ENDDO |
1109 |
C Calculate coriolis force on V |
1110 |
DO j=jMin,jMax |
1111 |
DO i=iMin,iMax |
1112 |
gvCD(i,j,k,bi,bj) = |
1113 |
& -0.5 _d 0 *( _fCori(i ,j,bi,bj) |
1114 |
& +_fCori(i,j-1,bi,bj) ) |
1115 |
& *uVelD(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj)*fvFac |
1116 |
ENDDO |
1117 |
ENDDO |
1118 |
#endif /* INCLUDE_CD_CODE */ |
1119 |
#ifndef INCLUDE_CD_CODE |
1120 |
C-- No CD scheme support |
1121 |
DO j=jMin,jMax |
1122 |
DO i=iMin,iMax |
1123 |
cf(i,j) = |
1124 |
& -0.5 _d 0 *(_fCori(i,j ,bi,bj)+_fCori(i,j-1,bi,bj)) |
1125 |
& *0.25 _d 0 |
1126 |
& *( uVel(i,j ,k,bi,bj)+uVel(i+1,j ,k,bi,bj) |
1127 |
& +uVel(i,j-1,k,bi,bj)+uVel(i+1,j-1,k,bi,bj) |
1128 |
& ) |
1129 |
ENDDO |
1130 |
ENDDO |
1131 |
DO j=jMin,jMax |
1132 |
DO i=iMin,iMax |
1133 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) |
1134 |
& +fvFac*cf(i,j) |
1135 |
ENDDO |
1136 |
ENDDO |
1137 |
#endif /* !INCLUDE_CD_CODE */ |
1138 |
|
1139 |
C-- Set du/dt on boundaries to zero |
1140 |
DO j=jMin,jMax |
1141 |
DO i=iMin,iMax |
1142 |
gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)*_maskW(i,j,k,bi,bj) |
1143 |
ENDDO |
1144 |
ENDDO |
1145 |
C-- Set dv/dt on boundaries to zero |
1146 |
DO j=jMin,jMax |
1147 |
DO i=iMin,iMax |
1148 |
gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)*_maskS(i,j,k,bi,bj) |
1149 |
ENDDO |
1150 |
ENDDO |
1151 |
|
1152 |
#ifdef INCLUDE_CD_CODE |
1153 |
C-- Save "previous time level" variables |
1154 |
DO j=1-OLy,sNy+OLy |
1155 |
DO i=1-OLx,sNx+OLx |
1156 |
uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj) |
1157 |
vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj) |
1158 |
ENDDO |
1159 |
ENDDO |
1160 |
#endif /* INCLUDE_CD_CODE */ |
1161 |
|
1162 |
RETURN |
1163 |
END |