/[MITgcm]/MITgcm/verification/tutorial_deep_convection/code/calc_gw.F
ViewVC logotype

Contents of /MITgcm/verification/tutorial_deep_convection/code/calc_gw.F

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


Revision 1.1 - (show annotations) (download)
Mon Sep 30 18:19:41 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64p
add 1rst version of isotropic 3-D Smagorinsky code (from Louis-Philippe),
  for now all in this exp. code dir.

1 C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.52 2013/07/28 21:07:03 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_MOM_COMMON
7 # include "MOM_COMMON_OPTIONS.h"
8 #endif
9 #define CALC_GW_NEW_THICK
10
11 CBOP
12 C !ROUTINE: CALC_GW
13 C !INTERFACE:
14 SUBROUTINE CALC_GW(
15 I bi, bj, KappaRU, KappaRV,
16 I myTime, myIter, myThid )
17 C !DESCRIPTION: \bv
18 C *==========================================================*
19 C | S/R CALC_GW
20 C | o Calculate vertical velocity tendency terms
21 C | ( Non-Hydrostatic only )
22 C *==========================================================*
23 C | In NH, the vertical momentum tendency must be
24 C | calculated explicitly and included as a source term
25 C | for a 3d pressure eqn. Calculate that term here.
26 C | This routine is not used in HYD calculations.
27 C *==========================================================*
28 C \ev
29
30 C !USES:
31 IMPLICIT NONE
32 C == Global variables ==
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "GRID.h"
37 #include "RESTART.h"
38 #include "SURFACE.h"
39 #include "DYNVARS.h"
40 #include "NH_VARS.h"
41 #ifdef ALLOW_ADDFLUID
42 # include "FFIELDS.h"
43 #endif
44 #ifdef ALLOW_MOM_COMMON
45 # include "MOM_VISC.h"
46 #endif
47
48 C !INPUT/OUTPUT PARAMETERS:
49 C == Routine arguments ==
50 C bi,bj :: current tile indices
51 C KappaRU :: vertical viscosity at U points
52 C KappaRV :: vertical viscosity at V points
53 C myTime :: Current time in simulation
54 C myIter :: Current iteration number in simulation
55 C myThid :: Thread number for this instance of the routine.
56 INTEGER bi,bj
57 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
58 _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59 _RL myTime
60 INTEGER myIter
61 INTEGER myThid
62
63 #ifdef ALLOW_NONHYDROSTATIC
64 #ifdef ALLOW_MOM_COMMON
65
66 C !LOCAL VARIABLES:
67 C == Local variables ==
68 C biharmonicVisc:: use horizontal biharmonic viscosity for vertical momentum
69 C iMin,iMax
70 C jMin,jMax
71 C xA :: W-Cell face area normal to X
72 C yA :: W-Cell face area normal to Y
73 C rThickC_W :: thickness (in r-units) of W-Cell at Western Edge
74 C rThickC_S :: thickness (in r-units) of W-Cell at Southern Edge
75 C rThickC_C :: thickness (in r-units) of W-Cell (centered on W pt)
76 C recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
77 C flx_NS :: vertical momentum flux, meridional direction
78 C flx_EW :: vertical momentum flux, zonal direction
79 C flxAdvUp :: vertical mom. advective flux, vertical direction (@ level k-1)
80 C flxDisUp :: vertical mom. dissipation flux, vertical direction (@ level k-1)
81 C flx_Dn :: vertical momentum flux, vertical direction (@ level k)
82 C gwDiss :: vertical momentum dissipation tendency
83 C gwAdd :: other tendencies (Coriolis, Metric-terms)
84 C gw_AB :: tendency increment from Adams-Bashforth
85 C del2w :: laplacian of wVel
86 C wFld :: local copy of wVel
87 C i,j,k :: Loop counters
88 LOGICAL biharmonicVisc
89 INTEGER iMin,iMax,jMin,jMax
90 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL rThickC_W (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL rThickC_S (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL rThickC_C (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97 _RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101 _RL gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL gw_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106 INTEGER i,j,k, km1, kp1
107 _RL mskM1, mskP1
108 _RL tmp_WbarZ
109 _RL uTrans, vTrans, rTrans
110 _RL viscLoc
111 PARAMETER( iMin = 1 , iMax = sNx )
112 PARAMETER( jMin = 1 , jMax = sNy )
113 CEOP
114 #ifdef ALLOW_DIAGNOSTICS
115 LOGICAL diagDiss, diagAdvec, diag_AB
116 LOGICAL DIAGNOSTICS_IS_ON
117 EXTERNAL DIAGNOSTICS_IS_ON
118 #endif /* ALLOW_DIAGNOSTICS */
119
120 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
121
122 #ifdef ALLOW_DIAGNOSTICS
123 IF ( useDiagnostics ) THEN
124 diagDiss = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid )
125 diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid )
126 diag_AB = DIAGNOSTICS_IS_ON( 'AB_gW ', myThid )
127 ELSE
128 diagDiss = .FALSE.
129 diagAdvec = .FALSE.
130 diag_AB = .FALSE.
131 ENDIF
132 #endif /* ALLOW_DIAGNOSTICS */
133
134 biharmonicVisc = viscA4W.NE.zeroRL
135 & .OR. ( useVariableVisc .AND. useBiharmonicVisc )
136
137 C-- Initialise gW to zero
138 DO k=1,Nr
139 DO j=1-OLy,sNy+OLy
140 DO i=1-OLx,sNx+OLx
141 gW(i,j,k,bi,bj) = 0.
142 ENDDO
143 ENDDO
144 ENDDO
145 C- Initialise gwDiss to zero
146 DO j=1-OLy,sNy+OLy
147 DO i=1-OLx,sNx+OLx
148 gwDiss(i,j) = 0.
149 ENDDO
150 ENDDO
151 IF (momViscosity) THEN
152 C- Initialize del2w to zero:
153 DO j=1-OLy,sNy+OLy
154 DO i=1-OLx,sNx+OLx
155 del2w(i,j) = 0. _d 0
156 ENDDO
157 ENDDO
158 ENDIF
159
160 C-- Boundaries condition at top (vertical advection of vertical momentum):
161 DO j=1-OLy,sNy+OLy
162 DO i=1-OLx,sNx+OLx
163 flxAdvUp(i,j) = 0.
164 c flxDisUp(i,j) = 0.
165 ENDDO
166 ENDDO
167
168 C--- Sweep down column
169 DO k=1,Nr
170 km1 = MAX( k-1, 1 )
171 kp1 = MIN( k+1,Nr )
172 mskM1 = 1.
173 mskP1 = 1.
174 IF ( k.EQ. 1 ) mskM1 = 0.
175 IF ( k.EQ.Nr ) mskP1 = 0.
176 IF ( k.GT.1 ) THEN
177 C-- Compute grid factor arround a W-point:
178 #ifdef CALC_GW_NEW_THICK
179 DO j=1-OLy,sNy+OLy
180 DO i=1-OLx,sNx+OLx
181 IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
182 & maskC(i,j, k ,bi,bj).EQ.0. ) THEN
183 recip_rThickC(i,j) = 0.
184 ELSE
185 C- valid in z & p coord.; also accurate if Interface @ middle between 2 centers
186 recip_rThickC(i,j) = 1. _d 0 /
187 & ( MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
188 & - MAX( R_low(i,j,bi,bj), rC(k) )
189 & )
190 ENDIF
191 ENDDO
192 ENDDO
193 IF (momViscosity) THEN
194 DO j=1-OLy,sNy+OLy
195 DO i=1-OLx,sNx+OLx
196 rThickC_C(i,j) = MAX( zeroRS,
197 & MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
198 & -MAX( R_low(i,j,bi,bj), rC(k) )
199 & )
200 ENDDO
201 ENDDO
202 DO j=1-OLy,sNy+OLy
203 DO i=1-OLx+1,sNx+OLx
204 rThickC_W(i,j) = MAX( zeroRS,
205 & MIN( rSurfW(i,j,bi,bj), rC(k-1) )
206 & -MAX( rLowW(i,j,bi,bj), rC(k) )
207 & )
208 C W-Cell Western face area:
209 xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
210 c & *deepFacF(k)
211 ENDDO
212 ENDDO
213 DO j=1-OLy+1,sNy+OLy
214 DO i=1-OLx,sNx+OLx
215 rThickC_S(i,j) = MAX( zeroRS,
216 & MIN( rSurfS(i,j,bi,bj), rC(k-1) )
217 & -MAX( rLowS(i,j,bi,bj), rC(k) )
218 & )
219 C W-Cell Southern face area:
220 yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
221 c & *deepFacF(k)
222 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
223 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
224 ENDDO
225 ENDDO
226 ENDIF
227 #else /* CALC_GW_NEW_THICK */
228 DO j=1-OLy,sNy+OLy
229 DO i=1-OLx,sNx+OLx
230 C- note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
231 IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
232 recip_rThickC(i,j) = 0.
233 ELSE
234 recip_rThickC(i,j) = 1. _d 0 /
235 & ( drF(k-1)*halfRS
236 & + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
237 & )
238 ENDIF
239 c IF (momViscosity) THEN
240 #ifdef NONLIN_FRSURF
241 rThickC_C(i,j) =
242 & drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
243 & + drF( k )*MIN( h0FacC(i,j,k ,bi,bj), halfRS )
244 #else
245 rThickC_C(i,j) =
246 & drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
247 & + drF( k )*MIN( _hFacC(i,j,k ,bi,bj), halfRS )
248 #endif
249 rThickC_W(i,j) =
250 & drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
251 & + drF( k )*MIN( _hFacW(i,j,k ,bi,bj), halfRS )
252 rThickC_S(i,j) =
253 & drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
254 & + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
255 C W-Cell Western face area:
256 xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
257 c & *deepFacF(k)
258 C W-Cell Southern face area:
259 yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
260 c & *deepFacF(k)
261 C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
262 C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
263 c ENDIF
264 ENDDO
265 ENDDO
266 #endif /* CALC_GW_NEW_THICK */
267 ELSEIF ( selectNHfreeSurf.GE.1 ) THEN
268 DO j=1-OLy,sNy+OLy
269 DO i=1-OLx,sNx+OLx
270 recip_rThickC(i,j) = recip_drC(k)
271 c rThickC_C(i,j) = drC(k)
272 c rThickC_W(i,j) = drC(k)
273 c rThickC_S(i,j) = drC(k)
274 c xA(i,j) = _dyG(i,j,bi,bj)*drC(k)
275 c yA(i,j) = _dxG(i,j,bi,bj)*drC(k)
276 ENDDO
277 ENDDO
278 ENDIF
279
280 C-- local copy of wVel:
281 DO j=1-OLy,sNy+OLy
282 DO i=1-OLx,sNx+OLx
283 wFld(i,j) = wVel(i,j,k,bi,bj)
284 ENDDO
285 ENDDO
286
287 C-- horizontal bi-harmonic dissipation
288 IF ( momViscosity .AND. k.GT.1 .AND. biharmonicVisc ) THEN
289
290 C- calculate the horizontal Laplacian of vertical flow
291 C Zonal flux d/dx W
292 IF ( useCubedSphereExchange ) THEN
293 C to compute d/dx(W), fill corners with appropriate values:
294 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
295 & wFld, bi,bj, myThid )
296 ENDIF
297 DO j=1-OLy,sNy+OLy
298 flx_EW(1-OLx,j)=0.
299 DO i=1-OLx+1,sNx+OLx
300 flx_EW(i,j) =
301 & ( wFld(i,j) - wFld(i-1,j) )
302 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
303 #ifdef COSINEMETH_III
304 & *sqCosFacU(j,bi,bj)
305 #endif
306 #ifdef ALLOW_OBCS
307 & *maskInW(i,j,bi,bj)
308 #endif
309 ENDDO
310 ENDDO
311
312 C Meridional flux d/dy W
313 IF ( useCubedSphereExchange ) THEN
314 C to compute d/dy(W), fill corners with appropriate values:
315 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
316 & wFld, bi,bj, myThid )
317 ENDIF
318 DO i=1-OLx,sNx+OLx
319 flx_NS(i,1-OLy)=0.
320 ENDDO
321 DO j=1-OLy+1,sNy+OLy
322 DO i=1-OLx,sNx+OLx
323 flx_NS(i,j) =
324 & ( wFld(i,j) - wFld(i,j-1) )
325 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
326 #ifdef ISOTROPIC_COS_SCALING
327 #ifdef COSINEMETH_III
328 & *sqCosFacV(j,bi,bj)
329 #endif
330 #endif
331 #ifdef ALLOW_OBCS
332 & *maskInS(i,j,bi,bj)
333 #endif
334 ENDDO
335 ENDDO
336
337 C del^2 W
338 C Divergence of horizontal fluxes
339 DO j=1-OLy,sNy+OLy-1
340 DO i=1-OLx,sNx+OLx-1
341 del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) )
342 & +( flx_NS(i,j+1)-flx_NS(i,j) )
343 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
344 & *recip_deepFac2F(k)
345 ENDDO
346 ENDDO
347 C end if biharmonic viscosity
348 ENDIF
349
350 IF ( momViscosity .AND. k.GT.1 ) THEN
351 C Viscous Flux on Western face
352 DO j=jMin,jMax
353 DO i=iMin,iMax+1
354 flx_EW(i,j)=
355 & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
356 & *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
357 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
358 & *cosFacU(j,bi,bj)
359 & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
360 & *(del2w(i,j)-del2w(i-1,j))
361 & *_recip_dxC(i,j,bi,bj)*xA(i,j)
362 #ifdef COSINEMETH_III
363 & *sqCosFacU(j,bi,bj)
364 #else
365 & *cosFacU(j,bi,bj)
366 #endif
367 ENDDO
368 ENDDO
369 C Viscous Flux on Southern face
370 DO j=jMin,jMax+1
371 DO i=iMin,iMax
372 flx_NS(i,j)=
373 & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
374 & *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
375 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
376 #ifdef ISOTROPIC_COS_SCALING
377 & *cosFacV(j,bi,bj)
378 #endif
379 & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
380 & *(del2w(i,j)-del2w(i,j-1))
381 & *_recip_dyC(i,j,bi,bj)*yA(i,j)
382 #ifdef ISOTROPIC_COS_SCALING
383 #ifdef COSINEMETH_III
384 & *sqCosFacV(j,bi,bj)
385 #else
386 & *cosFacV(j,bi,bj)
387 #endif
388 #endif
389 ENDDO
390 ENDDO
391 C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
392 DO j=jMin,jMax
393 DO i=iMin,iMax
394 C Interpolate vert viscosity to center of tracer-cell (level k):
395 viscLoc = ( KappaRU(i,j,k) +KappaRU(i+1,j,k)
396 & +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
397 & +KappaRV(i,j,k) +KappaRV(i,j+1,k)
398 & +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
399 & )*0.125 _d 0
400 flx_Dn(i,j) =
401 & - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1
402 & -wVel(i,j, k ,bi,bj) )*rkSign
403 & *recip_drF(k)*rA(i,j,bi,bj)
404 & *deepFac2C(k)*rhoFacC(k)
405 ENDDO
406 ENDDO
407 IF ( k.EQ.2 ) THEN
408 C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1)
409 DO j=jMin,jMax
410 DO i=iMin,iMax
411 C Interpolate horizontally (but not vertically) vert viscosity to center:
412 C Although background visc. might be defined at k=1, this is not
413 C generally true when using variable visc. (from vertical mixing scheme).
414 C Therefore, no vert. interp. and only horizontal interpolation.
415 viscLoc = ( KappaRU(i,j,k) + KappaRU(i+1,j,k)
416 & +KappaRV(i,j,k) + KappaRV(i,j+1,k)
417 & )*0.25 _d 0
418 flxDisUp(i,j) =
419 & - viscLoc*( wVel(i,j, k ,bi,bj)
420 & -wVel(i,j,k-1,bi,bj) )*rkSign
421 & *recip_drF(k-1)*rA(i,j,bi,bj)
422 & *deepFac2C(k-1)*rhoFacC(k-1)
423 C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero)
424 c flxDisUp(i,j) = 0.
425 ENDDO
426 ENDDO
427 ENDIF
428 C Tendency is minus divergence of viscous fluxes:
429 C anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
430 DO j=jMin,jMax
431 DO i=iMin,iMax
432 gwDiss(i,j) =
433 & -( ( flx_EW(i+1,j)-flx_EW(i,j) )
434 & + ( flx_NS(i,j+1)-flx_NS(i,j) )
435 & + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
436 & *recip_rhoFacF(k)
437 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
438 & *recip_deepFac2F(k)
439 C-- prepare for next level (k+1)
440 flxDisUp(i,j)=flx_Dn(i,j)
441 ENDDO
442 ENDDO
443 ENDIF
444
445 IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN
446 C- No-slip BCs impose a drag at walls...
447 CALL MOM_W_SIDEDRAG(
448 I bi,bj,k,
449 I wVel, del2w,
450 I rThickC_C, recip_rThickC,
451 I viscAh_W, viscA4_W,
452 O gwAdd,
453 I myThid )
454 DO j=jMin,jMax
455 DO i=iMin,iMax
456 gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
457 ENDDO
458 ENDDO
459 ENDIF
460
461 #ifdef ALLOW_SMAG_3D
462 IF ( useSmag3D .AND. k.GT.1 ) THEN
463 CALL MOM_W_SMAG_3D(
464 I str13(1-OLx,1-OLy,1,bi,bj), str23(1-OLx,1-OLy,1,bi,bj),
465 I str33(1-OLx,1-OLy,1,bi,bj),
466 I viscAh3d_00(1-OLx,1-OLy,1,bi,bj),
467 I viscAh3d_13(1-OLx,1-OLy,1,bi,bj),
468 I viscAh3d_23(1-OLx,1-OLy,1,bi,bj),
469 I rThickC_W, rThickC_S, rThickC_C, recip_rThickC,
470 O gwAdd,
471 I k, bi, bj, myThid )
472 DO j=jMin,jMax
473 DO i=iMin,iMax
474 gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
475 ENDDO
476 ENDDO
477 ENDIF
478 #endif /* ALLOW_SMAG_3D */
479
480 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
481
482 IF ( momAdvection ) THEN
483
484 IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
485 C Advective Flux on Western face
486 DO j=jMin,jMax
487 DO i=iMin,iMax+1
488 C transport through Western face area:
489 uTrans = (
490 & drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj)
491 & *rhoFacC(km1)*mskM1
492 & + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
493 & *rhoFacC(k)
494 & )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
495 flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL
496 c flx_EW(i,j)=
497 c & uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
498 ENDDO
499 ENDDO
500 C Advective Flux on Southern face
501 DO j=jMin,jMax+1
502 DO i=iMin,iMax
503 C transport through Southern face area:
504 vTrans = (
505 & drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj)
506 & *rhoFacC(km1)*mskM1
507 & +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
508 & *rhoFacC(k)
509 & )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
510 flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL
511 c flx_NS(i,j)=
512 c & vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
513 ENDDO
514 ENDDO
515 ENDIF
516 C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
517 c IF (.TRUE.) THEN
518 DO j=jMin,jMax
519 DO i=iMin,iMax
520 C NH in p-coord.: advect wSpeed [m/s] with rTrans
521 tmp_WbarZ = halfRL*
522 & ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k )
523 & +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 )
524 C transport through Lower face area:
525 rTrans = halfRL*
526 & ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
527 & +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
528 & *mskP1
529 & )*rA(i,j,bi,bj)
530 flx_Dn(i,j) = rTrans*tmp_WbarZ
531 ENDDO
532 ENDDO
533 c ENDIF
534 IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN
535 C Advective Flux on Upper face of W-Cell (= at surface)
536 DO j=jMin,jMax
537 DO i=iMin,iMax
538 tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k)
539 rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k)
540 & *rA(i,j,bi,bj)
541 flxAdvUp(i,j) = rTrans*tmp_WbarZ
542 c flxAdvUp(i,j) = 0.
543 ENDDO
544 ENDDO
545 ENDIF
546
547 IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN
548 C Tendency is minus divergence of advective fluxes:
549 C anelastic: all transports & advect. fluxes are scaled by rhoFac
550 DO j=jMin,jMax
551 DO i=iMin,iMax
552 C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero)
553 c IF (k.EQ.2) flxAdvUp(i,j) = 0.
554 gW(i,j,k,bi,bj) =
555 & -( ( flx_EW(i+1,j)-flx_EW(i,j) )
556 & + ( flx_NS(i,j+1)-flx_NS(i,j) )
557 & + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
558 & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
559 & *recip_deepFac2F(k)*recip_rhoFacF(k)
560 ENDDO
561 ENDDO
562 #ifdef ALLOW_ADDFLUID
563 IF ( selectAddFluid.GE.1 ) THEN
564 DO j=jMin,jMax
565 DO i=iMin,iMax
566 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
567 & + wVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0
568 & *( addMass(i,j,k,bi,bj)
569 & +addMass(i,j,km1,bi,bj)*mskM1 )
570 & *recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
571 & *recip_deepFac2F(k)*recip_rhoFacF(k)
572 ENDDO
573 ENDDO
574 ENDIF
575 #endif /* ALLOW_ADDFLUID */
576 ENDIF
577
578 DO j=jMin,jMax
579 DO i=iMin,iMax
580 C-- prepare for next level (k+1)
581 flxAdvUp(i,j)=flx_Dn(i,j)
582 ENDDO
583 ENDDO
584
585 c ELSE
586 C- if momAdvection / else
587 c DO j=1-OLy,sNy+OLy
588 c DO i=1-OLx,sNx+OLx
589 c gW(i,j,k,bi,bj) = 0. _d 0
590 c ENDDO
591 c ENDDO
592
593 C- endif momAdvection.
594 ENDIF
595
596 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
597
598 IF ( useNHMTerms .AND. k.GT.1 ) THEN
599 CALL MOM_W_METRIC_NH(
600 I bi,bj,k,
601 I uVel, vVel,
602 O gwAdd,
603 I myThid )
604 DO j=jMin,jMax
605 DO i=iMin,iMax
606 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
607 ENDDO
608 ENDDO
609 ENDIF
610 IF ( use3dCoriolis .AND. k.GT.1 ) THEN
611 CALL MOM_W_CORIOLIS_NH(
612 I bi,bj,k,
613 I uVel, vVel,
614 O gwAdd,
615 I myThid )
616 DO j=jMin,jMax
617 DO i=iMin,iMax
618 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
619 ENDDO
620 ENDDO
621 ENDIF
622
623 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
624
625 #ifdef ALLOW_DIAGNOSTICS
626 IF ( diagDiss ) THEN
627 CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ',
628 & k, 1, 2, bi,bj, myThid )
629 C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
630 C does it only if k=1 (never the case here)
631 c IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid)
632 ENDIF
633 IF ( diagAdvec ) THEN
634 CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec',
635 & k,Nr, 1, bi,bj, myThid )
636 c IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid)
637 ENDIF
638 #endif /* ALLOW_DIAGNOSTICS */
639
640 C-- Dissipation term inside the Adams-Bashforth:
641 IF ( momViscosity .AND. momDissip_In_AB) THEN
642 DO j=jMin,jMax
643 DO i=iMin,iMax
644 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
645 ENDDO
646 ENDDO
647 ENDIF
648
649 C- Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
650 C and save gW_[n] into gwNm1 for the next time step.
651 #ifdef ALLOW_ADAMSBASHFORTH_3
652 CALL ADAMS_BASHFORTH3(
653 I bi, bj, k, Nr,
654 U gW, gwNm, gw_AB,
655 I nHydStartAB, myIter, myThid )
656 #else /* ALLOW_ADAMSBASHFORTH_3 */
657 CALL ADAMS_BASHFORTH2(
658 I bi, bj, k, Nr,
659 U gW, gwNm1, gw_AB,
660 I nHydStartAB, myIter, myThid )
661 #endif /* ALLOW_ADAMSBASHFORTH_3 */
662 #ifdef ALLOW_DIAGNOSTICS
663 IF ( diag_AB ) THEN
664 CALL DIAGNOSTICS_FILL(gw_AB,'AB_gW ',k,1,2,bi,bj,myThid)
665 ENDIF
666 #endif /* ALLOW_DIAGNOSTICS */
667
668 C-- Dissipation term outside the Adams-Bashforth:
669 IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
670 DO j=jMin,jMax
671 DO i=iMin,iMax
672 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
673 ENDDO
674 ENDDO
675 ENDIF
676
677 C- end of the k loop
678 ENDDO
679
680 #ifdef ALLOW_DIAGNOSTICS
681 IF (useDiagnostics) THEN
682 CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid)
683 CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid)
684 ENDIF
685 #endif /* ALLOW_DIAGNOSTICS */
686
687 #endif /* ALLOW_MOM_COMMON */
688 #endif /* ALLOW_NONHYDROSTATIC */
689
690 RETURN
691 END

  ViewVC Help
Powered by ViewVC 1.1.22