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

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

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


Revision 1.5 - (show annotations) (download)
Fri Dec 3 00:48:57 2004 UTC (19 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint57m_post, checkpoint57g_pre, checkpoint57s_post, checkpoint58b_post, checkpoint57b_post, checkpoint57g_post, checkpoint57y_post, checkpoint57h_post, checkpoint57r_post, checkpoint57d_post, checkpoint57i_post, checkpoint59, checkpoint58, checkpoint57, checkpoint58f_post, checkpoint57n_post, checkpoint58d_post, checkpoint58a_post, checkpoint57z_post, checkpoint58y_post, checkpoint58t_post, checkpoint58m_post, checkpoint57l_post, checkpoint57t_post, checkpoint57v_post, checkpoint57f_post, checkpoint57a_post, checkpoint57h_pre, checkpoint58w_post, checkpoint57y_pre, checkpoint58o_post, checkpoint57c_post, checkpoint58p_post, checkpoint58q_post, checkpoint58e_post, checkpoint57c_pre, checkpoint58r_post, checkpoint58n_post, checkpoint57e_post, checkpoint59a, checkpoint57p_post, checkpint57u_post, checkpoint57q_post, eckpoint57e_pre, checkpoint58k_post, checkpoint58v_post, checkpoint58l_post, checkpoint57h_done, checkpoint57j_post, checkpoint57f_pre, checkpoint58g_post, checkpoint58x_post, checkpoint58h_post, checkpoint56c_post, checkpoint58j_post, checkpoint57a_pre, checkpoint57o_post, checkpoint57k_post, checkpoint57w_post, checkpoint58i_post, checkpoint57x_post, checkpoint58c_post, checkpoint58u_post, checkpoint58s_post
Changes since 1.4: +8 -0 lines
o OBCS as control variables
  - update ad_diff.list
  - remove balance of obcs controls from default
  - fix index bug nobcs in ctrl_init
  - fix dummy fields filen in ctrl_pack
  - add dummy weights for obcs

1
2 #include "CTRL_CPPOPTIONS.h"
3 #ifdef ALLOW_OBCS
4 # include "OBCS_OPTIONS.h"
5 #endif
6
7
8 subroutine ctrl_getobcss(
9 I mytime,
10 I myiter,
11 I mythid
12 & )
13
14 c ==================================================================
15 c SUBROUTINE ctrl_getobcss
16 c ==================================================================
17 c
18 c o Get southern 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 new flags: gebbie@mit.edu, 25 Jan 2003.
24 c
25 c ==================================================================
26 c SUBROUTINE ctrl_getobcss
27 c ==================================================================
28
29 implicit none
30
31 #ifdef ALLOW_OBCSS_CONTROL
32
33 c == global variables ==
34
35 #include "EEPARAMS.h"
36 #include "SIZE.h"
37 #include "PARAMS.h"
38 #include "GRID.h"
39 #include "OBCS.h"
40
41 #include "ctrl.h"
42 #include "ctrl_dummy.h"
43 #include "optim.h"
44
45 c == routine arguments ==
46
47 _RL mytime
48 integer myiter
49 integer mythid
50
51 c == local variables ==
52
53 integer bi,bj
54 integer i,j,k
55 integer itlo,ithi
56 integer jtlo,jthi
57 integer jmin,jmax
58 integer imin,imax
59 integer ilobcss
60 integer iobcs
61
62 _RL dummy
63 _RL obcssfac
64 logical obcssfirst
65 logical obcsschanged
66 integer obcsscount0
67 integer obcsscount1
68 integer jp1
69
70 cgg _RL maskxz (1-olx:snx+olx,nr,nsx,nsy)
71
72 logical doglobalread
73 logical ladinit
74
75 character*(80) fnameobcss
76
77 cgg( Variables for splitting barotropic/baroclinic vels.
78 _RL vbaro
79 _RL vtop
80 cgg)
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 jmin = 1-oly
95 jmax = sny+oly
96 imin = 1-olx
97 imax = snx+olx
98 jp1 = 1
99
100 cgg( Initialize variables for balancing volume flux.
101 vbaro = 0.d0
102 vtop = 0.d0
103 cgg)
104
105 c-- Now, read the control vector.
106 doglobalread = .false.
107 ladinit = .false.
108
109 if (optimcycle .ge. 0) then
110 ilobcss=ilnblnk( xx_obcss_file )
111 write(fnameobcss(1:80),'(2a,i10.10)')
112 & xx_obcss_file(1:ilobcss), '.', optimcycle
113 endif
114
115 c-- Get the counters, flags, and the interpolation factor.
116 call ctrl_get_gen_rec(
117 I xx_obcssstartdate, xx_obcssperiod,
118 O obcssfac, obcssfirst, obcsschanged,
119 O obcsscount0,obcsscount1,
120 I mytime, myiter, mythid )
121
122 do iobcs = 1,nobcs
123 if ( obcssfirst ) then
124 call active_read_xz_loc( fnameobcss, tmpfldxz,
125 & (obcsscount0-1)*nobcs+iobcs,
126 & doglobalread, ladinit, optimcycle,
127 & mythid, xx_obcss_dummy )
128
129 #ifdef ALLOW_CTRL_OBCS_BALANCE
130
131 if ( optimcycle .gt. 0) then
132 if (iobcs .eq. 3) then
133 cgg Special attention is needed for the normal velocity.
134 cgg For the north, this is the v velocity, iobcs = 4.
135 cgg This is done on a columnwise basis here.
136 do bj = jtlo,jthi
137 do bi = itlo, ithi
138 do i = imin,imax
139 j = OB_Js(I,bi,bj)
140
141 cgg The barotropic velocity is stored in the level 1.
142 vbaro = tmpfldxz(i,1,bi,bj)
143 tmpfldxz(i,1,bi,bj) = 0.d0
144 vtop = 0.d0
145
146 do k = 1,Nr
147 cgg If cells are not full, this should be modified with hFac.
148 cgg
149 cgg The xx field (tmpfldxz) does not contain the velocity at the
150 cgg surface level. This velocity is not independent; it must
151 cgg exactly balance the volume flux, since we are dealing with
152 cgg the baroclinic velocity structure..
153 vtop = tmpfldxz(i,k,bi,bj)*
154 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
155 cgg Add the barotropic velocity component.
156 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
157 tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
158 endif
159 enddo
160 cgg Compute the baroclinic velocity at level 1. Should balance flux.
161 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
162 & - vtop / delR(1)
163 enddo
164 enddo
165 enddo
166 endif
167
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 i = imin,imax
175 j = OB_Js(I,bi,bj)
176
177 cgg The barotropic velocity is stored in the level 1.
178 vbaro = tmpfldxz(i,1,bi,bj)
179 tmpfldxz(i,1,bi,bj) = 0.d0
180 vtop = 0.d0
181
182 do k = 1,Nr
183 cgg If cells are not full, this should be modified with hFac.
184 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 cgg exactly balance the volume flux, since we are dealing with
188 cgg the baroclinic velocity structure..
189 vtop = tmpfldxz(i,k,bi,bj)*
190 & maskW(i,j,k,bi,bj) * delR(k) + vtop
191 cgg Add the barotropic velocity component.
192 if (maskW(i,j,k,bi,bj) .ne. 0.) then
193 tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
194 endif
195 enddo
196 cgg Compute the baroclinic velocity at level 1. Should balance flux.
197 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
198 & - vtop / delR(1)
199 enddo
200 enddo
201 enddo
202 endif
203 endif
204
205 #endif /* ALLOW_CTRL_OBCS_BALANCE */
206
207 do bj = jtlo,jthi
208 do bi = itlo,ithi
209 do k = 1,nr
210 do i = imin,imax
211 xx_obcss1(i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
212 cgg & * maskxz (i,k,bi,bj)
213 enddo
214 enddo
215 enddo
216 enddo
217 endif
218
219 if ( (obcssfirst) .or. (obcsschanged)) then
220
221 do bj = jtlo,jthi
222 do bi = itlo,ithi
223 do k = 1,nr
224 do i = imin,imax
225 tmpfldxz(i,k,bi,bj) = xx_obcss1(i,k,bi,bj,iobcs)
226 enddo
227 enddo
228 enddo
229 enddo
230
231 call exf_swapffields_xz( tmpfldxz2, tmpfldxz, mythid)
232
233 do bj = jtlo,jthi
234 do bi = itlo,ithi
235 do k = 1,nr
236 do i = imin,imax
237 xx_obcss0(i,k,bi,bj,iobcs) = tmpfldxz2(i,k,bi,bj)
238 enddo
239 enddo
240 enddo
241 enddo
242
243 call active_read_xz_loc( fnameobcss, tmpfldxz,
244 & (obcsscount1-1)*nobcs+iobcs,
245 & doglobalread, ladinit, optimcycle,
246 & mythid, xx_obcss_dummy )
247
248 #ifdef ALLOW_CTRL_OBCS_BALANCE
249
250 if ( optimcycle .gt. 0) then
251 if (iobcs .eq. 3) then
252 cgg Special attention is needed for the normal velocity.
253 cgg For the north, this is the v velocity, iobcs = 4.
254 cgg This is done on a columnwise basis here.
255 do bj = jtlo,jthi
256 do bi = itlo, ithi
257 do i = imin,imax
258 j = OB_Js(I,bi,bj)
259
260 cgg The barotropic velocity is stored in the level 1.
261 vbaro = tmpfldxz(i,1,bi,bj)
262 tmpfldxz(i,1,bi,bj) = 0.d0
263 vtop = 0.d0
264
265 do k = 1,Nr
266 cgg If cells are not full, this should be modified with hFac.
267 cgg
268 cgg The xx field (tmpfldxz) does not contain the velocity at the
269 cgg surface level. This velocity is not independent; it must
270 cgg exactly balance the volume flux, since we are dealing with
271 cgg the baroclinic velocity structure..
272 vtop = tmpfldxz(i,k,bi,bj)*
273 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
274 cgg Add the barotropic velocity component.
275 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
276 tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
277 endif
278 enddo
279 cgg Compute the baroclinic velocity at level 1. Should balance flux.
280 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
281 & - vtop / delR(1)
282 enddo
283 enddo
284 enddo
285 endif
286
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 i = imin,imax
294 j = OB_Js(I,bi,bj)
295
296 cgg The barotropic velocity is stored in the level 1.
297 vbaro = tmpfldxz(i,1,bi,bj)
298 tmpfldxz(i,1,bi,bj) = 0.d0
299 vtop = 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 vtop = tmpfldxz(i,k,bi,bj)*
309 & maskW(i,j,k,bi,bj) * delR(k) + vtop
310 cgg Add the barotropic velocity component.
311 if (maskW(i,j,k,bi,bj) .ne. 0.) then
312 tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
313 endif
314 enddo
315 cgg Compute the baroclinic velocity at level 1. Should balance flux.
316 tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
317 & - vtop / delR(1)
318 enddo
319 enddo
320 enddo
321 endif
322 endif
323
324 #endif /* ALLOW_CTRL_OBCS_BALANCE */
325
326 do bj = jtlo,jthi
327 do bi = itlo,ithi
328 do k = 1,nr
329 do i = imin,imax
330 xx_obcss1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
331 cgg & * maskxz (i,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 i = 1,snx
344 j = OB_Js(I,bi,bj)
345 if (iobcs .EQ. 1) then
346 OBSt(i,k,bi,bj) = OBSt (i,k,bi,bj)
347 & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
348 & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
349 OBSt(i,k,bi,bj) = OBSt(i,k,bi,bj)
350 & *maskS(i,j+jp1,k,bi,bj)
351 else if (iobcs .EQ. 2) then
352 OBSs(i,k,bi,bj) = OBSs (i,k,bi,bj)
353 & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
354 & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
355 OBSs(i,k,bi,bj) = OBSs(i,k,bi,bj)
356 & *maskS(i,j+jp1,k,bi,bj)
357 else if (iobcs .EQ. 4) then
358 OBSu(i,k,bi,bj) = OBSu (i,k,bi,bj)
359 & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
360 & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
361 OBSu(i,k,bi,bj) = OBSu(i,k,bi,bj)
362 & *maskW(i,j,k,bi,bj)
363 else if (iobcs .EQ. 3) then
364 OBSv(i,k,bi,bj) = OBSv (i,k,bi,bj)
365 & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
366 & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
367 OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
368 & *maskS(i,j+jp1,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_OBCSS_CONTROL undefined */
379
380 c == routine arguments ==
381
382 _RL mytime
383 integer myiter
384 integer mythid
385
386 c-- CPP flag ALLOW_OBCSS_CONTROL undefined.
387
388 #endif /* ALLOW_OBCSS_CONTROL */
389
390 end
391

  ViewVC Help
Powered by ViewVC 1.1.22