C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:13:09:1999      *
C *                                                                            *
C *   REVISION    :  Removal of non-standard "<>" construct                    *
C *                  in formats                            JRH:17:05:2000      *
C *                  Minor mods. to format                 JRH:17:05:2000      *
C *                  Mod. to allow any year                JRH:07:09:2001      *
C *                  All necessary subroutines added to                        *
C *                  source code                           JRH:14:09:2001      *
C *                  Change of format for output of Z0     JRH:18:09:2001      *
C *                  Modified formats to cope with "bug" in                    *
C *                  latest g77 -- on reaching the end of                      *
C *                  a format statement, it does return to                     *
C *                  second outermost bracket BUT                              *
C *                  (incorrectly, I think) inserts a                          *
C *                  carriage return. The old statements                       *
C *                  have been commented out for retrieval                     *
C *                  later.                                JRH:01:11:2002      *
C *                  Addition of header re. software                           *
C *                  licensing                             JRH:23:12:2002      *
C *                                                                            *
C *   SOURCE      :  tda1537.f                                                 *
C *   ROUTINE NAME:  tda1537                                                   *
C *   TYPE        :  MAIN                                                      *
C *                                                                            *
C *   FUNCTION    :  Latest tidal analysis program (based on tda1154 and       *
C *                  tda1527). Major updates are:                              *
C *                                                                            *
C *                  1. It uses an input parameter file.                       *
C *                  2. It uses a "dot_tides" file to point to the directory   *
C *                     of the astronomical arguments, constituent lists and   *
C *                     harmonic constants for reference port.                 *
C *                  3. It has had potential Y2K problems fixed.               *
C *                  4. It has been changed to real*8.                         *
C *                  5. It uses SVD                                            *
C *                  6. It handles time zones properly                         *
C *                                                                            *
C *                  The location of the directory of astronomical arguments,  *
C *                  constituent lists and harmonic constants for the reference*
C *                  port is given in the file dot_tides in the local directory*
C *                  The value of the variable dot_tides should be selected    *
C *                  appropriately using the data statements supplied.         *
C *                                                                            *
C ******************************************************************************
C *                                                                            *
C *                                SOFTWARE LICENSING                          *
C *                                                                            *
C *                  Copyright (C) 2002 John Robert Hunter                     *
C *                                                                            *
C *                  This program is free software; you can redistribute       *
C *                  it and/or modify it under the terms of the GNU General    *
C *                  Public License as published by the Free Software          *
C *                  Foundation; either Version 2 of the license, or (at       *
C *                  your option) any later version.                           *
C *                                                                            *
C *                  This program is distributed in the hope that it will      *
C *                  be useful, but without any warranty; without even the     *
C *                  implied warranty of merchantability or fitness for a      *
C *                  particular purpose. See the GNU General Public License    *
C *                  for more details.                                         *
C *                                                                            *
C *                  A copy of the GNU General Public License is available     *
C *                  at http://www.gnu.org/copyleft/gpl.html or by writing     *
C *                  to the Free Software Foundation, Inc., 59 Temple Place    *
C *                  - Suite 330, Boston, MA 02111-1307, USA.                  *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      integer*4 pcon,pequ,pset,phead
C
      parameter(pcon=114,                 ! >/= Maximum no. of constituents used
     $          pequ=pcon*2+1,           ! >/= (maximum no. of constituents)*2+1
     $          pset=20,                 ! >/= Maximum number of input data sets
     $          phead=10)                   ! >/= Maximum number of header lines
C
      real*8 j1(pcon),j2(pcon)
      real*8 v1(pcon),v2(pcon)
      real*8 v0(pcon),vpv(pcon)
      real*8 j(pcon),sigma(pcon)
      real*8 hrel(pcon),grel(pcon)
C
      real*8 a(pequ,pequ),acov(pequ,pequ),x(pequ)
      real*8 usvd(pequ,pequ),vsvd(pequ,pequ),wsvd(pequ)
      real*8 e(pequ),tmpsvd(pequ)
      real*8 rmean(pset),dat(pset),res(pset),svar(pset)
      real*8 var(pset)
      real*8 sres(pset),r1(pset),r2(pset)
      real*8 h(pcon,pset),g(pcon,pset)
      real*8 b(pequ,pset),bcov(pequ,pset),c(pequ,pset)
      real*8 b1(pequ),c1(pequ)
C
      real*8 resiset
      real*8 jul_start,jul_year,year
      real*8 tz_t,tz_g,tz_out,del_tz
      real*8 t,tlast
      real*8 arg
      real*8 r,delta
      real*8 rc,rs
      real*8 w1,w2
      real*8 d2r
C
      real*4 sec
C
      integer*4 aisp(pcon),index_ref(pcon)
C
      integer*4 aispi
      integer*4 nhead,ihead
      integer*4 ncon,icon,nref,iref,isp
      integer*4 ndum
      integer*4 nequ,nequm1,iequ1,iequ2,ndat,idat,nset,iset
      integer*4 len,lenchar,lentid,lhschar
      integer*4 ilhs
      integer*4 indexi
      integer*4 iye,iye_start,iye_last,mon,idy,ihr,min
      integer*4 ic,ics,icr,is,iss,isr
      integer*4 ichar
      integer*4 job,info
      integer*4 ios,ifail
C
      logical logs,logr,lres,lyearchange,ltrep
      logical end_of_cons,end_of_file
      logical l_tz_t,l_tz_g
C
      character*120 head(phead)
      character*6 asplit,aref,atest,acon(pcon),arel(pcon)
C
      character*120 parameter_file,conslist_file,input_file
      character*120 hgout_file,resout_file,temp_file,port_ref
      character*120 conbuf,time_type,datbuf,tidal_dir
      character*80 formt
      character*40 refbuf
      character*8 dot_tides
      character*1 astar,tf
      character*1 slash
C
      data dot_tides/'.tides  '/                                          ! Unix
C     data dot_tides/'ARGS.DIR'/                                           ! DOS
C
      d2r=4.d0*datan(1.d0)/180.d0                           ! Degrees to radians
C
C     Find directory for astronomical arguments, constituent lists and 
C     harmonic constants for the reference port:
C
      open(unit=7,file=dot_tides,status='old',iostat=ios)
      call error_handler(ios,001)                              ! Error point 001
C
      read(7,8,iostat=ios) tidal_dir                        ! Required directory
    8 format(a120)
      call error_handler(ios,002)                              ! Error point 002
C
      close(unit=7)
C
      if(dot_tides.eq.'.tides  ') then                      ! Set path delimiter
        slash='/'                                                         ! Unix
      else
        slash=char(92)                                         ! (backslash) DOS
      endif
C
      lentid=lenchar(tidal_dir)
      tidal_dir=tidal_dir(1:lentid)//slash               ! Add appropriate slash
      lentid=lenchar(tidal_dir)
C
C     Read name of parameter file and open it:
C
      read (5,1,iostat=ios) parameter_file
    1 format(a120)
      call error_handler(ios,003)                              ! Error point 003
C
      open(unit=7,file=parameter_file,status='old',
     $     iostat=ios)
      call error_handler(ios,004)                              ! Error point 004
C
C     Output file of residuals?
C
      call areadfil3(7,'OPRESID',7,tf,ifail)
      call error_handler(ifail,005)                            ! Error point 005
C
      if(tf.ne.'t'.and.tf.ne.'T'.and.tf.ne.'f'.and.tf.ne.'F')
     $  call error_handler(1,006)                              ! Error point 006
      lres=(tf.eq.'t'.or.tf.eq.'T')
C
C     Filename for constituent lists:
C
      call areadfil3(7,'CONSLIST_FILE',13,conslist_file,ifail)
      call error_handler(ifail,007)                            ! Error point 007
C
      open(unit=11,file=tidal_dir(1:lentid)//conslist_file,
     $     status='old',iostat=ios)
      call error_handler(ios,008)                              ! Error point 008
C
C     Filename for input data:
C
      call areadfil3(7,'INPUT_FILE',10,input_file,ifail)
      call error_handler(ifail,009)                            ! Error point 009
C
      open(unit=13,file=input_file,status='old',
     $     iostat=ios)
      call error_handler(ios,010)                              ! Error point 010
C
C     Filename for (h,g) output data:
C
      call areadfil3(7,'HGOUT_FILE',10,hgout_file,ifail)
      call error_handler(ifail,011)                            ! Error point 011
C
      open(unit=8,file=hgout_file,status='unknown',
     $     iostat=ios)
      call error_handler(ios,012)                              ! Error point 012
C
      if(lres) then
C
C     Filename for residual output data:
C
        call areadfil3(7,'RESOUT_FILE',11,resout_file,ifail)
        call error_handler(ifail,013)                          ! Error point 013
C
        open(unit=10,file=resout_file,status='unknown',
     $       iostat=ios)
        call error_handler(ios,014)                            ! Error point 014
C
C     Filename for temporary storage (sequential access, binary):
C
        call areadfil3(7,'TEMP_FILE',9,temp_file,ifail)
        call error_handler(ifail,015)                          ! Error point 015
C
        open(unit=15,file=temp_file,status='unknown',
     $       access='sequential',form='unformatted',
     $       iostat=ios)
        call error_handler(ios,016)                            ! Error point 016
C
      endif
C
C     Header:
C
      call ireadfil2(7,'NHEAD',5,nhead,ifail)           ! Number of header lines
      call error_handler(ifail,017)                            ! Error point 017
C
      read(7,2) (head(ihead),ihead=1,nhead)
    2 format(a120)
C
C     Time zone of time series (if present):
C
      call dreadfil2(7,'TZ_T',4,tz_t,ifail)
      l_tz_t=(ifail.eq.0)
C
C     Type of time coordinate ("REP": time repeats every year;
C                              "NONREP": time is continuous from year to year):
C
      call areadfil3(7,'TIME_TYPE',9,time_type,ifail)
      call error_handler(ifail,018)                            ! Error point 018
C
      if(time_type.ne.'REP'.and.time_type.ne.'rep'.and.
     $   time_type.ne.'NONREP'.and.time_type.ne.'nonrep')
     $  call error_handler(1,019)                            ! Error point 019
      ltrep=(time_type.eq.'REP'.or.time_type.eq.'rep')
C
C     Input constituent lists:
C
      ncon=0
      nref=0
      end_of_cons=.false.
      end_of_file=.false.
C
      do while(ios.eq.0.and..not.end_of_cons.and..not.end_of_file)
C
        read(11,10,iostat=ios) conbuf
   10   format(a120)
C
        if(ios.eq.0) then                                         ! Read is O.K.
C
          if(conbuf(1:1).ne.'#') then                          ! Ignore comments
C
            len=lenchar(conbuf)
            if(len.le.6) then                        ! Only one constituent name
C
              ncon=ncon+1
              acon(ncon)=conbuf(1:6)
              call changecase(acon(ncon),acon(ncon),len,0) !Change to lower case
C
            else !Assume more than one constituent name, so related constituents
C
              end_of_cons=.true.
C
            endif
C
          endif
C
        else if(ios.gt.0) then                                           ! Error
C
          call error_handler(1,020)                            ! Error point 020
C
        else                                                ! Assume end of file
C
          end_of_file=.true.
C
        endif
C
      end do
C
      write(6,17) ncon
   17 format('No. of constituents = ',i4)
C
      nequ=2*ncon+1
      nequm1=nequ-1
C
      do icon=1,ncon                 ! Set index array for relating constutuents
        index_ref(icon)=0
      end do
C
      if(.not.end_of_file) then    ! There are some related constituents to read
C
        nref=0
        ios=0
C
        do while(ios.eq.0)
C
          if(nref.ne.0) read(11,10,iostat=ios) conbuf !Have already read 1st one
C
          if(ios.eq.0) then                                       ! Read is O.K.
C
            if(conbuf(1:1).ne.'#') then                        ! Ignore comments
C
              asplit=conbuf(1:6)
C
              ilhs=lhschar(conbuf(7:120))
              aref=conbuf(ilhs+6:ilhs+11)
C
              len=lenchar(asplit)                         ! Change to lower case
              call changecase(asplit,asplit,len,0)
              len=lenchar(aref)                           ! Change to lower case
              call changecase(aref,aref,len,0)
C
              do icon=1,ncon
                if(acon(icon).eq.asplit) isp=icon
                if(acon(icon).eq.aref) iref=icon
              end do
C
              index_ref(isp)=iref
              nref=nref+1
              aisp(nref)=isp
C
            endif
C
          else if(ios.gt.0) then                                         ! Error
C
            call error_handler(1,021)                          ! Error point 021
C
          endif
C
        end do
C
      else                                   ! There are no related constituents
C
        nref=0
C
      endif
C
      write(6,16) nref
   16 format('No. of related constituents = ',i4)
C
      l_tz_g=.false.
C
      if(nref.ne.0) then
C
C     Input reference data from port
C
        call areadfil3(7,'PORT_REF',8,port_ref,ifail)
        call error_handler(ifail,022)                          ! Error point 022
C
        len=lenchar(port_ref)
        call changecase(port_ref,port_ref,len,0)         ! Convert to lower case
C
        open(unit=7,file=tidal_dir(1:lentid)//port_ref,status='old',
     $       iostat=ios)
        call error_handler(ios,023)                            ! Error point 023
C
C     Note that constituent list may be in any order ....
C
        ndum=0
        ios=0
        do while(ios.eq.0)
          read(7,6,iostat=ios) refbuf
    6     format(a40)
          if(ios.gt.0) call error_handler(ios,024)             ! Error point 024
          if(ios.eq.0) then
            if(refbuf(1:1).eq.'#') then                              ! A comment
              if(refbuf(1:5).eq.'# TZ ') then              ! Found the time zone
                read(refbuf(6:40),*,iostat=ios) tz_g
                call error_handler(ios,025)                    ! Error point 025
                l_tz_g=.true.
              endif
            else                                            ! Harmonic constants
              if(refbuf(1:3).ne.'z0 '.and.refbuf(1:3).ne.'Z0 ') then    ! Not Z0
                len=index(refbuf,' ')-1  ! Find separator between name and (h,g)
                if(len.gt.6) len=6   ! Assume format is a6 and there is no space
                ndum=ndum+1
                arel(ndum)='      '                                 ! For safety
                write(formt,'("(a",i3,")")') len
                read(refbuf,formt,iostat=ios) arel(ndum)(1:len)
                call error_handler(ios,026)                    ! Error point 026
C
                do ichar=len+1,40                      ! Change any "*" to space
                  if(refbuf(ichar:ichar).eq.'*')
     $              refbuf(ichar:ichar)=' '
                end do
C
                read(refbuf(len+1:40),*,iostat=ios) h(ndum,1),g(ndum,1)
                call error_handler(ios,027)                    ! Error point 027
                len=lenchar(arel(ndum))
                call changecase(arel(ndum),arel(ndum),len,0)     ! Change to L.C
              endif
            endif
          endif
        end do
C
C     Set up array for related constituents:
C
        do icon=1,pcon
          hrel(icon)=0.d0
          grel(icon)=0.d0
        end do
C
        do iref=1,nref
C
          aispi=aisp(iref)
          asplit=acon(aispi)
          indexi=index_ref(aispi)
          aref=acon(indexi)
          logs=.false.
          logr=.false.
          icon=1
C
          do while((.not.logs.or..not.logr)
     $             .and.icon.le.ndum)
C
            atest=arel(icon)
C
            if(asplit.eq.atest) then
              hrel(aispi)=h(icon,1)       ! Found split const. - convert to rad.
              grel(aispi)=g(icon,1)
              logs=.true.
            endif
C
            if(aref.eq.atest) then
              hrel(indexi)=h(icon,1)       ! Found ref. const. - convert to rad.
              grel(indexi)=g(icon,1)
              logr=.true.
            endif
C
            icon=icon+1
C
          end do
C
          if(.not.logs.or..not.logr)
     $      call error_handler(ios,029)                        ! Error point 029
C
        end do
C
      endif
C
C     Sort out time zone:
C
      if(.not.l_tz_t.and..not.l_tz_g) then
        write(6,11)
   11   format('No time zone has been defined')
        call error_handler(ios,028)                            ! Error point 028
      else if(.not.l_tz_t.and.l_tz_g) then ! Assume time zone of related cons.
        tz_out=tz_g
        write(6,12) tz_g
   12   format('Time zone of input and output (',f6.2,
     $         ') assumed to be that of related constituents')
      else if(l_tz_t.and..not.l_tz_g) then     ! Assume time zone in param. file
        tz_out=tz_t
        write(6,13) tz_t
   13   format('Time zone of input and output (',f6.2,
     $         ') assumed to be that of parameter file')
      else
        tz_out=tz_t
        if(tz_t.ne.tz_g) then
          write(6,14)
   14     format('Different time zone between related constituents ',
     $         'and parameter file.')
        endif
        write(6,15) tz_t
   15   format('Time zone of output (',f6.2,
     $           ') selected to be that of parameter file')
      endif
C
C     Read number of data columns and starting year:
C
      read(13,*) nset,iye
      iye_start=iye                                     ! Remember year at start
      iye_last=iye
      call julday(iye,1,1,0,0,0.,jul_start,ifail)
      call error_handler(ifail,030)                            ! Error point 030
      jul_year=jul_start
C
C     Clear covariance matrices:
C
      do iset=1,pset
        var(iset)=0.d0
      end do
C
      do iequ1=1,pequ
C
        do iset=1,pset
          bcov(iequ1,iset)=0.d0
          b(iequ1,iset)=0.d0
        end do
C
        do iequ2=1,pequ
          acov(iequ1,iequ2)=0.d0
          a(iequ1,iequ2)=0.d0
        end do
C
      end do
C
      x(1)=1.d0
C
      ndat=0
      ios=0
C
      do while(ios.eq.0)
C
        read(13,4,iostat=ios) datbuf
    4   format(a120)
C
        if(ios.gt.0) call error_handler(ios,031)               ! Error point 031
C
        if(ios.eq.0) then
C
          read(datbuf,*,iostat=ios) t,(dat(iset),iset=1,nset)
C
          ndat=ndat+1
          if((ndat/500)*500.eq.ndat) then
            write(6,29) ndat
   29       format('Processing data point ',i5)
          endif
C
          if(ltrep) then
C
            if(ndat.eq.1) then
              lyearchange=.false.
            else
              lyearchange=(t.lt.tlast)
            endif
C
            if(lyearchange) then
              iye=iye+1                                            ! Update year
            endif
C
          else
C
C     Find year:
C
            call caldat(iye,mon,idy,ihr,min,sec,
     $                  dble(t)/86400.d0+jul_start,
     $                  ifail)
            call error_handler(ifail,032)                      ! Error point 032
C
            lyearchange=(iye.ne.iye_last)
C
          endif
C
          if(ndat.eq.1.or.lyearchange) then             ! Load astronomical data
C
C     Removed 7/9/2001:
C
C           if(iye.lt.1970.or.iye.gt.2030)
C    $        call error_handler(1,033)                        ! Error point 033
C
            iye_last=iye
C
            if(.not.ltrep) then

              call julday(iye,1,1,0,0,0.,jul_year,
     $                    ifail)
              call error_handler(ifail,034)                    ! Error point 034
C
            endif
C
            call ldast(acon,iye,ncon,sigma,j1,j2,v1,v2,v0,year,
     $                 tidal_dir,lentid,
     $                 ifail)
            call error_handler(ifail,035)                      ! Error point 035
C
          endif
C
C     Change time so that it is relative to start of current year, if necessary:
C
          if(.not.ltrep) t=t-sngl((jul_year-jul_start)*86400.d0)
C
C     Interpolate j and vpv:
C
          w2=(t-tz_out*3600.d0)/year                               ! Time in GMT
C         w2=t/year                                                 ! Local time
          w1=1.d0-w2
          do icon=1,ncon
            j(icon)=j1(icon)*w1+j2(icon)*w2
            vpv(icon)=v0(icon)+v1(icon)*w1+v2(icon)*w2
          end do
C
          do icon=1,ncon
            ic=2*icon
            is=ic+1
            arg=vpv(icon)+sigma(icon)*t
            x(ic)=j(icon)*dcos(arg*d2r)
            x(is)=j(icon)*dsin(arg*d2r)
          end do
C
          do iset=1,nset
            var(iset)=var(iset)+dat(iset)**2
          end do
C
          do iequ1=1,nequ
C
            do iequ2=iequ1,nequ
              acov(iequ1,iequ2)=acov(iequ1,iequ2)+x(iequ1)*x(iequ2)
            end do
C
            do iset=1,nset
              bcov(iequ1,iset)=bcov(iequ1,iset)+x(iequ1)*dat(iset)
            end do
          end do
C
C     Output X for temp. storage:
C
          if(lres) write(15)(x(iequ1),iequ1=1,nequ)
C
          tlast=t
C
        endif
C
      end do
C
      write(6,5) ndat
    5 format(/'No. of data points = ',i5/)
C
      do iequ1=1,nequm1                            ! Generate other half of ACOV
        do iequ2=iequ1+1,nequ
          acov(iequ2,iequ1)=acov(iequ1,iequ2)
        end do
      end do
C
C     Divide VAR, ACOV and BCOV by number of data points:
C
      do iset=1,nset
        var(iset)=var(iset)/dble(ndat)
      end do
C
      do iequ1=1,nequ
C
        do iequ2=1,nequ
          acov(iequ1,iequ2)=acov(iequ1,iequ2)/dble(ndat)
        end do
C
        do iset=1,nset
          bcov(iequ1,iset)=bcov(iequ1,iset)/dble(ndat)
        end do
C
      end do
C
C     Generate A and B:
C
      do iequ1=1,nequ
C
        logs=.false.
        if(iequ1.ne.1.and.nref.ne.0) then
          icon=iequ1/2                              ! Test for split constituent
          if(index_ref(icon).gt.0) logs=.true.
        endif
C
        do iequ2=1,nequ
          if(.not.logs) a(iequ1,iequ2)=acov(iequ1,iequ2)
        end do
C
        do iset=1,nset
          if(.not.logs) b(iequ1,iset)=bcov(iequ1,iset)
        end do
C
      end do
C
      if(nref.ne.0.) then
C
C     Deal with split constutuents:
C
        do icon=1,ncon
C
          if(index_ref(icon).ne.0) then
C
            r=hrel(icon)/hrel(index_ref(icon))
            delta=grel(icon)-grel(index_ref(icon))
C
C     Correct reference phases to final time zone:
C
            if(l_tz_g) then
              del_tz=(tz_out-tz_g)*3600.d0         ! Change in time zone (secs.)
              delta=delta+(sigma(icon)-sigma(index_ref(icon)))*del_tz
            endif
C
            ics=2*icon
            iss=ics+1
            icr=2*index_ref(icon)
            isr=icr+1
C
            rc=r*cos(delta*d2r)
            rs=r*sin(delta*d2r)
C
            a(ics,ics)=1.d0
            a(ics,icr)=-rc
            a(ics,isr)=rs
            a(iss,iss)=1.d0
            a(iss,icr)=-rs
            a(iss,isr)=-rc
C
          endif
C
        end do
C
      endif
C
C     Solve A*C=B, using SVD:
C
      job=11
C
      call dsvdc(a,pequ,nequ,nequ,
     $           wsvd,e,usvd,pequ,vsvd,pequ,
     $           tmpsvd,job,info)
C
      write(6,7) info
    7 format('Info from LINPACK dsvdc: ',i5)
C
      call error_handler(info,036)                             ! Error point 036
C
      write(6,3) (wsvd(iequ1),iequ1=1,nequ)
    3 format('Singular values : '/
     $       (5x,5e12.4))
C
C     Back substitute:
C
      do iset=1,nset
C
        do iequ1=1,nequ
          b1(iequ1)=b(iequ1,iset)
        end do
C
        call dsvbksb2(usvd,wsvd,vsvd,nequ,nequ,nequ,
     $                pequ,pequ,b1,c1,tmpsvd)
C
        do iequ1=1,nequ
          c(iequ1,iset)=c1(iequ1)
        end do
C
      end do
C
      do iset=1,nset
C
        rmean(iset)=c(1,iset)
C
        do icon=1,ncon
C
          ic=2*icon
          is=ic+1
          h(icon,iset)=sqrt(c(ic,iset)**2+c(is,iset)**2)
C
          if(c(is,iset).ne.0.d0.or.c(ic,iset).ne.0.d0) then
            arg=atan2(c(is,iset),c(ic,iset))/d2r
            if(arg.lt.0.d0) arg=arg+360.d0
          else
            arg=0.d0
          endif
C
          g(icon,iset)=arg
C
        end do
C
      end do
C
C     Output results:
C
      do ihead=1,nhead
        len=lenchar(head(ihead))
        write(formt,'("(a2,a",i3,")")') len
        write(8,formt) '# ',head(ihead)(1:len)
      end do
      write(8,30) '# TZ  ',tz_out,
     $            '#        ',(('   H      G   '),iset=1,nset)
C  30 format(a6,f6.2/
C    $       a9,(a14))
   30 format(a6,f6.2/
     $       a9,100(a14))
C
      write(8,19)(rmean(iset),iset=1,nset)
C  19 format('Z0      ',(f7.4,7x))
   19 format('Z0      ',100(f7.4,7x))
C
      do icon=1,ncon
        astar=' '
        if(index_ref(icon).gt.0) astar='*'
        len=lenchar(acon(icon))                          ! Convert to upper case
        call changecase(acon(icon),acon(icon),len,1)
        write(8,20) acon(icon),astar,(h(icon,iset),g(icon,iset),
     $              iset=1,nset)
C  20   format(a6,1x,a1,(1x,f5.3,2x,f5.1,1x))
   20   format(a6,1x,a1,100(1x,f5.3,2x,f5.1,1x))
      end do
C
C     Compute residual sum of squares;
C
      do iset=1,nset
C
        resiset=var(iset)
C
        do iequ1=1,nequ
          resiset=resiset-2.*c(iequ1,iset)*bcov(iequ1,iset)
          do iequ2=1,nequ
            resiset=resiset
     $               +c(iequ1,iset)*c(iequ2,iset)*acov(iequ1,iequ2)
          end do
        end do
C
        res(iset)=resiset
C
      end do
C
C     Compute statistics:
C
      do iset=1,nset
C
        var(iset)=var(iset)-bcov(1,iset)**2
        svar(iset)=sqrt(var(iset))
        sres(iset)=sqrt(abs(res(iset)))
        r1(iset)=res(iset)/var(iset)
        r2(iset)=sres(iset)/svar(iset)
C
      end do
C
      write(6,21)(iset,iset=1,nset)
C  21 format(//'STATISTICS'//'DATA SET NO :',19x,(i3,9x)/)
   21 format(//'STATISTICS'//'DATA SET NO :',19x,100(i3,9x)/)
C
      write(6,22)(var(iset),iset=1,nset)
C  22 format(/'VARIANCE ORIGINAL (VO)   ',(f12.3))
   22 format(/'VARIANCE ORIGINAL (VO)   ',100(f12.3))
C
      write(6,23)(res(iset),iset=1,nset)
C  23 format('VARIANCE RESIDUAL (VR)   ',(f12.5))
   23 format('VARIANCE RESIDUAL (VR)   ',100(f12.5))
C
      write(6,24)(r1(iset),iset=1,nset)
C  24 format('RATIO VR/VO',14x,(f12.4))
   24 format('RATIO VR/VO',14x,100(f12.4))
C
      write(6,25)(svar(iset),iset=1,nset)
C  25 format('SQRT(VO)',17x,(f12.3))
   25 format('SQRT(VO)',17x,100(f12.3))
C
      write(6,26)(sres(iset),iset=1,nset)
C  26 format('SQRT(VR)',17x,(f12.5))
   26 format('SQRT(VR)',17x,100(f12.5))
C
      write(6,27)(r2(iset),iset=1,nset)
C  27 format('RATIO SQRT(VR)/SQRT(VO)  ',(f12.4))
   27 format('RATIO SQRT(VR)/SQRT(VO)  ',100(f12.4))
C
      write(6,28)
   28 format(//)
C
C     Compute residuals:
C
      if(lres) then
C
        rewind 13
        rewind 15
        read(13,*)
        write(10,32) nset,iye_start
   32   format(1x,2i5)
C
        do idat=1,ndat
C
          read(13,*) t,(dat(iset),iset=1,nset)
          read(15) (x(iequ1),iequ1=1,nequ)
C
          do iset=1,nset
            do iequ1=1,nequ
              dat(iset)=dat(iset)-c(iequ1,iset)*x(iequ1)
            end do
          end do
C
          write(10,31) t,(dat(iset),iset=1,nset)
C  31     format(f10.0,(f10.3))
   31     format(f10.0,100(f10.3))
C
        end do
C
      endif
C
      stop
C
      end
C
      subroutine areadfil3(nin,anot,nchar,avar,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:08:09:1999      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1999      *
C *                                                                            *
C *   SOURCE      :  tda1537.f                                                 *
C *   ROUTINE NAME:  areadfil3                                                 *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Reads CHARACTER*120 variable from file based on keyword.  *
C *                  (Modified version of areadfil, without output to operator *
C *                   or log file)                                             *
C *                                                                            *
C *                  nin ..... File input device                               *
C *                  anot .... Keyword in file                                 *
C *                  nchar ... Number of characters in ANOT (max. 40)          *
C *                  avar .... Resultant CHARACTER*120 variable                *
C *                  ifail ... 0 for successful execution, otherwise 1         *
C *                                                                            *
C *                  (A modified version of areadfil2 for 120 characters.)     *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      integer*4 nin,nchar,ifail
      character*(*) anot
      character*120 avar
C
      integer*4 ios
      logical found
      character*161 buff
C
      rewind(nin)
      ios=0
      found=.false.
      do while(ios.eq.0.and..not.found)
        read(nin,1,iostat=ios) buff
    1   format(a160)
        if(buff(1:nchar).eq.anot.and.
     $     buff(nchar+1:nchar+1).eq.' ') then             ! Match has been found
          avar(1:120)=buff(nchar+2:nchar+121)
          found=.true.
        endif
      end do
      if(ios.eq.0) then
        ifail=0                                       ! Match found and no error
      else
        ifail=1                                                 ! No match found
      endif
      end
C
      subroutine ldast(acon,iye,ncon,sigma,j1,j2,v1,v2,v0,year,
     $                 tidal_dir,lentid,
     $                 ifail)
C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:10:09:1999      *
C *                                                                            *
C *   REVISION    :  Mod. to allow any year and both "98"                      *
C *                  and "1998" type formats               JRH:07:09:2001      *
C *                                                                            *
C *   SOURCE      :  tda1537.f                                                 *
C *   ROUTINE NAME:  ldast                                                     *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Loads astronomical data from disc (converts to radians    *
C *                  and radians/sec).                                         *
C *                                                                            *
C *                  NOTE this uses unit 9 for input.                          *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      integer pcon,ptry
C
      parameter(pcon=114,                    ! >/= Max. no. of constituents used
     $          ptry=3)                 ! No. of tries at different const. names
C
      real*8 sigma(pcon)
      real*8 j1(pcon),j2(pcon),v1(pcon),v2(pcon),v0(pcon)
      real*8 year
C
      integer*4 iye,ncon,lentid,ifail
C
      character*120 tidal_dir
      character*80 formt
      character*6 acon(pcon)
C
      real*8 jul_this,jul_next
      real*8 d2r
      real*8 vpv
      real*8 dum
C
      integer*4 icon,itry,nyear,iyear
      integer*4 len,lenchar
      integer*4 idum
      integer*4 ios
C
      logical first
      logical match,found
C
      character*6 atry(ptry,2)
C
      character*40 buf
      character*6 filnam
C
      data first/.true./
C
      data atry(1,1),atry(1,2)/'sig1  ','sigma1'/
      data atry(2,1),atry(2,2)/'the1  ','theta1'/
      data atry(3,1),atry(3,2)/'lam2  ','lamda2'/
C
      save first
C
      if(first) then
C
        d2r=datan(1.d0)*4.d0/180.d0                         ! Degrees to radians
C
      endif
C
      call julday(iye,1,1,0,0,0.,jul_this,ifail)
C
      if(ifail.ne.0) then
        ifail=1
        return
      endif
C
      call julday(iye+1,1,1,0,0,0.,jul_next,ifail)
C
      if(ifail.ne.0) then
        ifail=2
        return
      endif
C
      year=(jul_next-jul_this)*86400.d0      ! Length of current year in seconds
C
      write(6,4)
    4 format(/)
C
      do icon=1,ncon
C
        ios=0
        itry=0
        found=.false.
        do while(.not.found)
          len=lenchar(acon(icon))
          filnam(1:len)=acon(icon)(1:len)
          write(formt,'("(a27,a",i3,")")') len
          write(6,formt) 'Opening ast. arg. file for ',filnam(1:len)! Diagnostic
C
          open(unit=9,
     $         file=tidal_dir(1:lentid)//filnam(1:len)//'.con',
     $         status='old',iostat=ios)
C
          found=(ios.eq.0)                                 ! .true. if open O.K.
          match=(ios.eq.0)                   ! If open not O.K., try other names
C
          do while(.not.match.and.itry.lt.ptry)
            itry=itry+1
            if(acon(icon).eq.atry(itry,1)) then
              acon(icon)=atry(itry,2)
              match=.true.
            else if(acon(icon).eq.atry(itry,2)) then
              acon(icon)=atry(itry,1)
              match=.true.
            endif
          end do
C
          if(.not.found.and..not.match) then
            ifail=3
            return
          endif
        end do
C
        read(9,*,iostat=ios) dum
C
        if(ios.ne.0) then
          ifail=4
          return
        endif
C
        if(first) sigma(icon)=dum/3600.d0        ! Convert to degrees per second
C
        nyear=iye-1900                 ! Note this format only works around 2000
        if(nyear.ge.100) nyear=nyear-100
        iyear=9999
        do while(iyear.ne.nyear.and.iyear.ne.iye)                    ! &7/9/2001
C
          read(9,1,iostat=ios) buf
    1     format(a40)
C
          if(ios.ne.0) then
            ifail=5
            return
          endif
C
          read(buf,*,iostat=ios) iyear
          if(ios.ne.0) then
            ifail=6
            return
          endif
C
          if(iyear.eq.nyear.or.iyear.eq.iye) then                    ! &7/9/2001
C
            read(buf,*,iostat=ios) idum,j1(icon),v1(icon),vpv
            if(ios.ne.0) then
              ifail=7
              return
            endif
C
            v0(icon)=vpv-v1(icon)
C
            read(9,*,iostat=ios) idum,j2(icon),v2(icon)
            if(ios.ne.0) then
              ifail=8
              return
            endif
C
          endif
C
        end do
C
        close(unit=9)
C
      end do
C
      first=.false.
C
      return
C
      end
C
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MODELLING                         REF:JRH:10:01:1994      *
C *                                                                            *
C *   REVISION    :  All "1" array declarations changed to                     *
C *                  "*" (so that subscript check works)   JRH:10:01:1994      *
C *                                                                            *
C *   SOURCE      :  tad1537.f                                                 *
C *   ROUTINE NAME:  --                                                        *
C *   TYPE        :  SUBROUTINES                                               *
C *                                                                            *
C *   FUNCTION    :  Following are :                                           *
C *                                                                            *
C *                  LINPACK routines dsvdc and drot                           *
C *                  BLAS routines daxpy,ddot,dscal,dswap,dnrm2,drotg          *
C *                                                                            *
C ******************************************************************************
C
      subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
      integer ldx,n,p,ldu,ldv,job,info
      double precision x(ldx,*),s(*),e(*),u(ldu,*),v(ldv,*),work(*)
c
c
c     dsvdc is a subroutine to reduce a double precision nxp matrix x
c     by orthogonal transformations u and v to diagonal form.  the
c     diagonal elements s(i) are the singular values of x.  the
c     columns of u are the corresponding left singular vectors,
c     and the columns of v the right singular vectors.
c
c     on entry
c
c         x         double precision(ldx,p), where ldx.ge.n.
c                   x contains the matrix whose singular value
c                   decomposition is to be computed.  x is
c                   destroyed by dsvdc.
c
c         ldx       integer.
c                   ldx is the leading dimension of the array x.
c
c         n         integer.
c                   n is the number of rows of the matrix x.
c
c         p         integer.
c                   p is the number of columns of the matrix x.
c
c         ldu       integer.
c                   ldu is the leading dimension of the array u.
c                   (see below).
c
c         ldv       integer.
c                   ldv is the leading dimension of the array v.
c                   (see below).
c
c         work      double precision(n).
c                   work is a scratch array.
c
c         job       integer.
c                   job controls the computation of the singular
c                   vectors.  it has the decimal expansion ab
c                   with the following meaning
c
c                        a.eq.0    do not compute the left singular
c                                  vectors.
c                        a.eq.1    return the n left singular vectors
c                                  in u.
c                        a.ge.2    return the first min(n,p) singular
c                                  vectors in u.
c                        b.eq.0    do not compute the right singular
c                                  vectors.
c                        b.eq.1    return the right singular vectors
c                                  in v.
c
c     on return
c
c         s         double precision(mm), where mm=min(n+1,p).
c                   the first min(n,p) entries of s contain the
c                   singular values of x arranged in descending
c                   order of magnitude.
c
c         e         double precision(p),
c                   e ordinarily contains zeros.  however see the
c                   discussion of info for exceptions.
c
c         u         double precision(ldu,k), where ldu.ge.n.  if
c                                   joba.eq.1 then k.eq.n, if joba.ge.2
c                                   then k.eq.min(n,p).
c                   u contains the matrix of left singular vectors.
c                   u is not referenced if joba.eq.0.  if n.le.p
c                   or if joba.eq.2, then u may be identified with x
c                   in the subroutine call.
c
c         v         double precision(ldv,p), where ldv.ge.p.
c                   v contains the matrix of right singular vectors.
c                   v is not referenced if job.eq.0.  if p.le.n,
c                   then v may be identified with x in the
c                   subroutine call.
c
c         info      integer.
c                   the singular values (and their corresponding
c                   singular vectors) s(info+1),s(info+2),...,s(m)
c                   are correct (here m=min(n,p)).  thus if
c                   info.eq.0, all the singular values and their
c                   vectors are correct.  in any event, the matrix
c                   b = trans(u)*x*v is the bidiagonal matrix
c                   with the elements of s on its diagonal and the
c                   elements of e on its super-diagonal (trans(u)
c                   is the transpose of u).  thus the singular
c                   values of x and b are the same.
c
c     linpack. this version dated 08/14/78 .
c              correction made to shift 2/84.
c     g.w. stewart, university of maryland, argonne national lab.
c
c     dsvdc uses the following functions and subprograms.
c
c     external drot
c     blas daxpy,ddot,dscal,dswap,dnrm2,drotg
c     fortran dabs,dmax1,max0,min0,mod,dsqrt
c
c     internal variables
c
      integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,
     *        mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1
      double precision ddot,t
      double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn,
     *                 smm1,t1,test,ztest
      logical wantu,wantv
c
c
c     set the maximum number of iterations.
c
      maxit = 30
c
c     determine what is to be computed.
c
      wantu = .false.
      wantv = .false.
      jobu = mod(job,100)/10
      ncu = n
      if (jobu .gt. 1) ncu = min0(n,p)
      if (jobu .ne. 0) wantu = .true.
      if (mod(job,10) .ne. 0) wantv = .true.
c
c     reduce x to bidiagonal form, storing the diagonal elements
c     in s and the super-diagonal elements in e.
c
      info = 0
      nct = min0(n-1,p)
      nrt = max0(0,min0(p-2,n))
      lu = max0(nct,nrt)
      if (lu .lt. 1) go to 170
      do 160 l = 1, lu
         lp1 = l + 1
         if (l .gt. nct) go to 20
c
c           compute the transformation for the l-th column and
c           place the l-th diagonal in s(l).
c
            s(l) = dnrm2(n-l+1,x(l,l),1)
            if (s(l) .eq. 0.0d0) go to 10
               if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l))
               call dscal(n-l+1,1.0d0/s(l),x(l,l),1)
               x(l,l) = 1.0d0 + x(l,l)
   10       continue
            s(l) = -s(l)
   20    continue
         if (p .lt. lp1) go to 50
         do 40 j = lp1, p
            if (l .gt. nct) go to 30
            if (s(l) .eq. 0.0d0) go to 30
c
c              apply the transformation.
c
               t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l)
               call daxpy(n-l+1,t,x(l,l),1,x(l,j),1)
   30       continue
c
c           place the l-th row of x into  e for the
c           subsequent calculation of the row transformation.
c
            e(j) = x(l,j)
   40    continue
   50    continue
         if (.not.wantu .or. l .gt. nct) go to 70
c
c           place the transformation in u for subsequent back
c           multiplication.
c
            do 60 i = l, n
               u(i,l) = x(i,l)
   60       continue
   70    continue
         if (l .gt. nrt) go to 150
c
c           compute the l-th row transformation and place the
c           l-th super-diagonal in e(l).
c
            e(l) = dnrm2(p-l,e(lp1),1)
            if (e(l) .eq. 0.0d0) go to 80
               if (e(lp1) .ne. 0.0d0) e(l) = dsign(e(l),e(lp1))
               call dscal(p-l,1.0d0/e(l),e(lp1),1)
               e(lp1) = 1.0d0 + e(lp1)
   80       continue
            e(l) = -e(l)
            if (lp1 .gt. n .or. e(l) .eq. 0.0d0) go to 120
c
c              apply the transformation.
c
               do 90 i = lp1, n
                  work(i) = 0.0d0
   90          continue
               do 100 j = lp1, p
                  call daxpy(n-l,e(j),x(lp1,j),1,work(lp1),1)
  100          continue
               do 110 j = lp1, p
                  call daxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1)
  110          continue
  120       continue
            if (.not.wantv) go to 140
c
c              place the transformation in v for subsequent
c              back multiplication.
c
               do 130 i = lp1, p
                  v(i,l) = e(i)
  130          continue
  140       continue
  150    continue
  160 continue
  170 continue
c
c     set up the final bidiagonal matrix or order m.
c
      m = min0(p,n+1)
      nctp1 = nct + 1
      nrtp1 = nrt + 1
      if (nct .lt. p) s(nctp1) = x(nctp1,nctp1)
      if (n .lt. m) s(m) = 0.0d0
      if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m)
      e(m) = 0.0d0
c
c     if required, generate u.
c
      if (.not.wantu) go to 300
         if (ncu .lt. nctp1) go to 200
         do 190 j = nctp1, ncu
            do 180 i = 1, n
               u(i,j) = 0.0d0
  180       continue
            u(j,j) = 1.0d0
  190    continue
  200    continue
         if (nct .lt. 1) go to 290
         do 280 ll = 1, nct
            l = nct - ll + 1
            if (s(l) .eq. 0.0d0) go to 250
               lp1 = l + 1
               if (ncu .lt. lp1) go to 220
               do 210 j = lp1, ncu
                  t = -ddot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l)
                  call daxpy(n-l+1,t,u(l,l),1,u(l,j),1)
  210          continue
  220          continue
               call dscal(n-l+1,-1.0d0,u(l,l),1)
               u(l,l) = 1.0d0 + u(l,l)
               lm1 = l - 1
               if (lm1 .lt. 1) go to 240
               do 230 i = 1, lm1
                  u(i,l) = 0.0d0
  230          continue
  240          continue
            go to 270
  250       continue
               do 260 i = 1, n
                  u(i,l) = 0.0d0
  260          continue
               u(l,l) = 1.0d0
  270       continue
  280    continue
  290    continue
  300 continue
c
c     if it is required, generate v.
c
      if (.not.wantv) go to 350
         do 340 ll = 1, p
            l = p - ll + 1
            lp1 = l + 1
            if (l .gt. nrt) go to 320
            if (e(l) .eq. 0.0d0) go to 320
               do 310 j = lp1, p
                  t = -ddot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l)
                  call daxpy(p-l,t,v(lp1,l),1,v(lp1,j),1)
  310          continue
  320       continue
            do 330 i = 1, p
               v(i,l) = 0.0d0
  330       continue
            v(l,l) = 1.0d0
  340    continue
  350 continue
c
c     main iteration loop for the singular values.
c
      mm = m
      iter = 0
  360 continue
c
c        quit if all the singular values have been found.
c
c     ...exit
         if (m .eq. 0) go to 620
c
c        if too many iterations have been performed, set
c        flag and return.
c
         if (iter .lt. maxit) go to 370
            info = m
c     ......exit
            go to 620
  370    continue
c
c        this section of the program inspects for
c        negligible elements in the s and e arrays.  on
c        completion the variables kase and l are set as follows.
c
c           kase = 1     if s(m) and e(l-1) are negligible and l.lt.m
c           kase = 2     if s(l) is negligible and l.lt.m
c           kase = 3     if e(l-1) is negligible, l.lt.m, and
c                        s(l), ..., s(m) are not negligible (qr step).
c           kase = 4     if e(m-1) is negligible (convergence).
c
         do 390 ll = 1, m
            l = m - ll
c        ...exit
            if (l .eq. 0) go to 400
            test = dabs(s(l)) + dabs(s(l+1))
            ztest = test + dabs(e(l))
            if (ztest .ne. test) go to 380
               e(l) = 0.0d0
c        ......exit
               go to 400
  380       continue
  390    continue
  400    continue
         if (l .ne. m - 1) go to 410
            kase = 4
         go to 480
  410    continue
            lp1 = l + 1
            mp1 = m + 1
            do 430 lls = lp1, mp1
               ls = m - lls + lp1
c           ...exit
               if (ls .eq. l) go to 440
               test = 0.0d0
               if (ls .ne. m) test = test + dabs(e(ls))
               if (ls .ne. l + 1) test = test + dabs(e(ls-1))
               ztest = test + dabs(s(ls))
               if (ztest .ne. test) go to 420
                  s(ls) = 0.0d0
c           ......exit
                  go to 440
  420          continue
  430       continue
  440       continue
            if (ls .ne. l) go to 450
               kase = 3
            go to 470
  450       continue
            if (ls .ne. m) go to 460
               kase = 1
            go to 470
  460       continue
               kase = 2
               l = ls
  470       continue
  480    continue
         l = l + 1
c
c        perform the task indicated by kase.
c
         go to (490,520,540,570), kase
c
c        deflate negligible s(m).
c
  490    continue
            mm1 = m - 1
            f = e(m-1)
            e(m-1) = 0.0d0
            do 510 kk = l, mm1
               k = mm1 - kk + l
               t1 = s(k)
               call drotg(t1,f,cs,sn)
               s(k) = t1
               if (k .eq. l) go to 500
                  f = -sn*e(k-1)
                  e(k-1) = cs*e(k-1)
  500          continue
               if (wantv) call drot(p,v(1,k),1,v(1,m),1,cs,sn)
  510       continue
         go to 610
c
c        split at negligible s(l).
c
  520    continue
            f = e(l-1)
            e(l-1) = 0.0d0
            do 530 k = l, m
               t1 = s(k)
               call drotg(t1,f,cs,sn)
               s(k) = t1
               f = -sn*e(k)
               e(k) = cs*e(k)
               if (wantu) call drot(n,u(1,k),1,u(1,l-1),1,cs,sn)
  530       continue
         go to 610
c
c        perform one qr step.
c
  540    continue
c
c           calculate the shift.
c
            scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)),
     *                    dabs(s(l)),dabs(e(l)))
            sm = s(m)/scale
            smm1 = s(m-1)/scale
            emm1 = e(m-1)/scale
            sl = s(l)/scale
            el = e(l)/scale
            b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0
            c = (sm*emm1)**2
            shift = 0.0d0
            if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550
               shift = dsqrt(b**2+c)
               if (b .lt. 0.0d0) shift = -shift
               shift = c/(b + shift)
  550       continue
            f = (sl + sm)*(sl - sm) + shift
            g = sl*el
c
c           chase zeros.
c
            mm1 = m - 1
            do 560 k = l, mm1
               call drotg(f,g,cs,sn)
               if (k .ne. l) e(k-1) = f
               f = cs*s(k) + sn*e(k)
               e(k) = cs*e(k) - sn*s(k)
               g = sn*s(k+1)
               s(k+1) = cs*s(k+1)
               if (wantv) call drot(p,v(1,k),1,v(1,k+1),1,cs,sn)
               call drotg(f,g,cs,sn)
               s(k) = f
               f = cs*e(k) + sn*s(k+1)
               s(k+1) = -sn*e(k) + cs*s(k+1)
               g = sn*e(k+1)
               e(k+1) = cs*e(k+1)
               if (wantu .and. k .lt. n)
     *            call drot(n,u(1,k),1,u(1,k+1),1,cs,sn)
  560       continue
            e(m-1) = f
            iter = iter + 1
         go to 610
c
c        convergence.
c
  570    continue
c
c           make the singular value  positive.
c
            if (s(l) .ge. 0.0d0) go to 580
               s(l) = -s(l)
               if (wantv) call dscal(p,-1.0d0,v(1,l),1)
  580       continue
c
c           order the singular value.
c
  590       if (l .eq. mm) go to 600
c           ...exit
               if (s(l) .ge. s(l+1)) go to 600
               t = s(l)
               s(l) = s(l+1)
               s(l+1) = t
               if (wantv .and. l .lt. p)
     *            call dswap(p,v(1,l),1,v(1,l+1),1)
               if (wantu .and. l .lt. n)
     *            call dswap(n,u(1,l),1,u(1,l+1),1)
               l = l + 1
            go to 590
  600       continue
            iter = 0
            m = m - 1
  610    continue
      go to 360
  620 continue
      return
      end
C
      subroutine  drot (n,dx,incx,dy,incy,c,s)
c
c     applies a plane rotation.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dy(*),dtemp,c,s
      integer i,incx,incy,ix,iy,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c       code for unequal increments or equal increments not equal
c         to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = c*dx(ix) + s*dy(iy)
        dy(iy) = c*dy(iy) - s*dx(ix)
        dx(ix) = dtemp
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c       code for both increments equal to 1
c
   20 do 30 i = 1,n
        dtemp = c*dx(i) + s*dy(i)
        dy(i) = c*dy(i) - s*dx(i)
        dx(i) = dtemp
   30 continue
      return
      end
C
      subroutine daxpy(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dy(*),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
C
      double precision function ddot(n,dx,incx,dy,incy)
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dy(*),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
C
      subroutine  dscal(n,da,dx,incx)
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c     modified 3/93 to return if incx .le. 0.
c
      double precision da,dx(*)
      integer i,incx,m,mp1,n,nincx
c
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
C
      subroutine  dswap (n,dx,incx,dy,incy)
c
c     interchanges two vectors.
c     uses unrolled loops for increments equal one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(*),dy(*),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c       code for unequal increments or equal increments not equal
c         to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dx(ix)
        dx(ix) = dy(iy)
        dy(iy) = dtemp
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c       code for both increments equal to 1
c
c
c       clean-up loop
c
   20 m = mod(n,3)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dx(i)
        dx(i) = dy(i)
        dy(i) = dtemp
   30 continue
      if( n .lt. 3 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,3
        dtemp = dx(i)
        dx(i) = dy(i)
        dy(i) = dtemp
        dtemp = dx(i + 1)
        dx(i + 1) = dy(i + 1)
        dy(i + 1) = dtemp
        dtemp = dx(i + 2)
        dx(i + 2) = dy(i + 2)
        dy(i + 2) = dtemp
   50 continue
      return
      end
C
      double precision function dnrm2 ( n, dx, incx)
      integer i, incx, ix, j, n, next
      double precision   dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one
      data   zero, one /0.0d0, 1.0d0/
c
c     euclidean norm of the n-vector stored in dx() with storage
c     increment incx .
c     if    n .le. 0 return with result = 0.
c     if n .ge. 1 then incx must be .ge. 1
c
c           c.l.lawson, 1978 jan 08
c     modified to correct failure to update ix, 1/25/92.
c     modified 3/93 to return if incx .le. 0.
c
c     four phase method     using two built-in constants that are
c     hopefully applicable to all machines.
c         cutlo = maximum of  dsqrt(u/eps)  over all known machines.
c         cuthi = minimum of  dsqrt(v)      over all known machines.
c     where
c         eps = smallest no. such that eps + 1. .gt. 1.
c         u   = smallest positive no.   (underflow limit)
c         v   = largest  no.            (overflow  limit)
c
c     brief outline of algorithm..
c
c     phase 1    scans zero components.
c     move to phase 2 when a component is nonzero and .le. cutlo
c     move to phase 3 when a component is .gt. cutlo
c     move to phase 4 when a component is .ge. cuthi/m
c     where m = n for x() real and m = 2*n for complex.
c
c     values for cutlo and cuthi..
c     from the environmental parameters listed in the imsl converter
c     document the limiting values are as follows..
c     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
c                   univac and dec at 2**(-103)
c                   thus cutlo = 2**(-51) = 4.44089e-16
c     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
c                   thus cuthi = 2**(63.5) = 1.30438e19
c     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
c                   thus cutlo = 2**(-33.5) = 8.23181d-11
c     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
c     data cutlo, cuthi / 8.232d-11,  1.304d19 /
c     data cutlo, cuthi / 4.441e-16,  1.304e19 /
      data cutlo, cuthi / 8.232d-11,  1.304d19 /
c
      if(n .gt. 0 .and. incx.gt.0) go to 10
         dnrm2  = zero
         go to 300
c
   10 assign 30 to next
      sum = zero
      i = 1
      ix = 1
c                                                 begin main loop
   20    go to next,(30, 50, 70, 110)
   30 if( dabs(dx(i)) .gt. cutlo) go to 85
      assign 50 to next
      xmax = zero
c
c                        phase 1.  sum is zero
c
   50 if( dx(i) .eq. zero) go to 200
      if( dabs(dx(i)) .gt. cutlo) go to 85
c
c                                prepare for phase 2.
      assign 70 to next
      go to 105
c
c                                prepare for phase 4.
c
  100 continue
      ix = j
      assign 110 to next
      sum = (sum / dx(i)) / dx(i)
  105 xmax = dabs(dx(i))
      go to 115
c
c                   phase 2.  sum is small.
c                             scale to avoid destructive underflow.
c
   70 if( dabs(dx(i)) .gt. cutlo ) go to 75
c
c                     common code for phases 2 and 4.
c                     in phase 4 sum is large.  scale to avoid overflow.
c
  110 if( dabs(dx(i)) .le. xmax ) go to 115
         sum = one + sum * (xmax / dx(i))**2
         xmax = dabs(dx(i))
         go to 200
c
  115 sum = sum + (dx(i)/xmax)**2
      go to 200
c
c
c                  prepare for phase 3.
c
   75 sum = (sum * xmax) * xmax
c
c
c     for real or d.p. set hitest = cuthi/n
c     for complex      set hitest = cuthi/(2*n)
c
   85 hitest = cuthi/float( n )
c
c                   phase 3.  sum is mid-range.  no scaling.
c
      do 95 j = ix,n
      if(dabs(dx(i)) .ge. hitest) go to 100
         sum = sum + dx(i)**2
         i = i + incx
   95 continue
      dnrm2 = dsqrt( sum )
      go to 300
c
  200 continue
      ix = ix + 1
      i = i + incx
      if( ix .le. n ) go to 20
c
c              end of main loop.
c
c              compute square root and adjust for scaling.
c
      dnrm2 = xmax * dsqrt(sum)
  300 continue
      return
      end
C
      subroutine drotg(da,db,c,s)
c
c     construct givens plane rotation.
c     jack dongarra, linpack, 3/11/78.
c                    modified 9/27/86.
c
      double precision da,db,c,s,roe,scale,r,z
c
      roe = db
      if( dabs(da) .gt. dabs(db) ) roe = da
      scale = dabs(da) + dabs(db)
      if( scale .ne. 0.0d0 ) go to 10
         c = 1.0d0
         s = 0.0d0
         r = 0.0d0
         go to 20
   10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2)
      r = dsign(1.0d0,roe)*r
      c = da/r
      s = db/r
   20 z = s
      if( dabs(c) .gt. 0.0d0 .and. dabs(c) .le. s ) z = 1.0d0/c
      da = r
      db = z
      return
      end
C
      subroutine dsvbksb2(u,w,v,n,p,minnp,nmax,pmax,b,x,tmp)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MODELLING                         REF:JRH:11:01:1994      *
C *                                                                            *
C *   REVISION    :  Re-ordered declarations               JRH:17:05:2000      *
C *                                                                            *
C *   SOURCE      :  tda1537                                                   *
C *   ROUTINE NAME:  DSVBKSB2                                                  *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Singular value decomposition solution                     *
C *                  (from Numerical Recipes - D.P. version of SVBKSB,         *
C *                   modified for use with LINPACK)                           *
C *                                                                            *
C *                  U ..... Input N by MINNP matrix                           *
C *                  W ..... Input MINNP element weight vector                 *
C *                  V ..... Input P by MINNP matrix                           *
C *                          ( NOT the transpose of "V")                       *
C *                  N ..... Number of rows of U                               *
C *                  P ..... Number of rows of V                               *
C *                  MINNP . MIN(N,P)                                          *
C *                  NMAX .. Physical number of rows of U                      *
C *                  PMAX .. Physical number of rows of V                      *
C *                  B ..... Input N element vector (right hand side)          *
C *                  X ..... Output P element solution vector                  *
C *                  TMP ... Workspace MINNP element vector                    *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      integer*4 n,p,minnp,nmax,pmax
      real*8 u(nmax,*),w(*),v(pmax,*),b(*),x(*),tmp(*)
C
      real*8 s
      integer*4 i,j,jj
      do j=1,minnp
        s=0.
        if(w(j).ne.0.)then
          do i=1,n
            s=s+u(i,j)*b(i)
          end do
          s=s/w(j)
        endif
        tmp(j)=s
      end do
      do j=1,p
        s=0.
        do jj=1,minnp
          s=s+v(j,jj)*tmp(jj)
        end do
        x(j)=s
      end do
      return
      end
C
      subroutine caldat(iye,mon,idy,ihr,min,sec,julian,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:23:05:1991      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1991      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FVS                                                *
C *   ROUTINE NAME:  CALDAT                                                    *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Converts Julian Day to (year, month, day, hour, minute,   *
C *                  second) (based on Numerical Recipes)                      *
C *                                                                            *
C *                  IYE ........ Year                                         *
C *                  MON ........ Month                                        *
C *                  IDY ........ Day                                          *
C *                  IHR ........ Hour                                         *
C *                  MIN ........ Minute                                       *
C *                  SEC ........ Second                                       *
C *                  JULIAN ..... Julian Day (non-negative, but may be         *
C *                               non-integer)                                 *
C *                  IFAIL ...... 0 for successful execution, otherwise 1      *
C *                                                                            *
C ******************************************************************************
      real*8 julian
      real*4 sec
      integer*4 iye,mon,idy,ihr,min,ifail
      integer*4 igreg,jalpha,ja,jb,jc,jd,je,ijul
      parameter (igreg=2299161)               ! Cross-over to Gregorian Calendar
      if(julian.lt.0.d0) then                  ! Negative Julian Day not allowed
        ifail=1
        return
      endif
      ijul=idint(julian)                                   ! Integral Julian Day
      sec=sngl((julian-dble(ijul))*86400.d0)! Seconds from beginning of Jul. Day
      if(sec.ge.43200.) then                              ! In next calendar day
        ijul=ijul+1
        sec=sec-43200.                            ! Adjust from noon to midnight
      else                                                ! In same calendar day
        sec=sec+43200.                            ! Adjust from noon to midnight
      endif
      if(sec.ge.86400.) then              ! Final check to prevent time 24:00:00
        ijul=ijul+1
        sec=sec-86400.
      endif
      min=int(sec/60.)                  ! Integral minutes from beginning of day
      sec=sec-float(min*60)                   ! Seconds from beginning of minute
      ihr=min/60                          ! Integral hours from beginning of day
      min=min-ihr*60                   ! Integral minutes from beginning of hour
      if(ijul.ge.igreg)then                  ! Correction for Gregorian Calendar
        jalpha=idint((dble(ijul-1867216)-0.25d0)/36524.25d0)
        ja=ijul+1+jalpha-idint(0.25d0*dble(jalpha))
      else                                                       ! No correction
        ja=ijul
      endif
      jb=ja+1524
      jc=idint(6680.d0+(dble(jb-2439870)-122.1d0)/365.25d0)
      jd=365*jc+idint(0.25d0*dble(jc))
      je=idint(dble(jb-jd)/30.6001d0)
      idy=jb-jd-idint(30.6001d0*dble(je))
      mon=je-1
      if(mon.gt.12)mon=mon-12
      iye=jc-4715
      if(mon.gt.2)iye=iye-1
      if(iye.le.0)iye=iye-1
      ifail=0
      return
      end
C
      subroutine changecase(ain,aout,n,itype)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:03:12:1990      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1990      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FVS                                                *
C *   ROUTINE NAME:  CHANGECASE                                                *
C *   TYPE        :  MAIN                                                      *
C *                                                                            *
C *   FUNCTION    :  Changes case of character variable                        *
C *                                                                            *
C *                  AIN ..... input variable                                  *
C *                  AOUT .... output variable                                 *
C *                  N ....... number of characters in AIN and AOUT            *
C *                  ITYPE ... 0 for change to lower case                      *
C *                            1 for change to upper case                      *
C *                                                                            *
C ******************************************************************************
      character*(*) ain,aout
      if(itype.eq.0) then                                 ! Change to lower case
        do i=1,n
          idum=ichar(ain(i:i))
          if(idum.ge.65.and.idum.le.90) then
            idum=idum+32
            aout(i:i)=char(idum)
          else
            aout(i:i)=ain(i:i)
          endif
        end do
        return
      else if(itype.eq.1) then                            ! Change to upper case
        do i=1,n
          idum=ichar(ain(i:i))
          if(idum.ge.97.and.idum.le.122) then
            idum=idum-32
            aout(i:i)=char(idum)
          else
            aout(i:i)=ain(i:i)
          endif
        end do
        return
      else                                 ! ITYPE out of range ..... do nothing
        do i=1,n
          aout(i:i)=ain(i:i)
        end do
        return
      endif
      end
C
      subroutine dreadfil2(nin,anot,nchar,var,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:27:01:1995      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1995      *
C *                                                                            *
C *   SOURCE      :  forlib.f                                                  *
C *   ROUTINE NAME:  dreadfil2                                                 *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Reads REAL*8 variable from file based on keyword.         *
C *                  (Modified version of areadfil2)                           *
C *                                                                            *
C *                  nin ..... File input device                               *
C *                  anot .... Keyword in file                                 *
C *                  nchar ... Number of characters in ANOT (max. 40)          *
C *                  var ..... Resultant REAL*8 variable                       *
C *                  ifail ... 0 for successful execution, otherwise 1         *
C *                                                                            *
C ******************************************************************************
      real*8 var
      integer*4 nin,nchar,ifail
      character*(*) anot
C
      integer*4 ios
      logical found
      character*80 buff
C
      rewind(nin)
      ios=0
      found=.false.
      do while(ios.eq.0.and..not.found)
        read(nin,1,iostat=ios) buff
    1   format(a80)
        if(buff(1:nchar).eq.anot.and.
     $     buff(nchar+1:nchar+1).eq.' ') then             ! Match has been found
          read(buff(nchar+2:80),*,iostat=ios) var
          found=.true.
        endif
      end do
      if(ios.eq.0) then
        ifail=0                                       ! Match found and no error
      else
        ifail=1                                                 ! No match found
      endif
      end
C
      subroutine error_handler(ifail,error_point)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MODELLING                         REF:JRH:06:01:1995      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1995      *
C *                                                                            *
C *   SOURCE      :  forlib.f                                                  *
C *   ROUTINE NAME:  error_handler                                             *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Handles error (ifail.ne.0)                                *
C *                                                                            *
C ******************************************************************************
      integer*4 ifail,error_point
C
      if(ifail.ne.0) then
        write(6,1) ifail,error_point
    1   format(/' ifail returned as ',i5,' at error point ',i5,
     $          ' ..... program terminated'/)
        stop
      endif
      return
      end
C
      integer*4 function idint2(x)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:23:04:1991      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1991      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FVS                                                *
C *   ROUTINE NAME:  IDINT2                                                    *
C *   TYPE        :  FUNCTION                                                  *
C *                                                                            *
C *   FUNCTION    :  As INT, but trunctates towards -(infinity)                *
C                    This is D.P. version of INT2                              *
C *                                                                            *
C ******************************************************************************
      real*8 x
      idint2=idint(x)
      if(dble(idint2).eq.x) return
      if(x.lt.0.d0) idint2=idint2-1
      return
      end 
C
      subroutine ireadfil2(nin,anot,nchar,ivar,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:06:01:1995      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1995      *
C *                                                                            *
C *   SOURCE      :  forlib.f                                                  *
C *   ROUTINE NAME:  ireadfil2                                                 *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Reads INTEGER*4 variable from file based on keyword.      *
C *                  (Modified version of ireadfil, without output to operator *
C *                   or log file)                                             *
C *                                                                            *
C *                  nin ..... File input device                               *
C *                  anot .... Keyword in file                                 *
C *                  nchar ... Number of characters in ANOT (max. 40)          *
C *                  ivar .... Resultant INTEGER*4 variable                    *
C *                  ifail ... 0 for successful execution, otherwise 1         *
C *                                                                            *
C ******************************************************************************
      integer*4 nin,nchar,ivar,ifail
      character*(*) anot
C
      integer*4 ios
      logical found
      character*80 buff
C
      rewind(nin)
      ios=0
      found=.false.
      do while(ios.eq.0.and..not.found)
        read(nin,1,iostat=ios) buff
    1   format(a80)
        if(buff(1:nchar).eq.anot.and.
     $     buff(nchar+1:nchar+1).eq.' ') then             ! Match has been found
          read(buff(nchar+2:80),*,iostat=ios) ivar
          found=.true.
        endif
      end do
      if(ios.eq.0) then
        ifail=0                                       ! Match found and no error
      else
        ifail=1                                                 ! No match found
      endif
      end
C
      subroutine julday(iye,mon,idy,ihr,min,sec,julian,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:23:05:1991      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1991      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FVS                                                *
C *   ROUTINE NAME:  JULDAY                                                    *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Converts (year, month, day, hour, minute, second) to      *
C *                  Julian Day (based on Numerical Recipes)                   *
C *                                                                            *
C *                  IYE ........ Year                                         *
C *                  MON ........ Month                                        *
C *                  IDY ........ Day                                          *
C *                  IHR ........ Hour                                         *
C *                  MIN ........ Minute                                       *
C *                  SEC ........ Second                                       *
C *                  JULIAN ..... Julian Day (non-negative, but may be         *
C *                               non-integer)                                 *
C *                  IFAIL ...... 0 for successful execution, otherwise 1      *
C *                                                                            *
C ******************************************************************************
      real*8 julian
      real*4 sec
      integer*4 iye,mon,idy,ihr,min,ifail
      integer*4 iyyy,jy,jm,igreg,ja,ijul
      integer*4 idint2
      parameter (igreg=15+31*(10+12*1582))
C                              ..... Gregorian Calendar was adopted 15 Oct. 1582
      if(iye.eq.0.or.                                    ! There is no year zero
     $   iye.lt.-4713) then                    ! Julian Day must be non-neagtive
        ifail=1
        return
      endif
      if(iye.lt.0) then
        iyyy=iye+1
      else
        iyyy=iye
      endif
      if(mon.gt.2) then
        jy=iyyy
        jm=mon+1
      else
        jy=iyyy-1
        jm=mon+13
      endif
      ijul=idint2(365.25d0*dble(jy))+idint2(30.6001d0*dble(jm))
     $       +idy+1720995
      if(idy+31*(mon+12*iyyy).ge.igreg) then
C                                    ..... Test for change to Gregorian Calendar
        ja=idint(0.01d0*dble(jy))
        ijul=ijul+2-ja+idint(0.25d0*dble(ja))
      endif
      julian=dble(ijul)
     $         +dble(ihr)/24.d0+dble(min)/1440.d0+dble(sec)/86400.d0
     $         -0.5d0                         ! Correction from midnight to noon
      if(julian.lt.0.d0) then                  ! Julian Day must be non-negative
        ifail=1
        return
      endif
      ifail=0
      return
      end
C
      integer function lenchar(c)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Variable declarations                 JRH:30:09:1996      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  LENCHAR                                                   *
C *   TYPE        :  INTEGER FUNCTION                                          *
C *                                                                            *
C *   FUNCTION    :  Returns length of CHARACTER  variable (defined by         *
C *                  removing blank characters from right-hand side).          *
C *                                                                            *
C ******************************************************************************
      character*(*) c
C
      integer*4 itot,i                                                ! @30/9/96
C
      itot=len(c)
      do i=itot,1,-1
        if(c(i:i).ne.' ') go to 1
      end do
      lenchar=0                                           ! String is all blanks
      return
    1 lenchar=i
      return
      end
C
      integer function lhschar(c)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:07:09:1989      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1989      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  LHSCHAR                                                   *
C *   TYPE        :  INTEGER FUNCTION                                          *
C *                                                                            *
C *   FUNCTION    :  Returns index of first non-blank character of CHARACTER   *
C *                  variable (defined by removing blank characters from       *
C *                  left-hand side).                                          *
C *                                                                            *
C ******************************************************************************
      character*(*) c
      itot=len(c)
      do i=1,itot
        if(c(i:i).ne.' ') go to 1
      end do
      lhschar=0                                           ! String is all blanks
      return
    1 lhschar=i
      return
      end
C

