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

Diff of /MITgcm/pkg/ctrl/ctrl_getobcsw.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.1.2.2 by heimbach, Thu May 30 19:56:44 2002 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_getobcsw(
9         I                             mytime,
10         I                             myiter,
11         I                             mythid
12         &                           )
13    
14    c     ==================================================================
15    c     SUBROUTINE ctrl_getobcsw
16    c     ==================================================================
17    c
18    c     o Get western 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_getobcsw
25    c     ==================================================================
26    
27          implicit none
28    
29    #ifdef ALLOW_OBCSW_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 ilobcsw
58          integer iobcs
59          integer ip1
60    
61          _RL     dummy
62          _RL     obcswfac
63          logical obcswfirst
64          logical obcswchanged
65          integer obcswcount0
66          integer obcswcount1
67    
68    cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
69    
70          logical doglobalread
71          logical ladinit
72    
73          character*(80) fnameobcsw
74    
75    #ifdef BALANCE_CONTROL_VOLFLUX
76    cgg(  Variables for balancing volume flux.
77          _RL ubaro
78          _RL ubarocount
79    cgg)
80    #endif
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          ip1  = 1
99    
100    #ifdef BALANCE_CONTROL_VOLFLUX
101    cgg(  Initialize variables for balancing volume flux.
102          ubaro = 0.d0
103          ubarocount = 0.d0
104    cgg)
105    #endif
106    
107    c--   Now, read the control vector.
108          doglobalread = .false.
109          ladinit      = .false.
110    
111          if (optimcycle .ge. 0) then
112            ilobcsw=ilnblnk( xx_obcsw_file )
113            write(fnameobcsw(1:80),'(2a,i10.10)')
114         &       xx_obcsw_file(1:ilobcsw), '.', optimcycle
115          else
116            print*
117            print*,' ctrl_getobcsw: optimcycle not set correctly.'
118            print*,' ... stopped in ctrl_getobcsw.'
119          endif
120    
121    c--   Get the counters, flags, and the interpolation factor.
122          call ctrl_GetRec( 'xx_obcsw',
123         O                   obcswfac, obcswfirst, obcswchanged,
124         O                   obcswcount0,obcswcount1,
125         I                   mytime, myiter, mythid )
126    
127          do iobcs = 1,nobcs
128          
129    cgg        if ( (obcswfirst) .or. (obcswchanged) ) then
130    cgg          call active_read_yz( 'maskobcsw', maskyz,
131    cgg     &                         iobcs,
132    cgg     &                         doglobalread, ladinit, 0,
133    cgg     &                         mythid, dummy )
134    cgg        endif
135    
136            if ( obcswfirst ) then
137              call active_read_yz( fnameobcsw, tmpfldyz,
138         &                         (obcswcount0-1)*nobcs+iobcs,
139         &                         doglobalread, ladinit, optimcycle,
140         &                         mythid, xx_obcsw_dummy )
141    
142    #ifdef BALANCE_CONTROL_VOLFLUX
143              if ( optimcycle .gt. 0) then          
144                if (iobcs .eq. 3) then
145    cgg(        Make sure that the xx velocity field has a balanced net volume flux.
146    cgg         Find the barotropic flow normal to the boundary.
147    cgg         For the north, this is the v velocity, iobcs = 4.
148                  do bj = jtlo,jthi
149                    do bi = itlo, ithi
150                      do j = jmin,jmax
151                        i = OB_Iw(j,bi,bj)
152                        ubaro = 0.d0
153                        ubarocount = 0.d0
154                        do k = 1,nr
155    cgg(   If cells are not full, this should be modified with hFac.
156                          ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)
157         &                      * delZ(k) + ubaro
158                          ubarocount = maskW(i+ip1,j,k,bi,bj)
159         &                           * delZ(k) +ubarocount
160                        enddo
161                        if (ubarocount .eq. 0.) then
162                          ubaro = 0.
163                        else
164                          ubaro = ubaro / ubarocount
165                        endif
166                        do k = 1,nr
167                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro
168                        enddo
169                      enddo
170                    enddo
171                  enddo
172                endif
173              endif
174    cgg)
175    #endif
176    #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL
177              if (optimcycle .gt. 0) then
178                if ( iobcs .eq. 3) then
179                  do bj = jtlo,jthi
180                    do bi = itlo, ithi
181                      do k = 1,Nr
182                        do j = jmin,jmax
183                          i = OB_Iw(j,bi,bj)
184                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)
185         &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)
186                        enddo
187                      enddo
188                    enddo
189                  enddo
190                endif
191              endif
192    #endif
193    
194              do bj = jtlo,jthi
195                do bi = itlo,ithi
196                  do k = 1,nr
197                    do j = jmin,jmax
198                      xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
199                     enddo
200                  enddo
201                enddo
202              enddo
203            endif
204    
205            if ( (obcswfirst) .or. (obcswchanged)) then
206              
207    cgg(    This is a terribly long way to do it. However, the dimensions don't exactly
208    cgg     match up. I will blame Fortran for the ugliness.
209    
210              do bj = jtlo,jthi
211                do bi = itlo,ithi
212                  do k = 1,nr
213                    do j = jmin,jmax
214    
215                      tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
216                    enddo
217                  enddo
218                enddo
219              enddo
220    
221              call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
222    
223              do bj = jtlo,jthi
224                do bi = itlo,ithi
225                  do k = 1,nr
226                    do j = jmin,jmax
227                      xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)        
228                    enddo
229                  enddo
230                enddo
231              enddo
232    
233              call active_read_yz( fnameobcsw, tmpfldyz,
234         &                         (obcswcount1-1)*nobcs+iobcs,
235         &                         doglobalread, ladinit, optimcycle,
236         &                         mythid, xx_obcsw_dummy )
237    
238    #ifdef BALANCE_CONTROL_VOLFLUX
239              if (optimcycle .gt. 0) then
240                if (iobcs .eq. 3) then
241    cgg(        Make sure that the xx velocity field has a balanced net volume flux.
242    cgg         Find the barotropic flow normal to the boundary.
243    cgg         For the north, this is the v velocity, iobcs = 4.
244                  do bj = jtlo,jthi
245                    do bi = itlo, ithi
246                      do j = jmin,jmax
247                        i = OB_Iw(j,bi,bj)
248                        ubaro = 0.d0
249                        ubarocount = 0.d0
250                        do k = 1,nr
251    cgg(   If cells are not full, this should be modified with hFac.
252                          ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)
253         &                      * delZ(k) + ubaro
254                          ubarocount = maskW(i+ip1,j,k,bi,bj)
255         &                           * delZ(k) + ubarocount
256                        enddo
257                        if (ubarocount .eq. 0.) then
258                          ubaro = 0.
259                        else
260                          ubaro = ubaro / ubarocount
261                        endif
262                        do k = 1,nr
263                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro
264                        enddo
265                      enddo
266                    enddo
267                  enddo
268                endif
269              endif
270    cgg)
271    #endif
272    #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL
273              if (optimcycle .gt. 0) then
274                if ( iobcs .eq. 3) then
275                  do bj = jtlo,jthi
276                    do bi = itlo, ithi
277                      do k = 1,Nr
278                        do j = jmin,jmax
279                          i = OB_Iw(j,bi,bj)
280                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)
281         &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)
282                        enddo
283                      enddo
284                    enddo
285                  enddo
286                endif
287              endif
288    #endif
289    
290              do bj = jtlo,jthi
291                do bi = itlo,ithi
292                  do k = 1,nr
293                    do j = jmin,jmax
294                      xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
295                     enddo
296                  enddo
297                enddo
298              enddo
299            endif
300    
301    c--     Add control to model variable.
302            do bj = jtlo,jthi
303               do bi = itlo,ithi
304    c--        Calculate mask for tracer cells (0 => land, 1 => water).
305                  do k = 1,nr
306                     do j = 1,sny
307                        i = OB_Iw(j,bi,bj)
308                        if (iobcs .EQ. 1) then
309                           OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
310         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
311         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
312                           OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
313         &                      *maskW(i+ip1,j,k,bi,bj)
314                        else if (iobcs .EQ. 2) then
315                           OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
316         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
317         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
318                           OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
319         &                      *maskW(i+ip1,j,k,bi,bj)
320                        else if (iobcs .EQ. 3) then
321                           OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
322         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
323         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
324                           OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
325         &                      *maskW(i+ip1,j,k,bi,bj)
326                        else if (iobcs .EQ. 4) then
327                           OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
328         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
329         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
330                           OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
331         &                      *maskS(i,j,k,bi,bj)
332                        endif
333                     enddo
334                  enddo
335               enddo
336            enddo
337    
338    C--   End over iobcs loop
339          enddo
340    
341    #else /* ALLOW_OBCSW_CONTROL undefined */
342    
343    c     == routine arguments ==
344    
345          _RL     mytime
346          integer myiter
347          integer mythid
348    
349    c--   CPP flag ALLOW_OBCSW_CONTROL undefined.
350    
351    #endif /* ALLOW_OBCSW_CONTROL */
352    
353          end
354    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.2.2
  ViewVC Help
Powered by ViewVC 1.1.22