/[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.8 by jmc, Tue Oct 9 00:00:00 2007 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    #endif /* ALLOW_CTRL_OBCS_BALANCE */
206    
207          call active_read_yz( 'maskobcsw', maskyz,            do bj = jtlo,jthi
208       &                       iobcs,              do bi = itlo,ithi
209       &                       doglobalread, ladinit, 0,                do k = 1,nr
210       &                       mythid, dummy )                  do j = jmin,jmax
211                      xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
212          call active_read_yz( fnameobcsw, tmpfldyz,  cgg     &                                        *   maskyz (j,k,bi,bj)
213       &                       (obcswcount0-1)*nobcs+iobcs,                   enddo
214       &                       doglobalread, ladinit, optimcycle,                enddo
215       &                       mythid, xx_obcsw_dummy )              enddo
216              enddo
217          do bj = jtlo,jthi          endif
218            do bi = itlo,ithi  
219              do k = 1,nr          if ( (obcswfirst) .or. (obcswchanged)) then
220                do j = jmin,jmax  
221                  yz_obcs0(j,k,bi,bj)  = tmpfldyz (j,k,bi,bj)  cgg(    This is a terribly long way to do it. However, the dimensions do not exactly
222    cgg     match up. I will blame Fortran for the ugliness.
223    
224              do bj = jtlo,jthi
225                do bi = itlo,ithi
226                  do k = 1,nr
227                    do j = jmin,jmax
228                      tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)
229                    enddo
230                enddo                enddo
231              enddo              enddo
232            enddo            enddo
         enddo  
233    
234          call active_read_yz( fnameobcsw, tmpfldyz,            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)
235       &                       (obcswcount1-1)*nobcs+iobcs,  
236       &                       doglobalread, ladinit, optimcycle,            do bj = jtlo,jthi
237       &                       mythid, xx_obcsw_dummy )              do bi = itlo,ithi
238                  do k = 1,nr
239          do bj = jtlo,jthi                  do j = jmin,jmax
240            do bi = itlo,ithi                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
241              do k = 1,nr                  enddo
242                do j = jmin,jmax                enddo
243                  yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj)              enddo
244              enddo
245    
246              call active_read_yz( fnameobcsw, tmpfldyz,
247         &                         (obcswcount1-1)*nobcs+iobcs,
248         &                         doglobalread, ladinit, optimcycle,
249         &                         mythid, xx_obcsw_dummy )
250    
251    #ifdef ALLOW_CTRL_OBCS_BALANCE
252    
253              if ( optimcycle .gt. 0) then
254                if (iobcs .eq. 3) then
255    cgg         Special attention is needed for the normal velocity.
256    cgg         For the north, this is the v velocity, iobcs = 4.
257    cgg         This is done on a columnwise basis here.
258                  do bj = jtlo,jthi
259                    do bi = itlo, ithi
260                      do j = jmin,jmax
261                        i = OB_Iw(J,bi,bj)
262    
263    cgg         The barotropic velocity is stored in the level 1.
264                        ubaro = tmpfldyz(j,1,bi,bj)
265                        tmpfldyz(j,1,bi,bj) = 0.d0
266                        utop = 0.d0
267    
268                        do k = 1,Nr
269    cgg    If cells are not full, this should be modified with hFac.
270    cgg
271    cgg    The xx field (tmpfldxz) does not contain the velocity at the
272    cgg    surface level. This velocity is not independent; it must
273    cgg    exactly balance the volume flux, since we are dealing with
274    cgg    the baroclinic velocity structure..
275                          utop = tmpfldyz(j,k,bi,bj)*
276         &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop
277    cgg    Add the barotropic velocity component.
278                          if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
279                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
280                          endif
281                        enddo
282    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
283                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
284         &                                      - utop / delR(1)
285                      enddo
286                    enddo
287                  enddo
288                endif
289                if (iobcs .eq. 4) then
290    cgg         Special attention is needed for the normal velocity.
291    cgg         For the north, this is the v velocity, iobcs = 4.
292    cgg         This is done on a columnwise basis here.
293                  do bj = jtlo,jthi
294                    do bi = itlo, ithi
295                      do j = jmin,jmax
296                        i = OB_Iw(J,bi,bj)
297    
298    cgg         The barotropic velocity is stored in the level 1.
299                        ubaro = tmpfldyz(j,1,bi,bj)
300                        tmpfldyz(j,1,bi,bj) = 0.d0
301                        utop = 0.d0
302    
303                        do k = 1,Nr
304    cgg    If cells are not full, this should be modified with hFac.
305    cgg
306    cgg    The xx field (tmpfldxz) does not contain the velocity at the
307    cgg    surface level. This velocity is not independent; it must
308    cgg    exactly balance the volume flux, since we are dealing with
309    cgg    the baroclinic velocity structure..
310                          utop = tmpfldyz(j,k,bi,bj)*
311         &                maskS(i,j,k,bi,bj) * delR(k) + utop
312    cgg    Add the barotropic velocity component.
313                          if (maskS(i,j,k,bi,bj) .ne. 0.) then
314                            tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro
315                          endif
316                        enddo
317    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
318                        tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)
319         &                                      - utop / delR(1)
320                      enddo
321                    enddo
322                  enddo
323                endif
324              endif
325    
326    #endif /* ALLOW_CTRL_OBCS_BALANCE */
327    
328              do bj = jtlo,jthi
329                do bi = itlo,ithi
330                  do k = 1,nr
331                    do j = jmin,jmax
332                      xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
333    cgg     &                                        *   maskyz (j,k,bi,bj)
334                     enddo
335                enddo                enddo
336              enddo              enddo
337            enddo            enddo
338          enddo          endif
339    
340  c--     Add control to model variable.  c--     Add control to model variable.
341          do bj = jtlo,jthi          do bj = jtlo, jthi
342             do bi = itlo,ithi             do bi = itlo, ithi
343  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--        Calculate mask for tracer cells (0 => land, 1 => water).
344                do k = 1,nr                do k = 1,nr
345                   do j = 1,sny                   do j = 1,sny
346                        i = OB_Iw(j,bi,bj)
347                      if (iobcs .EQ. 1) then                      if (iobcs .EQ. 1) then
348                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
349       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
350       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
351                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
352       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
353                      else if (iobcs .EQ. 2) then                      else if (iobcs .EQ. 2) then
354                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
355       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
356       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
357                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
358       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
359                      else if (iobcs .EQ. 3) then                      else if (iobcs .EQ. 3) then
360                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
361       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
362       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
363                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
364       &                      *maskyz(j,k,bi,bj)       &                      *maskW(i+ip1,j,k,bi,bj)
365                      else if (iobcs .EQ. 4) then                      else if (iobcs .EQ. 4) then
366                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
367       &                      + obcswfac            *yz_obcs0(j,k,bi,bj)       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
368       &                      + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj)       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
369                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
370       &                      *maskyz(j,k,bi,bj)       &                      *maskS(i,j,k,bi,bj)
371                      endif                      endif
372                   enddo                   enddo
373                enddo                enddo
Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22