/[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.16 by edhill, Fri Jun 25 13:09:40 2004 UTC revision 1.23 by jmc, Tue Dec 13 23:16:52 2005 UTC
# Line 9  CBOP Line 9  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                                              
# Line 26  C     !USES: Line 26  C     !USES:
26        IMPLICIT NONE        IMPLICIT NONE
27  C     == Global variables ==  C     == Global variables ==
28  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.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 54  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  hFacWtmp        _RS  hFacWtmp
67        _RS  hFacStmp        _RS  hFacStmp
68          _RS  hFacCtmp
69          _RS  recip_hFacCtmp
70        _RL ab15,ab05        _RL ab15,ab05
71        _RL slipSideFac        _RL slipSideFac
72        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ        _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
# Line 67  C     I,J,K - Loop counters Line 75  C     I,J,K - Loop counters
75        PARAMETER(Half=0.5D0)        PARAMETER(Half=0.5D0)
76  CEOP  CEOP
77    
 ceh3 needs an IF ( useNONHYDROSTATIC ) THEN  
   
78        iMin = 1        iMin = 1
79        iMax = sNx        iMax = sNx
80        jMin = 1        jMin = 1
81        jMax = sNy        jMax = sNy
82    
83  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
84        ab15 =  1.5 _d 0 + abeps        IF (myIter .EQ. 0) THEN
85        ab05 = -0.5 _d 0 - abeps         ab15 =  1.0 _d 0
86           ab05 =  0.0 _d 0
87          ELSE
88           ab15 =  1.5 _d 0 + abeps
89           ab05 = -0.5 _d 0 - abeps
90          ENDIF
91    
92  C     Lateral friction (no-slip, free slip, or half slip):  C     Lateral friction (no-slip, free slip, or half slip):
93        IF ( no_slip_sides ) THEN        IF ( no_slip_sides ) THEN
# Line 84  C     Lateral friction (no-slip, free sl Line 95  C     Lateral friction (no-slip, free sl
95        ELSE        ELSE
96          slipSideFac =  1. _d 0          slipSideFac =  1. _d 0
97        ENDIF        ENDIF
98  CML   half slip was used before ; keep it for now, but half slip is  CML   half slip was used before ; keep the line for now, but half slip is
99  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.
100  C        slipSideFac = 0. _d 0  C        slipSideFac = 0. _d 0
101    
102        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
# Line 93  C        slipSideFac = 0. _d 0 Line 104  C        slipSideFac = 0. _d 0
104          DO K=1,Nr          DO K=1,Nr
105           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
106            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
107             gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)             gwNm1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
108             gW(i,j,k,bi,bj) = 0.             gW(i,j,k,bi,bj) = 0.
109            ENDDO            ENDDO
110           ENDDO           ENDDO
# Line 123  C Sweep down column Line 134  C Sweep down column
134            Kp1=Nr            Kp1=Nr
135            wOverRide=0.            wOverRide=0.
136           endif           endif
137    C     horizontal bi-harmonic dissipation
138             IF (momViscosity .AND. viscA4W.NE.0. ) THEN
139    C     calculate the horizontal Laplacian of vertical flow
140    C     Zonal flux d/dx W
141              DO j=1-Oly,sNy+Oly
142               fZon(1-Olx,j)=0.
143               DO i=1-Olx+1,sNx+Olx
144                fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
145         &           *_dyG(i,j,bi,bj)
146         &           *_recip_dxC(i,j,bi,bj)
147         &           *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
148    #ifdef COSINEMETH_III
149         &           *sqcosFacU(J,bi,bj)
150    #endif
151               ENDDO
152              ENDDO
153    C     Meridional flux d/dy W
154              DO i=1-Olx,sNx+Olx
155               fMer(I,1-Oly)=0.
156              ENDDO
157              DO j=1-Oly+1,sNy+Oly
158               DO i=1-Olx,sNx+Olx
159                fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
160         &           *_dxG(i,j,bi,bj)
161         &           *_recip_dyC(i,j,bi,bj)
162         &           *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
163    #ifdef ISOTROPIC_COS_SCALING
164    #ifdef COSINEMETH_III
165         &           *sqCosFacV(j,bi,bj)
166    #endif
167    #endif
168               ENDDO
169              ENDDO
170              
171    C     del^2 W
172    C     Difference of zonal fluxes ...
173              DO j=1-Oly,sNy+Oly
174               DO i=1-Olx,sNx+Olx-1
175                del2w(i,j)=fZon(i+1,j)-fZon(i,j)
176               ENDDO
177               del2w(sNx+Olx,j)=0.
178              ENDDO
179    
180    C     ... add difference of meridional fluxes and divide by volume
181              DO j=1-Oly,sNy+Oly-1
182               DO i=1-Olx,sNx+Olx
183    C     First compute the fraction of open water for the w-control volume
184    C     at the southern face
185                hFacCtmp=max(hFacC(I,J,K-1,bi,bj)-Half,0. _d 0)
186         &           +   min(hFacC(I,J,K  ,bi,bj),Half)
187                IF (hFacCtmp .GT. 0.) THEN
188                 recip_hFacCtmp = 1./hFacCtmp
189                ELSE
190                 recip_hFacCtmp = 0. _d 0
191                ENDIF
192                del2w(i,j)=recip_rA(i,j,bi,bj)
193         &           *recip_drC(k)*recip_hFacCtmp
194         &           *(
195         &           del2w(i,j)
196         &           +( fMer(i,j+1)-fMer(i,j) )
197         &           )
198               ENDDO
199              ENDDO
200    C-- No-slip BCs impose a drag at walls...
201    CML ************************************************************
202    CML   No-slip Boundary conditions for bi-harmonic dissipation
203    CML   need to be implemented here!
204    CML ************************************************************
205             ELSE
206    C-    Initialize del2w to zero:
207              DO j=1-Oly,sNy+Oly
208               DO i=1-Olx,sNx+Olx
209                del2w(i,j) = 0. _d 0
210               ENDDO
211              ENDDO
212             ENDIF
213    
214  C Flux on Southern face  C Flux on Southern face
215           DO J=jMin,jMax+1           DO J=jMin,jMax+1
216            DO I=iMin,iMax            DO I=iMin,iMax
# Line 135  C     at the southern face Line 223  C     at the southern face
223       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))       &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
224             Flx_NS(I,J,bi,bj)=             Flx_NS(I,J,bi,bj)=
225       &     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))
226       &    -viscAh*_recip_dyC(I,J,bi,bj)         &    -viscAhW*_recip_dyC(I,J,bi,bj)  
227       &       *(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))
228       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)       &        +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
229       &         *wVel(I,J,K,bi,bj))       &         *wVel(I,J,K,bi,bj))
230         &    +viscA4W*_recip_dyC(I,J,bi,bj)*(del2w(I,J)-del2w(I,J-1))
231    #ifdef ISOTROPIC_COS_SCALING
232    #ifdef COSINEMETH_III
233         &    *sqCosFacV(j,bi,bj)
234    #else
235         &    *CosFacV(j,bi,bj)
236    #endif
237    #endif
238  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
239  C     fraction of the w control volume and at the closed fraction of the  C     fraction of the w control volume and at the closed fraction of the
240  C     the control volume. A more compact but less intelligible version  C     the control volume. A more compact but less intelligible version
# Line 154  C     First compute the fraction of open Line 250  C     First compute the fraction of open
250  C     at the western face  C     at the western face
251             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)
252       &         +    min(hFacW(I,J,K  ,bi,bj),Half)       &         +    min(hFacW(I,J,K  ,bi,bj),Half)
253                   tmp_UbarZ=Half*(             tmp_UbarZ=Half*(
254       &         _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)
255       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))       &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
256             Flx_EW(I,J,bi,bj)=             Flx_EW(I,J,bi,bj)=
257       &     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))
258       &    -viscAh*_recip_dxC(I,J,bi,bj)       &    -viscAhW*_recip_dxC(I,J,bi,bj)
259       &      *(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))
260       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)       &       +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
261       &        *wVel(I,J,K,bi,bj) )       &        *wVel(I,J,K,bi,bj) )
262         &    +viscA4W*_recip_dxC(I,J,bi,bj)*(del2w(I,J)-del2w(I-1,J))
263    #ifdef COSINEMETH_III
264         &                *sqCosFacU(j,bi,bj)
265    #else
266         &                *CosFacU(j,bi,bj)
267    #endif
268  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
269  C     fraction of the w control volume and at the closed fraction of the  C     fraction of the w control volume and at the closed fraction of the
270  C     the control volume. A more compact but less intelligible version  C     the control volume. A more compact but less intelligible version
# Line 205  caja           and an hFacUS (above V po Line 307  caja           and an hFacUS (above V po
307         ENDDO         ENDDO
308        ENDDO        ENDDO
309    
310        
311        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
312         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
313          DO K=2,Nr          DO K=2,Nr
314           DO j=jMin,jMax           DO j=jMin,jMax
315            DO i=iMin,iMax            DO i=iMin,iMax
316             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
317       &     +deltatMom*( ab15*gW(i,j,k,bi,bj)       &     +deltatMom*nh_Am2*( ab15*gW(i,j,k,bi,bj)
318       &                 +ab05*gWNM1(i,j,k,bi,bj) )       &                 +ab05*gwNm1(i,j,k,bi,bj) )
319             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.
320            ENDDO            ENDDO
321           ENDDO           ENDDO
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22