/[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.2 by heimbach, Thu May 30 19:56:44 2002 UTC revision 1.9 by mlosch, Wed Jan 19 08:42:06 2011 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "CTRL_CPPOPTIONS.h"  #include "CTRL_CPPOPTIONS.h"
5  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 20  c       to dyn. fields Line 22  c       to dyn. fields
22  c  c
23  c     started: heimbach@mit.edu, 29-Aug-2001  c     started: heimbach@mit.edu, 29-Aug-2001
24  c  c
25    c     modified: gebbie@mit.edu, 18-Mar-2003
26  c     ==================================================================  c     ==================================================================
27  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
28  c     ==================================================================  c     ==================================================================
# Line 72  cgg      _RL maskyz   (1-oly:sny+oly,nr, Line 75  cgg      _RL maskyz   (1-oly:sny+oly,nr,
75    
76        character*(80) fnameobcsw        character*(80) fnameobcsw
77    
78  #ifdef BALANCE_CONTROL_VOLFLUX  cgg(  Variables for splitting barotropic/baroclinic vels.
 cgg(  Variables for balancing volume flux.  
79        _RL ubaro        _RL ubaro
80        _RL ubarocount        _RL utop
81  cgg)  cgg)
 #endif  
82    
83  c     == external functions ==  c     == external functions ==
84    
# Line 97  c     == end of interface == Line 98  c     == end of interface ==
98        imax = snx+olx        imax = snx+olx
99        ip1  = 1        ip1  = 1
100    
 #ifdef BALANCE_CONTROL_VOLFLUX  
101  cgg(  Initialize variables for balancing volume flux.  cgg(  Initialize variables for balancing volume flux.
102        ubaro = 0.d0        ubaro = 0.d0
103        ubarocount = 0.d0        utop = 0.d0
104  cgg)  cgg)
 #endif  
105    
106  c--   Now, read the control vector.  c--   Now, read the control vector.
107        doglobalread = .false.        doglobalread = .false.
# Line 110  c--   Now, read the control vector. Line 109  c--   Now, read the control vector.
109    
110        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
111          ilobcsw=ilnblnk( xx_obcsw_file )          ilobcsw=ilnblnk( xx_obcsw_file )
112          write(fnameobcsw(1:80),'(2a,i10.10)')          write(fnameobcsw(1:80),'(2a,i10.10)')
113       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle
       else  
         print*  
         print*,' ctrl_getobcsw: optimcycle not set correctly.'  
         print*,' ... stopped in ctrl_getobcsw.'  
114        endif        endif
115    
116  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
117        call ctrl_GetRec( 'xx_obcsw',        call ctrl_get_gen_rec(
118         I                   xx_obcswstartdate, xx_obcswperiod,
119       O                   obcswfac, obcswfirst, obcswchanged,       O                   obcswfac, obcswfirst, obcswchanged,
120       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
121       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
122    
123        do iobcs = 1,nobcs        do iobcs = 1,nobcs
         
 cgg        if ( (obcswfirst) .or. (obcswchanged) ) then  
 cgg          call active_read_yz( 'maskobcsw', maskyz,  
 cgg     &                         iobcs,  
 cgg     &                         doglobalread, ladinit, 0,  
 cgg     &                         mythid, dummy )  
 cgg        endif  
   
124          if ( obcswfirst ) then          if ( obcswfirst ) then
125            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
126       &                         (obcswcount0-1)*nobcs+iobcs,       &                         (obcswcount0-1)*nobcs+iobcs,
127       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
128       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
129    
130  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
131            if ( optimcycle .gt. 0) then            
132              if ( optimcycle .gt. 0) then
133              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
134  cgg(        Make sure that the xx velocity field has a balanced net volume flux.  cgg         Special attention is needed for the normal velocity.
 cgg         Find the barotropic flow normal to the boundary.  
135  cgg         For the north, this is the v velocity, iobcs = 4.  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                do bj = jtlo,jthi
138                  do bi = itlo, ithi                  do bi = itlo, ithi
139                    do j = jmin,jmax                    do j = jmin,jmax
140                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
141                      ubaro = 0.d0  
142                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
143                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
144  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
145                        ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)                      utop = 0.d0
146       &                      * delZ(k) + ubaro  
147                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
148       &                           * delZ(k) +ubarocount  cgg    If cells are not full, this should be modified with hFac.
149                      enddo  cgg
150                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
151                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
152                      else  cgg    exactly balance the volume flux, since we are dealing with
153                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
154                      endif                        utop = tmpfldyz(j,k,bi,bj)*
155                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
156                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
157                          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                      enddo
161    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
162                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
163         &                                      - utop / delR(1)
164                    enddo                    enddo
165                  enddo                  enddo
166                enddo                enddo
167              endif              endif
168            endif              if (iobcs .eq. 4) then
169  cgg)  cgg         Special attention is needed for the normal velocity.
170  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
171  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
172                do bj = jtlo,jthi                do bj = jtlo,jthi
173                  do bi = itlo, ithi                  do bi = itlo, ithi
174                    do k = 1,Nr                    do j = jmin,jmax
175                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
176                        i = OB_Iw(j,bi,bj)  
177                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
178       &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)                      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    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    cgg    exactly balance the volume flux, since we are dealing with
188    cgg    the baroclinic velocity structure..
189                          utop = tmpfldyz(j,k,bi,bj)*
190         &                maskS(i,j,k,bi,bj) * delR(k) + utop
191    cgg    Add the barotropic velocity component.
192                          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                      enddo
196    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
197                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
198         &                                      - utop / delR(1)
199                    enddo                    enddo
200                  enddo                  enddo
201                enddo                enddo
202              endif              endif
203            endif            endif
204  #endif  
205    #endif /* ALLOW_CTRL_OBCS_BALANCE */
206    
207            do bj = jtlo,jthi            do bj = jtlo,jthi
208              do bi = itlo,ithi              do bi = itlo,ithi
209                do k = 1,nr                do k = 1,nr
210                  do j = jmin,jmax                  do j = jmin,jmax
211                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
212    cgg     &                                        *   maskyz (j,k,bi,bj)
213                   enddo                   enddo
214                enddo                enddo
215              enddo              enddo
# Line 203  cgg) Line 217  cgg)
217          endif          endif
218    
219          if ( (obcswfirst) .or. (obcswchanged)) then          if ( (obcswfirst) .or. (obcswchanged)) then
             
 cgg(    This is a terribly long way to do it. However, the dimensions don't exactly  
 cgg     match up. I will blame Fortran for the ugliness.  
220    
221            do bj = jtlo,jthi            do bj = jtlo,jthi
222              do bi = itlo,ithi             do bi = itlo,ithi
223                do k = 1,nr              do k = 1,nr
224                  do j = jmin,jmax               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) = xx_obcsw1(j,k,bi,bj,iobcs)                tmpfldyz (j,k,bi,bj)       = 0. _d 0
227                  enddo               enddo
               enddo  
             enddo  
           enddo  
   
           call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)          
                 enddo  
               enddo  
228              enddo              enddo
229               enddo
230            enddo            enddo
231    
232            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
233       &                         (obcswcount1-1)*nobcs+iobcs,       &                         (obcswcount1-1)*nobcs+iobcs,
234       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
235       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
236    
237  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
238            if (optimcycle .gt. 0) then  
239              if ( optimcycle .gt. 0) then
240              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
241  cgg(        Make sure that the xx velocity field has a balanced net volume flux.  cgg         Special attention is needed for the normal velocity.
 cgg         Find the barotropic flow normal to the boundary.  
242  cgg         For the north, this is the v velocity, iobcs = 4.  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                do bj = jtlo,jthi
245                  do bi = itlo, ithi                  do bi = itlo, ithi
246                    do j = jmin,jmax                    do j = jmin,jmax
247                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
248                      ubaro = 0.d0  
249                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
250                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
251  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
252                        ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)                      utop = 0.d0
253       &                      * delZ(k) + ubaro  
254                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
255       &                           * delZ(k) + ubarocount  cgg    If cells are not full, this should be modified with hFac.
256                      enddo  cgg
257                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
258                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
259                      else  cgg    exactly balance the volume flux, since we are dealing with
260                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
261                      endif                        utop = tmpfldyz(j,k,bi,bj)*
262                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
263                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
264                          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                      enddo
268    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
269                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
270         &                                      - utop / delR(1)
271                    enddo                    enddo
272                  enddo                  enddo
273                enddo                enddo
274              endif              endif
275            endif              if (iobcs .eq. 4) then
276  cgg)  cgg         Special attention is needed for the normal velocity.
277  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
278  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
279                do bj = jtlo,jthi                do bj = jtlo,jthi
280                  do bi = itlo, ithi                  do bi = itlo, ithi
281                    do k = 1,Nr                    do j = jmin,jmax
282                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
283                        i = OB_Iw(j,bi,bj)  
284                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
285       &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)                      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    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    cgg    exactly balance the volume flux, since we are dealing with
295    cgg    the baroclinic velocity structure..
296                          utop = tmpfldyz(j,k,bi,bj)*
297         &                maskS(i,j,k,bi,bj) * delR(k) + utop
298    cgg    Add the barotropic velocity component.
299                          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                      enddo
303    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
304                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
305         &                                      - utop / delR(1)
306                    enddo                    enddo
307                  enddo                  enddo
308                enddo                enddo
309              endif              endif
310            endif            endif
311  #endif  
312    #endif /* ALLOW_CTRL_OBCS_BALANCE */
313    
314            do bj = jtlo,jthi            do bj = jtlo,jthi
315              do bi = itlo,ithi              do bi = itlo,ithi
316                do k = 1,nr                do k = 1,nr
317                  do j = jmin,jmax                  do j = jmin,jmax
318                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
319    cgg     &                                        *   maskyz (j,k,bi,bj)
320                   enddo                   enddo
321                enddo                enddo
322              enddo              enddo
# Line 299  cgg) Line 324  cgg)
324          endif          endif
325    
326  c--     Add control to model variable.  c--     Add control to model variable.
327          do bj = jtlo,jthi          do bj = jtlo, jthi
328             do bi = itlo,ithi             do bi = itlo, ithi
329  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
330                do k = 1,nr                do k = 1,nr
331                   do j = 1,sny                   do j = 1,sny
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22