/[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.24 by jmc, Thu Dec 15 02:06:06 2005 UTC revision 1.27 by baylor, Tue Jun 20 20:57:37 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           myTime, myIter, myThid )       I         KappaRU, KappaRV,
13         I         myTime, myIter, myThid )
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | S/R CALC_GW                                                C     | S/R CALC_GW
17  C     | o Calculate vert. velocity tendency terms ( NH, QH only )  C     | o Calculate vert. velocity tendency terms ( NH, QH only )
18  C     *==========================================================*  C     *==========================================================*
19  C     | In NH and QH, the vertical momentum tendency must be  C     | In NH and QH, 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
# Line 34  C     == Global variables == Line 35  C     == Global variables ==
35    
36  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38    C     KappaRU:: vertical viscosity at U points
39    C     KappaRV:: vertical viscosity at V points
40  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
41  C     myIter :: Current iteration number in simulation  C     myIter :: Current iteration number in simulation
42  C     myThid :: Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
43          _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
44          _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
45        _RL     myTime        _RL     myTime
46        INTEGER myIter        INTEGER myIter
47        INTEGER myThid        INTEGER myThid
# Line 67  C     i,j,k - Loop counters Line 72  C     i,j,k - Loop counters
72        _RS  hFacStmp        _RS  hFacStmp
73        _RS  hFacCtmp        _RS  hFacCtmp
74        _RS  recip_hFacCtmp        _RS  recip_hFacCtmp
       _RL ab15,ab05  
75        _RL slipSideFac        _RL slipSideFac
76        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
77          _RL vischere
78          _RL visc4here
79        _RL  Half        _RL  Half
80        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
81  CEOP  CEOP
82    
83    C Catch barotropic mode
84          IF ( Nr .LT. 2 ) RETURN
85    
86        iMin = 1        iMin = 1
87        iMax = sNx        iMax = sNx
88        jMin = 1        jMin = 1
89        jMax = sNy        jMax = sNy
90    
 C     Adams-Bashforth timestepping weights  
       IF (myIter .EQ. 0) THEN  
        ab15 =  1.0 _d 0  
        ab05 =  0.0 _d 0  
       ELSE  
        ab15 =  1.5 _d 0 + abeps  
        ab05 = -0.5 _d 0 - abeps  
       ENDIF  
   
91  C     Lateral friction (no-slip, free slip, or half slip):  C     Lateral friction (no-slip, free slip, or half slip):
92        IF ( no_slip_sides ) THEN        IF ( no_slip_sides ) THEN
93          slipSideFac = -1. _d 0          slipSideFac = -1. _d 0
94        ELSE        ELSE
95          slipSideFac =  1. _d 0          slipSideFac =  1. _d 0
96        ENDIF        ENDIF
97  CML   half slip was used before ; keep the line for now, but half slip is  CML   half slip was used before ; keep the line for now, but half slip is
98  CML   not used anywhere in the code as far as I can see.  CML   not used anywhere in the code as far as I can see.
99  C        slipSideFac = 0. _d 0  C        slipSideFac = 0. _d 0
100    
101    C For each tile
102        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
103         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
104    
105    C Initialise gW to zero
106          DO K=1,Nr          DO K=1,Nr
107           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
108            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
# Line 108  C        slipSideFac = 0. _d 0 Line 110  C        slipSideFac = 0. _d 0
110            ENDDO            ENDDO
111           ENDDO           ENDDO
112          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)  
113    
114  C Boundaries condition at top  C Boundaries condition at top
115          DO J=jMin,jMax          DO J=jMin,jMax
# Line 166  C     Meridional flux d/dy W Line 159  C     Meridional flux d/dy W
159  #endif  #endif
160             ENDDO             ENDDO
161            ENDDO            ENDDO
162              
163  C     del^2 W  C     del^2 W
164  C     Difference of zonal fluxes ...  C     Difference of zonal fluxes ...
165            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
# Line 181  C     ... add difference of meridional f Line 174  C     ... add difference of meridional f
174             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
175  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
176  C     at the southern face  C     at the southern face
177              hFacCtmp=max(hFacC(I,J,K-1,bi,bj)-Half,0. _d 0)              hFacCtmp=max( _hFacC(I,J,K-1,bi,bj)-Half,0. _d 0 )
178       &           +   min(hFacC(I,J,K  ,bi,bj),Half)       &           +   min( _hFacC(I,J,K  ,bi,bj)     ,Half )
179              IF (hFacCtmp .GT. 0.) THEN              IF (hFacCtmp .GT. 0.) THEN
180               recip_hFacCtmp = 1./hFacCtmp               recip_hFacCtmp = 1./hFacCtmp
181              ELSE              ELSE
# Line 215  C Flux on Southern face Line 208  C Flux on Southern face
208            DO I=iMin,iMax            DO I=iMin,iMax
209  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
210  C     at the southern face  C     at the southern face
211             hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)             hFacStmp=max(_hFacS(I,J,K-1,bi,bj)-Half,0. _d 0)
212       &         +    min(hFacS(I,J,K  ,bi,bj),Half)       &         +    min(_hFacS(I,J,K  ,bi,bj),Half)
213             tmp_VbarZ=Half*(             tmp_VbarZ=Half*(
214       &          _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)
215       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
216             Flx_NS(I,J,bi,bj)=             Flx_NS(I,J,bi,bj)=
217       &     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))
218       &    -viscAhW*_recip_dyC(I,J,bi,bj)         &    -viscAh_W(I,J,K,bi,bj)*_recip_dyC(I,J,bi,bj)
219       &       *(hFacStmp*(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))
220       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
221       &         *wVel(I,J,K,bi,bj))       &         *wVel(I,J,K,bi,bj))
222       &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))       &    +viscA4_W(I,J,K,bi,bj)
223         &      *_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))
224  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
225  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
226       &    *sqCosFacV(j,bi,bj)       &    *sqCosFacV(j,bi,bj)
# Line 247  C Flux on Western face Line 241  C Flux on Western face
241            DO I=iMin,iMax+1            DO I=iMin,iMax+1
242  C     First compute the fraction of open water for the w-control volume  C     First compute the fraction of open water for the w-control volume
243  C     at the western face  C     at the western face
244             hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)             hFacWtmp=max(_hFacW(I,J,K-1,bi,bj)-Half,0. _d 0)
245       &         +    min(hFacW(I,J,K  ,bi,bj),Half)       &         +    min(_hFacW(I,J,K  ,bi,bj),Half)
246             tmp_UbarZ=Half*(             tmp_UbarZ=Half*(
247       &         _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)
248       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
249             Flx_EW(I,J,bi,bj)=             Flx_EW(I,J,bi,bj)=
250       &     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))
251       &    -viscAhW*_recip_dxC(I,J,bi,bj)       &    -viscAh_W(I,J,K,bi,bj)*_recip_dxC(I,J,bi,bj)
252       &      *(hFacWtmp*(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))
253       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
254       &        *wVel(I,J,K,bi,bj) )       &        *wVel(I,J,K,bi,bj) )
255       &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))       &    +viscA4_W(I,J,K,bi,bj)
256         &      *_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))
257  #ifdef COSINEMETH_III  #ifdef COSINEMETH_III
258       &                *sqCosFacU(j,bi,bj)       &                *sqCosFacU(j,bi,bj)
259  #else  #else
260       &                *CosFacU(j,bi,bj)       &                *CosFacU(j,bi,bj)
261  #endif  #endif
262  C     The last term is the weighted average of the viscous stress at the open  C     The last term is the weighted average of the viscous stress at the open
# Line 275  CML     &       *wVel(I,J,K,bi,bi) + hFa Line 270  CML     &       *wVel(I,J,K,bi,bi) + hFa
270  C Flux on Lower face  C Flux on Lower face
271           DO J=jMin,jMax           DO J=jMin,jMax
272            DO I=iMin,iMax            DO I=iMin,iMax
273    C Interpolate vert viscosity to W points
274               vischere=0.125*(kappaRU(I,J,K)  +kappaRU(I+1,J,K)
275         &                    +kappaRU(I,J,Kp1)+kappaRU(I+1,J,Kp1)
276         &                    +kappaRV(I,J,K)  +kappaRV(I,J+1,K)
277         &                    +kappaRV(I,J,Kp1)+kappaRV(I,J+1,Kp1))
278             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
279             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)
280       &         +wOverRide*wVel(I,J,Kp1,bi,bj))       &         +wOverRide*wVel(I,J,Kp1,bi,bj))
281             Flx_Dn(I,J,bi,bj)=             Flx_Dn(I,J,bi,bj)=
282       &     tmp_WbarZ*tmp_WbarZ       &     tmp_WbarZ*tmp_WbarZ
283       &    -viscAr*recip_drF(K)       &    -vischere*recip_drF(K)
284       &       *( 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) )
285            ENDDO            ENDDO
286           ENDDO           ENDDO
# Line 297  C        Divergence of fluxes Line 297  C        Divergence of fluxes
297       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )       &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )
298       &       )       &       )
299  caja    * recip_hFacU(I,J,K,bi,bj)  caja    * recip_hFacU(I,J,K,bi,bj)
300  caja   NOTE:  This should be included    caja   NOTE:  This should be included
301  caja   but we need an hFacUW (above U points)  caja   but we need an hFacUW (above U points)
302  caja           and an hFacUS (above V points) too...  caja           and an hFacUS (above V points) too...
303            ENDDO            ENDDO
304           ENDDO           ENDDO
         ENDDO  
        ENDDO  
       ENDDO  
305    
306    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
307    
308        DO bj=myByLo(myThid),myByHi(myThid)  C-    Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
309         DO bi=myBxLo(myThid),myBxHi(myThid)  C     and save gW_[n] into gwNm1 for the next time step.
310          DO K=2,Nr  c#ifdef ALLOW_ADAMSBASHFORTH_3
311           DO j=jMin,jMax  c        CALL ADAMS_BASHFORTH3(
312            DO i=iMin,iMax  c    I                          bi, bj, k,
313             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)  c    U                          gW, gwNm,
314       &     +deltatMom*nh_Am2*( ab15*gW(i,j,k,bi,bj)  c    I                          momStartAB, myIter, myThid )
315       &                 +ab05*gwNm1(i,j,k,bi,bj) )  c#else /* ALLOW_ADAMSBASHFORTH_3 */
316             IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.           CALL ADAMS_BASHFORTH2(
317             gwNm1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)       I                          bi, bj, k,
318            ENDDO       U                          gW, gwNm1,
319           ENDDO       I                          myIter, myThid )
320          ENDDO  c#endif /* ALLOW_ADAMSBASHFORTH_3 */
        ENDDO  
       ENDDO  
321    
322  #ifdef ALLOW_OBCS  C-    end of the k loop
       IF (useOBCS) THEN  
 C--   This call is aesthetic: it makes the W field  
 C     consistent with the OBs but this has no algorithmic  
 C     impact. This is purely for diagnostic purposes.  
        DO bj=myByLo(myThid),myByHi(myThid)  
         DO bi=myBxLo(myThid),myBxHi(myThid)  
          DO K=1,Nr  
           CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )  
          ENDDO  
323          ENDDO          ENDDO
324    
325    C-    end of bi,bj loops
326         ENDDO         ENDDO
327        ENDIF        ENDDO
 #endif /* ALLOW_OBCS */  
328    
329  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
330    
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22