/[MITgcm]/MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcse.F
ViewVC logotype

Diff of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcse.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.5 by mmazloff, Wed Jan 19 21:08:32 2011 UTC revision 1.6 by mmazloff, Fri Mar 18 18:43:05 2011 UTC
# Line 66  c     == local variables == Line 66  c     == local variables ==
66        integer obcsecount0        integer obcsecount0
67        integer obcsecount1        integer obcsecount1
68        integer ip1        integer ip1
       integer nk,nz  
       _RL     tmpz (nr,nsx,nsy)  
       _RL     stmp  
69    
70  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
71          _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy)
72    
73        logical doglobalread        logical doglobalread
74        logical ladinit        logical ladinit
75    
76        character*(80) fnameobcse        character*(80) fnameobcse
77    
78  cgg(  Variables for splitting barotropic/baroclinic vels.  cmm( modes:
79        _RL ubaro        integer nk,nz
80        _RL utop        _RL     tmpz (nr,nsx,nsy)
81  cgg)        _RL     stmp
82    
83  c     == external functions ==  c     == external functions ==
84    
# Line 100  c     == end of interface == Line 98  c     == end of interface ==
98        imax = snx+olx        imax = snx+olx
99        ip1  = 0        ip1  = 0
100    
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       utop = 0.d0  
 cgg)  
   
101  c--   Now, read the control vector.  c--   Now, read the control vector.
102        doglobalread = .false.        doglobalread = .false.
103        ladinit      = .false.        ladinit      = .false.
104    
105        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
106          ilobcse=ilnblnk( xx_obcse_file )         ilobcse=ilnblnk( xx_obcse_file )
107          write(fnameobcse(1:80),'(2a,i10.10)')         write(fnameobcse(1:80),'(2a,i10.10)')
108       &       xx_obcse_file(1:ilobcse), '.', optimcycle       &      xx_obcse_file(1:ilobcse), '.', optimcycle
109        endif        endif
110    
111  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
# Line 124  c--   Get the counters, flags, and the i Line 117  c--   Get the counters, flags, and the i
117    
118        do iobcs = 1,nobcs        do iobcs = 1,nobcs
119    
120          if ( obcsefirst ) then         if ( obcsefirst ) then
121            call active_read_yz( fnameobcse, tmpfldyz,          call active_read_yz( fnameobcse, tmpfldyz,
122       &                         (obcsecount0-1)*nobcs+iobcs,       &                       (obcsecount0-1)*nobcs+iobcs,
123       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
124       &                         mythid, xx_obcse_dummy )       &                       mythid, xx_obcse_dummy )
125    
126            do bj = jtlo,jthi
127            if ( optimcycle .gt. 0) then           do bi = itlo,ithi
             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  
 cih    If open boundary.  
                      if ( OB_Ie(j,bi,bj) .ne. 0. ) then  
                     i = OB_Ie(J,bi,bj)  
 #ifdef ALLOW_OBCS_CONTROL_MODES  
 cih    Determine number of open vertical layers.  
                            nz = 0  
                            do k = 1,Nr  
                               nz = nz + maskW(i+ip1,j,k,bi,bj)  
                            end do  
 cih    Compute absolute velocities from the barotropic-baroclinic modes.  
 #ifdef ALLOW_CTRL_OBCS_BALANCE  
 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are  
 c     compatible - to ensure volume conservation can just set barotropic  
 c     mode amplitude to 0  
 c     however this means no inflow at every horizontal location....  
                    do k = 1,nr  
                      tmpfldyz(k,1,bi,bj)= 0.  
                    end do  
 #endif  
                            do k = 1,Nr  
                              if (k.le.nz) then  
                                 stmp = 0.  
                                 do nk = 1,nz  
                                  stmp = stmp +  
      &                            modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)  
                                 end do  
                                 tmpz(k,bi,bj) = stmp  
                              else  
                                 tmpz(k,bi,bj) = 0.  
                              end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacW(i+ip1,j,k,bi,bj)  
                            end do  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 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)  
 #endif  
 cih    End open boundary.  
                         end if  
                   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  
                      if ( OB_Ie(j,bi,bj) .ne. 0. ) then  
                     i = OB_Ie(J,bi,bj)  
128  #ifdef ALLOW_OBCS_CONTROL_MODES  #ifdef ALLOW_OBCS_CONTROL_MODES
129              if (iobcs .gt. 2) then
130               do j = jmin,jmax
131                i = OB_Ie(j,bi,bj)
132  cih    Determine number of open vertical layers.  cih    Determine number of open vertical layers.
133                             nz = 0              nz = 0
134                             do k = 1,Nr              do k = 1,Nr
135                                nz = nz + maskS(i,j,k,bi,bj)                if (iobcs .eq. 3) then
136                             end do                  nz = nz + maskW(i+ip1,j,k,bi,bj)
137  cih    Compute absolute velocities from the barotropic-baroclinic modes.                else
138  #ifdef ALLOW_CTRL_OBCS_BALANCE                  nz = nz + maskS(i,j,k,bi,bj)
139  CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are                endif
140  c     compatible - to ensure volume conservation can just set barotropic              end do
141  c     mode amplitude to 0  cih    Compute absolute velocities from the barotropic-baroclinic modes.
142  c     however this means no inflow at every horizontal location....              do k = 1,Nr
143                     do k = 1,nr               if (k.le.nz) then
144                       tmpfldyz(k,1,bi,bj)= 0.                stmp = 0.
145                     end do                do nk = 1,nz
146  #endif                 stmp = stmp +
147                             do k = 1,Nr       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
148                                if (k.le.nz) then                end do
149                                   stmp = 0.                 tmpz(k,bi,bj) = stmp
150                                   do nk = 1,nz               else
151                                    stmp = stmp +                tmpz(k,bi,bj) = 0.
152       &                             modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)               end if
                                  end do  
                                  tmpz(k,bi,bj) = stmp  
                               else  
                                  tmpz(k,bi,bj) = 0.  
                               end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacS(i,j,k,bi,bj)  
                            end do  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 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)  
 #endif  
 cih    End open boundary.  
                         end if  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcse1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
153              enddo              enddo
154                do k = 1,Nr
155                  if (iobcs .eq. 3) then
156                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
157         &            *recip_hFacW(i+ip1,j,k,bi,bj)
158                  else
159                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
160         &            *recip_hFacS(i,j,k,bi,bj)
161                  endif
162                end do
163               enddo
164              endif
165    #endif
166              do k = 1,nr
167               do j = jmin,jmax
168                xx_obcse1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
169    cgg   &                                        *   maskyz (j,k,bi,bj)
170               enddo
171            enddo            enddo
172          endif           enddo
173            enddo
174          if ( (obcsefirst) .or. (obcsechanged)) then         endif
175          
176            do bj = jtlo,jthi         if ( (obcsefirst) .or. (obcsechanged)) then
177             do bi = itlo,ithi          
178              do j = jmin,jmax          do bj = jtlo,jthi
179               do k = 1,nr           do bi = itlo,ithi
180                xx_obcse0(j,k,bi,bj,iobcs) = xx_obcse1(j,k,bi,bj,iobcs)            do j = jmin,jmax
181                tmpfldyz (j,k,bi,bj)       = 0. _d 0             do k = 1,nr
182               enddo              xx_obcse0(j,k,bi,bj,iobcs) = xx_obcse1(j,k,bi,bj,iobcs)
183              enddo              tmpfldyz (j,k,bi,bj)       = 0. _d 0
184             enddo             enddo
185            enddo            enddo
186             enddo
187            enddo
188            
189            call active_read_yz( fnameobcse, tmpfldyz,
190         &                       (obcsecount1-1)*nobcs+iobcs,
191         &                       doglobalread, ladinit, optimcycle,
192         &                       mythid, xx_obcse_dummy )
193    
194            call active_read_yz( fnameobcse, tmpfldyz,          do bj = jtlo,jthi
195       &                         (obcsecount1-1)*nobcs+iobcs,           do bi = itlo,ithi
      &                         doglobalread, ladinit, optimcycle,  
      &                         mythid, xx_obcse_dummy )  
   
           if ( optimcycle .gt. 0) then  
             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  
 cih    If open boundary.  
                      if ( OB_Ie(j,bi,bj) .ne. 0. ) then  
                     i = OB_Ie(J,bi,bj)  
 #ifdef ALLOW_OBCS_CONTROL_MODES  
 cih    Determine number of open vertical layers.  
                            nz = 0  
                            do k = 1,Nr  
                               nz = nz + maskW(i+ip1,j,k,bi,bj)  
                            end do  
 cih    Compute absolute velocities from the barotropic-baroclinic modes.  
 #ifdef ALLOW_CTRL_OBCS_BALANCE  
 CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are  
 c     compatible - to ensure volume conservation can just set barotropic  
 c     mode amplitude to 0  
 c     however this means no inflow at every horizontal location....  
                    do k = 1,nr  
                      tmpfldyz(k,1,bi,bj)= 0.  
                    end do  
 #endif  
                            do k = 1,Nr  
                              if (k.le.nz) then  
                                 stmp = 0.  
                                 do nk = 1,nz  
                                  stmp = stmp +  
      &                            modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)  
                                 end do  
                                 tmpz(k,bi,bj) = stmp  
                              else  
                                 tmpz(k,bi,bj) = 0.  
                              end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacW(i+ip1,j,k,bi,bj)  
                            end do  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 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)  
 #endif  
 cih    End open boundary.  
                         end if  
                   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  
                      if ( OB_Ie(j,bi,bj) .ne. 0. ) then  
                     i = OB_Ie(J,bi,bj)  
196  #ifdef ALLOW_OBCS_CONTROL_MODES  #ifdef ALLOW_OBCS_CONTROL_MODES
197              if (iobcs .gt. 2) then
198               do j = jmin,jmax
199                i = OB_Ie(j,bi,bj)
200  cih    Determine number of open vertical layers.  cih    Determine number of open vertical layers.
201                             nz = 0              nz = 0
202                             do k = 1,Nr              do k = 1,Nr
203                                nz = nz + maskS(i,j,k,bi,bj)                if (iobcs .eq. 3) then
204                             end do                  nz = nz + maskW(i+ip1,j,k,bi,bj)
205  cih    Compute absolute velocities from the barotropic-baroclinic modes.                else
206  #ifdef ALLOW_CTRL_OBCS_BALANCE                  nz = nz + maskS(i,j,k,bi,bj)
207  CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are                endif
208  c     compatible - to ensure volume conservation can just set barotropic              end do
209  c     mode amplitude to 0  cih    Compute absolute velocities from the barotropic-baroclinic modes.
210  c     however this means no inflow at every horizontal location....              do k = 1,Nr
211                     do k = 1,nr               if (k.le.nz) then
212                       tmpfldyz(k,1,bi,bj)= 0.                stmp = 0.
213                     end do                do nk = 1,nz
214  #endif                 stmp = stmp +
215                             do k = 1,Nr       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
216                                if (k.le.nz) then                end do
217                                   stmp = 0.                 tmpz(k,bi,bj) = stmp
218                                   do nk = 1,nz               else
219                                    stmp = stmp +                tmpz(k,bi,bj) = 0.
220       &                             modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)               endif
                                  end do  
                                  tmpz(k,bi,bj) = stmp  
                               else  
                                  tmpz(k,bi,bj) = 0.  
                               end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacS(i,j,k,bi,bj)  
                            end do  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 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)  
 #endif  
 cih    End open boundary.  
                         end if  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
221              enddo              enddo
222            enddo              do k = 1,Nr
223          endif                if (iobcs .eq. 3) then
224                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
225  c--     Add control to model variable.       &            *recip_hFacW(i+ip1,j,k,bi,bj)
226          do bj = jtlo,jthi                else
227             do bi = itlo,ithi                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
228  c--        Calculate mask for tracer cells (0 => land, 1 => water).       &            *recip_hFacS(i,j,k,bi,bj)
229                do k = 1,nr                endif
230                   do j = 1,sny              end do
                     i = OB_Ie(j,bi,bj)  
                     if (iobcs .EQ. 1) then  
                        OBEt(j,k,bi,bj) = OBEt (j,k,bi,bj)  
      &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)  
                        OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 2) then  
                        OBEs(j,k,bi,bj) = OBEs (j,k,bi,bj)  
      &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)  
                        OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 3) then  
                        OBEu(j,k,bi,bj) = OBEu (j,k,bi,bj)  
      &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)  
                        OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 4) then  
                        OBEv(j,k,bi,bj) = OBEv (j,k,bi,bj)  
      &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)  
                        OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)  
      &                      *maskS(i,j,k,bi,bj)  
                     endif  
                  enddo  
               enddo  
231             enddo             enddo
232              endif
233    #endif
234              do k = 1,nr
235               do j = jmin,jmax
236                xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
237    cgg   &                                        *   maskyz (j,k,bi,bj)
238               enddo
239              enddo
240             enddo
241          enddo          enddo
242           endif
243          
244    c--   Add control to model variable.
245           do bj = jtlo,jthi
246            do bi = itlo,ithi
247    c--   Calculate mask for tracer cells (0 => land, 1 => water).
248             do k = 1,nr
249              do j = 1,sny
250               i = OB_Ie(j,bi,bj)
251               if (iobcs .EQ. 1) then
252                OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
253         &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
254         &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
255                OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
256         &           *maskW(i+ip1,j,k,bi,bj)
257               else if (iobcs .EQ. 2) then
258                OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
259         &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
260         &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
261                OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
262         &           *maskW(i+ip1,j,k,bi,bj)
263               else if (iobcs .EQ. 3) then
264                OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
265         &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
266         &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
267                OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
268         &           *maskW(i+ip1,j,k,bi,bj)
269               else if (iobcs .EQ. 4) then
270                OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
271         &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
272         &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
273                OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
274         &           *maskS(i,j,k,bi,bj)
275               endif
276              enddo
277             enddo
278            enddo
279           enddo
280            
281  C--   End over iobcs loop  C--   End over iobcs loop
282        enddo        enddo
283    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22