/[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.2.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 15  c     ================================== Line 15  c     ==================================
15  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
16  c     ==================================================================  c     ==================================================================
17  c  c
18  c     o Get norhtern obc of the control vector and add it  c     o Get western obc of the control vector and add it
19  c       to dyn. fields  c       to dyn. fields
20  c  c
21  c     started: heimbach@mit.edu, 29-Aug-2001  c     started: heimbach@mit.edu, 29-Aug-2001
# Line 56  c     == local variables == Line 56  c     == local variables ==
56        integer imin,imax        integer imin,imax
57        integer ilobcsw        integer ilobcsw
58        integer iobcs        integer iobcs
59          integer ip1
60    
61        _RL     dummy        _RL     dummy
62        _RL     obcswfac        _RL     obcswfac
# Line 64  c     == local variables == Line 65  c     == local variables ==
65        integer obcswcount0        integer obcswcount0
66        integer obcswcount1        integer obcswcount1
67    
68        _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
69    
70        logical doglobalread        logical doglobalread
71        logical ladinit        logical ladinit
72    
73        character*(80) fnameobcsw        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 ==  c     == external functions ==
83    
# Line 89  c     == end of interface == Line 95  c     == end of interface ==
95        jmax = sny+oly        jmax = sny+oly
96        imin = 1-olx        imin = 1-olx
97        imax = snx+olx        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.  c--   Now, read the control vector.
108        doglobalread = .false.        doglobalread = .false.
# Line 111  c--   Get the counters, flags, and the i Line 125  c--   Get the counters, flags, and the i
125       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
126    
127        do iobcs = 1,nobcs        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          call active_read_yz( 'maskobcsw', maskyz,            do bj = jtlo,jthi
195       &                       iobcs,              do bi = itlo,ithi
196       &                       doglobalread, ladinit, 0,                do k = 1,nr
197       &                       mythid, dummy )                  do j = jmin,jmax
198                      xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
199          call active_read_yz( fnameobcsw, tmpfldyz,                   enddo
200       &                       (obcswcount0-1)*nobcs+iobcs,                enddo
201       &                       doglobalread, ladinit, optimcycle,              enddo
202       &                       mythid, xx_obcsw_dummy )            enddo
203            endif
204    
205          do bj = jtlo,jthi          if ( (obcswfirst) .or. (obcswchanged)) then
206            do bi = itlo,ithi            
207              do k = 1,nr  cgg(    This is a terribly long way to do it. However, the dimensions don't exactly
208                do j = jmin,jmax  cgg     match up. I will blame Fortran for the ugliness.
209                  yz_obcs0(j,k,bi,bj)  = tmpfldyz (j,k,bi,bj)  
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                enddo
218              enddo              enddo
219            enddo            enddo
         enddo  
220    
221          call active_read_yz( fnameobcsw, tmpfldyz,            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
      &                       (obcswcount1-1)*nobcs+iobcs,  
      &                       doglobalread, ladinit, optimcycle,  
      &                       mythid, xx_obcsw_dummy )  
222    
223          do bj = jtlo,jthi            do bj = jtlo,jthi
224            do bi = itlo,ithi              do bi = itlo,ithi
225              do k = 1,nr                do k = 1,nr
226                do j = jmin,jmax                  do j = jmin,jmax
227                  yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj)                    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                enddo
297              enddo              enddo
298            enddo            enddo
299          enddo          endif
300    
301  c--     Add control to model variable.  c--     Add control to model variable.
302          do bj = jtlo,jthi          do bj = jtlo,jthi
# Line 153  c--     Add control to model variable. Line 304  c--     Add control to model variable.
304  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
305                do k = 1,nr                do k = 1,nr
306                   do j = 1,sny                   do j = 1,sny
307                        i = OB_Iw(j,bi,bj)
308                      if (iobcs .EQ. 1) then                      if (iobcs .EQ. 1) then
309                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
310       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
311       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
312                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
313       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
314                      else if (iobcs .EQ. 2) then                      else if (iobcs .EQ. 2) then
315                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
316       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
317       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
318                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
319       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
320                      else if (iobcs .EQ. 3) then                      else if (iobcs .EQ. 3) then
321                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
322       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
323       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
324                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
325       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
326                      else if (iobcs .EQ. 4) then                      else if (iobcs .EQ. 4) then
327                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
328       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
329       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
330                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
331       &                      *maskyz(j,k,bi,bj)       &                      *maskS(i,j,k,bi,bj)
332                      endif                      endif
333                   enddo                   enddo
334                enddo                enddo
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.1.2.2
  ViewVC Help
Powered by ViewVC 1.1.22