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

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

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


Revision 1.7 - (hide annotations) (download)
Tue Oct 9 00:00:00 2007 UTC (16 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62q, checkpoint62p, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint60, checkpoint61, checkpoint62, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59k, checkpoint62b, checkpoint61f, checkpoint61n, checkpoint59j, checkpoint61q, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.6: +32 -30 lines
add missing cvs $Header:$ or $Name:$

1 jmc 1.7 C $Header: $
2     C $Name: $
3 heimbach 1.2
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    
70     cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
71    
72     logical doglobalread
73     logical ladinit
74    
75     character*(80) fnameobcse
76    
77     cgg( Variables for splitting barotropic/baroclinic vels.
78     _RL ubaro
79     _RL utop
80     cgg)
81    
82     c == external functions ==
83    
84     integer ilnblnk
85     external ilnblnk
86    
87    
88     c == end of interface ==
89    
90     jtlo = mybylo(mythid)
91     jthi = mybyhi(mythid)
92     itlo = mybxlo(mythid)
93     ithi = mybxhi(mythid)
94     jmin = 1-oly
95     jmax = sny+oly
96     imin = 1-olx
97     imax = snx+olx
98     ip1 = 0
99    
100     cgg( Initialize variables for balancing volume flux.
101     ubaro = 0.d0
102     utop = 0.d0
103     cgg)
104    
105     c-- Now, read the control vector.
106     doglobalread = .false.
107     ladinit = .false.
108    
109     if (optimcycle .ge. 0) then
110     ilobcse=ilnblnk( xx_obcse_file )
111 jmc 1.7 write(fnameobcse(1:80),'(2a,i10.10)')
112 heimbach 1.2 & xx_obcse_file(1:ilobcse), '.', optimcycle
113     endif
114    
115     c-- Get the counters, flags, and the interpolation factor.
116     call ctrl_get_gen_rec(
117     I xx_obcsestartdate, xx_obcseperiod,
118     O obcsefac, obcsefirst, obcsechanged,
119     O obcsecount0,obcsecount1,
120     I mytime, myiter, mythid )
121    
122     do iobcs = 1,nobcs
123    
124     if ( obcsefirst ) then
125 jmc 1.7 call active_read_yz( fnameobcse, tmpfldyz,
126 heimbach 1.2 & (obcsecount0-1)*nobcs+iobcs,
127     & doglobalread, ladinit, optimcycle,
128     & mythid, xx_obcse_dummy )
129    
130 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
131    
132 jmc 1.7 if ( optimcycle .gt. 0) then
133 heimbach 1.2 if (iobcs .eq. 3) then
134     cgg Special attention is needed for the normal velocity.
135     cgg For the north, this is the v velocity, iobcs = 4.
136     cgg This is done on a columnwise basis here.
137     do bj = jtlo,jthi
138     do bi = itlo, ithi
139     do j = jmin,jmax
140     i = OB_Ie(J,bi,bj)
141    
142     cgg The barotropic velocity is stored in the level 1.
143     ubaro = tmpfldyz(j,1,bi,bj)
144     tmpfldyz(j,1,bi,bj) = 0.d0
145     utop = 0.d0
146    
147     do k = 1,Nr
148     cgg If cells are not full, this should be modified with hFac.
149 jmc 1.7 cgg
150     cgg The xx field (tmpfldxz) does not contain the velocity at the
151     cgg surface level. This velocity is not independent; it must
152 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
153 jmc 1.7 cgg the baroclinic velocity structure..
154 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
155 heimbach 1.4 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
156 jmc 1.7 cgg Add the barotropic velocity component.
157 heimbach 1.2 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
158     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
159     endif
160     enddo
161     cgg Compute the baroclinic velocity at level 1. Should balance flux.
162 jmc 1.7 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
163 heimbach 1.4 & - utop / delR(1)
164 heimbach 1.2 enddo
165     enddo
166     enddo
167     endif
168     if (iobcs .eq. 4) then
169     cgg Special attention is needed for the normal velocity.
170     cgg For the north, this is the v velocity, iobcs = 4.
171     cgg This is done on a columnwise basis here.
172     do bj = jtlo,jthi
173     do bi = itlo, ithi
174     do j = jmin,jmax
175     i = OB_Ie(J,bi,bj)
176    
177     cgg The barotropic velocity is stored in the level 1.
178     ubaro = tmpfldyz(j,1,bi,bj)
179     tmpfldyz(j,1,bi,bj) = 0.d0
180     utop = 0.d0
181    
182     do k = 1,Nr
183     cgg If cells are not full, this should be modified with hFac.
184 jmc 1.7 cgg
185     cgg The xx field (tmpfldxz) does not contain the velocity at the
186     cgg surface level. This velocity is not independent; it must
187 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
188 jmc 1.7 cgg the baroclinic velocity structure..
189 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
190 heimbach 1.4 & maskS(i,j,k,bi,bj) * delR(k) + utop
191 jmc 1.7 cgg Add the barotropic velocity component.
192 heimbach 1.2 if (maskS(i,j,k,bi,bj) .ne. 0.) then
193     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
194     endif
195     enddo
196     cgg Compute the baroclinic velocity at level 1. Should balance flux.
197 jmc 1.7 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
198 heimbach 1.4 & - utop / delR(1)
199 heimbach 1.2 enddo
200     enddo
201     enddo
202     endif
203     endif
204    
205 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
206    
207 heimbach 1.2 do bj = jtlo,jthi
208     do bi = itlo,ithi
209     do k = 1,nr
210     do j = jmin,jmax
211     xx_obcse1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
212     cgg & * maskyz (j,k,bi,bj)
213     enddo
214     enddo
215     enddo
216     enddo
217     endif
218    
219     if ( (obcsefirst) .or. (obcsechanged)) then
220 jmc 1.7
221 heimbach 1.2 do bj = jtlo,jthi
222     do bi = itlo,ithi
223     do j = jmin,jmax
224     do k = 1,nr
225     tmpfldyz(j,k,bi,bj) = xx_obcse1(j,k,bi,bj,iobcs)
226     enddo
227     enddo
228     enddo
229     enddo
230    
231     call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
232    
233     do bj = jtlo,jthi
234     do bi = itlo,ithi
235     do j = jmin,jmax
236     do k = 1,nr
237     xx_obcse0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
238     enddo
239     enddo
240     enddo
241     enddo
242    
243 jmc 1.7 call active_read_yz( fnameobcse, tmpfldyz,
244 heimbach 1.2 & (obcsecount1-1)*nobcs+iobcs,
245     & doglobalread, ladinit, optimcycle,
246     & mythid, xx_obcse_dummy )
247    
248 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
249    
250 jmc 1.7 if ( optimcycle .gt. 0) then
251 heimbach 1.2 if (iobcs .eq. 3) then
252     cgg Special attention is needed for the normal velocity.
253     cgg For the north, this is the v velocity, iobcs = 4.
254     cgg This is done on a columnwise basis here.
255     do bj = jtlo,jthi
256     do bi = itlo, ithi
257     do j = jmin,jmax
258     i = OB_Ie(J,bi,bj)
259    
260     cgg The barotropic velocity is stored in the level 1.
261     ubaro = tmpfldyz(j,1,bi,bj)
262     tmpfldyz(j,1,bi,bj) = 0.d0
263     utop = 0.d0
264    
265     do k = 1,Nr
266     cgg If cells are not full, this should be modified with hFac.
267 jmc 1.7 cgg
268     cgg The xx field (tmpfldxz) does not contain the velocity at the
269     cgg surface level. This velocity is not independent; it must
270 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
271 jmc 1.7 cgg the baroclinic velocity structure..
272 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
273 heimbach 1.4 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
274 jmc 1.7 cgg Add the barotropic velocity component.
275 heimbach 1.2 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
276     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
277     endif
278     enddo
279     cgg Compute the baroclinic velocity at level 1. Should balance flux.
280 jmc 1.7 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
281 heimbach 1.4 & - utop / delR(1)
282 heimbach 1.2 enddo
283     enddo
284     enddo
285     endif
286     if (iobcs .eq. 4) then
287     cgg Special attention is needed for the normal velocity.
288     cgg For the north, this is the v velocity, iobcs = 4.
289     cgg This is done on a columnwise basis here.
290     do bj = jtlo,jthi
291     do bi = itlo, ithi
292     do j = jmin,jmax
293     i = OB_Ie(J,bi,bj)
294    
295     cgg The barotropic velocity is stored in the level 1.
296     ubaro = tmpfldyz(j,1,bi,bj)
297     tmpfldyz(j,1,bi,bj) = 0.d0
298     utop = 0.d0
299    
300     do k = 1,Nr
301     cgg If cells are not full, this should be modified with hFac.
302 jmc 1.7 cgg
303     cgg The xx field (tmpfldxz) does not contain the velocity at the
304     cgg surface level. This velocity is not independent; it must
305 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
306 jmc 1.7 cgg the baroclinic velocity structure..
307 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
308 heimbach 1.4 & maskS(i,j,k,bi,bj) * delR(k) + utop
309 jmc 1.7 cgg Add the barotropic velocity component.
310 heimbach 1.2 if (maskS(i,j,k,bi,bj) .ne. 0.) then
311     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
312     endif
313     enddo
314     cgg Compute the baroclinic velocity at level 1. Should balance flux.
315 jmc 1.7 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
316 heimbach 1.4 & - utop / delR(1)
317 heimbach 1.2 enddo
318     enddo
319     enddo
320     endif
321     endif
322    
323 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
324    
325 heimbach 1.2 do bj = jtlo,jthi
326     do bi = itlo,ithi
327     do k = 1,nr
328     do j = jmin,jmax
329     xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
330     cgg & * maskyz (j,k,bi,bj)
331     enddo
332     enddo
333     enddo
334     enddo
335     endif
336    
337     c-- Add control to model variable.
338     do bj = jtlo,jthi
339     do bi = itlo,ithi
340     c-- Calculate mask for tracer cells (0 => land, 1 => water).
341     do k = 1,nr
342     do j = 1,sny
343     i = OB_Ie(j,bi,bj)
344     if (iobcs .EQ. 1) then
345     OBEt(j,k,bi,bj) = OBEt (j,k,bi,bj)
346     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
347     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
348     OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
349     & *maskW(i+ip1,j,k,bi,bj)
350     else if (iobcs .EQ. 2) then
351     OBEs(j,k,bi,bj) = OBEs (j,k,bi,bj)
352     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
353     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
354     OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
355     & *maskW(i+ip1,j,k,bi,bj)
356     else if (iobcs .EQ. 3) then
357     OBEu(j,k,bi,bj) = OBEu (j,k,bi,bj)
358     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
359     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
360     OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
361     & *maskW(i+ip1,j,k,bi,bj)
362     else if (iobcs .EQ. 4) then
363     OBEv(j,k,bi,bj) = OBEv (j,k,bi,bj)
364     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
365     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
366     OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
367     & *maskS(i,j,k,bi,bj)
368     endif
369     enddo
370     enddo
371     enddo
372     enddo
373    
374     C-- End over iobcs loop
375     enddo
376    
377     #else /* ALLOW_OBCSE_CONTROL undefined */
378    
379     c == routine arguments ==
380    
381     _RL mytime
382     integer myiter
383     integer mythid
384    
385     c-- CPP flag ALLOW_OBCSE_CONTROL undefined.
386    
387     #endif /* ALLOW_OBCSE_CONTROL */
388    
389     end
390    

  ViewVC Help
Powered by ViewVC 1.1.22