/[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.2 by heimbach, Thu May 30 19:56:44 2002 UTC revision 1.12 by mmazloff, Wed Apr 20 19:14:07 2011 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 20  c       to dyn. fields Line 22  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 66  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  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_OBCS_CONTROL_MODES
80  cgg(  Variables for balancing volume flux.        integer nk,nz
81        _RL ubaro        _RL     tmpz (nr,nsx,nsy)
82        _RL ubarocount        _RL     stmp
 cgg)  
83  #endif  #endif
84    
85  c     == external functions ==  c     == external functions ==
# Line 97  c     == end of interface == Line 100  c     == end of interface ==
100        imax = snx+olx        imax = snx+olx
101        ip1  = 1        ip1  = 1
102    
 #ifdef BALANCE_CONTROL_VOLFLUX  
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       ubarocount = 0.d0  
 cgg)  
 #endif  
   
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
       else  
         print*  
         print*,' ctrl_getobcsw: optimcycle not set correctly.'  
         print*,' ... stopped in ctrl_getobcsw.'  
111        endif        endif
112    
113  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
114        call ctrl_GetRec( 'xx_obcsw',        call ctrl_get_gen_rec(
115         I                   xx_obcswstartdate, xx_obcswperiod,
116       O                   obcswfac, obcswfirst, obcswchanged,       O                   obcswfac, obcswfirst, obcswchanged,
117       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
118       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
119    
120        do iobcs = 1,nobcs        do iobcs = 1,nobcs
121                 if ( obcswfirst ) then
122  cgg        if ( (obcswfirst) .or. (obcswchanged) ) then          call active_read_yz( fnameobcsw, tmpfldyz,
123  cgg          call active_read_yz( 'maskobcsw', maskyz,       &                       (obcswcount0-1)*nobcs+iobcs,
124  cgg     &                         iobcs,       &                       doglobalread, ladinit, optimcycle,
125  cgg     &                         doglobalread, ladinit, 0,       &                       mythid, xx_obcsw_dummy )
 cgg     &                         mythid, dummy )  
 cgg        endif  
   
         if ( obcswfirst ) then  
           call active_read_yz( fnameobcsw, tmpfldyz,  
      &                         (obcswcount0-1)*nobcs+iobcs,  
      &                         doglobalread, ladinit, optimcycle,  
      &                         mythid, xx_obcsw_dummy )  
   
 #ifdef BALANCE_CONTROL_VOLFLUX  
           if ( optimcycle .gt. 0) then            
             if (iobcs .eq. 3) then  
 cgg(        Make sure that the xx velocity field has a balanced net volume flux.  
 cgg         Find the barotropic flow normal to the boundary.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(j,bi,bj)  
                     ubaro = 0.d0  
                     ubarocount = 0.d0  
                     do k = 1,nr  
 cgg(   If cells are not full, this should be modified with hFac.  
                       ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)  
      &                      * delZ(k) + ubaro  
                       ubarocount = maskW(i+ip1,j,k,bi,bj)  
      &                           * delZ(k) +ubarocount  
                     enddo  
                     if (ubarocount .eq. 0.) then  
                       ubaro = 0.  
                     else  
                       ubaro = ubaro / ubarocount  
                     endif  
                     do k = 1,nr  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 cgg)  
 #endif  
 #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do k = 1,Nr  
                     do j = jmin,jmax  
                       i = OB_Iw(j,bi,bj)  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  
      &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 #endif  
126    
127            do bj = jtlo,jthi          do bj = jtlo,jthi
128              do bi = itlo,ithi           do bi = itlo,ithi
129                do k = 1,nr  #ifdef ALLOW_OBCS_CONTROL_MODES
130                  do j = jmin,jmax            if (iobcs .gt. 2) then
131                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)             do j = jmin,jmax
132                   enddo              i = OB_Iw(j,bi,bj)
133                enddo  cih    Determine number of open vertical layers.
134                nz = 0
135                do k = 1,Nr
136                  if (iobcs .eq. 3) then
137                    nz = nz + maskW(i+ip1,j,k,bi,bj)
138                  else
139                    nz = nz + maskS(i,j,k,bi,bj)
140                  endif
141                end do
142    cih    Compute absolute velocities from the barotropic-baroclinic modes.
143                do k = 1,Nr
144                 if (k.le.nz) then
145                  stmp = 0.
146                  do nk = 1,nz
147                   stmp = stmp +
148         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
149                  end do
150                   tmpz(k,bi,bj) = stmp
151                 else
152                  tmpz(k,bi,bj) = 0.
153                 end if
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           endif
176    
177          if ( (obcswfirst) .or. (obcswchanged)) then         if ( (obcswfirst) .or. (obcswchanged)) then
178              
179  cgg(    This is a terribly long way to do it. However, the dimensions don't exactly          do bj = jtlo,jthi
180  cgg     match up. I will blame Fortran for the ugliness.           do bi = itlo,ithi
181              do k = 1,nr
182            do bj = jtlo,jthi             do j = jmin,jmax
183              do bi = itlo,ithi              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
184                do k = 1,nr              tmpfldyz (j,k,bi,bj)       = 0. _d 0
185                  do j = jmin,jmax             enddo
   
                   tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)  
                 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         &                       doglobalread, ladinit, optimcycle,
193         &                       mythid, xx_obcsw_dummy )
194    
195            do bj = jtlo,jthi          do bj = jtlo,jthi
196              do bi = itlo,ithi           do bi = itlo,ithi
197                do k = 1,nr  #ifdef ALLOW_OBCS_CONTROL_MODES
198                  do j = jmin,jmax            if (iobcs .gt. 2) then
199                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                     do j = jmin,jmax
200                  enddo              i = OB_Iw(j,bi,bj)
201                enddo  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 BALANCE_CONTROL_VOLFLUX                endif
231            if (optimcycle .gt. 0) then              end do
232              if (iobcs .eq. 3) then             enddo
 cgg(        Make sure that the xx velocity field has a balanced net volume flux.  
 cgg         Find the barotropic flow normal to the boundary.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(j,bi,bj)  
                     ubaro = 0.d0  
                     ubarocount = 0.d0  
                     do k = 1,nr  
 cgg(   If cells are not full, this should be modified with hFac.  
                       ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)  
      &                      * delZ(k) + ubaro  
                       ubarocount = maskW(i+ip1,j,k,bi,bj)  
      &                           * delZ(k) + ubarocount  
                     enddo  
                     if (ubarocount .eq. 0.) then  
                       ubaro = 0.  
                     else  
                       ubaro = ubaro / ubarocount  
                     endif  
                     do k = 1,nr  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 cgg)  
 #endif  
 #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do k = 1,Nr  
                     do j = jmin,jmax  
                       i = OB_Iw(j,bi,bj)  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  
      &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
233            endif            endif
234  #endif  #endif
235              do k = 1,nr
236            do bj = jtlo,jthi             do j = jmin,jmax
237              do bi = itlo,ithi              xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
238                do k = 1,nr  cgg   &                                        *   maskyz (j,k,bi,bj)
                 do j = jmin,jmax  
                   xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)  
                  enddo  
               enddo  
             enddo  
           enddo  
         endif  
   
 c--     Add control to model variable.  
         do bj = jtlo,jthi  
            do bi = itlo,ithi  
 c--        Calculate mask for tracer cells (0 => land, 1 => water).  
               do k = 1,nr  
                  do j = 1,sny  
                     i = OB_Iw(j,bi,bj)  
                     if (iobcs .EQ. 1) then  
                        OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 2) then  
                        OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 3) then  
                        OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 4) then  
                        OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)  
      &                      *maskS(i,j,k,bi,bj)  
                     endif  
                  enddo  
               enddo  
239             enddo             enddo
240              enddo
241             enddo
242          enddo          enddo
243           endif
244    
245    c--   Add control to model variable.
246           do bj = jtlo, jthi
247            do bi = itlo, ithi
248    c--   Calculate mask for tracer cells (0 => land, 1 => water).
249             do k = 1,nr
250              do j = 1,sny
251               i = OB_Iw(j,bi,bj)
252               if (iobcs .EQ. 1) then
253                OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
254         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
255         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
256                OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
257         &           *maskW(i+ip1,j,k,bi,bj)
258               else if (iobcs .EQ. 2) then
259                OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
260         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
261         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
262                OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
263         &           *maskW(i+ip1,j,k,bi,bj)
264               else if (iobcs .EQ. 3) then
265                OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
266         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
267         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
268                OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
269         &           *maskW(i+ip1,j,k,bi,bj)
270               else if (iobcs .EQ. 4) then
271                OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
272         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
273         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
274                OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
275         &           *maskS(i,j,k,bi,bj)
276               endif
277              enddo
278             enddo
279            enddo
280           enddo
281          
282  C--   End over iobcs loop  C--   End over iobcs loop
283        enddo        enddo
284    
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22