/[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.9 by mlosch, Wed Jan 19 08:42:06 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 15  c     ================================== Line 17  c     ==================================
17  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
18  c     ==================================================================  c     ==================================================================
19  c  c
20  c     o Get norhtern obc of the control vector and add it  c     o Get western obc of the control vector and add it
21  c       to dyn. fields  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 56  c     == local variables == Line 59  c     == local variables ==
59        integer imin,imax        integer imin,imax
60        integer ilobcsw        integer ilobcsw
61        integer iobcs        integer iobcs
62          integer ip1
63    
64        _RL     dummy        _RL     dummy
65        _RL     obcswfac        _RL     obcswfac
# Line 64  c     == local variables == Line 68  c     == local variables ==
68        integer obcswcount0        integer obcswcount0
69        integer obcswcount1        integer obcswcount1
70    
71        _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
72    
73        logical doglobalread        logical doglobalread
74        logical ladinit        logical ladinit
75    
76        character*(80) fnameobcsw        character*(80) fnameobcsw
77    
78          cgg(  Variables for splitting barotropic/baroclinic vels.
79          _RL ubaro
80          _RL utop
81    cgg)
82    
83  c     == external functions ==  c     == external functions ==
84    
# Line 89  c     == end of interface == Line 96  c     == end of interface ==
96        jmax = sny+oly        jmax = sny+oly
97        imin = 1-olx        imin = 1-olx
98        imax = snx+olx        imax = snx+olx
99          ip1  = 1
100    
101    cgg(  Initialize variables for balancing volume flux.
102          ubaro = 0.d0
103          utop = 0.d0
104    cgg)
105    
106  c--   Now, read the control vector.  c--   Now, read the control vector.
107        doglobalread = .false.        doglobalread = .false.
# Line 96  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
124            if ( obcswfirst ) then
125              call active_read_yz( fnameobcsw, tmpfldyz,
126         &                         (obcswcount0-1)*nobcs+iobcs,
127         &                         doglobalread, ladinit, optimcycle,
128         &                         mythid, xx_obcsw_dummy )
129    
130    #ifdef ALLOW_CTRL_OBCS_BALANCE
131    
132              if ( optimcycle .gt. 0) then
133                if (iobcs .eq. 3) then
134    cgg         Special attention is needed for the normal velocity.
135    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
138                    do bi = itlo, ithi
139                      do j = jmin,jmax
140                        i = OB_Iw(J,bi,bj)
141    
142    cgg         The barotropic velocity is stored in the level 1.
143                        ubaro = tmpfldyz(j,1,bi,bj)
144                        tmpfldyz(j,1,bi,bj) = 0.d0
145                        utop = 0.d0
146    
147                        do k = 1,Nr
148    cgg    If cells are not full, this should be modified with hFac.
149    cgg
150    cgg    The xx field (tmpfldxz) does not contain the velocity at the
151    cgg    surface level. This velocity is not independent; it must
152    cgg    exactly balance the volume flux, since we are dealing with
153    cgg    the baroclinic velocity structure..
154                          utop = tmpfldyz(j,k,bi,bj)*
155         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
156    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
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
165                    enddo
166                  enddo
167                endif
168                if (iobcs .eq. 4) then
169    cgg         Special attention is needed for the normal velocity.
170    cgg         For the north, this is the v velocity, iobcs = 4.
171    cgg         This is done on a columnwise basis here.
172                  do bj = jtlo,jthi
173                    do bi = itlo, ithi
174                      do j = jmin,jmax
175                        i = OB_Iw(J,bi,bj)
176    
177    cgg         The barotropic velocity is stored in the level 1.
178                        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
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
200                    enddo
201                  enddo
202                endif
203              endif
204    
205          call active_read_yz( 'maskobcsw', maskyz,  #endif /* ALLOW_CTRL_OBCS_BALANCE */
      &                       iobcs,  
      &                       doglobalread, ladinit, 0,  
      &                       mythid, dummy )  
   
         call active_read_yz( fnameobcsw, tmpfldyz,  
      &                       (obcswcount0-1)*nobcs+iobcs,  
      &                       doglobalread, ladinit, optimcycle,  
      &                       mythid, xx_obcsw_dummy )  
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                  yz_obcs0(j,k,bi,bj)  = 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
214                enddo                enddo
215              enddo              enddo
216            enddo            enddo
217          enddo          endif
218    
219          call active_read_yz( fnameobcsw, tmpfldyz,          if ( (obcswfirst) .or. (obcswchanged)) then
      &                       (obcswcount1-1)*nobcs+iobcs,  
      &                       doglobalread, ladinit, optimcycle,  
      &                       mythid, xx_obcsw_dummy )  
220    
221          do bj = jtlo,jthi            do bj = jtlo,jthi
222            do bi = itlo,ithi             do bi = itlo,ithi
223              do k = 1,nr              do k = 1,nr
224                do j = jmin,jmax               do j = jmin,jmax
225                  yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj)                xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
226                  tmpfldyz (j,k,bi,bj)       = 0. _d 0
227                 enddo
228                enddo
229               enddo
230              enddo
231    
232              call active_read_yz( fnameobcsw, tmpfldyz,
233         &                         (obcswcount1-1)*nobcs+iobcs,
234         &                         doglobalread, ladinit, optimcycle,
235         &                         mythid, xx_obcsw_dummy )
236    
237    #ifdef ALLOW_CTRL_OBCS_BALANCE
238    
239              if ( optimcycle .gt. 0) then
240                if (iobcs .eq. 3) then
241    cgg         Special attention is needed for the normal velocity.
242    cgg         For the north, this is the v velocity, iobcs = 4.
243    cgg         This is done on a columnwise basis here.
244                  do bj = jtlo,jthi
245                    do bi = itlo, ithi
246                      do j = jmin,jmax
247                        i = OB_Iw(J,bi,bj)
248    
249    cgg         The barotropic velocity is stored in the level 1.
250                        ubaro = tmpfldyz(j,1,bi,bj)
251                        tmpfldyz(j,1,bi,bj) = 0.d0
252                        utop = 0.d0
253    
254                        do k = 1,Nr
255    cgg    If cells are not full, this should be modified with hFac.
256    cgg
257    cgg    The xx field (tmpfldxz) does not contain the velocity at the
258    cgg    surface level. This velocity is not independent; it must
259    cgg    exactly balance the volume flux, since we are dealing with
260    cgg    the baroclinic velocity structure..
261                          utop = tmpfldyz(j,k,bi,bj)*
262         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
263    cgg    Add the barotropic velocity component.
264                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
265                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
266                          endif
267                        enddo
268    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
269                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
270         &                                      - utop / delR(1)
271                      enddo
272                    enddo
273                  enddo
274                endif
275                if (iobcs .eq. 4) then
276    cgg         Special attention is needed for the normal velocity.
277    cgg         For the north, this is the v velocity, iobcs = 4.
278    cgg         This is done on a columnwise basis here.
279                  do bj = jtlo,jthi
280                    do bi = itlo, ithi
281                      do j = jmin,jmax
282                        i = OB_Iw(J,bi,bj)
283    
284    cgg         The barotropic velocity is stored in the level 1.
285                        ubaro = tmpfldyz(j,1,bi,bj)
286                        tmpfldyz(j,1,bi,bj) = 0.d0
287                        utop = 0.d0
288    
289                        do k = 1,Nr
290    cgg    If cells are not full, this should be modified with hFac.
291    cgg
292    cgg    The xx field (tmpfldxz) does not contain the velocity at the
293    cgg    surface level. This velocity is not independent; it must
294    cgg    exactly balance the volume flux, since we are dealing with
295    cgg    the baroclinic velocity structure..
296                          utop = tmpfldyz(j,k,bi,bj)*
297         &                maskS(i,j,k,bi,bj) * delR(k) + utop
298    cgg    Add the barotropic velocity component.
299                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
300                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
301                          endif
302                        enddo
303    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
304                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
305         &                                      - utop / delR(1)
306                      enddo
307                    enddo
308                  enddo
309                endif
310              endif
311    
312    #endif /* ALLOW_CTRL_OBCS_BALANCE */
313    
314              do bj = jtlo,jthi
315                do bi = itlo,ithi
316                  do k = 1,nr
317                    do j = jmin,jmax
318                      xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
319    cgg     &                                        *   maskyz (j,k,bi,bj)
320                     enddo
321                enddo                enddo
322              enddo              enddo
323            enddo            enddo
324          enddo          endif
325    
326  c--     Add control to model variable.  c--     Add control to model variable.
327          do bj = jtlo,jthi          do bj = jtlo, jthi
328             do bi = itlo,ithi             do bi = itlo, ithi
329  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
330                do k = 1,nr                do k = 1,nr
331                   do j = 1,sny                   do j = 1,sny
332                        i = OB_Iw(j,bi,bj)
333                      if (iobcs .EQ. 1) then                      if (iobcs .EQ. 1) then
334                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
335       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
336       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
337                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
338       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
339                      else if (iobcs .EQ. 2) then                      else if (iobcs .EQ. 2) then
340                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
341       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
342       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
343                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
344       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
345                      else if (iobcs .EQ. 3) then                      else if (iobcs .EQ. 3) then
346                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu (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                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
350       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
351                      else if (iobcs .EQ. 4) then                      else if (iobcs .EQ. 4) then
352                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv (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                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
356       &                      *maskyz(j,k,bi,bj)       &                      *maskS(i,j,k,bi,bj)
357                      endif                      endif
358                   enddo                   enddo
359                enddo                enddo
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22