/[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.19 by adcroft, Tue May 31 14:49:38 2005 UTC revision 1.32 by jmc, Thu Jul 13 03:01:21 2006 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
3  C $Name$  C $Name$
4    
 #include "PACKAGES_CONFIG.h"  
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_GW(            SUBROUTINE CALC_GW(
11       I        myThid)       I               bi, bj, KappaRU, KappaRV,
12         I               myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | S/R CALC_GW                                                C     | S/R CALC_GW
16  C     | o Calculate vert. velocity tendency terms ( NH, QH only )  C     | o Calculate vertical velocity tendency terms
17    C     |   ( Non-Hydrostatic only )
18  C     *==========================================================*  C     *==========================================================*
19  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH, the vertical momentum tendency must be
20  C     | calculated explicitly and included as a source term  C     | calculated explicitly and included as a source term
21  C     | for a 3d pressure eqn. Calculate that term here.  C     | for a 3d pressure eqn. Calculate that term here.
22  C     | This routine is not used in HYD calculations.  C     | This routine is not used in HYD calculations.
23  C     *==========================================================*  C     *==========================================================*
24  C     \ev  C     \ev
25    
26  C     !USES:  C     !USES:
27        IMPLICIT NONE        IMPLICIT NONE
28  C     == Global variables ==  C     == Global variables ==
29  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
33  #include "GW.h"  #include "DYNVARS.h"
34  #include "CG3D.h"  #include "NH_VARS.h"
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38  C     myThid - Instance number for this innvocation of CALC_GW  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        INTEGER myThid
50    
51  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
52    
53  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
54  C     == Local variables ==  C     == Local variables ==
55  C     bi, bj,      :: Loop counters  C     iMin,iMax
56  C     iMin, iMax,  C     jMin,jMax
57  C     jMin, jMax  C     xA            :: W-Cell face area normal to X
58  C     flx_NS       :: Temp. used for fVol meridional terms.  C     yA            :: W-Cell face area normal to Y
59  C     flx_EW       :: Temp. used for fVol zonal terms.  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
60  C     flx_Up       :: Temp. used for fVol vertical terms.  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
61  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
62        INTEGER bi,bj,iMin,iMax,jMin,jMax  C     flx_NS        :: vertical momentum flux, meridional direction
63        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_EW        :: vertical momentum flux, zonal direction
64        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
65        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
66        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
67        _RL    fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     gwDiss        :: vertical momentum dissipation tendency
68        _RL    fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     i,j,k         :: Loop counters
69        _RL    del2w(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        INTEGER iMin,iMax,jMin,jMax
70  C     I,J,K - Loop counters        _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71        INTEGER i,j,k, kP1, kUp        _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        _RL  wOverride
85        _RS  hFacWtmp        _RL  tmp_WbarZ
86        _RS  hFacStmp        _RL  uTrans, vTrans, rTrans
87        _RS  hFacCtmp        _RL  viscLoc
88        _RS  recip_hFacCtmp        _RL  halfRL
89        _RL ab15,ab05        _RS  halfRS, zeroRS
90        _RL slipSideFac        PARAMETER( halfRL = 0.5D0 )
91        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
   
       _RL  Half  
       PARAMETER(Half=0.5D0)  
92  CEOP  CEOP
93    
94  ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  C Catch barotropic mode
95          IF ( Nr .LT. 2 ) RETURN
96    
97        iMin = 1        iMin = 1
98        iMax = sNx        iMax = sNx
99        jMin = 1        jMin = 1
100        jMax = sNy        jMax = sNy
101    
102  C     Adams-Bashforth timestepping weights  C--   Initialise gW to zero
103        ab15 =  1.5 _d 0 + abeps        DO k=1,Nr
104        ab05 = -0.5 _d 0 - abeps          DO j=1-OLy,sNy+OLy
105             DO i=1-OLx,sNx+OLx
 C     Lateral friction (no-slip, free slip, or half slip):  
       IF ( no_slip_sides ) THEN  
         slipSideFac = -1. _d 0  
       ELSE  
         slipSideFac =  1. _d 0  
       ENDIF  
 CML   half slip was used before ; keep the line for now, but half slip is  
 CML   not used anywhere in the code as far as I can see.  
 C        slipSideFac = 0. _d 0  
   
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO K=1,Nr  
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
            gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
106             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
           ENDDO  
107           ENDDO           ENDDO
108          ENDDO          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         ENDDO
115        ENDDO        ENDDO
116    
117  C Catch barotropic mode  C--   Boundaries condition at top
118        IF ( Nr .LT. 2 ) RETURN        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 For each tile  C---  Sweep down column
126        DO bj=myByLo(myThid),myByHi(myThid)        DO k=2,Nr
127         DO bi=myBxLo(myThid),myBxHi(myThid)          kp1=k+1
128            wOverRide=1.
129  C Boundaries condition at top          IF (k.EQ.Nr) THEN
130          DO J=jMin,jMax            kp1=Nr
131           DO I=iMin,iMax            wOverRide=0.
132            Flx_Dn(I,J,bi,bj)=0.          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           ENDDO
156          ENDDO          ENDDO
157    
158  C Sweep down column  C--   horizontal bi-harmonic dissipation
159          DO K=2,Nr          IF (momViscosity .AND. viscA4W.NE.0. ) THEN
160           Kp1=K+1  C-    calculate the horizontal Laplacian of vertical flow
          wOverRide=1.  
          if (K.EQ.Nr) then  
           Kp1=Nr  
           wOverRide=0.  
          endif  
 C     horizontal bi-harmonic dissipation  
          IF (momViscosity .AND. viscA4W.NE.0. ) THEN  
 C     calculate the horizontal Laplacian of vertical flow  
161  C     Zonal flux d/dx W  C     Zonal flux d/dx W
162            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
163             fZon(1-Olx,j)=0.             flx_EW(1-Olx,j)=0.
164             DO i=1-Olx+1,sNx+Olx             DO i=1-Olx+1,sNx+Olx
165              fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)              flx_EW(i,j) =
166       &           *_dyG(i,j,bi,bj)       &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
167       &           *_recip_dxC(i,j,bi,bj)       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
      &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))  
168  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
169       &           *sqcosFacU(J,bi,bj)       &              *sqcosFacU(j,bi,bj)
170  #endif  #endif
171             ENDDO             ENDDO
172            ENDDO            ENDDO
173  C     Meridional flux d/dy W  C     Meridional flux d/dy W
174            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
175             fMer(I,1-Oly)=0.             flx_NS(i,1-Oly)=0.
176            ENDDO            ENDDO
177            DO j=1-Oly+1,sNy+Oly            DO j=1-Oly+1,sNy+Oly
178             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
179              fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)              flx_NS(i,j) =
180       &           *_dxG(i,j,bi,bj)       &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
181       &           *_recip_dyC(i,j,bi,bj)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
      &           *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))  
182  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
183  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
184       &           *sqCosFacV(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
185  #endif  #endif
186  #endif  #endif
187             ENDDO             ENDDO
188            ENDDO            ENDDO
189              
190  C     del^2 W  C     del^2 W
191  C     Difference of zonal fluxes ...  C     Difference of zonal fluxes ...
192            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
193             DO i=1-Olx,sNx+Olx-1             DO i=1-Olx,sNx+Olx-1
194              del2w(i,j)=fZon(i+1,j)-fZon(i,j)              del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
195             ENDDO             ENDDO
196             del2w(sNx+Olx,j)=0.             del2w(sNx+Olx,j)=0.
197            ENDDO            ENDDO
# Line 174  C     Difference of zonal fluxes ... Line 199  C     Difference of zonal fluxes ...
199  C     ... add difference of meridional fluxes and divide by volume  C     ... add difference of meridional fluxes and divide by volume
200            DO j=1-Oly,sNy+Oly-1            DO j=1-Oly,sNy+Oly-1
201             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
202  C     First compute the fraction of open water for the w-control volume              del2w(i,j) = ( del2w(i,j)
203  C     at the southern face       &                    +(flx_NS(i,j+1)-flx_NS(i,j))
204              hFacCtmp=max(hFacC(I,J,K-1,bi,bj)-Half,0. _d 0)       &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
      &           +   min(hFacC(I,J,K  ,bi,bj),Half)  
             recip_hFacCtmp = 0. _d 0  
             IF (hFacCtmp .GT. 0.) THEN  
              recip_hFacCtmp = 1./hFacCtmp  
             ELSE  
              recip_hFacCtmp = 0. _d 0  
             ENDIF  
             del2w(i,j)=recip_rA(i,j,bi,bj)  
      &           *recip_drC(k)*recip_hFacCtmp  
      &           *(  
      &           del2w(i,j)  
      &           +( fMer(i,j+1)-fMer(i,j) )  
      &           )  
205             ENDDO             ENDDO
206            ENDDO            ENDDO
207  C-- No-slip BCs impose a drag at walls...  C-- No-slip BCs impose a drag at walls...
# Line 197  CML ************************************ Line 209  CML ************************************
209  CML   No-slip Boundary conditions for bi-harmonic dissipation  CML   No-slip Boundary conditions for bi-harmonic dissipation
210  CML   need to be implemented here!  CML   need to be implemented here!
211  CML ************************************************************  CML ************************************************************
212           ENDIF          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  C Flux on Southern face          IF (momViscosity) THEN
222           DO J=jMin,jMax+1  C Viscous Flux on Western face
223            DO I=iMin,iMax            DO j=jMin,jMax
224  C     First compute the fraction of open water for the w-control volume             DO i=iMin,iMax+1
225  C     at the southern face               flx_EW(i,j)=
226             hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
227       &         +    min(hFacS(I,J,K  ,bi,bj),Half)       &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
228             tmp_VbarZ=Half*(       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
229       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)  cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
230       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
231             Flx_NS(I,J,bi,bj)=       &              *(del2w(i,j)-del2w(i-1,j))
232       &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))       &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
233       &    -viscAhW*_recip_dyC(I,J,bi,bj)    cOld &              *_recip_dxC(i,j,bi,bj)*drC(k)
      &       *(hFacStmp*(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))  
      &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)  
      &         *wVel(I,J,K,bi,bj))  
      &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))  
 #ifdef ISOTROPIC_COS_SCALING  
234  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
235       &    *sqCosFacV(j,bi,bj)       &              *sqCosFacU(j,bi,bj)
236  #else  #else
237       &    *CosFacV(j,bi,bj)       &              *CosFacU(j,bi,bj)
238  #endif  #endif
239  #endif             ENDDO
 C     The last term is the weighted average of the viscous stress at the open  
 C     fraction of the w control volume and at the closed fraction of the  
 C     the control volume. A more compact but less intelligible version  
 C     of the last three lines is:  
 CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))  
 CML     &       *wVel(I,J,K,bi,bi) + hFacStmp*wVel(I,J-1,K,bi,bj) )  
240            ENDDO            ENDDO
241           ENDDO  C Viscous Flux on Southern face
242  C Flux on Western face            DO j=jMin,jMax+1
243           DO J=jMin,jMax             DO i=iMin,iMax
244            DO I=iMin,iMax+1               flx_NS(i,j)=
245  C     First compute the fraction of open water for the w-control volume       &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
246  C     at the western face       &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
247             hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
248       &         +    min(hFacW(I,J,K  ,bi,bj),Half)  cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
249                   tmp_UbarZ=Half*(       &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
250       &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)       &              *(del2w(i,j)-del2w(i,j-1))
251       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
252             Flx_EW(I,J,bi,bj)=  cOld &              *_recip_dyC(i,j,bi,bj)*drC(k)
253       &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))  #ifdef ISOTROPIC_COS_SCALING
      &    -viscAhW*_recip_dxC(I,J,bi,bj)  
      &      *(hFacWtmp*(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))  
      &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)  
      &        *wVel(I,J,K,bi,bj) )  
      &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))  
254  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
255       &                *sqCosFacU(j,bi,bj)       &              *sqCosFacV(j,bi,bj)
256  #else  #else
257       &                *CosFacU(j,bi,bj)       &              *CosFacV(j,bi,bj)
258    #endif
259  #endif  #endif
260  C     The last term is the weighted average of the viscous stress at the open             ENDDO
 C     fraction of the w control volume and at the closed fraction of the  
 C     the control volume. A more compact but less intelligible version  
 C     of the last three lines is:  
 CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))  
 CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )  
261            ENDDO            ENDDO
262           ENDDO  C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
263  C Flux on Lower face            DO j=jMin,jMax
264           DO J=jMin,jMax             DO i=iMin,iMax
265            DO I=iMin,iMax  C     Interpolate vert viscosity to center of tracer-cell (level k):
266             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)               viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
267             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)       &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
268       &         +wOverRide*wVel(I,J,Kp1,bi,bj))       &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
269             Flx_Dn(I,J,bi,bj)=       &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
270       &     tmp_WbarZ*tmp_WbarZ       &                 )*0.125 _d 0
271       &    -viscAr*recip_drF(K)               flx_Dn(i,j) =
272       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )       &          - 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            ENDDO
278           ENDDO  C     Tendency is minus divergence of viscous fluxes:
279  C        Divergence of fluxes            DO j=jMin,jMax
280           DO J=jMin,jMax             DO i=iMin,iMax
281            DO I=iMin,iMax               gwDiss(i,j) =
282             gW(I,J,K,bi,bj) = 0.       &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
283       &      -(       &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
284       &        +_recip_dxF(I,J,bi,bj)*(       &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
285       &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )       &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
286       &        +_recip_dyF(I,J,bi,bj)*(  cOld         gwDiss(i,j) =
287       &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )  cOld &        -(
288       &        +recip_drC(K)         *(  cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
289       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )  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  caja    * recip_hFacU(I,J,K,bi,bj)  c    &         )*recip_rThickC(i,j)
292  caja   NOTE:  This should be included    cOld &         )*recip_drC(k)
293  caja   but we need an hFacUW (above U points)  C--        prepare for next level (k+1)
294  caja           and an hFacUS (above V points) too...               flxDisUp(i,j)=flx_Dn(i,j)
295               ENDDO
296            ENDDO            ENDDO
297           ENDDO          ENDIF
         ENDDO  
        ENDDO  
       ENDDO  
298    
299                IF (no_slip_sides) THEN
300        DO bj=myByLo(myThid),myByHi(myThid)  C-     No-slip BCs impose a drag at walls...
301         DO bi=myBxLo(myThid),myBxHi(myThid)  c         CALL MOM_W_SIDEDRAG(
302          DO K=2,Nr  c    I        bi,bj,k,
303           DO j=jMin,jMax  c    O        gwAdd,
304            DO i=iMin,iMax  c    I        myThid)
305             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)  c         DO j=jMin,jMax
306       &     +deltatMom*nh_Am2*( ab15*gW(i,j,k,bi,bj)  c          DO i=iMin,iMax
307       &                 +ab05*gWNM1(i,j,k,bi,bj) )  c           gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
308             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.  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            ENDDO
328           ENDDO  C Advective Flux on Southern face
329          ENDDO            DO j=jMin,jMax+1
330         ENDDO             DO i=iMin,iMax
331        ENDDO  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  #ifdef ALLOW_OBCS          IF ( useNHMTerms ) THEN
374        IF (useOBCS) THEN            CALL MOM_W_METRIC_NH(
375  C--   This call is aesthetic: it makes the W field       I               bi,bj,k,
376  C     consistent with the OBs but this has no algorithmic       I               uVel, vVel,
377  C     impact. This is purely for diagnostic purposes.       O               gwAdd,
378         DO bj=myByLo(myThid),myByHi(myThid)       I               myThid )
379          DO bi=myBxLo(myThid),myBxHi(myThid)            DO j=jMin,jMax
380           DO K=1,Nr             DO i=iMin,iMax
381            CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )               gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
382           ENDDO             ENDDO
383          ENDDO            ENDDO
384         ENDDO          ENDIF
385        ENDIF          IF ( use3dCoriolis ) THEN
386  #endif /* ALLOW_OBCS */            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 */  #endif /* ALLOW_NONHYDROSTATIC */
436    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22