/[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.8 - (hide annotations) (download)
Tue Oct 9 00:00:00 2007 UTC (16 years, 11 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.7: +32 -30 lines
add missing cvs $Header:$ or $Name:$

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

  ViewVC Help
Powered by ViewVC 1.1.22