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

Diff of /MITgcm_contrib/SOSE/BoxAdj/code_ad/ctrl_getobcsn.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 67  c     == local variables == Line 67  c     == local variables ==
67        logical obcsnchanged        logical obcsnchanged
68        integer obcsncount0        integer obcsncount0
69        integer obcsncount1        integer obcsncount1
       integer nk,nz  
       _RL     tmpz (nr,nsx,nsy)  
       _RL     stmp  
70    
71  cgg      _RL maskxz   (1-olx:snx+olx,nr,nsx,nsy)  cgg      _RL maskxz   (1-olx:snx+olx,nr,nsx,nsy)
72          _RL tmpfldxz (1-olx:snx+olx,nr,nsx,nsy)
73    
74        logical doglobalread        logical doglobalread
75        logical ladinit        logical ladinit
76    
77        character*(80) fnameobcsn        character*(80) fnameobcsn
78    
79  cgg(  Variables for splitting barotropic/baroclinic vels.  cmm( modes:
80        _RL vbaro        integer nk,nz
81        _RL vtop        _RL     tmpz (nr,nsx,nsy)
82          _RL     stmp
83    
84  c     == external functions ==  c     == external functions ==
85    
# Line 105  cgg      imax = snx+olx Line 104  cgg      imax = snx+olx
104        imax = snx        imax = snx
105        jp1  = 0        jp1  = 0
106    
 cgg(  Initialize variables for balancing volume flux.  
       vbaro = 0.d0  
       vtop = 0.d0  
 cgg)  
   
107  c--   Now, read the control vector.  c--   Now, read the control vector.
108        doglobalread = .false.        doglobalread = .false.
109        ladinit      = .false.        ladinit      = .false.
# Line 128  c--   Get the counters, flags, and the i Line 122  c--   Get the counters, flags, and the i
122       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
123    
124        do iobcs = 1,nobcs        do iobcs = 1,nobcs
125          if ( obcsnfirst ) then         if ( obcsnfirst ) then
126            call active_read_xz( fnameobcsn, tmpfldxz,          call active_read_xz( fnameobcsn, tmpfldxz,
127       &                         (obcsncount0-1)*nobcs+iobcs,       &                       (obcsncount0-1)*nobcs+iobcs,
128       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
129       &                         mythid, xx_obcsn_dummy )       &                       mythid, xx_obcsn_dummy )
130    
131            if ( optimcycle .gt. 0) then          do bj = jtlo,jthi
132  c  normal velocity may require special treatment           do bi = itlo,ithi
             if (iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do i = imin,imax  
 cih    If open boundary.  
                         if ( OB_Jn(i,bi,bj) .ne. 0. ) then  
                            j = OB_Jn(i,bi,bj)  
133  #ifdef ALLOW_OBCS_CONTROL_MODES  #ifdef ALLOW_OBCS_CONTROL_MODES
134  c  ALLOW_OBCS_CONTROL_MODES means controls are in modes            if (iobcs .gt. 2) then
135               do i = imin,imax
136                j = OB_Jn(i,bi,bj)
137  cih    Determine number of open vertical layers.  cih    Determine number of open vertical layers.
138                             nz = 0              nz = 0
139                             do k = 1,Nr              do k = 1,Nr
140                                nz = nz + maskS(i,j+jp1,k,bi,bj)               if (iobcs .eq. 3) then
141                             end do                nz = nz + maskS(i,j+jp1,k,bi,bj)
142                 else
143                  nz = nz + maskW(i,j,k,bi,bj)
144                 endif
145                end do
146  cih    Compute absolute velocities from the barotropic-baroclinic modes.  cih    Compute absolute velocities from the barotropic-baroclinic modes.
147  #ifdef ALLOW_CTRL_OBCS_BALANCE              do k = 1,Nr
148  CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are               if (k.le.nz) then
149  c     compatible - to ensure volume conservation can just set barotropic                stmp = 0.
150  c     mode amplitude to 0                do nk = 1,nz
151  c     however this means no inflow at every horizontal location....                 stmp = stmp +
152                     do k = 1,nr       &         modesv(k,nk,nz)*tmpfldxz(i,nk,bi,bj)
153                       tmpfldxz(k,1,bi,bj)= 0.                end do
154                     end do                tmpz(k,bi,bj) = stmp
155  #endif               else
156                             do k = 1,Nr                tmpz(k,bi,bj) = 0.
157                                if (k.le.nz) then               end if
158                                   stmp = 0.              end do
159                                   do nk = 1,nz              do k = 1,Nr
160                                     stmp = stmp +               if (iobcs .eq. 3) then
161       &                              modesv(k,nk,nz)*tmpfldxz(i,nk,bi,bj)                tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)
162                                   end do       &         *recip_hFacS(i,j+jp1,k,bi,bj)
163                                   tmpz(k,bi,bj) = stmp               else
164                                else                tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)
                                  tmpz(k,bi,bj) = 0.  
                               end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacS(i,j+jp1,k,bi,bj)  
                            end do  
 c   if not controlling modes may need to balance obcs...  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 cgg         The barotropic velocity is stored in the level 1.  
                     vbaro = tmpfldxz(i,1,bi,bj)  
 cgg         Except for the special point which balances barotropic vol.flux.  
 cgg         Special column in the NW corner.  
                     if (ob_iw(j,bi,bj).eq.(i-1).and.  
      &                  ob_iw(j,bi,bj).ne. 0) then  
                         print*,'Apply shiftvel1 @ i,j'  
                         print*,shiftvel(1),i,j  
                       vbaro = shiftvel(1)  
                     endif  
                     tmpfldxz(i,1,bi,bj) = 0.d0  
                     vtop = 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..  
                       vtop = tmpfldxz(i,k,bi,bj)*  
      &                maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then  
                         tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro  
                       endif  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)  
      &                                      - vtop / delR(1)  
 #endif  
 cih    End if open boundary.  
                         end if  
 cih    End loop over x-points  
                      end do  
                   end do  
                end do  
 cih    End if iobcs = 3.  
             end if  
   
 cih    If tangential velocity    
             if (iobcs .eq. 4) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
 cih    Begin loop over x-points.  
                   do i = imin,imax  
 cih    If open boundary.  
                        if ( OB_Jn(i,bi,bj) .ne. 0. ) then  
                           j = OB_Jn(i,bi,bj)  
 #ifdef ALLOW_OBCS_CONTROL_MODES  
 cih    Determine number of open vertical layers.  
                           nz = 0  
                           do k = 1,Nr  
                              nz = nz + maskW(i,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  
                      tmpfldxz(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)*tmpfldxz(i,nk,bi,bj)  
                                 end do  
                                 tmpz(k,bi,bj) = stmp  
                              else  
                                 tmpz(k,bi,bj) = 0.  
                              end if  
                           end do  
                            do k = 1,Nr  
                         tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)  
165       &                           *recip_hFacW(i,j,k,bi,bj)       &                           *recip_hFacW(i,j,k,bi,bj)
166                             end do               endif
167  #elif defined (ALLOW_CTRL_OBCS_BALANCE)              end do
168  cgg         The barotropic velocity is stored in the level 1.             enddo
169                      vbaro = tmpfldxz(i,1,bi,bj)            endif
170  cgg         Except for the special point which balances barotropic vol.flux.  #endif
171  cgg         Special column in the NW corner.            do k = 1,nr
172                      tmpfldxz(i,1,bi,bj) = 0.d0             do i = imin,imax
173                      vtop = 0.d0              xx_obcsn1(i,k,bi,bj,iobcs)  = tmpfldxz (i,k,bi,bj)
174    cgg   &                                    *   maskxz (i,k,bi,bj)
175                      do k = 1,Nr             enddo
 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..  
                       vtop = tmpfldxz(i,k,bi,bj)*  
      &                maskW(i,j,k,bi,bj) * delR(k) + vtop  
 cgg    Add the barotropic velocity component.  
                       if (maskW(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)  
      &                                      - vtop / delR(1)  
 #endif /* ALLOW_OBCS_CONTROL_MODES */  
 cih    End open boundary.  
                        end if  
 cih    End loop over x-points  
                    end do  
                 end do  
              end do  
 cih    End if iobcs = 4.  
           end if  
 cih    End if optimcycle > 0 .    
          end if  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do i = imin,imax  
                   xx_obcsn1(i,k,bi,bj,iobcs)  = tmpfldxz (i,k,bi,bj)  
 cgg     &                                        *   maskxz (i,k,bi,bj)  
                 enddo  
               enddo  
             enddo  
176            enddo            enddo
177          endif           enddo
178            enddo
179          if ( (obcsnfirst) .or. (obcsnchanged)) then         endif
180    
181            do bj = jtlo,jthi         if ( (obcsnfirst) .or. (obcsnchanged)) then
182             do bi = itlo,ithi          
183              do k = 1,nr          do bj = jtlo,jthi
184               do i = imin,imax           do bi = itlo,ithi
185                xx_obcsn0(i,k,bi,bj,iobcs) = xx_obcsn1(i,k,bi,bj,iobcs)            do k = 1,nr
186                tmpfldxz (i,k,bi,bj)       = 0. _d 0             do i = imin,imax
187               enddo              xx_obcsn0(i,k,bi,bj,iobcs) = xx_obcsn1(i,k,bi,bj,iobcs)
188              enddo              tmpfldxz (i,k,bi,bj)       = 0. _d 0
189             enddo             enddo
190            enddo            enddo
191             enddo
192            enddo
193            
194            call active_read_xz( fnameobcsn, tmpfldxz,
195         &                       (obcsncount1-1)*nobcs+iobcs,
196         &                       doglobalread, ladinit, optimcycle,
197         &                       mythid, xx_obcsn_dummy )
198    
199            call active_read_xz( fnameobcsn, tmpfldxz,          do bj = jtlo,jthi
200       &                         (obcsncount1-1)*nobcs+iobcs,           do bi = itlo,ithi
      &                         doglobalread, ladinit, optimcycle,  
      &                         mythid, xx_obcsn_dummy )  
   
           if ( optimcycle .gt. 0) then  
 c  normal velocity may require special treatment  
             if (iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do i = imin,imax  
 cih    If open boundary.  
                         if ( OB_Jn(i,bi,bj) .ne. 0. ) then  
                            j = OB_Jn(i,bi,bj)  
201  #ifdef ALLOW_OBCS_CONTROL_MODES  #ifdef ALLOW_OBCS_CONTROL_MODES
202  c  ALLOW_OBCS_CONTROL_MODES means controls are in modes            if (iobcs .gt. 2) then
203               do i = imin,imax
204                j = OB_Jn(i,bi,bj)
205  cih    Determine number of open vertical layers.  cih    Determine number of open vertical layers.
206                             nz = 0              nz = 0
207                             do k = 1,Nr              do k = 1,Nr
208                                nz = nz + maskS(i,j+jp1,k,bi,bj)               if (iobcs .eq. 3) then
209                             end do                nz = nz + maskS(i,j+jp1,k,bi,bj)
210                 else
211                  nz = nz + maskW(i,j,k,bi,bj)
212                 endif
213                end do
214  cih    Compute absolute velocities from the barotropic-baroclinic modes.  cih    Compute absolute velocities from the barotropic-baroclinic modes.
215  #ifdef ALLOW_CTRL_OBCS_BALANCE              do k = 1,Nr
216  CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are               if (k.le.nz) then
217  c     compatible - to ensure volume conservation can just set barotropic                stmp = 0.
218  c     mode amplitude to 0                do nk = 1,nz
219  c     however this means no inflow at every horizontal location....                 stmp = stmp +
220                     do k = 1,nr       &         modesv(k,nk,nz)*tmpfldxz(i,nk,bi,bj)
221                       tmpfldxz(k,1,bi,bj)= 0.                end do
222                     end do                tmpz(k,bi,bj) = stmp
223  #endif               else
224                             do k = 1,Nr                tmpz(k,bi,bj) = 0.
225                                if (k.le.nz) then               end if
226                                   stmp = 0.              end do
227                                   do nk = 1,nz              do k = 1,Nr
228                                     stmp = stmp +               if (iobcs .eq. 3) then
229       &                              modesv(k,nk,nz)*tmpfldxz(i,nk,bi,bj)                tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)
230                                   end do       &         *recip_hFacS(i,j+jp1,k,bi,bj)
231                                   tmpz(k,bi,bj) = stmp               else
232                                else                tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)
                                  tmpz(k,bi,bj) = 0.  
                               end if  
                            end do  
                            do k = 1,Nr  
                         tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)  
      &                           *recip_hFacS(i,j+jp1,k,bi,bj)  
                            end do  
 c   if not controlling modes may need to balance obcs...  
 #elif defined (ALLOW_CTRL_OBCS_BALANCE)  
 cgg         The barotropic velocity is stored in the level 1.  
                     vbaro = tmpfldxz(i,1,bi,bj)  
 cgg         Except for the special point which balances barotropic vol.flux.  
 cgg         Special column in the NW corner.  
                     if (ob_iw(j,bi,bj).eq.(i-1).and.  
      &                  ob_iw(j,bi,bj).ne. 0) then  
                         print*,'Apply shiftvel1 @ i,j'  
                         print*,shiftvel(1),i,j  
                       vbaro = shiftvel(1)  
                     endif  
                     tmpfldxz(i,1,bi,bj) = 0.d0  
                     vtop = 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..  
                       vtop = tmpfldxz(i,k,bi,bj)*  
      &                maskS(i,j+jp1,k,bi,bj) * delR(k) + vtop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then  
                         tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro  
                       endif  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)  
      &                                      - vtop / delR(1)  
 #endif  
 cih    End if open boundary.  
                         end if  
 cih    End loop over x-points  
                      end do  
                   end do  
                end do  
 cih    End if iobcs = 3.  
             end if  
   
 cih    If tangential velocity    
             if (iobcs .eq. 4) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
 cih    Begin loop over x-points.  
                   do i = imin,imax  
 cih    If open boundary.  
                        if ( OB_Jn(i,bi,bj) .ne. 0. ) then  
                           j = OB_Jn(i,bi,bj)  
 #ifdef ALLOW_OBCS_CONTROL_MODES  
 cih    Determine number of open vertical layers.  
                           nz = 0  
                           do k = 1,Nr  
                              nz = nz + maskW(i,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  
                      tmpfldxz(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)*tmpfldxz(i,nk,bi,bj)  
                                 end do  
                                 tmpz(k,bi,bj) = stmp  
                              else  
                                 tmpz(k,bi,bj) = 0.  
                              end if  
                           end do  
                            do k = 1,Nr  
                         tmpfldxz(i,k,bi,bj) = tmpz(k,bi,bj)  
233       &                           *recip_hFacW(i,j,k,bi,bj)       &                           *recip_hFacW(i,j,k,bi,bj)
234                             end do               endif
235  #elif defined (ALLOW_CTRL_OBCS_BALANCE)              end do
 cgg         The barotropic velocity is stored in the level 1.  
                     vbaro = tmpfldxz(i,1,bi,bj)  
 cgg         Except for the special point which balances barotropic vol.flux.  
 cgg         Special column in the NW corner.  
                     tmpfldxz(i,1,bi,bj) = 0.d0  
                     vtop = 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..  
                       vtop = tmpfldxz(i,k,bi,bj)*  
      &                maskW(i,j,k,bi,bj) * delR(k) + vtop  
 cgg    Add the barotropic velocity component.  
                       if (maskW(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldxz(i,k,bi,bj) = tmpfldxz(i,k,bi,bj)+ vbaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldxz(i,1,bi,bj) = tmpfldxz(i,1,bi,bj)  
      &                                      - vtop / delR(1)  
 #endif /* ALLOW_OBCS_CONTROL_MODES */  
 cih    End open boundary.  
                        end if  
 cih    End loop over x-points  
                    end do  
                 end do  
              end do  
 cih    End if iobcs = 4.  
           end if  
 cih    End if optimcycle > 0 .    
          end if  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do i = imin,imax  
                   xx_obcsn1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)  
 cgg     &                                        *   maskxz (i,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 i = 1,snx  
                     j = OB_Jn(I,bi,bj)  
                     if (iobcs .EQ. 1) then  
                        OBNt(i,k,bi,bj) = OBNt (i,k,bi,bj)  
      &                 + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)  
                        OBNt(i,k,bi,bj) = OBNt(i,k,bi,bj)  
      &                         *maskS(i,j+jp1,k,bi,bj)  
 cgg     &                      *maskxz(i,k,bi,bj)  
                     else if (iobcs .EQ. 2) then  
                        OBNs(i,k,bi,bj) = OBNs (i,k,bi,bj)  
      &                 + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)  
                        OBNs(i,k,bi,bj) = OBNs(i,k,bi,bj)  
      &                         *maskS(i,j+jp1,k,bi,bj)  
 cgg     &                      *maskxz(i,k,bi,bj)  
                     else if (iobcs .EQ. 4) then  
                        OBNu(i,k,bi,bj) = OBNu (i,k,bi,bj)  
      &                 + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)  
                        OBNu(i,k,bi,bj) = OBNu(i,k,bi,bj)  
      &                         *maskW(i,j,k,bi,bj)  
 cgg     &                      *maskxz(i,k,bi,bj)  
                     else if (iobcs .EQ. 3) then  
                        OBNv(i,k,bi,bj) = OBNv (i,k,bi,bj)  
      &                 + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)  
                        OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)  
      &                         *maskS(i,j+jp1,k,bi,bj)  
 cgg     &                      *maskxz(i,k,bi,bj)  
                     endif  
                  enddo  
               enddo  
236             enddo             enddo
237              endif
238    #endif
239              do k = 1,nr
240               do i = imin,imax
241                xx_obcsn1 (i,k,bi,bj,iobcs) = tmpfldxz (i,k,bi,bj)
242    cgg   &                                        *   maskxz (i,k,bi,bj)
243               enddo
244              enddo
245             enddo
246          enddo          enddo
247              
248           endif
249    
250    c--   Add control to model variable.
251           do bj = jtlo,jthi
252            do bi = itlo,ithi
253    c--   Calculate mask for tracer cells (0 => land, 1 => water).
254             do k = 1,nr
255              do i = 1,snx
256               j = OB_Jn(I,bi,bj)
257               if (iobcs .EQ. 1) then
258                OBNt(i,k,bi,bj) = OBNt (i,k,bi,bj)
259         &           + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)
260         &           + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
261                OBNt(i,k,bi,bj) = OBNt(i,k,bi,bj)
262         &           *maskS(i,j+jp1,k,bi,bj)
263               else if (iobcs .EQ. 2) then
264                OBNs(i,k,bi,bj) = OBNs (i,k,bi,bj)
265         &           + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)
266         &           + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
267                OBNs(i,k,bi,bj) = OBNs(i,k,bi,bj)
268         &           *maskS(i,j+jp1,k,bi,bj)
269               else if (iobcs .EQ. 4) then
270                OBNu(i,k,bi,bj) = OBNu (i,k,bi,bj)
271         &           + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)
272         &           + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
273                OBNu(i,k,bi,bj) = OBNu(i,k,bi,bj)
274         &           *maskW(i,j,k,bi,bj)
275               else if (iobcs .EQ. 3) then
276                OBNv(i,k,bi,bj) = OBNv (i,k,bi,bj)
277         &           + obcsnfac            *xx_obcsn0(i,k,bi,bj,iobcs)
278         &           + (1. _d 0 - obcsnfac)*xx_obcsn1(i,k,bi,bj,iobcs)
279                OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
280         &           *maskS(i,j+jp1,k,bi,bj)
281               endif
282              enddo
283             enddo
284            enddo
285           enddo
286          
287  C--   End over iobcs loop  C--   End over iobcs loop
288        enddo        enddo
289    
290  #else /* ALLOW_OBCSS_CONTROL undefined */  #else /* ALLOW_OBCSN_CONTROL undefined */
291    
292  c     == routine arguments ==  c     == routine arguments ==
293    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22