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

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

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


Revision 1.9 - (hide annotations) (download)
Wed Jan 19 08:42:06 2011 UTC (14 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint62t, checkpoint62s, checkpoint62r
Changes since 1.8: +8 -22 lines
simplify code so that it does not need to use S/R exf_swapfields

1 mlosch 1.9 C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_getobcsw.F,v 1.8 2007/10/09 00:00:00 jmc Exp $
2 jmc 1.8 C $Name: $
3 heimbach 1.2
4     #include "CTRL_CPPOPTIONS.h"
5     #ifdef ALLOW_OBCS
6     # include "OBCS_OPTIONS.h"
7     #endif
8    
9    
10     subroutine ctrl_getobcsw(
11     I mytime,
12     I myiter,
13     I mythid
14     & )
15    
16     c ==================================================================
17     c SUBROUTINE ctrl_getobcsw
18     c ==================================================================
19     c
20     c o Get western obc of the control vector and add it
21     c to dyn. fields
22     c
23     c started: heimbach@mit.edu, 29-Aug-2001
24     c
25     c modified: gebbie@mit.edu, 18-Mar-2003
26     c ==================================================================
27     c SUBROUTINE ctrl_getobcsw
28     c ==================================================================
29    
30     implicit none
31    
32     #ifdef ALLOW_OBCSW_CONTROL
33    
34     c == global variables ==
35    
36     #include "EEPARAMS.h"
37     #include "SIZE.h"
38     #include "PARAMS.h"
39     #include "GRID.h"
40     #include "OBCS.h"
41    
42     #include "ctrl.h"
43     #include "ctrl_dummy.h"
44     #include "optim.h"
45    
46     c == routine arguments ==
47    
48     _RL mytime
49     integer myiter
50     integer mythid
51    
52     c == local variables ==
53    
54     integer bi,bj
55     integer i,j,k
56     integer itlo,ithi
57     integer jtlo,jthi
58     integer jmin,jmax
59     integer imin,imax
60     integer ilobcsw
61     integer iobcs
62     integer ip1
63    
64     _RL dummy
65     _RL obcswfac
66     logical obcswfirst
67     logical obcswchanged
68     integer obcswcount0
69     integer obcswcount1
70    
71     cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy)
72    
73     logical doglobalread
74     logical ladinit
75    
76     character*(80) fnameobcsw
77    
78     cgg( Variables for splitting barotropic/baroclinic vels.
79     _RL ubaro
80     _RL utop
81     cgg)
82    
83     c == external functions ==
84    
85     integer ilnblnk
86     external ilnblnk
87    
88    
89     c == end of interface ==
90    
91     jtlo = mybylo(mythid)
92     jthi = mybyhi(mythid)
93     itlo = mybxlo(mythid)
94     ithi = mybxhi(mythid)
95     jmin = 1-oly
96     jmax = sny+oly
97     imin = 1-olx
98     imax = snx+olx
99     ip1 = 1
100    
101     cgg( Initialize variables for balancing volume flux.
102     ubaro = 0.d0
103     utop = 0.d0
104     cgg)
105    
106     c-- Now, read the control vector.
107     doglobalread = .false.
108     ladinit = .false.
109    
110     if (optimcycle .ge. 0) then
111     ilobcsw=ilnblnk( xx_obcsw_file )
112 jmc 1.8 write(fnameobcsw(1:80),'(2a,i10.10)')
113 heimbach 1.2 & xx_obcsw_file(1:ilobcsw), '.', optimcycle
114     endif
115    
116     c-- Get the counters, flags, and the interpolation factor.
117     call ctrl_get_gen_rec(
118     I xx_obcswstartdate, xx_obcswperiod,
119     O obcswfac, obcswfirst, obcswchanged,
120     O obcswcount0,obcswcount1,
121     I mytime, myiter, mythid )
122    
123     do iobcs = 1,nobcs
124     if ( obcswfirst ) then
125 jmc 1.8 call active_read_yz( fnameobcsw, tmpfldyz,
126 heimbach 1.2 & (obcswcount0-1)*nobcs+iobcs,
127     & doglobalread, ladinit, optimcycle,
128     & mythid, xx_obcsw_dummy )
129    
130 mlosch 1.6 #ifdef ALLOW_CTRL_OBCS_BALANCE
131    
132 jmc 1.8 if ( optimcycle .gt. 0) then
133 heimbach 1.2 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 j = jmin,jmax
140     i = OB_Iw(J,bi,bj)
141    
142     cgg The barotropic velocity is stored in the level 1.
143     ubaro = tmpfldyz(j,1,bi,bj)
144     tmpfldyz(j,1,bi,bj) = 0.d0
145     utop = 0.d0
146    
147     do k = 1,Nr
148     cgg If cells are not full, this should be modified with hFac.
149 jmc 1.8 cgg
150     cgg The xx field (tmpfldxz) does not contain the velocity at the
151     cgg surface level. This velocity is not independent; it must
152 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
153 jmc 1.8 cgg the baroclinic velocity structure..
154 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
155 heimbach 1.5 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
156 jmc 1.8 cgg Add the barotropic velocity component.
157 heimbach 1.2 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
158     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
159     endif
160     enddo
161     cgg Compute the baroclinic velocity at level 1. Should balance flux.
162 jmc 1.8 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
163 heimbach 1.5 & - utop / delR(1)
164 heimbach 1.2 enddo
165     enddo
166     enddo
167     endif
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 j = jmin,jmax
175     i = OB_Iw(J,bi,bj)
176    
177     cgg The barotropic velocity is stored in the level 1.
178     ubaro = tmpfldyz(j,1,bi,bj)
179     tmpfldyz(j,1,bi,bj) = 0.d0
180     utop = 0.d0
181    
182     do k = 1,Nr
183     cgg If cells are not full, this should be modified with hFac.
184 jmc 1.8 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 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
188 jmc 1.8 cgg the baroclinic velocity structure..
189 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
190 heimbach 1.5 & maskS(i,j,k,bi,bj) * delR(k) + utop
191 jmc 1.8 cgg Add the barotropic velocity component.
192 heimbach 1.2 if (maskS(i,j,k,bi,bj) .ne. 0.) then
193     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
194     endif
195     enddo
196     cgg Compute the baroclinic velocity at level 1. Should balance flux.
197 jmc 1.8 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
198 heimbach 1.5 & - utop / delR(1)
199 heimbach 1.2 enddo
200     enddo
201     enddo
202     endif
203     endif
204    
205 mlosch 1.6 #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 j = jmin,jmax
211     xx_obcsw1(j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
212     cgg & * maskyz (j,k,bi,bj)
213     enddo
214     enddo
215     enddo
216     enddo
217     endif
218    
219     if ( (obcswfirst) .or. (obcswchanged)) then
220 jmc 1.8
221 heimbach 1.2 do bj = jtlo,jthi
222 mlosch 1.9 do bi = itlo,ithi
223     do k = 1,nr
224     do j = jmin,jmax
225     xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
226     tmpfldyz (j,k,bi,bj) = 0. _d 0
227     enddo
228 heimbach 1.2 enddo
229 mlosch 1.9 enddo
230 heimbach 1.2 enddo
231    
232 jmc 1.8 call active_read_yz( fnameobcsw, tmpfldyz,
233 heimbach 1.2 & (obcswcount1-1)*nobcs+iobcs,
234     & doglobalread, ladinit, optimcycle,
235     & mythid, xx_obcsw_dummy )
236    
237 mlosch 1.6 #ifdef ALLOW_CTRL_OBCS_BALANCE
238    
239 jmc 1.8 if ( optimcycle .gt. 0) then
240 heimbach 1.2 if (iobcs .eq. 3) then
241     cgg Special attention is needed for the normal velocity.
242     cgg For the north, this is the v velocity, iobcs = 4.
243     cgg This is done on a columnwise basis here.
244     do bj = jtlo,jthi
245     do bi = itlo, ithi
246     do j = jmin,jmax
247     i = OB_Iw(J,bi,bj)
248    
249     cgg The barotropic velocity is stored in the level 1.
250     ubaro = tmpfldyz(j,1,bi,bj)
251     tmpfldyz(j,1,bi,bj) = 0.d0
252     utop = 0.d0
253    
254     do k = 1,Nr
255     cgg If cells are not full, this should be modified with hFac.
256 jmc 1.8 cgg
257     cgg The xx field (tmpfldxz) does not contain the velocity at the
258     cgg surface level. This velocity is not independent; it must
259 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
260 jmc 1.8 cgg the baroclinic velocity structure..
261 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
262 heimbach 1.5 & maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
263 jmc 1.8 cgg Add the barotropic velocity component.
264 heimbach 1.2 if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
265     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
266     endif
267     enddo
268     cgg Compute the baroclinic velocity at level 1. Should balance flux.
269 jmc 1.8 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
270 heimbach 1.5 & - utop / delR(1)
271 heimbach 1.2 enddo
272     enddo
273     enddo
274     endif
275     if (iobcs .eq. 4) then
276     cgg Special attention is needed for the normal velocity.
277     cgg For the north, this is the v velocity, iobcs = 4.
278     cgg This is done on a columnwise basis here.
279     do bj = jtlo,jthi
280     do bi = itlo, ithi
281     do j = jmin,jmax
282     i = OB_Iw(J,bi,bj)
283    
284     cgg The barotropic velocity is stored in the level 1.
285     ubaro = tmpfldyz(j,1,bi,bj)
286     tmpfldyz(j,1,bi,bj) = 0.d0
287     utop = 0.d0
288    
289     do k = 1,Nr
290     cgg If cells are not full, this should be modified with hFac.
291 jmc 1.8 cgg
292     cgg The xx field (tmpfldxz) does not contain the velocity at the
293     cgg surface level. This velocity is not independent; it must
294 heimbach 1.2 cgg exactly balance the volume flux, since we are dealing with
295 jmc 1.8 cgg the baroclinic velocity structure..
296 heimbach 1.2 utop = tmpfldyz(j,k,bi,bj)*
297 heimbach 1.5 & maskS(i,j,k,bi,bj) * delR(k) + utop
298 jmc 1.8 cgg Add the barotropic velocity component.
299 heimbach 1.2 if (maskS(i,j,k,bi,bj) .ne. 0.) then
300     tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
301     endif
302     enddo
303     cgg Compute the baroclinic velocity at level 1. Should balance flux.
304 jmc 1.8 tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
305 heimbach 1.5 & - utop / delR(1)
306 heimbach 1.2 enddo
307     enddo
308     enddo
309     endif
310     endif
311    
312 mlosch 1.6 #endif /* ALLOW_CTRL_OBCS_BALANCE */
313    
314 heimbach 1.2 do bj = jtlo,jthi
315     do bi = itlo,ithi
316     do k = 1,nr
317     do j = jmin,jmax
318     xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
319     cgg & * maskyz (j,k,bi,bj)
320     enddo
321     enddo
322     enddo
323     enddo
324     endif
325    
326     c-- Add control to model variable.
327     do bj = jtlo, jthi
328     do bi = itlo, ithi
329     c-- Calculate mask for tracer cells (0 => land, 1 => water).
330     do k = 1,nr
331     do j = 1,sny
332     i = OB_Iw(j,bi,bj)
333     if (iobcs .EQ. 1) then
334     OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
335     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
336     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
337     OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
338     & *maskW(i+ip1,j,k,bi,bj)
339     else if (iobcs .EQ. 2) then
340     OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
341     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
342     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
343     OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
344     & *maskW(i+ip1,j,k,bi,bj)
345     else if (iobcs .EQ. 3) then
346     OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
347     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
348     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
349     OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
350     & *maskW(i+ip1,j,k,bi,bj)
351     else if (iobcs .EQ. 4) then
352     OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
353     & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs)
354     & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
355     OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
356     & *maskS(i,j,k,bi,bj)
357     endif
358     enddo
359     enddo
360     enddo
361     enddo
362    
363     C-- End over iobcs loop
364     enddo
365    
366     #else /* ALLOW_OBCSW_CONTROL undefined */
367    
368     c == routine arguments ==
369    
370     _RL mytime
371     integer myiter
372     integer mythid
373    
374     c-- CPP flag ALLOW_OBCSW_CONTROL undefined.
375    
376     #endif /* ALLOW_OBCSW_CONTROL */
377    
378     end
379    

  ViewVC Help
Powered by ViewVC 1.1.22