/[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.1 by heimbach, Tue Feb 5 20:23:58 2002 UTC revision 1.7 by heimbach, Mon May 14 22:02:33 2007 UTC
# Line 15  c     ================================== Line 15  c     ==================================
15  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
16  c     ==================================================================  c     ==================================================================
17  c  c
18  c     o Get norhtern obc of the control vector and add it  c     o Get western obc of the control vector and add it
19  c       to dyn. fields  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 56  c     == local variables == Line 57  c     == local variables ==
57        integer imin,imax        integer imin,imax
58        integer ilobcsw        integer ilobcsw
59        integer iobcs        integer iobcs
60          integer ip1
61    
62        _RL     dummy        _RL     dummy
63        _RL     obcswfac        _RL     obcswfac
# Line 64  c     == local variables == Line 66  c     == local variables ==
66        integer obcswcount0        integer obcswcount0
67        integer obcswcount1        integer obcswcount1
68    
69        _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
70    
71        logical doglobalread        logical doglobalread
72        logical ladinit        logical ladinit
73    
74        character*(80) fnameobcsw        character*(80) fnameobcsw
75    
76          cgg(  Variables for splitting barotropic/baroclinic vels.
77          _RL ubaro
78          _RL utop
79    cgg)
80    
81  c     == external functions ==  c     == external functions ==
82    
# Line 89  c     == end of interface == Line 94  c     == end of interface ==
94        jmax = sny+oly        jmax = sny+oly
95        imin = 1-olx        imin = 1-olx
96        imax = snx+olx        imax = snx+olx
97          ip1  = 1
98    
99    cgg(  Initialize variables for balancing volume flux.
100          ubaro = 0.d0
101          utop = 0.d0
102    cgg)
103    
104  c--   Now, read the control vector.  c--   Now, read the control vector.
105        doglobalread = .false.        doglobalread = .false.
# Line 98  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
122            if ( obcswfirst ) then
123              call active_read_yz( fnameobcsw, tmpfldyz,
124         &                         (obcswcount0-1)*nobcs+iobcs,
125         &                         doglobalread, ladinit, optimcycle,
126         &                         mythid, xx_obcsw_dummy )
127    
128    #ifdef ALLOW_CTRL_OBCS_BALANCE
129    
130              if ( optimcycle .gt. 0) then          
131                if (iobcs .eq. 3) then
132    cgg         Special attention is needed for the normal velocity.
133    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
136                    do bi = itlo, ithi
137                      do j = jmin,jmax
138                        i = OB_Iw(J,bi,bj)
139    
140    cgg         The barotropic velocity is stored in the level 1.
141                        ubaro = tmpfldyz(j,1,bi,bj)
142                        tmpfldyz(j,1,bi,bj) = 0.d0
143                        utop = 0.d0
144    
145                        do k = 1,Nr
146    cgg    If cells are not full, this should be modified with hFac.
147    cgg    
148    cgg    The xx field (tmpfldxz) does not contain the velocity at the
149    cgg    surface level. This velocity is not independent; it must
150    cgg    exactly balance the volume flux, since we are dealing with
151    cgg    the baroclinic velocity structure..
152                          utop = tmpfldyz(j,k,bi,bj)*
153         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
154    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
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
163                    enddo
164                  enddo
165                endif
166                if (iobcs .eq. 4) then
167    cgg         Special attention is needed for the normal velocity.
168    cgg         For the north, this is the v velocity, iobcs = 4.
169    cgg         This is done on a columnwise basis here.
170                  do bj = jtlo,jthi
171                    do bi = itlo, ithi
172                      do j = jmin,jmax
173                        i = OB_Iw(J,bi,bj)
174    
175    cgg         The barotropic velocity is stored in the level 1.
176                        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
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
198                    enddo
199                  enddo
200                endif
201              endif
202    
203    #endif /* ALLOW_CTRL_OBCS_BALANCE */
204    
205              do bj = jtlo,jthi
206                do bi = itlo,ithi
207                  do k = 1,nr
208                    do j = jmin,jmax
209                      xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
210    cgg     &                                        *   maskyz (j,k,bi,bj)
211                     enddo
212                  enddo
213                enddo
214              enddo
215            endif
216    
217            if ( (obcswfirst) .or. (obcswchanged)) then
218              
219    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.
221    
222              do bj = jtlo,jthi
223                do bi = itlo,ithi
224                  do k = 1,nr
225                    do j = jmin,jmax
226                      tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
227                    enddo
228                  enddo
229                enddo
230              enddo
231    
232              call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
233    
234          call active_read_yz( 'maskobcsw', maskyz,            do bj = jtlo,jthi
235       &                       iobcs,              do bi = itlo,ithi
236       &                       doglobalread, ladinit, 0,                do k = 1,nr
237       &                       mythid, dummy )                  do j = jmin,jmax
238                      xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
239          call active_read_yz( fnameobcsw, tmpfldyz,                  enddo
      &                       (obcswcount0-1)*nobcs+iobcs,  
      &                       doglobalread, ladinit, optimcycle,  
      &                       mythid, xx_obcsw_dummy )  
   
         do bj = jtlo,jthi  
           do bi = itlo,ithi  
             do k = 1,nr  
               do j = jmin,jmax  
                 yz_obcs0(j,k,bi,bj)  = tmpfldyz (j,k,bi,bj)  
240                enddo                enddo
241              enddo              enddo
242            enddo            enddo
         enddo  
243    
244          call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
245       &                       (obcswcount1-1)*nobcs+iobcs,       &                         (obcswcount1-1)*nobcs+iobcs,
246       &                       doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
247       &                       mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
248    
249          do bj = jtlo,jthi  #ifdef ALLOW_CTRL_OBCS_BALANCE
250            do bi = itlo,ithi  
251              do k = 1,nr            if ( optimcycle .gt. 0) then          
252                do j = jmin,jmax              if (iobcs .eq. 3) then
253                  yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj)  cgg         Special attention is needed for the normal velocity.
254    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
257                    do bi = itlo, ithi
258                      do j = jmin,jmax
259                        i = OB_Iw(J,bi,bj)
260    
261    cgg         The barotropic velocity is stored in the level 1.
262                        ubaro = tmpfldyz(j,1,bi,bj)
263                        tmpfldyz(j,1,bi,bj) = 0.d0
264                        utop = 0.d0
265    
266                        do k = 1,Nr
267    cgg    If cells are not full, this should be modified with hFac.
268    cgg    
269    cgg    The xx field (tmpfldxz) does not contain the velocity at the
270    cgg    surface level. This velocity is not independent; it must
271    cgg    exactly balance the volume flux, since we are dealing with
272    cgg    the baroclinic velocity structure..
273                          utop = tmpfldyz(j,k,bi,bj)*
274         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
275    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
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
284                    enddo
285                  enddo
286                endif
287                if (iobcs .eq. 4) then
288    cgg         Special attention is needed for the normal velocity.
289    cgg         For the north, this is the v velocity, iobcs = 4.
290    cgg         This is done on a columnwise basis here.
291                  do bj = jtlo,jthi
292                    do bi = itlo, ithi
293                      do j = jmin,jmax
294                        i = OB_Iw(J,bi,bj)
295    
296    cgg         The barotropic velocity is stored in the level 1.
297                        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
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
319                    enddo
320                  enddo
321                endif
322              endif
323    
324    #endif /* ALLOW_CTRL_OBCS_BALANCE */
325    
326              do bj = jtlo,jthi
327                do bi = itlo,ithi
328                  do k = 1,nr
329                    do j = jmin,jmax
330                      xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
331    cgg     &                                        *   maskyz (j,k,bi,bj)
332                     enddo
333                enddo                enddo
334              enddo              enddo
335            enddo            enddo
336          enddo          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
344                        i = OB_Iw(j,bi,bj)
345                      if (iobcs .EQ. 1) then                      if (iobcs .EQ. 1) then
346                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
347       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
348       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
349                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
350       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
351                      else if (iobcs .EQ. 2) then                      else if (iobcs .EQ. 2) then
352                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
353       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
354       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
355                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
356       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
357                      else if (iobcs .EQ. 3) then                      else if (iobcs .EQ. 3) then
358                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
359       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
360       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
361                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
362       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
363                      else if (iobcs .EQ. 4) then                      else if (iobcs .EQ. 4) then
364                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
365       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
366       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
367                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
368       &                      *maskyz(j,k,bi,bj)       &                      *maskS(i,j,k,bi,bj)
369                      endif                      endif
370                   enddo                   enddo
371                enddo                enddo
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22