/[MITgcm]/MITgcm/model/src/calc_gw.F
ViewVC logotype

Contents of /MITgcm/model/src/calc_gw.F

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


Revision 1.31 - (show annotations) (download)
Wed Jul 12 12:27:38 2006 UTC (17 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.30: +38 -38 lines
true flux-form, account for horizontal grid spacing.

1 C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.28 2006/07/07 20:09:37 jmc Exp $
2 C !DESCRIPTION: \bv
3 C $Name: $
4
5 #include "CPP_OPTIONS.h"
6
7 CBOP
8 C !ROUTINE: CALC_GW
9 C !INTERFACE:
10 SUBROUTINE CALC_GW(
11 I bi, bj, KappaRU, KappaRV,
12 I myTime, myIter, myThid )
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | S/R CALC_GW
16 C | o Calculate vertical velocity tendency terms
17 C | ( Non-Hydrostatic only )
18 C *==========================================================*
19 C | In NH, the vertical momentum tendency must be
20 C | calculated explicitly and included as a source term
21 C | for a 3d pressure eqn. Calculate that term here.
22 C | This routine is not used in HYD calculations.
23 C *==========================================================*
24 C \ev
25
26 C !USES:
27 IMPLICIT NONE
28 C == Global variables ==
29 #include "SIZE.h"
30 #include "EEPARAMS.h"
31 #include "PARAMS.h"
32 #include "GRID.h"
33 #include "DYNVARS.h"
34 #include "NH_VARS.h"
35
36 C !INPUT/OUTPUT PARAMETERS:
37 C == Routine arguments ==
38 C bi,bj :: current tile indices
39 C KappaRU :: vertical viscosity at U points
40 C KappaRV :: vertical viscosity at V points
41 C myTime :: Current time in simulation
42 C myIter :: Current iteration number in simulation
43 C myThid :: Thread number for this instance of the routine.
44 INTEGER bi,bj
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 myTime
48 INTEGER myIter
49 INTEGER myThid
50
51 #ifdef ALLOW_NONHYDROSTATIC
52
53 C !LOCAL VARIABLES:
54 C == Local variables ==
55 C iMin,iMax
56 C jMin,jMax
57 C xA :: W-Cell face area normal to X
58 C yA :: W-Cell face area normal to Y
59 C rThickC_W :: thickness (in r-units) of W-Cell at Western Edge
60 C rThickC_S :: thickness (in r-units) of W-Cell at Southern Edge
61 C recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
62 C flx_NS :: vertical momentum flux, meridional direction
63 C flx_EW :: vertical momentum flux, zonal direction
64 C flxAdvUp :: vertical mom. advective flux, vertical direction (@ level k-1)
65 C flxDisUp :: vertical mom. dissipation flux, vertical direction (@ level k-1)
66 C flx_Dn :: vertical momentum flux, vertical direction (@ level k)
67 C gwDiss :: vertical momentum dissipation tendency
68 C i,j,k :: Loop counters
69 INTEGER iMin,iMax,jMin,jMax
70 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL rThickC_W (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL rThickC_S (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 _RL flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80 _RL gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 _RL gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82 _RL del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 INTEGER i,j,k, kp1
84 _RL wOverride
85 _RL tmp_WbarZ
86 _RL uTrans, vTrans, rTrans
87 _RL viscLoc
88 _RL halfRL
89 _RS halfRS, zeroRS
90 PARAMETER( halfRL = 0.5D0 )
91 PARAMETER( halfRS = 0.5 , zeroRS = 0. )
92 CEOP
93
94 C Catch barotropic mode
95 IF ( Nr .LT. 2 ) RETURN
96
97 iMin = 1
98 iMax = sNx
99 jMin = 1
100 jMax = sNy
101
102 C-- Initialise gW to zero
103 DO k=1,Nr
104 DO j=1-OLy,sNy+OLy
105 DO i=1-OLx,sNx+OLx
106 gW(i,j,k,bi,bj) = 0.
107 ENDDO
108 ENDDO
109 ENDDO
110 C- Initialise gwDiss to zero
111 DO j=1-OLy,sNy+OLy
112 DO i=1-OLx,sNx+OLx
113 gwDiss(i,j) = 0.
114 ENDDO
115 ENDDO
116
117 C-- Boundaries condition at top
118 DO j=1-OLy,sNy+OLy
119 DO i=1-OLx,sNx+OLx
120 flxAdvUp(i,j) = 0.
121 flxDisUp(i,j) = 0.
122 ENDDO
123 ENDDO
124
125 C--- Sweep down column
126 DO k=2,Nr
127 kp1=k+1
128 wOverRide=1.
129 IF (k.EQ.Nr) THEN
130 kp1=Nr
131 wOverRide=0.
132 ENDIF
133 C-- Compute grid factor arround a W-point:
134 DO j=1-Oly,sNy+Oly
135 DO i=1-Olx,sNx+Olx
136 C- note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
137 rThickC_W(i,j) =
138 & drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
139 & + drF( k )*MIN( _hFacW(i,j,k ,bi,bj), halfRS )
140 rThickC_S(i,j) =
141 & drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
142 & + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
143 IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
144 recip_rThickC(i,j) = 0.
145 ELSE
146 recip_rThickC(i,j) = 1. _d 0 /
147 & ( drF(k-1)*halfRS +
148 & + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
149 & )
150 ENDIF
151 C W-Cell Western face area:
152 xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
153 C W-Cell Southern face area:
154 yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
155 ENDDO
156 ENDDO
157
158 C-- horizontal bi-harmonic dissipation
159 IF (momViscosity .AND. viscA4W.NE.0. ) THEN
160 C- calculate the horizontal Laplacian of vertical flow
161 C Zonal flux d/dx W
162 DO j=1-Oly,sNy+Oly
163 flx_EW(1-Olx,j)=0.
164 DO i=1-Olx+1,sNx+Olx
165 flx_EW(i,j) =
166 & (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
167 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
168 #ifdef COSINEMETH_III
169 & *sqcosFacU(j,bi,bj)
170 #endif
171 ENDDO
172 ENDDO
173 C Meridional flux d/dy W
174 DO i=1-Olx,sNx+Olx
175 flx_NS(i,1-Oly)=0.
176 ENDDO
177 DO j=1-Oly+1,sNy+Oly
178 DO i=1-Olx,sNx+Olx
179 flx_NS(i,j) =
180 & (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
181 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
182 #ifdef ISOTROPIC_COS_SCALING
183 #ifdef COSINEMETH_III
184 & *sqCosFacV(j,bi,bj)
185 #endif
186 #endif
187 ENDDO
188 ENDDO
189
190 C del^2 W
191 C Difference of zonal fluxes ...
192 DO j=1-Oly,sNy+Oly
193 DO i=1-Olx,sNx+Olx-1
194 del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
195 ENDDO
196 del2w(sNx+Olx,j)=0.
197 ENDDO
198
199 C ... add difference of meridional fluxes and divide by volume
200 DO j=1-Oly,sNy+Oly-1
201 DO i=1-Olx,sNx+Olx
202 del2w(i,j) = ( del2w(i,j)
203 & +(flx_NS(i,j+1)-flx_NS(i,j))
204 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
205 ENDDO
206 ENDDO
207 C-- No-slip BCs impose a drag at walls...
208 CML ************************************************************
209 CML No-slip Boundary conditions for bi-harmonic dissipation
210 CML need to be implemented here!
211 CML ************************************************************
212 ELSE
213 C- Initialize del2w to zero:
214 DO j=1-Oly,sNy+Oly
215 DO i=1-Olx,sNx+Olx
216 del2w(i,j) = 0. _d 0
217 ENDDO
218 ENDDO
219 ENDIF
220
221 IF (momViscosity) THEN
222 C Viscous Flux on Western face
223 DO j=jMin,jMax
224 DO i=iMin,iMax+1
225 flx_EW(i,j)=
226 & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
227 & *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
228 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
229 cOld & *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
230 & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
231 & *(del2w(i,j)-del2w(i-1,j))
232 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
233 cOld & *_recip_dxC(i,j,bi,bj)*drC(k)
234 #ifdef COSINEMETH_III
235 & *sqCosFacU(j,bi,bj)
236 #else
237 & *CosFacU(j,bi,bj)
238 #endif
239 ENDDO
240 ENDDO
241 C Viscous Flux on Southern face
242 DO j=jMin,jMax+1
243 DO i=iMin,iMax
244 flx_NS(i,j)=
245 & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
246 & *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
247 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
248 cOld & *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
249 & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
250 & *(del2w(i,j)-del2w(i,j-1))
251 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
252 cOld & *_recip_dyC(i,j,bi,bj)*drC(k)
253 #ifdef ISOTROPIC_COS_SCALING
254 #ifdef COSINEMETH_III
255 & *sqCosFacV(j,bi,bj)
256 #else
257 & *CosFacV(j,bi,bj)
258 #endif
259 #endif
260 ENDDO
261 ENDDO
262 C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
263 DO j=jMin,jMax
264 DO i=iMin,iMax
265 C Interpolate vert viscosity to center of tracer-cell (level k):
266 viscLoc = ( KappaRU(i,j,k) +KappaRU(i+1,j,k)
267 & +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
268 & +KappaRV(i,j,k) +KappaRV(i,j+1,k)
269 & +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
270 & )*0.125 _d 0
271 flx_Dn(i,j) =
272 & - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
273 & -wVel(i,j, k ,bi,bj) )*rkSign
274 & *recip_drF(k)*rA(i,j,bi,bj)
275 cOld & *recip_drF(k)
276 ENDDO
277 ENDDO
278 C Tendency is minus divergence of viscous fluxes:
279 DO j=jMin,jMax
280 DO i=iMin,iMax
281 gwDiss(i,j) =
282 & -( ( flx_EW(i+1,j)-flx_EW(i,j) )
283 & + ( flx_NS(i,j+1)-flx_NS(i,j) )
284 & + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
285 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
286 cOld gwDiss(i,j) =
287 cOld & -(
288 cOld & +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
289 cOld & +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
290 cOld & + ( flxDisUp(i,j)-flx_Dn(i,j) )
291 c & )*recip_rThickC(i,j)
292 cOld & )*recip_drC(k)
293 C-- prepare for next level (k+1)
294 flxDisUp(i,j)=flx_Dn(i,j)
295 ENDDO
296 ENDDO
297 ENDIF
298
299 IF (no_slip_sides) THEN
300 C- No-slip BCs impose a drag at walls...
301 c CALL MOM_W_SIDEDRAG(
302 c I bi,bj,k,
303 c O gwAdd,
304 c I myThid)
305 c DO j=jMin,jMax
306 c DO i=iMin,iMax
307 c gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
308 c ENDDO
309 c ENDDO
310 ENDIF
311
312 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313
314 IF ( momAdvection ) THEN
315 C Advective Flux on Western face
316 DO j=jMin,jMax
317 DO i=iMin,iMax+1
318 C transport through Western face area:
319 uTrans = (
320 & drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)
321 & + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
322 & )*halfRL*_dyG(i,j,bi,bj)
323 cOld & )*halfRL
324 flx_EW(i,j)=
325 & uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
326 ENDDO
327 ENDDO
328 C Advective Flux on Southern face
329 DO j=jMin,jMax+1
330 DO i=iMin,iMax
331 C transport through Southern face area:
332 vTrans = (
333 & drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)
334 & +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
335 & )*halfRL*_dxG(i,j,bi,bj)
336 cOld & )*halfRL
337 flx_NS(i,j)=
338 & vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
339 ENDDO
340 ENDDO
341 C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
342 DO j=jMin,jMax
343 DO i=iMin,iMax
344 tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)
345 & +wVel(i,j,kp1,bi,bj)*wOverRide )
346 C transport through Lower face area:
347 rTrans = tmp_WbarZ*rA(i,j,bi,bj)
348 flx_Dn(i,j) = rTrans*tmp_WbarZ
349 cOld flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ
350 ENDDO
351 ENDDO
352 C Tendency is minus divergence of advective fluxes:
353 DO j=jMin,jMax
354 DO i=iMin,iMax
355 gW(i,j,k,bi,bj) =
356 & -( ( flx_EW(i+1,j)-flx_EW(i,j) )
357 & + ( flx_NS(i,j+1)-flx_NS(i,j) )
358 & + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign
359 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
360 cOld gW(i,j,k,bi,bj) =
361 cOld & -(
362 cOld & +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
363 cOld & +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
364 cOld & + ( flxAdvUp(i,j)-flx_Dn(i,j) )
365 c & )*recip_rThickC(i,j)
366 cOld & )*recip_drC(k)
367 C-- prepare for next level (k+1)
368 flxAdvUp(i,j)=flx_Dn(i,j)
369 ENDDO
370 ENDDO
371 ENDIF
372
373 IF ( useNHMTerms ) THEN
374 CALL MOM_W_METRIC_NH(
375 I bi,bj,k,
376 I uVel, vVel,
377 O gwAdd,
378 I myThid )
379 DO j=jMin,jMax
380 DO i=iMin,iMax
381 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
382 ENDDO
383 ENDDO
384 ENDIF
385 IF ( useCoriolis ) THEN
386 CALL MOM_W_CORIOLIS_NH(
387 I bi,bj,k,
388 I uVel, vVel,
389 O gwAdd,
390 I myThid )
391 DO j=jMin,jMax
392 DO i=iMin,iMax
393 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
394 ENDDO
395 ENDDO
396 ENDIF
397
398 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
399
400 C-- Dissipation term inside the Adams-Bashforth:
401 IF ( momViscosity .AND. momDissip_In_AB) THEN
402 DO j=jMin,jMax
403 DO i=iMin,iMax
404 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
405 ENDDO
406 ENDDO
407 ENDIF
408
409 C- Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
410 C and save gW_[n] into gwNm1 for the next time step.
411 c#ifdef ALLOW_ADAMSBASHFORTH_3
412 c CALL ADAMS_BASHFORTH3(
413 c I bi, bj, k,
414 c U gW, gwNm,
415 c I momStartAB, myIter, myThid )
416 c#else /* ALLOW_ADAMSBASHFORTH_3 */
417 CALL ADAMS_BASHFORTH2(
418 I bi, bj, k,
419 U gW, gwNm1,
420 I myIter, myThid )
421 c#endif /* ALLOW_ADAMSBASHFORTH_3 */
422
423 C-- Dissipation term outside the Adams-Bashforth:
424 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
425 DO j=jMin,jMax
426 DO i=iMin,iMax
427 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
428 ENDDO
429 ENDDO
430 ENDIF
431
432 C- end of the k loop
433 ENDDO
434
435 #endif /* ALLOW_NONHYDROSTATIC */
436
437 RETURN
438 END

  ViewVC Help
Powered by ViewVC 1.1.22