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

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

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


Revision 1.2 - (show annotations) (download)
Tue Jun 24 16:07:06 2003 UTC (20 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51a_post, checkpoint51c_post, checkpoint51, checkpoint51b_post, checkpoint51b_pre
Changes since 1.1: +495 -0 lines
Merging for c51 vs. e34

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

  ViewVC Help
Powered by ViewVC 1.1.22