/[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.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 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          _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          cgg(  Variables for splitting barotropic/baroclinic vels.
80          _RL ubaro
81          _RL utop
82    cgg)
83    
84  c     == external functions ==  c     == external functions ==
85    
# Line 89  c     == end of interface == Line 97  c     == end of interface ==
97        jmax = sny+oly        jmax = sny+oly
98        imin = 1-olx        imin = 1-olx
99        imax = snx+olx        imax = snx+olx
100          ip1  = 1
101    
102    cgg(  Initialize variables for balancing volume flux.
103          ubaro = 0.d0
104          utop = 0.d0
105    cgg)
106    
107  c--   Now, read the control vector.  c--   Now, read the control vector.
108        doglobalread = .false.        doglobalread = .false.
# Line 96  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
127            if ( obcswfirst ) then
128              call active_read_yz( fnameobcsw, tmpfldyz,
129         &                         (obcswcount0-1)*nobcs+iobcs,
130         &                         doglobalread, ladinit, optimcycle,
131         &                         mythid, xx_obcsw_dummy )
132    
133    #ifdef ALLOW_CTRL_OBCS_BALANCE
134    
135              if ( optimcycle .gt. 0) then
136                if (iobcs .eq. 3) then
137    cgg         Special attention is needed for the normal velocity.
138    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
141                    do bi = itlo, ithi
142                      do j = jmin,jmax
143                        i = OB_Iw(J,bi,bj)
144    
145    cgg         The barotropic velocity is stored in the level 1.
146                        ubaro = tmpfldyz(j,1,bi,bj)
147                        tmpfldyz(j,1,bi,bj) = 0.d0
148                        utop = 0.d0
149    
150                        do k = 1,Nr
151    cgg    If cells are not full, this should be modified with hFac.
152    cgg
153    cgg    The xx field (tmpfldyz) does not contain the velocity at the
154    cgg    surface level. This velocity is not independent; it must
155    cgg    exactly balance the volume flux, since we are dealing with
156    cgg    the baroclinic velocity structure..
157                          utop = tmpfldyz(j,k,bi,bj)*
158         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
159    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
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
168                    enddo
169                  enddo
170                endif
171                if (iobcs .eq. 4) then
172    cgg         Special attention is needed for the normal velocity.
173    cgg         For the north, this is the v velocity, iobcs = 4.
174    cgg         This is done on a columnwise basis here.
175                  do bj = jtlo,jthi
176                    do bi = itlo, ithi
177                      do j = jmin,jmax
178                        i = OB_Iw(J,bi,bj)
179    
180    cgg         The barotropic velocity is stored in the level 1.
181                        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
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
203                    enddo
204                  enddo
205                endif
206              endif
207    
208          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 )  
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                  yz_obcs0(j,k,bi,bj)  = 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
217                enddo                enddo
218              enddo              enddo
219            enddo            enddo
220          enddo          endif
221    
222          call active_read_yz( fnameobcsw, tmpfldyz,          if ( (obcswfirst) .or. (obcswchanged)) then
      &                       (obcswcount1-1)*nobcs+iobcs,  
      &                       doglobalread, ladinit, optimcycle,  
      &                       mythid, xx_obcsw_dummy )  
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                  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)
229                enddo                tmpfldyz (j,k,bi,bj)       = 0. _d 0
230                 enddo
231              enddo              enddo
232               enddo
233            enddo            enddo
         enddo  
234    
235  c--     Add control to model variable.            call active_read_yz( fnameobcsw, tmpfldyz,
236          do bj = jtlo,jthi       &                         (obcswcount1-1)*nobcs+iobcs,
237             do bi = itlo,ithi       &                         doglobalread, ladinit, optimcycle,
238  c--        Calculate mask for tracer cells (0 => land, 1 => water).       &                         mythid, xx_obcsw_dummy )
239    
240    #ifdef ALLOW_CTRL_OBCS_BALANCE
241    
242              if ( optimcycle .gt. 0) then
243                if (iobcs .eq. 3) then
244    cgg         Special attention is needed for the normal velocity.
245    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
248                    do bi = itlo, ithi
249                      do j = jmin,jmax
250                        i = OB_Iw(J,bi,bj)
251    
252    cgg         The barotropic velocity is stored in the level 1.
253                        ubaro = tmpfldyz(j,1,bi,bj)
254                        tmpfldyz(j,1,bi,bj) = 0.d0
255                        utop = 0.d0
256    
257                        do k = 1,Nr
258    cgg    If cells are not full, this should be modified with hFac.
259    cgg
260    cgg    The xx field (tmpfldyz) does not contain the velocity at the
261    cgg    surface level. This velocity is not independent; it must
262    cgg    exactly balance the volume flux, since we are dealing with
263    cgg    the baroclinic velocity structure..
264                          utop = tmpfldyz(j,k,bi,bj)*
265         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
266    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
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
275                    enddo
276                  enddo
277                endif
278                if (iobcs .eq. 4) then
279    cgg         Special attention is needed for the normal velocity.
280    cgg         For the north, this is the v velocity, iobcs = 4.
281    cgg         This is done on a columnwise basis here.
282                  do bj = jtlo,jthi
283                    do bi = itlo, ithi
284                      do j = jmin,jmax
285                        i = OB_Iw(J,bi,bj)
286    
287    cgg         The barotropic velocity is stored in the level 1.
288                        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
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
310                    enddo
311                  enddo
312                endif
313              endif
314    
315    #endif /* ALLOW_CTRL_OBCS_BALANCE */
316    
317              do bj = jtlo,jthi
318                do bi = itlo,ithi
319                do k = 1,nr                do k = 1,nr
320                   do j = 1,sny                  do j = jmin,jmax
321                      if (iobcs .EQ. 1) then                    xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
322                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)  cgg     &                                        *   maskyz (j,k,bi,bj)
      &                      + obcswfac            *yz_obcs0(j,k,bi,bj)  
      &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)  
                        OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)  
      &                      *maskyz(j,k,bi,bj)  
                     else if (iobcs .EQ. 2) then  
                        OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)  
      &                      + obcswfac            *yz_obcs0(j,k,bi,bj)  
      &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)  
                        OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)  
      &                      *maskyz(j,k,bi,bj)  
                     else if (iobcs .EQ. 3) then  
                        OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)  
      &                      + obcswfac            *yz_obcs0(j,k,bi,bj)  
      &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)  
                        OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)  
      &                      *maskyz(j,k,bi,bj)  
                     else if (iobcs .EQ. 4) then  
                        OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)  
      &                      + obcswfac            *yz_obcs0(j,k,bi,bj)  
      &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)  
                        OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)  
      &                      *maskyz(j,k,bi,bj)  
                     endif  
323                   enddo                   enddo
324                enddo                enddo
325                enddo
326              enddo
327            endif
328    
329    c--   Add control to model variable.
330            do bj = jtlo, jthi
331             do bi = itlo, ithi
332    c--   Calculate mask for tracer cells (0 => land, 1 => water).
333              do k = 1,nr
334               do j = 1,sny
335                i = OB_Iw(j,bi,bj)
336                if (iobcs .EQ. 1) then
337                 OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
338         &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
339         &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
340                 OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
341         &            *maskW(i+ip1,j,k,bi,bj)
342                else if (iobcs .EQ. 2) then
343                 OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
344         &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
345         &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
346                 OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
347         &            *maskW(i+ip1,j,k,bi,bj)
348                else if (iobcs .EQ. 3) then
349                 OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
350         &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
351         &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
352                 OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
353         &            *maskW(i+ip1,j,k,bi,bj)
354                else if (iobcs .EQ. 4) then
355                 OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
356         &            + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
357         &            + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
358                 OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
359         &            *maskS(i,j,k,bi,bj)
360                endif
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.1  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22