/[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.2 - (show annotations) (download)
Wed Jan 19 12:39:13 2011 UTC (14 years, 11 months ago) by mmazloff
Branch: MAIN
Changes since 1.1: +17 -6 lines
account for partial cells

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 cnma Here we go into the modes' space
123 #ifdef CTRL_OBCS_MODES
124 cih
125 if ( optimcycle .ge. 0) then
126 cih If normal velocity
127 if (iobcs .eq. 3) then
128 cih Begin loop over y-points.
129 do bj = jtlo,jthi
130 do bi = itlo, ithi
131 do j = jmin,jmax
132 cih If open boundary.
133 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
134 i = OB_Iw(j,bi,bj)
135 cih Determine number of open vertical layers.
136 nz = 0
137 do k = 1,Nr
138 nz = nz + maskW(i+ip1,j,k,bi,bj)
139 end do
140 cih Compute absolute velocities from the barotropic-baroclinic modes.
141 do k = 1,Nr
142 if (k.le.nz) then
143 stmp = 0.
144 do nk = 1,nz
145 stmp = stmp +
146 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
147 end do
148 tmpz(k,bi,bj) = stmp
149 else
150 tmpz(k,bi,bj) = 0.
151 end if
152 end do
153 do k = 1,Nr
154 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
155 & *recip_hFacW(i+ip1,j,k,bi,bj)
156 end do
157 cih End open boundary.
158 end if
159 cih End loop over x-points
160 end do
161 end do
162 end do
163 cih End if iobcs = 3.
164 end if
165 cih
166 cih If tangential velocity
167 if (iobcs .eq. 4) then
168 cih Begin loop over y-points.
169 do bj = jtlo,jthi
170 do bi = itlo, ithi
171 do j = jmin,jmax
172 cih If open boundary.
173 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
174 i = OB_Iw(j,bi,bj)
175 cih Determine number of open vertical layers.
176 nz = 0
177 do k = 1,Nr
178 nz = nz + maskS(i,j,k,bi,bj)
179 end do
180 cih Compute absolute velocities from the barotropic-baroclinic modes.
181 do k = 1,Nr
182 if (k.le.nz) then
183 stmp = 0.
184 do nk = 1,nz
185 stmp = stmp +
186 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
187 end do
188 tmpz(k,bi,bj) = stmp
189 else
190 tmpz(k,bi,bj) = 0.
191 end if
192 end do
193 do k = 1,Nr
194 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
195 & *recip_hFacS(i,j,k,bi,bj)
196 end do
197 cih End open boundary.
198 end if
199 cih End loop over x-points
200 end do
201 end do
202 end do
203 cih End if iobcs = 4.
204 end if
205 cih End if optimcycle > 0 .
206 end if
207 cnma
208 #endif
209 cih
210 do bj = jtlo,jthi
211 do bi = itlo,ithi
212 do k = 1,nr
213 do j = jmin,jmax
214 xx_obcsw1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
215 enddo
216 enddo
217 enddo
218 enddo
219 cih End if obcswfirst
220 endif
221 cih
222 if ( (obcswfirst) .or. (obcswchanged) ) then
223 cih
224 do bj = jtlo,jthi
225 do bi = itlo,ithi
226 do k = 1,nr
227 do j = jmin,jmax
228 tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
229 enddo
230 enddo
231 enddo
232 enddo
233 cih
234 call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
235 cih
236 do bj = jtlo,jthi
237 do bi = itlo,ithi
238 do k = 1,nr
239 do j = jmin,jmax
240 xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
241 enddo
242 enddo
243 enddo
244 enddo
245 cih
246 call active_read_yz( fnameobcsw, tmpfldyz,
247 & (obcswcount1-1)*nobcs+iobcs,
248 & doglobalread, ladinit, optimcycle,
249 & mythid, xx_obcsw_dummy )
250 cnma
251 #ifdef CTRL_OBCS_MODES
252 cih
253 if ( optimcycle .ge. 0) then
254 cih If normal velocity
255 if (iobcs .eq. 3) then
256 cih Begin loop over y-points.
257 do bj = jtlo,jthi
258 do bi = itlo, ithi
259 do j = jmin,jmax
260 cih If open boundary.
261 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
262 i = OB_Iw(j,bi,bj)
263 cih Determine number of open vertical layers.
264 nz = 0
265 do k = 1,Nr
266 nz = nz + maskW(i+ip1,j,k,bi,bj)
267 end do
268 cih Compute absolute velocities from the barotropic-baroclinic modes.
269 do k = 1,Nr
270 if (k.le.nz) then
271 stmp = 0.
272 do nk = 1,nz
273 stmp = stmp +
274 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
275 end do
276 tmpz(k,bi,bj) = stmp
277 else
278 tmpz(k,bi,bj) = 0.
279 end if
280 end do
281 do k = 1,Nr
282 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
283 & *recip_hFacW(i+ip1,j,k,bi,bj)
284 end do
285 cih End open boundary.
286 end if
287 cih End loop over x-points
288 end do
289 end do
290 end do
291 cih End if iobcs = 3.
292 end if
293 cih
294 cih If tangential velocity
295 if (iobcs .eq. 4) then
296 cih Begin loop over y-points.
297 do bj = jtlo,jthi
298 do bi = itlo, ithi
299 do j = jmin,jmax
300 cih If open boundary.
301 if ( OB_Iw(j,bi,bj) .ne. 0. ) then
302 i = OB_Iw(j,bi,bj)
303 cih Determine number of open vertical layers.
304 nz = 0
305 do k = 1,Nr
306 nz = nz + maskS(i,j,k,bi,bj)
307 end do
308 cih Compute absolute velocities from the barotropic-baroclinic modes.
309 do k = 1,Nr
310 if (k.le.nz) then
311 stmp = 0.
312 do nk = 1,nz
313 stmp = stmp +
314 & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
315 end do
316 tmpz(k,bi,bj) = stmp
317 else
318 tmpz(k,bi,bj) = 0.
319 end if
320 end do
321 do k = 1,Nr
322 tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
323 & *recip_hFacS(i,j,k,bi,bj)
324 end do
325 cih End open boundary.
326 end if
327 cih End loop over x-points
328 end do
329 end do
330 end do
331 cih End if iobcs = 4.
332 end if
333 cih End if optimcycle > 0 .
334 end if
335 cnma
336 #endif
337 cih
338 do bj = jtlo,jthi
339 do bi = itlo,ithi
340 do k = 1,nr
341 do j = jmin,jmax
342 xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
343 enddo
344 enddo
345 enddo
346 enddo
347 cih End if obcswfirst .or. obcswchanged
348 endif
349
350 c-- Add control to model variable.
351 do bj = jtlo, jthi
352 do bi = itlo, ithi
353 c-- Calculate mask for tracer cells (0 => land, 1 => water).
354 do k = 1,nr
355 do j = 1,sny
356 i = OB_Iw(j,bi,bj)
357 if (iobcs .EQ. 1) then
358 OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
359 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
360 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
361 OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
362 & *maskW(i+ip1,j,k,bi,bj)
363 else if (iobcs .EQ. 2) then
364 OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
365 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
366 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
367 OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
368 & *maskW(i+ip1,j,k,bi,bj)
369 else if (iobcs .EQ. 3) then
370 OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
371 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
372 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
373 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
374 & *maskW(i+ip1,j,k,bi,bj)
375 else if (iobcs .EQ. 4) then
376 OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
377 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
378 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
379 OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
380 & *maskS(i,j,k,bi,bj)
381 endif
382 enddo
383 enddo
384 enddo
385 enddo
386
387 C-- End over iobcs loop
388 enddo
389
390 #else /* ALLOW_OBCSW_CONTROL undefined */
391
392 c == routine arguments ==
393
394 _RL mytime
395 integer myiter
396 integer mythid
397
398 c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
399
400 #endif /* ALLOW_OBCSW_CONTROL */
401
402 end
403

  ViewVC Help
Powered by ViewVC 1.1.22