/[MITgcm]/MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcsw.F
ViewVC logotype

Contents of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcsw.F

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


Revision 1.1 - (show annotations) (download)
Tue Jan 18 19:33:08 2011 UTC (14 years, 11 months ago) by mmazloff
Branch: MAIN
Adding in simple box model with code for controlling obcs with modal decomposition.

1
2 #include "CTRL_CPPOPTIONS.h"
3 #ifdef ALLOW_OBCS
4 # include "OBCS_OPTIONS.h"
5 #endif
6
7
8 subroutine ctrl_getobcsw(
9 I mytime,
10 I myiter,
11 I mythid
12 & )
13
14 c ==================================================================
15 c SUBROUTINE ctrl_getobcsw
16 c ==================================================================
17 c
18 c o Get western 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 modified: gebbie@mit.edu, 18-Mar-2003
24 c ==================================================================
25 c SUBROUTINE ctrl_getobcsw
26 c ==================================================================
27
28 implicit none
29
30 #ifdef ALLOW_OBCSW_CONTROL
31
32 c == global variables ==
33
34 #include "EEPARAMS.h"
35 #include "SIZE.h"
36 #include "PARAMS.h"
37 #include "GRID.h"
38 #include "OBCS.h"
39
40 #include "ctrl.h"
41 #include "ctrl_dummy.h"
42 #include "optim.h"
43
44 c == routine arguments ==
45
46 _RL mytime
47 integer myiter
48 integer mythid
49
50 c == local variables ==
51
52 integer bi,bj
53 integer i,j,k
54 integer itlo,ithi
55 integer jtlo,jthi
56 integer jmin,jmax
57 integer imin,imax
58 integer ilobcsw
59 integer iobcs
60 integer ip1
61
62 _RL dummy
63 _RL obcswfac
64 logical obcswfirst
65 logical obcswchanged
66 integer obcswcount0
67 integer obcswcount1
68 cih
69 Integer nk,nz
70 _RL tmpz (nr,nsx,nsy)
71 _RL stmp
72 character*(80) fnamein
73 cih
74 logical doglobalread
75 logical ladinit
76
77 character*(80) fnameobcsw
78
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 = 1
97
98
99 c-- Now, read the control vector.
100 doglobalread = .false.
101 ladinit = .false.
102
103 if (optimcycle .ge. 0) then
104 ilobcsw=ilnblnk( xx_obcsw_file )
105 write(fnameobcsw(1:80),'(2a,i10.10)')
106 & xx_obcsw_file(1:ilobcsw), '.', optimcycle
107 endif
108
109 c-- Get the counters, flags, and the interpolation factor.
110 call ctrl_get_gen_rec(
111 I xx_obcswstartdate, xx_obcswperiod,
112 O obcswfac, obcswfirst, obcswchanged,
113 O obcswcount0,obcswcount1,
114 I mytime, myiter, mythid )
115 cih
116 do iobcs = 1,nobcs
117 if ( obcswfirst ) then
118 call active_read_yz( fnameobcsw, tmpfldyz,
119 & (obcswcount0-1)*nobcs+iobcs,
120 & doglobalread, ladinit, optimcycle,
121 & mythid, xx_obcsw_dummy )
122 cih
123 if ( optimcycle .ge. 0) then
124 cih If normal velocity
125 if (iobcs .eq. 3) then
126 cih Begin loop over y-points.
127 do bj = jtlo,jthi
128 do bi = itlo, ithi
129 do j = jmin,jmax
130 cih If open boundary.
131 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
132 i = OB_Iw(j,bi,bj)
133 cih Determine number of open vertical layers.
134 nz = 0
135 do k = 1,Nr
136 nz = nz + maskW(i+ip1,j,k,bi,bj)
137 end do
138 cih Compute absolute velocities from the barotropic-baroclinic modes.
139 do k = 1,Nr
140 if (k.le.nz) then
141 stmp = 0.
142 do nk = 1,nz
143 stmp = stmp +
144 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
145 end do
146 tmpz(k,bi,bj) = stmp
147 else
148 tmpz(k,bi,bj) = 0.
149 end if
150 end do
151 do k = 1,Nr
152 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
153 end do
154 cih End open boundary.
155 end if
156 cih End loop over x-points
157 end do
158 end do
159 end do
160 cih End if iobcs = 3.
161 end if
162 cih
163 cih If tangential velocity
164 if (iobcs .eq. 4) then
165 cih Begin loop over y-points.
166 do bj = jtlo,jthi
167 do bi = itlo, ithi
168 do j = jmin,jmax
169 cih If open boundary.
170 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
171 i = OB_Iw(j,bi,bj)
172 cih Determine number of open vertical layers.
173 nz = 0
174 do k = 1,Nr
175 nz = nz + maskS(i,j,k,bi,bj)
176 end do
177 cih Compute absolute velocities from the barotropic-baroclinic modes.
178 do k = 1,Nr
179 if (k.le.nz) then
180 stmp = 0.
181 do nk = 1,nz
182 stmp = stmp +
183 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
184 end do
185 tmpz(k,bi,bj) = stmp
186 else
187 tmpz(k,bi,bj) = 0.
188 end if
189 end do
190 do k = 1,Nr
191 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
192 end do
193 cih End open boundary.
194 end if
195 cih End loop over x-points
196 end do
197 end do
198 end do
199 cih End if iobcs = 4.
200 end if
201 cih End if optimcycle > 0 .
202 end if
203 cih
204 do bj = jtlo,jthi
205 do bi = itlo,ithi
206 do k = 1,nr
207 do j = jmin,jmax
208 xx_obcsw1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
209 enddo
210 enddo
211 enddo
212 enddo
213 cih End if obcswfirst
214 endif
215 cih
216 if ( (obcswfirst) .or. (obcswchanged) ) then
217 cih
218 do bj = jtlo,jthi
219 do bi = itlo,ithi
220 do k = 1,nr
221 do j = jmin,jmax
222 tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
223 enddo
224 enddo
225 enddo
226 enddo
227 cih
228 call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
229 cih
230 do bj = jtlo,jthi
231 do bi = itlo,ithi
232 do k = 1,nr
233 do j = jmin,jmax
234 xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
235 enddo
236 enddo
237 enddo
238 enddo
239 cih
240 call active_read_yz( fnameobcsw, tmpfldyz,
241 & (obcswcount1-1)*nobcs+iobcs,
242 & doglobalread, ladinit, optimcycle,
243 & mythid, xx_obcsw_dummy )
244 cih
245 if ( optimcycle .ge. 0) then
246
247 cih If normal velocity
248 if (iobcs .eq. 3) then
249 cih Begin loop over y-points.
250 do bj = jtlo,jthi
251 do bi = itlo, ithi
252 do j = jmin,jmax
253 cih If open boundary.
254 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
255 i = OB_Iw(j,bi,bj)
256 cih Determine number of open vertical layers.
257 nz = 0
258 do k = 1,Nr
259 nz = nz + maskW(i+ip1,j,k,bi,bj)
260 end do
261 cih Compute absolute velocities from the barotropic-baroclinic modes.
262 do k = 1,Nr
263 if (k.le.nz) then
264 stmp = 0.
265 do nk = 1,nz
266 stmp = stmp +
267 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
268 end do
269 tmpz(k,bi,bj) = stmp
270 else
271 tmpz(k,bi,bj) = 0.
272 end if
273 end do
274 do k = 1,Nr
275 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
276 end do
277 cih End open boundary.
278 end if
279 cih End loop over x-points
280 end do
281 end do
282 end do
283 cih End if iobcs = 3.
284 end if
285 cih
286 cih If tangential velocity
287 if (iobcs .eq. 4) then
288 cih Begin loop over y-points.
289 do bj = jtlo,jthi
290 do bi = itlo, ithi
291 do j = jmin,jmax
292 cih If open boundary.
293 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
294 i = OB_Iw(j,bi,bj)
295 cih Determine number of open vertical layers.
296 nz = 0
297 do k = 1,Nr
298 nz = nz + maskS(i,j,k,bi,bj)
299 end do
300 cih Compute absolute velocities from the barotropic-baroclinic modes.
301 do k = 1,Nr
302 if (k.le.nz) then
303 stmp = 0.
304 do nk = 1,nz
305 stmp = stmp +
306 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
307 end do
308 tmpz(k,bi,bj) = stmp
309 else
310 tmpz(k,bi,bj) = 0.
311 end if
312 end do
313 do k = 1,Nr
314 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
315 end do
316 cih End open boundary.
317 end if
318 cih End loop over x-points
319 end do
320 end do
321 end do
322 cih End if iobcs = 4.
323 end if
324 cih End if optimcycle > 0 .
325 end if
326 cih
327 do bj = jtlo,jthi
328 do bi = itlo,ithi
329 do k = 1,nr
330 do j = jmin,jmax
331 xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
332 enddo
333 enddo
334 enddo
335 enddo
336 cih End if obcswfirst .or. obcswchanged
337 endif
338
339 c-- Add control to model variable.
340 do bj = jtlo, jthi
341 do bi = itlo, ithi
342 c-- Calculate mask for tracer cells (0 => land, 1 => water).
343 do k = 1,nr
344 do j = 1,sny
345 i = OB_Iw(j,bi,bj)
346 if (iobcs .EQ. 1) then
347 OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
348 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
349 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
350 OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
351 & *maskW(i+ip1,j,k,bi,bj)
352 else if (iobcs .EQ. 2) then
353 OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
354 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
355 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
356 OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
357 & *maskW(i+ip1,j,k,bi,bj)
358 else if (iobcs .EQ. 3) then
359 OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
360 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
361 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
362 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
363 & *maskW(i+ip1,j,k,bi,bj)
364 else if (iobcs .EQ. 4) then
365 OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
366 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
367 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
368 OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
369 & *maskS(i,j,k,bi,bj)
370 endif
371 enddo
372 enddo
373 enddo
374 enddo
375
376 C-- End over iobcs loop
377 enddo
378
379 #else /* ALLOW_OBCSW_CONTROL undefined */
380
381 c == routine arguments ==
382
383 _RL mytime
384 integer myiter
385 integer mythid
386
387 c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
388
389 #endif /* ALLOW_OBCSW_CONTROL */
390
391 end
392

  ViewVC Help
Powered by ViewVC 1.1.22