/[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.17 by mlosch, Fri Oct 1 16:15:29 2004 UTC revision 1.37 by jmc, Fri Oct 19 14:41:39 2007 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
 C     !DESCRIPTION: \bv  
2  C $Name$  C $Name$
3    
 #include "PACKAGES_CONFIG.h"  
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
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 "RESTART.h"
34  #include "CG3D.h"  #include "SURFACE.h"
35    #include "DYNVARS.h"
36    #include "NH_VARS.h"
37    
38  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
39  C     == Routine arguments ==  C     == Routine arguments ==
40  C     myThid - Instance number for this innvocation of CALC_GW  C     bi,bj   :: current tile indices
41    C     KappaRU :: vertical viscosity at U points
42    C     KappaRV :: vertical viscosity at V points
43    C     myTime  :: Current time in simulation
44    C     myIter  :: Current iteration number in simulation
45    C     myThid  :: Thread number for this instance of the routine.
46          INTEGER bi,bj
47          _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48          _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49          _RL     myTime
50          INTEGER myIter
51        INTEGER myThid        INTEGER myThid
52    
53  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
54    
55  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
56  C     == Local variables ==  C     == Local variables ==
57  C     bi, bj,      :: Loop counters  C     iMin,iMax
58  C     iMin, iMax,  C     jMin,jMax
59  C     jMin, jMax  C     xA            :: W-Cell face area normal to X
60  C     flx_NS       :: Temp. used for fVol meridional terms.  C     yA            :: W-Cell face area normal to Y
61  C     flx_EW       :: Temp. used for fVol zonal terms.  C     rMinW,rMaxW   :: column boundaries (r-units) at Western  Edge
62  C     flx_Up       :: Temp. used for fVol vertical terms.  C     rMinS,rMaxS   :: column boundaries (r-units) at Southern Edge
63  C     flx_Dn       :: Temp. used for fVol vertical terms.  C     rThickC_W     :: thickness (in r-units) of W-Cell at Western Edge
64        INTEGER bi,bj,iMin,iMax,jMin,jMax  C     rThickC_S     :: thickness (in r-units) of W-Cell at Southern Edge
65        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     rThickC_C     :: thickness (in r-units) of W-Cell (centered on W pt)
66        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     recip_rThickC :: reciprol thickness of W-Cell (centered on W-point)
67        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_NS        :: vertical momentum flux, meridional direction
68        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  C     flx_EW        :: vertical momentum flux, zonal direction
69  C     I,J,K - Loop counters  C     flxAdvUp      :: vertical mom. advective   flux, vertical direction (@ level k-1)
70        INTEGER i,j,k, kP1, kUp  C     flxDisUp      :: vertical mom. dissipation flux, vertical direction (@ level k-1)
71    C     flx_Dn        :: vertical momentum flux, vertical direction (@ level k)
72    C     gwDiss        :: vertical momentum dissipation tendency
73    C     i,j,k         :: Loop counters
74          INTEGER iMin,iMax,jMin,jMax
75          _RS    xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76          _RS    yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77          _RS    rMinW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78          _RS    rMaxW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79          _RS    rMinS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80          _RS    rMaxS (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL    rThickC_W    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL    rThickC_S    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL    rThickC_C    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL    recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86          _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL    flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL    flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RL    gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL    gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL    del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          INTEGER i,j,k, kp1
94        _RL  wOverride        _RL  wOverride
95        _RS  hFacWtmp        _RL  tmp_WbarZ
96        _RS  hFacStmp        _RL  uTrans, vTrans, rTrans
97        _RL ab15,ab05        _RL  viscLoc
98        _RL slipSideFac        _RL  halfRL
99        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RS  halfRS, zeroRS
100          PARAMETER( halfRL = 0.5D0 )
101        _RL  Half        PARAMETER( halfRS = 0.5 , zeroRS = 0. )
102        PARAMETER(Half=0.5D0)        PARAMETER( iMin = 1 , iMax = sNx )
103          PARAMETER( jMin = 1 , jMax = sNy )
104  CEOP  CEOP
105    
106  ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  C Catch barotropic mode
107          IF ( Nr .LT. 2 ) RETURN
108    
109        iMin = 1  C--   Initialise gW to zero
110        iMax = sNx        DO k=1,Nr
111        jMin = 1          DO j=1-OLy,sNy+OLy
112        jMax = sNy           DO i=1-OLx,sNx+OLx
   
 C     Adams-Bashforth timestepping weights  
       ab15 =  1.5 _d 0 + abeps  
       ab05 = -0.5 _d 0 - abeps  
   
 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 it 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)  
113             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
           ENDDO  
114           ENDDO           ENDDO
115          ENDDO          ENDDO
116          ENDDO
117    C-    Initialise gwDiss to zero
118          DO j=1-OLy,sNy+OLy
119           DO i=1-OLx,sNx+OLx
120             gwDiss(i,j) = 0.
121         ENDDO         ENDDO
122        ENDDO        ENDDO
123    
124  C Catch barotropic mode  C--   Boundaries condition at top
125        IF ( Nr .LT. 2 ) RETURN        DO j=1-OLy,sNy+OLy
126           DO i=1-OLx,sNx+OLx
127  C For each tile           flxAdvUp(i,j) = 0.
128        DO bj=myByLo(myThid),myByHi(myThid)           flxDisUp(i,j) = 0.
129         DO bi=myBxLo(myThid),myBxHi(myThid)         ENDDO
130          ENDDO
131  C Boundaries condition at top  C--   column boundaries :
132          DO J=jMin,jMax        IF (momViscosity) THEN
133           DO I=iMin,iMax          DO j=1-Oly,sNy+Oly
134            Flx_Dn(I,J,bi,bj)=0.           DO i=1-Olx+1,sNx+Olx
135               rMaxW(i,j) = MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
136               rMinW(i,j) = MAX(   R_low(i-1,j,bi,bj),   R_low(i,j,bi,bj) )
137             ENDDO
138            ENDDO
139            DO j=1-Oly+1,sNy+Oly
140             DO i=1-Olx,sNx+Olx
141               rMaxS(i,j) = MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
142               rMinS(i,j) = MAX(   R_low(i,j-1,bi,bj),   R_low(i,j,bi,bj) )
143           ENDDO           ENDDO
144          ENDDO          ENDDO
145          ENDIF
146    
147  C Sweep down column  C---  Sweep down column
148          DO K=2,Nr        DO k=2,Nr
149           Kp1=K+1          kp1=k+1
150           wOverRide=1.          wOverRide=1.
151           if (K.EQ.Nr) then          IF (k.EQ.Nr) THEN
152            Kp1=Nr            kp1=Nr
153            wOverRide=0.            wOverRide=0.
154           endif          ENDIF
155  C Flux on Southern face  C--   Compute grid factor arround a W-point:
156           DO J=jMin,jMax+1  #ifdef CALC_GW_NEW_THICK
157            DO I=iMin,iMax          DO j=1-Oly,sNy+Oly
158  C     First compute the fraction of open water for the w-control volume           DO i=1-Olx,sNx+Olx
159  C     at the southern face             IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR.
160             hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)       &          maskC(i,j, k ,bi,bj).EQ.0. ) THEN
161       &         +    min(hFacS(I,J,K  ,bi,bj),Half)               recip_rThickC(i,j) = 0.
162             tmp_VbarZ=Half*(             ELSE
163       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)  C-    valid in z & p coord.; also accurate if Interface @ middle between 2 centers
164       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))               recip_rThickC(i,j) = 1. _d 0 /
165             Flx_NS(I,J,bi,bj)=       &        (  MIN( Ro_surf(i,j,bi,bj),rC(k-1) )
166       &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))       &         - MAX( R_low(i,j,bi,bj),  rC(k)   )
167       &    -viscAhW*_recip_dyC(I,J,bi,bj)         &        )
168       &       *(hFacStmp*(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))             ENDIF
      &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)  
      &         *wVel(I,J,K,bi,bj))  
 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) )  
           ENDDO  
          ENDDO  
 C Flux on Western face  
          DO J=jMin,jMax  
           DO I=iMin,iMax+1  
 C     First compute the fraction of open water for the w-control volume  
 C     at the western face  
            hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)  
      &         +    min(hFacW(I,J,K  ,bi,bj),Half)  
                  tmp_UbarZ=Half*(  
      &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)  
      &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))  
            Flx_EW(I,J,bi,bj)=  
      &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))  
      &    -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) )  
 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 - hFacWtmp))  
 CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )  
           ENDDO  
          ENDDO  
 C Flux on Lower face  
          DO J=jMin,jMax  
           DO I=iMin,iMax  
            Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)  
            tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)  
      &         +wOverRide*wVel(I,J,Kp1,bi,bj))  
            Flx_Dn(I,J,bi,bj)=  
      &     tmp_WbarZ*tmp_WbarZ  
      &    -viscAr*recip_drF(K)  
      &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )  
           ENDDO  
          ENDDO  
 C        Divergence of fluxes  
          DO J=jMin,jMax  
           DO I=iMin,iMax  
            gW(I,J,K,bi,bj) = 0.  
      &      -(  
      &        +_recip_dxF(I,J,bi,bj)*(  
      &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )  
      &        +_recip_dyF(I,J,bi,bj)*(  
      &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )  
      &        +recip_drC(K)         *(  
      &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )  
      &       )  
 caja    * recip_hFacU(I,J,K,bi,bj)  
 caja   NOTE:  This should be included    
 caja   but we need an hFacUW (above U points)  
 caja           and an hFacUS (above V points) too...  
           ENDDO  
169           ENDDO           ENDDO
170          ENDDO          ENDDO
171         ENDDO          IF (momViscosity) THEN
172        ENDDO           DO j=1-Oly,sNy+Oly
173              DO i=1-Olx,sNx+Olx
174                   rThickC_C(i,j) = MAX( zeroRS,
175        DO bj=myByLo(myThid),myByHi(myThid)       &                           MIN( Ro_surf(i,j,bi,bj), rC(k-1) )
176         DO bi=myBxLo(myThid),myBxHi(myThid)       &                          -MAX(   R_low(i,j,bi,bj),  rC(k)  )
177          DO K=2,Nr       &                         )
          DO j=jMin,jMax  
           DO i=iMin,iMax  
            wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)  
      &     +deltatMom*( ab15*gW(i,j,k,bi,bj)  
      &                 +ab05*gWNM1(i,j,k,bi,bj) )  
            IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.  
178            ENDDO            ENDDO
179           ENDDO           ENDDO
180          ENDDO           DO j=1-Oly,sNy+Oly
181         ENDDO            DO i=1-Olx+1,sNx+Olx
182        ENDDO             rThickC_W(i,j) = MAX( zeroRS,
183         &                           MIN( rMaxW(i,j), rC(k-1) )
184  #ifdef ALLOW_OBCS       &                          -MAX( rMinW(i,j), rC(k)   )
185        IF (useOBCS) THEN       &                         )
186  C--   This call is aesthetic: it makes the W field  C     W-Cell Western face area:
187  C     consistent with the OBs but this has no algorithmic             xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
188  C     impact. This is purely for diagnostic purposes.  c    &                              *deepFacF(k)
189         DO bj=myByLo(myThid),myByHi(myThid)            ENDDO
190          DO bi=myBxLo(myThid),myBxHi(myThid)           ENDDO
191           DO K=1,Nr           DO j=1-Oly+1,sNy+Oly
192            CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )            DO i=1-Olx,sNx+Olx
193               rThickC_S(i,j) = MAX( zeroRS,
194         &                           MIN( rMaxS(i,j), rC(k-1) )
195         &                          -MAX( rMinS(i,j), rC(k)   )
196         &                         )
197    C     W-Cell Southern face area:
198               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
199    c    &                              *deepFacF(k)
200    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
201    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
202              ENDDO
203             ENDDO
204            ENDIF
205    #else /* CALC_GW_NEW_THICK */
206            DO j=1-Oly,sNy+Oly
207             DO i=1-Olx,sNx+Olx
208    C-    note: assume fluid @ smaller k than bottom: does not work in p-coordinate !
209               IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN
210                 recip_rThickC(i,j) = 0.
211               ELSE
212                 recip_rThickC(i,j) = 1. _d 0 /
213         &        ( drF(k-1)*halfRS
214         &        + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS )
215         &        )
216               ENDIF
217    c          IF (momViscosity) THEN
218    #ifdef NONLIN_FRSURF
219               rThickC_C(i,j) =
220         &          drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
221         &        + drF( k )*MIN( h0FacC(i,j,k  ,bi,bj), halfRS )
222    #else
223               rThickC_C(i,j) =
224         &          drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS )
225         &        + drF( k )*MIN( _hFacC(i,j,k  ,bi,bj), halfRS )
226    #endif
227               rThickC_W(i,j) =
228         &          drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS )
229         &        + drF( k )*MIN( _hFacW(i,j,k  ,bi,bj), halfRS )
230               rThickC_S(i,j) =
231         &          drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS )
232         &        + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS )
233    C     W-Cell Western face area:
234               xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j)
235    c    &                              *deepFacF(k)
236    C     W-Cell Southern face area:
237               yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j)
238    c    &                              *deepFacF(k)
239    C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC.
240    C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted)
241    c          ENDIF
242           ENDDO           ENDDO
243          ENDDO          ENDDO
244         ENDDO  #endif /* CALC_GW_NEW_THICK */
245        ENDIF  
246  #endif /* ALLOW_OBCS */  C--   horizontal bi-harmonic dissipation
247            IF (momViscosity .AND. viscA4W.NE.0. ) THEN
248    C-    calculate the horizontal Laplacian of vertical flow
249    C     Zonal flux d/dx W
250              DO j=1-Oly,sNy+Oly
251               flx_EW(1-Olx,j)=0.
252               DO i=1-Olx+1,sNx+Olx
253                flx_EW(i,j) =
254         &              (wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
255         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
256    #ifdef COSINEMETH_III
257         &              *sqCosFacU(j,bi,bj)
258    #endif
259               ENDDO
260              ENDDO
261    C     Meridional flux d/dy W
262              DO i=1-Olx,sNx+Olx
263               flx_NS(i,1-Oly)=0.
264              ENDDO
265              DO j=1-Oly+1,sNy+Oly
266               DO i=1-Olx,sNx+Olx
267                flx_NS(i,j) =
268         &              (wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
269         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
270    #ifdef ISOTROPIC_COS_SCALING
271    #ifdef COSINEMETH_III
272         &              *sqCosFacV(j,bi,bj)
273    #endif
274    #endif
275               ENDDO
276              ENDDO
277    
278    C     del^2 W
279    C     Difference of zonal fluxes ...
280              DO j=1-Oly,sNy+Oly
281               DO i=1-Olx,sNx+Olx-1
282                del2w(i,j)=flx_EW(i+1,j)-flx_EW(i,j)
283               ENDDO
284               del2w(sNx+Olx,j)=0.
285              ENDDO
286    
287    C     ... add difference of meridional fluxes and divide by volume
288              DO j=1-Oly,sNy+Oly-1
289               DO i=1-Olx,sNx+Olx
290                del2w(i,j) = ( del2w(i,j)
291         &                    +(flx_NS(i,j+1)-flx_NS(i,j))
292         &                   )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
293         &                    *recip_deepFac2F(k)
294               ENDDO
295              ENDDO
296    C-- No-slip BCs impose a drag at walls...
297    CML ************************************************************
298    CML   No-slip Boundary conditions for bi-harmonic dissipation
299    CML   need to be implemented here!
300    CML ************************************************************
301            ELSEIF (momViscosity) THEN
302    C-    Initialize del2w to zero:
303              DO j=1-Oly,sNy+Oly
304               DO i=1-Olx,sNx+Olx
305                del2w(i,j) = 0. _d 0
306               ENDDO
307              ENDDO
308            ENDIF
309    
310            IF (momViscosity) THEN
311    C Viscous Flux on Western face
312              DO j=jMin,jMax
313               DO i=iMin,iMax+1
314                 flx_EW(i,j)=
315         &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL
316         &              *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
317         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
318    cOld &              *_recip_dxC(i,j,bi,bj)*rThickC_W(i,j)
319         &              *cosFacU(j,bi,bj)
320         &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL
321         &              *(del2w(i,j)-del2w(i-1,j))
322         &              *_recip_dxC(i,j,bi,bj)*xA(i,j)
323    cOld &              *_recip_dxC(i,j,bi,bj)*drC(k)
324    #ifdef COSINEMETH_III
325         &              *sqCosFacU(j,bi,bj)
326    #else
327         &              *cosFacU(j,bi,bj)
328    #endif
329               ENDDO
330              ENDDO
331    C Viscous Flux on Southern face
332              DO j=jMin,jMax+1
333               DO i=iMin,iMax
334                 flx_NS(i,j)=
335         &       - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL
336         &              *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
337         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
338    cOld &              *_recip_dyC(i,j,bi,bj)*rThickC_S(i,j)
339    #ifdef ISOTROPIC_COS_SCALING
340         &              *cosFacV(j,bi,bj)
341    #endif
342         &       + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL
343         &              *(del2w(i,j)-del2w(i,j-1))
344         &              *_recip_dyC(i,j,bi,bj)*yA(i,j)
345    cOld &              *_recip_dyC(i,j,bi,bj)*drC(k)
346    #ifdef ISOTROPIC_COS_SCALING
347    #ifdef COSINEMETH_III
348         &              *sqCosFacV(j,bi,bj)
349    #else
350         &              *cosFacV(j,bi,bj)
351    #endif
352    #endif
353               ENDDO
354              ENDDO
355    C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k)
356              DO j=jMin,jMax
357               DO i=iMin,iMax
358    C     Interpolate vert viscosity to center of tracer-cell (level k):
359                 viscLoc = ( KappaRU(i,j,k)  +KappaRU(i+1,j,k)
360         &                  +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
361         &                  +KappaRV(i,j,k)  +KappaRV(i,j+1,k)
362         &                  +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
363         &                 )*0.125 _d 0
364                 flx_Dn(i,j) =
365         &          - viscLoc*( wVel(i,j,kp1,bi,bj)*wOverRide
366         &                     -wVel(i,j, k ,bi,bj) )*rkSign
367         &                   *recip_drF(k)*rA(i,j,bi,bj)
368         &                   *deepFac2C(k)*rhoFacC(k)
369    cOld &                   *recip_drF(k)
370               ENDDO
371              ENDDO
372    C     Tendency is minus divergence of viscous fluxes:
373    C     anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not
374              DO j=jMin,jMax
375               DO i=iMin,iMax
376                 gwDiss(i,j) =
377         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
378         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
379         &           + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign
380         &                                          *recip_rhoFacF(k)
381         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
382         &          *recip_deepFac2F(k)
383    cOld         gwDiss(i,j) =
384    cOld &        -(
385    cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
386    cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
387    cOld &          +                      ( flxDisUp(i,j)-flx_Dn(i,j) )
388    c    &         )*recip_rThickC(i,j)
389    cOld &         )*recip_drC(k)
390    C--        prepare for next level (k+1)
391                 flxDisUp(i,j)=flx_Dn(i,j)
392               ENDDO
393              ENDDO
394            ENDIF
395    
396            IF ( momViscosity .AND. no_slip_sides ) THEN
397    C-     No-slip BCs impose a drag at walls...
398              CALL MOM_W_SIDEDRAG(
399         I               bi,bj,k,
400         I               wVel, del2w,
401         I               rThickC_C, recip_rThickC,
402         I               viscAh_W, viscA4_W,
403         O               gwAdd,
404         I               myThid )
405              DO j=jMin,jMax
406               DO i=iMin,iMax
407                gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j)
408               ENDDO
409              ENDDO
410            ENDIF
411    
412    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
413    
414            IF ( momAdvection ) THEN
415    C Advective Flux on Western face
416              DO j=jMin,jMax
417               DO i=iMin,iMax+1
418    C     transport through Western face area:
419                 uTrans = (
420         &          drF(k-1)*_hFacW(i,j,k-1,bi,bj)*uVel(i,j,k-1,bi,bj)
421         &                  *rhoFacC(k-1)
422         &        + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj)
423         &                  *rhoFacC(k)
424         &                )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k)
425    cOld &                )*halfRL
426                 flx_EW(i,j)=
427         &         uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL
428               ENDDO
429              ENDDO
430    C Advective Flux on Southern face
431              DO j=jMin,jMax+1
432               DO i=iMin,iMax
433    C     transport through Southern face area:
434                 vTrans = (
435         &          drF(k-1)*_hFacS(i,j,k-1,bi,bj)*vVel(i,j,k-1,bi,bj)
436         &                  *rhoFacC(k-1)
437         &         +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj)
438         &                  *rhoFacC(k)
439         &                )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k)
440    cOld &                )*halfRL
441                 flx_NS(i,j)=
442         &         vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL
443               ENDDO
444              ENDDO
445    C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k)
446              DO j=jMin,jMax
447               DO i=iMin,iMax
448    C     NH in p-coord.: advect wSpeed [m/s] with rTrans
449                 tmp_WbarZ = halfRL*
450         &              ( wVel(i,j, k ,bi,bj)*rVel2wUnit(k)
451         &               +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*wOverRide )
452    C     transport through Lower face area:
453                 rTrans = halfRL*
454         &              ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k )
455         &               +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1)
456         &                                   *wOverRide
457         &              )*rA(i,j,bi,bj)
458                 flx_Dn(i,j) = rTrans*tmp_WbarZ
459    cOld         flx_Dn(i,j) = tmp_WbarZ*tmp_WbarZ
460               ENDDO
461              ENDDO
462    C     Tendency is minus divergence of advective fluxes:
463    C     anelastic: all transports & advect. fluxes are scaled by rhoFac
464              DO j=jMin,jMax
465               DO i=iMin,iMax
466                 gW(i,j,k,bi,bj) =
467         &        -(   ( flx_EW(i+1,j)-flx_EW(i,j) )
468         &           + ( flx_NS(i,j+1)-flx_NS(i,j) )
469         &           + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k)
470         &         )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j)
471         &          *recip_deepFac2F(k)*recip_rhoFacF(k)
472    cOld         gW(i,j,k,bi,bj) =
473    cOld &        -(
474    cOld &          +_recip_dxF(i,j,bi,bj)*( flx_EW(i+1,j)-flx_EW(i,j) )
475    cOld &          +_recip_dyF(i,j,bi,bj)*( flx_NS(i,j+1)-flx_NS(i,j) )
476    cOld &          +                      ( flxAdvUp(i,j)-flx_Dn(i,j) )
477    c    &         )*recip_rThickC(i,j)
478    cOld &         )*recip_drC(k)
479    C--          prepare for next level (k+1)
480                 flxAdvUp(i,j)=flx_Dn(i,j)
481               ENDDO
482              ENDDO
483            ENDIF
484    
485            IF ( useNHMTerms ) THEN
486              CALL MOM_W_METRIC_NH(
487         I               bi,bj,k,
488         I               uVel, vVel,
489         O               gwAdd,
490         I               myThid )
491              DO j=jMin,jMax
492               DO i=iMin,iMax
493                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
494               ENDDO
495              ENDDO
496            ENDIF
497            IF ( use3dCoriolis ) THEN
498              CALL MOM_W_CORIOLIS_NH(
499         I               bi,bj,k,
500         I               uVel, vVel,
501         O               gwAdd,
502         I               myThid )
503              DO j=jMin,jMax
504               DO i=iMin,iMax
505                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j)
506               ENDDO
507              ENDDO
508            ENDIF
509    
510    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
511    
512    C--   Dissipation term inside the Adams-Bashforth:
513            IF ( momViscosity .AND. momDissip_In_AB) THEN
514              DO j=jMin,jMax
515               DO i=iMin,iMax
516                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
517               ENDDO
518              ENDDO
519            ENDIF
520    
521    C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
522    C     and save gW_[n] into gwNm1 for the next time step.
523    c#ifdef ALLOW_ADAMSBASHFORTH_3
524    c       CALL ADAMS_BASHFORTH3(
525    c    I                         bi, bj, k,
526    c    U                         gW, gwNm,
527    c    I                         nHydStartAB, myIter, myThid )
528    c#else /* ALLOW_ADAMSBASHFORTH_3 */
529            CALL ADAMS_BASHFORTH2(
530         I                         bi, bj, k,
531         U                         gW, gwNm1,
532         I                         nHydStartAB, myIter, myThid )
533    c#endif /* ALLOW_ADAMSBASHFORTH_3 */
534    
535    C--   Dissipation term outside the Adams-Bashforth:
536            IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
537              DO j=jMin,jMax
538               DO i=iMin,iMax
539                 gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j)
540               ENDDO
541              ENDDO
542            ENDIF
543    
544    C-    end of the k loop
545          ENDDO
546    
547  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
548    
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22