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

Annotation 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.1 - (hide annotations) (download)
Thu May 5 18:17:26 2016 UTC (9 years, 2 months ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65w
rename files, add kate's change to obcs_balance

1 dgoldberg 1.1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_balance_flow.F,v 1.7 2012/09/17 22:03:18 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     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     #ifdef ALLOW_SHELFICE_REMESHING
311     IF ( conserve_ssh ) THEN
312     SEALEVEL = 0. _d 0
313     ETA = 0. _d 0
314     CALL SHELFICE_SEA_LEVEL_AVG( SEALEVEL,ETA, myThid )
315     C Restore the open ocean ETA sum to 0
316     IF (ETA .NE. 0. _d 0) THEN
317     inFlow = inFlow + ETA/deltaT
318     ENDIF
319     ENDIF
320     #endif /* */
321     #endif /* ALLOW_SHELFICE */
322     C KS16 end-----------------------------------------------------
323    
324     inFlow = inFlow / areaOB
325     ENDIF
326     IF ( OBCS_balanceFacE.GE.0. ) flowE = inFlow*OBCS_balanceFacE
327     IF ( OBCS_balanceFacW.GE.0. ) flowW = -inFlow*OBCS_balanceFacW
328     IF ( OBCS_balanceFacN.GE.0. ) flowN = inFlow*OBCS_balanceFacN
329     IF ( OBCS_balanceFacS.GE.0. ) flowS = -inFlow*OBCS_balanceFacS
330    
331     IF ( debugLevel.GE.debLevC .AND. areaOB.GT.0. ) THEN
332     WRITE(msgBuf,'(A,1P2E16.8)')
333     & 'OBCS_balance correction to OBEu,OBWu:', flowE, flowW
334     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
335     & SQUEEZE_RIGHT, myThid )
336     WRITE(msgBuf,'(A,1P2E16.8)')
337     & 'OBCS_balance correction to OBNv,OBSv:', flowN, flowS
338     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
339     & SQUEEZE_RIGHT, myThid )
340     ENDIF
341    
342     c IF ( .NOT.useOBCSbalance ) RETURN
343    
344     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
345     C-- Add correction:
346    
347     #ifdef ALLOW_OBCS_EAST
348     IF ( OBCS_balanceFacE.NE.0. ) THEN
349     DO bj=myByLo(myThid),myByHi(myThid)
350     DO bi=myBxLo(myThid),myBxHi(myThid)
351     IF ( tileHasOBE(bi,bj) ) THEN
352     DO k=1,Nr
353     DO j=1-OLy,sNy+OLy
354     iB = OB_Ie(j,bi,bj)
355     IF ( iB.NE.OB_indexNone ) THEN
356     OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
357     & + flowE*maskW(iB,j,k,bi,bj)
358     ENDIF
359     ENDDO
360     ENDDO
361     ENDIF
362     ENDDO
363     ENDDO
364     ENDIF
365     #endif /* ALLOW_OBCS_EAST */
366    
367     #ifdef ALLOW_OBCS_WEST
368     IF ( OBCS_balanceFacW.NE.0. ) THEN
369     DO bj=myByLo(myThid),myByHi(myThid)
370     DO bi=myBxLo(myThid),myBxHi(myThid)
371     IF ( tileHasOBW(bi,bj) ) THEN
372     DO k=1,Nr
373     DO j=1-OLy,sNy+OLy
374     iB = OB_Iw(j,bi,bj)
375     IF ( iB.NE.OB_indexNone ) THEN
376     OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
377     & + flowW*maskW(1+iB,j,k,bi,bj)
378     ENDIF
379     ENDDO
380     ENDDO
381     ENDIF
382     ENDDO
383     ENDDO
384     ENDIF
385     #endif /* ALLOW_OBCS_WEST */
386    
387     #ifdef ALLOW_OBCS_NORTH
388     C KS16. There was a bug here. Used to say OBCS_balanceFacW, not
389     C FacN!!!
390     IF ( OBCS_balanceFacN.NE.0. ) THEN
391     DO bj=myByLo(myThid),myByHi(myThid)
392     DO bi=myBxLo(myThid),myBxHi(myThid)
393     IF ( tileHasOBN(bi,bj) ) THEN
394     DO k=1,Nr
395     DO i=1-OLx,sNx+OLx
396     jB = OB_Jn(i,bi,bj)
397     IF ( jB.NE.OB_indexNone ) THEN
398     OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
399     & + flowN*maskS(i,jB,k,bi,bj)
400     ENDIF
401     ENDDO
402     ENDDO
403     ENDIF
404     ENDDO
405     ENDDO
406     ENDIF
407     #endif /* ALLOW_OBCS_NORTH */
408    
409     #ifdef ALLOW_OBCS_SOUTH
410     IF ( OBCS_balanceFacS.NE.0. ) THEN
411     DO bj=myByLo(myThid),myByHi(myThid)
412     DO bi=myBxLo(myThid),myBxHi(myThid)
413     IF ( tileHasOBS(bi,bj) ) THEN
414     DO k=1,Nr
415     DO i=1-OLx,sNx+OLx
416     jB = OB_Js(i,bi,bj)
417     IF ( jB.NE.OB_indexNone ) THEN
418     OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
419     & + flowS*maskS(i,1+jB,k,bi,bj)
420     ENDIF
421     ENDDO
422     ENDDO
423     ENDIF
424     ENDDO
425     ENDDO
426     ENDIF
427     #endif /* ALLOW_OBCS_SOUTH */
428    
429     #ifdef ALLOW_DEBUG
430     IF (debugMode) CALL DEBUG_LEAVE('OBCS_BALANCE_FLOW',myThid)
431     #endif
432    
433     #endif /* ALLOW_OBCS_BALANCE */
434     #endif /* ALLOW_OBCS */
435    
436     RETURN
437     END

  ViewVC Help
Powered by ViewVC 1.1.22