/[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.10 by mlosch, Mon Mar 7 09:24:10 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 66  c     == local variables == Line 69  c     == local variables ==
69        integer obcswcount1        integer obcswcount1
70    
71  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
72          _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy)
73    
74        logical doglobalread        logical doglobalread
75        logical ladinit        logical ladinit
76    
77        character*(80) fnameobcsw        character*(80) fnameobcsw
78    
79  #ifdef BALANCE_CONTROL_VOLFLUX  cgg(  Variables for splitting barotropic/baroclinic vels.
 cgg(  Variables for balancing volume flux.  
80        _RL ubaro        _RL ubaro
81        _RL ubarocount        _RL utop
82  cgg)  cgg)
 #endif  
83    
84  c     == external functions ==  c     == external functions ==
85    
# Line 97  c     == end of interface == Line 99  c     == end of interface ==
99        imax = snx+olx        imax = snx+olx
100        ip1  = 1        ip1  = 1
101    
 #ifdef BALANCE_CONTROL_VOLFLUX  
102  cgg(  Initialize variables for balancing volume flux.  cgg(  Initialize variables for balancing volume flux.
103        ubaro = 0.d0        ubaro = 0.d0
104        ubarocount = 0.d0        utop = 0.d0
105  cgg)  cgg)
 #endif  
106    
107  c--   Now, read the control vector.  c--   Now, read the control vector.
108        doglobalread = .false.        doglobalread = .false.
# Line 110  c--   Now, read the control vector. Line 110  c--   Now, read the control vector.
110    
111        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
112          ilobcsw=ilnblnk( xx_obcsw_file )          ilobcsw=ilnblnk( xx_obcsw_file )
113          write(fnameobcsw(1:80),'(2a,i10.10)')          write(fnameobcsw(1:80),'(2a,i10.10)')
114       &       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.'  
115        endif        endif
116    
117  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
118        call ctrl_GetRec( 'xx_obcsw',        call ctrl_get_gen_rec(
119         I                   xx_obcswstartdate, xx_obcswperiod,
120       O                   obcswfac, obcswfirst, obcswchanged,       O                   obcswfac, obcswfirst, obcswchanged,
121       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
122       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
123    
124    CML      print *,'ml-getobcs: ',myIter,obcswfirst,obcswchanged,
125    CML     &     obcswcount0,obcswcount1,obcswfac
126        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  
   
127          if ( obcswfirst ) then          if ( obcswfirst ) then
128            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
129       &                         (obcswcount0-1)*nobcs+iobcs,       &                         (obcswcount0-1)*nobcs+iobcs,
130       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
131       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
132    
133  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
134            if ( optimcycle .gt. 0) then            
135              if ( optimcycle .gt. 0) then
136              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
137  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.  
138  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
139    cgg         This is done on a columnwise basis here.
140                do bj = jtlo,jthi                do bj = jtlo,jthi
141                  do bi = itlo, ithi                  do bi = itlo, ithi
142                    do j = jmin,jmax                    do j = jmin,jmax
143                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
144                      ubaro = 0.d0  
145                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
146                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
147  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
148                        ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)                      utop = 0.d0
149       &                      * delZ(k) + ubaro  
150                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
151       &                           * delZ(k) +ubarocount  cgg    If cells are not full, this should be modified with hFac.
152                      enddo  cgg
153                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldyz) does not contain the velocity at the
154                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
155                      else  cgg    exactly balance the volume flux, since we are dealing with
156                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
157                      endif                        utop = tmpfldyz(j,k,bi,bj)*
158                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
159                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
160                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
161                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
162                          endif
163                      enddo                      enddo
164    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
165                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
166         &                                      - utop / delR(1)
167                    enddo                    enddo
168                  enddo                  enddo
169                enddo                enddo
170              endif              endif
171            endif              if (iobcs .eq. 4) then
172  cgg)  cgg         Special attention is needed for the normal velocity.
173  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
174  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
175                do bj = jtlo,jthi                do bj = jtlo,jthi
176                  do bi = itlo, ithi                  do bi = itlo, ithi
177                    do k = 1,Nr                    do j = jmin,jmax
178                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
179                        i = OB_Iw(j,bi,bj)  
180                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
181       &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
182                        tmpfldyz(j,1,bi,bj) = 0.d0
183                        utop = 0.d0
184    
185                        do k = 1,Nr
186    cgg    If cells are not full, this should be modified with hFac.
187    cgg
188    cgg    The xx field (tmpfldyz) does not contain the velocity at the
189    cgg    surface level. This velocity is not independent; it must
190    cgg    exactly balance the volume flux, since we are dealing with
191    cgg    the baroclinic velocity structure..
192                          utop = tmpfldyz(j,k,bi,bj)*
193         &                maskS(i,j,k,bi,bj) * delR(k) + utop
194    cgg    Add the barotropic velocity component.
195                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
196                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
197                          endif
198                      enddo                      enddo
199    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
200                      tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
201         &                                      - utop / delR(1)
202                    enddo                    enddo
203                  enddo                  enddo
204                enddo                enddo
205              endif              endif
206            endif            endif
207  #endif  
208    #endif /* ALLOW_CTRL_OBCS_BALANCE */
209    
210            do bj = jtlo,jthi            do bj = jtlo,jthi
211              do bi = itlo,ithi              do bi = itlo,ithi
212                do k = 1,nr                do k = 1,nr
213                  do j = jmin,jmax                  do j = jmin,jmax
214                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
215    cgg     &                                        *   maskyz (j,k,bi,bj)
216                   enddo                   enddo
217                enddo                enddo
218              enddo              enddo
# Line 203  cgg) Line 220  cgg)
220          endif          endif
221    
222          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.  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
   
                   tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)  
                 enddo  
               enddo  
             enddo  
           enddo  
   
           call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)  
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                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                        xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
229                  enddo                tmpfldyz (j,k,bi,bj)       = 0. _d 0
230                enddo               enddo
231              enddo              enddo
232               enddo
233            enddo            enddo
234    
235            call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
236       &                         (obcswcount1-1)*nobcs+iobcs,       &                         (obcswcount1-1)*nobcs+iobcs,
237       &                         doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
238       &                         mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
239    
240  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_CTRL_OBCS_BALANCE
241            if (optimcycle .gt. 0) then  
242              if ( optimcycle .gt. 0) then
243              if (iobcs .eq. 3) then              if (iobcs .eq. 3) then
244  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.  
245  cgg         For the north, this is the v velocity, iobcs = 4.  cgg         For the north, this is the v velocity, iobcs = 4.
246    cgg         This is done on a columnwise basis here.
247                do bj = jtlo,jthi                do bj = jtlo,jthi
248                  do bi = itlo, ithi                  do bi = itlo, ithi
249                    do j = jmin,jmax                    do j = jmin,jmax
250                      i = OB_Iw(j,bi,bj)                      i = OB_Iw(J,bi,bj)
251                      ubaro = 0.d0  
252                      ubarocount = 0.d0  cgg         The barotropic velocity is stored in the level 1.
253                      do k = 1,nr                      ubaro = tmpfldyz(j,1,bi,bj)
254  cgg(   If cells are not full, this should be modified with hFac.                      tmpfldyz(j,1,bi,bj) = 0.d0
255                        ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)                      utop = 0.d0
256       &                      * delZ(k) + ubaro  
257                        ubarocount = maskW(i+ip1,j,k,bi,bj)                      do k = 1,Nr
258       &                           * delZ(k) + ubarocount  cgg    If cells are not full, this should be modified with hFac.
259                      enddo  cgg
260                      if (ubarocount .eq. 0.) then  cgg    The xx field (tmpfldyz) does not contain the velocity at the
261                        ubaro = 0.  cgg    surface level. This velocity is not independent; it must
262                      else  cgg    exactly balance the volume flux, since we are dealing with
263                        ubaro = ubaro / ubarocount  cgg    the baroclinic velocity structure..
264                      endif                        utop = tmpfldyz(j,k,bi,bj)*
265                      do k = 1,nr       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
266                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  cgg    Add the barotropic velocity component.
267                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
268                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
269                          endif
270                      enddo                      enddo
271    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
272                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
273         &                                      - utop / delR(1)
274                    enddo                    enddo
275                  enddo                  enddo
276                enddo                enddo
277              endif              endif
278            endif              if (iobcs .eq. 4) then
279  cgg)  cgg         Special attention is needed for the normal velocity.
280  #endif  cgg         For the north, this is the v velocity, iobcs = 4.
281  #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  cgg         This is done on a columnwise basis here.
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
282                do bj = jtlo,jthi                do bj = jtlo,jthi
283                  do bi = itlo, ithi                  do bi = itlo, ithi
284                    do k = 1,Nr                    do j = jmin,jmax
285                      do j = jmin,jmax                      i = OB_Iw(J,bi,bj)
286                        i = OB_Iw(j,bi,bj)  
287                        tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  cgg         The barotropic velocity is stored in the level 1.
288       &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)                      ubaro = tmpfldyz(j,1,bi,bj)
289                        tmpfldyz(j,1,bi,bj) = 0.d0
290                        utop = 0.d0
291    
292                        do k = 1,Nr
293    cgg    If cells are not full, this should be modified with hFac.
294    cgg
295    cgg    The xx field (tmpfldyz) does not contain the velocity at the
296    cgg    surface level. This velocity is not independent; it must
297    cgg    exactly balance the volume flux, since we are dealing with
298    cgg    the baroclinic velocity structure..
299                          utop = tmpfldyz(j,k,bi,bj)*
300         &                maskS(i,j,k,bi,bj) * delR(k) + utop
301    cgg    Add the barotropic velocity component.
302                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
303                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
304                          endif
305                      enddo                      enddo
306    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
307                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
308         &                                      - utop / delR(1)
309                    enddo                    enddo
310                  enddo                  enddo
311                enddo                enddo
312              endif              endif
313            endif            endif
314  #endif  
315    #endif /* ALLOW_CTRL_OBCS_BALANCE */
316    
317            do bj = jtlo,jthi            do bj = jtlo,jthi
318              do bi = itlo,ithi              do bi = itlo,ithi
319                do k = 1,nr                do k = 1,nr
320                  do j = jmin,jmax                  do j = jmin,jmax
321                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
322    cgg     &                                        *   maskyz (j,k,bi,bj)
323                   enddo                   enddo
324                enddo                enddo
325              enddo              enddo
326            enddo            enddo
327          endif          endif
328    
329  c--     Add control to model variable.  c--   Add control to model variable.
330          do bj = jtlo,jthi          do bj = jtlo, jthi
331             do bi = itlo,ithi           do bi = itlo, ithi
332  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--   Calculate mask for tracer cells (0 => land, 1 => water).
333                do k = 1,nr            do k = 1,nr
334                   do j = 1,sny             do j = 1,sny
335                      i = OB_Iw(j,bi,bj)              i = OB_Iw(j,bi,bj)
336                      if (iobcs .EQ. 1) then              if (iobcs .EQ. 1) then
337                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)               OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
338       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
339       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
340                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)               OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
341       &                      *maskW(i+ip1,j,k,bi,bj)       &            *maskW(i+ip1,j,k,bi,bj)
342                      else if (iobcs .EQ. 2) then              else if (iobcs .EQ. 2) then
343                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)               OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
344       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
345       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
346                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)               OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
347       &                      *maskW(i+ip1,j,k,bi,bj)       &            *maskW(i+ip1,j,k,bi,bj)
348                      else if (iobcs .EQ. 3) then              else if (iobcs .EQ. 3) then
349                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)               OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
350       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
351       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
352                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)               OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
353       &                      *maskW(i+ip1,j,k,bi,bj)       &            *maskW(i+ip1,j,k,bi,bj)
354                      else if (iobcs .EQ. 4) then              else if (iobcs .EQ. 4) then
355                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)               OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
356       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
357       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
358                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)               OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
359       &                      *maskS(i,j,k,bi,bj)       &            *maskS(i,j,k,bi,bj)
360                      endif              endif
                  enddo  
               enddo  
361             enddo             enddo
362              enddo
363             enddo
364          enddo          enddo
365    
366  C--   End over iobcs loop  C--   End over iobcs loop
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22