/[MITgcm]/MITgcm/pkg/ctrl/ctrl_getobcss.F
ViewVC logotype

Diff of /MITgcm/pkg/ctrl/ctrl_getobcss.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_getobcss(
9         I                             mytime,
10         I                             myiter,
11         I                             mythid
12         &                           )
13    
14    c     ==================================================================
15    c     SUBROUTINE ctrl_getobcss
16    c     ==================================================================
17    c
18    c     o Get southern 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     new flags: gebbie@mit.edu, 25 Jan 2003.
24    c
25    c     ==================================================================
26    c     SUBROUTINE ctrl_getobcss
27    c     ==================================================================
28    
29          implicit none
30    
31    #ifdef ALLOW_OBCSS_CONTROL
32    
33    c     == global variables ==
34    
35    #include "EEPARAMS.h"
36    #include "SIZE.h"
37    #include "PARAMS.h"
38    #include "GRID.h"
39    #include "OBCS.h"
40    
41    #include "ctrl.h"
42    #include "ctrl_dummy.h"
43    #include "optim.h"
44    
45    c     == routine arguments ==
46    
47          _RL     mytime
48          integer myiter
49          integer mythid
50    
51    c     == local variables ==
52    
53          integer bi,bj
54          integer i,j,k
55          integer itlo,ithi
56          integer jtlo,jthi
57          integer jmin,jmax
58          integer imin,imax
59          integer ilobcss
60          integer iobcs
61    
62          _RL     dummy
63          _RL     obcssfac
64          logical obcssfirst
65          logical obcsschanged
66          integer obcsscount0
67          integer obcsscount1
68          integer jp1
69    
70    cgg      _RL maskxz   (1-olx:snx+olx,nr,nsx,nsy)
71    
72          logical doglobalread
73          logical ladinit
74    
75          character*(80) fnameobcss
76    
77    cgg(  Variables for splitting barotropic/baroclinic vels.
78          _RL vbaro
79          _RL vtop
80    cgg)
81    
82    c     == external functions ==
83    
84          integer  ilnblnk
85          external ilnblnk
86    
87    
88    c     == end of interface ==
89    
90          jtlo = mybylo(mythid)
91          jthi = mybyhi(mythid)
92          itlo = mybxlo(mythid)
93          ithi = mybxhi(mythid)
94          jmin = 1-oly
95          jmax = sny+oly
96          imin = 1-olx
97          imax = snx+olx
98          jp1  = 1
99    
100    cgg(  Initialize variables for balancing volume flux.
101          vbaro = 0.d0
102          vtop = 0.d0
103    cgg)
104    
105    c--   Now, read the control vector.
106          doglobalread = .false.
107          ladinit      = .false.
108    
109          if (optimcycle .ge. 0) then
110            ilobcss=ilnblnk( xx_obcss_file )
111            write(fnameobcss(1:80),'(2a,i10.10)')
112         &       xx_obcss_file(1:ilobcss), '.', optimcycle
113          endif
114    
115    c--   Get the counters, flags, and the interpolation factor.
116          call ctrl_get_gen_rec(
117         I                   xx_obcssstartdate, xx_obcssperiod,
118         O                   obcssfac, obcssfirst, obcsschanged,
119         O                   obcsscount0,obcsscount1,
120         I                   mytime, myiter, mythid )
121    
122          do iobcs = 1,nobcs
123            if ( obcssfirst ) then
124              call active_read_xz( fnameobcss, tmpfldxz,
125         &                         (obcsscount0-1)*nobcs+iobcs,
126         &                         doglobalread, ladinit, optimcycle,
127         &                         mythid, xx_obcss_dummy )
128    
129              if ( optimcycle .gt. 0) then          
130                if (iobcs .eq. 3) then
131    cgg         Special attention is needed for the normal velocity.
132    cgg         For the north, this is the v velocity, iobcs = 4.
133    cgg         This is done on a columnwise basis here.
134                  do bj = jtlo,jthi
135                    do bi = itlo, ithi
136                      do i = imin,imax
137                        j = OB_Js(I,bi,bj)
138    
139    cgg         The barotropic velocity is stored in the level 1.
140                        vbaro = tmpfldxz(i,1,bi,bj)
141                        tmpfldxz(i,1,bi,bj) = 0.d0
142                        vtop = 0.d0
143    
144                        do k = 1,Nr
145    cgg    If cells are not full, this should be modified with hFac.
146    cgg    
147    cgg    The xx field (tmpfldxz) does not contain the velocity at the
148    cgg    surface level. This velocity is not independent; it must
149    cgg    exactly balance the volume flux, since we are dealing with
150    cgg    the baroclinic velocity structure..
151                          vtop = tmpfldxz(i,k,bi,bj)*
152         &                maskS(i,j+jp1,k,bi,bj) * delZ(k) + vtop
153    cgg    Add the barotropic velocity component.
154                          if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
155                            tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
156                          endif
157                        enddo
158    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
159                        tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
160         &                                      - vtop / delZ(1)
161                      enddo
162                    enddo
163                  enddo
164                endif
165    
166                if (iobcs .eq. 4) then
167    cgg         Special attention is needed for the normal velocity.
168    cgg         For the north, this is the v velocity, iobcs = 4.
169    cgg         This is done on a columnwise basis here.
170                  do bj = jtlo,jthi
171                    do bi = itlo, ithi
172                      do i = imin,imax
173                        j = OB_Js(I,bi,bj)
174    
175    cgg         The barotropic velocity is stored in the level 1.
176                        vbaro = tmpfldxz(i,1,bi,bj)
177                        tmpfldxz(i,1,bi,bj) = 0.d0
178                        vtop = 0.d0
179    
180                        do k = 1,Nr
181    cgg    If cells are not full, this should be modified with hFac.
182    cgg    
183    cgg    The xx field (tmpfldxz) does not contain the velocity at the
184    cgg    surface level. This velocity is not independent; it must
185    cgg    exactly balance the volume flux, since we are dealing with
186    cgg    the baroclinic velocity structure..
187                          vtop = tmpfldxz(i,k,bi,bj)*
188         &                maskW(i,j,k,bi,bj) * delZ(k) + vtop
189    cgg    Add the barotropic velocity component.
190                          if (maskW(i,j,k,bi,bj) .ne. 0.) then
191                            tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
192                          endif
193                        enddo
194    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
195                        tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
196         &                                      - vtop / delZ(1)
197                      enddo
198                    enddo
199                  enddo
200                endif
201              endif
202    
203              do bj = jtlo,jthi
204                do bi = itlo,ithi
205                  do k = 1,nr
206                    do i = imin,imax
207                      xx_obcss1(i,k,bi,bj,iobcs)  = tmpfldxz (i,k,bi,bj)
208    cgg     &                                        *   maskxz (i,k,bi,bj)
209                    enddo
210                  enddo
211                enddo
212              enddo
213            endif
214    
215            if ( (obcssfirst) .or. (obcsschanged)) then
216              
217              do bj = jtlo,jthi
218                do bi = itlo,ithi
219                  do k = 1,nr
220                    do i = imin,imax
221                      tmpfldxz(i,k,bi,bj) = xx_obcss1(i,k,bi,bj,iobcs)
222                    enddo
223                  enddo
224                enddo
225              enddo
226    
227              call exf_swapffields_xz( tmpfldxz2, tmpfldxz, mythid)
228    
229              do bj = jtlo,jthi
230                do bi = itlo,ithi
231                  do k = 1,nr
232                    do i = imin,imax
233                      xx_obcss0(i,k,bi,bj,iobcs) = tmpfldxz2(i,k,bi,bj)
234                    enddo
235                  enddo
236                enddo
237              enddo
238    
239              call active_read_xz( fnameobcss, tmpfldxz,
240         &                         (obcsscount1-1)*nobcs+iobcs,
241         &                         doglobalread, ladinit, optimcycle,
242         &                         mythid, xx_obcss_dummy )
243    
244              if ( optimcycle .gt. 0) then          
245                if (iobcs .eq. 3) then
246    cgg         Special attention is needed for the normal velocity.
247    cgg         For the north, this is the v velocity, iobcs = 4.
248    cgg         This is done on a columnwise basis here.
249                  do bj = jtlo,jthi
250                    do bi = itlo, ithi
251                      do i = imin,imax
252                        j = OB_Js(I,bi,bj)
253    
254    cgg         The barotropic velocity is stored in the level 1.
255                        vbaro = tmpfldxz(i,1,bi,bj)
256                        tmpfldxz(i,1,bi,bj) = 0.d0
257                        vtop = 0.d0
258    
259                        do k = 1,Nr
260    cgg    If cells are not full, this should be modified with hFac.
261    cgg    
262    cgg    The xx field (tmpfldxz) does not contain the velocity at the
263    cgg    surface level. This velocity is not independent; it must
264    cgg    exactly balance the volume flux, since we are dealing with
265    cgg    the baroclinic velocity structure..
266                          vtop = tmpfldxz(i,k,bi,bj)*
267         &                maskS(i,j+jp1,k,bi,bj) * delZ(k) + vtop
268    cgg    Add the barotropic velocity component.
269                          if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
270                            tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
271                          endif
272                        enddo
273    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
274                        tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
275         &                                      - vtop / delZ(1)
276                      enddo
277                    enddo
278                  enddo
279                endif
280    
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 i = imin,imax
288                        j = OB_Js(I,bi,bj)
289    
290    cgg         The barotropic velocity is stored in the level 1.
291                        vbaro = tmpfldxz(i,1,bi,bj)
292                        tmpfldxz(i,1,bi,bj) = 0.d0
293                        vtop = 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                          vtop = tmpfldxz(i,k,bi,bj)*
303         &                maskW(i,j,k,bi,bj) * delZ(k) + vtop
304    cgg    Add the barotropic velocity component.
305                          if (maskW(i,j,k,bi,bj) .ne. 0.) then
306                            tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro
307                          endif
308                        enddo
309    cgg    Compute the baroclinic velocity at level 1. Should balance flux.
310                        tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)
311         &                                      - vtop / 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 i = imin,imax
322                      xx_obcss1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
323    cgg     &                                        *   maskxz (i,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 i = 1,snx
336                        j = OB_Js(I,bi,bj)
337                        if (iobcs .EQ. 1) then
338                           OBSt(i,k,bi,bj) = OBSt (i,k,bi,bj)
339         &                 + obcssfac            *xx_obcss0(i,k,bi,bj,iobcs)
340         &                 + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
341                           OBSt(i,k,bi,bj) = OBSt(i,k,bi,bj)
342         &                      *maskS(i,j+jp1,k,bi,bj)
343                        else if (iobcs .EQ. 2) then
344                           OBSs(i,k,bi,bj) = OBSs (i,k,bi,bj)
345         &                 + obcssfac            *xx_obcss0(i,k,bi,bj,iobcs)
346         &                 + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
347                           OBSs(i,k,bi,bj) = OBSs(i,k,bi,bj)
348         &                      *maskS(i,j+jp1,k,bi,bj)
349                        else if (iobcs .EQ. 4) then
350                           OBSu(i,k,bi,bj) = OBSu (i,k,bi,bj)
351         &                 + obcssfac            *xx_obcss0(i,k,bi,bj,iobcs)
352         &                 + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
353                           OBSu(i,k,bi,bj) = OBSu(i,k,bi,bj)
354         &                      *maskW(i,j,k,bi,bj)
355                        else if (iobcs .EQ. 3) then
356                           OBSv(i,k,bi,bj) = OBSv (i,k,bi,bj)
357         &                 + obcssfac            *xx_obcss0(i,k,bi,bj,iobcs)
358         &                 + (1. _d 0 - obcssfac)*xx_obcss1(i,k,bi,bj,iobcs)
359                           OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
360         &                      *maskS(i,j+jp1,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_OBCSS_CONTROL undefined */
371    
372    c     == routine arguments ==
373    
374          _RL     mytime
375          integer myiter
376          integer mythid
377    
378    c--   CPP flag ALLOW_OBCSS_CONTROL undefined.
379    
380    #endif /* ALLOW_OBCSS_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