program clhilo
      parameter(nimax=642,njmax=320,nijmax=nimax*njmax)

      real vals(20),ovals(20)
      character fmt*6,pfmt*8,qtitle*24
      character*23 udftdir/'/user/meteo/grads/udft/'/
      logical verb
      
      dimension fld_in(nijmax),xin(nimax),yin(njmax)
      data pfmt/'(      )'/

      verb=.false.
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
c      iunit_diag=6
      iunit_mm=13

      open (iunit_in,file=udftdir//'udf.clhilo.gfi',status='unknown',
     $     form='unformatted',err=810)
        
       if (ios .ne. 0) then
         print *,' Errore di apertura del file di output '
         print *,'udf.clhilo.gfo'
         call exit()
        endif

      open (iunit_out,file=udftdir//'udf.clhilo.gfo',status='unknown',
     $     form='unformatted',iostat=ios)

        if (ios .ne. 0) then
         print *,' Errore di apertura del file di output '
         print *,'udf.clhilo.gfo'
         call exit()
        endif

      open (iunit_diag,file=udftdir//'udf.clhilo.out',
     $     form='formatted',status='unknown',iostat=ios)

        if (ios .ne. 0) then
         print *,' Errore di apertura del file di output '
         print *,'udf.clhilo.out'
         call exit()
        endif

      open (iunit_mm,file=udftdir//'udf.clhilo.vals',
     $     form='formatted',status='unknown',iostat=ios)
        
       if (ios .ne. 0) then
         print *,' Errore di apertura del file di output '
         print *,'udf.clhilo.vals'
         call exit()
        endif

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))
      if(verb)
     $     write(iunit_diag,*) 'number of arguments = ',nargs
c         
c         stop if not one argument
c         
      if(nargs.ne.5) then
        ovals(1)=1.0
        write(*,12) 
 12     format(' ',/
     $       /' ','clhilo requires 5 arguments only!!',
     $       /' ',' ')
        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))
      niji=nii*nji
      xbegi=vals(8)
      dxi=vals(9)
      ybegi=vals(10)
      dyi=vals(11)
      if(verb) 
     $     write(iunit_diag,17)  itype,jtype,nii,nji,niji
   17 format(' ...itype,jtype,nii,nji,niji=',5i5)
      if(verb)
     $     write(iunit_diag,19)  nint(vals(6)),nint(vals(7))
   19 format(' ...vals(6),vals(7)=',2i5)

c         
c         barf if the array is not 2 dimensional
c     
      if(itype.ne.0.or.jtype.ne.1) then
        write(*,22) 
 22     format(' ',/
     $       /' ','clhilo deals only with horizontal fields:',
     $       /' ',' ')
        go to 999
      endif
      
c
c         grid record #2 is the data
c         
      call read_grid(iunit_in,fld_in,niji,istat)
      if(istat.ne.1) go to 820

      qtitle='input field from grads '
      if(verb)
     $     write(iunit_diag,103)  (fld_in(ii),ii=1,10)
  103 format(' ...input field is',/,2x,10e12.3)
c      call qprntu(fld_in,qtitle,undef,1,1,nii,nji,4,iunit_diag)

c     read i/j scaling records

      if(itype .ge. 0)  then
        call read_grid(iunit_in,xin,nii,istat)
      else
        do ii=1,nii
          xin(ii)=xbegi+(ii-1)*dxi
        enddo
      endif
      if(jtype .ge. 0)  then
        call read_grid(iunit_in,yin,nji,istat)
      else
        do jj=1,nji
          yin(jj)=ybegi+(jj-1)*dyi
        enddo
      endif
c         
c------------------------------------------------------------
c
c         read arguments after the grid
c
c------------------------------------------------------------
c         
c      nargs=nargs-1
c      do i=1,nargs
c        read (iunit_in,err=832) args(i)
c      end do
      
      read(iunit_in)  rmm
      read(iunit_in)  rad
      read(iunit_in)  cint
      read(iunit_in)  fmt
      maxmin=nint(rmm)
      if(verb)
     $     write(iunit_diag,1515)  rmm,rad,cint,maxmin,fmt
 1515 format(' ...rmm=',e12.5,'rad=',e12.5,'cint=',e12.5,
     1       ' maxmin=',i4,'  fmt=',a,'...')
      pfmt(2:7)=fmt

      nijo=niji
c
c     process data
c      
      call hilopr(fld_in,nii,nji,maxmin,iunit_mm,xin,yin,rad,cint,
     $     fmt,verb)
      close (unit=iunit_mm,status='keep')
       
      ovals(1)=1.0
      ovals(2)=0.0
      
      write(iunit_out) ovals
      write(iunit_out) vals
      if(verb)
     $     write(iunit_diag,153)  (fld_in(ii),ii=1,10)
  153 format(' ...output field is ',/,2x,10e12.3)
      if(verb)
     $     write(iunit_diag,157)  (vals(ii),ii=1,20)
  157 format(' ...output vals are ',2(/,2x,10e12.3))

      call write_grid(iunit_out,fld_in,nijo)

c     write i/j scaling records

      if(itype .ge. 0)  call write_grid(iunit_out,xin,nii)
      if(jtype .ge. 0)  call write_grid(iunit_out,yin,nji)
      stop
c         
c         error conditions 
c
 810  continue
      ovals(1) = 810.
      write(iunit_diag,*) ' '
      write(iunit_diag,*) 'clhilo error 810'
      write(iunit_diag,*) 'unable to open the grads input file to ',
     1                    'clhilo'
      write(iunit_diag,*) ' '
      go to 999

 820  continue
      ovals(1) = 820.
      write(iunit_diag,*) ' '
      write(iunit_diag,*) 'clhilo error 820'
      write(iunit_diag,*) 'the first argument to regrid is not a grid'
      write(iunit_diag,*) ' '
      go to 999

 832  continue
      ovals(1) = 832.
      write(iunit_diag,*) ' '
      write(iunit_diag,*) 'clhilo error 832'
      write(iunit_diag,*) 'arguments must be numbers' 
      write(iunit_diag,*) ' '
      go to 999
 999  continue
      write(iunit_out) ovals
      stop
      end
      
      subroutine read_grid(iunit,a,nij,istat)
      dimension a(nij)
      istat=1
      read(iunit,err=800) a
      go to 999
 800  continue
      istat=0
 999  continue
      return
      end
      
      subroutine write_grid(iunit,a,nij)
      dimension a(nij)
      write(iunit) a
      return
      end

c
c$$$  subprogram documentation block
c                .      .    .                                       .
c subprogram:    sindb        sine function from degrees input
c   prgmmr: s. j. lord       org: w/nmc22    date: 91-06-06
c
c abstract: sine function from degrees input.
c
c program history log:
c   91-06-06  s. j. lord
c   yy-mm-dd  modifier1   description of change
c   yy-mm-dd  modifier2   description of change
c
c usage:    call pgm-name(inarg1, inarg2, wrkarg, outarg1, ... )
c   input argument list:
c     inarg1   - generic description, including content, units,
c     inarg2   - type.  explain function if control variable.
c
c   output argument list:      (including work arrays)
c     wrkarg   - generic description, etc., as above.
c     outarg1  - explain completely if error return
c     errflag  - even if many lines are needed
c
c   input files:   (delete if no input files in subprogram)
c     ddname1  - generic name & content
c
c   output files:  (delete if no output files in subprogram)
c     ddname2  - generic name & content as above
c     ft06f001 - include if any printout
c
c remarks: list caveats, other helpful hints or information
c
c attributes:
c   language: indicate extensions, compiler options
c   machine:  nas, cyber, whatever
c
c$$$
      function sindb(arg)
c     degrad converts degrees to radians
      data degrad/0.017453/
      sindb=sin(arg*degrad)
      return
      end
c
c$$$  subprogram documentation block
c                .      .    .                                       .
c subprogram:    cosdb        cosine function from degrees input
c   prgmmr: s. j. lord       org: w/nmc22    date: 91-06-06
c
c abstract: returns cosine function from degrees input
c
c program history log:
c   91-06-06  s. j. lord
c   yy-mm-dd  modifier1   description of change
c   yy-mm-dd  modifier2   description of change
c
c usage:    call pgm-name(inarg1, inarg2, wrkarg, outarg1, ... )
c   input argument list:
c     inarg1   - generic description, including content, units,
c     inarg2   - type.  explain function if control variable.
c
c   output argument list:      (including work arrays)
c     wrkarg   - generic description, etc., as above.
c     outarg1  - explain completely if error return
c     errflag  - even if many lines are needed
c
c   input files:   (delete if no input files in subprogram)
c     ddname1  - generic name & content
c
c   output files:  (delete if no output files in subprogram)
c     ddname2  - generic name & content as above
c     ft06f001 - include if any printout
c
c remarks: list caveats, other helpful hints or information
c
c attributes:
c   language: indicate extensions, compiler options
c   machine:  nas, cyber, whatever
c
c$$$
      function cosdb(arg)
c     degrad converts degrees to radians
      data degrad/0.017453/
      cosdb=cos(arg*degrad)
      return
      end


      subroutine hilopr(a,ni,nj,mm,iunit_mm,xin,yin,rad,cint,
     $     fmt,verb)
      parameter (nmax=2000,nmin=nmax)
      dimension a(ni,nj),xin(ni),yin(nj)
      dimension xmax(nmax),ymax(nmax),xmin(nmin),ymin(nmin),
     1          valmax(nmax),valmin(nmin),ihisrt(nmax),ilosrt(nmax)
      character fmt*6,cout*10
      data rmis/-99999999.9/
      logical verb

      rad=rad*1000.
      nhi=0
      nlo=0
      nhimax=0
      nlomax=0
      if(mm .ge. 0)  then
      do j=2,nj-1
        do i=2,ni-1
          if(a(i,j) .ge. a(i+1,j  ) .and.
     1       a(i,j) .ge. a(i+1,j+1) .and.
     2       a(i,j) .ge. a(i  ,j+1) .and.
     3       a(i,j) .ge. a(i-1,j+1) .and.
     4       a(i,j) .gt. a(i-1,j  ) .and.
     5       a(i,j) .gt. a(i-1,j-1) .and.
     6       a(i,j) .gt. a(i  ,j-1) .and.
     7       a(i,j) .gt. a(i+1,j-1))  then
             if(nhi .lt. nmax)  then
               nhi=nhi+1
               dxinumc=a(i+1,j)-a(i-1,j)
               dxidenc=2.0*(a(i+1,j)-2.0*a(i,j)+a(i-1,j))
               if(abs(dxinumc) .gt. 0.5*abs(dxidenc))  then
                 dxi=0.0
                 ximax=xin(i)
               else
                 dxi=-dxinumc/dxidenc
                 ximax=xin(i)+dxi
               endif
               dyjnumc=a(i,j+1)-a(i,j-1)
               dyjdenc=2.0*(a(i,j+1)-2.0*a(i,j)+a(i,j-1))
               if(abs(dyjnumc) .gt. 0.5*abs(dyjdenc))  then
                 dyj=0.0
                 yjmax=yin(j)
               else
                 dyj=-dyjnumc/dyjdenc
                 yjmax=yin(j)+dyj
               endif
               dxinumu=a(i+1,j+1)-a(i-1,j+1)
               dxinumd=a(i+1,j-1)-a(i-1,j-1)
               dxidenu=2.0*(a(i+1,j+1)-2.0*a(i,j+1)+a(i-1,j+1))
               dxidend=2.0*(a(i+1,j-1)-2.0*a(i,j-1)+a(i-1,j-1))
               xintc=a(i,j)+0.5*dxi*(dxinumc+dxi*dxidenc)
               xintu=a(i,j+1)+0.5*dxi*(dxinumu+dxi*dxidenu)
               xintd=a(i,j-1)+0.5*dxi*(dxinumd+dxi*dxidend)
               valij=xintc+0.5*dyj*(xintu-xintd+dyj*(xintu-
     1                              2.0*xintc+xintd))
              
c               write(6,1211)  a(i-1,j+1),a(i,j+1),a(i+1,j+1)
c 1211          format(' ...a(i-1,j+1),a(i,j+1),a(i+1,j+1)=',3f12.4)
c               write(6,1212)  a(i-1,j),a(i,j),a(i+1,j)
c 1212          format(' ...a(i-1,j),a(i,j),a(i+1,j)=',6x,3f12.4)
c               write(6,1213)  a(i-1,j-1),a(i,j-1),a(i+1,j-1)
c 1213          format(' ...a(i-1,j-1),a(i,j-1),a(i+1,j-1)=',3f12.4)
c               write(6,1214)  ximax,yjmax,valij
c 1214          format(' ...ximax,yjmax,valij=',3f12.3)

               xmax(nhi)=ximax
               ymax(nhi)=yjmax
               valmax(nhi)=valij
c               write(6,1316)  nhi,xmax(nhi),ymax(nhi),valmax(nhi)
c 1316          format(' ...nhi,xmax,ymax,valmax=',f12.3)
             else
               write(6,1315)  nmax
 1315          format(' ******nhi exceeds nmax=',i4)
             endif
          endif
          
        enddo
      enddo
      nhimax=nhi

c     sort max values 
      call sortrl(valmax(1),ihisrt(1),nhimax)

      if(verb) then
        write(6,*)  '  #    maxval      index       x       y'
        do nhi=nhimax,1,-1
          write(6,5649)  nhi,valmax(nhi),
     $         ihisrt(nhi),xmax(ihisrt(nhi)),
     $         ymax(ihisrt(nhi))
 5649     format(' ',i4,f12.3,i4,2f12.3)
        enddo
      endif

      do nhi=nhimax,1,-1
        if(valmax(nhi) .ne. rmis)  then
          if(verb)  write(6,5678)  ihisrt(nhi),valmax(nhi)
 5678     format(' ...sorted max, ihisrt,valmax=',i4,f12.3)
          do ntest=nhi-1,1,-1
            if(valmax(ntest) .ne. rmis)  then
              dis=distsp(ymax(ihisrt(nhi)),xmax(ihisrt(nhi)),
     1             ymax(ihisrt(ntest)),xmax(ihisrt(ntest)))
              if(dis .lt. rad .and. abs(valmax(nhi)-
     1             valmax(ntest)) .le. cint)  then
                valmax(ntest)=rmis
                if(verb) write(6,5683)  ntest
 5683           format(' ......valmax for ntest=',i5,' is removed.')
              endif
            endif
          enddo
        endif
      enddo
      endif

      if(mm .le. 0)  then
      do j=2,nj-1
        do i=2,ni-1
          if(a(i,j) .le. a(i+1,j  ) .and.
     1       a(i,j) .le. a(i+1,j+1) .and.
     2       a(i,j) .le. a(i  ,j+1) .and.
     3       a(i,j) .le. a(i-1,j+1) .and.
     4       a(i,j) .lt. a(i-1,j  ) .and.
     5       a(i,j) .lt. a(i-1,j-1) .and.
     6       a(i,j) .lt. a(i  ,j-1) .and.
     7       a(i,j) .lt. a(i+1,j-1))  then
c             write(6,21)  i,j,a(i,j),xin(i),yin(j),fmt,fmt(1:1)
c   21        format(' ...min: i,j,a(i,j)=',2i4,e12.3,2f9.3,' fmt=',a,
c     1              '...fmt(1:1)=',a)

             if(nlo .lt. nmax)  then
               nlo=nlo+1
               dxinumc=a(i+1,j)-a(i-1,j)
               dxidenc=2.0*(a(i+1,j)-2.0*a(i,j)+a(i-1,j))
               if(abs(dxinumc) .gt. 0.5*abs(dxidenc))  then
                 dxi=0.0
                 ximin=xin(i)
               else
                 dxi=-dxinumc/dxidenc
                 ximin=xin(i)+dxi
               endif
               dyjnumc=a(i,j+1)-a(i,j-1)
               dyjdenc=2.0*(a(i,j+1)-2.0*a(i,j)+a(i,j-1))
               if(abs(dyjnumc) .gt. 0.5*abs(dyjdenc))  then
                 dyj=0.0
                 yjmin=yin(j)
               else
                 dyj=-dyjnumc/dyjdenc
                 yjmin=yin(j)+dyj
               endif
               dxinumu=a(i+1,j+1)-a(i-1,j+1)
               dxinumd=a(i+1,j-1)-a(i-1,j-1)
               dxidenu=2.0*(a(i+1,j+1)-2.0*a(i,j+1)+a(i-1,j+1))
               dxidend=2.0*(a(i+1,j-1)-2.0*a(i,j-1)+a(i-1,j-1))
               xintc=a(i,j)+0.5*dxi*(dxinumc+dxi*dxidenc)
               xintu=a(i,j+1)+0.5*dxi*(dxinumu+dxi*dxidenu)
               xintd=a(i,j-1)+0.5*dxi*(dxinumd+dxi*dxidend)
               valij=xintc+0.5*dyj*(xintu-xintd+dyj*(xintu-
     1                              2.0*xintc+xintd))
             
c               write(6,1211)  a(i-1,j+1),a(i,j+1),a(i+1,j+1)
c               write(6,1212)  a(i-1,j),a(i,j),a(i+1,j)
c               write(6,1213)  a(i-1,j-1),a(i,j-1),a(i+1,j-1)
c               write(6,1214)  ximin,yjmin,valij

               xmin(nlo)=ximin
               ymin(nlo)=yjmin
               valmin(nlo)=valij
             else
               write(6,1317)  nmax
 1317          format(' ******nlo exceeds nmax=',i4)
             endif 
          endif
          
        enddo
      enddo
      nlomax=nlo
      call sortrl(valmin(1),ilosrt(1),nlomax)

      if(verb) then
        write(6,*)  '  #    minval      index       x       y'
        do nlo=1,nlomax
          write(6,6649)  nlo,valmin(nlo),
     $         ilosrt(nlo),xmin(ilosrt(nlo)),
     $         ymin(ilosrt(nlo))
 6649     format(' ',i4,f12.3,i4,2f12.3)
        enddo
      endif
      do nlo=1,nlomax
        if(valmin(nlo) .ne. rmis)  then
          if(verb) write(6,6678)  nlo,ilosrt(nlo),valmin(nlo)
 6678     format(' ...sorted min, nlo,ilosrt,valmin=',2i4,f12.3)
          do ntest=nlo+1,nlomax
            if(valmin(ntest) .ne. rmis)  then
              dis=distsp(ymin(ilosrt(nlo)),xmin(ilosrt(nlo)),
     1                   ymin(ilosrt(ntest)),xmin(ilosrt(ntest)))
              if(dis .lt. rad .and. abs(valmin(nlo)-
     1           valmin(ntest)) .le. cint)  then
                 valmin(ntest)=rmis
                 if(verb) write(6,6683)  ntest
 6683            format(' ......valmin for ntest=',i5,' is removed.')
              endif
            endif
          enddo
        endif
      enddo
      endif
c--------------
c     for proximate highs and lows, we decide in favor of the lows
      do nhi=1,nhimax
        if(valmax(nhi) .ne. rmis)  then
          if(verb) write(6,7678)  ihisrt(nhi),valmax(nhi)
 7678     format(' ...sorted max, ihisrt,valmax=',i4,f12.3)
          do ntest=1,nlomax
          if(valmin(ntest) .ne. rmis)  then
            dis=distsp(ymax(ihisrt(nhi)),xmax(ihisrt(nhi)),
     1                 ymin(ilosrt(ntest)),xmin(ilosrt(ntest)))
            if(dis .lt. rad)  then
               valmax(nhi)=rmis
               if(verb) write(6,7683)  nhi
 7683          format(' ......valmax for nhi=',i5,' is removed.')
            endif
          endif
          enddo
        endif
      enddo
c------------------
      do nhi=nhimax,1,-1
        if(valmax(nhi) .ne. rmis)  then
          if(fmt(1:1) .eq. "f" .or. fmt(1:1) .eq. 'f')  then
            write(cout,'('//fmt//')')  valmax(nhi)
            write(iunit_mm,15)  xmax(ihisrt(nhi)),ymax(ihisrt(nhi)),cout
   15       format('max',2f9.3,1x,a) 
          else if (fmt(1:1) .eq. "i" .or. fmt(1:1) .eq. 'i')  then
            write(cout,'('//fmt//')')  nint(valmax(nhi))
            write(iunit_mm,15)  xmax(ihisrt(nhi)),ymax(ihisrt(nhi)),cout
          else
            if(verb) write(6,19)  fmt,fmt(1:1)
   19       format('  ...format not allowed, fmt,fmt(1:1)=',
     1             2(a,'...'))
          endif 
        endif
      enddo
      do nlo=1,nlomax
        if(valmin(nlo) .ne. rmis)  then
          if(fmt(1:1) .eq. "f" .or. fmt(1:1) .eq. 'f')  then
            write(cout,'('//fmt//')')  valmin(nlo)
            write(iunit_mm,25)  xmin(ilosrt(nlo)),ymin(ilosrt(nlo)),cout
   25       format('min',2f9.3,1x,a) 
          else if (fmt(1:1) .eq. "i" .or. fmt(1:1) .eq. 'i')  then
            write(cout,'('//fmt//')')  nint(valmin(nlo))
            write(iunit_mm,25)  xmin(ilosrt(nlo)),ymin(ilosrt(nlo)),cout
          else
            if(verb) write(6,29)  fmt,fmt(1:1)
   29       format('  ...format not allowed, fmt,fmt(1:1)=',
     1             2(a,'...'))
          endif 
        endif
      enddo
      return
      end

c$$$  subprogram documentation block
c                .      .    .                                       .
c subprogram:    sortrl      sorts real numbers
c   prgmmr: s. j. lord       org: w/nmc22    date: 91-06-04
c
c abstract: sorts real numbers.  output array is the index of
c   the input values that are sorted.
c
c program history log:
c   91-06-04  s. j. lord (modified from ncar code)
c
c usage:    call sortrl(a,la,nl)
c   input argument list:
c     a        - array of elements to be sorted.
c     nl       - number of elements to be sorted.
c
c   output argument list:      (including work arrays)
c     la       - integer array containing the index of the sorted
c              - elements.  sorting is from small to large.  e.g.
c              - la(1) contains the index of the smallest element in
c              - array.  la(nl) contains the index of the largest.
c
c
c remarks: none
c
c attributes:
c   language: fortran 77
c   machine:  any
c
c$$$
      subroutine sortrl(a,la,nl)
c
c          entry sortrl(a,la,nl)    sort up real numbers
c  ** revised (6/13/84) for the use in vax-11
c     arguments ...
c       a    input array of nl elements to be sorted or re-ordered
c       la   output array of nl elements in which the original location
c            of the sorted elements of a are saved, or
c            input array to specify the reordering of array a by sorted
c       nl   the number of elements to be treated
c
      save
c
      dimension a(nl),la(nl),ls1(64),ls2(64)
      data nsx/64/,warned/-1./
c
c     set the original order in la
c
    5 do 6 l=1,nl
    6 la(l)=l
c
c     separate negatives from positives
c
   10 l = 0
      m = nl + 1
   12 l = l + 1
      if(l.ge.m) go to 19
   13 if(a(l)) 12,15,15
   15 m = m - 1
      if(l.ge.m) go to 19
   16 if(a(m)) 18,15,15
   18 az = a(m)
      a(m) = a(l)
      a(l) = az
      lz = la(m)
      la(m) = la(l)
      la(l) = lz
      go to 12
   19 l = l - 1
      lnegx=l
 2019 continue
c
c     note that min and max for interval (1,nl) have not been determined
c
      ls1(1) = 0
      l2 = nl + 1
      ns = 1
c
c     step up
c
   20 ls1(ns) = ls1(ns) + 1
      ls2(ns) = l
      ns = ns + 1
      if(ns.gt.nsx) go to 80
      l1 = l + 1
      ls1(ns) = l1
      l2 = l2 - 1
      go to 40
c
c     step down
c
   30 ns=ns-1
      if (ns.le.0) go to 90
      l1 = ls1(ns)
      l2 = ls2(ns)
   40 if(l2.le.l1) go to 30
c
c     find max and min of the interval (l1,l2)
c
   50 if (a(l1)-a(l2) .le. 0) go to 52
      an = a(l2)
      ln = l2
      ax = a(l1)
      lx = l1
      go to 54
   52 an = a(l1)
      ln = l1
      ax = a(l2)
      lx = l2
   54 l1a = l1 + 1
      l2a = l2 - 1
      if(l1a.gt.l2a) go to 60
c
      do 58 l=l1a,l2a
      if (a(l)-ax .gt. 0) go to 56
      if (a(l)-an .ge. 0) go to 58
      an = a(l)
      ln = l
      go to 58
   56 ax = a(l)
      lx = l
   58 continue
c
c     if all elements are equal (an=ax), step down
c
   60 if (an .eq. ax)  go to 30
c
c     place min at l1, and max at l2
c     if either ln=l2 or lx=l1, first exchange l1 and l2
c
      if(ln.eq.l2.or.lx.eq.l1) go to 62
      go to 64
   62 az=a(l1)
      a(l1)=a(l2)
      a(l2)=az
      lz=la(l1)
      la(l1)=la(l2)
      la(l2)=lz
c
c     min to l1, if ln is not at either end
c
   64 if(ln.eq.l1.or.ln.eq.l2) go to 66
      a(ln)=a(l1)
      a(l1)=an
      lz=la(ln)
      la(ln)=la(l1)
      la(l1)=lz
c
c     max to l2, if lx is not at either end
c
   66 if(lx.eq.l2.or.lx.eq.l1) go to 68
      a(lx)=a(l2)
      a(l2)=ax
      lz=la(lx)
      la(lx)=la(l2)
      la(l2)=lz
c
c     if only three elements in (l1,l2), step down.
c
   68 if(l1a.ge.l2a) go to 30
c
c     set a criterion to split the interval (l1a,l2a)
c     ac is an approximate arithmetic average of ax and an,
c     provided that ax is greater than an.  (it is the case, here)
c     ** if a is distributed exponentially, geometric mean may
c        be more efficient
c
      ac = (ax+an)/2
c
c     min at l1 and max at l2 are outside the interval
c
      l = l1
      m = l2
   72 l = l + 1
      if(l.ge.m) go to 78
   73 if (a(l)-ac .le. 0) go to 72
   75 m = m - 1
      if(l.ge.m) go to 78
   76 if (a(m)-ac .gt. 0) go to 75
      az = a(m)
      a(m) = a(l)
      a(l) = az
      lz = la(m)
      la(m) = la(l)
      la(l) = lz
      go to 72
c
c     since 75 is entered only if 73 is false, 75 is not tentative
c     but 72 is tentative, and must be corrected if no false 76 occurs
c
   78 l = l - 1
      go to 20
   80 write(6,85) nsx
   85 format('0 === sorting incomplete. split exceeded',i3,' ===',/)
   90 return
      end
c
c$$$  subprogram documentation block
c                .      .    .                                       .
c subprogram:    distsp      distance on great circle
c   prgmmr: s. j. lord       org: w/nmc22    date: 91-06-06
c
c abstract: calculates distance on great circle between two lat/lon
c   points.
c
c program history log:
c   91-06-06  s. j. lord
c   yy-mm-dd  modifier1   description of change
c   yy-mm-dd  modifier2   description of change
c
c usage:    dxy=distsp(dlat1,dlon1,dlat2,dlon2)
c   input argument list:
c     dlat1    - latitude of point 1 (-90<=lat<=90)
c     dlon1    - longitude of point 1 (-180 to 180 or 0 to 360)
c     dlat2    - latitude of point 2 (-90<=lat<=90)
c     dlon1    - longitude of point 2
c
c
c remarks: distance is in meters
c
c attributes:
c   language: indicate extensions, compiler options
c   machine:  nas, cyber, whatever
c
c$$$
      function distsp(dlat1,dlon1,dlat2,dlon2)
      data rearth/6.37e6/
c
      xxd=cosdb(dlon1-dlon2)*cosdb(dlat1)*cosdb(dlat2)+
     1    sindb(dlat1)*sindb(dlat2)
c
      xxm=amin1(1.0,amax1(-1.0,xxd))
c
      distsp=acos(xxm)*rearth
      return
      end