CTITLE 'REPORT_ERROR'
 
      subroutine report_error(type,ierr,action,file, kill, prog)
 
 
*     Routine to Report that an error has occurred and possibly
*     stop the program (depending on the value of Kill: no kill if
*     0, otherwize kill)
 
*         ierr      - IOSTAT error
*         kill      - Kill program if value is not zero.
*         len_type  - Used length of type
*         len_action    - Used length of action
*         len_file  - Used length of file
*         len_prog  - Used length of prog
*         trimlen   - Get length of string
 
      integer*4 ierr, kill, len_type, len_action, len_file, len_prog,
     .    trimlen
 
*             type  - Type of error
*   action          - Thing being done when error occurred
*   file            - Name of file being acted on
*   prog            - Program/subroutine name of calling
*                   - routine.
 
      character*(*) type, action, file, prog
      
*   Additional variables needed for use of report_stat with this
*   routine

*   rcpar   - Get the name of the program calling this routine
*   len_mod - Length of the module name
*   jerr    - IOSTAT error creating message line
*   message - The concatinated error message
*   module  - name of program

      integer*4 rcpar , len_mod, jerr
      character*128 module
      character*128 message
      
 
****  If error is not zero, print message
      if( ierr.ne.0 ) then          
          call check_ascii( type )
          call check_ascii( action )
          call check_ascii( file )
          call check_ascii( prog )

          len_type = max(1,trimlen(type))
          len_action = max(1,trimlen(action))
          len_file = max(1,trimlen(file))
          len_prog = max(1,trimlen(prog))
 
          if( len_file+len_action+len_type+len_prog.lt.55 ) then
              write(*,100) type(1:len_type), ierr,
     .                     action(1:len_action),
     .                     file(1:len_file), prog(1:len_prog)
  100         format(1x,a,' error ',i6,' occurred ',a,'ing ',a,
     .                    ' in ',a)
          else
              write(*,120) type(1:len_type), ierr,
     .                     action(1:len_action),
     .                     file(1:len_file), prog(1:len_prog)
  120         format(1x,a,' error ',i6,' occurred ',a,'ing',/,
     .                  a,/,' in ',a)
          end if

****      Set up the report_stat string
          len_mod = rcpar(0, module)
          if( len_mod.le.0 ) then
              module = 'MAIN'
              len_mod = 4
          end if
          write(message,140, iostat=jerr) type(1:len_type), 
     .                                action(1:len_action)
  140     format(a,' error ',a,'ing file')
          if( kill.ne.0 ) then
              call report_stat('fatal',module, prog, file, 
     .                         message,ierr)
          else
              call report_stat('warning',module, prog,file, 
     .                          message,ierr)
          end if
 
****      Now see if we should kill
          if( kill.ne.0 ) then
              write(*,200) prog(1:len_prog)
  200         format(' Program terminating in routine ',a)
              stop ' Stopped from report_error'
          end if
      end if
 
***** Thats all
      return
      end
 
CTITLE SUB_CHAR
 
      subroutine sub_char ( buffer, old, new )
 
 
*     Routine to search string buffer and replace all occurrences
*     of string OLD with string NEW
 
*   end     - End position of moving buffer
*   i,j     - Loop counter
*   len_buffer  - length of buffer (full dimensioned length)
*   len_old - Length of old string
*   len_new - Length of new string
*   len_used    - Used portion of buffer
*   offset  - Difference in length between new and old string
*   pos     - Position of old in string
*   start   - start position for moving buffer
*   trimlen - Length of used portion of string
 
      integer*4 end, i, len_buffer, len_old, len_new, len_used,
     .    offset, pos, start, trimlen
 
*   finished    - Indicates we have run out of string
 
      logical finished
 
*   buffer  - Buffer to be manipulated
*   old     - String to be replaced
*   new     - string to replace old with
 
 
      character*(*) buffer, old, new
 
***** Start by get sizes of strings
 
      len_old = len( old )
      len_new = len( new )
      len_buffer = len( buffer )
 
*     Check to make sure old is not a substring of new.  If it is the
*     algorithm used here will not work
 
      pos = index ( new, old )
      if( pos.gt.0 ) RETURN
 
      finished = .false.
 
      do while ( .not. finished )
 
          pos = index( buffer, old )
 
*                             ! Make substiution
          if( pos.gt.0 ) then
 
*             Now copy string after pos
              len_used = trimlen( buffer )
 
              offset = len_new - len_old
              start  = pos - offset + 1
              end    = min(len_used,len_buffer+offset)
 
*                                             ! Move string starting at top
              if( len_new.gt.len_old ) then
                  do i = end, start, -1
                      buffer(i+offset:i+offset) = buffer(i:i)
                  end do
              end if
 
*                                             ! Move string starting a bottom
              if( len_new.lt.len_old ) then
                  do i = start, end
                      buffer(i+offset:i+offset) = buffer(i:i)
                  end do
*                                             ! Clear the characters left at
                  buffer(end+offset+1:) = ' '
*                                             ! the end
              end if
 
****          Now put in new string
              buffer(pos:pos+len_new-1) = new
 
*                     ! we have finished
          ELSE
 
              finished = .true.
 
          ENDIF
*                     ! Searching through string
      END DO
 
****  Thats all
      return
      end
 
CTITLE 'TRIMLEN'
 
 
      integer*4 function trimlen(string)
 

* MOD TAH 050129: Adding counting null's as well as blanks at the 
*     end of strings.  Useful for strings returned from C that have
*     nulls at the end.
 
*     Routine to return the length of the used portion of string.
*     Length is with trailing blanks removed.
 
*         len_string    - declared length of string
*         i             - Loop counter
 
      integer*4 len_string, i
 
*       blanks          - Indicates blanks are still being
*                       - found
 
      logical blanks
 
*             string    - the string whose length we want
 
      character*(*) string
 
***** Get full length of string, and work backwards
 
      len_string = LEN(string)
      i = len_string
 
      blanks = .true.
 
*     Scan string from end to front
      do while ( blanks )
          if( i.gt.0 ) then
* MOD TAH 050129: Add or check on null
              if( string(i:i).eq.' ' .or.
     .            string(i:i).eq.char(0) ) then
                  i = i - 1
              else
                  blanks = .false.
              end if
          else
*                                     ! Force exit at end of string
              blanks = .false.
          end if
      end do
 
***** Save length of string
      trimlen = i
 
***** Thats all
      return
      end
 
 
CTITLE 'YMDHMS_TO_JD'
 
      SUBROUTINE YMDHMS_to_JD ( date, seconds, epoch )
*     -------------------------
*    .,      26JUL85                <871210.1317>
 
*     Author: T.Herring            3:53 PM  THU., 10  JULY, 1986
*
*$S "YMDHMS_to_JD"
*-----------------------------------------------------------------------
*     Routine to convert a calender date with hours, minutes and
*     seconds to a Modified Julian date. The calender date is
*     ordered as year, month, day of month, hours, and
*     minutes.  These values are stored in a single I*2 array.
*
*     If the year is greater than 1000 then the it is assumed to
*     contain the centuries.
*
*     This routine is only valid for dates after 1600 Jan 0.
*
*     CALLING SEQUENCE:
*     =================
*     CALL ymdhms_to_JD ( date, seconds, epoch)
*
*     WHERE:
*     date    is a 5 element array containing year, month, day of month,
*             hour, and minutes.
*             (I*4 5 element array INPUT)
*     seconds is the seconds part of the epoch
*             (REAL*8 INPUT)
*     epoch   is the JD with fractional days corresponding to date and
*             seconds.
*             (REAL*8 OUTPUT)
*
*----------------------------------------------------------------------
*$E
 
 
*         day           - day of month
*         date(5)     - 5 element date array with year, month, day,
*               - hours and minutes.
 
*   days_to_month(12)   - number of days from start of year to
*               - each month
 
*   leap_days   - the number of leap days we need to include
 
*   month       - month of year
*   year        - the years since 1900.0
 
*   years_from_1600 - number of years since 1600 Jan 0.
 
*   days_from_1600  - number of days since 1600 Jan 0.
 
      integer*4 day, date(5), days_to_month(12), leap_days, month,
     .    year, years_from_1600, days_from_1600
 
*       epoch   - The JD corresponding to date and seconds
*   fraction    - Fraction of a day.  We compute this separately
*               - to avoiod some rounding error when added to MJD
*   mjd         - the computes MJD with fractional days
*   seconds     - the seconds part of the epoch.
 
      real*8 epoch, fraction, mjd, seconds
 
*       leap_year   - Indicates that this is a leap year
 
      logical leap_year
 
      data  days_to_month /   0,  31,  59,  90, 120, 151,
     .                      181, 212, 243, 273, 304, 334 /
 
***** START, Make sure year is from 1900
 
      year      = date(1)
      month     = date(2)
      day       = date(3)

* Y2K mod to map small number years to resonabale values 
      if( year.lt.50 ) then
          year = year + 2000
      else if ( year.lt. 200 ) then
          year = year + 1900
      end if
 
***** Compute number of years from 1600
 
      years_from_1600 = year - 1600
 
***** Now compute number of leap days upto the start of the year
*     we are in (i.e., if the year we are in is a leap year do not
*     add the leap day yet)
 
      leap_days =   (years_from_1600 -   1)/  4
     .            - (years_from_1600 +  99)/100
     .            + (years_from_1600 + 399)/400  + 1
 
      if( years_from_1600.eq.0 ) leap_days = leap_days - 1
 
***** Now see if we are in leap year
 
      leap_year = .false.
      if(   mod(years_from_1600,  4).eq.0     .and.
     .     (mod(years_from_1600,100).ne.0.or.
     .      mod(years_from_1600,400).eq.0)       ) leap_year = .true.
 
***** Now compute number of days sinec 1600
 
      days_from_1600 = years_from_1600*365  + leap_days +
     .                 days_to_month(month) + day
 
***** Add extra day if we are after Februrary and this is a leap year
      if( month.gt.2 .and. leap_year ) then
          days_from_1600 = days_from_1600 + 1
      end if
 
***** Compute the mjd and add in the fraction of a day part
*     The MJD of 1600 Jan 0 is -94554
 
      fraction  = seconds/86400.d0  + date(5)/1440.d0 + date(4)/24.d0
 
      mjd   = -94554.d0 + days_from_1600 + fraction
      epoch = mjd + 2 400 000.5d0
 
***** THATS ALL
      RETURN
      end
 
CTITLE 'CASEFOLD'
 
      subroutine casefold( string )
 
 
*     Routine to convert a string to upper case.
 
*         i     - Loop counter
*   len_string  - used length of string
*   trimlen     - returns string length
 
      integer*4 i, len_string, trimlen
 
*             string    - Used portion of string
 
      character*(*) string
 
***** Get length of string
      len_string = trimlen(string)
 
*     Now loop and convert
      do i = 1, len_string
          if( string(i:i).ge.'a' .and. string(i:i).le.'z') then
              string(i:i) = CHAR( ICHAR(string(i:i))-32 )
          end if
      end do
 
****  Thats all
      return
      end


CTITLE 'caseunfold'
 
      subroutine caseunfold( string )
 
 
*     Routine to convert a string to lower case.
 
*         i     - Loop counter
*   len_string  - used length of string
*   trimlen     - returns string length
 
      integer*4 i, len_string, trimlen
 
*             string    - Used portion of string
 
      character*(*) string
 
***** Get length of string
      len_string = trimlen(string)
 
*     Now loop and convert
      do i = 1, len_string
          if( string(i:i).ge.'A' .and. string(i:i).le.'Z') then
              string(i:i) = CHAR( ICHAR(string(i:i))+32 )
          end if
      end do
 
****  Thats all
      return
      end
 

 
CTITLE CHECK_ASCII

      subroutine check_ascii(string)

*     Routine to check a string for non-ascii characters.  Any 
*     non-ascii characters are replaced with the @ character

* PASSED Variables

* string  - The string to checked.  (ichar should be between 32-127
*      for each character in the string.  If non-ascii characters
*      found, all charcaters after 32 are set to blanks.

      character*(*) string

* LOCAL variables

*   i  - loop counter
*   ascii_ichar   -  Value of ascii character (nulls are set to blank)

       integer*4 i, ascii_ichar

*   non_ascii - Set true if non-ascii characters are found

       logical non_ascii

*****  Set non_ascii false to start
       non_ascii = .false.
       do i = 1, len(string)
           ascii_ichar = ichar(string(i:i))
          if( ascii_ichar.eq.0 )  then
              string(i:i) = ' '
              ascii_ichar = 32
           else if ( ascii_ichar.lt.32 .or.
     .               ascii_ichar.ge.128 ) then
              string(i:i) = '@'
              non_ascii = .true.
           end if
       end do

*****  If we had non-ascii characters then clear the end of the
*      string
       if( non_ascii .and. len(string).gt.32 ) then
           string(32:)  = ' '
       end if

*****  Thats all
       return
       end

CTITLE 'JD_TO_YMDHMS'
 
      SUBROUTINE JD_to_YMDHMS ( epoch, date, seconds )
*     ------------------------
*    .,      10JUL86                <871210.1317>
 
*     Author: T.Herring           11:28 AM  THU., 10  JULY, 1986
*
*$S "JD_to_YMDHMS"
*----------------------------------------------------------------------
*     Routine to convert a Julian date (with fractions of a day)
*     or a full Julian date, to a calender data with hours and minutes,
*     and floating-point seconds (Real*8).
*
*     NOTE: If a full Julian date is used the resolution of the seconds
*     will only be about 10 microseconds,  a MJD should yield a resolution
*     of about 0.1 microseconds in the seconds.
*
*     This routine is only valid for dates after 1600 Jan 0 (MJD -94554)
*
*     CALLING SEQUENCE:
*     =================
*     CALL JD_to_YMDHMS ( epoch, date, seconds )
*
*     WHERE:
*     epoch   is a modified or full julian date with possibly a fractional
*             part of a day.  The two type of input are destinguished by
*             the full Julian date being greater than 2,000,000.
*             (REAL*8 INPUT)
*     date    is an array containing the calender date with (full) year,
*             month of year, day of month, hours, minutes. Note the year
*             will be returned with the centuries added.
*             (I*4 5 element array OUTPUT)
*     seconds is the floating point seconds part of the MJD
*             (REAL*8 OUTPUT)
*
*----------------------------------------------------------------------
*
*$E
 
 
 
*         century     - Number of century from 1600 Jan 0.
*         date(5)     - the calender date corresponding to
*                     - 'epoch' resolution to the minute
*   day         - day of month
*   day_of_year - Number of days from start of year
*   days_to_month(13)   - number of days to start of each month
*               - in a non-leap-year
 
*   month       - month of year
 
*   year        - years since start of century
*   years_from_1600 - Number of years since 1600 Jan 0.
 
*   days_from_1600  - Number of days elapsed since 1600 Jan 0.
*               - (MJD -94554.0 Julian date 2305447.0)
      integer*4 century, date(5), day, day_of_year, days_to_month(13),
     .    month, year, years_from_1600, days_from_1600
 
*      epoch    - the julian date or modified julian date
*               - to be converted to calender date
*   fraction    - the fraction of a day part of MJD
 
*   mjd         - epoch converted to a MJD
*   mjd_day     - the whole number days in the mjd
 
*   seconds     - the seconds part of the MJD (<60)
 
      real*8 epoch, fraction, mjd, mjd_day, seconds
 
*       leap_year   - Indicates that this days is a leap year
 
      logical leap_year
 
      data days_to_month /   0,  31,  59,  90, 120, 151, 181,
     .                     212, 243, 273, 304, 334, 365       /
 
 
***** START, convert full jd to mjd if we have too.
*                                      ! epoch must be Full julian date
      if( epoch.gt.2000000.d0 ) then
*                                      ! convert to MJD.
 
          mjd = epoch - 2 400 000.5d0
      else
          mjd = epoch
      end if
 
***** Remove the fractional part of the mjd
 
      mjd_day  = aint ( mjd )
      fraction =  mod ( mjd,1.d0 )
 
      if( mjd.lt.0 .and. fraction.ne.0.d0 ) then
          mjd_day  = mjd_day - 1
          fraction = fraction + 1.d0
      end if
 
***** Now convert MJD (even day) to date (year, month, day )
*     Get number of days since 1600.
 
      days_from_1600 = mjd_day - (-94554.0d0)
      years_from_1600 = days_from_1600/365.d0
 
***** Now compute day_of_year and years_from_1600 accounting for
*     leap years
 
*                        ! Just to get us into loop
      day_of_year = 0
      do while ( day_of_year.le.0 )
 
          century = years_from_1600/100
 
          day_of_year =  days_from_1600 - years_from_1600*365
     .                   - (years_from_1600 -   1)/  4
     .                   + (years_from_1600 +  99)/100
     .                   - (years_from_1600 + 399)/400 - 1
 
*         If we are 1600 then add one day
          if( years_from_1600.eq.0 ) then
              day_of_year = day_of_year + 1
          end if
 
*         See if the leap days have taken us to a earlier year
          if( day_of_year.le.0 ) then
              years_from_1600 = years_from_1600 - 1
          end if
      end do
 
***** We now have number of days from start of year and the year
*     Convert years back to start of century
 
      year    = mod( years_from_1600, 100)
 
***** See if this is a leap year
 
      leap_year = .false.
*                             ! we are at beginning of century
      if( year.eq.0 ) then
          if( mod(century,4).eq.0 ) leap_year = .true.
      else
          if( mod(year,4)   .eq.0 ) leap_year = .true.
      end if
 
***** If the day of year is less than 60 then the leap years do no not
*     matter,  if we are greater than or equal to 60, need to account
*     for the leap years
 
*                                     ! Dont worry about leap years
      IF( day_of_year.lt.60 ) THEN
          if( day_of_year.le.31 ) then
              month = 1
              day   = day_of_year
*                 ! we are in February
          else
              month = 2
              day   = day_of_year - 31
          end if
*                                     ! Need to account for leap years
      ELSE
          if( leap_year .and. day_of_year.eq.60 ) then
              month  = 2
              day    = 29
          else
              if( leap_year ) day_of_year = day_of_year - 1
 
*****         Now find month
              month = 2
              do while ( day_of_year.gt. days_to_month(month) )
                  month = month + 1
              end do
              month = month - 1
              day   = day_of_year - days_to_month(month)
          end if
      END IF
 
***** Now save the date
 
      date(1) = years_from_1600 + 1600
      date(2) = month
      date(3) = day
 
***** Now convert the fraction of a day to hours, minutes and seconds
 
      date(4) = fraction*24.d0
      date(5) = fraction*1440.d0 - date(4)*60.d0
 
      seconds = 86400.d0*fraction - date(4)*3600.d0
     .                            - date(5)*60.d0
 
***** Thats all
      RETURN
      end
 
CTITLE 'PROPER_RUNSTRING'
 
      subroutine proper_runstring( help_file, program, kill)
 
 
*     Routine to ouput the runstring for nutcorr
*     To use this subroutine the user needs to set the enviroment
*     variable HELP_DIR. This directory is then preprended to
*     the file name
 
*           ierr        - IOSTAT error (opening help file)
*   trimlen             - Length of string
*   kill                - indicates iof program should be killed.
*                         0 means do not stop program
*                        >0 kill program
*                        <0 kill program but list help in 24 line
*                           segment (with q to quit)
*   len_dir             - Length of directory name
*   cnt                 - Keeps track of number of lines written
*                        so that user can be quizzed to continue if
*                        kill set < 0
 
      integer*4 ierr, trimlen, kill, len_dir, cnt
 
*            line       - Line read from input.
 
      character*100 line
 
*             help_file - Name of help file
*       program         - Name of program.
 
      character*(*) help_file, program

*     Full_help_file  - Help file name with directory added.

      character*256 Full_help_file

*     ans - Answer to query for more lines
      character*4 ans
 
***** open the help file.  First construct name
      call getenv('HELP_DIR',Full_help_file)
      len_dir = trimlen(Full_help_file)
      if( len_dir.gt.0 ) then
           Full_help_file(len_dir+1:) = help_file
      else
           Full_help_file = help_file
      end if
 
      open(999, file=full_help_file, iostat=ierr, status='old')
*     NOD TAH 980813: If file not found tell user that to do
      if( ierr.ne.0 ) then
          write(*,110)
 110      format(/,'**ERROR** Opening help file.  Check that',
     .             ' enviroment variable HELP_DIR points to',/,
     .             '          directory with help files.')
      end if
      call report_error('IOSTAT',ierr,'open',full_help_file, 0,
     .                  program)
 
*                             ! Loop until EOF
      cnt = 0
      do while ( ierr.eq.0 )
          read(999,'(a)', iostat=ierr) line
          cnt = cnt + 1
          if( ierr.eq.0 ) write(*,'(a)') line(1:max(trimlen(line),1))
          if( mod(cnt,24).eq.0 .and. kill.lt.0 ) then

*             see if we should continue list
              write(*,120)
 120          format('More? (q to quit) ',$)
              read(*,'(a)') ans
              if( ans(1:1).eq.'q' .or. ans(1:1).eq.'Q' ) then
                 ierr = -1
              end if
          end if
      end do
 
***** close the input
      close(999)
 
****  Now see if we should stop
      if( kill.ne.0 ) then
          write(*,200) program(1:max(trimlen(program),1))
  200     format(' PROGRAM ',a,' terminating')
          stop ' Incomplete runstring'
      end if
 
****  Thats all
      return
      end
 
CTITLE RCPAR
 
      integer*4 function rcpar( iel, arg )

 
*     Routine to emulate RCPAR using the igetarg UNIX subroutine
*     Modified to use getarg
 
*         iel       - Element of runstring to get
*       igetarg     - UNIX routine to read runstring
*       len_arg     - Length of arg.
*       trimlen     - Get length of string
*       offset      - Offset to be applied to the passsed element
*                    (0 is assumed to program name, 1 first argument)
 
      integer*4 iel, len_arg, trimlen, offset
 
*             arg   - Arg of runstring
 
      character*(*) arg
      character*4 test_arg
      
      data offset / -1 /
 
****  Get length of argument and runstring
* MOD TAH 010610: To see where the count starts for getarg
      if( offset.lt.0 ) then
          call getarg(0, test_arg)
	  len_arg = trimlen(test_arg)
	  if( len_arg.eq.0 ) then
	      offset = 1
	  else
	      offset = 0
	  end if
      end if
      
      len_arg = LEN(arg)
      call getarg( iel+offset, arg )
      rcpar = trimlen( arg )
 
***** Thats all
      return
      end
 
CTITLE report_stat
     
      subroutine report_stat( typep, modulep, routinep, filep,
     .            messagep, ierr )
 
*     Routine to report on the status on gamit and globk processing.
*     The routine reports status messages which indicate where in the
*     processing we are; warnings which are are not fatal to the processing
*     but could be used as diagnostic of later problems; and fatal errors
*     which cause the program to stop running.
 
*     For all reports, the messsages are appended to files named
*     'module' .status, .warn or .fatal depending on the type of error
*     These files are opened with append status and closed before the
*     routine ends (the later is to ensure that the UNIX file buffers
*     flushed so that the user can see the current status).
 
*     Standard examples of the use of this routine would be:
*     STATUS:
*     call report_stat('STATUS','MODEL','model',' ','Beginning process',0)
*     WARNING:
*     call report_stat('WARNING','MODEL','wrtout',' ',
*    .                  'Site radius too large', 0)
*     call report_stat('WARNING','GLOBK','update_aprioris',line,
*    .                   'IOSTAT error reading string ',ierr)
*     FATAL:
*     call report_stat('FATAL','SOLVE','invert',' ',
*    .                   'Normal equations non-positive definite',0)
*     call report_stat('FATAL','GLOBK','read_glb_mar',glb_mar_file,
*    .                   'IOSTAT error opening file ',ierr)
*     CLEAR:
*     call report_stat('CLEAR','GLOBK',' ',' ',' ',0)

* ROUTINES used from HP1000 library:
* trimlen    - Length of string
* casefold   - Converts string to upper case
* caseunfold - Converts string to lower case
* systime  - Returns system time (available on both HP and Sun)

*
* PASSED Variables
 
*   ierr    - Error number associated with report.  This value is
*           - printed for warnings and fatal errors when non-zero.
 
      integer*4 ierr
 
*   type    - Type of message:  Only three types are allowed, any
*           - other type will cause this routine to print a warning
*           - message.  The types are (only first four characters
*           - checked:
*           - status  -- General reporting of status for progress
*           -       information
*           - warning -- Non-fatal error in the processing, could be
*           -       diagonistic of later fatal errors
*           - fatal   -- Fatal error, If ierr=0 then only the message
*           -       is printed, otherwize the error number is
*           -       printed with the message.
*           - clear   -- Removes the status, warning and fatal files
*                   for a given module name.
 
*   module  - Name of module (main program) (e.g., model, solve)
*   subroutine  - Name of the subroutine from which the report is
*           - made.
*   file    - Name of the file or string being manipulated.  This is
*             printed whenever a non-zero length string is passed.
*             It is included with (Name xxxxx) at the end of the line.
*   message - Message to be printed
                                              
      character*(*) typep, modulep, routinep, filep, messagep
      character*(256) type, module, routine, file, message

* LOCAL VARIABLES
 
*   rep_date(5) - Yr, mon, day, hr, min that the report is made
*   trimlen     - Library routine to returned the non-blank length
*               - of a string
 
*   len_mod     - Length of the module string
*   len_sub     - Length of the subroutine string
*   len_mes     - Length of the message
*   len_fil     - Length of the file name
*   len_typ     - Length of the status type (variable type)
 
*   jerr        - IOSTAT error writing to the status file
*   kerr        - IOSTAT error writing to stdout
*   cerr        - IOSTAT error closing the status file
 
*   rep_unit        - Unit number assigned to the report file.  If we have
*               - trouble with the opening the unit, rep_unit is set
*               - 6 so that output will goto to stdout.

*   i           - Counter used to find / or start of string in the module
*                 name

      integer*4 rep_date(5), trimlen, len_mod, len_sub, len_mes,
     .    len_fil, len_typ, jerr, kerr, cerr, rep_unit, i
 
*   rep_sec     - Seconds that the report is made
 
      real*8 rep_sec
 
*   rep_file        - Name of the report file to which the messages are
*               - written.  Generated from the module name with an an
*               - extent based on the type.  (NOTE: Module is upper
*               - case when name generated).
*   mod_name - Names of the module for creating the file name
*                 (Removes the full information if present)
 
 
      character*256 rep_file, mod_name
      character*6 prog_name
         
***** Make a copy of all passed variables.  Stops us overwritting
*     a fixed string.
      type    = typep
      module  = modulep
      routine = routinep
      file    = filep
      message = messagep

****  Start: First get the system time so that we know when the report
*     was made
 
      call systime( rep_date, rep_sec )

*     Remove the 19-- portiion of the date so that HP and Sun look the
*     same
      rep_date(1) = rep_date(1) - int(rep_date(1)/100)*100
 
*     See what type of report we are to make (this saves
*     unnecessary opening of files).  Casefold the type incase user
*     calls with the wrong case.  Module and subroutine are set case for
*     "nice" looking output.
 
      call casefold( type )
      call casefold( module )
      call caseunfold( routine )
 
*     Get the lengths of the strings.  (1 is minumum length to stop fortran
*     error is string(1:0) is attempted.)

      if( ichar(module(1:1)) .eq. 0 ) module(1:1)  = ' '
      if( ichar(routine(1:1)).eq. 0 ) routine(1:1) = ' '
      if( ichar(message(1:1)).eq. 0 ) message(1:1) = ' '
      if( ichar(file(1:1))   .eq. 0 ) file(1:1)    = ' '
      if( ichar(type(1:1))   .eq. 0 ) type(1:1)    = ' '

*     Make sure all the strings contain ascii characters
      call check_ascii( module )
      call check_ascii( routine )
      call check_ascii( message )
      call check_ascii( file )
      call check_ascii( type )
 
      len_mod = max(1,trimlen(module))
      len_sub = max(1,trimlen(routine))
      len_mes = max(1,trimlen(message))
      len_fil = max(1,trimlen(file))
      len_typ = max(1,trimlen(type))
*
*     Now strip out the module name from the end the module string.
*     (Look for /'s in the name in case full name passed)
      i = len_mod
      do while ( i.gt.1 .and. module(i:i).ne.'/' )
          i = i - 1
      end do
      if( module(i:i).eq.'/' ) i = i + 1
      mod_name = module(i:len_mod) 
*     The following will result in the report files being named 'GAMIT'
*     for the main GAMIT batch sequence, but individual names for the
*     non-batch modules (MAKEXP, MAKEX, FIXDRV, and utilities) and
*     for the /kf routines except for AUTCLN.  Currently, the GAMIT 
*     batch script erases the 'GAMIT' files but all the other modules
*     erase their files upon executation with a 'CLEAR' call of report_stat.

      if( mod_name(1:3).eq.'ARC'.or. 
     .    mod_name(1:6).eq.'YAWTAB'.or.
     .    mod_name(1:6).eq.'OCTTAB'.or.
     .    mod_name(1:5).eq.'MODEL'.or. 
     .    mod_name(1:6).eq.'AUTECL' .or.
     .    mod_name(1:6).eq.'AUTCLN'.or.
     .    mod_name(1:5).eq.'CFMRG'.or.
     .    mod_name(1:5).eq.'SOLVE'.or.
     .    mod_name(1:6).eq.'SCANDD' ) then
           prog_name = 'GAMIT' 
      else
           prog_name = mod_name(1:6)
      endif
      
*     Set the error flags for errors in this routine to zero.  At end we
*     check these are still zero.
      jerr = 0
      kerr = 0
      cerr = 0
 
*****  Report the appropriate type of messages.
*     STATUS report
      if( type(1:4).eq.'STAT' ) then
 
*         Generate the name the file to write the status do and get
*         an available unit number to open the file with status append.
          call rep_open( prog_name, '.status', rep_file, rep_unit)
 
*         Now write the status message
          if( len_fil.le.1 ) then
              write(rep_unit, 120, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                    routine(1:len_sub),
     .            message(1:len_mes)
              write(*, 120, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes)
 
 120          format('STATUS :',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .            a,'/',a,': ',a)
           else
*             if file name passed print this as well)           
              write(rep_unit, 140, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil)
              write(*, 140, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil)
 
 140          format('STATUS :',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .            a,'/',a,': ',a,' (Name ',a,')')
           end if
           
****  Warning messages that are not fatal (can be for information or to
*     report problems with file access and/or string reading.
      else if ( type(1:4).eq.'WARN' ) then
 
*         Generate file name and open
          call rep_open( prog_name, '.warning', rep_file, rep_unit)
 
*         If the error flag is zero, assume the file name is not
*         needed
          if( ierr.eq. 0 ) then

*             If the file name is zero or 1 charcater long, just
*             print the message, otherwize output the file name
*             as well.
              if( len_fil.le.1 ) then
                  write(rep_unit, 220, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
                  write(*, 220, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
 
 220              format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .                a,'/',a,': ',a)
               else
*                 Print the file name as well.
                  write(rep_unit, 230, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
                  write(*, 230, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
 
 230              format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .                a,'/',a,': ',a,' (Name ',a,')')
              end if
 
          else
*             Report the messsage and the file name and the error number
              write(rep_unit, 240, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
              write(*, 240, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
 240          format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .            a,'/',a,': ',a,1x,a,' ERROR ',i5)
 
          end if
 
***** FATAL error messages:  These are processed the same as warnings
*     but will cause the program to stop running.
      else if ( type(1:4).eq.'FATA' ) then
 
*         Generate file name and open
          call rep_open( prog_name, '.fatal', rep_file, rep_unit)
 
*         If the error flag is zero, assume the file name is not
*         needed
          if( ierr.eq. 0 ) then

*             Again check to see if the file name has been passed
              if( len_fil.le.1 ) then
                  write(rep_unit, 320, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
                  write(*, 320, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
 
 320              format('FATAL  :',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .                a,'/',a,': ',a)
              else

*                 Print the file name as well
                  write(rep_unit, 330, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
                  write(*, 330, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
 
 330              format('FATAL  :',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .                a,'/',a,': ',a,' (Name ',a,')')
              end if         
          else
*             Report the messsage and the file name and the error number
              write(rep_unit, 340, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
              write(*, 340, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
 340          format('FATAL  :',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .            a,'/',a,': ',a,1x,a,' ERROR ',i5)  
*             if the error is a 101, it's like this is a g77 compiler with a 99-unit-number limit
              if( ierr.eq.101 ) then   
                write(rep_unit,341,iostat=jerr) 
                write(*       ,341,iostat=jerr)
 341            format('FATAL  : Error 101 suggests you are using'
     .             ,' a g77 compiler with Fortran unit numbers limited'
     .                ,' to 99; see GAMIT installation README')
              endif
          end if
 
*         For the fatal errors we close the rep_file and stop the program
*         running
          close( rep_unit, iostat=cerr )
          stop 'FATAL Error: Stop from report_stat'
 
***** See if we clearing the files
      else if( type(1:4).eq.'CLEA' ) then

*         Open each of the types of files and delete on closing
          if( len_mod.gt.1 ) then
              call rep_open( prog_name, '.status', rep_file, rep_unit)
              close(rep_unit, iostat=cerr, status='delete')
              call rep_open( prog_name, '.warning', rep_file, rep_unit)
              close(rep_unit, iostat=cerr, status='delete')
              call rep_open( prog_name, '.fatal', rep_file, rep_unit)
              close(rep_unit, iostat=cerr, status='delete')
          end if
      else
 
*         An unknown type of report has been called.  First warn the user
*         that an unknown type has been attempted, and then print out the
*         message information as if it were a warning.
*         Generate file name and open
          call rep_open( prog_name, '.warning', rep_file, rep_unit)
          write(rep_unit, 420, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod), routine(1:len_sub),
     .                type(1:len_typ)
          write(*, 420, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod), routine(1:len_sub),
     .                type(1:len_typ)
 420      format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x,
     .            a,'/',a,': Unknown type of status. Error in call to ',
     .           'report_stat -- ',a)
 
*         Now print the standard warning format.
          if( ierr.eq. 0 ) then

*             See if the file name has been passed.
              if( len_fil.le.1 ) then
                  write(rep_unit, 220, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
                  write(*, 220, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes)
               else
*                 Print the file name as well.
                  write(rep_unit, 230, iostat=jerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
                  write(*, 230, iostat=kerr) rep_date, rep_sec,
     .                    module(1:len_mod),
     .                    routine(1:len_sub),
     .                    message(1:len_mes),
     .                    file(1:len_fil)
              end if
          else
*             Report the messsage and the file name and the error number
              write(rep_unit, 240, iostat=jerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
              write(*, 240, iostat=kerr) rep_date, rep_sec,
     .                module(1:len_mod),
     .                routine(1:len_sub),
     .                message(1:len_mes),
     .                file(1:len_fil),  ierr
          end if
 
      end if
 
****  Now clean up before exiting. Close the report file so that it will
*     be updated and check if any errors occurred during the writing of
*     the messages
 
      close( rep_unit, iostat=cerr )
      if( jerr.ne.0 ) then
          write(*,520) jerr, rep_file(1:max(1,trimlen(rep_file)))
 520      format('WARNING: IOSTAT error ',i4,' occurred in ',
     .        'report_stat writing the status file ', a)
      end if
      if( kerr.ne.0 ) then
          write(*,540) kerr
 540      format('WARNING: IOSTAT error ',i4,' occurred in ',
     .            'report_stat writing the status report to stdout')
      end if
 
      if( cerr.ne.0 ) then
          write(*,560) cerr, rep_file(1:max(1,trimlen(rep_file)))
 560      format('WARNING: IOSTAT error ',i4,' occurred in ',
     .            'report_stat closing the status file ', a)
      end if
 
****  Thats all
      return
      end
 
CTITLE REP_OPEN
 
      subroutine rep_open( module, extent, rep_file, rep_unit )
 
*     This routine will open the status/warning/fatal file with access
*     append.  The unit number is searched to find an available unit
*     and if no unit can be found the unit is returned as 6 so that
*     stdout will be used.
 
* PASSED VARIABLES
 
*   rep_unit    - Unit number for the file.  Returned by this
*               - routine
 
      integer*4 rep_unit
 
*   module      - Name of the module calling report_stat
*   extent      - extent to be given to file name (status, warning
*               - or fatal).
*   rep_file        - Name given to output file (returned).
 
      character*(*) module, extent, rep_file
 
* LOCAL VARIABLES
 
*   ierr        - IOSTAT error when file opened.
*   trimlen     - Length of non-blank portion of string
 
      integer*4 ierr, trimlen
 
*   unit_open   - Set true by inquire if unit is already in use.
 
      logical unit_open

      character*8 access

 
****  First generate the file name
 
      rep_file = module(1:max(1,trimlen(module))) // extent
 
*     Now find an available unit.  Units 60-100 are searched
      unit_open = .true.
      rep_unit  =    60
      do while ( unit_open .and. rep_unit.lt.100 )
          inquire( unit=rep_unit, opened=unit_open)
 
*         If the unit is in use, try the next one
          if( unit_open ) rep_unit = rep_unit + 1
      end do
 
***   If the unit is still open, then we have a major problem.   We can
*     find a unit number for output so warning user and set the unit
*     to 6
      if( unit_open ) then
          write(*,120) extent(2:)
 120      format('WARNING: Cannot find available unit number for ',a,
     .            ' reports.  Output to stdout')
          rep_unit = 6
      else
 
*         Found unit number, so open file
          access = 'append'
          open(rep_unit, file=rep_file, status='unknown', 
     .         access=access, iostat=ierr)
          if( ierr.ne.0 ) then
 
*             Error opening the output file.  Tell user of problem
              write(*,140) ierr, extent(2:),
     .                rep_file(1:max(1,trimlen(rep_file))),
     .                rep_unit
 140          format('WARNING: IOSTAT error ',i4,' occurred opening ',
     .               a,' file ',a,' on unit ',i3)
              rep_unit = 6
          end if
      end if
 
***** Thats all
      return
      end
 


CTITLE 'dvdot'
 
      subroutine dvdot(r, v1, inc1, v2,inc2, iter)
c
c     Routine for a double precision dot product.  Result returned
c     in r

c
*   inc1,inc2       - increment on V1, and v2
*   i               - loop counter
*   i1,i2           - position in V1 and v2
*   iter            - number of values to be searched

      integer*4 inc1,inc2, i, i1,i2, iter
*
c
c   v1(inc1,1),v2(inc2,1)     ! vectors to be dotted 
c   r   -  result of dot product
      real*8 v1(1), v2(1), r    
*
c
c
****  Loop swapping values
      i1 = 1
      i2 = 1
      r = 0.d0
      do i = 1, iter
         r = r + v1(i1)*v2(i2)
         i1 = i1 + inc1
         i2 = i2 + inc2
      end do
c
***** Thats all
      return
      END
 
ctitle systime
 
      subroutine systime( date, secs)
 
*     This routine will return the current system date and seconds.
*     The date is a five element array with year, month, day, hour, minute.
*
*   Note: time is returned to nearest second at the moment.
 
* Variables
 
*   date(5)     - Year, month, day, hour, min, sec
 
      integer*4 date(5)
 
*   secs        - Seconds part of the time.
 
      real*8 secs
 
* Functions used:
*     IDATE       ! Sun routine to return date as day, month, year
*   ITIME       ! Returns time in hr, min, seconds
 
* LOCAL VARIABLES
 
*   jdate(3)    - Return from idate
*   jtime(3)    - Return from itime.
 
      integer*4 jdate(3), jtime(3)

c      external idate
 
****  Get the date
      call idate(jdate)
      date(1) = jdate(3)
      date(2) = jdate(2)
      date(3) = jdate(1)

cifc If you are using the intel fortran compiler prior to version 8 use the idate calls below      
cifc      call idate(jdate(1),jdate(2),jdate(3))
cifc      date(1) = jdate(3)
cifc      date(2) = jdate(1)
cifc      date(3) = jdate(2)
 
****  Get the hours, min and seconds
      call itime(jtime)
      date(4) = jtime(1)
      date(5) = jtime(2)
 
      secs = jtime(3)
 
****  Thats all
      return
      end
 

      subroutine XYZ_to_GEOD(rot_mat, site_pos, geod_pos )
 
 
*     Routine to compute the geodetic latitudes, longitude and
*     elipsoidal height from the cartesian XYZ coordinates of
*     the site.  The algorithm used is adopted from:
*     Heiskanen W A and H. Moritz, Physical Geodesy, W.H.Freeman
*     and company, San Francisco, 364, 1967. Chapter 5.3.
*
*     The routine also returns the tranformation between XYZ and
*     local NEU system.
*
*     RESTRICTION: User must supply the XYZ coordinates.  If these
*     are not available the can be computed from geod_pos using:
*
*     X = (N+h) cos(phi) cos(lambda)
*     Y = (N+h) cos(phi) sin(lambda)
*     Z = [(1-e**2)N+h]  sin(phi)
*
*     where  e**2 = 2f-f**2; f is flattening
*            N**2 = a**/[1-(e*sin(phi))**2]
*     N is East-West Radius of curvature.
*
      include 'const_param.h'
 
*   dlatdN      - Derivative of latitudes with respect to the
*               - (horizontal) North component of position
*               - This derivative should not be confused with
*               - dNdlat (same symbol N but different meanings)
*   dNdlat      - the derivative of the curvature (rad_curve)
*               - with respect to latitudes.
*   dXdlat      - Derivative of the X coordinate wrt latitudes
*   dYdlat      - Derivative of the Y coordinate wrt latitudes
*   dZdlat      - Derivative of the Z coordinate wrt latitudes
 
*   eccsq       - eccentricity squared computed from flattening
*   equ_rad     - radial distance of site from rotation axis.
*   geod_pos(3) - the geodetic coordinates.  The order of the
*               - values is:
*               - (1) - Geodetic co-latitudes (rads)
*               - (2) - Geodetic longitude (positive east) (rad)
*               - (3) - ellipsoidal height (m)
*   lat_i       - approximate latitudes used in iteration
*   lat_p       - approximate latitudes from previous iteration
*   long        - longitude of the site.
*   h_i         - approximate height used in iteration
*   h_p         - approximate height from previous iteration
*   rad_lat     - radius to be used in computing the geodetic
*               - latitudes
*   rad_curve    - radius of curvature at the site (N above)
*   rot_mat(3,3)    - Transformation matrix XYZ and NEU where
*               - N is -colatitudes, E is longitude and U is
*               - along ellipsoidal height
*   site_pos(3) - XYZ coordinates of the site (m)
 
*   tolerance   - toleranace for convergence on geodetic
*               - coordinates.
 
      real*8 dlatdN, dNdlat, dXdlat, dYdlat, dZdlat, eccsq, equ_rad,
     .    geod_pos(3), lat_i, lat_p, long, h_i, h_p, rad_lat,
     .    rad_curve, rot_mat(3,3), site_pos(3), tolerance

*   niter       - Number of iterations.  EXITS if more than 50 
*                 iterations are required.  (Removes problem if
*                 coordinates are NaN in IEEE floating point)

      integer*4 niter
 
*   converged   - Indicate values have converged.
 
      logical converged
 
      data
*                                     ! Converge to 0.1 mm
     .    tolerance  / 0.0001d0 /
 
***** Start, get the latitudes and height by iteration.
 
      equ_rad = sqrt( site_pos(1)**2 + site_pos(2)**2 )
      eccsq   = 2.d0*earth_flat - earth_flat**2
 
*                                            ! Set previous iteration values
      lat_p = atan2( site_pos(3), equ_rad)
      h_p   = 0.d0
 
      converged = .false.
      niter = 0
 
      do while ( .not. converged )
 
*         Get radius of curvature using previous latitudes estimate
          rad_curve = earth_rad /
     .               sqrt(1.d0 - eccsq*sin(lat_p)**2 )
          rad_lat  = equ_rad *
     .               ( 1.d0 - eccsq*rad_curve/(rad_curve+h_p) )
 
          lat_i = atan2( site_pos(3), rad_lat)
 
*                                          ! Use cos lat formula
          if( abs(lat_i).lt. pi/4 ) then
               h_i   = equ_rad/cos(lat_i) - rad_curve
*                                           ! Use sin lat formula
           else
               h_i   = site_pos(3)/sin(lat_i) - (1.d0-eccsq)*rad_curve
           end if
 
 
*         Check for convergence
          if( abs(h_i-h_p)              .lt. tolerance .and.
*                                                             ! Converged
     .        abs(lat_i-lat_p)*rad_curve.lt. tolerance ) then
 
              converged = .true.
          end if

*         Check for two many iterations
          niter = niter + 1
          if( niter.gt.50 ) then
              write(*,'('' XYZ_to_GEOD ERROR: Failure to converge'')')
              converged = .true.
          end if
 
*         Save the latest values
          h_p   = h_i
          lat_p = lat_i
 
*                     ! iterating for latitudes and height
      end do
 
 
***** Save the final values
      long = atan2( site_pos(2),site_pos(1) )
 
*                                     ! colatitudes
      geod_pos(1) = pi/2.d0 - lat_i
*                                     ! Add 2*pi
      if( long.lt.0 ) then
          geod_pos(2) = 2*pi + long
      else
          geod_pos(2) = long
      end if
 
      geod_pos(3) = h_i
 
***** Now do the transformation between XYZ and NEU
 
*     Compute derivative of the radius of curvature with to latitudes
* MOD TAH 980216:  Fixed error in dNdLat (derivative in error by factor
*     of two).  Long standing bug pointed out by Dan Leibach at CfA. 
      dNdlat = earth_rad*eccsq *sin(lat_i)*cos(lat_i) /
     .         sqrt( (1.d0-eccsq* sin(lat_i)**2)**3 )
 
*     Now do NORTH component  (First do derivate wrt latitudes, then compute
*     the Northing derivative)
 
      dXdlat  = (dNdlat * cos(lat_i) -
     .                (rad_curve+h_i)* sin(lat_i) )*cos(long)
      dYdlat  = (dNdlat * cos(lat_i) -
     .                (rad_curve+h_i)* sin(lat_i) )*sin(long)
      dZdlat  = (1.d0-eccsq)*dNdlat         *sin(lat_i) +
     .             ((1.d0-eccsq)*rad_curve + h_i)*cos(lat_i)
 
      dlatdN  = sqrt(dXdlat**2 + dYdlat**2 + dZdlat**2 )
 
*     Now do rotation matrix
 
      rot_mat(1,1) = dXdlat / dlatdN
      rot_mat(1,2) = dYdlat / dlatdN
      rot_mat(1,3) = dZdlat / dlatdN
 
 
*     Now do EAST component
      rot_mat(2,1) = -sin(long)
      rot_mat(2,2) =  cos(long)
      rot_mat(2,3) =  0.d0
 
*     Now do UP component
 
      rot_mat(3,1) = cos(lat_i) * cos(long)
      rot_mat(3,2) = cos(lat_i) * sin(long)
      rot_mat(3,3) = sin(lat_i)
 
***** Thats all
      return
      end
 
 
      subroutine GEOD_to_XYZ(geod_pos, site_pos )
 
 
*     Routine to convert the geodetic co-latitudes, longitude and
*     elipsoidal height to the cartesian XYZ coordinates of
*     the site.  The algorithm used is adopted from:
*     Heiskanen W A and H. Moritz, Physical Geodesy, W.H.Freeman
*     and company, San Francisco, 364, 1967. Chapter 5.3.
*
*     The routine also returns the tranformation between XYZ and
*     local NEU system.
*
*     The following formula is used for the conversion:
*     NOTE: The colatiude is passed.

*     X = (N+h) cos(phi) cos(lambda)
*     Y = (N+h) cos(phi) sin(lambda)
*     Z = [(1-e**2)N+h]  sin(phi)
*
*     where  e**2 = 2f-f**2; f is flattening
*            N**2 = a**/[1-(e*sin(phi))**2]
*     N is East-West Radius of curvature.
*
      include 'const_param.h'
 
*   geod_pos(3) - the geodetic coordinates.  The order of the
*               - values is:
*               - (1) - Geodetic co-latitudes (rads)
*               - (2) - Geodetic longitude (positive east) (rad)
*               - (3) - ellipsoidal height (m)
*   rad_curve    - radius of curvature at the site (N above)
*   site_pos(3) - XYZ coordinates of the site (m)
*   eccsq       - Eccentricity
 
      real*8  geod_pos(3),  rad_curve, site_pos(3), eccsq

***** Start, get the latitudes and height by iteration.
 
      eccsq   = 2.d0*earth_flat - earth_flat**2
 
      rad_curve = earth_rad /
     .            sqrt(1.d0 - eccsq*cos(geod_pos(1))**2 )

      site_pos(1) = (rad_curve+geod_pos(3))*sin(geod_pos(1)) *
     .               cos(geod_pos(2))
      site_pos(2) = (rad_curve+geod_pos(3))*sin(geod_pos(1)) *
     .               sin(geod_pos(2))
      site_pos(3) = ((1.d0-eccsq)*rad_curve+geod_pos(3))*
     .               cos(geod_pos(1))

 
***** Thats all
      return
      end
 

 

 
 
