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 |