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

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

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


Revision 1.5 - (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.4: +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_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
69 cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
70
71 logical doglobalread
72 logical ladinit
73
74 character*(80) fnameobcsw
75
76 cgg( Variables for splitting barotropic/baroclinic vels.
77 _RL ubaro
78 _RL utop
79 cgg)
80
81 c == external functions ==
82
83 integer ilnblnk
84 external ilnblnk
85
86
87 c == end of interface ==
88
89 jtlo = mybylo(mythid)
90 jthi = mybyhi(mythid)
91 itlo = mybxlo(mythid)
92 ithi = mybxhi(mythid)
93 jmin = 1-oly
94 jmax = sny+oly
95 imin = 1-olx
96 imax = snx+olx
97 ip1 = 1
98
99 cgg( Initialize variables for balancing volume flux.
100 ubaro = 0.d0
101 utop = 0.d0
102 cgg)
103
104 c-- Now, read the control vector.
105 doglobalread = .false.
106 ladinit = .false.
107
108 if (optimcycle .ge. 0) then
109 ilobcsw=ilnblnk( xx_obcsw_file )
110 write(fnameobcsw(1:80),'(2a,i10.10)')
111 & xx_obcsw_file(1:ilobcsw), '.', optimcycle
112 endif
113
114 c-- Get the counters, flags, and the interpolation factor.
115 call ctrl_get_gen_rec(
116 I xx_obcswstartdate, xx_obcswperiod,
117 O obcswfac, obcswfirst, obcswchanged,
118 O obcswcount0,obcswcount1,
119 I mytime, myiter, mythid )
120
121 do iobcs = 1,nobcs
122 if ( obcswfirst ) then
123 call active_read_yz_loc( fnameobcsw, tmpfldyz,
124 & (obcswcount0-1)*nobcs+iobcs,
125 & doglobalread, ladinit, optimcycle,
126 & mythid, xx_obcsw_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_Iw(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_Iw(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_obcsw1(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 ( (obcswfirst) .or. (obcswchanged)) then
214
215 cgg( This is a terribly long way to do it. However, the dimensions do not exactly
216 cgg match up. I will blame Fortran for the ugliness.
217
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
228 call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
229
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
240 call active_read_yz_loc( fnameobcsw, tmpfldyz,
241 & (obcswcount1-1)*nobcs+iobcs,
242 & doglobalread, ladinit, optimcycle,
243 & mythid, xx_obcsw_dummy )
244
245 if ( optimcycle .gt. 0) then
246 if (iobcs .eq. 3) then
247 cgg Special attention is needed for the normal velocity.
248 cgg For the north, this is the v velocity, iobcs = 4.
249 cgg This is done on a columnwise basis here.
250 do bj = jtlo,jthi
251 do bi = itlo, ithi
252 do j = jmin,jmax
253 i = OB_Iw(J,bi,bj)
254
255 cgg The barotropic velocity is stored in the level 1.
256 ubaro = tmpfldyz(j,1,bi,bj)
257 tmpfldyz(j,1,bi,bj) = 0.d0
258 utop = 0.d0
259
260 do k = 1,Nr
261 cgg If cells are not full, this should be modified with hFac.
262 cgg
263 cgg The xx field (tmpfldxz) does not contain the velocity at the
264 cgg surface level. This velocity is not independent; it must
265 cgg exactly balance the volume flux, since we are dealing with
266 cgg the baroclinic velocity structure..
267 utop = tmpfldyz(j,k,bi,bj)*
268 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
269 cgg Add the barotropic velocity component.
270 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
271 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
272 endif
273 enddo
274 cgg Compute the baroclinic velocity at level 1. Should balance flux.
275 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
276 & - utop / delR(1)
277 enddo
278 enddo
279 enddo
280 endif
281 if (iobcs .eq. 4) then
282 cgg Special attention is needed for the normal velocity.
283 cgg For the north, this is the v velocity, iobcs = 4.
284 cgg This is done on a columnwise basis here.
285 do bj = jtlo,jthi
286 do bi = itlo, ithi
287 do j = jmin,jmax
288 i = OB_Iw(J,bi,bj)
289
290 cgg The barotropic velocity is stored in the level 1.
291 ubaro = tmpfldyz(j,1,bi,bj)
292 tmpfldyz(j,1,bi,bj) = 0.d0
293 utop = 0.d0
294
295 do k = 1,Nr
296 cgg If cells are not full, this should be modified with hFac.
297 cgg
298 cgg The xx field (tmpfldxz) does not contain the velocity at the
299 cgg surface level. This velocity is not independent; it must
300 cgg exactly balance the volume flux, since we are dealing with
301 cgg the baroclinic velocity structure..
302 utop = tmpfldyz(j,k,bi,bj)*
303 & maskS(i,j,k,bi,bj) * delR(k) + utop
304 cgg Add the barotropic velocity component.
305 if (maskS(i,j,k,bi,bj) .ne. 0.) then
306 tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
307 endif
308 enddo
309 cgg Compute the baroclinic velocity at level 1. Should balance flux.
310 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
311 & - utop / delR(1)
312 enddo
313 enddo
314 enddo
315 endif
316 endif
317
318 do bj = jtlo,jthi
319 do bi = itlo,ithi
320 do k = 1,nr
321 do j = jmin,jmax
322 xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
323 cgg & * maskyz (j,k,bi,bj)
324 enddo
325 enddo
326 enddo
327 enddo
328 endif
329
330 c-- Add control to model variable.
331 do bj = jtlo, jthi
332 do bi = itlo, ithi
333 c-- Calculate mask for tracer cells (0 => land, 1 => water).
334 do k = 1,nr
335 do j = 1,sny
336 i = OB_Iw(j,bi,bj)
337 if (iobcs .EQ. 1) then
338 OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
339 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
340 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
341 OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
342 & *maskW(i+ip1,j,k,bi,bj)
343 else if (iobcs .EQ. 2) then
344 OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
345 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
346 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
347 OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
348 & *maskW(i+ip1,j,k,bi,bj)
349 else if (iobcs .EQ. 3) then
350 OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
351 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
352 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
353 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
354 & *maskW(i+ip1,j,k,bi,bj)
355 else if (iobcs .EQ. 4) then
356 OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
357 & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
358 & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
359 OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
360 & *maskS(i,j,k,bi,bj)
361 endif
362 enddo
363 enddo
364 enddo
365 enddo
366
367 C-- End over iobcs loop
368 enddo
369
370 #else /* ALLOW_OBCSW_CONTROL undefined */
371
372 c == routine arguments ==
373
374 _RL mytime
375 integer myiter
376 integer mythid
377
378 c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
379
380 #endif /* ALLOW_OBCSW_CONTROL */
381
382 end
383

  ViewVC Help
Powered by ViewVC 1.1.22