/[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 by heimbach, Tue Feb 5 20:23:58 2002 UTC revision 1.2 by heimbach, Tue Jun 24 16:07:06 2003 UTC
# Line 0  Line 1 
1    
2    #include "CTRL_CPPOPTIONS.h"
3    #ifdef ALLOW_OBCS
4    # include "OBCS_OPTIONS.h"
5    #endif
6    
7    
8          subroutine ctrl_getobcsw(
9         I                             mytime,
10         I                             myiter,
11         I                             mythid
12         &                           )
13    
14    c     ==================================================================
15    c     SUBROUTINE ctrl_getobcsw
16    c     ==================================================================
17    c
18    c     o Get western obc of the control vector and add it
19    c       to dyn. fields
20    c
21    c     started: heimbach@mit.edu, 29-Aug-2001
22    c
23    c     modified: gebbie@mit.edu, 18-Mar-2003
24    c     ==================================================================
25    c     SUBROUTINE ctrl_getobcsw
26    c     ==================================================================
27    
28          implicit none
29    
30    #ifdef ALLOW_OBCSW_CONTROL
31    
32    c     == global variables ==
33    
34    #include "EEPARAMS.h"
35    #include "SIZE.h"
36    #include "PARAMS.h"
37    #include "GRID.h"
38    #include "OBCS.h"
39    
40    #include "ctrl.h"
41    #include "ctrl_dummy.h"
42    #include "optim.h"
43    
44    c     == routine arguments ==
45    
46          _RL     mytime
47          integer myiter
48          integer mythid
49    
50    c     == local variables ==
51    
52          integer bi,bj
53          integer i,j,k
54          integer itlo,ithi
55          integer jtlo,jthi
56          integer jmin,jmax
57          integer imin,imax
58          integer ilobcsw
59          integer iobcs
60          integer ip1
61    
62          _RL     dummy
63          _RL     obcswfac
64          logical obcswfirst
65          logical obcswchanged
66          integer obcswcount0
67          integer obcswcount1
68    
69    cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
70    
71          logical doglobalread
72          logical ladinit
73    
74          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 ==
82    
83          integer  ilnblnk
84          external ilnblnk
85    
86    
87    c     == end of interface ==
88    
89          jtlo = mybylo(mythid)
90          jthi = mybyhi(mythid)
91          itlo = mybxlo(mythid)
92          ithi = mybxhi(mythid)
93          jmin = 1-oly
94          jmax = sny+oly
95          imin = 1-olx
96          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.
105          doglobalread = .false.
106          ladinit      = .false.
107    
108          if (optimcycle .ge. 0) then
109            ilobcsw=ilnblnk( xx_obcsw_file )
110            write(fnameobcsw(1:80),'(2a,i10.10)')
111         &       xx_obcsw_file(1:ilobcsw), '.', optimcycle
112          endif
113    
114    c--   Get the counters, flags, and the interpolation factor.
115          call ctrl_get_gen_rec(
116         I                   xx_obcswstartdate, xx_obcswperiod,
117         O                   obcswfac, obcswfirst, obcswchanged,
118         O                   obcswcount0,obcswcount1,
119         I                   mytime, myiter, mythid )
120    
121          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              do bj = jtlo,jthi
231                do bi = itlo,ithi
232                  do k = 1,nr
233                    do j = jmin,jmax
234                      xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)
235                    enddo
236                  enddo
237                enddo
238              enddo
239    
240              call active_read_yz( fnameobcsw, tmpfldyz,
241         &                         (obcswcount1-1)*nobcs+iobcs,
242         &                         doglobalread, ladinit, optimcycle,
243         &                         mythid, xx_obcsw_dummy )
244    
245              if ( optimcycle .gt. 0) then          
246                if (iobcs .eq. 3) then
247    cgg         Special attention is needed for the normal velocity.
248    cgg         For the north, this is the v velocity, iobcs = 4.
249    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
326                enddo
327              enddo
328            endif
329    
330    c--     Add control to model variable.
331            do bj = jtlo, jthi
332               do bi = itlo, ithi
333    c--        Calculate mask for tracer cells (0 => land, 1 => water).
334                  do k = 1,nr
335                     do j = 1,sny
336                        i = OB_Iw(j,bi,bj)
337                        if (iobcs .EQ. 1) then
338                           OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
339         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
340         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
341                           OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
342         &                      *maskW(i+ip1,j,k,bi,bj)
343                        else if (iobcs .EQ. 2) then
344                           OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
345         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
346         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
347                           OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
348         &                      *maskW(i+ip1,j,k,bi,bj)
349                        else if (iobcs .EQ. 3) then
350                           OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
351         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
352         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
353                           OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
354         &                      *maskW(i+ip1,j,k,bi,bj)
355                        else if (iobcs .EQ. 4) then
356                           OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
357         &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
358         &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
359                           OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
360         &                      *maskS(i,j,k,bi,bj)
361                        endif
362                     enddo
363                  enddo
364               enddo
365            enddo
366    
367    C--   End over iobcs loop
368          enddo
369    
370    #else /* ALLOW_OBCSW_CONTROL undefined */
371    
372    c     == routine arguments ==
373    
374          _RL     mytime
375          integer myiter
376          integer mythid
377    
378    c--   CPP flag ALLOW_OBCSW_CONTROL undefined.
379    
380    #endif /* ALLOW_OBCSW_CONTROL */
381    
382          end
383    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22