program regrid
ccc      parameter(nimax=642,njmax=320,n1max=nimax*2)
ccc      parameter(nimax=751,njmax=801,n1max=nimax*2)
      parameter(nimax=5000,njmax=5000,n1max=nimax*2)

      real vals(20),ovals(20),args(8)

      dimension xin(n1max),yin(n1max),xout(n1max),yout(n1max)
      dimension gxout(n1max),gyout(n1max)

      dimension dumx(n1max)

      dimension incard(80)

      dimension fld_in(nimax,njmax),area_in(nimax,njmax),
     $     fld_out(nimax,njmax),dum1(nimax,njmax)

      character card*80,qtitle*24,regrid_method*40
      character incard*1

      logical cyclicxi,wrapxi,
     $     diag_out,vote,bessel,
     $     npole_point,spole_point

      equivalence (incard(1),card)
C         
C         gfi - input file to the program with the GrADS data object
C         
C         gfo - output file to be read by GrADS
C
      iunit_in=8
      iunit_out=10
      iunit_diag=12
     
      open (iunit_in,file='udf.regrid.gfi',status='unknown',
     $     form='unformatted',err=810)

      open (iunit_out,file='udf.regrid.gfo',status='unknown',
     $     form='unformatted')

      open (iunit_diag,file='udf.regrid.out',
     $     form='formatted',status='unknown')
C         
C         constants
C
      pi=3.141592654
      deg2rad=pi/180.0
      rad2deg=180.0/pi
C         
C         defaults
C
      diag_out=.false.
cc      diag_out=.true.
      
C         
C         user defined beginning lat/lon
C
      rlatbeg_set=-999.0
      rlonbeg_set=-999.0
C         
C         pole points
C         
      spole_point=.false.
      spole_point=.false.

C         
C         the first record contains number of args
C         and other parameters for future implementations
C

      read (iunit_in) vals
      nargs=nint(vals(1))
      write(iunit_diag,*) 'number of arguments = ',nargs
C         
C         stop if only one argument
C         
      if(nargs.le.1) then
        write(*,12) 
 12     format(' ',/
     $       /' ','regrid requires 2 or more arguements:',
     $       /' ',' ')
        go to 999
      endif
C         
C         the second record is the field ("expr" in the user defined
C         function table) and has is input using four reads
C         
C         grid record 1
C
      read(iunit_in,err=820) vals
C         
C         grid parameters
C
      undef=vals(1)
      itype=nint(vals(2))
      jtype=nint(vals(3))
      nii=nint(vals(4))
      nji=nint(vals(5))
      ilineari=nint(vals(6))
      jlineari=nint(vals(7))
      xbeg=vals(8)
      
      write(*,*) 'itype=',itype
      write(*,*) 'jtype=',jtype
      write(*,*) 'nii=',nii
      write(*,*) 'nji=',nji
      write(*,*) 'ilineari=',ilineari
      write(*,*) 'jlineari=',jlineari
      write(*,*) 'xbeg=',xbeg
      
C         
C         check the input grid type
C         
C         abort if not 2-d x,y
C         
      if(itype.ne.0.or.jtype.ne.1) go to 830
C         
C         get grid increments if uniform
C         
      if(ilineari.eq.1) then
        dx=vals(9)
      else
        dx=-999
      endif
      ybeg=vals(10)

      if(jlineari.eq.1) then
        dy=vals(11)
      else
        dy=-999
      endif
C         
C         only get time values if a one of the dimensions is time
C     
      if(itype.eq.3.or.jtype.eq.3) then
        yrbeg=vals(12)
        mobeg=vals(13)
        dabeg=vals(14)
        hrbeg=vals(15)
        mnbeg=vals(16)
        dmn=vals(17)
        dmo=vals(18)
      endif
C         
      write(iunit_diag,102) undef,itype,jtype,nii,nji,
     $     ilineari,jlineari,xbeg,dx,ybeg,dy
 102  format(' ',
     $     /,10x,'          undefined value = ',e13.5,
     $     /,10x,'         x dimension type = ',i3,
     $     /,10x,'         y dimension type = ',i3,
     $     /,10x,'         x dimension size = ',i5,
     $     /,10x,'         y dimension size = ',i5,
     $     /,10x,' uniform grid spacing in x ? ',i2,
     $     /,10x,' uniform grid spacing in y ? ',i2,
     $     /,10x,'   x world coor beginning  = ',g12.4,
     $     /,10x,'                        dx = ',g12.4,
     $     /,10x,'   y world coor beginning  = ',g12.4,
     $     /,10x,'                        dy = ',g12.4)


      if(itype.eq.3.or.jtype.eq.3) then
        write(iunit_diag,104) 
     $       yrbeg,mobeg,dabeg,hrbeg,mnbeg,dmn,dmo
 104    format(' ',
     $       /,10x,'                year start = ',g12.4,
     $       /,10x,'               month start = ',g12.4,
     $       /,10x,'                 day start = ',g12.4,
     $       /,10x,'                hour start = ',g12.4,
     $       /,10x,'              minute start = ',g12.4,
     $       /,10x,'           minute increment = ',g12.4,
     $       /,10x,'            month increment = ',g12.4)
      endif

C
C         grid record #2 the data
C         
      call read_grid(iunit_in,fld_in,nii,nji,istat)
      if(istat.ne.1) go to 820

      qtitle='input field from GrADS '
      call qprntu(fld_in,qtitle,undef,1,1,nii,nji,4,iunit_diag)
C         
C         grid record #3 and 4 -- grid-to-world coordinate maps
C         of the input grid
C         
      read(iunit_in,err=820) (xin(i),i=1,nii)
      write(iunit_diag,*) (xin(i),i=1,nii)
      read(iunit_in,err=820) (yin(j),j=1,nji)
      write(iunit_diag,*) (yin(j),j=1,nji)
C         
C------------------------------------------------------------
C         
C         input grid parameters
C         check for cyclic continuity in x 
C         
C------------------------------------------------------------
C
      cyclicxi=.false.
      wrapxi=.false.
C
C         make sure the x dimension is longitude and is uniform
C         
      if(itype.eq.0.and.ilineari.eq.1) then
        xlen=(nii-1)*dx
        if(xlen.ge.360.0) then
          wrapxi=.true.
          cyclicxi=.true.
          niifix=360.0/dx
          xlen=360.0
        endif

        xlen=nii*dx
        if(xlen.eq.360.0) then
          wrapxi=.false.
          cyclicxi=.true.
        endif

      endif

      write(iunit_diag,*) 'cyclicxi: ',cyclicxi,wrapxi
C         
C------------------------------------------------------------
C
C         if wrapped in x, then trim the grid
C
C         dump the input data into a dummy array
C         and then load the trimmed grid into the original field array
C
C------------------------------------------------------------
C         
      if(wrapxi) then
        call vload2(dum1,fld_in,nii,nji)
        call trim_grid(dum1,nii,niifix,nji,fld_in)
C         
C         set the input size to the new trimmed x dimension
C         
        nii=niifix
      endif
C         
C------------------------------------------------------------
C         
C         set up the input grid
C         
C------------------------------------------------------------
C
C         boundaries of input grid boxes
C         
      rlonbegi=xin(1)
      rlonendi=xin(nii)
      if(cyclicxi) rlonedi=xin(1)+360.0

      rlatbegi=yin(1)
      rlatendi=yin(nji)

      niip1=nii+1
      njip1=nji+1
C         
C         x in
C         
      do i=2,nii
        dumx(i)=0.5*(xin(i-1)+xin(i))
      end do

      dumx(1)=xin(1)-0.5*(xin(2)-xin(1))
      dumx(niip1)=xin(nii)+0.5*(xin(nii)-xin(nii-1))

      do i=1,niip1
        xin(i)=dumx(i)
      end do
C         
C         y in
C
      do j=2,nji
        dumx(j)=0.5*(yin(j-1)+yin(j))
      end do

      dumx(1)=yin(1)-0.5*(yin(2)-yin(1))
      dumx(njip1)=yin(nji)+0.5*(yin(nji)-yin(nji-1))
C         
C         make sure the input grid does not extend beyond the poles
C
      do j=1,njip1
        yin(j)=dumx(j)
        if(yin(j).lt.-90.0) yin(j)=-90.0
        if(yin(j).gt.90.0) yin(j)=90.0
      end do

      write(iunit_diag,*) 'input grid x boundaries'
      write(iunit_diag,*) (xin(i),i=1,niip1)
      write(iunit_diag,*) 'input grid y boundaries'
      write(iunit_diag,*) (yin(j),j=1,njip1)
C         
C------------------------------------------------------------
C
C         read arguments after the grid
C
C------------------------------------------------------------
C         
      nargs=nargs-1
      do i=1,nargs
        read (iunit_in,err=832) args(i)
      end do
C         
C         if one argument, then assume uniform lat lon with equal dx
C         
      if(nargs.eq.1) then
        iout_grid_type=1
        iregrid_method=1
        dxout=args(1)
        dyout=args(1)
C         
C         if two arguments, dx=dy and the second selects regrid method
C         
      else if(nargs.eq.2) then
        iout_grid_type=1
        iregrid_method=1
        dxout=args(1)
        dyout=dxout
        iregrid_type=nint(args(2))
        
C         
C         method
C         
        bessel=.false.
        vote=.false.
C         
C         check for valid type
C         
        if(iregrid_type.lt.1.or.iregrid_type.ge.5) go to 836
        
        if(iregrid_type.eq.1.or.iregrid_type.eq.11) then
          iregrid_method=1
        else if(iregrid_type.eq.2.or.iregrid_type.eq.12) then
          iregrid_method=1
          vote=.true.
        else if(iregrid_type.eq.3.or.iregrid_type.eq.13) then
          iregrid_method=2
        else if(iregrid_type.eq.4.or.iregrid_type.eq.14) then
          iregrid_method=2
          bessel=.true.
        end if
C         
C         three or more arguments
C         
      else if(nargs.ge.3) then
C         
C         the third argument indicates the type of interpolation
C         and output grid
C         
C         1 - box ave - uniform lat/lon
C         2 - vote box ave ; uniform lat/lon
C         3 - bilinear interpolation ; uniform lat/lon
C         4 - bessel interpolation ; uniform lat/lon
C         
C         11 - box ave ; gaussian grid
C         12 - vote box ave ; gaussian
C         13 - bilinear interp ; gaussian
C         14 - bessel interp ; gaussian
C         
        iregrid_type=nint(args(3))
C         
C         check for valid type
C         
        if(iregrid_type.lt.1.or.iregrid_type.gt.15.or.
     $       (iregrid_type.ge.5.and.iregrid_type.le.10)) go to 836
        
        
        if(iregrid_type.le.10) then
          
          dxout=args(1)
          dyout=args(2)
          ilinearo=1
          jlinearo=1
          iout_grid_type=1
          
        else if(iregrid_type.le.20) then
          
          niog=nint(args(1))
          njog=nint(args(2))
          nwaves=niog/3-1
          
C         
C         calculate the gaussian latitudes and 
C         the grid spacing for the gaussian longitudes
C         
          call gauss_lat_nmc(dumx,njog)
          
          ilinearo=0
          jlinearo=0
          dxout=360.0/float(niog)
          dyoutg=180.0/float(njog)
          iout_grid_type=2
          
        endif
C         
C         method
C         
        bessel=.false.
        vote=.false.
        
        if(iregrid_type.eq.1.or.iregrid_type.eq.11) then
          iregrid_method=1
        else if(iregrid_type.eq.2.or.iregrid_type.eq.12) then
          iregrid_method=1
          vote=.true.
        else if(iregrid_type.eq.3.or.iregrid_type.eq.13) then
          iregrid_method=2
        else if(iregrid_type.eq.4.or.iregrid_type.eq.14) then
          iregrid_method=2
          bessel=.true.
        end if

        if(vote) then
C         
C         voting defaults --
C         50% of grid box must be covered regardless of the number
C         of candidates
C         
          rmin_vote_max=0.5
          rmin_vote_min=0.5

        endif
        
      endif
C         
C         SPECIAL OPTIONS
C         
C         case 1 - voting options
C         case 2 - user defined beginning lat/lon
C
      if(nargs.gt.3) then
C         
C         voting
C         
        if(vote) then
          
          if (nargs.eq.4) then
C         
C         min and max limits are identical
C         
            if(args(4).lt.0.0.and.args(4).gt.1.0) then
              go to 840
            else
              rmin_vote_max=args(4)
              rmin_vote_min=args(4)
            endif
            
          else if(nargs.ge.5) then
C         
C         different limits
C         
            if(args(4).lt.0.0.or.args(4).gt.1.0.or.
     $           args(5).lt.0.0.or.args(5).gt.1.0) then
              go to 840
            else
              rmin_vote_max=args(4)
              rmin_vote_min=args(5)
            endif

C         
C         set the lat lon bounds with the args 6 and 7
C         if voting is being controlled
C
            if(dxout.lt.0.0.or.dyout.lt.0.0) then
          
              if(dxout.lt.0.0) then
                dxout=abs(dxout)
                rlonbeg_set=args(6)
              endif
            
              if(dyout.lt.0.0) then
                dyout=abs(dyout)
                rlatbeg_set=args(7)
              
              endif
              
            endif

          endif
C         
C         overide of lat/lon bound setting for uniform lat/lon
C         output grids
C         
          
        else if(iregrid_type.le.10) then

          if(dxout.lt.0.0.or.dyout.lt.0.0) then
          
            if(dxout.lt.0.0) then
              dxout=abs(dxout)
              rlonbeg_set=args(4)
            endif
            
            if(dyout.lt.0.0) then
              dyout=abs(dyout)
              rlatbeg_set=args(5)
              
            endif
          else
C         
C         dlon/dlat greater than 0.0, poorly formed user-defined
C         beginning lat/lon point
C         
            go to 842
            
          endif

          
C         
C         too many arguments ; abort
C         
        else
          
          go to  838
          
        endif

      endif

C
C------------------------------------------------------------
C         
C         set up the output grid based on the arguments
C
C------------------------------------------------------------
C
C         iout_grid_type 1 - uniform lat/lon in both directions
C
      if(iout_grid_type.eq.1) then

        rlonbego=ifix(rlonbegi/dxout)*dxout
        rlonendo=ifix(rlonendi/dxout)*dxout
        if(cyclicxi) rlonendo=rlonbego+(360.0-dxout)
        rlatbego=ifix(rlatbegi/dyout)*dyout
        rlatendo=ifix(rlatendi/dyout)*dyout
C         
C         check for a global input gaussian grid
C         
        yleni=yin(njip1)-yin(1)
        if(yleni+dyout.ge.180.0) then
          rlatbego=-90.0
          rlatendo=90.0
        endif

        nio=nint((rlonendo-rlonbego)/dxout)+1
        njo=nint((rlatendo-rlatbego)/dyout)+1

C         
C         make sure we have at least a 2x2 output grid
C 
        if(nio.lt.2) then
          nio=2
          rlonendo=rlonbego+2.0*dxout
        endif
        if(njo.lt.2) then
          njo=2
          rlatendo=rlonbego+2.0*dyout
        endif
C         
C         user defined beginning longitude
C         
        if(rlonbeg_set.ne.-999.0) then
          rlonbego=rlonbeg_set
C         
C         cylic continuity check
C
          if(cyclicxi) then
            nio=nint(360.0/dxout)
          else
            nio=nint((rlonendi-rlonbego)/dxout)+1
          endif
          if(nio.lt.2) nio=2
          rlonendo=rlonbego+(nio-1)*dxout
        endif

        if(rlatbeg_set.ne.-999.0) then
          rlatbego=rlatbeg_set
          njo=nint((rlatendi-rlatbego)/dyout)+1
          if(njo.lt.2) njo=2
          rlatendo=rlatbego+(njo-1)*dyout
C         
C         bounds check
C         
          if(rlatendo.gt.90.0) then
            rlatendo=rlatendo-dyout
            njo=njo-1
          endif

        endif

        if(diag_out) then
          write(iunit_diag,*) ' '
          write(iunit_diag,*) 'iout_grid_type=1'
          write(iunit_diag,*) 
     $         'i= ',rlonbegi,rlonendi,rlatbegi,rlatendi
          write(iunit_diag,*) 
     $         'o= ',rlonbego,rlonendo,rlatbego,rlatendo
          write(iunit_diag,*) 'nio,njo = ',nio,njo
        endif

C         
C         type 2 - gaussian in latitude
C
      else if(iout_grid_type.eq.2) then

        rlonbego=ifix(rlonbegi/dxout)*dxout
        rlonendo=ifix(rlonendi/dxout)*dxout
        if(cyclicxi) rlonendo=rlonbego+(360.0-dxout)
        nio=nint((rlonendo-rlonbego))/dxout+1
C         
C         limits of the input grid (handles gaussian --> gaussian)
C         
        xbegi=yin(1)
        xendi=yin(njip1)

        jj=0
        do j=1,njog
          if(dumx(j).ge.xbegi.and.dumx(j).le.xendi) then
            jj=jj+1
            yout(jj)=dumx(j)
          endif
        end do

        njo=jj
        rlatbego=yout(1)
        rlatendo=yout(njo)

        if(diag_out) then
          write(iunit_diag,*) ' '
          write(iunit_diag,*) 'iout_grid_type=2'
          write(iunit_diag,*) 'dxout = ',dxout
          write(iunit_diag,*)
     $         'i= ',rlonbegi,rlonendi,rlatbegi,rlatendi
          write(iunit_diag,*) 
     $         'o= ',rlonbego,rlonendo,rlatbego,rlatendo
          write(iunit_diag,*) 'nio, njo = ',nio,njo
        endif
 
      else

C         
C         invalid grid type; abort
C
        go to 834

      endif
C         
C         output grid characteristics to GrADS
C         

      if(iregrid_method.eq.1.and..not.vote) 
     $     regrid_method='box averaging                  '
      if(iregrid_method.eq.1.and.vote)
     $     regrid_method='box averaging with VOTING      '
      if(iregrid_method.eq.2.and..not.bessel)
     $     regrid_method='bilinear interpolation         '
      if(iregrid_method.eq.2.and.bessel)
     $     regrid_method='bessel interpolation           '

      if(iout_grid_type.eq.1) then
        write(*,'(a)')  ' '
        write(*,'(a)')
     $       'the output grid is UNIFORM lat/lon:'
        write(*,'(a,f5.2,a,f5.2,a)')
     $       'dx = ',dxout,' deg and dy = ',dyout,' deg'
        write(*,'(a,i3,a,i3)')
     $       '# points in i(lon) = ',nio,
     $       '  # points j(lat) = ',njo
        write(*,'(a,f7.2,a,f7.2,a,f6.2,a,f6.2)') 
     $       'lon extent = ',rlonbego,' to ',rlonendo,
     $       ' lat extend = ',rlatbego,' to ',rlatendo
      endif

      if(iout_grid_type.eq.2) then
        write(*,'(a)')  ' '
        write(*,'(a,i3,a)')
     $       'the output grid is ~ T',nwaves,' GAUSSIAN:'
        write(*,'(a,f7.2,a,f5.2,a)')
     $       'dx = ',dxout,' deg and dy ~ ',dyoutg,' deg'
        write(*,'(a,i3,a,i3)')
     $       '# points in i(lon) = ',nio,'  # points j(lat) = ',njo
        write(*,'(a,f7.2,a,f7.2,a,f6.2,a,f6.2)') 
     $       'lon extent = ',rlonbego,' to ',rlonendo,
     $       '   lat extend = ',rlatbego,' to ',rlatendo
      endif
      write(*,'(a,a)') 
     $     'regrid method is: ',regrid_method
      if(vote) then
        write(*,'(a,f4.2,a,f4.2)')
     $       'vote parameters:  max fract area = ',rmin_vote_max,
     $       '  min frac area = ',rmin_vote_min
      endif
      write(*,'(a)')  ' '

C
C         boundaries of the output boxes
C         
      niop1=nio+1
      njop1=njo+1

C
C         X
C
      do i=1,niop1
        xout(i)=rlonbego+(i-1)*dxout-0.5*dxout
      end do
C         
C         y
C         
C         make sure the output grid does not extend beyond the poles
C         
      if(iout_grid_type.eq.1) then

        do j=1,njop1
          yout(j)=rlatbego+(j-1)*dyout-0.5*dyout
          if(yout(j).lt.-90.0) yout(j)=-90.0
          if(yout(j).gt.90.0) yout(j)=90.0
        end do
C         
C         check for pole points
C         
        if(rlatendo.eq.90.0) npole_point=.true.
        if(rlatbego.eq.-90.0) spole_point=.true.

      else if(iout_grid_type.eq.2) then

        dumx(1)=yout(1)-0.5*(yout(2)-yout(1))
        dumx(njop1)=yout(njo)+0.5*(yout(njo)-yout(njo-1))
        do j=2,njo
          dumx(j)=0.5*(yout(j-1)+yout(j))
        end do

        do j=1,njop1
          yout(j)=dumx(j)
        end do

      endif

C         
C------------------------------------------------------------
C
C         input-output grid box relationship 
C
C------------------------------------------------------------
C         
C         make sure longitudes of the input/output grids
C         are in positive deg
C         
      if(xin(1).lt.0.0) then
        do i=1,niip1
          xin(i)=xin(i)+360.0
        enddo
      endif
         
      if(xout(1).lt.0.0) then
        do i=1,niop1
          xout(i)=xout(i)+360.0
        enddo
      endif
C         
C         calculate the location of the output grid box boundaries 
C         w.r.t. the input grid box boundaries
C
      call in_out_boundaries(xin,yin,nii,nji,
     $     xout,yout,nio,njo,cyclicxi,
     $     niip1,njip1,niop1,njop1,
     $     iunit_diag,
     $     gxout,gyout)
      
      write(iunit_diag,*) 'output grid x boundaries'
      do i=1,niop1
        write(iunit_diag,*) ' i = ',i,' xout = ',xout(i),
     $       ' gxout = ',gxout(i)
      end do

      write(iunit_diag,*) 'output grid y boundaries'
      do j=1,njop1
        write(iunit_diag,*) ' j = ',j,' yout = ',yout(j),
     $       ' gyout = ',gyout(j)
      end do
C         
C         calculate sfc area of each grid box of input grid
C         
      call sfc_area(fld_in,xin,yin,undef,nii,nji,
     $     area_in,iunit_diag)

cc      qtitle='grid box area on unit s'
cc      call qprntn(area_in,qtitle,1,1,nii,nji,4,iunit_diag)
C         
C------------------------------------------------------------
C         
C         do the regrid
C
C------------------------------------------------------------
C
C         box averaging or "clumping" with a "voting" option
C         where the output grid equals the value of the
C         input grid which accounts for the most area.  voting
C         is used for discontinuos data such as soil type
C

      if(iregrid_method.eq.1) then

        call box_ave(fld_in,area_in,undef,gxout,gyout,
     $       nii,nji,nio,njo,fld_out,iunit_diag,vote,istat,
     $       rmin_vote_max,rmin_vote_min)
C         
C         FNOC bilinear/bessel interpolation
C         
      else if(iregrid_method.eq.2) then

        call bssl_interp(fld_in,undef,gxout,gyout,
     $       nii,nji,nio,njo,fld_out,iunit_diag,
     $       cyclicxi,spole_point,npole_point,bessel,istat)

      endif
C         
C         check for pole points
C         
      write(iunit_diag,'(6e13.5)') (fld_out(i,1),i=1,nio)
      if(spole_point.or.npole_point) then
        call fix_poles(fld_out,nio,njo,undef,
     $       spole_point,npole_point)
      endif
      write(iunit_diag,'(6e13.5)') (fld_out(i,1),i=1,nio)

      qtitle='fld_out                '
      call qprntu(fld_out,qtitle,undef,1,1,nio,njo,4,iunit_diag)
C
C------------------------------------------------------------
C
C         write out return info for GrADS
C
C------------------------------------------------------------
C
      ovals(1) = 0.0
      write(iunit_out) ovals
C         
C         modify the appropriate input grid parameters for 
C         defining the output grid for GrADS
C

      vals(4)=float(nio)
      vals(5)=float(njo)

      vals(6)=float(ilinearo)
      vals(7)=float(jlinearo)

      vals(8)=rlonbego
      vals(9)=dxout
      if(iout_grid_type.eq.1) then
        vals(10)=rlatbego
        vals(11)=dyout
      else if(iout_grid_type.eq.2) then
        vals(10)=0.0
        vals(11)=0.0
      endif

C         

      write(iunit_out) vals

      call write_grid(iunit_out,fld_out,nio,njo)
C         
C         define the lat/lon of the output grid points
C         
      do i=1,nio
        xout(i)=rlonbego+(i-1)*dxout
      end do
      write(iunit_out) (xout(i),i=1,nio)

      if(iout_grid_type.eq.1) then

        do j=1,njo
          dumx(j)=(yout(j)+yout(j+1))*0.5
          yout(j)=rlatbego+(j-1)*dyout
        end do

      else if(iout_grid_type.eq.2) then

        do j=1,njo
          dumx(j)=(yout(j)+yout(j+1))*0.5
        end do
        do j=1,njo
          yout(j)=dumx(j)
        end do

      endif

      write(iunit_out) (yout(j),j=1,njo)
c
      go to 999
C         
C         error conditions 
C
 810  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 810'
      write(*,*) 'unable to open the GrADS input file to regrid'
      write(*,*) ' '
      go to 999

 820  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 820'
      write(*,*) 'the first arguement to regrid is not a grid'
      write(*,*) ' '
      go to 999
      
 830  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 830'
      write(*,*) 'the grid must be 2-d and x-y (lon-lat)' 
      write(*,*) ' '
      go to 999

 832  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 832'
      write(*,*) 'arguments must be numbers' 
      write(*,*) ' '
      go to 999

 834  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 834'
      write(*,*) 'invalid output grid iout_grid_type  = ',
     $     iout_grid_type 
      write(*,*) ' '
      go to 999

 836  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 836'
      write(*,*) 'invalid interpolation option',
     $     ' iregrid_type = ',iregrid_type 
      write(*,*) ' '
      go to 999

 838  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 838'
      write(*,*) '4 or more arguments, not a special case'
      write(*,*) ' '
      go to 999

 840  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 840'
      write(*,*) 'rmin_vote outside the interval [1,0]'
      write(*,*) ' '
      go to 999

 842  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 842'
      write(*,*) 'dlon/dlat greater than zero when attempting'
      write(*,*) 'to set beginning lat/lon of the output grid'
      write(*,*) ' '
      go to 999


 999  continue
      stop
      end