/[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.15 by mlosch, Tue May 25 07:48:22 2004 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        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"
29  #include "DYNVARS.h"  #include "DYNVARS.h"
 #include "FFIELDS.h"  
30  #include "EEPARAMS.h"  #include "EEPARAMS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #include "GRID.h"  #include "GRID.h"
 #include "CG2D.h"  
33  #include "GW.h"  #include "GW.h"
34  #include "CG3D.h"  #include "CG3D.h"
35    
36    C     !INPUT/OUTPUT PARAMETERS:
37  C     == Routine arguments ==  C     == Routine arguments ==
38  C     myThid - Instance number for this innvocation of CALC_GW  C     myThid - Instance number for this innvocation of CALC_GW
39        INTEGER myThid        INTEGER myThid
40    
41  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
42    
43    C     !LOCAL VARIABLES:
44  C     == Local variables ==  C     == Local variables ==
45    C     bi, bj,      :: Loop counters
46    C     iMin, iMax,
47    C     jMin, jMax
48    C     flx_NS       :: Temp. used for fVol meridional terms.
49    C     flx_EW       :: Temp. used for fVol zonal terms.
50    C     flx_Up       :: Temp. used for fVol vertical terms.
51    C     flx_Dn       :: Temp. used for fVol vertical terms.
52        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)  
   
53        _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)
54        _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)
55        _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)
# Line 42  C     == Local variables == Line 57  C     == Local variables ==
57  C     I,J,K - Loop counters  C     I,J,K - Loop counters
58        INTEGER i,j,k, kP1, kUp        INTEGER i,j,k, kP1, kUp
59        _RL  wOverride        _RL  wOverride
60        _RS  hFacROpen        _RS  hFacWtmp
61        _RS  hFacRClosed        _RS  hFacStmp
62        _RL ab15,ab05        _RL ab15,ab05
63          _RL slipSideFac
64        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
65    
66        _RL  Half        _RL  Half
67        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
68    CEOP
69    
70    ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
71    
72  #define I0 1        iMin = 1
73  #define In sNx        iMax = sNx
74  #define J0 1        jMin = 1
75  #define Jn sNy        jMax = sNy
76    
77  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
78        ab15=1.5+abeps        ab15 =  1.5 _d 0 + abeps
79        ab05=-0.5-abeps        ab05 = -0.5 _d 0 - abeps
80    
81    C     Lateral friction (no-slip, free slip, or half slip):
82          IF ( no_slip_sides ) THEN
83            slipSideFac = -1. _d 0
84          ELSE
85            slipSideFac =  1. _d 0
86          ENDIF
87    CML   half slip was used before ; keep it for now, but half slip is
88    CML   not used anywhere in the code as far as I can see
89    C        slipSideFac = 0. _d 0
90    
91        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
92         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
93          DO K=1,Nr          DO K=1,Nr
94           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
95            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
96               gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
97             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
98            ENDDO            ENDDO
99           ENDDO           ENDDO
# Line 79  C For each tile Line 109  C For each tile
109         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
110    
111  C Boundaries condition at top  C Boundaries condition at top
112          DO J=J0,Jn          DO J=jMin,jMax
113           DO I=I0,In           DO I=iMin,iMax
114            Flx_Dn(I,J,bi,bj)=0.            Flx_Dn(I,J,bi,bj)=0.
115           ENDDO           ENDDO
116          ENDDO          ENDDO
# Line 94  C Sweep down column Line 124  C Sweep down column
124            wOverRide=0.            wOverRide=0.
125           endif           endif
126  C Flux on Southern face  C Flux on Southern face
127           DO J=J0,Jn+1           DO J=jMin,jMax+1
128            DO I=I0,In            DO I=iMin,iMax
129    C     First compute the fraction of open water for the w-control volume
130    C     at the southern face
131               hFacStmp=max(hFacS(I,J,K-1,bi,bj)-Half,0)
132         &         +    min(hFacS(I,J,K  ,bi,bj),Half)
133             tmp_VbarZ=Half*(             tmp_VbarZ=Half*(
134       &          _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)
135       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
136             Flx_NS(I,J,bi,bj)=             Flx_NS(I,J,bi,bj)=
137       &     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))
138       &    -viscAh*_recip_dyC(I,J,bi,bj)*(       &    -viscAh*_recip_dyC(I,J,bi,bj)  
139       &                   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))
140         &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
141         &         *wVel(I,J,K,bi,bj))
142    C     The last term is the weighted average of the viscous stress at the open
143    C     fraction of the w control volume and at the closed fraction of the
144    C     the control volume. A more compact but less intelligible version
145    C     of the last three lines is:
146    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
147    CML     &       *wVel(I,J,K,bi,bi) + hFacStmp*wVel(I,J-1,K,bi,bj) )
148            ENDDO            ENDDO
149           ENDDO           ENDDO
150  C Flux on Western face  C Flux on Western face
151           DO J=J0,Jn           DO J=jMin,jMax
152            DO I=I0,In+1            DO I=iMin,iMax+1
153             tmp_UbarZ=Half*(  C     First compute the fraction of open water for the w-control volume
154    C     at the western face
155               hFacWtmp=max(hFacW(I,J,K-1,bi,bj)-Half,0)
156         &         +    min(hFacW(I,J,K  ,bi,bj),Half)
157                     tmp_UbarZ=Half*(
158       &         _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)
159       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
160             Flx_EW(I,J,bi,bj)=             Flx_EW(I,J,bi,bj)=
161       &     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))
162       &    -viscAh*_recip_dxC(I,J,bi,bj)*(       &    -viscAh*_recip_dxC(I,J,bi,bj)
163       &                   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))
164         &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
165         &        *wVel(I,J,K,bi,bj) )
166    C     The last term is the weighted average of the viscous stress at the open
167    C     fraction of the w control volume and at the closed fraction of the
168    C     the control volume. A more compact but less intelligible version
169    C     of the last three lines is:
170    CML     &       *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))
171    CML     &       *wVel(I,J,K,bi,bi) + hFacWtmp*wVel(I-1,J,K,bi,bj) )
172            ENDDO            ENDDO
173           ENDDO           ENDDO
174  C Flux on Lower face  C Flux on Lower face
175           DO J=J0,Jn           DO J=jMin,jMax
176            DO I=I0,In            DO I=iMin,iMax
177             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)             Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
178             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))             tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)
179         &         +wOverRide*wVel(I,J,Kp1,bi,bj))
180             Flx_Dn(I,J,bi,bj)=             Flx_Dn(I,J,bi,bj)=
181       &     tmp_WbarZ*tmp_WbarZ       &     tmp_WbarZ*tmp_WbarZ
182       &    -viscAr*recip_drF(K)*( wVel(I,J,K,bi,bj)       &    -viscAr*recip_drF(K)
183       &                -wOverRide*wVel(I,J,Kp1,bi,bj) )       &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
184            ENDDO            ENDDO
185           ENDDO           ENDDO
186  C        Divergence of fluxes  C        Divergence of fluxes
187           DO J=J0,Jn           DO J=jMin,jMax
188            DO I=I0,In            DO I=iMin,iMax
189             gW(I,J,K,bi,bj) = 0.             gW(I,J,K,bi,bj) = 0.
190       &      -(       &      -(
191       &        +_recip_dxF(I,J,bi,bj)*(       &        +_recip_dxF(I,J,bi,bj)*(
# Line 154  caja           and an hFacUS (above V po Line 209  caja           and an hFacUS (above V po
209        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
210         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
211          DO K=2,Nr          DO K=2,Nr
212           DO j=J0,Jn           DO j=jMin,jMax
213            DO i=I0,In            DO i=iMin,iMax
214             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
215       &     +deltatMom*( ab15*gW(i,j,k,bi,bj)       &     +deltatMom*( ab15*gW(i,j,k,bi,bj)
216       &                 +ab05*gWNM1(i,j,k,bi,bj) )       &                 +ab05*gWNM1(i,j,k,bi,bj) )
# Line 165  caja           and an hFacUS (above V po Line 220  caja           and an hFacUS (above V po
220          ENDDO          ENDDO
221         ENDDO         ENDDO
222        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  
223    
224  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
225        IF (useOBCS) THEN        IF (useOBCS) THEN
# Line 196  C     impact. This is purely for diagnos Line 240  C     impact. This is purely for diagnos
240    
241        RETURN        RETURN
242        END        END
   
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22