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

Annotation of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcse.F

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


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

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

  ViewVC Help
Powered by ViewVC 1.1.22