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

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

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

revision 1.32 by jmc, Thu Jul 13 03:01:21 2006 UTC revision 1.36 by jmc, Mon Mar 12 23:48:24 2007 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
 C     !DESCRIPTION: \bv  
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #define CALC_GW_NEW_THICK
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
# Line 30  C     == Global variables == Line 30  C     == Global variables ==
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
33    #include "SURFACE.h"
34  #include "DYNVARS.h"  #include "DYNVARS.h"
35  #include "NH_VARS.h"  #include "NH_VARS.h"
36    
# Line 56  C     iMin,iMax Line 57  C     iMin,iMax
57  C     jMin,jMax  C     jMin,jMax
58  C     xA            :: W-Cell face area normal to X  C     xA            :: W-Cell face area normal to X
59  C     yA            :: W-Cell face area normal to Y  C     yA            :: W-Cell face area normal to Y
60    C     rMinW,rMaxW   :: column boundaries (r-units) at Western  Edge
61    C     rMinS,rMaxS   :: column boundaries (r-units) at Southern Edge
62  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
63  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
64    C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
65  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
66  C     flx_NS        :: vertical momentum flux, meridional direction  C     flx_NS        :: vertical momentum flux, meridional direction
67  C     flx_EW        :: vertical momentum flux, zonal direction  C     flx_EW        :: vertical momentum flux, zonal direction
# Line 69  C     i,j,k         :: Loop counters Line 73  C     i,j,k         :: Loop counters
73        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
74        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76          _RS    rMinW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RS    rMaxW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RS    rMinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RS    rMaxS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 89  C     i,j,k         :: Loop counters Line 98  C     i,j,k         :: Loop counters
98        _RS  halfRS, zeroRS        _RS  halfRS, zeroRS
99        PARAMETER( halfRL = 0.5D0 )        PARAMETER( halfRL = 0.5D0 )
100        PARAMETER( halfRS = 0.5 , zeroRS = 0. )        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
101          PARAMETER( iMin = 1 , iMax = sNx )
102          PARAMETER( jMin = 1 , jMax = sNy )
103  CEOP  CEOP
104    
105  C Catch barotropic mode  C Catch barotropic mode
106        IF ( Nr .LT. 2 ) RETURN        IF ( Nr .LT. 2 ) RETURN
107    
       iMin = 1  
       iMax = sNx  
       jMin = 1  
       jMax = sNy  
   
108  C--   Initialise gW to zero  C--   Initialise gW to zero
109        DO k=1,Nr        DO k=1,Nr
110          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 121  C--   Boundaries condition at top Line 127  C--   Boundaries condition at top
127           flxDisUp(i,j) = 0.           flxDisUp(i,j) = 0.
128         ENDDO         ENDDO
129        ENDDO        ENDDO
130    C--   column boundaries :
131          IF (momViscosity) THEN
132            DO j=1-Oly,sNy+Oly
133             DO i=1-Olx+1,sNx+Olx
134               rMaxW(i,j) = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
135               rMinW(i,j) = MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
136             ENDDO
137            ENDDO
138            DO j=1-Oly+1,sNy+Oly
139             DO i=1-Olx,sNx+Olx
140               rMaxS(i,j) = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
141               rMinS(i,j) = MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
142             ENDDO
143            ENDDO
144          ENDIF
145    
146  C---  Sweep down column  C---  Sweep down column
147        DO k=2,Nr        DO k=2,Nr
# Line 131  C---  Sweep down column Line 152  C---  Sweep down column
152            wOverRide=0.            wOverRide=0.
153          ENDIF          ENDIF
154  C--   Compute grid factor arround a W-point:  C--   Compute grid factor arround a W-point:
155    #ifdef CALC_GW_NEW_THICK
156            DO j=1-Oly,sNy+Oly
157             DO i=1-Olx,sNx+Olx
158               IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
159         &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
160                 recip_rThickC(i,j) = 0.
161               ELSE
162    C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
163                 recip_rThickC(i,j) = 1. _d 0 /
164         &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
165         &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
166         &        )
167               ENDIF
168             ENDDO
169            ENDDO
170            IF (momViscosity) THEN
171             DO j=1-Oly,sNy+Oly
172              DO i=1-Olx,sNx+Olx
173               rThickC_C(i,j) = MAX( zeroRS,
174         &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
175         &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
176         &                         )
177              ENDDO
178             ENDDO
179             DO j=1-Oly,sNy+Oly
180              DO i=1-Olx+1,sNx+Olx
181               rThickC_W(i,j) = MAX( zeroRS,
182         &                           MIN( rMaxW(i,j), rC(k-1) )
183         &                          -MAX( rMinW(i,j), rC(k)   )
184         &                         )
185    C     W-Cell Western face area:
186               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
187    c    &                              *deepFacF(k)
188              ENDDO
189             ENDDO
190             DO j=1-Oly+1,sNy+Oly
191              DO i=1-Olx,sNx+Olx
192               rThickC_S(i,j) = MAX( zeroRS,
193         &                           MIN( rMaxS(i,j), rC(k-1) )
194         &                          -MAX( rMinS(i,j), rC(k)   )
195         &                         )
196    C     W-Cell Southern face area:
197               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
198    c    &                              *deepFacF(k)
199    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
200    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
201              ENDDO
202             ENDDO
203            ENDIF
204    #else /* CALC_GW_NEW_THICK */
205          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
206           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
207  C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !  C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
            rThickC_W(i,j) =  
      &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )  
      &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )  
            rThickC_S(i,j) =  
      &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )  
      &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )  
208             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN             IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
209               recip_rThickC(i,j) = 0.               recip_rThickC(i,j) = 0.
210             ELSE             ELSE
211               recip_rThickC(i,j) = 1. _d 0 /               recip_rThickC(i,j) = 1. _d 0 /
212       &        ( drF(k-1)*halfRS +       &        ( drF(k-1)*halfRS
213       &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )       &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
214       &        )       &        )
215             ENDIF             ENDIF
216    c          IF (momViscosity) THEN
217    #ifdef NONLIN_FRSURF
218               rThickC_C(i,j) =
219         &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
220         &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
221    #else
222               rThickC_C(i,j) =
223         &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
224         &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
225    #endif
226               rThickC_W(i,j) =
227         &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
228         &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
229               rThickC_S(i,j) =
230         &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
231         &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
232  C     W-Cell Western face area:  C     W-Cell Western face area:
233             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
234    c    &                              *deepFacF(k)
235  C     W-Cell Southern face area:  C     W-Cell Southern face area:
236             yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)             yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
237    c    &                              *deepFacF(k)
238    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
239    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
240    c          ENDIF
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243    #endif /* CALC_GW_NEW_THICK */
244    
245  C--   horizontal bi-harmonic dissipation  C--   horizontal bi-harmonic dissipation
246          IF (momViscosity .AND. viscA4W.NE.0. ) THEN          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
# Line 166  C     Zonal flux d/dx W Line 253  C     Zonal flux d/dx W
253       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
254       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
255  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
256       &              *sqcosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
257  #endif  #endif
258             ENDDO             ENDDO
259            ENDDO            ENDDO
# Line 202  C     ... add difference of meridional f Line 289  C     ... add difference of meridional f
289              del2w(i,j) = ( del2w(i,j)              del2w(i,j) = ( del2w(i,j)
290       &                    +(flx_NS(i,j+1)-flx_NS(i,j))       &                    +(flx_NS(i,j+1)-flx_NS(i,j))
291       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
292         &                    *recip_deepFac2F(k)
293             ENDDO             ENDDO
294            ENDDO            ENDDO
295  C-- No-slip BCs impose a drag at walls...  C-- No-slip BCs impose a drag at walls...
# Line 209  CML ************************************ Line 297  CML ************************************
297  CML   No-slip Boundary conditions for bi-harmonic dissipation  CML   No-slip Boundary conditions for bi-harmonic dissipation
298  CML   need to be implemented here!  CML   need to be implemented here!
299  CML ************************************************************  CML ************************************************************
300          ELSE          ELSEIF (momViscosity) THEN
301  C-    Initialize del2w to zero:  C-    Initialize del2w to zero:
302            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
303             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
# Line 227  C Viscous Flux on Western face Line 315  C Viscous Flux on Western face
315       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
316       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
317  cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)  cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
318         &              *cosFacU(j,bi,bj)
319       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
320       &              *(del2w(i,j)-del2w(i-1,j))       &              *(del2w(i,j)-del2w(i-1,j))
321       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
# Line 234  cOld &              *_recip_dxC(i,j,bi,b Line 323  cOld &              *_recip_dxC(i,j,bi,b
323  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
324       &              *sqCosFacU(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
325  #else  #else
326       &              *CosFacU(j,bi,bj)       &              *cosFacU(j,bi,bj)
327  #endif  #endif
328             ENDDO             ENDDO
329            ENDDO            ENDDO
# Line 246  C Viscous Flux on Southern face Line 335  C Viscous Flux on Southern face
335       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
336       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
337  cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)  cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
338    #ifdef ISOTROPIC_COS_SCALING
339         &              *cosFacV(j,bi,bj)
340    #endif
341       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
342       &              *(del2w(i,j)-del2w(i,j-1))       &              *(del2w(i,j)-del2w(i,j-1))
343       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
# Line 254  cOld &              *_recip_dyC(i,j,bi,b Line 346  cOld &              *_recip_dyC(i,j,bi,b
346  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
347       &              *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
348  #else  #else
349       &              *CosFacV(j,bi,bj)       &              *cosFacV(j,bi,bj)
350  #endif  #endif
351  #endif  #endif
352             ENDDO             ENDDO
# Line 272  C     Interpolate vert viscosity to cent Line 364  C     Interpolate vert viscosity to cent
364       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide       &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
365       &                     -wVel(i,j, k ,bi,bj) )*rkSign       &                     -wVel(i,j, k ,bi,bj) )*rkSign
366       &                   *recip_drF(k)*rA(i,j,bi,bj)       &                   *recip_drF(k)*rA(i,j,bi,bj)
367         &                   *deepFac2C(k)*rhoFacC(k)
368  cOld &                   *recip_drF(k)  cOld &                   *recip_drF(k)
369             ENDDO             ENDDO
370            ENDDO            ENDDO
371  C     Tendency is minus divergence of viscous fluxes:  C     Tendency is minus divergence of viscous fluxes:
372    C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
373            DO j=jMin,jMax            DO j=jMin,jMax
374             DO i=iMin,iMax             DO i=iMin,iMax
375               gwDiss(i,j) =               gwDiss(i,j) =
376       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
377       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
378       &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign       &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
379         &                                          *recip_rhoFacF(k)
380       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
381         &          *recip_deepFac2F(k)
382  cOld         gwDiss(i,j) =  cOld         gwDiss(i,j) =
383  cOld &        -(  cOld &        -(
384  cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )  cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
# Line 296  C--        prepare for next level (k+1) Line 392  C--        prepare for next level (k+1)
392            ENDDO            ENDDO
393          ENDIF          ENDIF
394    
395          IF (no_slip_sides) THEN          IF ( momViscosity .AND. no_slip_sides ) THEN
396  C-     No-slip BCs impose a drag at walls...  C-     No-slip BCs impose a drag at walls...
397  c         CALL MOM_W_SIDEDRAG(            CALL MOM_W_SIDEDRAG(
398  c    I        bi,bj,k,       I               bi,bj,k,
399  c    O        gwAdd,       I               wVel, del2w,
400  c    I        myThid)       I               rThickC_C, recip_rThickC,
401  c         DO j=jMin,jMax       I               viscAh_W, viscA4_W,
402  c          DO i=iMin,iMax       O               gwAdd,
403  c           gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)       I               myThid )
404  c          ENDDO            DO j=jMin,jMax
405  c         ENDDO             DO i=iMin,iMax
406                gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
407               ENDDO
408              ENDDO
409          ENDIF          ENDIF
410    
411  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 318  C Advective Flux on Western face Line 417  C Advective Flux on Western face
417  C     transport through Western face area:  C     transport through Western face area:
418               uTrans = (               uTrans = (
419       &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)       &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)
420         &                  *rhoFacC(k-1)
421       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)       &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
422       &                )*halfRL*_dyG(i,j,bi,bj)       &                  *rhoFacC(k)
423         &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
424  cOld &                )*halfRL  cOld &                )*halfRL
425               flx_EW(i,j)=               flx_EW(i,j)=
426       &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL       &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
# Line 331  C Advective Flux on Southern face Line 432  C Advective Flux on Southern face
432  C     transport through Southern face area:  C     transport through Southern face area:
433               vTrans = (               vTrans = (
434       &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)       &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)
435         &                  *rhoFacC(k-1)
436       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)       &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
437       &                )*halfRL*_dxG(i,j,bi,bj)       &                  *rhoFacC(k)
438         &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
439  cOld &                )*halfRL  cOld &                )*halfRL
440               flx_NS(i,j)=               flx_NS(i,j)=
441       &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL       &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
# Line 341  cOld &                )*halfRL Line 444  cOld &                )*halfRL
444  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)  C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
445            DO j=jMin,jMax            DO j=jMin,jMax
446             DO i=iMin,iMax             DO i=iMin,iMax
447               tmp_WbarZ = halfRL*( wVel(i,j, k ,bi,bj)  C     NH in p-coord.: advect wSpeed [m/s] with rTrans
448       &                           +wVel(i,j,kp1,bi,bj)*wOverRide )               tmp_WbarZ = halfRL*
449         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)
450         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*wOverRide )
451  C     transport through Lower face area:  C     transport through Lower face area:
452               rTrans = tmp_WbarZ*rA(i,j,bi,bj)               rTrans = halfRL*
453         &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
454         &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
455         &                                   *wOverRide
456         &              )*rA(i,j,bi,bj)
457               flx_Dn(i,j) = rTrans*tmp_WbarZ               flx_Dn(i,j) = rTrans*tmp_WbarZ
458  cOld         flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ  cOld         flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ
459             ENDDO             ENDDO
460            ENDDO            ENDDO
461  C     Tendency is minus divergence of advective fluxes:  C     Tendency is minus divergence of advective fluxes:
462    C     anelastic: all transports & advect. fluxes are scaled by rhoFac
463            DO j=jMin,jMax            DO j=jMin,jMax
464             DO i=iMin,iMax             DO i=iMin,iMax
465               gW(i,j,k,bi,bj) =               gW(i,j,k,bi,bj) =
466       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
467       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
468       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign       &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
469       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
470         &          *recip_deepFac2F(k)*recip_rhoFacF(k)
471  cOld         gW(i,j,k,bi,bj) =  cOld         gW(i,j,k,bi,bj) =
472  cOld &        -(  cOld &        -(
473  cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )  cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
Legend:
Removed from v.1.32  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22