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

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

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


Revision 1.5 - (hide annotations) (download)
Fri Dec 3 00:48:57 2004 UTC (19 years, 6 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 heimbach 1.2
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 heimbach 1.3 call active_read_xz_loc( fnameobcss, tmpfldxz,
125 heimbach 1.2 & (obcsscount0-1)*nobcs+iobcs,
126     & doglobalread, ladinit, optimcycle,
127     & mythid, xx_obcss_dummy )
128    
129 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
130    
131 heimbach 1.2 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 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
155 heimbach 1.2 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 heimbach 1.4 & - vtop / delR(1)
163 heimbach 1.2 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 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
191 heimbach 1.2 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 heimbach 1.4 & - vtop / delR(1)
199 heimbach 1.2 enddo
200     enddo
201     enddo
202     endif
203     endif
204    
205 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
206    
207 heimbach 1.2 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 heimbach 1.3 call active_read_xz_loc( fnameobcss, tmpfldxz,
244 heimbach 1.2 & (obcsscount1-1)*nobcs+iobcs,
245     & doglobalread, ladinit, optimcycle,
246     & mythid, xx_obcss_dummy )
247    
248 mlosch 1.5 #ifdef ALLOW_CTRL_OBCS_BALANCE
249    
250 heimbach 1.2 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 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
274 heimbach 1.2 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 heimbach 1.4 & - vtop / delR(1)
282 heimbach 1.2 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 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
310 heimbach 1.2 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 heimbach 1.4 & - vtop / delR(1)
318 heimbach 1.2 enddo
319     enddo
320     enddo
321     endif
322     endif
323    
324 mlosch 1.5 #endif /* ALLOW_CTRL_OBCS_BALANCE */
325    
326 heimbach 1.2 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