/[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.7 by heimbach, Mon May 14 22:02:33 2007 UTC
# Line 20  c       to dyn. fields Line 20  c       to dyn. fields
20  c  c
21  c     started: heimbach@mit.edu, 29-Aug-2001  c     started: heimbach@mit.edu, 29-Aug-2001
22  c  c
23    c     modified: gebbie@mit.edu, 18-Mar-2003
24  c     ==================================================================  c     ==================================================================
25  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
26  c     ==================================================================  c     ==================================================================
# Line 72  cgg      _RL maskyz   (1-oly:sny+oly,nr, Line 73  cgg      _RL maskyz   (1-oly:sny+oly,nr,
73    
74        character*(80) fnameobcsw        character*(80) fnameobcsw
75    
76  #ifdef BALANCE_CONTROL_VOLFLUX  cgg(  Variables for splitting barotropic/baroclinic vels.
 cgg(  Variables for balancing volume flux.  
77        _RL ubaro        _RL ubaro
78        _RL ubarocount        _RL utop
79  cgg)  cgg)
 #endif  
80    
81  c     == external functions ==  c     == external functions ==
82    
# Line 97  c     == end of interface == Line 96  c     == end of interface ==
96        imax = snx+olx        imax = snx+olx
97        ip1  = 1        ip1  = 1
98    
 #ifdef BALANCE_CONTROL_VOLFLUX  
99  cgg(  Initialize variables for balancing volume flux.  cgg(  Initialize variables for balancing volume flux.
100        ubaro = 0.d0        ubaro = 0.d0
101        ubarocount = 0.d0        utop = 0.d0
102  cgg)  cgg)
 #endif  
103    
104  c--   Now, read the control vector.  c--   Now, read the control vector.
105        doglobalread = .false.        doglobalread = .false.
# Line 112  c--   Now, read the control vector. Line 109  c--   Now, read the control vector.
109          ilobcsw=ilnblnk( xx_obcsw_file )          ilobcsw=ilnblnk( xx_obcsw_file )
110          write(fnameobcsw(1:80),'(2a,i10.10)')          write(fnameobcsw(1:80),'(2a,i10.10)')
111       &       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.'  
112        endif        endif
113    
114  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
115        call ctrl_GetRec( 'xx_obcsw',        call ctrl_get_gen_rec(
116         I                   xx_obcswstartdate, xx_obcswperiod,
117       O                   obcswfac, obcswfirst, obcswchanged,       O                   obcswfac, obcswfirst, obcswchanged,
118       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
119       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
120    
121        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  
   
122          if ( obcswfirst ) then          if ( obcswfirst ) then
123            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
124       &                         (obcswcount0-1)*nobcs+iobcs,       &                         (obcswcount0-1)*nobcs+iobcs,
125       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
126       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
127    
128  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
129    
130            if ( optimcycle .gt. 0) then                      if ( optimcycle .gt. 0) then          
131              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
132  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.  
133  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
134    cgg         This is done on a columnwise basis here.
135                do bj = jtlo,jthi                do bj = jtlo,jthi
136                  do bi = itlo, ithi                  do bi = itlo, ithi
137                    do j = jmin,jmax                    do j = jmin,jmax
138                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
139                      ubaro = 0.d0  
140                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
141                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
142  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
143                        ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)                      utop = 0.d0
144       &                      * delZ(k) + ubaro  
145                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
146       &                           * delZ(k) +ubarocount  cgg    If cells are not full, this should be modified with hFac.
147                      enddo  cgg    
148                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
149                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
150                      else  cgg    exactly balance the volume flux, since we are dealing with
151                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
152                      endif                        utop = tmpfldyz(j,k,bi,bj)*
153                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
154                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
155                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
156                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
157                          endif
158                      enddo                      enddo
159    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
160                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
161         &                                      - utop / delR(1)
162                    enddo                    enddo
163                  enddo                  enddo
164                enddo                enddo
165              endif              endif
166            endif              if (iobcs .eq. 4) then
167  cgg)  cgg         Special attention is needed for the normal velocity.
168  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
169  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
170                do bj = jtlo,jthi                do bj = jtlo,jthi
171                  do bi = itlo, ithi                  do bi = itlo, ithi
172                    do k = 1,Nr                    do j = jmin,jmax
173                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
174                        i = OB_Iw(j,bi,bj)  
175                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
176       &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
177                        tmpfldyz(j,1,bi,bj) = 0.d0
178                        utop = 0.d0
179    
180                        do k = 1,Nr
181    cgg    If cells are not full, this should be modified with hFac.
182    cgg    
183    cgg    The xx field (tmpfldxz) does not contain the velocity at the
184    cgg    surface level. This velocity is not independent; it must
185    cgg    exactly balance the volume flux, since we are dealing with
186    cgg    the baroclinic velocity structure..
187                          utop = tmpfldyz(j,k,bi,bj)*
188         &                maskS(i,j,k,bi,bj) * delR(k) + utop
189    cgg    Add the barotropic velocity component.
190                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
191                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
192                          endif
193                      enddo                      enddo
194    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
195                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
196         &                                      - utop / delR(1)
197                    enddo                    enddo
198                  enddo                  enddo
199                enddo                enddo
200              endif              endif
201            endif            endif
202  #endif  
203    #endif /* ALLOW_CTRL_OBCS_BALANCE */
204    
205            do bj = jtlo,jthi            do bj = jtlo,jthi
206              do bi = itlo,ithi              do bi = itlo,ithi
207                do k = 1,nr                do k = 1,nr
208                  do j = jmin,jmax                  do j = jmin,jmax
209                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
210    cgg     &                                        *   maskyz (j,k,bi,bj)
211                   enddo                   enddo
212                enddo                enddo
213              enddo              enddo
# Line 204  cgg) Line 216  cgg)
216    
217          if ( (obcswfirst) .or. (obcswchanged)) then          if ( (obcswfirst) .or. (obcswchanged)) then
218                        
219  cgg(    This is a terribly long way to do it. However, the dimensions don't exactly  cgg(    This is a terribly long way to do it. However, the dimensions do not exactly
220  cgg     match up. I will blame Fortran for the ugliness.  cgg     match up. I will blame Fortran for the ugliness.
221    
222            do bj = jtlo,jthi            do bj = jtlo,jthi
223              do bi = itlo,ithi              do bi = itlo,ithi
224                do k = 1,nr                do k = 1,nr
225                  do j = jmin,jmax                  do j = jmin,jmax
   
226                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
227                  enddo                  enddo
228                enddo                enddo
# Line 224  cgg     match up. I will blame Fortran f Line 235  cgg     match up. I will blame Fortran f
235              do bi = itlo,ithi              do bi = itlo,ithi
236                do k = 1,nr                do k = 1,nr
237                  do j = jmin,jmax                  do j = jmin,jmax
238                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                            xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
239                  enddo                  enddo
240                enddo                enddo
241              enddo              enddo
# Line 235  cgg     match up. I will blame Fortran f Line 246  cgg     match up. I will blame Fortran f
246       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
247       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
248    
249  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
250            if (optimcycle .gt. 0) then  
251              if ( optimcycle .gt. 0) then          
252              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
253  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.  
254  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
255    cgg         This is done on a columnwise basis here.
256                do bj = jtlo,jthi                do bj = jtlo,jthi
257                  do bi = itlo, ithi                  do bi = itlo, ithi
258                    do j = jmin,jmax                    do j = jmin,jmax
259                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
260                      ubaro = 0.d0  
261                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
262                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
263  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
264                        ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)                      utop = 0.d0
265       &                      * delZ(k) + ubaro  
266                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
267       &                           * delZ(k) + ubarocount  cgg    If cells are not full, this should be modified with hFac.
268                      enddo  cgg    
269                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
270                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
271                      else  cgg    exactly balance the volume flux, since we are dealing with
272                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
273                      endif                        utop = tmpfldyz(j,k,bi,bj)*
274                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
275                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
276                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
277                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
278                          endif
279                      enddo                      enddo
280    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
281                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
282         &                                      - utop / delR(1)
283                    enddo                    enddo
284                  enddo                  enddo
285                enddo                enddo
286              endif              endif
287            endif              if (iobcs .eq. 4) then
288  cgg)  cgg         Special attention is needed for the normal velocity.
289  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
290  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
291                do bj = jtlo,jthi                do bj = jtlo,jthi
292                  do bi = itlo, ithi                  do bi = itlo, ithi
293                    do k = 1,Nr                    do j = jmin,jmax
294                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
295                        i = OB_Iw(j,bi,bj)  
296                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
297       &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
298                        tmpfldyz(j,1,bi,bj) = 0.d0
299                        utop = 0.d0
300    
301                        do k = 1,Nr
302    cgg    If cells are not full, this should be modified with hFac.
303    cgg    
304    cgg    The xx field (tmpfldxz) does not contain the velocity at the
305    cgg    surface level. This velocity is not independent; it must
306    cgg    exactly balance the volume flux, since we are dealing with
307    cgg    the baroclinic velocity structure..
308                          utop = tmpfldyz(j,k,bi,bj)*
309         &                maskS(i,j,k,bi,bj) * delR(k) + utop
310    cgg    Add the barotropic velocity component.
311                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
312                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
313                          endif
314                      enddo                      enddo
315    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
316                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
317         &                                      - utop / delR(1)
318                    enddo                    enddo
319                  enddo                  enddo
320                enddo                enddo
321              endif              endif
322            endif            endif
323  #endif  
324    #endif /* ALLOW_CTRL_OBCS_BALANCE */
325    
326            do bj = jtlo,jthi            do bj = jtlo,jthi
327              do bi = itlo,ithi              do bi = itlo,ithi
328                do k = 1,nr                do k = 1,nr
329                  do j = jmin,jmax                  do j = jmin,jmax
330                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
331    cgg     &                                        *   maskyz (j,k,bi,bj)
332                   enddo                   enddo
333                enddo                enddo
334              enddo              enddo
# Line 299  cgg) Line 336  cgg)
336          endif          endif
337    
338  c--     Add control to model variable.  c--     Add control to model variable.
339          do bj = jtlo,jthi          do bj = jtlo, jthi
340             do bi = itlo,ithi             do bi = itlo, ithi
341  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
342                do k = 1,nr                do k = 1,nr
343                   do j = 1,sny                   do j = 1,sny
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22