/[MITgcm]/MITgcm_contrib/verification_other/shelfice_remeshing/code/obcs_balance_flow.F
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/shelfice_remeshing/code/obcs_balance_flow.F

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


Revision 1.2 - (show annotations) (download)
Wed Jul 6 18:03:40 2016 UTC (9 years ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65y, checkpoint67a, checkpoint67b, checkpoint67d, HEAD
Changes since 1.1: +1 -3 lines
new files for vertical remeshing test only

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

  ViewVC Help
Powered by ViewVC 1.1.22