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

Contents of /MITgcm/pkg/ctrl/ctrl_getobcse.F

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


Revision 1.4 - (show annotations) (download)
Sun Oct 26 00:58:03 2003 UTC (20 years, 6 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
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 call active_read_yz_loc( fnameobcse, tmpfldyz,
124 & (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 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
152 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 & - utop / delR(1)
160 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 & maskS(i,j,k,bi,bj) * delR(k) + utop
187 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 & - utop / delR(1)
195 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 call active_read_yz_loc( fnameobcse, tmpfldyz,
238 & (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 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
266 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 & - utop / delR(1)
274 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 & maskS(i,j,k,bi,bj) * delR(k) + utop
301 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 & - utop / delR(1)
309 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