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

Annotation of /MITgcm/pkg/ctrl/ctrl_getobcsn.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, 7 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_getobcsn(
9     I mytime,
10     I myiter,
11     I mythid
12     & )
13    
14     c ==================================================================
15     c SUBROUTINE ctrl_getobcsn
16     c ==================================================================
17     c
18     c o Get northern 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_getobcsn
26     c ==================================================================
27    
28     implicit none
29    
30     #ifdef ALLOW_OBCSN_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 ilobcsn
59     integer iobcs
60     integer jp1
61    
62     _RL dummy
63     _RL obcsnfac
64     logical obcsnfirst
65     logical obcsnchanged
66     integer obcsncount0
67     integer obcsncount1
68    
69     cgg _RL maskxz (1-olx:snx+olx,nr,nsx,nsy)
70    
71     logical doglobalread
72     logical ladinit
73    
74     character*(80) fnameobcsn
75    
76     cgg( Variables for splitting barotropic/baroclinic vels.
77     _RL vbaro
78     _RL vtop
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     cgg jmin = 1-oly
93     cgg jmax = sny+oly
94     cgg imin = 1-olx
95     cgg imax = snx+olx
96    
97     jmin = 1
98     jmax = sny
99     imin = 1
100     imax = snx
101     jp1 = 0
102    
103     cgg( Initialize variables for balancing volume flux.
104     vbaro = 0.d0
105     vtop = 0.d0
106     cgg)
107    
108     c-- Now, read the control vector.
109     doglobalread = .false.
110     ladinit = .false.
111    
112     if (optimcycle .ge. 0) then
113     ilobcsn=ilnblnk( xx_obcsn_file )
114     write(fnameobcsn(1:80),'(2a,i10.10)')
115     & xx_obcsn_file(1:ilobcsn), '.', optimcycle
116     endif
117    
118     c-- Get the counters, flags, and the interpolation factor.
119     call ctrl_get_gen_rec(
120     I xx_obcsnstartdate, xx_obcsnperiod,
121     O obcsnfac, obcsnfirst, obcsnchanged,
122     O obcsncount0,obcsncount1,
123     I mytime, myiter, mythid )
124    
125     do iobcs = 1,nobcs
126     if ( obcsnfirst ) then
127 heimbach 1.3 call active_read_xz_loc( fnameobcsn, tmpfldxz,
128 heimbach 1.2 & (obcsncount0-1)*nobcs+iobcs,
129     & doglobalread, ladinit, optimcycle,
130     & mythid, xx_obcsn_dummy )
131    
132     if ( optimcycle .gt. 0) then
133     if (iobcs .eq. 3) then
134     cgg Special attention is needed for the normal velocity.
135     cgg For the north, this is the v velocity, iobcs = 4.
136     cgg This is done on a columnwise basis here.
137     do bj = jtlo,jthi
138     do bi = itlo, ithi
139     do i = imin,imax
140    
141     cgg The barotropic velocity is stored in the level 1.
142     vbaro = tmpfldxz(i,1,bi,bj)
143     cgg Except for the special point which balances barotropic vol.flux.
144     cgg Special column in the NW corner.
145     j = OB_Jn(I,bi,bj)
146     if (ob_iw(j,bi,bj).eq.(i-1).and.
147     & ob_iw(j,bi,bj).ne. 0) then
148     print*,'Apply shiftvel1 @ i,j'
149     print*,shiftvel(1),i,j
150     vbaro = shiftvel(1)
151     endif
152     tmpfldxz(i,1,bi,bj) = 0.d0
153     vtop = 0.d0
154    
155     do k = 1,Nr
156     cgg If cells are not full, this should be modified with hFac.
157     cgg
158     cgg The xx field (tmpfldxz) does not contain the velocity at the
159     cgg surface level. This velocity is not independent; it must
160     cgg exactly balance the volume flux, since we are dealing with
161     cgg the baroclinic velocity structure..
162     vtop = tmpfldxz(i,k,bi,bj)*
163 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
164 heimbach 1.2 cgg Add the barotropic velocity component.
165     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
166     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
167     endif
168     enddo
169     cgg Compute the baroclinic velocity at level 1. Should balance flux.
170     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
171 heimbach 1.4 & - vtop / delR(1)
172 heimbach 1.2 enddo
173     enddo
174     enddo
175     endif
176    
177     if (iobcs .eq. 4) then
178     cgg Special attention is needed for the normal velocity.
179     cgg For the north, this is the v velocity, iobcs = 4.
180     cgg This is done on a columnwise basis here.
181     do bj = jtlo,jthi
182     do bi = itlo, ithi
183     do i = imin,imax
184    
185     cgg The barotropic velocity is stored in the level 1.
186     vbaro = tmpfldxz(i,1,bi,bj)
187     cgg Except for the special point which balances barotropic vol.flux.
188     cgg Special column in the NW corner.
189     j = OB_Jn(I,bi,bj)
190     tmpfldxz(i,1,bi,bj) = 0.d0
191     vtop = 0.d0
192    
193     do k = 1,Nr
194     cgg If cells are not full, this should be modified with hFac.
195     cgg
196     cgg The xx field (tmpfldxz) does not contain the velocity at the
197     cgg surface level. This velocity is not independent; it must
198     cgg exactly balance the volume flux, since we are dealing with
199     cgg the baroclinic velocity structure..
200     vtop = tmpfldxz(i,k,bi,bj)*
201 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
202 heimbach 1.2 cgg Add the barotropic velocity component.
203     if (maskW(i,j,k,bi,bj) .ne. 0.) then
204     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
205     endif
206     enddo
207     cgg Compute the baroclinic velocity at level 1. Should balance flux.
208     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
209 heimbach 1.4 & - vtop / delR(1)
210 heimbach 1.2 enddo
211     enddo
212     enddo
213     endif
214    
215     endif
216    
217     do bj = jtlo,jthi
218     do bi = itlo,ithi
219     do k = 1,nr
220     do i = imin,imax
221     xx_obcsn1(i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
222     cgg & * maskxz (i,k,bi,bj)
223     enddo
224     enddo
225     enddo
226     enddo
227     endif
228    
229     if ( (obcsnfirst) .or. (obcsnchanged)) then
230    
231     do bj = jtlo,jthi
232     do bi = itlo,ithi
233     do k = 1,nr
234     do i = imin,imax
235     tmpfldxz(i,k,bi,bj) = xx_obcsn1(i,k,bi,bj,iobcs)
236     enddo
237     enddo
238     enddo
239     enddo
240    
241     call exf_swapffields_xz( tmpfldxz2, tmpfldxz, mythid)
242    
243     do bj = jtlo,jthi
244     do bi = itlo,ithi
245     do k = 1,nr
246     do i = imin,imax
247     xx_obcsn0(i,k,bi,bj,iobcs) = tmpfldxz2(i,k,bi,bj)
248     enddo
249     enddo
250     enddo
251     enddo
252    
253 heimbach 1.3 call active_read_xz_loc( fnameobcsn, tmpfldxz,
254 heimbach 1.2 & (obcsncount1-1)*nobcs+iobcs,
255     & doglobalread, ladinit, optimcycle,
256     & mythid, xx_obcsn_dummy )
257    
258     if ( optimcycle .gt. 0) then
259     if (iobcs .eq. 3) then
260     cgg Special attention is needed for the normal velocity.
261     cgg For the north, this is the v velocity, iobcs = 3.
262     cgg This is done on a columnwise basis here.
263     do bj = jtlo,jthi
264     do bi = itlo, ithi
265     do i = imin,imax
266    
267     cgg The barotropic velocity is stored in the level 1.
268     vbaro = tmpfldxz(i,1,bi,bj)
269     cgg Except for the special point which balances barotropic vol.flux.
270     cgg Special column in the NW corner.
271     j = OB_Jn(I,bi,bj)
272     if (ob_iw(j,bi,bj).eq.(i-1).and.
273     & ob_iw(j,bi,bj).ne. 0) then
274     print*,'correct vbaro',vbaro
275     print*,'Apply shiftvel2 @ i,j'
276     print*,shiftvel(2),i,j
277     vbaro = shiftvel(2)
278     endif
279     tmpfldxz(i,1,bi,bj) = 0.d0
280     vtop = 0.d0
281    
282     do k = 1,Nr
283     cgg If cells are not full, this should be modified with hFac.
284     cgg
285     cgg The xx field (tmpfldxz) does not contain the velocity at the
286     cgg surface level. This velocity is not independent; it must
287     cgg exactly balance the volume flux, since we are dealing with
288     cgg the baroclinic velocity structure..
289     vtop = tmpfldxz(i,k,bi,bj)*
290 heimbach 1.4 & maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop
291 heimbach 1.2 cgg Add the barotropic velocity component.
292     if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
293     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
294     endif
295     enddo
296     cgg Compute the baroclinic velocity at level 1. Should balance flux.
297     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
298 heimbach 1.4 & - vtop / delR(1)
299 heimbach 1.2 enddo
300     enddo
301     enddo
302     endif
303     if (iobcs .eq. 4) then
304     cgg Special attention is needed for the normal velocity.
305     cgg For the north, this is the v velocity, iobcs = 3.
306     cgg This is done on a columnwise basis here.
307     do bj = jtlo,jthi
308     do bi = itlo, ithi
309     do i = imin,imax
310    
311     cgg The barotropic velocity is stored in the level 1.
312     vbaro = tmpfldxz(i,1,bi,bj)
313     cgg Except for the special point which balances barotropic vol.flux.
314     cgg Special column in the NW corner.
315     j = OB_Jn(I,bi,bj)
316     tmpfldxz(i,1,bi,bj) = 0.d0
317     vtop = 0.d0
318    
319     do k = 1,Nr
320     cgg If cells are not full, this should be modified with hFac.
321     cgg
322     cgg The xx field (tmpfldxz) does not contain the velocity at the
323     cgg surface level. This velocity is not independent; it must
324     cgg exactly balance the volume flux, since we are dealing with
325     cgg the baroclinic velocity structure..
326     vtop = tmpfldxz(i,k,bi,bj)*
327 heimbach 1.4 & maskW(i,j,k,bi,bj) * delR(k) + vtop
328 heimbach 1.2 cgg Add the barotropic velocity component.
329     if (maskW(i,j,k,bi,bj) .ne. 0.) then
330     tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
331     endif
332     enddo
333     cgg Compute the baroclinic velocity at level 1. Should balance flux.
334     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
335 heimbach 1.4 & - vtop / delR(1)
336 heimbach 1.2 enddo
337     enddo
338     enddo
339     endif
340     endif
341    
342     do bj = jtlo,jthi
343     do bi = itlo,ithi
344     do k = 1,nr
345     do i = imin,imax
346     xx_obcsn1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
347     cgg & * maskxz (i,k,bi,bj)
348     enddo
349     enddo
350     enddo
351     enddo
352    
353     endif
354    
355     c-- Add control to model variable.
356     do bj = jtlo,jthi
357     do bi = itlo,ithi
358     c-- Calculate mask for tracer cells (0 => land, 1 => water).
359     do k = 1,nr
360     do i = 1,snx
361     j = OB_Jn(I,bi,bj)
362     if (iobcs .EQ. 1) then
363     OBNt(i,k,bi,bj) = OBNt (i,k,bi,bj)
364     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
365     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
366     OBNt(i,k,bi,bj) = OBNt(i,k,bi,bj)
367     & *maskS(i,j+jp1,k,bi,bj)
368     cgg & *maskxz(i,k,bi,bj)
369     else if (iobcs .EQ. 2) then
370     OBNs(i,k,bi,bj) = OBNs (i,k,bi,bj)
371     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
372     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
373     OBNs(i,k,bi,bj) = OBNs(i,k,bi,bj)
374     & *maskS(i,j+jp1,k,bi,bj)
375     cgg & *maskxz(i,k,bi,bj)
376     else if (iobcs .EQ. 4) then
377     OBNu(i,k,bi,bj) = OBNu (i,k,bi,bj)
378     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
379     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
380     OBNu(i,k,bi,bj) = OBNu(i,k,bi,bj)
381     & *maskW(i,j,k,bi,bj)
382     cgg & *maskxz(i,k,bi,bj)
383     else if (iobcs .EQ. 3) then
384     OBNv(i,k,bi,bj) = OBNv (i,k,bi,bj)
385     & + obcsnfac *xx_obcsn0(i,k,bi,bj,iobcs)
386     & + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
387     OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
388     & *maskS(i,j+jp1,k,bi,bj)
389     cgg & *maskxz(i,k,bi,bj)
390     endif
391     enddo
392     enddo
393     enddo
394     enddo
395    
396     C-- End over iobcs loop
397     enddo
398    
399     #else /* ALLOW_OBCSN_CONTROL undefined */
400    
401     c == routine arguments ==
402    
403     _RL mytime
404     integer myiter
405     integer mythid
406    
407     c-- CPP flag ALLOW_OBCSN_CONTROL undefined.
408    
409     #endif /* ALLOW_OBCSN_CONTROL */
410    
411     end
412    
413    
414    
415    
416    

  ViewVC Help
Powered by ViewVC 1.1.22