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

Annotation of /MITgcm/pkg/ctrl/ctrl_getobcsn.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: +38 -36 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_getobcsn(
11     I mytime,
12     I myiter,
13     I mythid
14     & )
15    
16     c ==================================================================
17     c SUBROUTINE ctrl_getobcsn
18     c ==================================================================
19     c
20     c o Get northern 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_getobcsn
28     c ==================================================================
29    
30     implicit none
31    
32     #ifdef ALLOW_OBCSN_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 ilobcsn
61     integer iobcs
62     integer jp1
63    
64     _RL dummy
65     _RL obcsnfac
66     logical obcsnfirst
67     logical obcsnchanged
68     integer obcsncount0
69     integer obcsncount1
70    
71     cgg _RL maskxz (1-olx:snx+olx,nr,nsx,nsy)
72    
73     logical doglobalread
74     logical ladinit
75    
76     character*(80) fnameobcsn
77    
78     cgg( Variables for splitting barotropic/baroclinic vels.
79     _RL vbaro
80     _RL vtop
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     cgg jmin = 1-oly
95     cgg jmax = sny+oly
96     cgg imin = 1-olx
97     cgg imax = snx+olx
98    
99     jmin = 1
100     jmax = sny
101     imin = 1
102     imax = snx
103     jp1 = 0
104    
105     cgg( Initialize variables for balancing volume flux.
106     vbaro = 0.d0
107     vtop = 0.d0
108     cgg)
109    
110     c-- Now, read the control vector.
111     doglobalread = .false.
112     ladinit = .false.
113    
114     if (optimcycle .ge. 0) then
115     ilobcsn=ilnblnk( xx_obcsn_file )
116 jmc 1.7 write(fnameobcsn(1:80),'(2a,i10.10)')
117 heimbach 1.2 & xx_obcsn_file(1:ilobcsn), '.', optimcycle
118     endif
119    
120     c-- Get the counters, flags, and the interpolation factor.
121     call ctrl_get_gen_rec(
122     I xx_obcsnstartdate, xx_obcsnperiod,
123     O obcsnfac, obcsnfirst, obcsnchanged,
124     O obcsncount0,obcsncount1,
125     I mytime, myiter, mythid )
126    
127     do iobcs = 1,nobcs
128     if ( obcsnfirst ) then
129 jmc 1.7 call active_read_xz( fnameobcsn, tmpfldxz,
130 heimbach 1.2 & (obcsncount0-1)*nobcs+iobcs,
131     & doglobalread, ladinit, optimcycle,
132     & mythid, xx_obcsn_dummy )
133    
134 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
135    
136 jmc 1.7 if ( optimcycle .gt. 0) then
137 heimbach 1.2 if (iobcs .eq. 3) then
138     cgg Special attention is needed for the normal velocity.
139     cgg For the north, this is the v velocity, iobcs = 4.
140     cgg This is done on a columnwise basis here.
141     do bj = jtlo,jthi
142     do bi = itlo, ithi
143     do i = imin,imax
144    
145     cgg The barotropic velocity is stored in the level 1.
146     vbaro = tmpfldxz(i,1,bi,bj)
147     cgg Except for the special point which balances barotropic vol.flux.
148 jmc 1.7 cgg Special column in the NW corner.
149 heimbach 1.2 j = OB_Jn(I,bi,bj)
150     if (ob_iw(j,bi,bj).eq.(i-1).and.
151     & ob_iw(j,bi,bj).ne. 0) then
152     print*,'Apply shiftvel1 @ i,j'
153     print*,shiftvel(1),i,j
154     vbaro = shiftvel(1)
155     endif
156     tmpfldxz(i,1,bi,bj) = 0.d0
157     vtop = 0.d0
158    
159     do k = 1,Nr
160     cgg If cells are not full, this should be modified with hFac.
161 jmc 1.7 cgg
162     cgg The xx field (tmpfldxz) does not contain the velocity at the
163     cgg surface level. This velocity is not independent; it must
164 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
165 jmc 1.7 cgg the baroclinic velocity structure..
166 heimbach 1.2 vtop = tmpfldxz(i,k,bi,bj)*
167 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
168 jmc 1.7 cgg Add the barotropic velocity component.
169 heimbach 1.2 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
170     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
171     endif
172     enddo
173     cgg Compute the baroclinic velocity at level 1. Should balance flux.
174 jmc 1.7 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
175 heimbach 1.4 & - vtop / delR(1)
176 heimbach 1.2 enddo
177     enddo
178     enddo
179     endif
180    
181     if (iobcs .eq. 4) then
182     cgg Special attention is needed for the normal velocity.
183     cgg For the north, this is the v velocity, iobcs = 4.
184     cgg This is done on a columnwise basis here.
185     do bj = jtlo,jthi
186     do bi = itlo, ithi
187     do i = imin,imax
188    
189     cgg The barotropic velocity is stored in the level 1.
190     vbaro = tmpfldxz(i,1,bi,bj)
191     cgg Except for the special point which balances barotropic vol.flux.
192 jmc 1.7 cgg Special column in the NW corner.
193 heimbach 1.2 j = OB_Jn(I,bi,bj)
194     tmpfldxz(i,1,bi,bj) = 0.d0
195     vtop = 0.d0
196    
197     do k = 1,Nr
198     cgg If cells are not full, this should be modified with hFac.
199 jmc 1.7 cgg
200     cgg The xx field (tmpfldxz) does not contain the velocity at the
201     cgg surface level. This velocity is not independent; it must
202 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
203 jmc 1.7 cgg the baroclinic velocity structure..
204 heimbach 1.2 vtop = tmpfldxz(i,k,bi,bj)*
205 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
206 jmc 1.7 cgg Add the barotropic velocity component.
207 heimbach 1.2 if (maskW(i,j,k,bi,bj) .ne. 0.) then
208     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
209     endif
210     enddo
211     cgg Compute the baroclinic velocity at level 1. Should balance flux.
212 jmc 1.7 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
213 heimbach 1.4 & - vtop / delR(1)
214 heimbach 1.2 enddo
215     enddo
216     enddo
217     endif
218    
219     endif
220    
221 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
222    
223 heimbach 1.2 do bj = jtlo,jthi
224     do bi = itlo,ithi
225     do k = 1,nr
226 jmc 1.7 do i = imin,imax
227 heimbach 1.2 xx_obcsn1(i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
228     cgg & * maskxz (i,k,bi,bj)
229     enddo
230     enddo
231     enddo
232     enddo
233     endif
234    
235     if ( (obcsnfirst) .or. (obcsnchanged)) then
236 jmc 1.7
237 heimbach 1.2 do bj = jtlo,jthi
238     do bi = itlo,ithi
239     do k = 1,nr
240     do i = imin,imax
241     tmpfldxz(i,k,bi,bj) = xx_obcsn1(i,k,bi,bj,iobcs)
242     enddo
243     enddo
244     enddo
245     enddo
246    
247     call exf_swapffields_xz( tmpfldxz2, tmpfldxz, mythid)
248    
249     do bj = jtlo,jthi
250     do bi = itlo,ithi
251     do k = 1,nr
252     do i = imin,imax
253     xx_obcsn0(i,k,bi,bj,iobcs) = tmpfldxz2(i,k,bi,bj)
254     enddo
255     enddo
256     enddo
257     enddo
258    
259 jmc 1.7 call active_read_xz( fnameobcsn, tmpfldxz,
260 heimbach 1.2 & (obcsncount1-1)*nobcs+iobcs,
261     & doglobalread, ladinit, optimcycle,
262     & mythid, xx_obcsn_dummy )
263    
264 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
265    
266 jmc 1.7 if ( optimcycle .gt. 0) then
267 heimbach 1.2 if (iobcs .eq. 3) then
268     cgg Special attention is needed for the normal velocity.
269     cgg For the north, this is the v velocity, iobcs = 3.
270     cgg This is done on a columnwise basis here.
271     do bj = jtlo,jthi
272     do bi = itlo, ithi
273     do i = imin,imax
274    
275     cgg The barotropic velocity is stored in the level 1.
276     vbaro = tmpfldxz(i,1,bi,bj)
277     cgg Except for the special point which balances barotropic vol.flux.
278 jmc 1.7 cgg Special column in the NW corner.
279 heimbach 1.2 j = OB_Jn(I,bi,bj)
280     if (ob_iw(j,bi,bj).eq.(i-1).and.
281     & ob_iw(j,bi,bj).ne. 0) then
282     print*,'correct vbaro',vbaro
283     print*,'Apply shiftvel2 @ i,j'
284     print*,shiftvel(2),i,j
285     vbaro = shiftvel(2)
286     endif
287     tmpfldxz(i,1,bi,bj) = 0.d0
288     vtop = 0.d0
289    
290     do k = 1,Nr
291     cgg If cells are not full, this should be modified with hFac.
292 jmc 1.7 cgg
293     cgg The xx field (tmpfldxz) does not contain the velocity at the
294     cgg surface level. This velocity is not independent; it must
295 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
296 jmc 1.7 cgg the baroclinic velocity structure..
297 heimbach 1.2 vtop = tmpfldxz(i,k,bi,bj)*
298 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
299 jmc 1.7 cgg Add the barotropic velocity component.
300 heimbach 1.2 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
301     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
302     endif
303     enddo
304     cgg Compute the baroclinic velocity at level 1. Should balance flux.
305 jmc 1.7 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
306 heimbach 1.4 & - vtop / delR(1)
307 heimbach 1.2 enddo
308     enddo
309     enddo
310     endif
311     if (iobcs .eq. 4) then
312     cgg Special attention is needed for the normal velocity.
313     cgg For the north, this is the v velocity, iobcs = 3.
314     cgg This is done on a columnwise basis here.
315     do bj = jtlo,jthi
316     do bi = itlo, ithi
317     do i = imin,imax
318    
319     cgg The barotropic velocity is stored in the level 1.
320     vbaro = tmpfldxz(i,1,bi,bj)
321     cgg Except for the special point which balances barotropic vol.flux.
322 jmc 1.7 cgg Special column in the NW corner.
323 heimbach 1.2 j = OB_Jn(I,bi,bj)
324     tmpfldxz(i,1,bi,bj) = 0.d0
325     vtop = 0.d0
326    
327     do k = 1,Nr
328     cgg If cells are not full, this should be modified with hFac.
329 jmc 1.7 cgg
330     cgg The xx field (tmpfldxz) does not contain the velocity at the
331     cgg surface level. This velocity is not independent; it must
332 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
333 jmc 1.7 cgg the baroclinic velocity structure..
334 heimbach 1.2 vtop = tmpfldxz(i,k,bi,bj)*
335 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
336 jmc 1.7 cgg Add the barotropic velocity component.
337 heimbach 1.2 if (maskW(i,j,k,bi,bj) .ne. 0.) then
338     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
339     endif
340     enddo
341     cgg Compute the baroclinic velocity at level 1. Should balance flux.
342 jmc 1.7 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
343 heimbach 1.4 & - vtop / delR(1)
344 heimbach 1.2 enddo
345     enddo
346     enddo
347     endif
348     endif
349    
350 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
351    
352 heimbach 1.2 do bj = jtlo,jthi
353     do bi = itlo,ithi
354     do k = 1,nr
355     do i = imin,imax
356     xx_obcsn1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
357     cgg & * maskxz (i,k,bi,bj)
358     enddo
359     enddo
360     enddo
361     enddo
362    
363     endif
364    
365     c-- Add control to model variable.
366     do bj = jtlo,jthi
367     do bi = itlo,ithi
368     c-- Calculate mask for tracer cells (0 => land, 1 => water).
369     do k = 1,nr
370     do i = 1,snx
371     j = OB_Jn(I,bi,bj)
372     if (iobcs .EQ. 1) then
373     OBNt(i,k,bi,bj) = OBNt (i,k,bi,bj)
374     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
375     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
376     OBNt(i,k,bi,bj) = OBNt(i,k,bi,bj)
377     & *maskS(i,j+jp1,k,bi,bj)
378     cgg & *maskxz(i,k,bi,bj)
379     else if (iobcs .EQ. 2) then
380     OBNs(i,k,bi,bj) = OBNs (i,k,bi,bj)
381     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
382     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
383     OBNs(i,k,bi,bj) = OBNs(i,k,bi,bj)
384     & *maskS(i,j+jp1,k,bi,bj)
385     cgg & *maskxz(i,k,bi,bj)
386     else if (iobcs .EQ. 4) then
387     OBNu(i,k,bi,bj) = OBNu (i,k,bi,bj)
388     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
389     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
390     OBNu(i,k,bi,bj) = OBNu(i,k,bi,bj)
391     & *maskW(i,j,k,bi,bj)
392     cgg & *maskxz(i,k,bi,bj)
393     else if (iobcs .EQ. 3) then
394     OBNv(i,k,bi,bj) = OBNv (i,k,bi,bj)
395     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
396     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
397     OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
398     & *maskS(i,j+jp1,k,bi,bj)
399     cgg & *maskxz(i,k,bi,bj)
400     endif
401     enddo
402     enddo
403     enddo
404     enddo
405    
406     C-- End over iobcs loop
407     enddo
408 jmc 1.7
409 heimbach 1.2 #else /* ALLOW_OBCSN_CONTROL undefined */
410    
411     c == routine arguments ==
412    
413     _RL mytime
414     integer myiter
415     integer mythid
416    
417     c-- CPP flag ALLOW_OBCSN_CONTROL undefined.
418    
419     #endif /* ALLOW_OBCSN_CONTROL */
420    
421     end
422    
423    
424    
425    
426    

  ViewVC Help
Powered by ViewVC 1.1.22