/[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.5 by adcroft, Fri Feb 2 21:04:47 2001 UTC revision 1.22 by jmc, Tue Nov 8 02:14:10 2005 UTC
# Line 1  Line 1 
1    C $Header$
2    C     !DESCRIPTION: \bv
3    C $Name$
4    
5    #include "PACKAGES_CONFIG.h"
6  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
7    
8    CBOP
9    C     !ROUTINE: CALC_GW
10    C     !INTERFACE:
11        SUBROUTINE CALC_GW(            SUBROUTINE CALC_GW(    
12       I        myThid)       I           myTime, myIter, myThid )
13  C     /==========================================================\  C     !DESCRIPTION: \bv
14  C     | S/R CALC_GW                                              |  C     *==========================================================*
15  C     \==========================================================/  C     | S/R CALC_GW                                              
16        IMPLICIT NONE  C     | o Calculate vert. velocity tendency terms ( NH, QH only )
17    C     *==========================================================*
18    C     | In NH and QH, the vertical momentum tendency must be
19    C     | calculated explicitly and included as a source term
20    C     | for a 3d pressure eqn. Calculate that term here.
21    C     | This routine is not used in HYD calculations.
22    C     *==========================================================*
23    C     \ev
24    
25    C     !USES:
26          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 "CG2D.h"  #include "DYNVARS.h"
33  #include "GW.h"  #include "NH_VARS.h"
 #include "CG3D.h"  
34    
35    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
45    
46    C     !LOCAL VARIABLES:
47  C     == Local variables ==  C     == Local variables ==
48    C     bi, bj,      :: Loop counters
49    C     iMin, iMax,
50    C     jMin, jMax
51    C     flx_NS       :: Temp. used for fVol meridional terms.
52    C     flx_EW       :: Temp. used for fVol zonal terms.
53    C     flx_Up       :: Temp. used for fVol vertical terms.
54    C     flx_Dn       :: Temp. used for fVol vertical terms.
55        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
       _RL      aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      cF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL      pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL    fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL    fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
   
56        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
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          _RS  hFacCtmp
69          _RS  recip_hFacCtmp
70        _RL ab15,ab05        _RL ab15,ab05
71          _RL slipSideFac
72        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
73    
74        _RL  Half        _RL  Half
75        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
76    CEOP
77    
78  #define I0 1        iMin = 1
79  #define In sNx        iMax = sNx
80  #define J0 1        jMin = 1
81  #define Jn sNy        jMax = sNy
82    
83  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
84        ab15=1.5+abeps        ab15 =  1.5 _d 0 + abeps
85        ab05=-0.5-abeps        ab05 = -0.5 _d 0 - abeps
86    
87    C     Lateral friction (no-slip, free slip, or half slip):
88          IF ( no_slip_sides ) THEN
89            slipSideFac = -1. _d 0
90          ELSE
91            slipSideFac =  1. _d 0
92          ENDIF
93    CML   half slip was used before ; keep the line for now, but half slip is
94    CML   not used anywhere in the code as far as I can see.
95    C        slipSideFac = 0. _d 0
96    
97        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
98         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
99          DO K=1,Nr          DO K=1,Nr
100           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
101            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
102               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
# Line 79  C For each tile Line 115  C For each tile
115         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
116    
117  C Boundaries condition at top  C Boundaries condition at top
118          DO J=J0,Jn          DO J=jMin,jMax
119           DO I=I0,In           DO I=iMin,iMax
120            Flx_Dn(I,J,bi,bj)=0.            Flx_Dn(I,J,bi,bj)=0.
121           ENDDO           ENDDO
122          ENDDO          ENDDO
# Line 93  C Sweep down column Line 129  C Sweep down column
129            Kp1=Nr            Kp1=Nr
130            wOverRide=0.            wOverRide=0.
131           endif           endif
132    C     horizontal bi-harmonic dissipation
133             IF (momViscosity .AND. viscA4W.NE.0. ) THEN
134    C     calculate the horizontal Laplacian of vertical flow
135    C     Zonal flux d/dx W
136              DO j=1-Oly,sNy+Oly
137               fZon(1-Olx,j)=0.
138               DO i=1-Olx+1,sNx+Olx
139                fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
140         &           *_dyG(i,j,bi,bj)
141         &           *_recip_dxC(i,j,bi,bj)
142         &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
143    #ifdef COSINEMETH_III
144         &           *sqcosFacU(J,bi,bj)
145    #endif
146               ENDDO
147              ENDDO
148    C     Meridional flux d/dy W
149              DO i=1-Olx,sNx+Olx
150               fMer(I,1-Oly)=0.
151              ENDDO
152              DO j=1-Oly+1,sNy+Oly
153               DO i=1-Olx,sNx+Olx
154                fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
155         &           *_dxG(i,j,bi,bj)
156         &           *_recip_dyC(i,j,bi,bj)
157         &           *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
158    #ifdef ISOTROPIC_COS_SCALING
159    #ifdef COSINEMETH_III
160         &           *sqCosFacV(j,bi,bj)
161    #endif
162    #endif
163               ENDDO
164              ENDDO
165              
166    C     del^2 W
167    C     Difference of zonal fluxes ...
168              DO j=1-Oly,sNy+Oly
169               DO i=1-Olx,sNx+Olx-1
170                del2w(i,j)=fZon(i+1,j)-fZon(i,j)
171               ENDDO
172               del2w(sNx+Olx,j)=0.
173              ENDDO
174    
175    C     ... add difference of meridional fluxes and divide by volume
176              DO j=1-Oly,sNy+Oly-1
177               DO i=1-Olx,sNx+Olx
178    C     First compute the fraction of open water for the w-control volume
179    C     at the southern face
180                hFacCtmp=max(hFacC(I,J,K-1,bi,bj)-Half,0. _d 0)
181         &           +   min(hFacC(I,J,K  ,bi,bj),Half)
182                IF (hFacCtmp .GT. 0.) THEN
183                 recip_hFacCtmp = 1./hFacCtmp
184                ELSE
185                 recip_hFacCtmp = 0. _d 0
186                ENDIF
187                del2w(i,j)=recip_rA(i,j,bi,bj)
188         &           *recip_drC(k)*recip_hFacCtmp
189         &           *(
190         &           del2w(i,j)
191         &           +( fMer(i,j+1)-fMer(i,j) )
192         &           )
193               ENDDO
194              ENDDO
195    C-- No-slip BCs impose a drag at walls...
196    CML ************************************************************
197    CML   No-slip Boundary conditions for bi-harmonic dissipation
198    CML   need to be implemented here!
199    CML ************************************************************
200             ELSE
201    C-    Initialize del2w to zero:
202              DO j=1-Oly,sNy+Oly
203               DO i=1-Olx,sNx+Olx
204                del2w(i,j) = 0. _d 0
205               ENDDO
206              ENDDO
207             ENDIF
208    
209  C Flux on Southern face  C Flux on Southern face
210           DO J=J0,Jn+1           DO J=jMin,jMax+1
211            DO I=I0,In            DO I=iMin,iMax
212    C     First compute the fraction of open water for the w-control volume
213    C     at the southern face
214               hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)
215         &         +    min(hFacS(I,J,K  ,bi,bj),Half)
216             tmp_VbarZ=Half*(             tmp_VbarZ=Half*(
217       &          _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)
218       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
219             Flx_NS(I,J,bi,bj)=             Flx_NS(I,J,bi,bj)=
220       &     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))
221       &    -viscAh*_recip_dyC(I,J,bi,bj)*(       &    -viscAhW*_recip_dyC(I,J,bi,bj)  
222       &                   wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj) )       &       *(hFacStmp*(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))
223         &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
224         &         *wVel(I,J,K,bi,bj))
225         &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))
226    #ifdef ISOTROPIC_COS_SCALING
227    #ifdef COSINEMETH_III
228         &    *sqCosFacV(j,bi,bj)
229    #else
230         &    *CosFacV(j,bi,bj)
231    #endif
232    #endif
233    C     The last term is the weighted average of the viscous stress at the open
234    C     fraction of the w control volume and at the closed fraction of the
235    C     the control volume. A more compact but less intelligible version
236    C     of the last three lines is:
237    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
238    CML     &       *wVel(I,J,K,bi,bi) + hFacStmp*wVel(I,J-1,K,bi,bj) )
239            ENDDO            ENDDO
240           ENDDO           ENDDO
241  C Flux on Western face  C Flux on Western face
242           DO J=J0,Jn           DO J=jMin,jMax
243            DO I=I0,In+1            DO I=iMin,iMax+1
244    C     First compute the fraction of open water for the w-control volume
245    C     at the western face
246               hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)
247         &         +    min(hFacW(I,J,K  ,bi,bj),Half)
248             tmp_UbarZ=Half*(             tmp_UbarZ=Half*(
249       &         _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)
250       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
251             Flx_EW(I,J,bi,bj)=             Flx_EW(I,J,bi,bj)=
252       &     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))
253       &    -viscAh*_recip_dxC(I,J,bi,bj)*(       &    -viscAhW*_recip_dxC(I,J,bi,bj)
254       &                   wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj) )       &      *(hFacWtmp*(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
255         &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
256         &        *wVel(I,J,K,bi,bj) )
257         &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))
258    #ifdef COSINEMETH_III
259         &                *sqCosFacU(j,bi,bj)
260    #else
261         &                *CosFacU(j,bi,bj)
262    #endif
263    C     The last term is the weighted average of the viscous stress at the open
264    C     fraction of the w control volume and at the closed fraction of the
265    C     the control volume. A more compact but less intelligible version
266    C     of the last three lines is:
267    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))
268    CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )
269            ENDDO            ENDDO
270           ENDDO           ENDDO
271  C Flux on Lower face  C Flux on Lower face
272           DO J=J0,Jn           DO J=jMin,jMax
273            DO I=I0,In            DO I=iMin,iMax
274             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
275             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)
276         &         +wOverRide*wVel(I,J,Kp1,bi,bj))
277             Flx_Dn(I,J,bi,bj)=             Flx_Dn(I,J,bi,bj)=
278       &     tmp_WbarZ*tmp_WbarZ       &     tmp_WbarZ*tmp_WbarZ
279       &    -viscAr*recip_drF(K)*( wVel(I,J,K,bi,bj)       &    -viscAr*recip_drF(K)
280       &                -wOverRide*wVel(I,J,Kp1,bi,bj) )       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
281            ENDDO            ENDDO
282           ENDDO           ENDDO
283  C        Divergence of fluxes  C        Divergence of fluxes
284           DO J=J0,Jn           DO J=jMin,jMax
285            DO I=I0,In            DO I=iMin,iMax
286             gW(I,J,K,bi,bj) = 0.             gW(I,J,K,bi,bj) = 0.
287       &      -(       &      -(
288       &        +_recip_dxF(I,J,bi,bj)*(       &        +_recip_dxF(I,J,bi,bj)*(
# Line 150  caja           and an hFacUS (above V po Line 302  caja           and an hFacUS (above V po
302         ENDDO         ENDDO
303        ENDDO        ENDDO
304    
305        
306        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
307         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
308          DO K=2,Nr          DO K=2,Nr
309           DO j=J0,Jn           DO j=jMin,jMax
310            DO i=I0,In            DO i=iMin,iMax
311             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
312       &     +deltatMom*( ab15*gW(i,j,k,bi,bj)       &     +deltatMom*nh_Am2*( ab15*gW(i,j,k,bi,bj)
313       &                 +ab05*gWNM1(i,j,k,bi,bj) )       &                 +ab05*gwNm1(i,j,k,bi,bj) )
314             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.
315            ENDDO            ENDDO
316           ENDDO           ENDDO
317          ENDDO          ENDDO
318         ENDDO         ENDDO
319        ENDDO        ENDDO
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO K=1,Nr  
          DO j=J0,Jn  
           DO i=I0,In  
            gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
320    
321  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
322        IF (useOBCS) THEN        IF (useOBCS) THEN
# Line 196  C     impact. This is purely for diagnos Line 337  C     impact. This is purely for diagnos
337    
338        RETURN        RETURN
339        END        END
   
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22