subroutine in_out_boundaries(xin,yin,nii,nji,
     $     xout,yout,nio,njo,cyclicxi,
     $     niip1,njip1,niop1,njop1,
     $     iunit_diag,
     $     gxout,gyout)
C         
C  Purpose:
C
C         calculate the location of grid box boundaries of
C         an "output" grid w.r.t. an "input" grid
C         
C         used in a 2-D regriding process, i.e.,
C         input --> output
C         
C  Input variables:
C         
C         xin - longitudes of the input grid
C         yin - latitudes of the input grid
C         xout - longtitudes of the output grid
C         yout - latitudes of the output grid
C         
C         cyclicxi - flag whether the input grid is cyclically continuous
C         in x
C         
C         nii - size of the x dimension of the input grid
C         nji - size of the y dimension of the input grid
C         nio - size of the x dimension of the output grid
C         njo - size of the y dimension of the output grid
C
C         niip1= nii+1, etc.
C
C         iunit_diag - unit number to write diagnositics
C
C  Output variables:         
C
C         gxout - x location of the output grid in input grid units
C         gyout - y location of the output grid in input grid units
C         

      dimension xin(*),yin(*),xout(*),yout(*),
     $     gxout(*),gyout(*)

      logical cyclicxi

C         
C         locate the output grid w.r.t. input grid in x
C         ALLOW FOR CYCLIC CONTINUITY IN X!!
C         
      ibeg=1
      iend=niip1
      
      do i=1,niop1
        
        x0=xout(i)
        
        do ii=ibeg,iend
C         
C         point is before start of input grid 
C         
          x1=xin(ii)
          x2=xin(niop1)

          if(x0.lt.x1) then
            if(cyclicxi) then
              ic=nii
              do while(x0.lt.x1)
                x2=x1
                x1=xin(ic)-360.0
                ic=ic-1
              enddo
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ic)+dx0
            else
              gxout(i)=0.0
            endif
            
            ibeg=1
            go to 50
            
C         
C         point is within the input grid
C         
          else if(ii.ge.1.and.ii.le.nii) then
            
            x2=xin(ii+1)
            
cccc            write(iunit_diag,*) ibeg,i,ii,x0,x1,x2
            
            if(x0.ge.x1.and.x0.lt.x2) then
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ii-1)+dx0
              ibeg=ii
              go to 50
            endif
           
C         
C         point is beyond input grid
C         
          else if(x0.ge.x2) then
            ic=2
            if(cyclicxi) then
              x1=xin(niip1)
              x2=xin(ic)+360.0

cccc              write(iunit_diag,*) 'after',i,ic,x0,x1,x2
              do while(x0.lt.x1.or.x0.ge.x2)
                ic=ic+1
                x1=x2
                x2=xin(ic)+360.0
              end do
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ic-2)+dx0
            else
              gxout(i)=float(nii)
            endif
            
            ibeg=niip1
            go to 50
            
          endif
          
        end do
        
 50     continue
        
      end do
C         
C         locate the output grid w.r.t. input grid in y
C         NO CYCLIC CONTINUITY!!
C         
      jbeg=1
      jend=njip1
      
      do j=1,njop1
        
        y0=yout(j)
        
        do jj=jbeg,jend
C         
C         point is before start of input grid 
C         
          y1=yin(jj)
          
          if(y0.lt.y1) then

            gyout(j)=0.0
            jbeg=1
cccc            write(iunit_diag,*) 'before y ',j,jj,gyout(j)
            go to 60
C         
C         point is within the input grid
C         
          else if(jj.ge.1.and.jj.le.nji) then
            
            y2=yin(jj+1)
            
cccc            write(iunit_diag,*) 'jbeg... ',jbeg,j,jj,y0,y1,y2
            
            if(y0.ge.y1.and.y0.lt.y2) then
              dy0=(y0-y1)/(y2-y1)
              gyout(j)=float(jj-1)+dy0
              jbeg=jj
              go to 60
            endif
C         
C         point is beyond input grid
C         
          else if(y0.ge.y2) then
            gyout(j)=float(nji)
            jbeg=njip1
cccc            write(iunit_diag,*) 'after y ',j,jj,gyout(j)
            go to 60
          endif
          
        end do
        
 60     continue
        
      end do

      return
      end