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

Diff of /MITgcm/pkg/ctrl/ctrl_getobcse.F

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

revision 1.1 by heimbach, Tue Feb 5 20:23:58 2002 UTC revision 1.2 by heimbach, Tue Jun 24 16:07:06 2003 UTC
# Line 0  Line 1 
1    
2    #include "CTRL_CPPOPTIONS.h"
3    #ifdef ALLOW_OBCS
4    # include "OBCS_OPTIONS.h"
5    #endif
6    
7    
8          subroutine ctrl_getobcse(
9         I                             mytime,
10         I                             myiter,
11         I                             mythid
12         &                           )
13    
14    c     ==================================================================
15    c     SUBROUTINE ctrl_getobcse
16    c     ==================================================================
17    c
18    c     o Get eastern 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     ==================================================================
24    c     SUBROUTINE ctrl_getobcse
25    c     ==================================================================
26    
27          implicit none
28    
29    #ifdef ALLOW_OBCSE_CONTROL
30    
31    c     == global variables ==
32    
33    #include "EEPARAMS.h"
34    #include "SIZE.h"
35    #include "PARAMS.h"
36    #include "GRID.h"
37    #include "OBCS.h"
38    
39    #include "ctrl.h"
40    #include "ctrl_dummy.h"
41    #include "optim.h"
42    
43    c     == routine arguments ==
44    
45          _RL     mytime
46          integer myiter
47          integer mythid
48    
49    c     == local variables ==
50    
51          integer bi,bj
52          integer i,j,k
53          integer itlo,ithi
54          integer jtlo,jthi
55          integer jmin,jmax
56          integer imin,imax
57          integer ilobcse
58          integer iobcs
59    
60          _RL     dummy
61          _RL     obcsefac
62          logical obcsefirst
63          logical obcsechanged
64          integer obcsecount0
65          integer obcsecount1
66          integer ip1
67    
68    cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
69    
70          logical doglobalread
71          logical ladinit
72    
73          character*(80) fnameobcse
74    
75    cgg(  Variables for splitting barotropic/baroclinic vels.
76          _RL ubaro
77          _RL utop
78    cgg)
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          jmin = 1-oly
93          jmax = sny+oly
94          imin = 1-olx
95          imax = snx+olx
96          ip1  = 0
97    
98    cgg(  Initialize variables for balancing volume flux.
99          ubaro = 0.d0
100          utop = 0.d0
101    cgg)
102    
103    c--   Now, read the control vector.
104          doglobalread = .false.
105          ladinit      = .false.
106    
107          if (optimcycle .ge. 0) then
108            ilobcse=ilnblnk( xx_obcse_file )
109            write(fnameobcse(1:80),'(2a,i10.10)')
110         &       xx_obcse_file(1:ilobcse), '.', optimcycle
111          endif
112    
113    c--   Get the counters, flags, and the interpolation factor.
114          call ctrl_get_gen_rec(
115         I                   xx_obcsestartdate, xx_obcseperiod,
116         O                   obcsefac, obcsefirst, obcsechanged,
117         O                   obcsecount0,obcsecount1,
118         I                   mytime, myiter, mythid )
119    
120          do iobcs = 1,nobcs
121    
122            if ( obcsefirst ) then
123              call active_read_yz( fnameobcse, tmpfldyz,
124         &                         (obcsecount0-1)*nobcs+iobcs,
125         &                         doglobalread, ladinit, optimcycle,
126         &                         mythid, xx_obcse_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_Ie(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) * delZ(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 / delZ(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_Ie(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) * delZ(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 / delZ(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_obcse1(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 ( (obcsefirst) .or. (obcsechanged)) then
214              
215              do bj = jtlo,jthi
216                do bi = itlo,ithi
217                  do j = jmin,jmax
218                    do k = 1,nr
219                      tmpfldyz(j,k,bi,bj) = xx_obcse1(j,k,bi,bj,iobcs)
220                    enddo
221                  enddo
222                enddo
223              enddo
224    
225              call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
226    
227              do bj = jtlo,jthi
228                do bi = itlo,ithi
229                  do j = jmin,jmax
230                    do k = 1,nr
231                      xx_obcse0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
232                    enddo
233                  enddo
234                enddo
235              enddo
236    
237              call active_read_yz( fnameobcse, tmpfldyz,
238         &                         (obcsecount1-1)*nobcs+iobcs,
239         &                         doglobalread, ladinit, optimcycle,
240         &                         mythid, xx_obcse_dummy )
241    
242              if ( optimcycle .gt. 0) then          
243                if (iobcs .eq. 3) then
244    cgg         Special attention is needed for the normal velocity.
245    cgg         For the north, this is the v velocity, iobcs = 4.
246    cgg         This is done on a columnwise basis here.
247                  do bj = jtlo,jthi
248                    do bi = itlo, ithi
249                      do j = jmin,jmax
250                        i = OB_Ie(J,bi,bj)
251    
252    cgg         The barotropic velocity is stored in the level 1.
253                        ubaro = tmpfldyz(j,1,bi,bj)
254                        tmpfldyz(j,1,bi,bj) = 0.d0
255                        utop = 0.d0
256    
257                        do k = 1,Nr
258    cgg    If cells are not full, this should be modified with hFac.
259    cgg    
260    cgg    The xx field (tmpfldxz) does not contain the velocity at the
261    cgg    surface level. This velocity is not independent; it must
262    cgg    exactly balance the volume flux, since we are dealing with
263    cgg    the baroclinic velocity structure..
264                          utop = tmpfldyz(j,k,bi,bj)*
265         &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop
266    cgg    Add the barotropic velocity component.
267                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
268                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
269                          endif
270                        enddo
271    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
272                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
273         &                                      - utop / delZ(1)
274                      enddo
275                    enddo
276                  enddo
277                endif
278                if (iobcs .eq. 4) then
279    cgg         Special attention is needed for the normal velocity.
280    cgg         For the north, this is the v velocity, iobcs = 4.
281    cgg         This is done on a columnwise basis here.
282                  do bj = jtlo,jthi
283                    do bi = itlo, ithi
284                      do j = jmin,jmax
285                        i = OB_Ie(J,bi,bj)
286    
287    cgg         The barotropic velocity is stored in the level 1.
288                        ubaro = tmpfldyz(j,1,bi,bj)
289                        tmpfldyz(j,1,bi,bj) = 0.d0
290                        utop = 0.d0
291    
292                        do k = 1,Nr
293    cgg    If cells are not full, this should be modified with hFac.
294    cgg    
295    cgg    The xx field (tmpfldxz) does not contain the velocity at the
296    cgg    surface level. This velocity is not independent; it must
297    cgg    exactly balance the volume flux, since we are dealing with
298    cgg    the baroclinic velocity structure..
299                          utop = tmpfldyz(j,k,bi,bj)*
300         &                maskS(i,j,k,bi,bj) * delZ(k) + utop
301    cgg    Add the barotropic velocity component.
302                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
303                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
304                          endif
305                        enddo
306    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
307                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
308         &                                      - utop / delZ(1)
309                      enddo
310                    enddo
311                  enddo
312                endif
313              endif
314    
315              do bj = jtlo,jthi
316                do bi = itlo,ithi
317                  do k = 1,nr
318                    do j = jmin,jmax
319                      xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
320    cgg     &                                        *   maskyz (j,k,bi,bj)
321                     enddo
322                  enddo
323                enddo
324              enddo
325            endif
326    
327    c--     Add control to model variable.
328            do bj = jtlo,jthi
329               do bi = itlo,ithi
330    c--        Calculate mask for tracer cells (0 => land, 1 => water).
331                  do k = 1,nr
332                     do j = 1,sny
333                        i = OB_Ie(j,bi,bj)
334                        if (iobcs .EQ. 1) then
335                           OBEt(j,k,bi,bj) = OBEt (j,k,bi,bj)
336         &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
337         &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
338                           OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
339         &                      *maskW(i+ip1,j,k,bi,bj)
340                        else if (iobcs .EQ. 2) then
341                           OBEs(j,k,bi,bj) = OBEs (j,k,bi,bj)
342         &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
343         &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
344                           OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
345         &                      *maskW(i+ip1,j,k,bi,bj)
346                        else if (iobcs .EQ. 3) then
347                           OBEu(j,k,bi,bj) = OBEu (j,k,bi,bj)
348         &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
349         &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
350                           OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
351         &                      *maskW(i+ip1,j,k,bi,bj)
352                        else if (iobcs .EQ. 4) then
353                           OBEv(j,k,bi,bj) = OBEv (j,k,bi,bj)
354         &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
355         &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
356                           OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
357         &                      *maskS(i,j,k,bi,bj)
358                        endif
359                     enddo
360                  enddo
361               enddo
362            enddo
363    
364    C--   End over iobcs loop
365          enddo
366    
367    #else /* ALLOW_OBCSE_CONTROL undefined */
368    
369    c     == routine arguments ==
370    
371          _RL     mytime
372          integer myiter
373          integer mythid
374    
375    c--   CPP flag ALLOW_OBCSE_CONTROL undefined.
376    
377    #endif /* ALLOW_OBCSE_CONTROL */
378    
379          end
380    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22