/[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.8 by jmc, Tue Oct 9 00:00:00 2007 UTC revision 1.13 by jmc, Tue May 24 20:48:28 2011 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  # include "OBCS_OPTIONS.h"  # include "OBCS_OPTIONS.h"
7  #endif  #endif
8    
   
9        subroutine ctrl_getobcsw(        subroutine ctrl_getobcsw(
10       I                             mytime,       I                             mytime,
11       I                             myiter,       I                             myiter,
# Line 29  c     ================================== Line 28  c     ==================================
28    
29        implicit none        implicit none
30    
 #ifdef ALLOW_OBCSW_CONTROL  
   
31  c     == global variables ==  c     == global variables ==
32    #ifdef ALLOW_OBCSW_CONTROL
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36  #include "GRID.h"  #include "GRID.h"
37  #include "OBCS.h"  c#include "OBCS_PARAMS.h"
38    #include "OBCS_GRID.h"
39    #include "OBCS_FIELDS.h"
40    
41  #include "ctrl.h"  #include "ctrl.h"
42  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
43  #include "optim.h"  #include "optim.h"
44    #endif /* ALLOW_OBCSW_CONTROL */
45    
46  c     == routine arguments ==  c     == routine arguments ==
   
47        _RL     mytime        _RL     mytime
48        integer myiter        integer myiter
49        integer mythid        integer mythid
50    
51    #ifdef ALLOW_OBCSW_CONTROL
52  c     == local variables ==  c     == local variables ==
53    
54        integer bi,bj        integer bi,bj
# Line 69  c     == local variables == Line 69  c     == local variables ==
69        integer obcswcount1        integer obcswcount1
70    
71  cgg      _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.  #ifdef ALLOW_OBCS_CONTROL_MODES
80        _RL ubaro        integer nk,nz
81        _RL utop        _RL     tmpz (nr,nsx,nsy)
82  cgg)        _RL     stmp
83    #endif
84    
85  c     == external functions ==  c     == external functions ==
86    
# Line 98  c     == end of interface == Line 100  c     == end of interface ==
100        imax = snx+olx        imax = snx+olx
101        ip1  = 1        ip1  = 1
102    
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       utop = 0.d0  
 cgg)  
   
103  c--   Now, read the control vector.  c--   Now, read the control vector.
104        doglobalread = .false.        doglobalread = .false.
105        ladinit      = .false.        ladinit      = .false.
106    
107        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
108          ilobcsw=ilnblnk( xx_obcsw_file )         ilobcsw=ilnblnk( xx_obcsw_file )
109          write(fnameobcsw(1:80),'(2a,i10.10)')         write(fnameobcsw(1:80),'(2a,i10.10)')
110       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle       &      xx_obcsw_file(1:ilobcsw), '.', optimcycle
111        endif        endif
112    
113  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
# Line 121  c--   Get the counters, flags, and the i Line 118  c--   Get the counters, flags, and the i
118       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
119    
120        do iobcs = 1,nobcs        do iobcs = 1,nobcs
121          if ( obcswfirst ) then         if ( obcswfirst ) then
122            call active_read_yz( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
123       &                         (obcswcount0-1)*nobcs+iobcs,       &                       (obcswcount0-1)*nobcs+iobcs,
124       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
125       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
126    
127  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
128             do bi = itlo,ithi
129            if ( optimcycle .gt. 0) then  #ifdef ALLOW_OBCS_CONTROL_MODES
130              if (iobcs .eq. 3) then            if (iobcs .gt. 2) then
131  cgg         Special attention is needed for the normal velocity.             do j = jmin,jmax
132  cgg         For the north, this is the v velocity, iobcs = 4.              i = OB_Iw(j,bi,bj)
133  cgg         This is done on a columnwise basis here.  cih    Determine number of open vertical layers.
134                do bj = jtlo,jthi              nz = 0
135                  do bi = itlo, ithi              do k = 1,Nr
136                    do j = jmin,jmax                if (iobcs .eq. 3) then
137                      i = OB_Iw(J,bi,bj)                  nz = nz + maskW(i+ip1,j,k,bi,bj)
138                  else
139  cgg         The barotropic velocity is stored in the level 1.                  nz = nz + maskS(i,j,k,bi,bj)
140                      ubaro = tmpfldyz(j,1,bi,bj)                endif
141                      tmpfldyz(j,1,bi,bj) = 0.d0              end do
142                      utop = 0.d0  cih    Compute absolute velocities from the barotropic-baroclinic modes.
143                do k = 1,Nr
144                      do k = 1,Nr               if (k.le.nz) then
145  cgg    If cells are not full, this should be modified with hFac.                stmp = 0.
146  cgg                do nk = 1,nz
147  cgg    The xx field (tmpfldxz) does not contain the velocity at the                 stmp = stmp +
148  cgg    surface level. This velocity is not independent; it must       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
149  cgg    exactly balance the volume flux, since we are dealing with                end do
150  cgg    the baroclinic velocity structure..                 tmpz(k,bi,bj) = stmp
151                        utop = tmpfldyz(j,k,bi,bj)*               else
152       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop                tmpz(k,bi,bj) = 0.
153  cgg    Add the barotropic velocity component.               end if
                       if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                   tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
             if (iobcs .eq. 4) then  
 cgg         Special attention is needed for the normal velocity.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
 cgg         This is done on a columnwise basis here.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(J,bi,bj)  
   
 cgg         The barotropic velocity is stored in the level 1.  
                     ubaro = tmpfldyz(j,1,bi,bj)  
                     tmpfldyz(j,1,bi,bj) = 0.d0  
                     utop = 0.d0  
   
                     do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg  
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    surface level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                       utop = tmpfldyz(j,k,bi,bj)*  
      &                maskS(i,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                   tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
   
 #endif /* ALLOW_CTRL_OBCS_BALANCE */  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
154              enddo              enddo
155                do k = 1,Nr
156                  if (iobcs .eq. 3) then
157                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
158         &            *recip_hFacW(i+ip1,j,k,bi,bj)
159                  else
160                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
161         &            *recip_hFacS(i,j,k,bi,bj)
162                  endif
163                end do
164               enddo
165              endif
166    #endif
167              do k = 1,nr
168               do j = jmin,jmax
169                xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
170    cgg   &                                        *   maskyz (j,k,bi,bj)
171               enddo
172            enddo            enddo
173          endif           enddo
174            enddo
175          if ( (obcswfirst) .or. (obcswchanged)) then         endif
176    
177  cgg(    This is a terribly long way to do it. However, the dimensions do not exactly         if ( (obcswfirst) .or. (obcswchanged)) then
 cgg     match up. I will blame Fortran for the ugliness.  
178    
179            do bj = jtlo,jthi          do bj = jtlo,jthi
180              do bi = itlo,ithi           do bi = itlo,ithi
181                do k = 1,nr            do k = 1,nr
182                  do j = jmin,jmax             do j = jmin,jmax
183                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
184                  enddo              tmpfldyz (j,k,bi,bj)       = 0. _d 0
185                enddo             enddo
             enddo  
186            enddo            enddo
187             enddo
188            enddo
189    
190            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)          call active_read_yz( fnameobcsw, tmpfldyz,
191         &                       (obcswcount1-1)*nobcs+iobcs,
192            do bj = jtlo,jthi       &                       doglobalread, ladinit, optimcycle,
193              do bi = itlo,ithi       &                       mythid, xx_obcsw_dummy )
194                do k = 1,nr  
195                  do j = jmin,jmax          do bj = jtlo,jthi
196                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)           do bi = itlo,ithi
197                  enddo  #ifdef ALLOW_OBCS_CONTROL_MODES
198                enddo            if (iobcs .gt. 2) then
199               do j = jmin,jmax
200                i = OB_Iw(j,bi,bj)
201    cih    Determine number of open vertical layers.
202                nz = 0
203                do k = 1,Nr
204                  if (iobcs .eq. 3) then
205                    nz = nz + maskW(i+ip1,j,k,bi,bj)
206                  else
207                    nz = nz + maskS(i,j,k,bi,bj)
208                  endif
209                end do
210    cih    Compute absolute velocities from the barotropic-baroclinic modes.
211                do k = 1,Nr
212                 if (k.le.nz) then
213                  stmp = 0.
214                  do nk = 1,nz
215                   stmp = stmp +
216         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
217                  end do
218                   tmpz(k,bi,bj) = stmp
219                 else
220                   tmpz(k,bi,bj) = 0.
221                 end if
222              enddo              enddo
223            enddo              do k = 1,Nr
224                  if (iobcs .eq. 3) then
225            call active_read_yz( fnameobcsw, tmpfldyz,                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
226       &                         (obcswcount1-1)*nobcs+iobcs,       &            *recip_hFacW(i+ip1,j,k,bi,bj)
227       &                         doglobalread, ladinit, optimcycle,                else
228       &                         mythid, xx_obcsw_dummy )                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
229         &            *recip_hFacS(i,j,k,bi,bj)
230  #ifdef ALLOW_CTRL_OBCS_BALANCE                endif
231                end do
232            if ( optimcycle .gt. 0) then             enddo
             if (iobcs .eq. 3) then  
 cgg         Special attention is needed for the normal velocity.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
 cgg         This is done on a columnwise basis here.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(J,bi,bj)  
   
 cgg         The barotropic velocity is stored in the level 1.  
                     ubaro = tmpfldyz(j,1,bi,bj)  
                     tmpfldyz(j,1,bi,bj) = 0.d0  
                     utop = 0.d0  
   
                     do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg  
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    surface level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                       utop = tmpfldyz(j,k,bi,bj)*  
      &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
             if (iobcs .eq. 4) then  
 cgg         Special attention is needed for the normal velocity.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
 cgg         This is done on a columnwise basis here.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(J,bi,bj)  
   
 cgg         The barotropic velocity is stored in the level 1.  
                     ubaro = tmpfldyz(j,1,bi,bj)  
                     tmpfldyz(j,1,bi,bj) = 0.d0  
                     utop = 0.d0  
   
                     do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg  
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    surface level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                       utop = tmpfldyz(j,k,bi,bj)*  
      &                maskS(i,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
233            endif            endif
234    #endif
235  #endif /* ALLOW_CTRL_OBCS_BALANCE */            do k = 1,nr
236               do j = jmin,jmax
237            do bj = jtlo,jthi              xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
238              do bi = itlo,ithi  cgg   &                                        *   maskyz (j,k,bi,bj)
239                do k = 1,nr             enddo
                 do j = jmin,jmax  
                   xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
             enddo  
240            enddo            enddo
241          endif           enddo
242            enddo
243           endif
244    
245  c--     Add control to model variable.  c--   Add control to model variable.
246          do bj = jtlo, jthi         do bj = jtlo, jthi
247             do bi = itlo, ithi          do bi = itlo, ithi
248  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--   Calculate mask for tracer cells (0 => land, 1 => water).
249                do k = 1,nr           do k = 1,nr
250                   do j = 1,sny            do j = 1,sny
251                      i = OB_Iw(j,bi,bj)             i = OB_Iw(j,bi,bj)
252                      if (iobcs .EQ. 1) then             if (iobcs .EQ. 1) then
253                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
254       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
255       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
256                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
257       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
258                      else if (iobcs .EQ. 2) then             else if (iobcs .EQ. 2) then
259                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
260       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
261       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
262                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
263       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
264                      else if (iobcs .EQ. 3) then             else if (iobcs .EQ. 3) then
265                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
266       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
267       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
268                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
269       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
270                      else if (iobcs .EQ. 4) then             else if (iobcs .EQ. 4) then
271                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
272       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
273       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
274                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
275       &                      *maskS(i,j,k,bi,bj)       &           *maskS(i,j,k,bi,bj)
276                      endif             endif
277                   enddo            enddo
278                enddo           enddo
            enddo  
279          enddo          enddo
280           enddo
281    
282  C--   End over iobcs loop  C--   End over iobcs loop
283        enddo        enddo
284    
 #else /* ALLOW_OBCSW_CONTROL undefined */  
   
 c     == routine arguments ==  
   
       _RL     mytime  
       integer myiter  
       integer mythid  
   
 c--   CPP flag ALLOW_OBCSW_CONTROL undefined.  
   
285  #endif /* ALLOW_OBCSW_CONTROL */  #endif /* ALLOW_OBCSW_CONTROL */
286    
287          return
288        end        end
   
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22