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

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

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


Revision 1.7 - (hide annotations) (download)
Mon May 14 22:02:33 2007 UTC (18 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59c, checkpoint59b, checkpoint59h
Changes since 1.6: +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_getobcsw(
9     I mytime,
10     I myiter,
11     I mythid
12     & )
13    
14     c ==================================================================
15     c SUBROUTINE ctrl_getobcsw
16     c ==================================================================
17     c
18     c o Get western 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 modified: gebbie@mit.edu, 18-Mar-2003
24     c ==================================================================
25     c SUBROUTINE ctrl_getobcsw
26     c ==================================================================
27    
28     implicit none
29    
30     #ifdef ALLOW_OBCSW_CONTROL
31    
32     c == global variables ==
33    
34     #include "EEPARAMS.h"
35     #include "SIZE.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38     #include "OBCS.h"
39    
40     #include "ctrl.h"
41     #include "ctrl_dummy.h"
42     #include "optim.h"
43    
44     c == routine arguments ==
45    
46     _RL mytime
47     integer myiter
48     integer mythid
49    
50     c == local variables ==
51    
52     integer bi,bj
53     integer i,j,k
54     integer itlo,ithi
55     integer jtlo,jthi
56     integer jmin,jmax
57     integer imin,imax
58     integer ilobcsw
59     integer iobcs
60     integer ip1
61    
62     _RL dummy
63     _RL obcswfac
64     logical obcswfirst
65     logical obcswchanged
66     integer obcswcount0
67     integer obcswcount1
68    
69     cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
70    
71     logical doglobalread
72     logical ladinit
73    
74     character*(80) fnameobcsw
75    
76     cgg( Variables for splitting barotropic/baroclinic vels.
77     _RL ubaro
78     _RL utop
79     cgg)
80    
81     c == external functions ==
82    
83     integer ilnblnk
84     external ilnblnk
85    
86    
87     c == end of interface ==
88    
89     jtlo = mybylo(mythid)
90     jthi = mybyhi(mythid)
91     itlo = mybxlo(mythid)
92     ithi = mybxhi(mythid)
93     jmin = 1-oly
94     jmax = sny+oly
95     imin = 1-olx
96     imax = snx+olx
97     ip1 = 1
98    
99     cgg( Initialize variables for balancing volume flux.
100     ubaro = 0.d0
101     utop = 0.d0
102     cgg)
103    
104     c-- Now, read the control vector.
105     doglobalread = .false.
106     ladinit = .false.
107    
108     if (optimcycle .ge. 0) then
109     ilobcsw=ilnblnk( xx_obcsw_file )
110     write(fnameobcsw(1:80),'(2a,i10.10)')
111     & xx_obcsw_file(1:ilobcsw), '.', optimcycle
112     endif
113    
114     c-- Get the counters, flags, and the interpolation factor.
115     call ctrl_get_gen_rec(
116     I xx_obcswstartdate, xx_obcswperiod,
117     O obcswfac, obcswfirst, obcswchanged,
118     O obcswcount0,obcswcount1,
119     I mytime, myiter, mythid )
120    
121     do iobcs = 1,nobcs
122     if ( obcswfirst ) then
123 heimbach 1.7 call active_read_yz( fnameobcsw, tmpfldyz,
124 heimbach 1.2 & (obcswcount0-1)*nobcs+iobcs,
125     & doglobalread, ladinit, optimcycle,
126     & mythid, xx_obcsw_dummy )
127    
128 mlosch 1.6 #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_Iw(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.5 & 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.5 & - 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_Iw(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.5 & 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.5 & - utop / delR(1)
197 heimbach 1.2 enddo
198     enddo
199     enddo
200     endif
201     endif
202    
203 mlosch 1.6 #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_obcsw1(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 ( (obcswfirst) .or. (obcswchanged)) then
218    
219 edhill 1.4 cgg( This is a terribly long way to do it. However, the dimensions do not exactly
220 heimbach 1.2 cgg match up. I will blame Fortran for the ugliness.
221    
222     do bj = jtlo,jthi
223     do bi = itlo,ithi
224     do k = 1,nr
225     do j = jmin,jmax
226     tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
227     enddo
228     enddo
229     enddo
230     enddo
231    
232     call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
233    
234     do bj = jtlo,jthi
235     do bi = itlo,ithi
236     do k = 1,nr
237     do j = jmin,jmax
238     xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
239     enddo
240     enddo
241     enddo
242     enddo
243    
244 heimbach 1.7 call active_read_yz( fnameobcsw, tmpfldyz,
245 heimbach 1.2 & (obcswcount1-1)*nobcs+iobcs,
246     & doglobalread, ladinit, optimcycle,
247     & mythid, xx_obcsw_dummy )
248    
249 mlosch 1.6 #ifdef ALLOW_CTRL_OBCS_BALANCE
250    
251 heimbach 1.2 if ( optimcycle .gt. 0) then
252     if (iobcs .eq. 3) then
253     cgg Special attention is needed for the normal velocity.
254     cgg For the north, this is the v velocity, iobcs = 4.
255     cgg This is done on a columnwise basis here.
256     do bj = jtlo,jthi
257     do bi = itlo, ithi
258     do j = jmin,jmax
259     i = OB_Iw(J,bi,bj)
260    
261     cgg The barotropic velocity is stored in the level 1.
262     ubaro = tmpfldyz(j,1,bi,bj)
263     tmpfldyz(j,1,bi,bj) = 0.d0
264     utop = 0.d0
265    
266     do k = 1,Nr
267     cgg If cells are not full, this should be modified with hFac.
268     cgg
269     cgg The xx field (tmpfldxz) does not contain the velocity at the
270     cgg surface level. This velocity is not independent; it must
271     cgg exactly balance the volume flux, since we are dealing with
272     cgg the baroclinic velocity structure..
273     utop = tmpfldyz(j,k,bi,bj)*
274 heimbach 1.5 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
275 heimbach 1.2 cgg Add the barotropic velocity component.
276     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
277     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
278     endif
279     enddo
280     cgg Compute the baroclinic velocity at level 1. Should balance flux.
281     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
282 heimbach 1.5 & - utop / delR(1)
283 heimbach 1.2 enddo
284     enddo
285     enddo
286     endif
287     if (iobcs .eq. 4) then
288     cgg Special attention is needed for the normal velocity.
289     cgg For the north, this is the v velocity, iobcs = 4.
290     cgg This is done on a columnwise basis here.
291     do bj = jtlo,jthi
292     do bi = itlo, ithi
293     do j = jmin,jmax
294     i = OB_Iw(J,bi,bj)
295    
296     cgg The barotropic velocity is stored in the level 1.
297     ubaro = tmpfldyz(j,1,bi,bj)
298     tmpfldyz(j,1,bi,bj) = 0.d0
299     utop = 0.d0
300    
301     do k = 1,Nr
302     cgg If cells are not full, this should be modified with hFac.
303     cgg
304     cgg The xx field (tmpfldxz) does not contain the velocity at the
305     cgg surface level. This velocity is not independent; it must
306     cgg exactly balance the volume flux, since we are dealing with
307     cgg the baroclinic velocity structure..
308     utop = tmpfldyz(j,k,bi,bj)*
309 heimbach 1.5 & maskS(i,j,k,bi,bj) * delR(k) + utop
310 heimbach 1.2 cgg Add the barotropic velocity component.
311     if (maskS(i,j,k,bi,bj) .ne. 0.) then
312     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
313     endif
314     enddo
315     cgg Compute the baroclinic velocity at level 1. Should balance flux.
316     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
317 heimbach 1.5 & - utop / delR(1)
318 heimbach 1.2 enddo
319     enddo
320     enddo
321     endif
322     endif
323    
324 mlosch 1.6 #endif /* ALLOW_CTRL_OBCS_BALANCE */
325    
326 heimbach 1.2 do bj = jtlo,jthi
327     do bi = itlo,ithi
328     do k = 1,nr
329     do j = jmin,jmax
330     xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
331     cgg & * maskyz (j,k,bi,bj)
332     enddo
333     enddo
334     enddo
335     enddo
336     endif
337    
338     c-- Add control to model variable.
339     do bj = jtlo, jthi
340     do bi = itlo, ithi
341     c-- Calculate mask for tracer cells (0 => land, 1 => water).
342     do k = 1,nr
343     do j = 1,sny
344     i = OB_Iw(j,bi,bj)
345     if (iobcs .EQ. 1) then
346     OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
347     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
348     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
349     OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
350     & *maskW(i+ip1,j,k,bi,bj)
351     else if (iobcs .EQ. 2) then
352     OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
353     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
354     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
355     OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
356     & *maskW(i+ip1,j,k,bi,bj)
357     else if (iobcs .EQ. 3) then
358     OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
359     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
360     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
361     OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
362     & *maskW(i+ip1,j,k,bi,bj)
363     else if (iobcs .EQ. 4) then
364     OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
365     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
366     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
367     OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
368     & *maskS(i,j,k,bi,bj)
369     endif
370     enddo
371     enddo
372     enddo
373     enddo
374    
375     C-- End over iobcs loop
376     enddo
377    
378     #else /* ALLOW_OBCSW_CONTROL undefined */
379    
380     c == routine arguments ==
381    
382     _RL mytime
383     integer myiter
384     integer mythid
385    
386     c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
387    
388     #endif /* ALLOW_OBCSW_CONTROL */
389    
390     end
391    

  ViewVC Help
Powered by ViewVC 1.1.22