/[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.6 - (hide annotations) (download)
Mon May 14 22:02:33 2007 UTC (17 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59c, checkpoint59b, checkpoint59h
Changes since 1.5: +2 -2 lines
Cleanup suggested by M. Mazloff (remove _loc)

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

  ViewVC Help
Powered by ViewVC 1.1.22