/[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.1.2.3 by heimbach, Thu Jun 19 15:18:48 2003 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    
 #ifdef BALANCE_CONTROL_VOLFLUX  
128            if ( optimcycle .gt. 0) then                      if ( optimcycle .gt. 0) then          
129              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
130  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.  
131  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
132    cgg         This is done on a columnwise basis here.
133                do bj = jtlo,jthi                do bj = jtlo,jthi
134                  do bi = itlo, ithi                  do bi = itlo, ithi
135                    do j = jmin,jmax                    do j = jmin,jmax
136                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
137                      ubaro = 0.d0  
138                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
139                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
140  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
141                        ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)                      utop = 0.d0
142       &                      * delZ(k) + ubaro  
143                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
144       &                           * delZ(k) +ubarocount  cgg    If cells are not full, this should be modified with hFac.
145                      enddo  cgg    
146                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
147                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
148                      else  cgg    exactly balance the volume flux, since we are dealing with
149                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
150                      endif                        utop = tmpfldyz(j,k,bi,bj)*
151                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop
152                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
153                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
154                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
155                          endif
156                      enddo                      enddo
157    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
158                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
159         &                                      - utop / delZ(1)
160                    enddo                    enddo
161                  enddo                  enddo
162                enddo                enddo
163              endif              endif
164            endif              if (iobcs .eq. 4) then
165  cgg)  cgg         Special attention is needed for the normal velocity.
166  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
167  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
168                do bj = jtlo,jthi                do bj = jtlo,jthi
169                  do bi = itlo, ithi                  do bi = itlo, ithi
170                    do k = 1,Nr                    do j = jmin,jmax
171                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
172                        i = OB_Iw(j,bi,bj)  
173                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
174       &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
175                        tmpfldyz(j,1,bi,bj) = 0.d0
176                        utop = 0.d0
177    
178                        do k = 1,Nr
179    cgg    If cells are not full, this should be modified with hFac.
180    cgg    
181    cgg    The xx field (tmpfldxz) does not contain the velocity at the
182    cgg    surface level. This velocity is not independent; it must
183    cgg    exactly balance the volume flux, since we are dealing with
184    cgg    the baroclinic velocity structure..
185                          utop = tmpfldyz(j,k,bi,bj)*
186         &                maskS(i,j,k,bi,bj) * delZ(k) + utop
187    cgg    Add the barotropic velocity component.
188                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
189                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
190                          endif
191                      enddo                      enddo
192    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
193                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
194         &                                      - utop / delZ(1)
195                    enddo                    enddo
196                  enddo                  enddo
197                enddo                enddo
198              endif              endif
199            endif            endif
 #endif  
200    
201            do bj = jtlo,jthi            do bj = jtlo,jthi
202              do bi = itlo,ithi              do bi = itlo,ithi
203                do k = 1,nr                do k = 1,nr
204                  do j = jmin,jmax                  do j = jmin,jmax
205                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
206    cgg     &                                        *   maskyz (j,k,bi,bj)
207                   enddo                   enddo
208                enddo                enddo
209              enddo              enddo
# Line 211  cgg     match up. I will blame Fortran f Line 219  cgg     match up. I will blame Fortran f
219              do bi = itlo,ithi              do bi = itlo,ithi
220                do k = 1,nr                do k = 1,nr
221                  do j = jmin,jmax                  do j = jmin,jmax
   
222                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
223                  enddo                  enddo
224                enddo                enddo
# Line 224  cgg     match up. I will blame Fortran f Line 231  cgg     match up. I will blame Fortran f
231              do bi = itlo,ithi              do bi = itlo,ithi
232                do k = 1,nr                do k = 1,nr
233                  do j = jmin,jmax                  do j = jmin,jmax
234                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                            xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
235                  enddo                  enddo
236                enddo                enddo
237              enddo              enddo
# Line 235  cgg     match up. I will blame Fortran f Line 242  cgg     match up. I will blame Fortran f
242       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
243       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
244    
245  #ifdef BALANCE_CONTROL_VOLFLUX            if ( optimcycle .gt. 0) then          
           if (optimcycle .gt. 0) then  
246              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
247  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.  
248  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
249    cgg         This is done on a columnwise basis here.
250                do bj = jtlo,jthi                do bj = jtlo,jthi
251                  do bi = itlo, ithi                  do bi = itlo, ithi
252                    do j = jmin,jmax                    do j = jmin,jmax
253                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
254                      ubaro = 0.d0  
255                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
256                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
257  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
258                        ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)                      utop = 0.d0
259       &                      * delZ(k) + ubaro  
260                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
261       &                           * delZ(k) + ubarocount  cgg    If cells are not full, this should be modified with hFac.
262                      enddo  cgg    
263                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldxz) does not contain the velocity at the
264                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
265                      else  cgg    exactly balance the volume flux, since we are dealing with
266                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
267                      endif                        utop = tmpfldyz(j,k,bi,bj)*
268                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop
269                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
270                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
271                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
272                          endif
273                      enddo                      enddo
274    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
275                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
276         &                                      - utop / delZ(1)
277                    enddo                    enddo
278                  enddo                  enddo
279                enddo                enddo
280              endif              endif
281            endif              if (iobcs .eq. 4) then
282  cgg)  cgg         Special attention is needed for the normal velocity.
283  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
284  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
285                do bj = jtlo,jthi                do bj = jtlo,jthi
286                  do bi = itlo, ithi                  do bi = itlo, ithi
287                    do k = 1,Nr                    do j = jmin,jmax
288                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
289                        i = OB_Iw(j,bi,bj)  
290                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
291       &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
292                        tmpfldyz(j,1,bi,bj) = 0.d0
293                        utop = 0.d0
294    
295                        do k = 1,Nr
296    cgg    If cells are not full, this should be modified with hFac.
297    cgg    
298    cgg    The xx field (tmpfldxz) does not contain the velocity at the
299    cgg    surface level. This velocity is not independent; it must
300    cgg    exactly balance the volume flux, since we are dealing with
301    cgg    the baroclinic velocity structure..
302                          utop = tmpfldyz(j,k,bi,bj)*
303         &                maskS(i,j,k,bi,bj) * delZ(k) + utop
304    cgg    Add the barotropic velocity component.
305                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
306                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
307                          endif
308                      enddo                      enddo
309    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
310                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
311         &                                      - utop / delZ(1)
312                    enddo                    enddo
313                  enddo                  enddo
314                enddo                enddo
315              endif              endif
316            endif            endif
 #endif  
317    
318            do bj = jtlo,jthi            do bj = jtlo,jthi
319              do bi = itlo,ithi              do bi = itlo,ithi
320                do k = 1,nr                do k = 1,nr
321                  do j = jmin,jmax                  do j = jmin,jmax
322                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
323    cgg     &                                        *   maskyz (j,k,bi,bj)
324                   enddo                   enddo
325                enddo                enddo
326              enddo              enddo
# Line 299  cgg) Line 328  cgg)
328          endif          endif
329    
330  c--     Add control to model variable.  c--     Add control to model variable.
331          do bj = jtlo,jthi          do bj = jtlo, jthi
332             do bi = itlo,ithi             do bi = itlo, ithi
333  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
334                do k = 1,nr                do k = 1,nr
335                   do j = 1,sny                   do j = 1,sny
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.1.2.3
  ViewVC Help
Powered by ViewVC 1.1.22