/[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.12 by jmc, Mon Apr 5 21:46:18 2004 UTC revision 1.26 by heimbach, Wed Jun 7 01:55:12 2006 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
3  C $Name$  C $Name$
4    
5  #include "PACKAGES_CONFIG.h"  c #include "PACKAGES_CONFIG.h"
6  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
7    
8  CBOP  CBOP
9  C     !ROUTINE: CALC_GW  C     !ROUTINE: CALC_GW
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE CALC_GW(            SUBROUTINE CALC_GW(
12       I        myThid)       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 vert. velocity tendency terms ( NH, QH only )
17  C     *==========================================================*  C     *==========================================================*
18  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH and QH, the vertical momentum tendency must be
19  C     | calculated explicitly and included as a source term  C     | calculated explicitly and included as a source term
20  C     | for a 3d pressure eqn. Calculate that term here.  C     | for a 3d pressure eqn. Calculate that term here.
21  C     | This routine is not used in HYD calculations.  C     | This routine is not used in HYD calculations.
22  C     *==========================================================*  C     *==========================================================*
23  C     \ev  C     \ev
24    
25  C     !USES:  C     !USES:
26        IMPLICIT NONE        IMPLICIT NONE
27  C     == Global variables ==  C     == Global variables ==
28  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
29  #include "EEPARAMS.h"  #include "EEPARAMS.h"
30  #include "PARAMS.h"  #include "PARAMS.h"
31  #include "GRID.h"  #include "GRID.h"
32  #include "GW.h"  #include "DYNVARS.h"
33  #include "CG3D.h"  #include "NH_VARS.h"
34    
35  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
36  C     == Routine arguments ==  C     == Routine arguments ==
37  C     myThid - Instance number for this innvocation of CALC_GW  C     myTime :: Current time in simulation
38    C     myIter :: Current iteration number in simulation
39    C     myThid :: Thread number for this instance of the routine.
40          _RL     myTime
41          INTEGER myIter
42        INTEGER myThid        INTEGER myThid
43    
44  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
# Line 55  C     flx_Dn       :: Temp. used for fVo Line 57  C     flx_Dn       :: Temp. used for fVo
57        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
59        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
60  C     I,J,K - Loop counters        _RL    fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61        INTEGER i,j,k, kP1, kUp        _RL    fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62          _RL    del2w(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63    C     i,j,k - Loop counters
64          INTEGER i,j,k, kP1
65        _RL  wOverride        _RL  wOverride
66        _RS  hFacROpen        _RS  hFacWtmp
67        _RS  hFacRClosed        _RS  hFacStmp
68        _RL ab15,ab05        _RS  hFacCtmp
69          _RS  recip_hFacCtmp
70        _RL slipSideFac        _RL slipSideFac
71        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
72    
73        _RL  Half        _RL  Half
74        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
   
 #define I0 1  
 #define In sNx  
 #define J0 1  
 #define Jn sNy  
75  CEOP  CEOP
76    
77  ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  C Catch barotropic mode
78          IF ( Nr .LT. 2 ) RETURN
79    
80  C     Adams-Bashforth timestepping weights        iMin = 1
81        ab15 =  1.5 _d 0 + abeps        iMax = sNx
82        ab05 = -0.5 _d 0 - abeps        jMin = 1
83          jMax = sNy
84    
85  C     Lateral friction (no-slip, free slip, or half slip):  C     Lateral friction (no-slip, free slip, or half slip):
86        IF ( no_slip_sides ) THEN        IF ( no_slip_sides ) THEN
87          slipSideFac = -Half          slipSideFac = -1. _d 0
88        ELSE        ELSE
89          slipSideFac =  Half          slipSideFac =  1. _d 0
90        ENDIF        ENDIF
91  C-    half slip was used before ; keep it for now.  CML   half slip was used before ; keep the line for now, but half slip is
92          slipSideFac = 0. _d 0  CML   not used anywhere in the code as far as I can see.
93    C        slipSideFac = 0. _d 0
94    
95    C For each tile
96        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
97         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
98    
99    C Initialise gW to zero
100          DO K=1,Nr          DO K=1,Nr
101           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
102            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
            gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
103             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
104            ENDDO            ENDDO
105           ENDDO           ENDDO
106          ENDDO          ENDDO
        ENDDO  
       ENDDO  
   
 C Catch barotropic mode  
       IF ( Nr .LT. 2 ) RETURN  
   
 C For each tile  
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
107    
108  C Boundaries condition at top  C Boundaries condition at top
109          DO J=J0,Jn          DO J=jMin,jMax
110           DO I=I0,In           DO I=iMin,iMax
111            Flx_Dn(I,J,bi,bj)=0.            Flx_Dn(I,J,bi,bj)=0.
112           ENDDO           ENDDO
113          ENDDO          ENDDO
# Line 123  C Sweep down column Line 120  C Sweep down column
120            Kp1=Nr            Kp1=Nr
121            wOverRide=0.            wOverRide=0.
122           endif           endif
123    C     horizontal bi-harmonic dissipation
124             IF (momViscosity .AND. viscA4W.NE.0. ) THEN
125    C     calculate the horizontal Laplacian of vertical flow
126    C     Zonal flux d/dx W
127              DO j=1-Oly,sNy+Oly
128               fZon(1-Olx,j)=0.
129               DO i=1-Olx+1,sNx+Olx
130                fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
131         &           *_dyG(i,j,bi,bj)
132         &           *_recip_dxC(i,j,bi,bj)
133         &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
134    #ifdef COSINEMETH_III
135         &           *sqcosFacU(J,bi,bj)
136    #endif
137               ENDDO
138              ENDDO
139    C     Meridional flux d/dy W
140              DO i=1-Olx,sNx+Olx
141               fMer(I,1-Oly)=0.
142              ENDDO
143              DO j=1-Oly+1,sNy+Oly
144               DO i=1-Olx,sNx+Olx
145                fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
146         &           *_dxG(i,j,bi,bj)
147         &           *_recip_dyC(i,j,bi,bj)
148         &           *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
149    #ifdef ISOTROPIC_COS_SCALING
150    #ifdef COSINEMETH_III
151         &           *sqCosFacV(j,bi,bj)
152    #endif
153    #endif
154               ENDDO
155              ENDDO
156    
157    C     del^2 W
158    C     Difference of zonal fluxes ...
159              DO j=1-Oly,sNy+Oly
160               DO i=1-Olx,sNx+Olx-1
161                del2w(i,j)=fZon(i+1,j)-fZon(i,j)
162               ENDDO
163               del2w(sNx+Olx,j)=0.
164              ENDDO
165    
166    C     ... add difference of meridional fluxes and divide by volume
167              DO j=1-Oly,sNy+Oly-1
168               DO i=1-Olx,sNx+Olx
169    C     First compute the fraction of open water for the w-control volume
170    C     at the southern face
171                hFacCtmp=max( _hFacC(I,J,K-1,bi,bj)-Half,0. _d 0 )
172         &           +   min( _hFacC(I,J,K  ,bi,bj)     ,Half )
173                IF (hFacCtmp .GT. 0.) THEN
174                 recip_hFacCtmp = 1./hFacCtmp
175                ELSE
176                 recip_hFacCtmp = 0. _d 0
177                ENDIF
178                del2w(i,j)=recip_rA(i,j,bi,bj)
179         &           *recip_drC(k)*recip_hFacCtmp
180         &           *(
181         &           del2w(i,j)
182         &           +( fMer(i,j+1)-fMer(i,j) )
183         &           )
184               ENDDO
185              ENDDO
186    C-- No-slip BCs impose a drag at walls...
187    CML ************************************************************
188    CML   No-slip Boundary conditions for bi-harmonic dissipation
189    CML   need to be implemented here!
190    CML ************************************************************
191             ELSE
192    C-    Initialize del2w to zero:
193              DO j=1-Oly,sNy+Oly
194               DO i=1-Olx,sNx+Olx
195                del2w(i,j) = 0. _d 0
196               ENDDO
197              ENDDO
198             ENDIF
199    
200  C Flux on Southern face  C Flux on Southern face
201           DO J=J0,Jn+1           DO J=jMin,jMax+1
202            DO I=I0,In            DO I=iMin,iMax
203    C     First compute the fraction of open water for the w-control volume
204    C     at the southern face
205               hFacStmp=max(_hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)
206         &         +    min(_hFacS(I,J,K  ,bi,bj),Half)
207             tmp_VbarZ=Half*(             tmp_VbarZ=Half*(
208       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)       &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)
209       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
210             Flx_NS(I,J,bi,bj)=             Flx_NS(I,J,bi,bj)=
211       &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))       &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))
212       &    -viscAh*_recip_dyC(I,J,bi,bj)         &    -viscAhW*_recip_dyC(I,J,bi,bj)
213       &      *(1. _d 0 + slipSideFac*       &       *(hFacStmp*(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))
214       &         (maskS(I,J,K-1,bi,bj)+maskS(I,J,K,bi,bj)-2. _d 0))       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
215       &                   *(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))       &         *wVel(I,J,K,bi,bj))
216         &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))
217    #ifdef ISOTROPIC_COS_SCALING
218    #ifdef COSINEMETH_III
219         &    *sqCosFacV(j,bi,bj)
220    #else
221         &    *CosFacV(j,bi,bj)
222    #endif
223    #endif
224    C     The last term is the weighted average of the viscous stress at the open
225    C     fraction of the w control volume and at the closed fraction of the
226    C     the control volume. A more compact but less intelligible version
227    C     of the last three lines is:
228    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
229    CML     &       *wVel(I,J,K,bi,bi) + hFacStmp*wVel(I,J-1,K,bi,bj) )
230            ENDDO            ENDDO
231           ENDDO           ENDDO
232  C Flux on Western face  C Flux on Western face
233           DO J=J0,Jn           DO J=jMin,jMax
234            DO I=I0,In+1            DO I=iMin,iMax+1
235    C     First compute the fraction of open water for the w-control volume
236    C     at the western face
237               hFacWtmp=max(_hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)
238         &         +    min(_hFacW(I,J,K  ,bi,bj),Half)
239             tmp_UbarZ=Half*(             tmp_UbarZ=Half*(
240       &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)       &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)
241       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
242             Flx_EW(I,J,bi,bj)=             Flx_EW(I,J,bi,bj)=
243       &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))       &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))
244       &    -viscAh*_recip_dxC(I,J,bi,bj)       &    -viscAhW*_recip_dxC(I,J,bi,bj)
245       &      *(1. _d 0 + slipSideFac*       &      *(hFacWtmp*(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
246       &         (maskW(I,J,K-1,bi,bj)+maskW(I,J,K,bi,bj)-2. _d 0))       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
247       &                   *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))       &        *wVel(I,J,K,bi,bj) )
248         &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))
249    #ifdef COSINEMETH_III
250         &                *sqCosFacU(j,bi,bj)
251    #else
252         &                *CosFacU(j,bi,bj)
253    #endif
254    C     The last term is the weighted average of the viscous stress at the open
255    C     fraction of the w control volume and at the closed fraction of the
256    C     the control volume. A more compact but less intelligible version
257    C     of the last three lines is:
258    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))
259    CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )
260            ENDDO            ENDDO
261           ENDDO           ENDDO
262  C Flux on Lower face  C Flux on Lower face
263           DO J=J0,Jn           DO J=jMin,jMax
264            DO I=I0,In            DO I=iMin,iMax
265             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
266             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)
267         &         +wOverRide*wVel(I,J,Kp1,bi,bj))
268             Flx_Dn(I,J,bi,bj)=             Flx_Dn(I,J,bi,bj)=
269       &     tmp_WbarZ*tmp_WbarZ       &     tmp_WbarZ*tmp_WbarZ
270       &    -viscAr*recip_drF(K)       &    -viscAr*recip_drF(K)
271       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
272            ENDDO            ENDDO
273           ENDDO           ENDDO
274  C        Divergence of fluxes  C        Divergence of fluxes
275           DO J=J0,Jn           DO J=jMin,jMax
276            DO I=I0,In            DO I=iMin,iMax
277             gW(I,J,K,bi,bj) = 0.             gW(I,J,K,bi,bj) = 0.
278       &      -(       &      -(
279       &        +_recip_dxF(I,J,bi,bj)*(       &        +_recip_dxF(I,J,bi,bj)*(
# Line 175  C        Divergence of fluxes Line 284  C        Divergence of fluxes
284       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )
285       &       )       &       )
286  caja    * recip_hFacU(I,J,K,bi,bj)  caja    * recip_hFacU(I,J,K,bi,bj)
287  caja   NOTE:  This should be included    caja   NOTE:  This should be included
288  caja   but we need an hFacUW (above U points)  caja   but we need an hFacUW (above U points)
289  caja           and an hFacUS (above V points) too...  caja           and an hFacUS (above V points) too...
290            ENDDO            ENDDO
291           ENDDO           ENDDO
         ENDDO  
        ENDDO  
       ENDDO  
292    
293        C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO K=2,Nr  
          DO j=J0,Jn  
           DO i=I0,In  
            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.  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
294    
295  #ifdef ALLOW_OBCS  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
296        IF (useOBCS) THEN  C     and save gW_[n] into gwNm1 for the next time step.
297  C--   This call is aesthetic: it makes the W field  c#ifdef ALLOW_ADAMSBASHFORTH_3
298  C     consistent with the OBs but this has no algorithmic  c        CALL ADAMS_BASHFORTH3(
299  C     impact. This is purely for diagnostic purposes.  c    I                          bi, bj, k,
300         DO bj=myByLo(myThid),myByHi(myThid)  c    U                          gW, gwNm,
301          DO bi=myBxLo(myThid),myBxHi(myThid)  c    I                          momStartAB, myIter, myThid )
302           DO K=1,Nr  c#else /* ALLOW_ADAMSBASHFORTH_3 */
303            CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )           CALL ADAMS_BASHFORTH2(
304           ENDDO       I                          bi, bj, k,
305         U                          gW, gwNm1,
306         I                          myIter, myThid )
307    c#endif /* ALLOW_ADAMSBASHFORTH_3 */
308    
309    C-    end of the k loop
310          ENDDO          ENDDO
311    
312    C-    end of bi,bj loops
313         ENDDO         ENDDO
314        ENDIF        ENDDO
 #endif /* ALLOW_OBCS */  
315    
316  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
317    
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22