/[MITgcm]/MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcsw.F
ViewVC logotype

Contents of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcsw.F

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


Revision 1.5 - (show annotations) (download)
Wed Jan 19 21:08:32 2011 UTC (14 years, 11 months ago) by mmazloff
Branch: MAIN
Changes since 1.4: +234 -114 lines
 Bringing code up to date and ensuring backward compatibility

1 C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_getobcsw.F,v 1.9 2011/01/19 08:42:06 mlosch Exp $
2 C $Name: $
3
4 #include "CTRL_CPPOPTIONS.h"
5 #ifdef ALLOW_OBCS
6 # include "OBCS_OPTIONS.h"
7 #endif
8
9
10 subroutine ctrl_getobcsw(
11 I mytime,
12 I myiter,
13 I mythid
14 & )
15
16 c ==================================================================
17 c SUBROUTINE ctrl_getobcsw
18 c ==================================================================
19 c
20 c o Get western obc of the control vector and add it
21 c to dyn. fields
22 c
23 c started: heimbach@mit.edu, 29-Aug-2001
24 c
25 c modified: gebbie@mit.edu, 18-Mar-2003
26 c ==================================================================
27 c SUBROUTINE ctrl_getobcsw
28 c ==================================================================
29
30 implicit none
31
32 #ifdef ALLOW_OBCSW_CONTROL
33
34 c == global variables ==
35
36 #include "EEPARAMS.h"
37 #include "SIZE.h"
38 #include "PARAMS.h"
39 #include "GRID.h"
40 #include "OBCS.h"
41
42 #include "ctrl.h"
43 #include "ctrl_dummy.h"
44 #include "optim.h"
45
46 c == routine arguments ==
47
48 _RL mytime
49 integer myiter
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 ilobcsw
61 integer iobcs
62 integer ip1
63
64 _RL dummy
65 _RL obcswfac
66 logical obcswfirst
67 logical obcswchanged
68 integer obcswcount0
69 integer obcswcount1
70 integer nk,nz
71 _RL tmpz (nr,nsx,nsy)
72 _RL stmp
73
74 cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
75
76 logical doglobalread
77 logical ladinit
78
79 character*(80) fnameobcsw
80
81 cgg( Variables for splitting barotropic/baroclinic vels.
82 _RL ubaro
83 _RL utop
84 cgg)
85
86 c == external functions ==
87
88 integer ilnblnk
89 external ilnblnk
90
91
92 c == end of interface ==
93
94 jtlo = mybylo(mythid)
95 jthi = mybyhi(mythid)
96 itlo = mybxlo(mythid)
97 ithi = mybxhi(mythid)
98 jmin = 1-oly
99 jmax = sny+oly
100 imin = 1-olx
101 imax = snx+olx
102 ip1 = 1
103
104 cgg( Initialize variables for balancing volume flux.
105 ubaro = 0.d0
106 utop = 0.d0
107 cgg)
108
109 c-- Now, read the control vector.
110 doglobalread = .false.
111 ladinit = .false.
112
113 if (optimcycle .ge. 0) then
114 ilobcsw=ilnblnk( xx_obcsw_file )
115 write(fnameobcsw(1:80),'(2a,i10.10)')
116 & xx_obcsw_file(1:ilobcsw), '.', optimcycle
117 endif
118
119 c-- Get the counters, flags, and the interpolation factor.
120 call ctrl_get_gen_rec(
121 I xx_obcswstartdate, xx_obcswperiod,
122 O obcswfac, obcswfirst, obcswchanged,
123 O obcswcount0,obcswcount1,
124 I mytime, myiter, mythid )
125
126 do iobcs = 1,nobcs
127 if ( obcswfirst ) then
128 call active_read_yz( fnameobcsw, tmpfldyz,
129 & (obcswcount0-1)*nobcs+iobcs,
130 & doglobalread, ladinit, optimcycle,
131 & mythid, xx_obcsw_dummy )
132
133
134 if ( optimcycle .gt. 0) then
135 if (iobcs .eq. 3) then
136 cgg Special attention is needed for the normal velocity.
137 cgg For the north, this is the v velocity, iobcs = 4.
138 cgg This is done on a columnwise basis here.
139 do bj = jtlo,jthi
140 do bi = itlo, ithi
141 do j = jmin,jmax
142 cih If open boundary.
143 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
144 i = OB_Iw(J,bi,bj)
145 #ifdef ALLOW_OBCS_CONTROL_MODES
146 cih Determine number of open vertical layers.
147 nz = 0
148 do k = 1,Nr
149 nz = nz + maskW(i+ip1,j,k,bi,bj)
150 end do
151 cih Compute absolute velocities from the barotropic-baroclinic modes.
152 #ifdef ALLOW_CTRL_OBCS_BALANCE
153 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are
154 c compatible - to ensure volume conservation can just set barotropic
155 c mode amplitude to 0
156 c however this means no inflow at every horizontal location....
157 do k = 1,nr
158 tmpfldyz(k,1,bi,bj)= 0.
159 end do
160 #endif
161 do k = 1,Nr
162 if (k.le.nz) then
163 stmp = 0.
164 do nk = 1,nz
165 stmp = stmp +
166 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
167 end do
168 tmpz(k,bi,bj) = stmp
169 else
170 tmpz(k,bi,bj) = 0.
171 end if
172 end do
173 do k = 1,Nr
174 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
175 & *recip_hFacW(i+ip1,j,k,bi,bj)
176 end do
177 #elif defined (ALLOW_CTRL_OBCS_BALANCE)
178 cgg The barotropic velocity is stored in the level 1.
179 ubaro = tmpfldyz(j,1,bi,bj)
180 tmpfldyz(j,1,bi,bj) = 0.d0
181 utop = 0.d0
182
183 do k = 1,Nr
184 cgg If cells are not full, this should be modified with hFac.
185 cgg
186 cgg The xx field (tmpfldxz) does not contain the velocity at the
187 cgg surface level. This velocity is not independent; it must
188 cgg exactly balance the volume flux, since we are dealing with
189 cgg the baroclinic velocity structure..
190 utop = tmpfldyz(j,k,bi,bj)*
191 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
192 cgg Add the barotropic velocity component.
193 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
194 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
195 endif
196 enddo
197 cgg Compute the baroclinic velocity at level 1. Should balance flux.
198 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
199 & - utop / delR(1)
200 #endif
201 cih End open boundary.
202 end if
203 enddo
204 enddo
205 enddo
206 endif
207 if (iobcs .eq. 4) then
208 cgg Special attention is needed for the normal velocity.
209 cgg For the north, this is the v velocity, iobcs = 4.
210 cgg This is done on a columnwise basis here.
211 do bj = jtlo,jthi
212 do bi = itlo, ithi
213 do j = jmin,jmax
214 cih If open boundary.
215 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
216 i = OB_Iw(J,bi,bj)
217 #ifdef ALLOW_OBCS_CONTROL_MODES
218 cih Determine number of open vertical layers.
219 nz = 0
220 do k = 1,Nr
221 nz = nz + maskS(i,j,k,bi,bj)
222 end do
223 cih Compute absolute velocities from the barotropic-baroclinic modes.
224 #ifdef ALLOW_CTRL_OBCS_BALANCE
225 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are
226 c compatible - to ensure volume conservation can just set barotropic
227 c mode amplitude to 0
228 c however this means no inflow at every horizontal location....
229 do k = 1,nr
230 tmpfldyz(k,1,bi,bj)= 0.
231 end do
232 #endif
233 do k = 1,Nr
234 if (k.le.nz) then
235 stmp = 0.
236 do nk = 1,nz
237 stmp = stmp +
238 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
239 end do
240 tmpz(k,bi,bj) = stmp
241 else
242 tmpz(k,bi,bj) = 0.
243 end if
244 end do
245 do k = 1,Nr
246 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
247 & *recip_hFacS(i,j,k,bi,bj)
248 end do
249 #elif defined (ALLOW_CTRL_OBCS_BALANCE)
250 cgg The barotropic velocity is stored in the level 1.
251 ubaro = tmpfldyz(j,1,bi,bj)
252 tmpfldyz(j,1,bi,bj) = 0.d0
253 utop = 0.d0
254
255 do k = 1,Nr
256 cgg If cells are not full, this should be modified with hFac.
257 cgg
258 cgg The xx field (tmpfldxz) does not contain the velocity at the
259 cgg surface level. This velocity is not independent; it must
260 cgg exactly balance the volume flux, since we are dealing with
261 cgg the baroclinic velocity structure..
262 utop = tmpfldyz(j,k,bi,bj)*
263 & maskS(i,j,k,bi,bj) * delR(k) + utop
264 cgg Add the barotropic velocity component.
265 if (maskS(i,j,k,bi,bj) .ne. 0.) then
266 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
267 endif
268 enddo
269 cgg Compute the baroclinic velocity at level 1. Should balance flux.
270 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
271 & - utop / delR(1)
272 #endif
273 cih End if iobcs = 4.
274 end if
275 enddo
276 enddo
277 enddo
278 endif
279 endif
280
281 do bj = jtlo,jthi
282 do bi = itlo,ithi
283 do k = 1,nr
284 do j = jmin,jmax
285 xx_obcsw1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
286 cgg & * maskyz (j,k,bi,bj)
287 enddo
288 enddo
289 enddo
290 enddo
291 endif
292
293 if ( (obcswfirst) .or. (obcswchanged)) then
294
295 do bj = jtlo,jthi
296 do bi = itlo,ithi
297 do k = 1,nr
298 do j = jmin,jmax
299 xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
300 tmpfldyz (j,k,bi,bj) = 0. _d 0
301 enddo
302 enddo
303 enddo
304 enddo
305
306 call active_read_yz( fnameobcsw, tmpfldyz,
307 & (obcswcount1-1)*nobcs+iobcs,
308 & doglobalread, ladinit, optimcycle,
309 & mythid, xx_obcsw_dummy )
310
311 if ( optimcycle .gt. 0) then
312 if (iobcs .eq. 3) then
313 cgg Special attention is needed for the normal velocity.
314 cgg For the north, this is the v velocity, iobcs = 4.
315 cgg This is done on a columnwise basis here.
316 do bj = jtlo,jthi
317 do bi = itlo, ithi
318 do j = jmin,jmax
319 cih If open boundary.
320 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
321 i = OB_Iw(J,bi,bj)
322 #ifdef ALLOW_OBCS_CONTROL_MODES
323 cih Determine number of open vertical layers.
324 nz = 0
325 do k = 1,Nr
326 nz = nz + maskW(i+ip1,j,k,bi,bj)
327 end do
328 cih Compute absolute velocities from the barotropic-baroclinic modes.
329 #ifdef ALLOW_CTRL_OBCS_BALANCE
330 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are
331 c compatible - to ensure volume conservation can just set barotropic
332 c mode amplitude to 0
333 c however this means no inflow at every horizontal location....
334 do k = 1,nr
335 tmpfldyz(k,1,bi,bj)= 0.
336 end do
337 #endif
338 do k = 1,Nr
339 if (k.le.nz) then
340 stmp = 0.
341 do nk = 1,nz
342 stmp = stmp +
343 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
344 end do
345 tmpz(k,bi,bj) = stmp
346 else
347 tmpz(k,bi,bj) = 0.
348 end if
349 end do
350 do k = 1,Nr
351 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
352 & *recip_hFacW(i+ip1,j,k,bi,bj)
353 end do
354 #elif defined (ALLOW_CTRL_OBCS_BALANCE)
355 cgg The barotropic velocity is stored in the level 1.
356 ubaro = tmpfldyz(j,1,bi,bj)
357 tmpfldyz(j,1,bi,bj) = 0.d0
358 utop = 0.d0
359
360 do k = 1,Nr
361 cgg If cells are not full, this should be modified with hFac.
362 cgg
363 cgg The xx field (tmpfldxz) does not contain the velocity at the
364 cgg surface level. This velocity is not independent; it must
365 cgg exactly balance the volume flux, since we are dealing with
366 cgg the baroclinic velocity structure..
367 utop = tmpfldyz(j,k,bi,bj)*
368 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
369 cgg Add the barotropic velocity component.
370 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
371 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
372 endif
373 enddo
374 cgg Compute the baroclinic velocity at level 1. Should balance flux.
375 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
376 & - utop / delR(1)
377 #endif
378 cih End open boundary.
379 end if
380 enddo
381 enddo
382 enddo
383 endif
384 if (iobcs .eq. 4) then
385 cgg Special attention is needed for the normal velocity.
386 cgg For the north, this is the v velocity, iobcs = 4.
387 cgg This is done on a columnwise basis here.
388 do bj = jtlo,jthi
389 do bi = itlo, ithi
390 do j = jmin,jmax
391 cih If open boundary.
392 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
393 i = OB_Iw(J,bi,bj)
394 #ifdef ALLOW_OBCS_CONTROL_MODES
395 cih Determine number of open vertical layers.
396 nz = 0
397 do k = 1,Nr
398 nz = nz + maskS(i,j,k,bi,bj)
399 end do
400 cih Compute absolute velocities from the barotropic-baroclinic modes.
401 #ifdef ALLOW_CTRL_OBCS_BALANCE
402 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are
403 c compatible - to ensure volume conservation can just set barotropic
404 c mode amplitude to 0
405 c however this means no inflow at every horizontal location....
406 do k = 1,nr
407 tmpfldyz(k,1,bi,bj)= 0.
408 end do
409 #endif
410 do k = 1,Nr
411 if (k.le.nz) then
412 stmp = 0.
413 do nk = 1,nz
414 stmp = stmp +
415 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
416 end do
417 tmpz(k,bi,bj) = stmp
418 else
419 tmpz(k,bi,bj) = 0.
420 end if
421 end do
422 do k = 1,Nr
423 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
424 & *recip_hFacS(i,j,k,bi,bj)
425 end do
426 #elif defined (ALLOW_CTRL_OBCS_BALANCE)
427 cgg The barotropic velocity is stored in the level 1.
428 ubaro = tmpfldyz(j,1,bi,bj)
429 tmpfldyz(j,1,bi,bj) = 0.d0
430 utop = 0.d0
431
432 do k = 1,Nr
433 cgg If cells are not full, this should be modified with hFac.
434 cgg
435 cgg The xx field (tmpfldxz) does not contain the velocity at the
436 cgg surface level. This velocity is not independent; it must
437 cgg exactly balance the volume flux, since we are dealing with
438 cgg the baroclinic velocity structure..
439 utop = tmpfldyz(j,k,bi,bj)*
440 & maskS(i,j,k,bi,bj) * delR(k) + utop
441 cgg Add the barotropic velocity component.
442 if (maskS(i,j,k,bi,bj) .ne. 0.) then
443 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
444 endif
445 enddo
446 cgg Compute the baroclinic velocity at level 1. Should balance flux.
447 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
448 & - utop / delR(1)
449 #endif
450 cih End if iobcs = 4.
451 end if
452 enddo
453 enddo
454 enddo
455 endif
456 endif
457
458 do bj = jtlo,jthi
459 do bi = itlo,ithi
460 do k = 1,nr
461 do j = jmin,jmax
462 xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
463 cgg & * maskyz (j,k,bi,bj)
464 enddo
465 enddo
466 enddo
467 enddo
468 endif
469
470 c-- Add control to model variable.
471 do bj = jtlo, jthi
472 do bi = itlo, ithi
473 c-- Calculate mask for tracer cells (0 => land, 1 => water).
474 do k = 1,nr
475 do j = 1,sny
476 i = OB_Iw(j,bi,bj)
477 if (iobcs .EQ. 1) then
478 OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
479 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
480 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
481 OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
482 & *maskW(i+ip1,j,k,bi,bj)
483 else if (iobcs .EQ. 2) then
484 OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
485 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
486 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
487 OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
488 & *maskW(i+ip1,j,k,bi,bj)
489 else if (iobcs .EQ. 3) then
490 OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
491 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
492 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
493 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
494 & *maskW(i+ip1,j,k,bi,bj)
495 else if (iobcs .EQ. 4) then
496 OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
497 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
498 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
499 OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
500 & *maskS(i,j,k,bi,bj)
501 endif
502 enddo
503 enddo
504 enddo
505 enddo
506
507 C-- End over iobcs loop
508 enddo
509
510 #else /* ALLOW_OBCSW_CONTROL undefined */
511
512 c == routine arguments ==
513
514 _RL mytime
515 integer myiter
516 integer mythid
517
518 c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
519
520 #endif /* ALLOW_OBCSW_CONTROL */
521
522 end
523

  ViewVC Help
Powered by ViewVC 1.1.22