C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Rationalise CHARACTER variables       JRH:10:01:1986      *
C *                  AISP() and AISPI declared INTEGER     JRH:10:01:1986      *
C *                  Upper to lower case in operator O/P   JRH:10:01:1986      *
C *                  Output device nos. changed            JRH:10:01:1986      *
C *                  PRTNEW.DAT renamed PORT.DAT           JRH:10:01:1986      *
C *                  DOPEN replaced by VAX version         JRH:10:01:1986      *
C *                  Mods. to OPEN                         JRH:10:01:1986      *
C *                  NAG name for single precision         JRH:10:01:1986      *
C *                  Modification to estimate of original                      *
C *                  variance                              JRH:10:01:1986      *
C *                  Data files converted to FORTRAN type                      *
C *                  with extra space at start of each                         *
C *                  record                                JRH:13:01:1986      *
C *                  DOPEN modified                        JRH:13:01:1986      *
C *                  Code for MATTID included              JRH:13:01:1986      *
C *                  Checks on array bounds included       JRH:13:01:1986      *
C *                  Statistics output also to disc file   JRH:13:01:1986      *
C *                  Division of arrays VAR, ACOV and BCOV                     *
C *                  by no. of data points to prevent                          *
C *                  overflows                             JRH:15:01:1986      *
C *                  Minor mod.                            JRH:15:01:1986      *
C *                  VAR and RESK to DOUBLE PRECISION to                       *
C *                  improve estimate of residual variance JRH:15:01:1986      *
C *                  Mod. for use of MATTID with different                     *
C *                  number of constituents without                            *
C *                  re-linking                            JRH:12:05:1989      *
C *                  Further mods. as above                JRH:15:05:1989      *
C *                  Minor mod. involving ATAN2            JRH:19:05:1989      *
C *                  Transfer to Sun plus associated mods. JRH:03:12:1990      *
C *                  Mod. for constituent names :                              *
C *                    (1) 6 character                                         *
C *                    (2) Equivalence of SIGMA1 and SIG1,                     *
C *                                       THETA1 and THE1,                     *
C *                                       LAMDA2 and LAM2                      *
C *                    (3) Possiblity of lower case names  JRH:03:12:1990      *
C *                  Inclusion of log file etc.            JRH:03:12:1990      *
C *                  Removal of subroutines already in                         *
C *                  FORLIB                                JRH:03:12:1990      *
C *                  Continuation of above mods.           JRH:04:12:1990      *
C *                  Mods. for 114 constituents and                            *
C *                  1970-2030                             JRH:04:12:1990      *
C *                  Change of astronomical argument                           *
C *                  filenames from *.dat to *.con         JRH:04:12:1990      *
C *                  Change to read data in free format    JRH:04:07:1994      *
C *                  Correction of spelling of tlst to                         *
C *                  tlast                                 JRH:05:07:1994      *
C *                  pset etc. increased to 10             JRH:14:03:1997      *
C *                  pset etc. increased to 20             JRH:01:05:1998      *
C *                  Other minor mods.                     JRH:01:05:1998      *
C *                  Correction so that year at head of                        *
C *                  residual file is O.K.                 JRH:01:05:1998      *
C *                  Mod. to output format for mean        JRH:31:08:1999      *
C *                  Minor mods. to formats                JRH:14:09:2001      *
C *                  Renamed tda1154.f                     JRH:14:09:2001      *
C *                  All necessary subroutines added to                        *
C *                  source code                           JRH:14:09:2001      *
C *                  Addition of header re. software                           *
C *                  licensing                             JRH:23:12:2002      *
C *                                                                            *
C *   SOURCE      :  tda1154.f                                                 *
C *   ROUTINE NAME:  tda1154                                                   *
C *   TYPE        :  MAIN                                                      *
C *                                                                            *
C *   FUNCTION    :  Tidal analysis program                                    *
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 ******************************************************************************
      integer pcon,pequ,pset                                          ! @3/12/90 
      parameter(pcon=114,        ! >/= Maximum no. of constituents used &4/12/90
     $          pequ=229,       ! >/= (maximum no. of constituents)*2+1 &4/12/90
     $          pset=20)       ! >/= Maximum number of input data sets # &1/5/98
      dimension j1(pcon),vpv(pcon),v1(pcon),j2(pcon),v2(pcon),v0(pcon)
      dimension j(pcon),sigma(pcon),index(pcon)                       ! &10/1/86
      dimension hrel(pcon),grel(pcon)                                 ! &10/1/86
      integer aisp(pcon),aispi                               ! @10/1/86 &4/12/90
      character*6 asplit,aref,atest,acon(pcon),arel(pcon)    ! @10/1/86 &4/12/90
      dimension a(pequ,pequ),acov(pequ,pequ),x(pequ)                  ! &4/12/90
C     dimension aa(pequ,pequ),bb(pequ,pset),wkspce(pequ)          ! NAG &4/12/90
      dimension rmean(pset),dat(pset),res(pset),svar(pset)            ! &15/1/86
      double precision var(pset)                                      ! @15/1/86
      dimension sres(pset),r1(pset),r2(pset)
      dimension h(pcon,pset),g(pcon,pset)
      dimension b(pequ,pset),bcov(pequ,pset),c(pequ,pset)             ! &4/12/90
      dimension amat(pequ,pequ),bmat(pequ,pset),cmat(pequ,pset)! MATTID @12/5/89
      double precision resk                                           ! @15/1/86
      real j1,j2,j
      logical logs,logr,lres,lyearchange                              ! &3/12/90
      logical yesno                                                   ! @3/12/90
      logical sing                                             ! MATTID @13/1/86
      character*80 fmt                                              ! @14/9/2001
      character*80 head1,head2,head3,port,port1                       ! @10/1/86
      character*14 ahg                                                ! @10/1/86
      character*1 astar                                      ! @10/1/86 &4/12/90
      data maxcon,maxequ,maxset/pcon,pequ,pset/
      data nopin,nptin,nastin,nconin,ndatin,ntemp     ! Set I/O devices &3/12/90
     $    /5    ,7    ,9     ,11    ,13    ,15   /                    ! &3/12/90
      data nopout,nhgo1,nhgo2,nrsout,nlog                             ! &3/12/90
     $    /6    ,8     ,10   ,12    ,14   /                           ! &10/1/86
      data rpt/600000./    ! Set interval (secs.) for interp. of J, VPV &3/12/90
      data ahg/'   H      G   '/
      call dskopen(nopin,nopout,nlog,nlog,0,                          ! @3/12/90
     $             'Filename for log file ?                 ',        ! @3/12/90
     $             'unknown   ')                                      ! @3/12/90
      f2=3.14159265/180.
      lres=yesno(nopin,nopout,nlog,                                   ! @3/12/90
     $           'O/P file of residuals (Y or N) ?        ')          ! @3/12/90
      open(unit=nptin,file='port.dat',status='old')                   ! &10/1/86
      call dskopen(nopin,nopout,nlog,nconin,0,                        ! @3/12/90
     $             'Filename for constituent lists ?        ',        ! @3/12/90
     $             'unknown   ')                                      ! @3/12/90
      call dskopen(nopin,nopout,nlog,ndatin,0,                        ! @3/12/90
     $             'Filename for input data ?               ',        ! @3/12/90
     $             'unknown   ')                                      ! @3/12/90
      call dskopen(nopin,nopout,nlog,nhgo1,0,                         ! @3/12/90
     $             'Filename for output (1) of H and G ?    ',        ! @3/12/90
     $             'unknown   ')                                      ! @3/12/90
      call dskopen(nopin,nopout,nlog,nhgo2,0,                         ! @3/12/90
     $             'Filename for output (2) of H and G ?    ',        ! @3/12/90
     $             'unknown   ')                                      ! @3/12/90
      if(lres) call dskopen(nopin,nopout,nlog,ntemp,1,       ! @3/12/90 &4/12/90
     $                      'Filename for temporary storage ?        ',!@3/12/90
     $                      'unknown   ') ! Unform. sequ. acc. @3/12/90 &4/12/90
      if(lres) call dskopen(nopin,nopout,nlog,nrsout,0,               ! @3/12/90
     $                      'Filename for output of residuals ?      ',!@3/12/90
     $                      'unknown   ')                             ! @3/12/90
      write(nopout,52)
   52 format(' Input header (3 lines) for O/P files :')               ! &10/1/86
      read(nopin,53) head1,head2,head3
   53 format((a80))                                                   ! &10/1/86
      write(nlog,17) head1,head2,head3                                ! @3/12/90
   17 format((1x,a80))                                     ! @3/12/90 &14/9/2001
      read(nconin,6) ncon                      ! Input constituent list &3/12/90
    6 format(1x,i5)                                        ! &10/1/86 &14/9/2001
      read(nconin,7)(acon(i),i=1,ncon)
    7 format((1x,10(a6)))                         ! &13/1/86 &4/12/90 &14/9/2001
      do i=1,ncon                                ! Change to lower case @4/12/90
        len=lenchar(acon(i))                                          ! @4/12/90
        call changecase(acon(i),acon(i),len,0)                        ! @4/12/90
      end do                                                          ! @4/12/90
      write(nopout,26) ncon
      write(nlog,26) ncon                                             ! @3/12/90
   26 format(' No. of constituents = ',i4)                            ! &10/1/86
      nequ=2*ncon+1
      nequm1=nequ-1
      do k=1,ncon           ! Set index array for relating constutuents &3/12/90
        index(k)=0                                                    ! &3/12/90
      end do                                                          ! @3/12/90
      i=1                         ! Input data for related constituents &3/12/90
   12 read(nconin,9,end=10) asplit,aref
    9 format((1x,2(a6)))                          ! &10/1/86 &4/12/90 &14/9/2001
      len=lenchar(asplit)                        ! Change to lower case @4/12/90
      call changecase(asplit,asplit,len,0)                            ! @4/12/90
      len=lenchar(aref)                          ! Change to lower case @4/12/90
      call changecase(aref,aref,len,0)                                ! @4/12/90
      do k=1,ncon                                                     ! &3/12/90
        if(acon(k).eq.asplit) isp=k
        if(acon(k).eq.aref) iref=k                                    ! &3/12/90
      end do                                                          ! @3/12/90
      index(isp)=iref
      aisp(i)=isp
      i=i+1
      go to 12
   10 nref=i-1
      write(nopout,25) nref
      write(nlog,25) nref                                             ! @3/12/90
   25 format(' No. of related constituents = ',i4)                    ! &10/1/86
      if(nref.ne.0) then                                              ! &3/12/90
C     Input reference data from port                                  ! &3/12/90
C     - note that constituent list may be in any order ....           ! &3/12/90
   34   write(nopout,13)
        write(nlog,13)                                                ! @3/12/90
   13   format(' Input reference port : ')                 ! &10/1/86 &14/9/2001
        read(nopin,14) port
   14   format(a80)                                                   ! &10/1/86
        write(nlog,11) port                                           ! @3/12/90
   11   format(1x,a80)                                     ! @3/12/90 &14/9/2001
        len=lenchar(port)                                             ! @3/12/90
        call changecase(port,port,len,0)        ! Convert to lower case @3/12/90
   22   read(nptin,95) port1                                          ! &10/1/86
   95   format(1x,a80)                                     ! &10/1/86 &14/9/2001
        len=lenchar(port1)                                            ! @3/12/90
        call changecase(port1,port1,len,0)      ! Convert to lower case @3/12/90
        if(port1(1:5).eq.' stop') then                       ! &10/1/86 &3/12/90
          write(nopout,16)
          write(nlog,16)                                              ! @3/12/90
   16     format(' Port not listed')                                  ! &10/1/86
          rewind nptin
          go to 34
        endif                                                         ! @3/12/90
   15   if(port.ne.port1) then                               ! &10/1/86 &3/12/90
   18     read(nptin,95) port1                                        ! &10/1/86
          read(nptin,95) port1                                        ! &10/1/86
          read(nptin,20) ndum
   20     format(1x,i5)                                    ! &10/1/86 &14/9/2001
          nrec=1+3*((ndum-1)/10+1)
          do i=1,nrec                                                 ! &3/12/90
            read(nptin,95) port1                             ! &10/1/86 &3/12/90
          end do                                                      ! @3/12/90
          go to 22
        endif                                                         ! @3/12/90
        read(nptin,95) port1                                          ! &10/1/86
        read(nptin,95) port1                                          ! &10/1/86
        read(nptin,20) ndum
        read(nptin,7)(arel(i),i=1,ndum)
        do i=1,ndum                             ! Convert to lower case @4/12/90
          len=lenchar(arel(i))                                        ! @4/12/90
          call changecase(arel(i),arel(i),len,0)                      ! @4/12/90
        end do                                                        ! @4/12/90
        read(nptin,95) port1                                          ! &10/1/86
        read(nptin,23)(h(i,1),i=1,ndum)
   23   format((1x,10f10.3))                               ! &10/1/86 &14/9/2001
        read(nptin,24)(g(i,1),i=1,ndum)
   24   format((1x,10f10.1))                               ! &10/1/86 &14/9/2001
        do k=1,maxcon           ! Set up array for related constituents &3/12/90
          hrel(k)=0.
          grel(k)=0.                                                  ! &3/12/90
        end do                                                        ! @3/12/90
        do i=1,nref                                                   ! &3/12/90
          aispi=aisp(i)
          asplit=acon(aispi)
          indexi=index(aispi)
          aref=acon(indexi)
          logs=.false.
          logr=.false.
          do k=1,ndum                                                 ! &3/12/90
            atest=arel(k)
            if(asplit.eq.atest) then                                  ! &3/12/90
              hrel(aispi)=h(k,1) ! Found split const. - convert to rad. &3/12/90
              grel(aispi)=g(k,1)*f2
              logs=.true.
            endif                                                     ! @3/12/90
            if(aref.eq.atest) then                                    ! &3/12/90
              hrel(indexi)=h(k,1) ! Found ref. const. - convert to rad. &3/12/90
              grel(indexi)=g(k,1)*f2
              logr=.true.
            endif                                                     ! @3/12/90
            if(logs.and.logr) go to 32
C                           ..... Jump out of loop if both const. found &3/12/90
          end do                                                      ! @3/12/90
          write(nopout,33)                      ! Constituent not found &3/12/90
          write(nlog,33)                                              ! @3/12/90
   33     format(' Required related constituent',                     ! &3/12/90
     $           ' not found in port data')                  ! &10/1/86 &3/12/90
          rewind nptin
          go to 34
   32     continue
        end do                                                        ! @3/12/90
      endif
      read(ndatin,*) nset,iye                                          ! &4/7/94
    1 format(1x,2i5)                                        ! &10/1/86 14/9/2001
      iye_start=iye                             ! Remember year at start @1/5/98
      do k=1,maxset                         ! Clear covariance matrices &3/12/90
        var(k)=0.                                                     ! &3/12/90
      end do                                                          ! @3/12/90
      do i=1,maxequ                                                   ! &3/12/90
        do k=1,maxset                                                 ! &3/12/90
          bcov(i,k)=0.
          b(i,k)=0.                                                   ! &3/12/90
        end do                                                        ! @3/12/90
        do k=1,maxequ                                                 ! &3/12/90
          acov(i,k)=0.
          a(i,k)=0.                                                   ! &3/12/90
        end do                                                        ! @3/12/90
      end do                                                          ! @3/12/90
      x(1)=1.
      i=1
      write(nopout,97)                                                ! @15/1/86
      write(nlog,97)                                                  ! @3/12/90
   97 format(/)                                                       ! @15/1/86
    4 read(ndatin,*,end=3) t,(dat(k),k=1,nset)                         ! &4/7/94
      if((i/100)*100.eq.i) then                                       ! &3/12/90
        write(nopout,93) i                                            ! &3/12/90
   93   format(' Processing data point ',i5)                 ! &15/1/86 &4/12/90
        write(nlog,19) i                                              ! @3/12/90
   19   format(' Processing data point ',i5)                          ! @3/12/90
      endif                                                           ! @3/12/90
      if(i.eq.1) then                                                 ! @3/12/90
        lyearchange=.false.                                           ! @3/12/90
      else                                                            ! @3/12/90
        lyearchange=(t.lt.tlast)                              ! @3/12/90 &5/7/94
      endif                                                           ! @3/12/90
      if(lyearchange) then                                            ! &3/12/90
        iye=iye+1                                         ! Update year &3/12/90
      endif                                                           ! @3/12/90
      if(i.eq.1.or.lyearchange) then           ! Load astronomical data @3/12/90
        if(iye.lt.1970.or.iye.gt.2030) then                           ! &4/12/90
          write(nopout,85)
          write(nlog,85)                                              ! &3/12/90
   85     format(/' Year out of range 1970-2030',                     ! &4/12/90
     $            ' ..... program terminated'/)              ! &10/1/86 &3/12/90
          stop
        endif                                                         ! @3/12/90
        call ldast(acon,iye,ncon,nastin,nopout,nlog,sigma,
     $             j1,j2,v1,v2,v0,year,iint)                 ! &3/12/90 &4/12/90
      endif                                                           ! @3/12/90
      idum=t/rpt+1.    ! Astron. data loaded - interp. J, VPV regularly &3/12/90
      if(iint.ne.idum) then                                           ! &4/12/90
        iint=idum
        ridum=idum
        tdum=(ridum-0.5)*rpt
        do k=1,ncon                                                   ! &4/12/90
          j(k)=(j1(k)*(year-tdum)+j2(k)*tdum)/year
          vpv(k)=v0(k)+(v1(k)*(year-tdum)+v2(k)*tdum)/year            ! &4/12/90
        end do                                                        ! @4/12/90
      endif                                                           ! @4/12/90
      do k=1,ncon                                                     ! &4/12/90
        ic=2*k
        is=ic+1
        arg=vpv(k)+sigma(k)*t
        x(ic)=j(k)*cos(arg)
        x(is)=j(k)*sin(arg)                                           ! &4/12/90
      end do                                                          ! @4/12/90
      do k=1,nset                                                     ! &4/12/90
        var(k)=var(k)+dat(k)**2                                       ! &4/12/90
      end do                                                          ! @4/12/90
      do k=1,nequ                                                     ! &4/12/90
        do kk=k,nequ                                                  ! &4/12/90
          acov(k,kk)=acov(k,kk)+x(k)*x(kk)                            ! &4/12/90
        end do                                                        ! @4/12/90
        do ks=1,nset                                                  ! &4/12/90
          bcov(k,ks)=bcov(k,ks)+x(k)*dat(ks)                          ! &4/12/90
        end do                                                        ! @4/12/90
      end do                                                          ! @4/12/90
      if(lres) write(ntemp)(x(k),k=1,nequ) ! Output X for temp. storage &4/12/90
      i=i+1
      tlast=t
      goto 4
    3 ndat=i-1
      rn=ndat                                                         ! @15/1/86
      do k=1,nequm1                       ! Generate other half of ACOV &4/12/90
        kp1=k+1
        do kk=kp1,nequ                                                ! &4/12/90
          acov(kk,k)=acov(k,kk)                                       ! &4/12/90
        end do                                                        ! @4/12/90
      end do                                                          ! @4/12/90
      do k=1,nset !Divide VAR, ACOV & BCOV by no. of data pts. @15/1/86 &4/12/90
        var(k)=var(k)/rn                                              ! @15/1/86
      end do                                                          ! @15/1/86
      do k=1,nequ                                                     ! @15/1/86
        do kk=1,nequ                                                  ! @15/1/86
          acov(k,kk)=acov(k,kk)/rn                                    ! @15/1/86
        end do                                                        ! @15/1/86
        do ks=1,nset                                                  ! @15/1/86
          bcov(k,ks)=bcov(k,ks)/rn                                    ! @15/1/86
        end do                                                        ! @15/1/86
      end do                                                          ! @15/1/86
      do k=1,nequ                                    ! Generate A and B &4/12/90
        logs=.false.
        if(k.ne.1.and.nref.ne.0) then                                 ! &4/12/90
          icon=k/2                         ! Test for split constituent &4/12/90
          if(index(icon).gt.0) logs=.true.
        endif                                                         ! @4/12/90
        do kk=1,nequ                                                  ! &4/12/90
          if(.not.logs) a(k,kk)=acov(k,kk)
        end do                                                        ! &4/12/90
        do ks=1,nset                                                  ! &4/12/90
          if(.not.logs) b(k,ks)=bcov(k,ks)
        end do                                                        ! &4/12/90
      end do                                                          ! &4/12/90
      write(nopout,5) ndat
      write(nlog,5) ndat                                              ! @4/12/90
    5 format(/' No. of data points = ',i5/)                  ! &10/1/86 &4/12/90
      if(nref.ne.0.) then                                             ! &4/12/90
        do k=1,ncon                                ! Split constituents &4/12/90
          if(index(k).ne.0) then                                      ! &4/12/90
            indexk=index(k)
            r=hrel(k)/hrel(indexk)
            delta=grel(k)-grel(indexk)
            ics=2*k
            iss=ics+1
            icr=2*indexk
            isr=icr+1
            rc=r*cos(delta)
            rs=r*sin(delta)
            a(ics,ics)=1.
            a(ics,icr)=-rc
            a(ics,isr)=rs
            a(iss,iss)=1.
            a(iss,icr)=-rs
            a(iss,isr)=-rc
          endif                                                       ! @4/12/90
        end do                                                        ! @4/12/90
      endif                                                           ! &4/12/90
C...............................................................................
C     Solve A*C=B :
C...............................................................................
C     NAG section :                                                     &4/12/90
C...............................................................................
C     ifail=0                                                     ! NAG &13/1/86
C     call f04aee(a,maxequ,b,maxequ,nequ,nset,c,maxequ,wkspce,aa, ! NAG &13/1/86
C    $            maxequ,bb,maxequ,ifail)                         ! NAG &13/1/86
C     if(ifail.eq.0) go to 48                                     ! NAG &13/1/86
C...............................................................................
C     MATTID section :                                                  &4/12/90
C...............................................................................
      call aconv(a,amat,maxequ,maxequ,nequ,nequ,nequ,nequ)     ! MATTID @12/5/89
      call aconv(b,bmat,maxequ,maxset,nequ,nset,nequ,nset)     ! MATTID @12/5/89
C     Note that this matrix inversion routine (solvm) puts     ! MATTID &13/1/86
C     result in BMAT                                           ! MATTID &15/5/89
      call solvmx(amat,bmat,nequ,nset,sing,det,dtnorm)         ! MATTID &15/5/89
      call mxcopy(cmat,bmat,nequ,nset)                         ! MATTID &13/1/86
      call aconv(cmat,c,nequ,nset,maxequ,maxset,nequ,nset)     ! MATTID @12/5/89
      if(.not.sing) go to 48                                   ! MATTID &13/1/86
C...............................................................................
      write(nopout,49)                         ! Matrix solution failed &4/12/90
   49 format(' Matrix solution failed')                               ! &10/1/86
      stop
C...............................................................................
   48 do i=1,nset                                                     ! &4/12/90
        rmean(i)=c(1,i)
        do k=1,ncon                                                   ! &4/12/90
          ic=2*k
          is=ic+1
          h(k,i)=sqrt(c(ic,i)**2+c(is,i)**2)
          if(c(is,i).ne.0..or.c(ic,i).ne.0.) then                     ! @19/5/89
            arg=atan2(c(is,i),c(ic,i))/f2
            if(arg.lt.0.) arg=arg+360.
          else                                                        ! @19/5/89
            arg=0.                                                    ! @19/5/89
          endif                                                       ! @19/5/89
          g(k,i)=arg
        end do                                                        ! @4/12/90
      end do                                                          ! @4/12/90
      write(fmt,54) nset                                            ! @14/9/2001
   54 format('(1x,a80/1x,a80/1x,a80///91x,',i5,'a14)')              ! @14/9/2001
      write(nhgo1,fmt) head1,head2,head3,(ahg,i=1,nset)             ! &14/9/2001
C        ..... Output results (1)                            ! &10/1/86 &4/12/90
C                                                               ..... &14/9/2001
      write(fmt,55) nset                                            ! @14/9/2001
   55 format('(/''  MEAN   '',',i5,'(f6.3,8x))')                    ! @14/9/2001
      write(nhgo1,fmt)(rmean(i),i=1,nset)                           ! &14/9/2001
      do k=1,ncon                                                     ! &4/12/90
        astar=' '
        if(index(k).gt.0) astar='*'
        len=lenchar(acon(k))                    ! Convert to upper case @4/12/90
        call changecase(acon(k),acon(k),len,1)                        ! @4/12/90
        write(fmt,57) nset                                          ! @14/9/2001
   57   format('(1x,a6,1x,a1,',i5,'(1x,f5.3,2x,f5.1,1x))')          ! @14/9/2001
        write(nhgo1,fmt) acon(k),astar,(h(k,i),g(k,i),i=1,nset)       ! &4/12/90
C                                                               ..... &14/9/2001
      end do                                                          ! @4/12/90
      do i=1,nset                                  ! Output results (2) &4/12/90
        write(nhgo2,58) head1,head2,head3,ncon
   58   format(1x,a80/1x,a80/1x,a80/1x,i5)                 ! &10/1/86 &14/9/2001
        write(nhgo2,59)(acon(k),k=1,ncon)
   59   format((1x,10(a6)))                       ! &10/1/86 &4/12/90 &14/9/2001
        write(nhgo2,90) rmean(i)
   90   format(1x,f10.3)                                   ! &10/1/86 &14/9/2001
        write(nhgo2,60)(h(k,i),k=1,ncon)
   60   format((1x,10f10.3))                               ! &10/1/86 &14/9/2001
        write(nhgo2,61)(g(k,i),k=1,ncon)
   61   format((1x,10f10.1))                               ! &10/1/86 &14/9/2001
      end do                                                          ! @4/12/90
      do k=1,nset                     ! Compute residual sum of squares &4/12/90
        resk=var(k)
        do i=1,nequ                                                   ! &4/12/90
          resk=resk-2.*c(i,k)*bcov(i,k)
          do ii=1,nequ                                                ! &4/12/90
            resk=resk+c(i,k)*c(ii,k)*acov(i,ii)                       ! &4/12/90
          end do                                                      ! @4/12/90
        end do                                                        ! @4/12/90
        res(k)=resk                                                   ! &4/12/90
      end do                                                          ! @4/12/90
      do k=1,nset                                  ! Compute statistics &4/12/90
        var(k)=var(k)-bcov(1,k)**2                                    ! &15/1/86
        svar(k)=sqrt(var(k))
        sres(k)=sqrt(abs(res(k)))
        r1(k)=res(k)/var(k)
        r2(k)=sres(k)/svar(k)                                         ! &4/12/90
      end do                                                          ! @4/12/90
      write(fmt,69) nset                                            ! @14/9/2001
   69 format('(//'' STATISTICS''//'' DATA SET NO :'',19x,
     $      ',i5,'(i3,9x)/)')                                       ! @14/9/2001
      write(nopout,fmt)(k,k=1,nset)                ! &4/12/90 &1/5/98 &14/9/2001
      write(nlog,fmt)(k,k=1,nset)                          ! @4/12/90 &14/9/2001
      write(fmt,70) nset                                            ! @14/9/2001
   70 format('(/'' VARIANCE ORIGINAL (VO)   '',',i5,'f12.3)')       ! @14/9/2001
      write(nopout,fmt)(var(k),k=1,nset)
      write(nlog,fmt)(var(k),k=1,nset)                     ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(var(k),k=1,nset)                    ! @13/1/86 &14/9/2001
      write(fmt,71) nset                                            ! @14/9/2001
   71 format('('' VARIANCE RESIDUAL (VR)   '',',i5,'f12.5)')        ! @14/9/2001
      write(nopout,fmt)(res(k),k=1,nset)                            ! &14/9/2001
      write(nlog,fmt)(res(k),k=1,nset)                     ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(res(k),k=1,nset)                    ! @13/1/86 &14/9/2001
      write(fmt,72) nset                                            ! @14/9/2001
   72 format('('' RATIO VR/VO'',14x,',i5,'f12.4)')                  ! @14/9/2001
      write(nopout,fmt)(r1(k),k=1,nset)                             ! &14/9/2001
      write(nlog,fmt)(r1(k),k=1,nset)                      ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(r1(k),k=1,nset)                     ! @13/1/86 &14/9/2001
      write(fmt,73) nset                                            ! @14/9/2001
   73 format('('' SQRT(VO)'',17x,',i5,'f12.3)')                     ! @14/9/2001
      write(nopout,fmt)(svar(k),k=1,nset)                           ! &14/9/2001
      write(nlog,fmt)(svar(k),k=1,nset)                    ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(svar(k),k=1,nset)                   ! @13/1/86 &14/9/2001
      write(fmt,74) nset                                            ! @14/9/2001
   74 format('('' SQRT(VR)'',17x,',i5,'f12.5)')                     ! @14/9/2001
      write(nopout,fmt)(sres(k),k=1,nset)                           ! &14/9/2001
      write(nlog,fmt)(sres(k),k=1,nset)                    ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(sres(k),k=1,nset)                   ! @13/1/86 &14/9/2001
      write(fmt,75) nset                                            ! @14/9/2001
   75 format('('' RATIO SQRT(VR)/SQRT(VO)  '',',i5,'f12.4)')        ! @14/9/2001
      write(nopout,fmt)(r2(k),k=1,nset)                             ! &14/9/2001
      write(nlog,fmt)(r2(k),k=1,nset)                      ! @4/12/90 &14/9/2001
      write(nhgo1,fmt)(r2(k),k=1,nset)                     ! @13/1/86 &14/9/2001
      write(nopout,91)
      write(nlog,91)                                                  ! @4/12/90
      write(nhgo1,91)                                                 ! @13/1/86
   91 format(//)
      if(lres) then                                 ! Compute residuals &4/12/90
        rewind ndatin
        rewind ntemp
        read(ndatin,*)                                                 ! &4/7/94
        write(nrsout,1) nset,iye_start                                 ! &1/5/98
        do i=1,ndat                                                   ! &4/12/90
          read(ndatin,*) t,(dat(k),k=1,nset)                           ! &4/7/94
          read(ntemp) (x(k),k=1,nequ)
          do k=1,nset                                                 ! &4/12/90
            do kk=1,nequ                                              ! &4/12/90
              dat(k)=dat(k)-c(kk,k)*x(kk)
            end do                                                    ! @4/12/90
          end do                                                      ! @4/12/90
          write(fmt,2) nset                                         ! @14/9/2001
    2     format('(1x,f10.0,',i5,'f10.3)')                          ! @14/9/2001
          write(nrsout,fmt) t,(dat(k),k=1,nset)                     ! &14/9/2001
        end do                                                        ! @4/12/90
      endif                                                           ! @4/12/90
      stop
      end
C
      subroutine ldast(acon,iye,ncon,nastin,nopout,nlog,sigma,
     $                 j1,j2,v1,v2,v0,year,iint)                      ! &4/12/90
C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Rationalise character variables       JRH:10:01:1986      *
C *                  Data files converted to FORTRAN type                      *
C *                  with extra space at start of each                         *
C *                  record                                JRH:13:01:1986      *
C *                  Minor mod.                            JRH:14:01:1986      *
C *                  Transfer to Sun plus associated mods. JRH:04:12:1990      *
C *                  Mods. for 114 constituents and                            *
C *                  1970-2030                             JRH:04:12:1990      *
C *                  Minor mods. to formats                JRH:14:09:2001      *
C *                                                                            *
C *   SOURCE      :  TDA1154.FVS                                               *
C *   ROUTINE NAME:  LDAST                                                     *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Loads astronomical data from disc (converts to radians    *
C *                  and radians/sec).                                         *
C *                                                                            *
C ******************************************************************************
      integer pcon,ptry                                               ! @4/12/90
      parameter(pcon=114,  ! >/= Max. no. of constituents used &13/1/86 &4/12/90
     $          ptry=3)      ! No. of tries at different const. names # @4/12/90
      dimension sigma(pcon),j1(pcon),j2(pcon),v1(pcon)                ! &10/1/86
      dimension v2(pcon),v0(pcon)
      real j1,j2
      logical first
      logical ltry(ptry)                                              ! @4/12/90
      character*6 acon(pcon)                                 ! @10/1/86 &4/12/90
      character*6 atry(ptry,2)                                        ! @4/12/90
      character*10 filnam                                             ! @4/12/90
      data first/.true./
      data ntry/ptry/                                                 ! @4/12/90
      data atry(1,1),atry(1,2)/'sig1  ','sigma1'/                     ! @4/12/90
      data atry(2,1),atry(2,2)/'the1  ','theta1'/                     ! @4/12/90
      data atry(3,1),atry(3,2)/'lam2  ','lamda2'/                     ! @4/12/90
      save first                                                      ! @4/12/90
      if(first) then                                                  ! &4/12/90
        f1=3.14159265*2./(3600.*360.)
        f2=3.14159265/180.
      endif                                                           ! @4/12/90
      iint=0                                                          ! &4/12/90
      year=31536000.
      if((iye/4)*4.eq.iye) year=year+86400.
      nyear=iye-1969
      write(6,9)                                                      ! @15/1/86
    9 format(/)                                                       ! @15/1/86
      do i=1,ncon                                                     ! &4/12/90
        do k=1,ntry                                                   ! @4/12/90
          ltry(k)=.true.                                              ! @4/12/90
        end do                                                        ! @4/12/90
    1   len=lenchar(acon(i))                                          ! @4/12/90
        filnam(1:len)=acon(i)                                         ! @4/12/90
        filnam(len+1:len+4)='.con'                                    ! @4/12/90
        write(nopout,8) filnam(1:len+4)                      ! &15/1/86 &4/12/90
        write(nlog,8) filnam(1:len+4)                                 ! @4/12/90
    8   format(' Opening file ',a10)                         ! &15/1/86 &4/12/90
        open(unit=nastin,file=filnam(1:len+4),
     $       status='old',iostat=ios)                        ! &10/1/86 &4/12/90
        if(ios.gt.0) then                          ! Error opening file @4/12/90
C     Try other possible names :                                        @4/12/90
          do k=1,ntry                                                 ! @4/12/90
            if(acon(i).eq.atry(k,1).and.ltry(k)) then                 ! @4/12/90
              acon(i)=atry(k,2)                                       ! @4/12/90
              ltry(k)=.false.                                         ! @4/12/90
              go to 1                                                 ! @4/12/90
            endif                                                     ! @4/12/90
            if(acon(i).eq.atry(k,2).and.ltry(k)) then                 ! @4/12/90
              acon(i)=atry(k,1)                                       ! @4/12/90
              ltry(k)=.false.                                         ! @4/12/90
              go to 1                                                 ! @4/12/90
            endif                                                     ! @4/12/90
          end do                                                      ! @4/12/90
          write(nopout,4)                                             ! @4/12/90
          write(nlog,4)                                               ! @4/12/90
    4     format(/' Could not find constituent ',                     ! @4/12/90
     $            ' ..... program terminated'/)                       ! @4/12/90
          stop                                                        ! @4/12/90
        endif                                                         ! @4/12/90
        read(nastin,2) dum
    2   format(1x,f14.10)                                  ! &13/1/86 &14/9/2001
        if(first) sigma(i)=dum*f1                                     ! &4/12/90
    3   do ii=1,62                                                    ! &4/12/90
          read(nastin,5) dj,dv,dvpv
    5     format(3x,f6.3,f8.2,f6.1)                                   ! &13/1/86
          if(ii.eq.nyear) then                                        ! &4/12/90
            j1(i)=dj
            v1(i)=dv*f2
            v0(i)=(dvpv-dv)*f2
          else                                                        ! @4/12/90
            if(ii.eq.nyear+1) then                                    ! &4/12/90
              j2(i)=dj
              v2(i)=dv*f2
            endif                                                     ! @4/12/90
          endif                                                       ! @4/12/90
        end do                                                        ! @4/12/90
        close(unit=nastin)                                            ! &4/12/90
      end do                                                          ! @4/12/90
      first=.false.
      return
      end
C
      subroutine aconv(ain,aout,min,nin,mout,nout,mtrans,ntrans)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:12:05:1989      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1989      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  ACONV                                                     *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Converts array ain(min,nin) to aout(mout,nout) between    *
C *                  the limits (mtrans,ntrans)                                *
C *                                                                            *
C ******************************************************************************
      dimension ain(min,nin),aout(mout,nout)
      if(mtrans.gt.min.or.mtrans.gt.mout.or.
     $   ntrans.gt.nin.or.ntrans.gt.nout) stop                 ! Error condition
      do i=1,mout
        do j=1,nout
          aout(i,j)=0.                                   ! Initially clear array
        end do
      end do
      do i=1,mtrans
        do j=1,ntrans
          aout(i,j)=ain(i,j)
        end do
      end do
      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 dskopen(nin,nout,nlog,ndev,nbt,anot,stat)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:12:01:1987      *
C *                                                                            *
C *   REVISION    :  Inclusion of binary file option       JRH:21:01:1987      *
C *                  Minor mod.                            JRH:24:01:1987      *
C *                  Mod. to allow direct access for                           *
C *                  binary file                           JRH:18:02:1987      *
C *                  Mod. to include INQUIRE               JRH:31:03:1987      *
C *                  Mod. so that files are always SHARED  JRH:09:11:1988      *
C *                  Removal of "shared" option in "open"  JRH:06:07:1989      *
C *                  "recordsize" changed to "recl"        JRH:12:07:1989      *
C *                  Mods. for change from words to bytes                      *
C *                  in "open"                             JRH:12:07:1989      *
C *                  Mod. to INQUIRE option, since UNIX                        *
C *                  direct-access binary files do not store                   *
C *                  record length. It is now assumed that,                    *
C *                  when this option is used, the file is                     *
C *                  DIRECT ACCESS BINARY and that first 4                     *
C *                  bytes hold (no. bytes)/4 (ie. no. of                      *
C *                  VMS words).                           JRH:19:07:1989      *
C *                  Mod. (for SUN) so that record length                      *
C *                  omitted for sequential direct-access                      *
C *                  file                                  JRH:22:08:1989      *
C *                  Expansion of FILNAM to 30 characters  JRH:22:09:1989      *
C *                  Expansion of FILNAM to 40 characters  JRH:02:05:1990      *
C *                  Minor modification to formats         JRH:15:05:2000      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  DSKOPEN                                                   *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Interactive routine for opening disc file.                *
C *                  Modified version of DOPEN with O/P to log file.           *
C *                                                                            *
C *                  NOTE : Log file must be first file to be openned.         *
C *                                                                            *
C *                    NIN ..... Operator input LU (5)                         *
C *                    NOUT .... Operator output LU (6)                        *
C *                    NLOG .... Log file LU                                   *
C *                    NDEV .... Device LU to be openned                       *
C *                    NBT  .... For ASCII file, 0                             *
C *                              For binary sequential access file, no. of     *
C *                              bytes per record                              *
C *                              For binary direct access file, -(no. of       *
C *                              bytes per record)                             *
C *                    ANOT .... Annotation to appear on operator's terminal   *
C *                              (40 chars.)                                   *
C *                    STAT .... (a) 'INQUIRE' .... File is understood to      *
C *                                  exist and to be binary, and no. of bytes  *
C *                                  per record is returned in NBT             *
C *                              (b) otherwise, STAT is status of file         *
C *                                                                            *
C ******************************************************************************
      character*40 anot
      character*10 stat,status,acc                                    ! &31/3/87
      character*40 filnam                                     ! &22/9/89 &2/5/90
      write(nout,1) anot,' '                                          ! &24/1/87
    1 format(/1x,a40,a1)                                              ! &24/1/87
      read(nin,2) filnam
    2 format(a40)                                             ! &22/9/89 &2/5/90
      if(nbt.eq.0) then                                ! ASCII @21/1/87 &12/7/89
        open(unit=ndev,file=filnam,status=stat)               ! &9/11/88 &6/7/89
      else                                                     ! Binary @21/1/87
        if(nbt.gt.0) then                  ! Sequential access @31/3/87 &12/7/89
          acc='SEQUENTIAL'                                            ! @31/3/87
        else                                            ! Direct access @31/3/87
          acc='DIRECT'                                                ! @31/3/87
        endif                                                         ! @31/3/87
        if(stat.eq.'INQUIRE') then                                    ! @31/3/87
          status='OLD'                                                ! @31/3/87
C         inquire(file=filnam,recl=nbyte) ! Find record length @31/3/87 &19/7/89
C     (Vax Note : record length is apparently given in bytes if file is not 
C     openned)                                                        ! &12/7/89
          open(unit=ndev,file=filnam,status=status,                   ! @19/7/89
     $         form='UNFORMATTED',recl=4,                             ! @19/7/89
     $         access='DIRECT')                                       ! @19/7/89
          read(ndev,rec=1) imax                                       ! @19/7/89
          nbyte=imax*4                                                ! @19/7/89
          close(unit=ndev)                                            ! @19/7/89
          nbt=nbyte                                          ! @31/3/87 &12/7/89
        else                                                          ! @31/3/87
          status=stat                                                 ! @31/3/87
          nbyte=iabs(nbt)                ! Record length given @31/3/87 &12/7/89
        endif                                                         ! @31/3/87
        if(acc.eq.'SEQUENTIAL') then                                  ! @22/8/89
          open(unit=ndev,file=filnam,status=status,                   ! @22/8/89
     $         form='UNFORMATTED',                                    ! @22/8/89
     $         access=acc)                                            ! @22/8/89
        else                                                          ! @22/8/89
          open(unit=ndev,file=filnam,status=status,                   ! @31/3/87
     $         form='UNFORMATTED',recl=nbyte,                ! @31/3/87 &12/7/89
     $         access=acc)                           ! &31/3/87 &9/11/88 &6/7/89
        endif                                                         ! @22/8/89
      endif                                                           ! @21/1/87
      write(nlog,3) anot,filnam                                       ! &19/2/87
    3 format(/1x,a40,1x,a40)                         ! &19/2/87 &22/9/89 &2/5/90
      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
      subroutine mxcopy(a,b,m,n)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MATHS                             REF:JRH:14:01:1986      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1986      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  MXCOPY                                                    *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Forms matrix equality A=B.                                *
C *                  A and B both have M rows and N cols. (M*N).               *
C *                                                                            *
C ******************************************************************************
      dimension a(m,n),b(m,n)
      do 1 i=1,m
      do 1 j=1,n
      a(i,j)=b(i,j)
  1   continue
      return
      end
C
      subroutine mxgaus(a,ind,n,sing,det)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MATHS                             REF:JRH:14:01:1986      *
C *                                                                            *
C *   REVISION    :  Increase of array dimensions          JRH:14:09:2001      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  MXGAUS                                                    *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Performs an in-situ reduction of an N*N matrix A into     *
C *                  quasi-triangular matrices L (a strictly lower triangle)   *
C *                  and U (an upper triangle) by the Gauss-Doolittle method.  *
C *                  The pivots are along the diagonal of U such that :        *
C *                    A=(L+I)*U                                               *
C *                  No permutations of rows or columns of A are performed and *
C *                  the pivotal row subscripts are recorded in the non-local  *
C *                  vector IND(N). Because of this, this procedure is not     *
C *                  recommended for use separately, but is called by SOLVMX.  *
C *                  SING is assigned .TRUE. if A is ill-conditioned,          *
C *                  otherwise .FALSE..                                        *
C *                  DET is determinant of A.                                  *
C *                                                                            *
C *                  Ref.: First Course in Numerical Analysis, Ralston,        *
C *                        Ch. 9, P. 394.                                      *
C *                                                                            *
C ******************************************************************************
      dimension a(n,n),ind(n)
C     NOTE :
C     ******
C     ARRAY D SHOULD HAVE DIMENSION AT LEAST N (OR "M" IN SOLVM) :
C
      double precision d(229),dspace,x
      integer r,rr
      logical sing
      det=0.
      nrc=2
      do 99 r=1,n
      do 1 k=1,n
      d(k)=dble(a(k,r))
  1   continue
      if (r.eq.1) goto 7
      rr=r-1
      do 2 j=1,rr
      jj=ind(j)
      dspace=d(jj)
      a(j,r)=sngl(dspace)
      d(jj)=d(j)
      jj=j+1
      do 3 i=jj,n
      d(i)=d(i)-dble(a(i,j))*dspace
  3   continue
  2   continue
  7   dspace=d(r)
      if (n.ne.r) goto 10
      if (dspace.eq.0.d0) goto 81
      ind(n)=n
      a(n,n)=sngl(dspace)
      goto 70
 10   ii=r
      rr=r+1
      do 4 i=rr,n
      if (dabs(dspace)-dabs(d(i))) 5,4,4
  5   dspace=d(i)
      ii=i
  4   continue
      if (dspace.eq.0.d0) goto 81
      ind(r)=ii
      if (r.ne.ii) nrc=nrc+1
      a(r,r)=sngl(dspace)
      jj=ind(r)
      d(jj)=d(r)
      do 6 i=rr,n
      x=d(i)/dspace
      a(i,r)=sngl(x)
  6   continue
 99   continue
C     THIS METHOD CAUSES THE LOWER TRIANGLE TO BE PERMUTATED WITH
C     RESPECT TO THE PIVOTAL ROW SUBCRIPTS.
C     REARRANGING LOWER TRIANGLE...
 70   if (n.eq.2) goto 80
      nn=n-1
      do 8 i=2,nn
      ii=i-1
      iii=ind(i)
      do 9 j=1,ii
      y=a(i,j)
      a(i,j)=a(iii,j)
      a(iii,j)=y
  9   continue
  8   continue
 80   sing=.false.
C     CALCULATING DETERMINANT OF A
C     DET A=PRODUCT OF DIAGONAL TERMS OF UPPER TRI-ANGLE
      dspace=dble(a(1,1))
      do 16 i=2,n
      dspace=dspace*dble(a(i,i))
 16   continue
      det=sngl(dspace)
C     CALCULATING SIGN OF DET A
      i=nrc-2*(nrc/2)
      if (i.ne.0) det=-det
      return
81    sing=.true.
      return
      end
C
      subroutine solvmx(a,b,m,n,sing,det,dtnorm)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MATHS                             REF:JRH:14:01:1986      *
C *                                                                            *
C *   REVISION    :  Removal of spurious characters in                         *
C *                  header                                JRH:18:12:1989      *
C *                  Increase of array dimensions          JRH:14:09:2001      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  SOLVMX                                                    *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Solves the matrix equation AX=B writing the result in B.  *
C *                  A is destroyed, being left in the Gauss-Doolittle form.   *
C *                  A has M rows and M cols. (M*M).                           *
C *                  B has M rows and N cols. (M*N).                           *
C *                  A is M*M, B is M*N.                                       *
C *                  SING is assigned .TRUE. if A is ill-conditioned,          *
C *                  otherwise .FALSE..                                        *
C *                  DET is determinant of A.                                  *
C *                  DTNORM is the normalised determinant of A.                *
C *                  Calls MXGAUS.                                             *
C *                                                                            *
C ******************************************************************************
C     NOTE :
C     ******
C     ARRAYS IND AND D SHOULD HAVE DIMENSION AT LEAST M :
C
      dimension a(m,m),b(m,n),ind(229)
      double precision d(229),dspace
      logical sing
C     CALCULATING FACTOR FOR NORMALISING DETERMINANT OF A
      do 13 i=1,m
      dspace=0.d0
      do 14 j=1,m
      dspace=dspace+dble(a(i,j))**2
 14   continue
      d(i)=dspace
 13   continue
      dspace=d(1)
      do 15 i=2,m
      dspace=dspace*d(i)
 15   continue
      d(1)=dsqrt(dspace)
C     L,U TRANSFORMATION OF A
      call mxgaus(a,ind,m,sing,det)
      if (sing) goto 80
C     NORMALISED DETERMINANT
      dspace=dble(det)
      dtnorm=sngl(dspace/d(1))
C     PERMUTATING THE ROWS OF B
      mm=m-1
      do 1 i=1,mm
      ii=ind(i)
      do 9 j=1,n
      x=b(i,j)
      b(i,j)=b(ii,j)
      b(ii,j)=x
  9   continue
  1   continue
C     CALCULATING Y AND INSERTING IT IN B. (I.E.LY=B)
      do 10 k=1,n
      do 4 i=1,m
      dspace=0.d0
      if (i.eq.1) goto 3
      ii=i-1
      do 2 j=1,ii
      dspace=dspace+dble(a(i,j))*d(j)
  2   continue
  3   d(i)=dble(b(i,k))-dspace
  4   continue
      do 11 i=1,m
      b(i,k)=sngl(d(i))
 11   continue
 10   continue
C     CALCULATING X AND INSERTING IT IN B. (I.E. UX=Y)
      do 12 k=1,n
      do 7 i=1,m
      ii=m+1-i
      dspace=0.d0
      if (ii.eq.m) goto 6
      jj=ii+1
      do 5 j=jj,m
      dspace=dspace+dble(a(ii,j))*d(j)
  5   continue
  6   d(ii)=(dble(b(ii,k))-dspace)/dble(a(ii,ii))
  7   continue
      do 8 i=1,m
      b(i,k)=sngl(d(i))
  8   continue
 12   continue
 80   return
      end
C
      logical function yesno(nin,nout,nlog,anot)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:27:08:1986      *
C *                                                                            *
C *   REVISION    :  Minor mod. to output                  JRH:19:02:1987      *
C *                  Minor mod. to formats                 JRH:15:05:2000      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  YESNO                                                     *
C *   TYPE        :  LOGICAL FUNCTION                                          *
C *                                                                            *
C *   FUNCTION    :  Returns .TRUE. or .FALSE. depending on operator reply     *
C *                  (Y or N, respectively).                                   *
C *                                                                            *
C *                  NIN ..... Operator input device                           *
C *                  NOUT .... Operator output device                          *
C *                  NLOG .... Log device                                      *
C *                  ANOT .... Anotation (40 chars.)                           *
C *                                                                            *
C ******************************************************************************
      character*40 anot
      character*1 adum
    4 write(nout,1) anot,' '                                          ! &19/2/87
      write(nlog,1) anot,' '                                          ! &19/2/87
    1 format(/1x,a40,a1)                                              ! &19/2/87
      read(nin,2) adum
    2 format(a1)
      write(nlog,3) adum
    3 format(1x,a1)
      if(adum.ne.'Y'.and.adum.ne.'N'.and.
     $   adum.ne.'y'.and.adum.ne.'n') go to 4
      if(adum.eq.'Y'.or.adum.eq.'y') then
        yesno=.true.
      else
        yesno=.false.
      endif
      return
      end
C
