/[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.28 - (hide annotations) (download)
Fri Jul 7 20:09:37 2006 UTC (17 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.27: +123 -123 lines
- take bi,bj loops outside of calc_gw (since KappaRU,V are local array)
- start cleaning-up but leave the bugs (so that results don't change)

1 jmc 1.28 C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.27 2006/06/20 20:57:37 baylor 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.11
100 jmc 1.28 C-- Initialise gW to zero
101     DO k=1,Nr
102     DO j=1-OLy,sNy+OLy
103     DO i=1-OLx,sNx+OLx
104 adcroft 1.1 gW(i,j,k,bi,bj) = 0.
105     ENDDO
106     ENDDO
107 jmc 1.28 ENDDO
108 adcroft 1.1
109 jmc 1.28 C-- Boundaries condition at top
110     DO j=jMin,jMax
111     DO i=iMin,iMax
112     flx_Up(i,j)=0.
113     ENDDO
114     ENDDO
115 adcroft 1.1
116 jmc 1.28 C--- Sweep down column
117     DO k=2,Nr
118     kp1=k+1
119     wOverRide=1.
120     IF (k.EQ.Nr) THEN
121     kp1=Nr
122 adcroft 1.1 wOverRide=0.
123 jmc 1.28 ENDIF
124     C-- horizontal bi-harmonic dissipation
125     IF (momViscosity .AND. viscA4W.NE.0. ) THEN
126     C- calculate the horizontal Laplacian of vertical flow
127 mlosch 1.18 C Zonal flux d/dx W
128     DO j=1-Oly,sNy+Oly
129     fZon(1-Olx,j)=0.
130     DO i=1-Olx+1,sNx+Olx
131 jmc 1.28 C- Problem here: drF(k)*_hFacC & fZon are not at the same Horiz.&Vert. location
132 mlosch 1.18 fZon(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
133     & *_dyG(i,j,bi,bj)
134     & *_recip_dxC(i,j,bi,bj)
135     & *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
136     #ifdef COSINEMETH_III
137 jmc 1.28 & *sqcosFacU(j,bi,bj)
138 mlosch 1.18 #endif
139     ENDDO
140 jmc 1.28 ENDDO
141 mlosch 1.18 C Meridional flux d/dy W
142     DO i=1-Olx,sNx+Olx
143 jmc 1.28 fMer(i,1-Oly)=0.
144 mlosch 1.18 ENDDO
145     DO j=1-Oly+1,sNy+Oly
146     DO i=1-Olx,sNx+Olx
147 jmc 1.28 C- Problem here: drF(k)*_hFacC & fMer are not at the same Horiz.&Vert. location
148 mlosch 1.18 fMer(i,j) = drF(k)*_hFacC(i,j,k,bi,bj)
149     & *_dxG(i,j,bi,bj)
150     & *_recip_dyC(i,j,bi,bj)
151     & *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
152     #ifdef ISOTROPIC_COS_SCALING
153     #ifdef COSINEMETH_III
154     & *sqCosFacV(j,bi,bj)
155     #endif
156     #endif
157     ENDDO
158     ENDDO
159 jmc 1.25
160 mlosch 1.18 C del^2 W
161     C Difference of zonal fluxes ...
162     DO j=1-Oly,sNy+Oly
163     DO i=1-Olx,sNx+Olx-1
164     del2w(i,j)=fZon(i+1,j)-fZon(i,j)
165     ENDDO
166     del2w(sNx+Olx,j)=0.
167     ENDDO
168    
169     C ... add difference of meridional fluxes and divide by volume
170     DO j=1-Oly,sNy+Oly-1
171     DO i=1-Olx,sNx+Olx
172     C First compute the fraction of open water for the w-control volume
173     C at the southern face
174 jmc 1.28 hFacCtmp=max( _hFacC(i,j,k-1,bi,bj)-Half,0. _d 0 )
175     & + min( _hFacC(i,j,k ,bi,bj) ,Half )
176 mlosch 1.18 IF (hFacCtmp .GT. 0.) THEN
177     recip_hFacCtmp = 1./hFacCtmp
178     ELSE
179     recip_hFacCtmp = 0. _d 0
180     ENDIF
181     del2w(i,j)=recip_rA(i,j,bi,bj)
182     & *recip_drC(k)*recip_hFacCtmp
183     & *(
184     & del2w(i,j)
185     & +( fMer(i,j+1)-fMer(i,j) )
186     & )
187     ENDDO
188     ENDDO
189     C-- No-slip BCs impose a drag at walls...
190     CML ************************************************************
191     CML No-slip Boundary conditions for bi-harmonic dissipation
192     CML need to be implemented here!
193     CML ************************************************************
194 jmc 1.28 ELSE
195 jmc 1.21 C- Initialize del2w to zero:
196     DO j=1-Oly,sNy+Oly
197     DO i=1-Olx,sNx+Olx
198     del2w(i,j) = 0. _d 0
199     ENDDO
200     ENDDO
201 jmc 1.28 ENDIF
202 mlosch 1.18
203 adcroft 1.1 C Flux on Southern face
204 jmc 1.28 DO j=jMin,jMax+1
205     DO i=iMin,iMax
206 mlosch 1.15 C First compute the fraction of open water for the w-control volume
207     C at the southern face
208 jmc 1.28 hFacStmp=max(_hFacS(i,j,k-1,bi,bj)-Half,0. _d 0)
209     & + min(_hFacS(i,j,k ,bi,bj),Half)
210 adcroft 1.1 tmp_VbarZ=Half*(
211 jmc 1.28 & _hFacS(i,j,k-1,bi,bj)*vVel( i ,j,k-1,bi,bj)
212     & +_hFacS(i,j, k ,bi,bj)*vVel( i ,j, k ,bi,bj))
213     flx_NS(i,j)=
214     & tmp_VbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))
215     & -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*Half
216     & *_recip_dyC(i,j,bi,bj)
217     & *(hFacStmp*(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj))
218     C- Problem here: No-slip bc CANNOT be written in term of a flux
219 mlosch 1.15 & +(1. _d 0 - hFacStmp)*(1. _d 0 - slipSideFac)
220 jmc 1.28 & *wVel(i,j,k,bi,bj))
221     & +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*Half
222     & *_recip_dyC(i,j,bi,bj)*(del2w(i,j)-del2w(i,j-1))
223 mlosch 1.18 #ifdef ISOTROPIC_COS_SCALING
224     #ifdef COSINEMETH_III
225     & *sqCosFacV(j,bi,bj)
226     #else
227     & *CosFacV(j,bi,bj)
228     #endif
229     #endif
230 mlosch 1.15 C The last term is the weighted average of the viscous stress at the open
231     C fraction of the w control volume and at the closed fraction of the
232     C the control volume. A more compact but less intelligible version
233     C of the last three lines is:
234     CML & *( (1 _d 0 - slipSideFac*(1 _d 0 - hFacStmp))
235 jmc 1.28 CML & *wVel(i,j,k,bi,bi) + hFacStmp*wVel(i,j-1,k,bi,bj) )
236 adcroft 1.1 ENDDO
237 jmc 1.28 ENDDO
238 adcroft 1.1 C Flux on Western face
239 jmc 1.28 DO j=jMin,jMax
240     DO i=iMin,iMax+1
241 mlosch 1.15 C First compute the fraction of open water for the w-control volume
242     C at the western face
243 jmc 1.28 hFacWtmp=max(_hFacW(i,j,k-1,bi,bj)-Half,0. _d 0)
244     & + min(_hFacW(i,j,k ,bi,bj),Half)
245 jmc 1.21 tmp_UbarZ=Half*(
246 jmc 1.28 & _hFacW(i,j,k-1,bi,bj)*uVel( i ,j,k-1,bi,bj)
247     & +_hFacW(i,j, k ,bi,bj)*uVel( i ,j, k ,bi,bj))
248     flx_EW(i,j)=
249     & tmp_UbarZ*Half*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))
250     & -(viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*Half
251     & *_recip_dxC(i,j,bi,bj)
252     & *(hFacWtmp*(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj))
253     C- Problem here: No-slip bc CANNOT be written in term of a flux
254 mlosch 1.15 & +(1 _d 0 - hFacWtmp)*(1 _d 0 - slipSideFac)
255 jmc 1.28 & *wVel(i,j,k,bi,bj) )
256     & +(viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*Half
257     & *_recip_dxC(i,j,bi,bj)*(del2w(i,j)-del2w(i-1,j))
258 mlosch 1.18 #ifdef COSINEMETH_III
259     & *sqCosFacU(j,bi,bj)
260 jmc 1.25 #else
261 mlosch 1.18 & *CosFacU(j,bi,bj)
262     #endif
263 mlosch 1.15 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 jmc 1.28 CML & *wVel(i,j,k,bi,bi) + hFacWtmp*wVel(i-1,j,k,bi,bj) )
269 adcroft 1.1 ENDDO
270 jmc 1.28 ENDDO
271 adcroft 1.1 C Flux on Lower face
272 jmc 1.28 DO j=jMin,jMax
273     DO i=iMin,iMax
274 baylor 1.27 C Interpolate vert viscosity to W points
275 jmc 1.28 viscLoc = ( KappaRU(i,j,k) +KappaRU(i+1,j,k)
276     & +KappaRU(i,j,kp1)+KappaRU(i+1,j,kp1)
277     & +KappaRV(i,j,k) +KappaRV(i,j+1,k)
278     & +KappaRV(i,j,kp1)+KappaRV(i,j+1,kp1)
279     & )*0.125 _d 0
280     tmp_WbarZ = Half*( wVel(i,j, k ,bi,bj)
281     & +wVel(i,j,kp1,bi,bj)*wOverRide )
282     flx_Dn(i,j)=
283     & tmp_WbarZ*tmp_WbarZ
284     & -viscLoc*recip_drF(k)
285     & *( wVel(i,j, k ,bi,bj)
286     & -wVel(i,j,kp1,bi,bj)*wOverRide )
287 adcroft 1.1 ENDDO
288 jmc 1.28 ENDDO
289     C Divergence of fluxes
290     DO j=jMin,jMax
291     DO i=iMin,iMax
292     gW(i,j,k,bi,bj) = 0.
293 adcroft 1.1 & -(
294 jmc 1.28 & +_recip_dxF(i,j,bi,bj)*(
295     & flx_EW(i+1,j)-flx_EW(i,j) )
296     & +_recip_dyF(i,j,bi,bj)*(
297     & flx_NS(i,j+1)-flx_NS(i,j) )
298     & +recip_drC(k) *(
299     & flx_Up(i,j) -flx_Dn(i,j) )
300 adcroft 1.1 & )
301 jmc 1.28 caja * recip_hFacU(i,j,k,bi,bj)
302 jmc 1.25 caja NOTE: This should be included
303 adcroft 1.1 caja but we need an hFacUW (above U points)
304     caja and an hFacUS (above V points) too...
305 jmc 1.28 C-- prepare for next level (k+1)
306     flx_Up(i,j)=flx_Dn(i,j)
307 adcroft 1.1 ENDDO
308 jmc 1.28 ENDDO
309 adcroft 1.1
310 jmc 1.25 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
311    
312     C- Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights)
313     C and save gW_[n] into gwNm1 for the next time step.
314     c#ifdef ALLOW_ADAMSBASHFORTH_3
315 jmc 1.28 c CALL ADAMS_BASHFORTH3(
316     c I bi, bj, k,
317     c U gW, gwNm,
318     c I momStartAB, myIter, myThid )
319 jmc 1.25 c#else /* ALLOW_ADAMSBASHFORTH_3 */
320 jmc 1.28 CALL ADAMS_BASHFORTH2(
321     I bi, bj, k,
322     U gW, gwNm1,
323     I myIter, myThid )
324 jmc 1.25 c#endif /* ALLOW_ADAMSBASHFORTH_3 */
325 jmc 1.21
326 jmc 1.25 C- end of the k loop
327 adcroft 1.4 ENDDO
328    
329 adcroft 1.1 #endif /* ALLOW_NONHYDROSTATIC */
330    
331     RETURN
332     END

  ViewVC Help
Powered by ViewVC 1.1.22