/[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.4 - (hide annotations) (download)
Sun Oct 26 00:58:03 2003 UTC (20 years, 7 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint52l_pre, hrcube4, hrcube5, checkpoint52d_pre, checkpoint56b_post, checkpoint52j_pre, checkpoint54d_post, checkpoint54e_post, checkpoint52l_post, checkpoint52k_post, checkpoint55, checkpoint54, checkpoint56, checkpoint53, checkpoint52, checkpoint52f_post, checkpoint54f_post, checkpoint51t_post, checkpoint55i_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55c_post, checkpoint52e_pre, checkpoint52e_post, checkpoint53d_post, checkpoint52b_pre, checkpoint54b_post, checkpoint52m_post, checkpoint55g_post, checkpoint51q_post, checkpoint52b_post, checkpoint52c_post, checkpoint52f_pre, checkpoint55d_post, checkpoint54a_pre, checkpoint53c_post, checkpoint55d_pre, checkpoint55j_post, checkpoint54a_post, checkpoint55h_post, checkpoint51r_post, checkpoint55b_post, checkpoint53a_post, checkpoint55f_post, checkpoint52d_post, checkpoint53g_post, checkpoint52a_pre, checkpoint52i_post, checkpoint52h_pre, checkpoint56a_post, checkpoint53f_post, checkpoint52j_post, branch-netcdf, checkpoint52n_post, checkpoint53b_pre, checkpoint55a_post, checkpoint51o_post, checkpoint53b_post, checkpoint52a_post, ecco_c52_e35, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-nonh, netcdf-sm0
Changes since 1.3: +8 -8 lines
Replacing delZ by delR in pkg/ctrl/

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.3 call active_read_yz_loc( fnameobcse, tmpfldyz,
124 heimbach 1.2 & (obcsecount0-1)*nobcs+iobcs,
125     & doglobalread, ladinit, optimcycle,
126     & mythid, xx_obcse_dummy )
127    
128     if ( optimcycle .gt. 0) then
129     if (iobcs .eq. 3) then
130     cgg Special attention is needed for the normal velocity.
131     cgg For the north, this is the v velocity, iobcs = 4.
132     cgg This is done on a columnwise basis here.
133     do bj = jtlo,jthi
134     do bi = itlo, ithi
135     do j = jmin,jmax
136     i = OB_Ie(J,bi,bj)
137    
138     cgg The barotropic velocity is stored in the level 1.
139     ubaro = tmpfldyz(j,1,bi,bj)
140     tmpfldyz(j,1,bi,bj) = 0.d0
141     utop = 0.d0
142    
143     do k = 1,Nr
144     cgg If cells are not full, this should be modified with hFac.
145     cgg
146     cgg The xx field (tmpfldxz) does not contain the velocity at the
147     cgg surface level. This velocity is not independent; it must
148     cgg exactly balance the volume flux, since we are dealing with
149     cgg the baroclinic velocity structure..
150     utop = tmpfldyz(j,k,bi,bj)*
151 heimbach 1.4 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
152 heimbach 1.2 cgg Add the barotropic velocity component.
153     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
154     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
155     endif
156     enddo
157     cgg Compute the baroclinic velocity at level 1. Should balance flux.
158     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
159 heimbach 1.4 & - utop / delR(1)
160 heimbach 1.2 enddo
161     enddo
162     enddo
163     endif
164     if (iobcs .eq. 4) then
165     cgg Special attention is needed for the normal velocity.
166     cgg For the north, this is the v velocity, iobcs = 4.
167     cgg This is done on a columnwise basis here.
168     do bj = jtlo,jthi
169     do bi = itlo, ithi
170     do j = jmin,jmax
171     i = OB_Ie(J,bi,bj)
172    
173     cgg The barotropic velocity is stored in the level 1.
174     ubaro = tmpfldyz(j,1,bi,bj)
175     tmpfldyz(j,1,bi,bj) = 0.d0
176     utop = 0.d0
177    
178     do k = 1,Nr
179     cgg If cells are not full, this should be modified with hFac.
180     cgg
181     cgg The xx field (tmpfldxz) does not contain the velocity at the
182     cgg surface level. This velocity is not independent; it must
183     cgg exactly balance the volume flux, since we are dealing with
184     cgg the baroclinic velocity structure..
185     utop = tmpfldyz(j,k,bi,bj)*
186 heimbach 1.4 & maskS(i,j,k,bi,bj) * delR(k) + utop
187 heimbach 1.2 cgg Add the barotropic velocity component.
188     if (maskS(i,j,k,bi,bj) .ne. 0.) then
189     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
190     endif
191     enddo
192     cgg Compute the baroclinic velocity at level 1. Should balance flux.
193     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
194 heimbach 1.4 & - utop / delR(1)
195 heimbach 1.2 enddo
196     enddo
197     enddo
198     endif
199     endif
200    
201     do bj = jtlo,jthi
202     do bi = itlo,ithi
203     do k = 1,nr
204     do j = jmin,jmax
205     xx_obcse1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
206     cgg & * maskyz (j,k,bi,bj)
207     enddo
208     enddo
209     enddo
210     enddo
211     endif
212    
213     if ( (obcsefirst) .or. (obcsechanged)) then
214    
215     do bj = jtlo,jthi
216     do bi = itlo,ithi
217     do j = jmin,jmax
218     do k = 1,nr
219     tmpfldyz(j,k,bi,bj) = xx_obcse1(j,k,bi,bj,iobcs)
220     enddo
221     enddo
222     enddo
223     enddo
224    
225     call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
226    
227     do bj = jtlo,jthi
228     do bi = itlo,ithi
229     do j = jmin,jmax
230     do k = 1,nr
231     xx_obcse0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
232     enddo
233     enddo
234     enddo
235     enddo
236    
237 heimbach 1.3 call active_read_yz_loc( fnameobcse, tmpfldyz,
238 heimbach 1.2 & (obcsecount1-1)*nobcs+iobcs,
239     & doglobalread, ladinit, optimcycle,
240     & mythid, xx_obcse_dummy )
241    
242     if ( optimcycle .gt. 0) then
243     if (iobcs .eq. 3) then
244     cgg Special attention is needed for the normal velocity.
245     cgg For the north, this is the v velocity, iobcs = 4.
246     cgg This is done on a columnwise basis here.
247     do bj = jtlo,jthi
248     do bi = itlo, ithi
249     do j = jmin,jmax
250     i = OB_Ie(J,bi,bj)
251    
252     cgg The barotropic velocity is stored in the level 1.
253     ubaro = tmpfldyz(j,1,bi,bj)
254     tmpfldyz(j,1,bi,bj) = 0.d0
255     utop = 0.d0
256    
257     do k = 1,Nr
258     cgg If cells are not full, this should be modified with hFac.
259     cgg
260     cgg The xx field (tmpfldxz) does not contain the velocity at the
261     cgg surface level. This velocity is not independent; it must
262     cgg exactly balance the volume flux, since we are dealing with
263     cgg the baroclinic velocity structure..
264     utop = tmpfldyz(j,k,bi,bj)*
265 heimbach 1.4 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
266 heimbach 1.2 cgg Add the barotropic velocity component.
267     if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
268     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
269     endif
270     enddo
271     cgg Compute the baroclinic velocity at level 1. Should balance flux.
272     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
273 heimbach 1.4 & - utop / delR(1)
274 heimbach 1.2 enddo
275     enddo
276     enddo
277     endif
278     if (iobcs .eq. 4) then
279     cgg Special attention is needed for the normal velocity.
280     cgg For the north, this is the v velocity, iobcs = 4.
281     cgg This is done on a columnwise basis here.
282     do bj = jtlo,jthi
283     do bi = itlo, ithi
284     do j = jmin,jmax
285     i = OB_Ie(J,bi,bj)
286    
287     cgg The barotropic velocity is stored in the level 1.
288     ubaro = tmpfldyz(j,1,bi,bj)
289     tmpfldyz(j,1,bi,bj) = 0.d0
290     utop = 0.d0
291    
292     do k = 1,Nr
293     cgg If cells are not full, this should be modified with hFac.
294     cgg
295     cgg The xx field (tmpfldxz) does not contain the velocity at the
296     cgg surface level. This velocity is not independent; it must
297     cgg exactly balance the volume flux, since we are dealing with
298     cgg the baroclinic velocity structure..
299     utop = tmpfldyz(j,k,bi,bj)*
300 heimbach 1.4 & maskS(i,j,k,bi,bj) * delR(k) + utop
301 heimbach 1.2 cgg Add the barotropic velocity component.
302     if (maskS(i,j,k,bi,bj) .ne. 0.) then
303     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
304     endif
305     enddo
306     cgg Compute the baroclinic velocity at level 1. Should balance flux.
307     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
308 heimbach 1.4 & - utop / delR(1)
309 heimbach 1.2 enddo
310     enddo
311     enddo
312     endif
313     endif
314    
315     do bj = jtlo,jthi
316     do bi = itlo,ithi
317     do k = 1,nr
318     do j = jmin,jmax
319     xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
320     cgg & * maskyz (j,k,bi,bj)
321     enddo
322     enddo
323     enddo
324     enddo
325     endif
326    
327     c-- Add control to model variable.
328     do bj = jtlo,jthi
329     do bi = itlo,ithi
330     c-- Calculate mask for tracer cells (0 => land, 1 => water).
331     do k = 1,nr
332     do j = 1,sny
333     i = OB_Ie(j,bi,bj)
334     if (iobcs .EQ. 1) then
335     OBEt(j,k,bi,bj) = OBEt (j,k,bi,bj)
336     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
337     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
338     OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
339     & *maskW(i+ip1,j,k,bi,bj)
340     else if (iobcs .EQ. 2) then
341     OBEs(j,k,bi,bj) = OBEs (j,k,bi,bj)
342     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
343     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
344     OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
345     & *maskW(i+ip1,j,k,bi,bj)
346     else if (iobcs .EQ. 3) then
347     OBEu(j,k,bi,bj) = OBEu (j,k,bi,bj)
348     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
349     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
350     OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
351     & *maskW(i+ip1,j,k,bi,bj)
352     else if (iobcs .EQ. 4) then
353     OBEv(j,k,bi,bj) = OBEv (j,k,bi,bj)
354     & + obcsefac *xx_obcse0(j,k,bi,bj,iobcs)
355     & + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
356     OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
357     & *maskS(i,j,k,bi,bj)
358     endif
359     enddo
360     enddo
361     enddo
362     enddo
363    
364     C-- End over iobcs loop
365     enddo
366    
367     #else /* ALLOW_OBCSE_CONTROL undefined */
368    
369     c == routine arguments ==
370    
371     _RL mytime
372     integer myiter
373     integer mythid
374    
375     c-- CPP flag ALLOW_OBCSE_CONTROL undefined.
376    
377     #endif /* ALLOW_OBCSE_CONTROL */
378    
379     end
380    

  ViewVC Help
Powered by ViewVC 1.1.22