C Copyright(c) 1997, Space Science and Engineering Center, UW-Madison
C Refer to "McIDAS Software Acquisition and Distribution Policies"
C in the file  mcidas/data/license.txt

C
C *** $Id: adirserv.f,v 1.33 1998/07/16 19:23:26 krauss Rel $ ***

c ADIRSERV - ADDE image file directory block server

c Remarks:
c      ADIRSERV is intended to be executed via client applications
c      which request image file directory blocks. All communication
c      is via a control block embedded within the client/server
c      transaction protocol. The client interface to this module is
c      through the subroutine MCADIR.


      IMPLICIT NONE

      ! symbolic constants and shared data

      character*8   servername
      parameter (servername = 'ADIRSERV')

      ! external functions

      integer index 
      INTEGER M0CMDPARSE
      integer m0istraceset
      integer m0splice
      integer m0sxdatasetinfo
      integer m0sxgetclientrequest
      integer m0sxresolv

      ! local variables

      CHARACTER*500 REQST               !  full request after resolver
      CHARACTER*500 group               !  full request after resolver
      CHARACTER*500 dataset             !  full request after resolver
      CHARACTER*500 type                !  full request after resolver
      CHARACTER*500 info                !  full request after resolver
      CHARACTER*500 comment             !  full request after resolver
      CHARACTER*120 A                   !  request, or part of it from J
      CHARACTER*72  B                   !  return message for bad transaction
      CHARACTER*32 NAME                 !  name of group/dataset
      CHARACTER*24 PROGRAM              !  constructed server name for exec()
      CHARACTER*12 KIND                 !  data type from resolver

      CHARACTER*500 TMPREQ
      integer LEN_TRIM
      integer LENREQ

      INTEGER I
      INTEGER ISTAT
      INTEGER J(64)                     !  transaction control block
      INTEGER KREQ(125)                 !  REQST as integers, for read
      INTEGER LENCMD
      integer min_range			! minimum area number in dataset
      integer max_range			! maximum area number in dataset
      integer rc
      integer realflag			! needed by m0sxdatasetinfo

      EQUIVALENCE ( KREQ, REQST)
C
C - - ZERO UNUSED FIELDS IN CONTROL BLOCK
      CALL ZEROW(64,J)
C
C - - GET McIDAS ENVIRONMENT
      CALL INITBLOK(I)

      IF(I.NE.0) then

	 b =
     &  'Image directory server unable to initialize McIDAS environment'

         call movcw(b,j(44))
	 j(43) = -117

	 goto 999

      endif
C
C - - READ REQUEST HEADER
      istat = m0sxgetclientrequest (j, reqst)

C     Move beginning of reqst back to J for accounting
      istat = m0istraceset (reqst)
      call mctrace(1,servername, reqst)

      A=REQST
      CALL MOVCW(A,J(11))
C-------Check authorization
      CALL M0VSERV(J)

      istat = m0sxdatasetinfo(REQST,group,dataset,type,kind,info,
     &		comment,min_range,max_range,realflag)
c
c - tack info field back on the end of the request
c
      LENREQ = LEN_TRIM(REQST)
      TMPREQ = REQST(1:LENREQ)//' "'//info
      REQST  = TMPREQ
c
c
c
c--- fail if m0sxdatasetinfo did not succeed!
c
      IF (ISTAT .LT. 0) THEN

	  write(b,FMT='(a27,1x,i4)') 'm0sxdatasetinfo returns: ',
     &						istat
           call mctrace(1,servername, b)

	  rc = m0splice(reqst,name)

	  if (rc .lt. 0) then
	     b = 'Invalid dataset name'//reqst
	     call movcw(b,j(44))
	     j(43) = -118
	     goto 999
	  endif

	  b =
     &  'Image directory server unable to resolve this dataset: '//name

         call mctrace(1,servername,b)

          call movcw(b,j(44))
	  j(43) = -118

	  goto 999

      ENDIF

C
C---------If area, do by internal call; else do by exec
      IF(KIND.EQ.'AREA') THEN
C     - - Initialize accounting
          CALL M0SXLOGI(J)
C     - - PARSE COMMAND
          CALL M0CMDPUT(M0CMDPARSE(REQST,LENCMD))
C==========   CALL ACTUAL SERVER   ==========
          CALL   ADIR(J,min_range,max_range)
c----
c   if this program is calling a Global Imaging (GI) server
c   then you have to do some funny things:  For one, the GI server
c     expects that the REQST string has been mangled by m0sxresolv,
c     so a call to that function is necessary.  Secondly, the subserver
c     program is all upcase, so don't do the mclocase that is done
c     otherwise...

      elseif (kind.eq.'GI') then

	  i = index(KIND,' ')-1
	  program = kind(1:i)//'ADIR'
          KIND='AREA'
          istat = m0sxresolv(REQST,kind)

	  if (istat .lt. 0) then
	     b = 'Invalid dataset name'//reqst
	     call movcw(b,j(44))
	     j(43) = -118
	     goto 999
	  endif

          CALL   M0SUBSERV( J, PROGRAM, REQST)
C----------this return indicates an error
          J(43) = -114
          B='ADIRSERV failed on exec of '//PROGRAM
          CALL MOVCW(B,J(44))
      ELSE
C-------------this will be an exec()
          I=INDEX(KIND,' ')-1
          PROGRAM=KIND(1:I)//'ADIR'
	  call mclocase(program)
          CALL M0SUBSERV( J, PROGRAM, REQST )
C-------------if we are still in control, exec() failed
          J(43)=-114
          B='ADIRSERV failed on exec of '//PROGRAM
          CALL MOVCW(B,J(44))
      ENDIF

 999  CALL M0SXDONE(J)

      END

*$ Name:
*$	adir - Gets area directories matching search criteria
*$
*$ Interface:
*$	subroutine
*$	adir(integer j, integer min_areanumber, integer max_areanumber)
*$
*$ Input:
*$	j		- communications block
*$	min_areanumber	- minimum area number in dataset
*$	max_areanumber	- maximum area number in dataset
*$
*$ Input and Output:
*$	none
*$
*$ Output:
*$	none
*$
*$ Remarks:
*$	260 bytes are sent for each area directory entry.
*$	The first 4 are the absolute number of the entry on the server
*$	The next 4 are the relative number within the cylinder
*$	The rest are words 2 thru 64 of the area directory entry
*$	Word 1 is not sent, since it will always be 0 for a valid entry
*$	and invalid entries just aren't sent.

      subroutine adir(j,min_areanumber,max_areanumber)

      IMPLICIT NONE

c --- include areaparm.inc to get the value MAX_AUXBLOCK_SIZE
      include 'areaparm.inc'

      INTEGER J(64)			! communication block
      INTEGER min_areanumber		! min area number in dataset
      INTEGER max_areanumber		! max area number in dataset

      ! symbolic constants and shared data

      character*8   servername
      parameter (servername = 'ADIRSERV')

      INTEGER MAXBND			! max # of bands
      PARAMETER (MAXBND = 32)

      INTEGER MAXCARDS			! max # of comment cards
      PARAMETER (MAXCARDS = 500)

      INTEGER MAXNAM			! max # of expanded key names
      PARAMETER (MAXNAM = 9)

      INTEGER MAXTYP			! max # of keys per area
      PARAMETER (MAXTYP = 20)


      INTEGER   search_list(MAX_AREA_NUMBER)
      INTEGER   work_array(MAX_AREA_NUMBER)
      COMMON/STORE/search_list

      ! external functions

      CHARACTER*4   CLIT
      CHARACTER*12  CFG
      CHARACTER*12  CFI
      CHARACTER*12  CFU

      integer all_check
      INTEGER ARANAM
      INTEGER ICGET
      INTEGER ISCHAR
      INTEGER KB1OPT
      INTEGER KBPREP
      INTEGER LEN_TRIM
      INTEGER LIT
      INTEGER LLTORA
      INTEGER LUC
      integer lwfile
      INTEGER m0specs
      INTEGER MCARGIHR
      INTEGER MCARGINT
      INTEGER MCARGSTR
      integer mcargiyd
      integer mcdaytimetosec
      integer mciydtocyd
      INTEGER NV1SAE
      INTEGER NVSET
      INTEGER VOLNAM

      LOGICAL BTEST

      ! local variables

      CHARACTER*72  area_path
      CHARACTER*12  area_file
      CHARACTER*256 ARGDUM
      CHARACTER*12  CBLOCK		! block type to transfer
      CHARACTER*80  CTEMP
      CHARACTER*12  EXPTXT
      CHARACTER*12  GETAUX
      CHARACTER*12  KEYEXP(MAXNAM)
      CHARACTER*4   KEYNAM(MAXNAM)

      DOUBLE PRECISION AZIMUTH
      DOUBLE PRECISION RANGE1
      DOUBLE PRECISION RANGE2

      INTEGER ALL
      INTEGER area_dir(64)
      INTEGER area_index
      integer area_number
      integer area_date              ! area directory(5) in ccyyddd format
      INTEGER AUXBUF( MAX_AUXBLOCK_SIZE/4 ) ! buffer to hold AUX block
      INTEGER AUXFLAG                   ! indicates AUX block transfer
      INTEGER BBAND
      INTEGER BDAY
      integer begin_list
      INTEGER BTIM
      INTEGER CARDS(MAXCARDS*20)
      INTEGER EBAND
      INTEGER EDAY
      integer end_list
      INTEGER ETIM
      integer file_ok
      integer filter_index
      integer good_searchs
      INTEGER IBAND
      INTEGER IBUF(64)
      INTEGER IELE
      INTEGER IKEY
      INTEGER ILIN
      INTEGER INICK
      INTEGER IRES
      INTEGER ISTAT
      INTEGER KEYS(MAXTYP)
      INTEGER LEN_CTEMP
      INTEGER min_relativenum
      INTEGER max_relativenum
      INTEGER total_areas
      INTEGER NA
      INTEGER NABS
      INTEGER NAUX			! byte size of AUX block
      INTEGER NB
      INTEGER NBYTES
      INTEGER NCARD
      INTEGER NINT
      INTEGER OLDUC
      INTEGER RELATV
      INTEGER search_ok
      INTEGER SS1
      INTEGER SS2
      INTEGER STYPE
      INTEGER temp_areanum
      INTEGER temp_time
      INTEGER TIME_CHECK


      REAL DG2RAD
      REAL PI
      REAL XDUM
      REAL XELE
      REAL XLAT
      REAL XLAT2
      REAL XLIN
      REAL XLON
      REAL XLON2
      REAL XLNRES
      REAL XLTRES

C------ poke uc to suppress all messages to stdout
 
      olduc = luc(-31)
      call puc(0, -31)

      NA=0

      data keynam /'ALB ',   'BRIT',       'DBZ ',        'RAD ',
     &             'RAW ',   'TE64',       'TEMP',        'VIP ',
     &             'REFL'/

      data keyexp /'ALBEDO', 'BRIGHTNESS', 'DECIBELS',    'RADIANCE',
     &             'RAW',    'TE64',       'TEMPERATURE', 'VIP LEVEL',
     &             'REFLECTIVITY'/

      pi = 4.0 * atan(1.0)
      dg2rad = pi / 180.0

      ALL = lit('ALL ')
      call swbyt4(ALL, 1)

C------ parse command line
 
C------ beginning and ending satellite ids
 
      istat = mcargint(0, 'SS', 1, -1,  0, -1, ss1, argdum)
      istat = mcargint(0, 'SS', 2, ss1, 0, -1, ss2, argdum)

C------ day

      istat = mcargiyd(0, 'DAY', 1, -1,   0, -1, bday, argdum)
      istat = mcargiyd(0, 'DAY', 2, bday, 0, -1, eday, argdum)
 
C------ time
 
      istat = mcargihr(0, 'TIM.E', 1,  -1, 0, -1,btim, argdum)
      istat = mcargihr(0, 'TIM.E', 2,btim, 0, -1,etim, argdum)

C------ band

      istat = mcargint(0, 'BAN.D', 1,   -1, 0, -1, bband, argdum)
      istat = mcargint(0, 'BAN.D', 2,bband, 0, -1, eband, argdum)

C------ max and min relative numbers

      istat = mcargint(0, ' ', 2, 0, 0, -1, min_relativenum, argdum)
      istat = mcargint(0, ' ', 3, 0, 0, -1, max_relativenum, argdum)

      NABS=min_areanumber

      if (max_areanumber .lt. min_areanumber) 
     &     max_areanumber=min_areanumber

      if (max_relativenum .lt. min_relativenum) 
     &    max_relativenum=min_relativenum

C------ block

      AUXFLAG = 0
      istat = mcargstr(0, 'BLOCK', 1, ' ', cblock)
      if( cblock(1:3).eq.'AUX' ) AUXFLAG = 1

C------ Case of specifying a subgroup

      if (min_relativenum .eq. ALL) then
         all_check=1
         relatv = 0
         min_relativenum = 0

C------ Case of positive numbers

      elseif (min_relativenum .gt. 0) then
         relatv = 0
         total_areas = max_areanumber-min_areanumber+1

         if (max_relativenum .gt. total_areas) then
              max_relativenum = total_areas
         endif

         max_areanumber = max_relativenum + min_areanumber - 1
         min_areanumber = min_relativenum + min_areanumber - 1

         if (min_relativenum .gt. max_areanumber) goto 101

C------- Case of time ordered (0 or negative)

      elseif (min_relativenum .le. 0) then
         relatv = 1

C-------Test to see if -negative number is to large, if set number
C-------equal to the number of areas minus 1

         if (min_relativenum .le. -1 * (max_areanumber
     &      -min_areanumber)) then
                 min_relativenum = -1 * (max_areanumber
     &          -min_areanumber)
         endif

      endif

C-------Test bounds against legal area numbers

           if (relatv .eq. 0) then
              if (((min_areanumber .le. 0) .or.
     &             (min_areanumber .gt. MAX_AREA_NUMBER)) .or.
     &            ((max_areanumber .le. 0) .or.
     &             (max_areanumber .gt. MAX_AREA_NUMBER))) goto 102
           endif

C-------- Filter out unwanted areas and sort newest to oldest

      good_searchs = 0
      file_ok = 0


      do 52 area_index=min_areanumber,max_areanumber

           istat = aranam(area_index,area_file)
           istat = volnam(area_file,area_path)

           if (lwfile(area_path) .eq. 1) then
                file_ok = file_ok + 1

                if (file_ok .gt. MAX_AREA_NUMBER) then
                     goto 102
                else
                     call readd(area_index,area_dir)
                endif


C-------- Convert area date + time to total seconds

                if (relatv .eq. 1) then
                     istat = mciydtocyd(area_dir(4),area_date)
                     if (istat .lt. 0) goto 52

                     istat = mcdaytimetosec(area_date, area_dir(5),
     &                                      time_check)
                     if (istat .lt. 0) goto 52
                endif
C-------Test against search conditions

                search_ok=m0specs(ss1,ss2,bday,eday,btim,etim,
     &                          bband,eband,area_dir)
                if (search_ok .eq. 0) then
                     good_searchs = good_searchs + 1
                     search_list(good_searchs)=area_index

C------- Add time to list of times

                     if (relatv .eq. 1) then
                          work_array(good_searchs) = time_check
                     end if

                endif
           endif

52       continue

         if (good_searchs .eq. 0) goto 101

C-------Time sort

           if (relatv .eq. 1) then
                do 61 area_index=1,good_searchs
                     do 51 filter_index=area_index,good_searchs
                         if (work_array(area_index) .ge.
     &                       work_array(filter_index))  goto 51

                         temp_time=work_array(area_index)
                         work_array(area_index)=work_array(filter_index)
                         work_array(filter_index)=temp_time

                         temp_areanum=search_list(area_index)
                         search_list(area_index)=
     &                           search_list(filter_index)
                         search_list(filter_index)=temp_areanum
 51                   continue
 61             continue
           end if

C------Determine the loop bounds for listing for time and non-based
C------requests

           if (relatv .eq. 1) then
               begin_list =  1
               if (good_searchs .le. -min_relativenum) then
                    end_list = good_searchs
               else
                    end_list = -min_relativenum + 1
               endif
           else
               begin_list = 1
               end_list = good_searchs
           endif

C-----Send response.  For each area, directory and comments, if any

      do 10 area_index=begin_list,end_list
           area_number=search_list(area_index)

C-------Get area directory entry

           call readd(area_number,area_dir)

C-------Get comment cards

           ncard=0
           istat = ICGET(0,CARDS)
  5        istat = ICGET(area_number,CARDS(1+20*NCARD))

c
c--- save some room for the cards that will follow
c
           if ((istat .eq. 0) .and. (ncard .lt. maxcards-100)) then
               ncard=ncard+1
               goto 5
           endif
c
c--- attach some cards that the client applications might need
c--- we are putting these values here because there is no room
c--- in the area directory and the comment card section is
c--- currently the only expandable structure we have to pass
c--- variable length information.
c
c
c--- lat and lon at center of image
c

           istat = mcargstr(0, 'AUX', 1, 'N', getaux)

           if (getaux(1:1) .eq. 'Y') then

                ilin = area_dir(6) + ((area_dir(9)/2)  * area_dir(12))
                iele = area_dir(7) + ((area_dir(10)/2) * area_dir(13))

                xlin = ilin
                xele = iele

                istat = nvset('AREA', area_number)

          if (istat .eq. 0) then

                istat = nv1sae(xlin, xele, xdum, xlat, xlon, xdum)

                if (istat .eq. 0) then

	            ctemp = 'Center latitude = '  // cfg(xlat)
	            call movcw(ctemp, CARDS(1+20*ncard))
	            ncard = ncard + 1

	            ctemp = 'Center longitude = ' // cfg(xlon)
	            call movcw(ctemp, CARDS(1+20*ncard))
	            ncard = ncard + 1

                endif

          endif	! end of center lat/lon

c
c--- resolution of image
c
c
c--- only do this operation if the center lat/lon was successful
c
         if (istat .eq. 0) then
c
c--- check for pole as center point
c

            if(nint(abs(xlat)) .eq. 90) then
               ilin = area_dir(6) + area_dir(12)*area_dir(9)/4
               iele = area_dir(7) + area_dir(13)*area_dir(10)/4
               xlin = ilin
               xele = iele
               istat = nv1sae(xlin, xele, 0.0, xlat, xlon, xdum)
            endif

            ilin = ilin + area_dir(12)
            xlin = ilin
            xele = iele
            istat = nv1sae(xlin, xele, 0.0, xlat2, xlon2, xdum)

	    istat = lltora(xlat,xlon,xlat2,xlon2,range1,azimuth)

	    xlat = xlat2
	    xlon = xlon2

            iele = iele + area_dir(13)
            xlin = ilin
            xele = iele
            istat = nv1sae(xlin, xele, 0.0, xlat2, xlon2, xdum)

	    istat = lltora(xlat,xlon,xlat2,xlon2,range2,azimuth)

	    xltres = real((range1 + range2) / 2.0)
	    xlnres = xltres

c
c--- latitude resolution
c
            ctemp = 'Computed Latitude resolution (km) = '//
     &               cfg( real(range1) )
            call bsquez(ctemp)
            call movcw(ctemp, CARDS(1+20*ncard))
            ncard = ncard + 1

            ires = nint(xltres)

            if(ires .eq. 0) then
               ctemp = 'Latitude resolution (km) = ' // cfg(xltres)
            else
               ctemp = 'Latitude resolution (km) = ' // cfi(ires)
            endif

            call bsquez(ctemp)
            call movcw(ctemp, CARDS(1+20*ncard))
            ncard = ncard + 1

c
c--- longitude resolution
c
            ctemp = 'Computed Longitude resolution (km) = '//
     &               cfg( real(range2) )
            call bsquez(ctemp)
            call movcw(ctemp, CARDS(1+20*ncard))
            ncard = ncard + 1

            ires = nint(xlnres)

            if(ires .eq. 0) then
               ctemp = 'Longitude resolution (km) = ' // cfg(xlnres)
            else
               ctemp = 'Longitude resolution (km) = ' // cfi(ires)
            endif

            call bsquez(ctemp)
            call movcw(ctemp, CARDS(1+20*ncard))
            ncard = ncard + 1

         endif	! end of resolution

c
c--- valid calibration types
c

         istat = kbprep(1, area_dir(52))

         if (istat .eq. 0) then
c
c--- get the valid types from kb1opt
c

	    stype = area_dir(52)
	    if (area_dir(57) .ne. 0) stype = area_dir(57)

            ibuf(1)  = area_dir(3)
            ibuf(2)  = area_dir(4)
            ibuf(3)  = area_dir(5)
            ibuf(17) = area_number
            ibuf(37) = stype

            do 300 iband=1,MAXBND
c
c--- work only with valid bands
c
	       if (btest(area_dir(19), iband-1)) then

                  ibuf(4) = iband

                  istat = kb1opt('KEYS', ibuf, keys)

                  if (istat .eq. 0) then

                     do 310 ikey=2,keys(1)+1
c
c--- find an expanded name for the
c--- key if it exists
c
                        exptxt = ' '
                        do 320 inick=1,MAXNAM

		           if (clit(keys(ikey)) .eq. keynam(inick)) then
                              exptxt = keyexp(inick)
                           endif

 320                    continue

                        ctemp = 'Valid calibration unit for band '
     &                       // cfu(iband) // ' = ' // clit(keys(ikey))
			len_ctemp = len_trim( ctemp )
			ctemp = ctemp(1:len_ctemp)//'   "'//exptxt//'"'
                        call movcw(ctemp, CARDS(1+20*ncard))
                        ncard = ncard + 1

 310                 continue

	          endif

               endif

 300        continue

         endif	! end of valid cal types

      endif	! end of aux information

      IF(NCARD.GE.MAXCARDS-2) NCARD=0
      NBYTES=260+80*NCARD
      area_dir(64)=NCARD

C------AUX block transfer
      NAUX = 0
      if ((area_dir(60) .ne. 0) .and. (AUXFLAG .eq. 1)) then
         CALL MCTRACE(1, servername, 'sending the AUX block')
         NAUX = area_dir(34) - area_dir(60)
         IF( NAUX.GT.MAX_AUXBLOCK_SIZE ) GOTO 103
         CALL ARAGET( area_number, area_dir(60), NAUX, AUXBUF )
      endif

C------Account for bytes sent
      J(41)=J(41)+NBYTES+NAUX
      NB=NBYTES+NAUX

C------Record header = length
      CALL SWBYT4(NB,1)
      CALL M0SXSEND(4,NB)
C-------Transform to network format
      area_dir(1)=area_number-NABS+1
      CALL SWBYT4(area_dir(1),20)
      IF(ISCHAR(area_dir(21)).EQ.0) CALL SWBYT4(area_dir(21),1)
      CALL SWBYT4(area_dir(22),3)
      CALL SWBYT4(area_dir(33),19)
      CALL SWBYT4(area_dir(54),1)
      CALL SWBYT4(area_dir(59),6)
C------Transmit
      CALL SWBYT4(area_number,1)
      CALL M0SXSEND(4,area_number)
      CALL M0SXSEND(256,area_dir)
      IF(NCARD.NE.0) CALL M0SXSEND(NCARD*80,CARDS)
      IF(NAUX.NE.0)  CALL M0SXSEND(NAUX,AUXBUF)
 10   CONTINUE
      call puc(olduc, -31)
      RETURN

 101  area_path='No images satisfy the selection criteria'
      J(43)=-51
      CALL MOVCW(area_path,J(44))
      call puc(olduc, -31)
      RETURN

 102  area_path='Maximum number of images for a group exceeded'
      J(43)=-52
      CALL MOVCW(area_path,J(44))
      call puc(olduc, -31)
      RETURN

 103  area_path='Maximum size of AUX block buffer exceeded'
      J(43)=-54
      CALL MOVCW(area_path,J(44))
      call puc(olduc, -31)
      RETURN

      END

*$ Name:
*$      lltora - converts lat/lon to range and azimuth
*$
*$ Interface:
*$      integer function
*$      lltora(float rrlat, float rrlon, float pplat, float pplon,
*$             double range, double azimuth)
*$
*$ Input:
*$      rrlat    - reference latitude
*$      rrlon    - reference longitude
*$      pplat    - measured latitude
*$      pplon    - measured longitude
*$
*$ Input and Output:
*$      none
*$
*$ Output:
*$      range   - distance from reference to measured in km
*$      azimuth - direction from reference to measured in degrees
*$
*$ Return values:
*$       0      - success
*$
*$ Remarks:
*$      This function assumes a west positive convention
*$
*$ Categories:
*$      NAVIGATION, CONVERTER

      INTEGER FUNCTION LLTORA(RRLAT,RRLON,PPLAT,PPLON,RANGE,AZMUTH)

      IMPLICIT NONE

      ! external functions

      ! local variables

      DOUBLE PRECISION	ANG
      DOUBLE PRECISION	ARCL
      DOUBLE PRECISION	AZMUTH
      DOUBLE PRECISION	CPLAT
      DOUBLE PRECISION	CPLON
      DOUBLE PRECISION	CRLAT
      DOUBLE PRECISION	CRLON
      DOUBLE PRECISION	DIST
      DOUBLE PRECISION	PLAT
      DOUBLE PRECISION	PLON
      DOUBLE PRECISION	RANGE
      DOUBLE PRECISION	RLAT
      DOUBLE PRECISION	RLON
      DOUBLE PRECISION	SRLAT
      DOUBLE PRECISION	SRLON
      DOUBLE PRECISION	SPLAT
      DOUBLE PRECISION	SPLON
      DOUBLE PRECISION	XX
      DOUBLE PRECISION	YY
      DOUBLE PRECISION	ZZ

      REAL	COSAZ
      REAL	COSC
      REAL	R
      REAL	PPLAT
      REAL	PPLON
      REAL	RRLAT
      REAL	RRLON
      REAL	TEST
      REAL	Z

      DATA Z/0.017453292/,R/6371./

      RLAT = RRLAT*Z
      RLON = RRLON*Z
      PLAT = PPLAT*Z
      PLON = PPLON*Z

      CRLAT = COS(RLAT)
      CRLON = COS(RLON)
      SRLAT = SIN(RLAT)
      SRLON = SIN(RLON)
      CPLAT = COS(PLAT)
      CPLON = COS(PLON)
      SPLAT = SIN(PLAT)
      SPLON = SIN(PLON)

      XX = CPLAT*CPLON-CRLAT*CRLON
      YY = CPLAT*SPLON-CRLAT*SRLON
      ZZ = SPLAT-SRLAT

      DIST = SQRT(XX*XX+YY*YY+ZZ*ZZ)

      ARCL = 2.*ASIN(DIST/2.)*R
      ANG = 0

      IF(DIST.NE.0.0) ANG = SIN(RLON-PLON)*COS(PLAT)/SIN(DIST)

C     WRITE(*,*) 'LLTORA',ANG

      IF(ABS(ANG).GT.1.0) ANG = DSIGN(1.D0,ANG)
      ANG = ASIN(ANG)/Z
      TEST = SIN(DIST)*CRLAT

      COSC = SRLAT*SPLAT+CRLAT*CPLAT*COS(PLON-RLON)
      COSAZ = 0.0

      IF(TEST.NE.0.0) COSAZ = (SPLAT-COSC*SRLAT)/TEST

C     WRITE(*,*) 'COSAZ',COSAZ

      IF(COSAZ.LT.0.0) ANG = 180.-ANG

      AZMUTH = ANG
      RANGE = ARCL

      LLTORA = 0

      RETURN
      END
