/[MITgcm]/MITgcm/pkg/ctrl/ctrl_obcsvol.F
ViewVC logotype

Annotation of /MITgcm/pkg/ctrl/ctrl_obcsvol.F

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


Revision 1.9 - (hide annotations) (download)
Mon May 14 22:02:34 2007 UTC (17 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59c, checkpoint59b, checkpoint59h
Changes since 1.8: +8 -8 lines
Cleanup suggested by M. Mazloff (remove _loc)

1 heimbach 1.2
2     #include "CTRL_CPPOPTIONS.h"
3     #ifdef ALLOW_OBCS
4     # include "OBCS_OPTIONS.h"
5     #endif
6    
7     subroutine ctrl_obcsvol(
8     I mytime,
9     I myiter,
10     I mythid
11     & )
12    
13     c ==================================================================
14     c SUBROUTINE ctrl_obcsvol
15     c ==================================================================
16     c
17     c o volumetrically balance the control vector contribution.
18     c o Assume the calendar is identical
19     c for all open boundaries. Need to save the barotropic adjustment
20     c velocity so it can be used in all ctrl_getobcs files.
21     c ==================================================================
22     c SUBROUTINE ctrl_obcsvol
23     c ==================================================================
24    
25     implicit none
26    
27     c == global variables ==
28    
29     #include "EEPARAMS.h"
30     #include "SIZE.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33     #include "DYNVARS.h"
34     #ifdef ALLOW_OBCS
35     # include "OBCS.h"
36     #endif
37    
38     #include "ctrl.h"
39     #include "ctrl_dummy.h"
40     #include "optim.h"
41    
42     c == routine arguments ==
43    
44     integer myiter
45     _RL mytime
46     integer mythid
47    
48 mlosch 1.8 #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL
49 heimbach 1.2 c == local variables ==
50    
51     integer bi,bj
52     integer i,j,k
53     integer itlo,ithi
54     integer jtlo,jthi
55     integer jmin,jmax
56     integer imin,imax
57     integer irec
58     integer il
59     integer iobcs
60     integer ip1
61     integer jp1
62     integer nrec
63     integer ilfld
64     integer igg
65    
66     _RL sumvol
67     _RL sumarea
68     _RL tmpflux
69     _RL tmparea
70     _RL dummy
71     _RL gg
72     _RL tmpx
73     _RL tmpy
74     character*(80) fnamefldn
75     character*(80) fnameflds
76     character*(80) fnamefldw
77     character*(80) fnameflde
78    
79     logical doglobalread
80     logical ladinit
81     logical obcsnfirst, obcsnchanged
82     integer obcsncount0, obcsncount1
83     _RL obcsnfac
84    
85     #ifdef ECCO_VERBOSE
86     character*(MAX_LEN_MBUF) msgbuf
87     #endif
88    
89     c == external functions ==
90    
91     integer ilnblnk
92     external ilnblnk
93    
94     c == end of interface ==
95    
96     jtlo = mybylo(mythid)
97     jthi = mybyhi(mythid)
98     itlo = mybxlo(mythid)
99     ithi = mybxhi(mythid)
100     jmin = 1
101     jmax = sny
102     imin = 1
103     imax = snx
104    
105     c-- Read tiled data.
106     doglobalread = .false.
107     ladinit = .false.
108    
109     cgg Assume the number of records is the same for
110     cgg all boundaries.
111    
112     tmpflux= 0. d 0
113     tmparea= 0. d 0
114     sumarea= 0. d 0
115     sumvol = 0. d 0
116    
117     #ifdef ECCO_VERBOSE
118     _BEGIN_MASTER( mythid )
119     write(msgbuf,'(a)') ' '
120     call print_message( msgbuf, standardmessageunit,
121     & SQUEEZE_RIGHT , mythid)
122     write(msgbuf,'(a)') ' '
123     call print_message( msgbuf, standardmessageunit,
124     & SQUEEZE_RIGHT , mythid)
125     write(msgbuf,'(a,i9.8)')
126     & ' ctrl_obcsvol: number of records to process: ',nrec
127     call print_message( msgbuf, standardmessageunit,
128     & SQUEEZE_RIGHT , mythid)
129     write(msgbuf,'(a)') ' '
130     call print_message( msgbuf, standardmessageunit,
131     & SQUEEZE_RIGHT , mythid)
132     _END_MASTER( mythid )
133     #endif
134    
135     c-- Get the counters, flags, and the interpolation factor.
136     call ctrl_GetRec( 'xx_obcsn',
137     O obcsnfac, obcsnfirst, obcsnchanged,
138     O obcsncount0,obcsncount1,
139     I mytime, myiter, mythid )
140    
141     if (optimcycle .ge. 0) then
142     c
143     ilfld=ilnblnk( xx_obcsn_file )
144     write(fnamefldn(1:80),'(2a,i10.10)')
145     & xx_obcsn_file(1:ilfld),'.', optimcycle
146     ilfld=ilnblnk( xx_obcss_file )
147     write(fnameflds(1:80),'(2a,i10.10)')
148     & xx_obcss_file(1:ilfld),'.',optimcycle
149     ilfld=ilnblnk( xx_obcsw_file )
150     write(fnamefldw(1:80),'(2a,i10.10)')
151     & xx_obcsw_file(1:ilfld),'.',optimcycle
152     ilfld=ilnblnk( xx_obcse_file )
153     write(fnameflde(1:80),'(2a,i10.10)')
154     & xx_obcse_file(1:ilfld),'.',optimcycle
155     c
156     endif
157    
158     c-- Loop over records. For north boundary, we only need V velocity.
159    
160     if ( obcsnfirst ) then
161    
162     shiftvel(1) = 0. d0
163     shiftvel(2) = 0. d0
164    
165     #ifdef ALLOW_OBCSN_CONTROL
166     jp1 = 0
167    
168 heimbach 1.9 call active_read_xz(fnamefldn,tmpfldxz,
169 heimbach 1.3 & (obcsncount0-1)*nobcs+4, doglobalread,
170 heimbach 1.2 & ladinit, optimcycle, mythid
171     & , xx_obcsn_dummy )
172    
173     c-- Loop over this thread's tiles.
174     do bj = jtlo,jthi
175     do bi = itlo,ithi
176    
177     tmpflux = 0. d0
178     tmparea = 0. d0
179    
180     do k = 1, Nr
181     do i = imin,imax
182     j = Ob_Jn(I,bi,bj)
183     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
184 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
185 heimbach 1.2 if(ob_iw(j,bi,bj).ne.i .and. ob_ie(j,bi,bj).ne.i) then
186     tmpx = tmpfldxz(i,k,bi,bj)
187     cgg -- Positive is flux in.
188 heimbach 1.5 tmpflux = tmpflux - tmpx* delR(k)*dxg(i,j+jp1,bi,bj)
189     tmparea = tmparea + delR(k) * dxg(i,j+jp1,bi,bj)
190 heimbach 1.2 endif
191     endif
192     enddo
193     enddo
194    
195     sumarea = sumarea+ tmparea
196     sumvol = sumvol + tmpflux
197     enddo
198     enddo
199     #endif
200    
201     #ifdef ALLOW_OBCSS_CONTROL
202     jp1 = 1
203    
204 heimbach 1.9 call active_read_xz(fnameflds,tmpfldxz,
205 heimbach 1.3 & (obcsncount0-1)*nobcs+4, doglobalread,
206 heimbach 1.2 & ladinit, optimcycle, mythid
207     & , xx_obcss_dummy )
208    
209     c-- Loop over this thread's tiles.
210     do bj = jtlo,jthi
211     do bi = itlo,ithi
212    
213     tmpflux = 0. d 0
214     tmparea = 0. d 0
215    
216     do k = 1, Nr
217     do i = imin,imax
218     j = Ob_Js(I,bi,bj)
219     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
220 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
221 heimbach 1.2 if (ob_iw(j,bi,bj).ne.i .and.ob_ie(j,bi,bj).ne.i) then
222     tmpx = tmpfldxz(i,k,bi,bj)
223     cgg -- Positive is flux in.
224 heimbach 1.5 tmpflux = tmpflux + tmpx* delR(k)*dxg(i,j+jp1,bi,bj)
225     tmparea = tmparea + delR(k) * dxg(i,j+jp1,bi,bj)
226 heimbach 1.2 endif
227     endif
228     enddo
229     enddo
230     sumarea = sumarea+ tmparea
231     sumvol = sumvol + tmpflux
232     enddo
233     enddo
234     #endif
235    
236     #ifdef ALLOW_OBCSW_CONTROL
237     ip1 = 1
238    
239 heimbach 1.9 call active_read_yz( fnamefldw, tmpfldyz,
240 heimbach 1.2 & (obcsncount0-1)*nobcs+3, doglobalread,
241     & ladinit, optimcycle, mythid
242     & , xx_obcsw_dummy )
243    
244     c-- Loop over this thread's tiles.
245     do bj = jtlo,jthi
246     do bi = itlo,ithi
247    
248     c-- Determine the weights to be used.
249     tmpflux = 0. d 0
250     tmparea = 0. d 0
251    
252     do k = 1, Nr
253     do j = jmin,jmax
254     i = ob_iw(j,bi,bj)
255     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
256 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
257 heimbach 1.2 if (ob_jn(i,bi,bj).ne.j .and. ob_js(i,bi,bj).ne.j)then
258     tmpy = tmpfldyz(j,k,bi,bj)
259     cgg -- Positive is flux in.
260 heimbach 1.5 tmpflux = tmpflux + tmpy* delR(k)*dyg(i+ip1,j,bi,bj)
261     tmparea = tmparea + delR(k)*dyg(i+ip1,j,bi,bj)
262 heimbach 1.2 endif
263     endif
264     enddo
265     enddo
266     sumarea =sumarea + tmparea
267     sumvol = sumvol + tmpflux
268     enddo
269     enddo
270     #endif
271    
272     #ifdef ALLOW_OBCSE_CONTROL
273     ip1 = 0
274    
275 heimbach 1.9 call active_read_yz( fnameflde, tmpfldyz,
276 heimbach 1.2 (obcsncount0-1)*nobcs+3, doglobalread,
277     ladinit, optimcycle, mythid
278     & , xx_obcse_dummy )
279    
280     c-- Loop over this thread's tiles.
281     do bj = jtlo,jthi
282     do bi = itlo,ithi
283    
284     c-- Determine the weights to be used.
285     tmpflux = 0. d 0
286     tmparea = 0. d 0
287    
288     do k = 1, Nr
289     do j = jmin,jmax
290     i = ob_ie(j,bi,bj)
291     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
292 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
293 heimbach 1.2 if (ob_jn(i,bi,bj).ne.j .and. ob_js(i,bi,bj).ne.j)then
294     tmpy = tmpfldyz(j,k,bi,bj)
295     cgg -- Positive is flux in.
296 heimbach 1.5 tmpflux = tmpflux - tmpy* delR(k)*dyg(i+ip1,j,bi,bj)
297     tmparea = tmparea + delR(k) *dyg(i+ip1,j,bi,bj)
298 heimbach 1.2 endif
299     endif
300     enddo
301     enddo
302     sumarea = sumarea+ tmparea
303     sumvol = sumvol + tmpflux
304     enddo
305     enddo
306     #endif
307    
308     c-- Do the global summation.
309     _GLOBAL_SUM_R8( sumvol, mythid )
310     _GLOBAL_SUM_R8( sumarea,mythid )
311    
312     shiftvel(2) = sumvol /sumarea
313     endif
314     cgg End of the obcsnfirst loop.
315    
316     if ( ( obcsnfirst) .or. (obcsnchanged)) then
317    
318     cgg Swap the value.
319     shiftvel(1) = shiftvel(2)
320    
321     sumvol = 0. d0
322     sumarea= 0. d0
323    
324     #ifdef ALLOW_OBCSN_CONTROL
325     jp1 = 0
326    
327 heimbach 1.9 call active_read_xz(fnamefldn,tmpfldxz,
328 heimbach 1.3 & (obcsncount1-1)*nobcs+4, doglobalread,
329 heimbach 1.2 & ladinit, optimcycle, mythid
330     & , xx_obcsn_dummy )
331     c-- Loop over this thread's tiles.
332     do bj = jtlo,jthi
333     do bi = itlo,ithi
334    
335     tmpflux = 0. d0
336     tmparea = 0. d0
337    
338     do k = 1, Nr
339     do i = imin,imax
340     j = Ob_Jn(I,bi,bj)
341     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
342 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
343 heimbach 1.2 if (ob_iw(j,bi,bj).ne.i .and. ob_ie(j,bi,bj).ne.i)then
344     tmpx = tmpfldxz(i,k,bi,bj)
345     cgg -- Positive is flux in.
346 heimbach 1.5 tmpflux = tmpflux - tmpx* delR(k)*dxg(i,j+jp1,bi,bj)
347     tmparea = tmparea + delR(k) * dxg(i,j+jp1,bi,bj)
348 heimbach 1.2 endif
349     endif
350     enddo
351     enddo
352    
353     sumarea = sumarea+ tmparea
354     sumvol = sumvol + tmpflux
355     enddo
356     enddo
357     #endif
358    
359     #ifdef ALLOW_OBCSS_CONTROL
360     jp1 = 1
361    
362 heimbach 1.9 call active_read_xz(fnameflds,tmpfldxz,
363 heimbach 1.3 & (obcsncount1-1)*nobcs+4, doglobalread,
364 heimbach 1.2 & ladinit, optimcycle, mythid
365     & , xx_obcss_dummy )
366    
367     c-- Loop over this thread's tiles.
368     do bj = jtlo,jthi
369     do bi = itlo,ithi
370    
371     tmpflux = 0. d 0
372     tmparea = 0. d 0
373    
374     do k = 1, Nr
375     do i = imin,imax
376     j = Ob_Js(I,bi,bj)
377     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
378 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
379 heimbach 1.2 if (ob_iw(j,bi,bj).ne.i .and. ob_ie(j,bi,bj).ne.i)then
380     tmpx = tmpfldxz(i,k,bi,bj)
381     cgg -- Positive is flux in.
382 heimbach 1.5 tmpflux = tmpflux + tmpx* delR(k)*dxg(i,j+jp1,bi,bj)
383     tmparea = tmparea + delR(k) * dxg(i,j+jp1,bi,bj)
384 heimbach 1.2 endif
385     endif
386     enddo
387     enddo
388     sumarea = sumarea+ tmparea
389     sumvol = sumvol + tmpflux
390     enddo
391     enddo
392     #endif
393    
394     #ifdef ALLOW_OBCSW_CONTROL
395     ip1 = 1
396    
397 heimbach 1.9 call active_read_yz( fnamefldw, tmpfldyz,
398 heimbach 1.2 & (obcsncount1-1)*nobcs+3, doglobalread,
399     & ladinit, optimcycle, mythid
400     & , xx_obcsw_dummy )
401    
402     c-- Loop over this thread's tiles.
403     do bj = jtlo,jthi
404     do bi = itlo,ithi
405    
406     c-- Determine the weights to be used.
407     tmpflux = 0. d 0
408     tmparea = 0. d 0
409    
410     do k = 1, Nr
411     do j = jmin,jmax
412     i = ob_iw(j,bi,bj)
413     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
414 edhill 1.7 cgg -- Do not let corners contribute.
415 heimbach 1.2 if (ob_jn(i,bi,bj).ne.j .and. ob_js(i,bi,bj).ne.j)then
416     tmpy = tmpfldyz(j,k,bi,bj)
417     cgg -- Positive is flux in.
418 heimbach 1.5 tmpflux = tmpflux + tmpy* delR(k)*dyg(i+ip1,j,bi,bj)
419     tmparea = tmparea + delR(k)*dyg(i+ip1,j,bi,bj)
420 heimbach 1.2 endif
421     endif
422     enddo
423     enddo
424     sumarea =sumarea + tmparea
425     sumvol = sumvol + tmpflux
426     enddo
427     enddo
428     #endif
429    
430     #ifdef ALLOW_OBCSE_CONTROL
431     ip1 = 0
432    
433 heimbach 1.9 call active_read_yz( fnameflde, tmpfldyz,
434 heimbach 1.2 (obcsncount1-1)*nobcs+3, doglobalread,
435     ladinit, optimcycle, mythid
436     & , xx_obcse_dummy )
437    
438     c-- Loop over this thread's tiles.
439     do bj = jtlo,jthi
440     do bi = itlo,ithi
441    
442     c-- Determine the weights to be used.
443     tmpflux = 0. d 0
444     tmparea = 0. d 0
445    
446     do k = 1, Nr
447     do j = jmin,jmax
448     i = ob_ie(j,bi,bj)
449     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
450 edhill 1.4 cgg -- Do not let the corners contribute to the volume flux.
451 heimbach 1.2 if (ob_jn(i,bi,bj).ne.j .and. ob_js(i,bi,bj).ne.j)then
452     tmpy = tmpfldyz(j,k,bi,bj)
453     cgg -- Positive is flux in.
454 heimbach 1.5 tmpflux = tmpflux - tmpy* delR(k)*dyg(i+ip1,j,bi,bj)
455     tmparea = tmparea + delR(k) *dyg(i+ip1,j,bi,bj)
456 heimbach 1.2 endif
457     endif
458     enddo
459     enddo
460     sumarea = sumarea+ tmparea
461     sumvol = sumvol + tmpflux
462     enddo
463     enddo
464     #endif
465    
466     c-- Do the global summation.
467     _GLOBAL_SUM_R8( sumvol, mythid )
468     _GLOBAL_SUM_R8( sumarea,mythid )
469    
470     shiftvel(2) = sumvol /sumarea
471     endif
472     cgg End of the obcsnfirst, obcsnchanged loop.
473    
474 mlosch 1.8 #endif /* BALANCE_CONTROL_VOLFLUX_GLOBAL */
475 heimbach 1.2
476     return
477     end
478    
479    
480    
481    
482    
483    
484    

  ViewVC Help
Powered by ViewVC 1.1.22