/[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.8 by jmc, Tue Oct 9 00:00:00 2007 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
220              
221  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
222  cgg     match up. I will blame Fortran for the ugliness.  cgg     match up. I will blame Fortran for the ugliness.
223    
224            do bj = jtlo,jthi            do bj = jtlo,jthi
225              do bi = itlo,ithi              do bi = itlo,ithi
226                do k = 1,nr                do k = 1,nr
227                  do j = jmin,jmax                  do j = jmin,jmax
   
228                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
229                  enddo                  enddo
230                enddo                enddo
# Line 224  cgg     match up. I will blame Fortran f Line 237  cgg     match up. I will blame Fortran f
237              do bi = itlo,ithi              do bi = itlo,ithi
238                do k = 1,nr                do k = 1,nr
239                  do j = jmin,jmax                  do j = jmin,jmax
240                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                            xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
241                  enddo                  enddo
242                enddo                enddo
243              enddo              enddo
244            enddo            enddo
245    
246            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
247       &                         (obcswcount1-1)*nobcs+iobcs,       &                         (obcswcount1-1)*nobcs+iobcs,
248       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
249       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
250    
251  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
252            if (optimcycle .gt. 0) then  
253              if ( optimcycle .gt. 0) then
254              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
255  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.  
256  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
257    cgg         This is done on a columnwise basis here.
258                do bj = jtlo,jthi                do bj = jtlo,jthi
259                  do bi = itlo, ithi                  do bi = itlo, ithi
260                    do j = jmin,jmax                    do j = jmin,jmax
261                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
262                      ubaro = 0.d0  
263                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
264                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
265  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
266                        ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)                      utop = 0.d0
267       &                      * delZ(k) + ubaro  
268                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
269       &                           * delZ(k) + ubarocount  cgg    If cells are not full, this should be modified with hFac.
270                      enddo  cgg
271                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
272                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
273                      else  cgg    exactly balance the volume flux, since we are dealing with
274                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
275                      endif                        utop = tmpfldyz(j,k,bi,bj)*
276                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
277                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
278                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
279                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
280                          endif
281                      enddo                      enddo
282    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
283                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
284         &                                      - utop / delR(1)
285                    enddo                    enddo
286                  enddo                  enddo
287                enddo                enddo
288              endif              endif
289            endif              if (iobcs .eq. 4) then
290  cgg)  cgg         Special attention is needed for the normal velocity.
291  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
292  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
293                do bj = jtlo,jthi                do bj = jtlo,jthi
294                  do bi = itlo, ithi                  do bi = itlo, ithi
295                    do k = 1,Nr                    do j = jmin,jmax
296                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
297                        i = OB_Iw(j,bi,bj)  
298                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
299       &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
300                        tmpfldyz(j,1,bi,bj) = 0.d0
301                        utop = 0.d0
302    
303                        do k = 1,Nr
304    cgg    If cells are not full, this should be modified with hFac.
305    cgg
306    cgg    The xx field (tmpfldxz) does not contain the velocity at the
307    cgg    surface level. This velocity is not independent; it must
308    cgg    exactly balance the volume flux, since we are dealing with
309    cgg    the baroclinic velocity structure..
310                          utop = tmpfldyz(j,k,bi,bj)*
311         &                maskS(i,j,k,bi,bj) * delR(k) + utop
312    cgg    Add the barotropic velocity component.
313                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
314                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
315                          endif
316                      enddo                      enddo
317    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
318                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
319         &                                      - utop / delR(1)
320                    enddo                    enddo
321                  enddo                  enddo
322                enddo                enddo
323              endif              endif
324            endif            endif
325  #endif  
326    #endif /* ALLOW_CTRL_OBCS_BALANCE */
327    
328            do bj = jtlo,jthi            do bj = jtlo,jthi
329              do bi = itlo,ithi              do bi = itlo,ithi
330                do k = 1,nr                do k = 1,nr
331                  do j = jmin,jmax                  do j = jmin,jmax
332                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
333    cgg     &                                        *   maskyz (j,k,bi,bj)
334                   enddo                   enddo
335                enddo                enddo
336              enddo              enddo
# Line 299  cgg) Line 338  cgg)
338          endif          endif
339    
340  c--     Add control to model variable.  c--     Add control to model variable.
341          do bj = jtlo,jthi          do bj = jtlo, jthi
342             do bi = itlo,ithi             do bi = itlo, ithi
343  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
344                do k = 1,nr                do k = 1,nr
345                   do j = 1,sny                   do j = 1,sny
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22