/[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.2 by heimbach, Tue Jun 24 16:07:06 2003 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              if ( optimcycle .gt. 0) then          
129                if (iobcs .eq. 3) then
130    cgg         Special attention is needed for the normal velocity.
131    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
134                    do bi = itlo, ithi
135                      do j = jmin,jmax
136                        i = OB_Iw(J,bi,bj)
137    
138    cgg         The barotropic velocity is stored in the level 1.
139                        ubaro = tmpfldyz(j,1,bi,bj)
140                        tmpfldyz(j,1,bi,bj) = 0.d0
141                        utop = 0.d0
142    
143                        do k = 1,Nr
144    cgg    If cells are not full, this should be modified with hFac.
145    cgg    
146    cgg    The xx field (tmpfldxz) does not contain the velocity at the
147    cgg    surface level. This velocity is not independent; it must
148    cgg    exactly balance the volume flux, since we are dealing with
149    cgg    the baroclinic velocity structure..
150                          utop = tmpfldyz(j,k,bi,bj)*
151         &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop
152    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
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
161                    enddo
162                  enddo
163                endif
164                if (iobcs .eq. 4) then
165    cgg         Special attention is needed for the normal velocity.
166    cgg         For the north, this is the v velocity, iobcs = 4.
167    cgg         This is done on a columnwise basis here.
168                  do bj = jtlo,jthi
169                    do bi = itlo, ithi
170                      do j = jmin,jmax
171                        i = OB_Iw(J,bi,bj)
172    
173    cgg         The barotropic velocity is stored in the level 1.
174                        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
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
196                    enddo
197                  enddo
198                endif
199              endif
200    
201              do bj = jtlo,jthi
202                do bi = itlo,ithi
203                  do k = 1,nr
204                    do j = jmin,jmax
205                      xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
206    cgg     &                                        *   maskyz (j,k,bi,bj)
207                     enddo
208                  enddo
209                enddo
210              enddo
211            endif
212    
213            if ( (obcswfirst) .or. (obcswchanged)) then
214              
215    cgg(    This is a terribly long way to do it. However, the dimensions don't exactly
216    cgg     match up. I will blame Fortran for the ugliness.
217    
218              do bj = jtlo,jthi
219                do bi = itlo,ithi
220                  do k = 1,nr
221                    do j = jmin,jmax
222                      tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
223                    enddo
224                  enddo
225                enddo
226              enddo
227    
228              call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
229    
230          call active_read_yz( 'maskobcsw', maskyz,            do bj = jtlo,jthi
231       &                       iobcs,              do bi = itlo,ithi
232       &                       doglobalread, ladinit, 0,                do k = 1,nr
233       &                       mythid, dummy )                  do j = jmin,jmax
234                      xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
235          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)  
236                enddo                enddo
237              enddo              enddo
238            enddo            enddo
         enddo  
239    
240          call active_read_yz( fnameobcsw, tmpfldyz,            call active_read_yz( fnameobcsw, tmpfldyz,
241       &                       (obcswcount1-1)*nobcs+iobcs,       &                         (obcswcount1-1)*nobcs+iobcs,
242       &                       doglobalread, ladinit, optimcycle,       &                         doglobalread, ladinit, optimcycle,
243       &                       mythid, xx_obcsw_dummy )       &                         mythid, xx_obcsw_dummy )
244    
245          do bj = jtlo,jthi            if ( optimcycle .gt. 0) then          
246            do bi = itlo,ithi              if (iobcs .eq. 3) then
247              do k = 1,nr  cgg         Special attention is needed for the normal velocity.
248                do j = jmin,jmax  cgg         For the north, this is the v velocity, iobcs = 4.
249                  yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj)  cgg         This is done on a columnwise basis here.
250                  do bj = jtlo,jthi
251                    do bi = itlo, ithi
252                      do j = jmin,jmax
253                        i = OB_Iw(J,bi,bj)
254    
255    cgg         The barotropic velocity is stored in the level 1.
256                        ubaro = tmpfldyz(j,1,bi,bj)
257                        tmpfldyz(j,1,bi,bj) = 0.d0
258                        utop = 0.d0
259    
260                        do k = 1,Nr
261    cgg    If cells are not full, this should be modified with hFac.
262    cgg    
263    cgg    The xx field (tmpfldxz) does not contain the velocity at the
264    cgg    surface level. This velocity is not independent; it must
265    cgg    exactly balance the volume flux, since we are dealing with
266    cgg    the baroclinic velocity structure..
267                          utop = tmpfldyz(j,k,bi,bj)*
268         &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop
269    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
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
278                    enddo
279                  enddo
280                endif
281                if (iobcs .eq. 4) then
282    cgg         Special attention is needed for the normal velocity.
283    cgg         For the north, this is the v velocity, iobcs = 4.
284    cgg         This is done on a columnwise basis here.
285                  do bj = jtlo,jthi
286                    do bi = itlo, ithi
287                      do j = jmin,jmax
288                        i = OB_Iw(J,bi,bj)
289    
290    cgg         The barotropic velocity is stored in the level 1.
291                        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
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
313                    enddo
314                  enddo
315                endif
316              endif
317    
318              do bj = jtlo,jthi
319                do bi = itlo,ithi
320                  do k = 1,nr
321                    do j = jmin,jmax
322                      xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
323    cgg     &                                        *   maskyz (j,k,bi,bj)
324                     enddo
325                enddo                enddo
326              enddo              enddo
327            enddo            enddo
328          enddo          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
336                        i = OB_Iw(j,bi,bj)
337                      if (iobcs .EQ. 1) then                      if (iobcs .EQ. 1) then
338                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
339       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
340       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
341                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
342       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
343                      else if (iobcs .EQ. 2) then                      else if (iobcs .EQ. 2) then
344                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
345       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
346       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
347                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
348       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
349                      else if (iobcs .EQ. 3) then                      else if (iobcs .EQ. 3) then
350                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
351       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
352       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
353                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
354       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
355                      else if (iobcs .EQ. 4) then                      else if (iobcs .EQ. 4) then
356                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
357       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
358       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
359                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
360       &                      *maskyz(j,k,bi,bj)       &                      *maskS(i,j,k,bi,bj)
361                      endif                      endif
362                   enddo                   enddo
363                enddo                enddo
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22