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

Contents of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcss.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: +261 -140 lines
 Bringing code up to date and ensuring backward compatibility

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

  ViewVC Help
Powered by ViewVC 1.1.22