/[MITgcm]/MITgcm_contrib/dgoldberg/depth_control_no_nsa/code/obcs_balance_flow.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/depth_control_no_nsa/code/obcs_balance_flow.F

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


Revision 1.1 - (show annotations) (download)
Thu Dec 7 23:21:13 2017 UTC (6 years, 5 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
test case for depth control w/out cg2d_nsa

1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_balance_flow.F,v 1.8 2016/07/06 22:48:56 jmc Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7 CBOP
8 C !ROUTINE: OBCS_BALANCE_FLOW
9
10 C !INTERFACE:
11 SUBROUTINE OBCS_BALANCE_FLOW( myTime, myIter, myThid )
12
13 C !DESCRIPTION:
14 C *==========================================================*
15 C | SUBROUTINE OBCS_BALANCE_FLOW
16 C | o Modify OB normal flow to ensure no net inflow
17 C *==========================================================*
18
19 C !USES:
20 IMPLICIT NONE
21
22 C === Global variables ===
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "OBCS_PARAMS.h"
28 #include "OBCS_GRID.h"
29 #include "OBCS_FIELDS.h"
30
31 C !INPUT/OUTPUT PARAMETERS:
32 _RL myTime
33 INTEGER myIter
34 INTEGER myThid
35 CEOP
36
37 #ifdef ALLOW_OBCS
38 #ifdef ALLOW_OBCS_BALANCE
39
40 C !FUNCTIONS:
41
42 C !LOCAL VARIABLES:
43 C bi, bj :: tile indices
44 C i,j,k :: loop indices
45 C iB, jB :: local index of open boundary
46 C msgBuf :: Informational/error message buffer
47 INTEGER bi, bj
48 INTEGER i, j, k, iB, jB
49 CHARACTER*(MAX_LEN_MBUF) msgBuf
50 _RL areaOB, areaE, areaW, areaN, areaS, tmpA
51 _RL inFlow, flowE, flowW, flowN, flowS
52 _RL tileArea(nSx,nSy)
53 _RL tileFlow(nSx,nSy)
54 _RL tileAreaOB(nSx,nSy)
55 _RL tileInFlow(nSx,nSy)
56 LOGICAL flag
57
58 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
59 C-- Old method (OBCS_balanceFac < 0): balance each OB separately
60 C-- New method applies to all OB with BCS_balanceFac >=0 :
61 C ensure that the net inflow through all OB is balanced by correcting
62 C each OB normal flow with a uniform velocity, using the corresponding
63 C weight factor OBCS_balanceFac.
64 C e.g., OBCS_balanceFac_E,W,N,S= 1, -1, 2, 0 :
65 C => correct Western OBWu (by substracting a uniform velocity) to ensure
66 C zero net transport through Western OB
67 C => correct Eastern and Northern normal flow (twice larger Northern
68 C velocity correction than Eastern correction) to ensure that
69 C the total inflow through E+N+S OB is balanced
70
71 #ifdef ALLOW_DEBUG
72 IF (debugMode) CALL DEBUG_ENTER('OBCS_BALANCE_FLOW',myThid)
73 #endif
74
75 C-- Integrate the transport through each OB
76 flag = .FALSE.
77 areaOB = 0. _d 0
78 inFlow = 0. _d 0
79 DO bj=myByLo(myThid),myByHi(myThid)
80 DO bi=myBxLo(myThid),myBxHi(myThid)
81 tileAreaOB(bi,bj) = 0.
82 tileInFlow(bi,bj) = 0.
83 ENDDO
84 ENDDO
85
86 #ifdef ALLOW_OBCS_EAST
87 areaE = 0. _d 0
88 flowE = 0. _d 0
89 flag = flag .OR. ( OBCS_balanceFacE.GT.0. )
90 DO bj=myByLo(myThid),myByHi(myThid)
91 DO bi=myBxLo(myThid),myBxHi(myThid)
92 tileArea(bi,bj) = 0.
93 tileFlow(bi,bj) = 0.
94 IF ( tileHasOBE(bi,bj) ) THEN
95 DO k=1,Nr
96 DO j=1,sNy
97 iB = OB_Ie(j,bi,bj)
98 C- If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
99 C communicates with tile interior (sNx+1) rather than with halo region (i=1)
100 IF ( iB.NE.OB_indexNone .AND. iB.GT.1 ) THEN
101 tmpA = drF(k)*_hFacW(iB,j,k,bi,bj)*dyG(iB,j,bi,bj)
102 & *maskInW(iB,j,bi,bj)
103 tileArea(bi,bj) = tileArea(bi,bj) + tmpA
104 tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBEu(j,k,bi,bj)
105 ENDIF
106 ENDDO
107 ENDDO
108 IF ( OBCS_balanceFacE.GE.0. ) THEN
109 tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
110 tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
111 & + tileArea(bi,bj)*OBCS_balanceFacE
112 ENDIF
113 areaE = areaE + tileArea(bi,bj)
114 flowE = flowE + tileFlow(bi,bj)
115 ENDIF
116 ENDDO
117 ENDDO
118 c WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
119 c & 'OBCS_balance it,areaE,flowE=', myIter, areaE, flowE
120 IF ( OBCS_balanceFacE.LT.0. ) THEN
121 CALL GLOBAL_SUM_TILE_RL( tileArea, areaE, myThid )
122 IF ( areaE.GT.0. ) THEN
123 CALL GLOBAL_SUM_TILE_RL( tileFlow, flowE, myThid )
124 IF ( debugLevel.GE.debLevC ) THEN
125 WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
126 & myIter, ' ) correct OBEu:', flowE, -flowE/areaE
127 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
128 & SQUEEZE_RIGHT, myThid )
129 ENDIF
130 flowE = -flowE/areaE
131 ENDIF
132 ENDIF
133 #endif /* ALLOW_OBCS_EAST */
134
135 #ifdef ALLOW_OBCS_WEST
136 areaW = 0. _d 0
137 flowW = 0. _d 0
138 flag = flag .OR. ( OBCS_balanceFacW.GT.0. )
139 DO bj=myByLo(myThid),myByHi(myThid)
140 DO bi=myBxLo(myThid),myBxHi(myThid)
141 tileArea(bi,bj) = 0.
142 tileFlow(bi,bj) = 0.
143 IF ( tileHasOBW(bi,bj) ) THEN
144 DO k=1,Nr
145 DO j=1,sNy
146 iB = OB_Iw(j,bi,bj)
147 C- If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
148 C communicates with tile interior (i=0) rather than with halo region (i=sNx)
149 IF ( iB.NE.OB_indexNone .AND. iB.LT.sNx ) THEN
150 tmpA = drF(k)*_hFacW(1+iB,j,k,bi,bj)*dyG(1+iB,j,bi,bj)
151 & *maskInW(1+iB,j,bi,bj)
152 tileArea(bi,bj) = tileArea(bi,bj) + tmpA
153 tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBWu(j,k,bi,bj)
154 ENDIF
155 ENDDO
156 ENDDO
157 IF ( OBCS_balanceFacW.GE.0. ) THEN
158 tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
159 tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
160 & + tileArea(bi,bj)*OBCS_balanceFacW
161 ENDIF
162 areaW = areaW + tileArea(bi,bj)
163 flowW = flowW + tileFlow(bi,bj)
164 ENDIF
165 ENDDO
166 ENDDO
167 c WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
168 c & 'OBCS_balance it,areaW,flowW=', myIter, areaW, flowW
169 IF ( OBCS_balanceFacW.LT.0. ) THEN
170 CALL GLOBAL_SUM_TILE_RL( tileArea, areaW, myThid )
171 IF ( areaW.GT.0. ) THEN
172 CALL GLOBAL_SUM_TILE_RL( tileFlow, flowW, myThid )
173 IF ( debugLevel.GE.debLevC ) THEN
174 WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
175 & myIter, ' ) correct OBWu:', flowW, -flowW/areaW
176 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
177 & SQUEEZE_RIGHT, myThid )
178 ENDIF
179 flowW = -flowW/areaW
180 ENDIF
181 ENDIF
182 #endif /* ALLOW_OBCS_WEST */
183
184 #ifdef ALLOW_OBCS_NORTH
185 areaN = 0. _d 0
186 flowN = 0. _d 0
187 flag = flag .OR. ( OBCS_balanceFacN.GT.0. )
188 DO bj=myByLo(myThid),myByHi(myThid)
189 DO bi=myBxLo(myThid),myBxHi(myThid)
190 tileArea(bi,bj) = 0.
191 tileFlow(bi,bj) = 0.
192 IF ( tileHasOBN(bi,bj) ) THEN
193 DO k=1,Nr
194 DO i=1,sNx
195 jB = OB_Jn(i,bi,bj)
196 C- If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
197 C communicates with tile interior (sNy+1) rather than with halo region (j=1)
198 IF ( jB.NE.OB_indexNone .AND. jB.GT.1 ) THEN
199 tmpA = drF(k)*_hFacS(i,jB,k,bi,bj)*dxG(i,jB,bi,bj)
200 & *maskInS(i,jB,bi,bj)
201 tileArea(bi,bj) = tileArea(bi,bj) + tmpA
202 tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBNv(i,k,bi,bj)
203 ENDIF
204 ENDDO
205 ENDDO
206 IF ( OBCS_balanceFacN.GE.0. ) THEN
207 tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
208 tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
209 & + tileArea(bi,bj)*OBCS_balanceFacN
210 ENDIF
211 areaN = areaN + tileArea(bi,bj)
212 flowN = flowN + tileFlow(bi,bj)
213 ENDIF
214 ENDDO
215 ENDDO
216 c WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
217 c & 'OBCS_balance it,areaN,flowN=', myIter, areaN, flowN
218 IF ( OBCS_balanceFacN.LT.0. ) THEN
219 CALL GLOBAL_SUM_TILE_RL( tileArea, areaN, myThid )
220 IF ( areaN.GT.0. ) THEN
221 CALL GLOBAL_SUM_TILE_RL( tileFlow, flowN, myThid )
222 IF ( debugLevel.GE.debLevC ) THEN
223 WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
224 & myIter, ' ) correct OBNv:', flowN, -flowN/areaN
225 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
226 & SQUEEZE_RIGHT, myThid )
227 ENDIF
228 flowN = -flowN/areaN
229 ENDIF
230 ENDIF
231 #endif /* ALLOW_OBCS_NORTH */
232
233 #ifdef ALLOW_OBCS_SOUTH
234 areaS = 0. _d 0
235 flowS = 0. _d 0
236 flag = flag .OR. ( OBCS_balanceFacS.GT.0. )
237 DO bj=myByLo(myThid),myByHi(myThid)
238 DO bi=myBxLo(myThid),myBxHi(myThid)
239 tileArea(bi,bj) = 0.
240 tileFlow(bi,bj) = 0.
241 IF ( tileHasOBS(bi,bj) ) THEN
242 DO k=1,Nr
243 DO i=1,sNx
244 jB = OB_Js(i,bi,bj)
245 C- If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
246 C communicates with tile interior (j=0) rather than with halo region (j=sNy)
247 IF ( jB.NE.OB_indexNone .AND. jB.LT.sNy ) THEN
248 tmpA = drF(k)*_hFacS(i,1+jB,k,bi,bj)*dxG(i,1+jB,bi,bj)
249 & *maskInS(i,1+jB,bi,bj)
250 tileArea(bi,bj) = tileArea(bi,bj) + tmpA
251 tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBSv(i,k,bi,bj)
252 ENDIF
253 ENDDO
254 ENDDO
255 IF ( OBCS_balanceFacS.GE.0. ) THEN
256 tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
257 tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
258 & + tileArea(bi,bj)*OBCS_balanceFacS
259 ENDIF
260 areaS = areaS + tileArea(bi,bj)
261 flowS = flowS + tileFlow(bi,bj)
262 ENDIF
263 ENDDO
264 ENDDO
265 c WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
266 c & 'OBCS_balance it,areaS,flowS=', myIter, areaS, flowS
267 IF ( OBCS_balanceFacS.LT.0. ) THEN
268 CALL GLOBAL_SUM_TILE_RL( tileArea, areaS, myThid )
269 IF ( areaS.GT.0. ) THEN
270 CALL GLOBAL_SUM_TILE_RL( tileFlow, flowS, myThid )
271 IF ( debugLevel.GE.debLevC ) THEN
272 WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
273 & myIter, ' ) correct OBSv:', flowS, -flowS/areaS
274 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
275 & SQUEEZE_RIGHT, myThid )
276 ENDIF
277 flowS = -flowS/areaS
278 ENDIF
279 ENDIF
280 #endif /* ALLOW_OBCS_SOUTH */
281
282 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
283 C-- Calculate a unique velocity correction for all (OBCS_balanceFac>0) OB
284 C and correct each OB using corresponding weight factor OBCS_balanceFac
285
286 IF ( flag ) CALL GLOBAL_SUM_TILE_RL( tileAreaOB, areaOB, myThid )
287 IF ( areaOB.GT.0. ) THEN
288 CALL GLOBAL_SUM_TILE_RL( tileInFlow, inFlow, myThid )
289 IF ( debugLevel.GE.debLevC ) THEN
290 WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
291 & myIter, ' ) correct for inFlow:', inFlow, inFlow/areaOB
292 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
293 & SQUEEZE_RIGHT, myThid )
294 ENDIF
295 inFlow = inFlow / areaOB
296 ENDIF
297 IF ( OBCS_balanceFacE.GE.0. ) flowE = inFlow*OBCS_balanceFacE
298 IF ( OBCS_balanceFacW.GE.0. ) flowW = -inFlow*OBCS_balanceFacW
299 IF ( OBCS_balanceFacN.GE.0. ) flowN = inFlow*OBCS_balanceFacN
300 IF ( OBCS_balanceFacS.GE.0. ) flowS = -inFlow*OBCS_balanceFacS
301
302 IF ( debugLevel.GE.debLevC .AND. areaOB.GT.0. ) THEN
303 WRITE(msgBuf,'(A,1P2E16.8)')
304 & 'OBCS_balance correction to OBEu,OBWu:', flowE, flowW
305 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
306 & SQUEEZE_RIGHT, myThid )
307 WRITE(msgBuf,'(A,1P2E16.8)')
308 & 'OBCS_balance correction to OBNv,OBSv:', flowN, flowS
309 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
310 & SQUEEZE_RIGHT, myThid )
311 ENDIF
312
313 c IF ( .NOT.useOBCSbalance ) RETURN
314
315 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
316 C-- Add correction:
317
318 #ifdef ALLOW_OBCS_EAST
319 IF ( OBCS_balanceFacE.NE.0. ) THEN
320 DO bj=myByLo(myThid),myByHi(myThid)
321 DO bi=myBxLo(myThid),myBxHi(myThid)
322 IF ( tileHasOBE(bi,bj) ) THEN
323 DO k=1,Nr
324 DO j=1-OLy,sNy+OLy
325 iB = OB_Ie(j,bi,bj)
326 IF ( iB.NE.OB_indexNone ) THEN
327 OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
328 & + flowE*maskW(iB,j,k,bi,bj)
329 ENDIF
330 ENDDO
331 ENDDO
332 ENDIF
333 ENDDO
334 ENDDO
335 ENDIF
336 #endif /* ALLOW_OBCS_EAST */
337
338 #ifdef ALLOW_OBCS_WEST
339 IF ( OBCS_balanceFacW.NE.0. ) THEN
340 DO bj=myByLo(myThid),myByHi(myThid)
341 DO bi=myBxLo(myThid),myBxHi(myThid)
342 IF ( tileHasOBW(bi,bj) ) THEN
343 DO k=1,Nr
344 DO j=1-OLy,sNy+OLy
345 iB = OB_Iw(j,bi,bj)
346 IF ( iB.NE.OB_indexNone ) THEN
347 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
348 & + flowW*maskW(1+iB,j,k,bi,bj)
349 ENDIF
350 ENDDO
351 ENDDO
352 ENDIF
353 ENDDO
354 ENDDO
355 ENDIF
356 #endif /* ALLOW_OBCS_WEST */
357
358 #ifdef ALLOW_OBCS_NORTH
359 IF ( OBCS_balanceFacN.NE.0. ) THEN
360 DO bj=myByLo(myThid),myByHi(myThid)
361 DO bi=myBxLo(myThid),myBxHi(myThid)
362 IF ( tileHasOBN(bi,bj) ) THEN
363 DO k=1,Nr
364 DO i=1-OLx,sNx+OLx
365 jB = OB_Jn(i,bi,bj)
366 IF ( jB.NE.OB_indexNone ) THEN
367 OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
368 & + flowN*maskS(i,jB,k,bi,bj)
369 ENDIF
370 ENDDO
371 ENDDO
372 ENDIF
373 ENDDO
374 ENDDO
375 ENDIF
376 #endif /* ALLOW_OBCS_NORTH */
377
378 #ifdef ALLOW_OBCS_SOUTH
379 IF ( OBCS_balanceFacS.NE.0. ) THEN
380 DO bj=myByLo(myThid),myByHi(myThid)
381 DO bi=myBxLo(myThid),myBxHi(myThid)
382 IF ( tileHasOBS(bi,bj) ) THEN
383 DO k=1,Nr
384 DO i=1-OLx,sNx+OLx
385 jB = OB_Js(i,bi,bj)
386 IF ( jB.NE.OB_indexNone ) THEN
387 OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
388 & + flowS*maskS(i,1+jB,k,bi,bj)
389 ENDIF
390 ENDDO
391 ENDDO
392 ENDIF
393 ENDDO
394 ENDDO
395 ENDIF
396 #endif /* ALLOW_OBCS_SOUTH */
397
398 #ifdef ALLOW_DEBUG
399 IF (debugMode) CALL DEBUG_LEAVE('OBCS_BALANCE_FLOW',myThid)
400 #endif
401
402 #endif /* ALLOW_OBCS_BALANCE */
403 #endif /* ALLOW_OBCS */
404
405 RETURN
406 END

  ViewVC Help
Powered by ViewVC 1.1.22