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

Annotation 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.2 - (hide annotations) (download)
Tue Nov 5 13:53:26 2013 UTC (10 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
- move to pkg/mom_common and model/src (previously in tutorial_deep_convection
  code) 2nd version of isotropic 3-D Smagorinsky code interface: strain and
  viscosity are locally declared in dynmics.F and pass as argument to CALC_GW;
  ensure that all field value that are used are set.

1 jmc 1.2 C $Header: /u/gcmpack/MITgcm/verification/tutorial_deep_convection/code/calc_gw.F,v 1.1 2013/09/30 18:19:41 jmc Exp $
2 jmc 1.1 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