/[MITgcm]/MITgcm/model/src/calc_gw.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_gw.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.29 - (hide annotations) (download)
Mon Jul 10 15:17:19 2006 UTC (17 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.28: +3 -1 lines
need to fix the side-drag implementation; for now, just impose free-slip

1 jmc 1.29 C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.28 2006/07/07 20:09:37 jmc Exp $
2 cnh 1.9 C !DESCRIPTION: \bv
3 jmc 1.7 C $Name: $
4 edhill 1.10
5 adcroft 1.1 #include "CPP_OPTIONS.h"
6    
7 cnh 1.9 CBOP
8     C !ROUTINE: CALC_GW
9     C !INTERFACE:
10 jmc 1.25 SUBROUTINE CALC_GW(
11 jmc 1.28 I bi, bj, KappaRU, KappaRV,
12     I myTime, myIter, myThid )
13 cnh 1.9 C !DESCRIPTION: \bv
14     C *==========================================================*
15 jmc 1.25 C | S/R CALC_GW
16     C | o Calculate vert. velocity tendency terms ( NH, QH only )
17 cnh 1.9 C *==========================================================*
18     C | In NH and QH, the vertical momentum tendency must be
19 jmc 1.25 C | calculated explicitly and included as a source term
20 cnh 1.9 C | for a 3d pressure eqn. Calculate that term here.
21     C | This routine is not used in HYD calculations.
22     C *==========================================================*
23 jmc 1.25 C \ev
24 cnh 1.9
25     C !USES:
26 adcroft 1.1 IMPLICIT NONE
27     C == Global variables ==
28     #include "SIZE.h"
29     #include "EEPARAMS.h"
30     #include "PARAMS.h"
31     #include "GRID.h"
32 jmc 1.22 #include "DYNVARS.h"
33     #include "NH_VARS.h"
34 adcroft 1.1
35 cnh 1.9 C !INPUT/OUTPUT PARAMETERS:
36 adcroft 1.1 C == Routine arguments ==
37 jmc 1.28 C bi,bj :: current tile indices
38     C KappaRU :: vertical viscosity at U points
39     C KappaRV :: vertical viscosity at V points
40     C myTime :: Current time in simulation
41     C myIter :: Current iteration number in simulation
42     C myThid :: Thread number for this instance of the routine.
43     INTEGER bi,bj
44 baylor 1.27 _RL KappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
45     _RL KappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
46 jmc 1.20 _RL myTime
47     INTEGER myIter
48 adcroft 1.1 INTEGER myThid
49    
50 adcroft 1.3 #ifdef ALLOW_NONHYDROSTATIC
51    
52 cnh 1.9 C !LOCAL VARIABLES:
53 adcroft 1.1 C == Local variables ==
54 cnh 1.9 C iMin, iMax,
55     C jMin, jMax
56     C flx_NS :: Temp. used for fVol meridional terms.
57     C flx_EW :: Temp. used for fVol zonal terms.
58     C flx_Up :: Temp. used for fVol vertical terms.
59     C flx_Dn :: Temp. used for fVol vertical terms.
60 jmc 1.28 INTEGER iMin,iMax,jMin,jMax
61     _RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62     _RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63     _RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64     _RL flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 jmc 1.21 C i,j,k - Loop counters
69 jmc 1.28 INTEGER i,j,k, kp1
70 adcroft 1.1 _RL wOverride
71 mlosch 1.15 _RS hFacWtmp
72     _RS hFacStmp
73 mlosch 1.18 _RS hFacCtmp
74     _RS recip_hFacCtmp
75 jmc 1.28 _RL slipSideFac
76     _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
77     _RL viscLoc
78 adcroft 1.1 _RL Half
79     PARAMETER(Half=0.5D0)
80 cnh 1.9 CEOP
81 edhill 1.10
82 jmc 1.25 C Catch barotropic mode
83     IF ( Nr .LT. 2 ) RETURN
84    
85 mlosch 1.14 iMin = 1
86     iMax = sNx
87     jMin = 1
88     jMax = sNy
89    
90 jmc 1.11 C Lateral friction (no-slip, free slip, or half slip):
91     IF ( no_slip_sides ) THEN
92 mlosch 1.15 slipSideFac = -1. _d 0
93 jmc 1.11 ELSE
94 jmc 1.25 slipSideFac = 1. _d 0
95 jmc 1.11 ENDIF
96 mlosch 1.18 CML half slip was used before ; keep the line for now, but half slip is
97     CML not used anywhere in the code as far as I can see.
98 mlosch 1.14 C slipSideFac = 0. _d 0
99 jmc 1.29 C- need to fix the side-drag implementation; for now, always impose free-slip
100     slipSideFac = 1. _d 0
101 jmc 1.11
102 jmc 1.28 C-- Initialise gW to zero
103     DO k=1,Nr
104     DO j=1-OLy,sNy+OLy
105     DO i=1-OLx,sNx+OLx
106 adcroft 1.1 gW(i,j,k,bi,bj) = 0.
107     ENDDO
108     ENDDO
109 jmc 1.28 ENDDO
110 adcroft 1.1
111 jmc 1.28 C-- Boundaries condition at top
112     DO j=jMin,jMax
113     DO i=iMin,iMax
114     flx_Up(i,j)=0.
115     ENDDO
116     ENDDO
117 adcroft 1.1
118 jmc 1.28 C--- Sweep down column
119     DO k=2,Nr
120     kp1=k+1
121     wOverRide=1.
122     IF (k.EQ.Nr) THEN
123     kp1=Nr
124 adcroft 1.1 wOverRide=0.
125 jmc 1.28 ENDIF
126     C-- horizontal bi-harmonic dissipation
127     IF (momViscosity .AND. viscA4W.NE.0. ) THEN
128     C- calculate the horizontal Laplacian of vertical flow
129 mlosch 1.18 C Zonal flux d/dx W
130     DO j=1-Oly,sNy+Oly
131     fZon(1-Olx,j)=0.
132     DO i=1-Olx+1,sNx+Olx
133 jmc 1.28 C- Problem here: drF(k)*_hFacC & fZon are not at the same Horiz.&Vert. location
134 mlosch 1.18 fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
135     & *_dyG(i,j,bi,bj)
136     & *_recip_dxC(i,j,bi,bj)
137     & *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
138     #ifdef COSINEMETH_III
139 jmc 1.28 & *sqcosFacU(j,bi,bj)
140 mlosch 1.18 #endif
141     ENDDO
142 jmc 1.28 ENDDO
143 mlosch 1.18 C Meridional flux d/dy W
144     DO i=1-Olx,sNx+Olx
145 jmc 1.28 fMer(i,1-Oly)=0.
146 mlosch 1.18 ENDDO
147     DO j=1-Oly+1,sNy+Oly
148     DO i=1-Olx,sNx+Olx
149 jmc 1.28 C- Problem here: drF(k)*_hFacC & fMer are not at the same Horiz.&Vert. location
150 mlosch 1.18 fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
151     & *_dxG(i,j,bi,bj)
152     & *_recip_dyC(i,j,bi,bj)
153     & *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
154     #ifdef ISOTROPIC_COS_SCALING
155     #ifdef COSINEMETH_III
156     & *sqCosFacV(j,bi,bj)
157     #endif
158     #endif
159     ENDDO
160     ENDDO
161 jmc 1.25
162 mlosch 1.18 C del^2 W
163     C Difference of zonal fluxes ...
164     DO j=1-Oly,sNy+Oly
165     DO i=1-Olx,sNx+Olx-1
166     del2w(i,j)=fZon(i+1,j)-fZon(i,j)
167     ENDDO
168     del2w(sNx+Olx,j)=0.
169     ENDDO
170    
171     C ... add difference of meridional fluxes and divide by volume
172     DO j=1-Oly,sNy+Oly-1
173     DO i=1-Olx,sNx+Olx
174     C First compute the fraction of open water for the w-control volume
175     C at the southern face
176 jmc 1.28 hFacCtmp=max( _hFacC(i,j,k-1,bi,bj)-Half,0. _d 0 )
177     & + min( _hFacC(i,j,k ,bi,bj) ,Half )
178 mlosch 1.18 IF (hFacCtmp .GT. 0.) THEN
179     recip_hFacCtmp = 1./hFacCtmp
180     ELSE
181     recip_hFacCtmp = 0. _d 0
182     ENDIF
183     del2w(i,j)=recip_rA(i,j,bi,bj)
184     & *recip_drC(k)*recip_hFacCtmp
185     & *(
186     & del2w(i,j)
187     & +( fMer(i,j+1)-fMer(i,j) )
188     & )
189     ENDDO
190     ENDDO
191     C-- No-slip BCs impose a drag at walls...
192     CML ************************************************************
193     CML No-slip Boundary conditions for bi-harmonic dissipation
194     CML need to be implemented here!
195     CML ************************************************************
196 jmc 1.28 ELSE
197 jmc 1.21 C- Initialize del2w to zero:
198     DO j=1-Oly,sNy+Oly
199     DO i=1-Olx,sNx+Olx
200     del2w(i,j) = 0. _d 0
201     ENDDO
202     ENDDO
203 jmc 1.28 ENDIF
204 mlosch 1.18
205 adcroft 1.1 C Flux on Southern face
206 jmc 1.28 DO j=jMin,jMax+1
207     DO i=iMin,iMax
208 mlosch 1.15 C First compute the fraction of open water for the w-control volume
209     C at the southern face
210 jmc 1.28 hFacStmp=max(_hFacS(i,j,k-1,bi,bj)-Half,0. _d 0)
211     & + min(_hFacS(i,j,k ,bi,bj),Half)
212 adcroft 1.1 tmp_VbarZ=Half*(
213 jmc 1.28 & _hFacS(i,j,k-1,bi,bj)*vVel( i ,j,k-1,bi,bj)
214     & +_hFacS(i,j, k ,bi,bj)*vVel( i ,j, k ,bi,bj))
215     flx_NS(i,j)=
216     & tmp_VbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))
217     & -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*Half
218     & *_recip_dyC(i,j,bi,bj)
219     & *(hFacStmp*(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
220     C- Problem here: No-slip bc CANNOT be written in term of a flux
221 mlosch 1.15 & +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
222 jmc 1.28 & *wVel(i,j,k,bi,bj))
223     & +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*Half
224     & *_recip_dyC(i,j,bi,bj)*(del2w(i,j)-del2w(i,j-1))
225 mlosch 1.18 #ifdef ISOTROPIC_COS_SCALING
226     #ifdef COSINEMETH_III
227     & *sqCosFacV(j,bi,bj)
228     #else
229     & *CosFacV(j,bi,bj)
230     #endif
231     #endif
232 mlosch 1.15 C The last term is the weighted average of the viscous stress at the open
233     C fraction of the w control volume and at the closed fraction of the
234     C the control volume. A more compact but less intelligible version
235     C of the last three lines is:
236     CML & *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
237 jmc 1.28 CML & *wVel(i,j,k,bi,bi) + hFacStmp*wVel(i,j-1,k,bi,bj) )
238 adcroft 1.1 ENDDO
239 jmc 1.28 ENDDO
240 adcroft 1.1 C Flux on Western face
241 jmc 1.28 DO j=jMin,jMax
242     DO i=iMin,iMax+1
243 mlosch 1.15 C First compute the fraction of open water for the w-control volume
244     C at the western face
245 jmc 1.28 hFacWtmp=max(_hFacW(i,j,k-1,bi,bj)-Half,0. _d 0)
246     & + min(_hFacW(i,j,k ,bi,bj),Half)
247 jmc 1.21 tmp_UbarZ=Half*(
248 jmc 1.28 & _hFacW(i,j,k-1,bi,bj)*uVel( i ,j,k-1,bi,bj)
249     & +_hFacW(i,j, k ,bi,bj)*uVel( i ,j, k ,bi,bj))
250     flx_EW(i,j)=
251     & tmp_UbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))
252     & -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*Half
253     & *_recip_dxC(i,j,bi,bj)
254     & *(hFacWtmp*(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
255     C- Problem here: No-slip bc CANNOT be written in term of a flux
256 mlosch 1.15 & +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
257 jmc 1.28 & *wVel(i,j,k,bi,bj) )
258     & +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*Half
259     & *_recip_dxC(i,j,bi,bj)*(del2w(i,j)-del2w(i-1,j))
260 mlosch 1.18 #ifdef COSINEMETH_III
261     & *sqCosFacU(j,bi,bj)
262 jmc 1.25 #else
263 mlosch 1.18 & *CosFacU(j,bi,bj)
264     #endif
265 mlosch 1.15 C The last term is the weighted average of the viscous stress at the open
266     C fraction of the w control volume and at the closed fraction of the
267     C the control volume. A more compact but less intelligible version
268     C of the last three lines is:
269     CML & *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacWtmp))
270 jmc 1.28 CML & *wVel(i,j,k,bi,bi) + hFacWtmp*wVel(i-1,j,k,bi,bj) )
271 adcroft 1.1 ENDDO
272 jmc 1.28 ENDDO
273 adcroft 1.1 C Flux on Lower face
274 jmc 1.28 DO j=jMin,jMax
275     DO i=iMin,iMax
276 baylor 1.27 C Interpolate vert viscosity to W points
277 jmc 1.28 viscLoc = ( KappaRU(i,j,k) +KappaRU(i+1,j,k)
278     & +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
279     & +KappaRV(i,j,k) +KappaRV(i,j+1,k)
280     & +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
281     & )*0.125 _d 0
282     tmp_WbarZ = Half*( wVel(i,j, k ,bi,bj)
283     & +wVel(i,j,kp1,bi,bj)*wOverRide )
284     flx_Dn(i,j)=
285     & tmp_WbarZ*tmp_WbarZ
286     & -viscLoc*recip_drF(k)
287     & *( wVel(i,j, k ,bi,bj)
288     & -wVel(i,j,kp1,bi,bj)*wOverRide )
289 adcroft 1.1 ENDDO
290 jmc 1.28 ENDDO
291     C Divergence of fluxes
292     DO j=jMin,jMax
293     DO i=iMin,iMax
294     gW(i,j,k,bi,bj) = 0.
295 adcroft 1.1 & -(
296 jmc 1.28 & +_recip_dxF(i,j,bi,bj)*(
297     & flx_EW(i+1,j)-flx_EW(i,j) )
298     & +_recip_dyF(i,j,bi,bj)*(
299     & flx_NS(i,j+1)-flx_NS(i,j) )
300     & +recip_drC(k) *(
301     & flx_Up(i,j) -flx_Dn(i,j) )
302 adcroft 1.1 & )
303 jmc 1.28 caja * recip_hFacU(i,j,k,bi,bj)
304 jmc 1.25 caja NOTE: This should be included
305 adcroft 1.1 caja but we need an hFacUW (above U points)
306     caja and an hFacUS (above V points) too...
307 jmc 1.28 C-- prepare for next level (k+1)
308     flx_Up(i,j)=flx_Dn(i,j)
309 adcroft 1.1 ENDDO
310 jmc 1.28 ENDDO
311 adcroft 1.1
312 jmc 1.25 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313    
314     C- Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
315     C and save gW_[n] into gwNm1 for the next time step.
316     c#ifdef ALLOW_ADAMSBASHFORTH_3
317 jmc 1.28 c CALL ADAMS_BASHFORTH3(
318     c I bi, bj, k,
319     c U gW, gwNm,
320     c I momStartAB, myIter, myThid )
321 jmc 1.25 c#else /* ALLOW_ADAMSBASHFORTH_3 */
322 jmc 1.28 CALL ADAMS_BASHFORTH2(
323     I bi, bj, k,
324     U gW, gwNm1,
325     I myIter, myThid )
326 jmc 1.25 c#endif /* ALLOW_ADAMSBASHFORTH_3 */
327 jmc 1.21
328 jmc 1.25 C- end of the k loop
329 adcroft 1.4 ENDDO
330    
331 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
332    
333     RETURN
334     END

  ViewVC Help
Powered by ViewVC 1.1.22