      SUBROUTINE ld_sites(success)
c============================================================================
c PROCEDURE NAME: 	_ld_sites
c PROCEDURE TYPE: 	_subroutine
c VERSION DATE:		_23.II.98
c
c PURPOSE:
c
c CALLING PARAMETERS:
c   Name	Type	i/o	units	Description
c   --------	----	---	-----	--------------
c   success      l4      o              successful completion flag
c
c EXTERNAL ROUTINES REQUIRED:
c
c NOTES & CAVEATS:
c
c REVISION HISTORY:
c
c============================================================================
      IMPLICIT NONE
c----------------------------------------------------------------------------
c Compile parameters
c
      INCLUDE "ovpindex.prm"

c----------------------------------------------------------------------------
c Calling parameters
c
      LOGICAL success

c----------------------------------------------------------------------------
c Common blocks
c
      INCLUDE "sites.cmn"
      INCLUDE "ovpindex.cmn"

      REAL*4 v(max_sites)
      INTEGER*4 m(max_sites)
      INTEGER*2 order(max_sites)
      COMMON /workspace/ v, order, m

c----------------------------------------------------------------------------
c Local variables
c
      INTEGER*4 lun, i, j, k, l
      INTEGER*4 sl_len, str_1st_blank
      CHARACTER*1 dum
      REAL*4 temp
      REAL*8 radeg      /57.295779513082323/    !deg-rad conversion

c============================================================================

c-- Write to the log file.

      WRITE(loglun,'(A)')'LD_SITES ...'

c-- Open the site list file.

      CALL get_lun(lun)
      sl_len= str_1st_blank(slfile)-1
      OPEN(UNIT=lun, file=slfile(1:sl_len), 
     &       form='formatted', status='old', err=900)


c-- Read in the sites. (Toss header record first)

      READ(lun,1) dum

      i=1
  100 CONTINUE
        READ(lun,2,end=101,err=901) siteno(i), lat_site(i), lon_site(i)
        write(6,*) siteno(i), lat_site(i), lon_site(i)

        IF(i .EQ. max_sites) GOTO 902
        i= i + 1
        GOTO 100

  101 CONTINUE
      n_sites= i - 1

      CLOSE(lun)
      CALL unget_lun(lun)


c-- Build arrays of coarse-filtering lats and lons, lat_min and lat_max,
c     and lon_range.  Using end-points for the lat speeds up the coarse
c     filtering on lats, while using the center and range for the lons
c     makes it easier to handle the transition at the dateline.

      write(6,*) acceptang
      DO j=1, n_sites
        lat_min(j)= lat_site(j) - acceptang
        lat_max(j)= lat_site(j) + acceptang

        temp= 90.0 - abs(lat_site(j))
        IF(temp .GT. acceptang) THEN
          temp= asin(sin(acceptang/radeg)/sin(temp/radeg))*radeg
        ELSE
          temp= 185.
        ENDIF
        lon_range(j)= temp

      END DO

      

c-- Sort by latitude. (Inefficient, but easy-to-code insertion sort)

c--- Build index array

      order(1)= 1
      DO j=2, n_sites 
        k=1
  201   CONTINUE
        IF(lat_site(j) .LT. lat_site(order(k))) THEN
          DO l=j, k+1, -1 
            order(l)= order(l-1)
          END DO
          order(k)=j
        ELSEIF (k .EQ. j-1) THEN
          order(j)=j 
        ELSE
          k=k+1
          GOTO 201
        END IF
      END DO

c--- Sort the actual data into their own arrays

c---- lats.

      DO i= 1, n_sites
        v(i)= lat_site(i)
      END DO

         DO i=1, n_sites
           lat_site(i)=v(order(i))
         END DO

c---- minimum lats.

      DO i= 1, n_sites
        v(i)= lat_min(i)
      END DO

         DO i=1, n_sites
           lat_min(i)=v(order(i))
         END DO


c---- maximum lats.

      DO i= 1, n_sites
        v(i)= lat_max(i)
      END DO

         DO i=1, n_sites
           lat_max(i)=v(order(i))
         END DO


c---- lons.
     
      DO i= 1, n_sites
        v(i)= lon_site(i)
      END DO

         DO i=1, n_sites
           lon_site(i)=v(order(i))
         END DO
     
c---- lons ranges.
     
      DO i= 1, n_sites
        v(i)= lon_range(i)
      END DO

         DO i=1, n_sites
           lon_range(i)=v(order(i))
         END DO


c---- site ids.

      DO i= 1, n_sites
        m(i)= siteno(i)
      END DO

         DO i=1, n_sites
           siteno(i)=m(order(i))
         END DO





c--If debugging is turned on, write the sorted site list to a file

      IF(debuglevel .GT. 2) THEN
        CALL get_lun(lun)
        OPEN(unit=lun, file='OUTPUT.sitelist',form='formatted')
        WRITE(lun,3)
        DO i=1, n_sites
          WRITE(lun,4) siteno(i), lat_site(i), lat_min(i), lat_max(i),
     &                            lon_site(i), lon_range(i)
        END DO
        CLOSE(lun)
        CALL unget_lun(lun)
      END IF

      success= .TRUE.

      WRITE(loglun,10) n_sites

      RETURN

  900 CONTINUE
      WRITE(6,*)'Could not open file SITELIST for reading.'
      success= .FALSE.
      CALL unget_lun(lun)
      RETURN

  901 CONTINUE
      WRITE(6,*)'Eror while reading SITELIST file at record ', i
      success= .FALSE.
      CLOSE(lun)
      CALL unget_lun(lun)
      RETURN

  902 CONTINUE
      WRITE(6,*)'Attempt to read more sites than <max_sites>.'
      WRITE(6,*)'-- Change parameter in ovpass.prm and recompile.'
      success= .FALSE.
      CLOSE(lun)
      CALL unget_lun(lun)
      RETURN

    1 FORMAT(A)
    2 FORMAT(I3,4X,F6.2,F8.2)
    3 FORMAT('Diagnostic site list (should be sorted by lat.)'/
     &   'Site ID',T20,'Latitude',T60,'Longitude'/
     &   '-------',T20,'--------',T60,'---------'/
     &   T20,'Nominal  Minimum Maximum',T60,'Nominal  Width'/
     &   '#')
    4 FORMAT(I5,T20,3F10.2,T60,2F10.2)
   10 FORMAT('Number of sites to process:',t30,i7)

      END
