/[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.4 - (hide 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.3: +8 -8 lines
Replacing delZ by delR in pkg/ctrl/

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     if ( optimcycle .gt. 0) then
130     if (iobcs .eq. 3) then
131     cgg Special attention is needed for the normal velocity.
132     cgg For the north, this is the v velocity, iobcs = 4.
133     cgg This is done on a columnwise basis here.
134     do bj = jtlo,jthi
135     do bi = itlo, ithi
136     do i = imin,imax
137     j = OB_Js(I,bi,bj)
138    
139     cgg The barotropic velocity is stored in the level 1.
140     vbaro = tmpfldxz(i,1,bi,bj)
141     tmpfldxz(i,1,bi,bj) = 0.d0
142     vtop = 0.d0
143    
144     do k = 1,Nr
145     cgg If cells are not full, this should be modified with hFac.
146     cgg
147     cgg The xx field (tmpfldxz) does not contain the velocity at the
148     cgg surface level. This velocity is not independent; it must
149     cgg exactly balance the volume flux, since we are dealing with
150     cgg the baroclinic velocity structure..
151     vtop = tmpfldxz(i,k,bi,bj)*
152 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
153 heimbach 1.2 cgg Add the barotropic velocity component.
154     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
155     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
156     endif
157     enddo
158     cgg Compute the baroclinic velocity at level 1. Should balance flux.
159     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
160 heimbach 1.4 & - vtop / delR(1)
161 heimbach 1.2 enddo
162     enddo
163     enddo
164     endif
165    
166     if (iobcs .eq. 4) then
167     cgg Special attention is needed for the normal velocity.
168     cgg For the north, this is the v velocity, iobcs = 4.
169     cgg This is done on a columnwise basis here.
170     do bj = jtlo,jthi
171     do bi = itlo, ithi
172     do i = imin,imax
173     j = OB_Js(I,bi,bj)
174    
175     cgg The barotropic velocity is stored in the level 1.
176     vbaro = tmpfldxz(i,1,bi,bj)
177     tmpfldxz(i,1,bi,bj) = 0.d0
178     vtop = 0.d0
179    
180     do k = 1,Nr
181     cgg If cells are not full, this should be modified with hFac.
182     cgg
183     cgg The xx field (tmpfldxz) does not contain the velocity at the
184     cgg surface level. This velocity is not independent; it must
185     cgg exactly balance the volume flux, since we are dealing with
186     cgg the baroclinic velocity structure..
187     vtop = tmpfldxz(i,k,bi,bj)*
188 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
189 heimbach 1.2 cgg Add the barotropic velocity component.
190     if (maskW(i,j,k,bi,bj) .ne. 0.) then
191     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
192     endif
193     enddo
194     cgg Compute the baroclinic velocity at level 1. Should balance flux.
195     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
196 heimbach 1.4 & - vtop / delR(1)
197 heimbach 1.2 enddo
198     enddo
199     enddo
200     endif
201     endif
202    
203     do bj = jtlo,jthi
204     do bi = itlo,ithi
205     do k = 1,nr
206     do i = imin,imax
207     xx_obcss1(i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
208     cgg & * maskxz (i,k,bi,bj)
209     enddo
210     enddo
211     enddo
212     enddo
213     endif
214    
215     if ( (obcssfirst) .or. (obcsschanged)) then
216    
217     do bj = jtlo,jthi
218     do bi = itlo,ithi
219     do k = 1,nr
220     do i = imin,imax
221     tmpfldxz(i,k,bi,bj) = xx_obcss1(i,k,bi,bj,iobcs)
222     enddo
223     enddo
224     enddo
225     enddo
226    
227     call exf_swapffields_xz( tmpfldxz2, tmpfldxz, mythid)
228    
229     do bj = jtlo,jthi
230     do bi = itlo,ithi
231     do k = 1,nr
232     do i = imin,imax
233     xx_obcss0(i,k,bi,bj,iobcs) = tmpfldxz2(i,k,bi,bj)
234     enddo
235     enddo
236     enddo
237     enddo
238    
239 heimbach 1.3 call active_read_xz_loc( fnameobcss, tmpfldxz,
240 heimbach 1.2 & (obcsscount1-1)*nobcs+iobcs,
241     & doglobalread, ladinit, optimcycle,
242     & mythid, xx_obcss_dummy )
243    
244     if ( optimcycle .gt. 0) then
245     if (iobcs .eq. 3) then
246     cgg Special attention is needed for the normal velocity.
247     cgg For the north, this is the v velocity, iobcs = 4.
248     cgg This is done on a columnwise basis here.
249     do bj = jtlo,jthi
250     do bi = itlo, ithi
251     do i = imin,imax
252     j = OB_Js(I,bi,bj)
253    
254     cgg The barotropic velocity is stored in the level 1.
255     vbaro = tmpfldxz(i,1,bi,bj)
256     tmpfldxz(i,1,bi,bj) = 0.d0
257     vtop = 0.d0
258    
259     do k = 1,Nr
260     cgg If cells are not full, this should be modified with hFac.
261     cgg
262     cgg The xx field (tmpfldxz) does not contain the velocity at the
263     cgg surface level. This velocity is not independent; it must
264     cgg exactly balance the volume flux, since we are dealing with
265     cgg the baroclinic velocity structure..
266     vtop = tmpfldxz(i,k,bi,bj)*
267 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
268 heimbach 1.2 cgg Add the barotropic velocity component.
269     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
270     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
271     endif
272     enddo
273     cgg Compute the baroclinic velocity at level 1. Should balance flux.
274     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
275 heimbach 1.4 & - vtop / delR(1)
276 heimbach 1.2 enddo
277     enddo
278     enddo
279     endif
280    
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 i = imin,imax
288     j = OB_Js(I,bi,bj)
289    
290     cgg The barotropic velocity is stored in the level 1.
291     vbaro = tmpfldxz(i,1,bi,bj)
292     tmpfldxz(i,1,bi,bj) = 0.d0
293     vtop = 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     vtop = tmpfldxz(i,k,bi,bj)*
303 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
304 heimbach 1.2 cgg Add the barotropic velocity component.
305     if (maskW(i,j,k,bi,bj) .ne. 0.) then
306     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
307     endif
308     enddo
309     cgg Compute the baroclinic velocity at level 1. Should balance flux.
310     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
311 heimbach 1.4 & - vtop / delR(1)
312 heimbach 1.2 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 i = imin,imax
322     xx_obcss1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
323     cgg & * maskxz (i,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 i = 1,snx
336     j = OB_Js(I,bi,bj)
337     if (iobcs .EQ. 1) then
338     OBSt(i,k,bi,bj) = OBSt (i,k,bi,bj)
339     & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
340     & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
341     OBSt(i,k,bi,bj) = OBSt(i,k,bi,bj)
342     & *maskS(i,j+jp1,k,bi,bj)
343     else if (iobcs .EQ. 2) then
344     OBSs(i,k,bi,bj) = OBSs (i,k,bi,bj)
345     & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
346     & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
347     OBSs(i,k,bi,bj) = OBSs(i,k,bi,bj)
348     & *maskS(i,j+jp1,k,bi,bj)
349     else if (iobcs .EQ. 4) then
350     OBSu(i,k,bi,bj) = OBSu (i,k,bi,bj)
351     & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
352     & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
353     OBSu(i,k,bi,bj) = OBSu(i,k,bi,bj)
354     & *maskW(i,j,k,bi,bj)
355     else if (iobcs .EQ. 3) then
356     OBSv(i,k,bi,bj) = OBSv (i,k,bi,bj)
357     & + obcssfac *xx_obcss0(i,k,bi,bj,iobcs)
358     & + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
359     OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
360     & *maskS(i,j+jp1,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_OBCSS_CONTROL undefined */
371    
372     c == routine arguments ==
373    
374     _RL mytime
375     integer myiter
376     integer mythid
377    
378     c-- CPP flag ALLOW_OBCSS_CONTROL undefined.
379    
380     #endif /* ALLOW_OBCSS_CONTROL */
381    
382     end
383    

  ViewVC Help
Powered by ViewVC 1.1.22