C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Rationalise CHARACTER variables       JRH:09:01:1986      *
C *                  Inclusion of diagnostic option        JRH:09:01:1986      *
C *                  Upper to lower case in operator O/P   JRH:09:01:1986      *
C *                  Disc channel 6 to 8                   JRH:09:01:1986      *
C *                  PRTNEW.DAT renamed PORT.DAT           JRH:10:01:1986      *
C *                  Inclusion of DOPEN                    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. to diagnostic mode         JRH:13:01:1986      *
C *                  Minor mod.                            JRH:14:01:1986      *
C *                  "Hidden" facility to convert feet to                      *
C *                  metres                                JRH:14:01:1986      *
C *                  Minor mod.                            JRH:15:01:1986      *
C *                  Mod. for HW and LW with other than                        *
C *                  semi-diurnal tides                    JRH:18:05:1989      *
C *                  Transfer to Sun plus associated mods. JRH:29:11: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:29:11:1990      *
C *                  Continuation of above mods.           JRH:30:11:1990      *
C *                  Inclusion of log file etc.            JRH:30:11:1990      *
C *                  Removal of subroutines already in                         *
C *                  FORLIB                                JRH:30:11:1990      *
C *                  Continuation of above mods.           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 *                  Move of label 49 to correct position  JRH:08:04:1995      *   
C *                  Initialisation of itest and n removed                     *
C *                  from data statements                  JRH:25:02:1998      *
C *                  Minor mods. to formats                JRH:17:05:2000      *
C *                  Renamed tpd1150.f                     JRH:14:09:2001      *
C *                  All necessary subroutines added to                        *
C *                  source code                           JRH:14:09:2001      *
C *                  Mod. to change diagnostic output from                     *
C *                  one month to one year                 JRH:14:09:2001      *
C *                  Addition of header re. software                           *
C *                  licensing                             JRH:23:12:2002      *
C *                                                                            *
C *   SOURCE      :  tpd1150.f                                                 *
C *   ROUTINE NAME:  tpd1150                                                   *
C *   TYPE        :  MAIN                                                      *
C *                                                                            *
C *   FUNCTION    :  Tide prediction 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,ptry                                               ! @4/12/90
      parameter(pcon=114,        ! >/= Maximum no. of constituents used @4/12/90
     $          ptry=3)      ! No. of tries at different const. names # @4/12/90
      dimension h(pcon),g(pcon),j1(pcon),vpv(pcon),v1(pcon)           ! &4/12/90
      dimension j2(pcon),v2(pcon),v0(pcon),j(pcon)                    ! &4/12/90
      dimension pred1(0:24)                                           ! &18/5/89
      dimension ihrt(2,2,2),mint(2,2,2),pred(2,2),tlst(2)
      dimension ihrhat(2),minhat(2),ihrlat(2),minlat(2)
      dimension sigma(pcon)                                           ! &4/12/90
      dimension sn(2001)                              ! Dimensions N+1 &29/11/90
      dimension iarray(5)                                            ! @30/11/90
      character*80 head1,head2                                         ! @9/1/86
      character*80 port,port1                                          ! @9/1/86
      character*40 aarray(1)                                         ! @30/11/90
      character*10 filnam                                             ! @3/12/90
      character*6 am2,acon(pcon)                              ! @9/1/86 &4/12/90
      character*6 atry(ptry,2)                                        ! @4/12/90
      character*5 amnhat,amnlat,amon(12)                               ! @9/1/86
      character*1 atype,alen,units,awtype                             ! @18/5/89
      logical ltry(ptry)                                              ! @4/12/90
      common sigma,sn,itest,rn,n2,n4,rn2,rn4,f2,f4,ncon,n
      real j1,j2,j
      data (amon(i),i=1,12)/'  JAN','  FEB','MARCH','APRIL','  MAY',   ! &9/1/86
     1' JUNE',' JULY','  AUG',' SEPT','  OCT','  NOV','  DEC'/         ! &9/1/86
      data am2/'m2    '/                                             ! &29/11/90
      data iint/0/
C     data itest/0/                                                   ! &25/2/98
C     data n/2000/                                                    ! &25/2/98
      data rpt/600000./     ! Set interval for interpolation of J, VPV &29/11/90
      data nopin,nptin,nastin,nopout,nprout,nlog
     $    /5    ,7    ,9     ,6     ,8     ,10  / !Set devices &9/1/86 &30/11/90
      data limit,eps/10,4.e-7/        ! Set values for Newton's method &29/11/90
      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
      itest=0                                                         ! @25/2/98
      n=2000                                                          ! @25/2/98
      f2m=12.*2.54/100.                                ! Feet to metres &13/1/86
      call dskopen(nopin,nopout,nlog,nlog,0,                         ! @30/11/90
     $             'Filename for log file ?                 ',       ! @30/11/90
     $             'unknown   ')                                     ! @30/11/90
      call dskopen(nopin,nopout,nlog,nprout,0,                       ! @30/11/90
     $             'Filename for output   ?                 ',       ! @30/11/90
     $             'unknown   ')                                     ! @30/11/90
      open(unit=nptin,file='port.dat',status='old')         ! &10/1/86 &30/11/90
      f1=3.14159265*2./(3600.*360.)
      rn=n                                                    ! Set SN &29/11/90
      f4=3.14159265/(2*rn)
      np1=n+1
      do i=1,np1                                                     ! &29/11/90
        ri=i-1
        sn(i)=sin(ri*f4)                                             ! &29/11/90
      end do                                                         ! @29/11/90
   36 write(nopout,30)
      write(nlog,30)                                                 ! @30/11/90
   30 format(' Input port : ')                                         ! &9/1/86
      read(nopin,7) port
    7 format(a80)                                                     ! &13/1/86
      write(nlog,6) port                                             ! @30/11/90
    6 format(1x,a80)                                                 ! @30/11/90
C     Scan port for "!" (symbol for conv. ft. to m.)                  ! @14/1/86
      units='M'                                                       ! @14/1/86
      do i=1,80                                                       ! @14/1/86
        if(port(i:i).eq.'!') then                                     ! @14/1/86
          units='F'                                                   ! @14/1/86
          port(i:i)=' '                                               ! @14/1/86
        endif                                                         ! @14/1/86
      end do                                                          ! @14/1/86
      call changecase(port,port,80,0)           ! Convert to lower case @3/12/90
   33 read(nptin,95) port1                                            ! &13/1/86
   95 format(1x,a80)                                                  ! @13/1/86
      call changecase(port1,port1,80,0)         ! Convert to lower case @3/12/90
      if(port1(1:5).eq.' stop') then                ! &9/1/86 &29/11/90 &3/12/90
        write(nopout,35)
        write(nlog,35)                                               ! @30/11/90
   35   format(' Port not listed')                                     ! &9/1/86
        rewind(nptin)
        go to 36
      else if(port.ne.port1) then                            ! &9/1/86 &29/11/90
        read(nptin,95) port1                                ! &13/1/86 &29/11/90
        read(nptin,95) port1                                          ! &13/1/86
        read(nptin,81) ncon
        nrec=1+3*((ncon-1)/10+1)
        do i=1,nrec                                                  ! &29/11/90
          read(nptin,95) port1                              ! &13/1/86 &29/11/90
        end do                                                       ! @29/11/90
        go to 33
      endif                                                          ! @29/11/90
      read(nptin,95) head1               ! Input header lines @10/1/86 &29/11/90
      read(nptin,95) head2                                            ! @10/1/86
      read(nptin,81) ncon  ! Input no. of constituents and const. list &29/11/90
   81 format(1x,i5)                                                   ! &10/1/86
      read (nptin,82)(acon(i),i=1,ncon)
   82 format((1x,10(a6)))                                   ! &10/1/86 &29/11/90
      do i=1,ncon                              ! Convert to lower case @29/11/90
        len=lenchar(acon(i))                                         ! @29/11/90
        call changecase(acon(i),acon(i),len,0)                        ! @3/12/90
      end do                                                         ! @29/11/90
      do i=1,ncon                                ! Find M2 constituent &29/11/90
        if(acon(i).eq.am2) then                                      ! &29/11/90
          im2=i                                                      ! @29/11/90
          go to 89                                                   ! @29/11/90
        endif                                                        ! @29/11/90
      end do                                                         ! @29/11/90
      im2=0
   89 continue
      read(nptin,2) pmean            ! Input mean value of predictions &29/11/90
    2 format(1x,f10.3)                                                ! &10/1/86
      read(nptin,3) (h(i),i=1,ncon)                       ! Input H, G &29/11/90
    3 format((1x,10f10.3))                                            ! &10/1/86
      read(nptin,4) (g(i),i=1,ncon)
    4 format((1x,10f10.1))                                            ! &10/1/86
      if(units.eq.'F') then  ! Convert from ft. to m. if nec. @14/1/86 &29/11/90
        pmean=pmean*f2m                                               ! @14/1/86
        do i=1,ncon                                                   ! @14/1/86
          h(i)=h(i)*f2m                                               ! @14/1/86
        end do                                                        ! @14/1/86
      endif                                                           ! @14/1/86
   41 call iread(nopin,nopout,nlog,                                  ! @30/11/90
     $           'Input range of years (first,last)       ',' ',     ! @30/11/90
     $           iarray,0,0,                                         ! @30/11/90
     $           'x,2i5                                   ')         ! @30/11/90
      call iread(nopin,nopout,nlog,                                  ! @30/11/90
     $           '(1970-2030) :                           ','$',      ! @4/12/90
     $           iarray,1,2,                                         ! @30/11/90
     $           'x,2i5                                   ')         ! @30/11/90
      iyfst=iarray(1)                                                ! @30/11/90
      iylst=iarray(2)                                                ! @30/11/90
      if(iyfst.gt.iylst.or.iyfst.lt.1970.or.iylst.gt.2030) then      ! &4/12 /90
        write(nopout,40)
        write(nlog,40)                                               ! @30/11/90
   40   format(' Input error')                                         ! &9/1/86
        go to 41
      endif                                                          ! @29/11/90
   39 call aread(nopin,nopout,nlog,                                  ! @30/11/90
     $           'Type H for hourly predictions,          ',' ',     ! @30/11/90
     $           aarray,0,0,                                         ! @30/11/90
     $           'x,a40                                   ')         ! @30/11/90
      call aread(nopin,nopout,nlog,                                  ! @30/11/90
     $           '     W for HW and LW,                   ',' ',     ! @30/11/90
     $           aarray,0,0,                                         ! @30/11/90
     $           'x,a40                                   ')         ! @30/11/90
      call aread(nopin,nopout,nlog,                                  ! @30/11/90
     $           '     D for diagnostic output :          ','$',     ! @30/11/90
     $           aarray,1,1,                                         ! @30/11/90
     $           'x,a40                                   ')         ! @30/11/90
      atype=aarray(1)(1:1)                                           ! @30/11/90
      call changecase(atype,atype,1,0)          ! Convert to lower case @3/12/90
      if(atype.ne.'h'.and.atype.ne.'w'.and.atype.ne.'d') then  !&9/1/86 &3/12/90
        write(nopout,40)
        write(nlog,40)                                               ! @30/11/90
        go to 39
      endif                                                          ! @29/11/90
      if(atype.ne.'h'.and.atype.ne.'d') then                  ! &9/1/86 @3/12/90
   46   call aread(nopin,nopout,nlog,                                ! @30/11/90
     $             'Type S for short output,                ',' ',   ! @30/11/90
     $             aarray,0,0,                                       ! @30/11/90
     $             'x,a40                                   ')       ! @30/11/90
        call aread(nopin,nopout,nlog,                                ! @30/11/90
     $             '     L for long output :                ','$',   ! @30/11/90
     $             aarray,1,1,                                       ! @30/11/90
     $             'x,a40                                   ')       ! @30/11/90
        alen=aarray(1)(1:1)                                          ! @30/11/90
        call changecase(alen,alen,1,0)          ! Convert to lower case @3/12/90
        if(alen.ne.'s'.and.alen.ne.'l') then                 ! &18/5/89 &3/12/90
          write(nopout,40)
          write(nlog,40)                                             ! @30/11/90
          go to 46
        endif                                                        ! @29/11/90
   97   call aread(nopin,nopout,nlog,                                ! @30/11/90
     $             'Type S for semi-diurnal tide,           ',' ',   ! @30/11/90
     $             aarray,0,0,                                       ! @30/11/90
     $             'x,a40                                   ')       ! @30/11/90
        call aread(nopin,nopout,nlog,                                ! @30/11/90
     $             '     M for mixed or diurnal tide :      ','$',   ! @30/11/90
     $             aarray,1,1,                                       ! @30/11/90
     $             'x,a40                                   ')       ! @30/11/90
        awtype=aarray(1)(1:1)                                        ! @30/11/90
        call changecase(awtype,awtype,1,0)      ! Convert to lower case @3/12/90
        if(awtype.ne.'s'.and.awtype.ne.'m') then             ! @18/5/89 &3/12/90
          write(nopout,40)                                            ! @18/5/89
          write(nlog,40)                                             ! @30/11/90
          go to 97                                                    ! @18/5/89
        endif                                                         ! @18/5/89
      endif                                                          ! @29/11/90
C     Main loop :
      if(atype.eq.'d') iylst=iyfst       ! Force only one year @13/1/86 &3/12/90
      do iyed=iyfst,iylst                         ! Year loop &13/1/86 &29/11/90
        iye=iyed
        write(nopout,90) iye
        write(nlog,90) iye                                           ! @30/11/90
   90   format(/' Starting year',i5/)                        ! &15/1/86 &4/12/90
        nyear=iye-1969                        ! Input SIGMA, J, V, VPV &29/11/90
        do i=1,ncon                                                  ! &29/11/90
          do k=1,ntry                                                 ! @4/12/90
            ltry(k)=.true.                                            ! @4/12/90
          end do                                                      ! @4/12/90
    1     len=lenchar(acon(i))                                        ! @3/12/90
          filnam(1:len)=acon(i)                                       ! @3/12/90
          filnam(len+1:len+4)='.con'                                  ! @3/12/90
          write(nopout,91) filnam(1:len+4)                   ! &14/1/86 &3/12/90
   91     format(' Opening file ',a10)                       ! &15/1/86 &3/12/90
          write(nlog,9) filnam(1:len+4)                               ! @3/12/90
    9     format(' Opening file ',a10)                                ! &3/12/90
          open(unit=nastin,file=filnam(1:len+4),
     $         status='old',iostat=ios)                       ! &9/1/86 &3/12/90
          if(ios.gt.0) then                       ! Error opening file @29/11/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,5)                                          ! @30/11/90
            write(nlog,5)                                            ! @30/11/90
    5       format(/' Could not find constituent ',                  ! @30/11/90
     $              ' ..... program terminated'/)                    ! @30/11/90
            stop                                                     ! @30/11/90
          endif                                                      ! @29/11/90
          read(nastin,85) dum
   85     format(1x,f14.10)                                           ! &13/1/86
          if(iyed.eq.iyfst) then                                     ! &29/11/90
            sigma(i)=dum*f1                                ! Set SIGMA &29/11/90
          endif                                                      ! @29/11/90
          do ii=1,62                                                  ! &4/12/90
            read(nastin,51) dj,dv,dvpv
   51       format(3x,f6.3,f8.2,f6.1)                                 ! &13/1/86
            if(ii.eq.nyear) then                                     ! &29/11/90
              j1(i)=dj
              v1(i)=dv
              vpv(i)=dvpv
            else if(ii.eq.nyear+1) then                              ! &29/11/90
              j2(i)=dj
              v2(i)=dv
            endif                                                    ! &29/11/90
          end do                                                     ! @29/11/90
          close(unit=nastin)                                         ! &29/11/90
        end do                                                       ! @29/11/90
        do i=1,ncon                                                  ! &29/11/90
          v0(i)=vpv(i)-v1(i)                                         ! &29/11/90
        end do                                                       ! @29/11/90
        year=31536000.
        if((iye/4)*4.eq.iye) year=year+86400.
        if(atype.ne.'w') then   ! Hourly or Diagnostic predict. &9/1/86 &3/12/90
          if(atype.eq.'d') then                               ! @9/1/86 &3/12/90
C     Following removed, 14/9/2001:
C           call iread(nopin,nopout,nlog,                            ! @30/11/90
C    $                 'Input month number :                    ','$',!@30/11/90
C    $                 iarray,1,1,                                   ! @30/11/90
C    $                 'x,2i5                                   ')   ! @30/11/90
C           idiag=iarray(1)                                          ! @30/11/90
C           write(nprout,93) iye,amon(idiag)                           ! @9/1/86
C  93       format('     2',i5,a5)                                    ! @10/1/86
            write(nprout,93) iye                                    ! @14/9/2001
   93       format('     2',i5)                                     ! @14/9/2001
          endif                                                        ! @9/1/86
          t=0.
    8     call etd(iye,mon,idy,ihr,min,sec,t)
          if(idy.eq.1.and.atype.ne.'d') then                  ! &9/1/86 &3/12/90
            write(nprout,10) head1
   10       format(1h1///////////1h ,a80)                             ! &13/1/86
            write(nprout,11) head2
   11       format(/1h ,a80)                                          ! &13/1/86
            write(nprout,12) amon(mon),iye
   12       format(/61x,a5,i5)
            write(nprout,13)
   13       format(/'  DATE',58x,'TIME')                              ! &3/12/90
            write(nprout,14)
   14       format('       0000 0100 0200 0300 0400 0500 0600 0700', ! &30/11/90
     $             ' 0800 0900 1000 1100 1200 1300 1400 1500',       ! &30/11/90
     $             ' 1600 1700 1800 1900 2000 2100 2200 2300'/)      ! &30/11/90
          endif                                                      ! @30/11/90
          do i=1,24                                                  ! &30/11/90
            ri=i
            tim=t+(ri-1.)*3600.
            idum=tim/rpt+1.  ! Mod. to interp. J, VPV on regular basis &30/11/90
            if(iint.ne.idum) then                                    ! &30/11/90
              iint=idum
              ridum=idum
              tdum=(ridum-0.5)*rpt
              do k=1,ncon                                            ! &30/11/90
                j(k)=(j1(k)*(year-tdum)+j2(k)*tdum)/year
                vpv(k)=v0(k)+(v1(k)*(year-tdum)+v2(k)*tdum)/year     ! &30/11/90
              end do                                                 ! @30/11/90
            endif                                                    ! @30/11/90
            pred1(i)=tide(tim,h,g,j,vpv)+pmean                       ! &30/11/90
          end do                                                     ! @30/11/90
          if(atype.eq.'h') then                               ! @9/1/86 &3/12/90
            write(nprout,26) idy,(pred1(i),i=1,24)
   26       format(1x,i4,1h ,24f5.2)                                  ! &10/1/86
          else                                                         ! @9/1/86
C           if(mon.eq.idiag) then                           ! @9/1/86 &14/9/2001
            do i=1,24                                                  ! @9/1/86
              write(nprout,94) t+float(i-1)*3600.,pred1(i),pred1(i)*2. ! @9/1/86
   94         format(1x,f10.0,2f10.3)                                 ! @10/1/86
            end do                                                     ! @9/1/86
C           endif                                                      ! @9/1/86
          endif                                             ! @9/1/86 &14/9/2001
          t=t+86400.
          if(mon.ne.12.or.idy.ne.31) go to 8
        else                                   ! HW and LW predictions @30/11/90
          if(awtype.eq.'s') then           ! Semi-diurnal tide @18/5/89 &3/12/90
            if(im2.eq.0) then             ! Abort if no M2 constituent &30/11/90
              write(nopout,109)                                      ! @30/11/90
              write(nlog,109)                                        ! @30/11/90
  109         format(/' M2 absent from data -',                      ! @30/11/90
     $                ' use hourly predictions instead',             ! @30/11/90
     $                ' ..... program terminated'/)                  ! @30/11/90
              stop                                                   ! @30/11/90
            endif                                                    ! @30/11/90
            iyp=iye           ! Set initial values, using M2 tide only &30/11/90
            monp=0
            idp=1
            dt=44714.
            hthat=pmean
            htlat=pmean
            thetah=g(im2)-vpv(im2)
            if(thetah.lt.0.) thetah=thetah+360.
            thetal=thetah-180.
            if(thetal.lt.0.) thetal=thetal+360.
            t=thetal
            if(thetah.lt.thetal) t=thetah
            tl=t-180.
            t=t*dt/360.
            tl=tl*dt/360.
            tlst(1)=tl
            tlst(2)=tl
   54       do i=1,limit                             ! Newton's method &30/11/90
              idum=t/rpt+1.  ! Mod. to interp. J, VPV on regular basis &30/11/90
              if(iint.ne.idum) then                                  ! &30/11/90
                iint=idum
                ridum=idum
                tdum=(ridum-0.5)*rpt
                do k=1,ncon                                          ! &30/11/90
                  j(k)=(j1(k)*(year-tdum)+j2(k)*tdum)/year
                  vpv(k)=v0(k)+(v1(k)*(year-tdum)+v2(k)*tdum)/year   ! &30/11/90
                end do                                               ! @30/11/90
              endif                                                  ! @30/11/90
              diff=tided(t,h,g,j,vpv,1)                              ! &30/11/90
              if(abs(diff).lt.eps) go to 56
              t=t-diff/tided(t,h,g,j,vpv,2)                          ! &30/11/90
            end do                                                   ! @30/11/90
            write(nopout,80)                     ! Failure to converge &30/11/90
            write(nlog,80)                       ! Failure to converge &30/11/90
   80       format(/' Failure to converge',                          ! &30/11/90
     $              ' ..... program terminated'/)                    ! @30/11/90
            stop
   56       if(t.lt.0.) then                                         ! &30/11/90
              tlst(1)=t
              tlst(2)=t
              t=t+dt/2.
              go to 54
            endif                                                    ! @30/11/90
            call etd(iye,mon,idy,ihr,min,sec,t)                      ! @30/11/90
            if(monp.ne.0) then                                       ! &30/11/90
              if(idy.eq.idp) go to 59
              if(alen.ne.'s') then                 ! Output day &9/1/86 &3/12/90
                nform=1
                do ii=1,2                                            ! &30/11/90
                  do jj=1,2                                          ! &30/11/90
                    if(ihrt(ii,jj,1).eq.-1) nform=2*ii+jj-1          ! &30/11/90
                  end do                                             ! &30/11/90
                end do                                               ! &30/11/90
                if(nform.eq.1)                                       ! &30/11/90
     $            write(nprout,61) idp,((ihrt(ii,jj,1),ihrt(ii,jj,2),! &30/11/90
     $                                   mint(ii,jj,1),mint(ii,jj,2),! &30/11/90
     $                                   pred(ii,jj),jj=1,2),ii=1,2) ! &30/11/90
   61           format(i5,2(2x,4i1,f6.2),1x,2(2x,4i1,f6.2))
                if(nform.eq.2)                                       ! &30/11/90
     $            write(nprout,62) idp,                              ! &30/11/90
     $              ihrt(1,2,1),ihrt(1,2,2),mint(1,2,1),mint(1,2,2), ! &30/11/90
     $                pred(1,2),                                     ! &30/11/90
     $              ihrt(2,1,1),ihrt(2,1,2),mint(2,1,1),mint(2,1,2), ! &30/11/90
     $                pred(2,1),                                     ! &30/11/90
     $              ihrt(2,2,1),ihrt(2,2,2),mint(2,2,1),mint(2,2,2), ! &30/11/90
     $                pred(2,2)                                      ! &30/11/90
   62             format(i5,14x,4i1,f6.2,1x,2(2x,4i1,f6.2))
                if(nform.eq.3)                                       ! &30/11/90
     $            write(nprout,63) idp,                              ! &30/11/90
     $              ihrt(1,1,1),ihrt(1,1,2),mint(1,1,1),mint(1,1,2), ! &30/11/90
     $                pred(1,1),                                     ! &30/11/90
     $              ihrt(2,1,1),ihrt(2,1,2),mint(2,1,1),mint(2,1,2), ! &30/11/90
     $                pred(2,1),                                     ! &30/11/90
     $              ihrt(2,2,1),ihrt(2,2,2),mint(2,2,1),mint(2,2,2), ! &30/11/90
     $                pred(2,2)                                      ! &30/11/90
   63             format(i5,2x,4i1,f6.2,13x,2(2x,4i1,f6.2))
                if(nform.eq.4)                                       ! &30/11/90
     $            write(nprout,64) idp,                              ! &30/11/90
     $              ihrt(1,1,1),ihrt(1,1,2),mint(1,1,1),mint(1,1,2), ! &30/11/90
     $                pred(1,1),                                     ! &30/11/90
     $              ihrt(1,2,1),ihrt(1,2,2),mint(1,2,1),mint(1,2,2), ! &30/11/90
     $                pred(1,2),                                     ! &30/11/90
     $              ihrt(2,2,1),ihrt(2,2,2),mint(2,2,1),mint(2,2,2), ! &30/11/90
     $                pred(2,2)                                      ! &30/11/90
   64             format(i5,2(2x,4i1,f6.2),15x,4i1,f6.2)
                if(nform.eq.5)                                       ! &30/11/90
     $            write(nprout,65) idp,                              ! &30/11/90
     $              ihrt(1,1,1),ihrt(1,1,2),mint(1,1,1),mint(1,1,2), ! &30/11/90
     $                pred(1,1),                                     ! &30/11/90
     $              ihrt(1,2,1),ihrt(1,2,2),mint(1,2,1),mint(1,2,2), ! &30/11/90
     $                pred(1,2),                                     ! &30/11/90
     $              ihrt(2,1,1),ihrt(2,1,2),mint(2,1,1),mint(2,1,2), ! &30/11/90
     $                pred(2,1)                                      ! &30/11/90
   65             format(i5,2(2x,4i1,f6.2),3x,4i1,f6.2,12x)
              endif                                                  ! @30/11/90
              idp=idy                                      ! Reset day &30/11/90
            endif                                                    ! @30/11/90
            do ii=1,2                                                ! &30/11/90
              do jj=1,2                                              ! &30/11/90
                ihrt(ii,jj,1)=-1                                     ! &30/11/90
              end do                                                 ! @30/11/90
            end do                                                   ! @30/11/90
            icount=0 !Mod. to count no. of estims. of HW or LW per day &30/11/90
            if(mon.ne.monp) then                                     ! &30/11/90
              if(iye.ne.iyp) then                       ! Test for end &30/11/90
                amnhat=amon(monhat)               ! Output HAT and LAT &30/11/90
                amnlat=amon(monlat)
                write(nprout,68) hthat,ihrhat(1),ihrhat(2),minhat(1),! &30/11/90
     $                           minhat(2),idyhat,amnhat,iyp,        ! &30/11/90
     $                           htlat,ihrlat(1),ihrlat(2),minlat(1),! &30/11/90
     $                           minlat(2),idylat,amnlat,iyp         ! &30/11/90
   68           format('1',////' HAT = ',f5.2,' at ',4i1,i3,1x,a5,i5 ! &30/11/90
     $                       //' LAT = ',f5.2,' at ',4i1,i3,1x,a5,i5)! &30/11/90
                go to 49
              endif                                                  ! @30/11/90
              if(alen.ne.'s') then               ! Update month &9/1/86 &4/12/90
                write(nprout,69) head1(1:55),head1(56:80)              ! &9/1/86
   69           format(1h1///1x,a55/1x,a25)                            ! &9/1/86
                write(nprout,70) head2(1:55),head2(56:80)              ! &9/1/86
   70           format(/1x,a55/1x,a25)                                 ! &9/1/86
                write(nprout,71) amon(mon),iye
   71           format(/22x,a5,i5/13x,'HIGH WATER',15x,'LOW WATER'/  ! &30/11/90
     $                 '  DATE',2(2(' TIME   HT. '),1x)/)            ! &30/11/90
              endif                                                  ! @30/11/90
              monp=mon                                               ! &30/11/90
            endif                                                    ! @30/11/90
C     Mod. to count number of estimates of HW or LW per day :
   59       icount=icount+1
            if(icount.gt.4) then         ! Abort if too many estimates &30/11/90
              write(nopout,107)                                      ! @30/11/90
              write(nlog,107)                                        ! @30/11/90
  107         format(/' More than 4 estimates of HW or LW in one day'/!@30/11/90
     $                ' Tides are probably diurnal -',               ! @30/11/90
     $                ' use hourly predictions instead',             ! @30/11/90
     $                ' ..... program terminated'/)                  ! @30/11/90
              stop                                                   ! @30/11/90
            endif                                                    ! @30/11/90
            height=tide(t,h,g,j,vpv)+pmean
            itype=1                        ! Find if high or low water &30/11/90
            if(height.lt.pmean) itype=2
            i=ihr/12+1
            ii=ihr/10
            ihrt(itype,i,1)=ii
            ihrt(itype,i,2)=ihr-10*ii
            ii=min/10
            mint(itype,i,1)=ii
            mint(itype,i,2)=min-10*ii
            pred(itype,i)=height
            if(height.gt.hthat) then                                 ! &30/11/90
              hthat=height
              monhat=mon
              idyhat=idy
              ihrhat(1)=ihrt(itype,i,1)
              ihrhat(2)=ihrt(itype,i,2)
              minhat(1)=mint(itype,i,1)
              minhat(2)=mint(itype,i,2)
            else if(height.lt.htlat) then                            ! &30/11/90
              htlat=height
              monlat=mon
              idylat=idy
              ihrlat(1)=ihrt(itype,i,1)
              ihrlat(2)=ihrt(itype,i,2)
              minlat(1)=mint(itype,i,1)
              minlat(2)=mint(itype,i,2)
            endif                                                    ! @30/11/90
            tlst(itype)=t
            ii=3-itype
            t=tlst(ii)+dt
            go to 54
          else         ! Diurnal or mixed tide (section added 18/5/89) @30/11/90
            t=0.
            hthat=pmean
            htlat=pmean
            pred1(0)=tided(0.,h,g,j1,vpv,1)
   99       call etd(iye,mon,idy,ihr,min,sec,t)
            if(alen.eq.'l') then                                      ! &3/12/90
              if(idy.eq.1) then
                write(nprout,101) head1
  101           format(1h1,a80)
                write(nprout,11) head2
                write(nprout,102) amon(mon),iye
  102           format(/1x,a5,i5)
                nline=4
              endif
              if(nline.le.47) then !Assume no more than 12 tides / day &30/11/90
                write(nprout,110) idy,amon(mon)
  110           format(/1x,i2,1x,a5)
                nline=nline+2
              else
                write(nprout,111) idy,amon(mon)
  111           format(/'1',i2,1x,a5)
                nline=2
              endif
            endif
            do i=1,24
              tim=t+float(i)*3600.
              idum=tim/rpt+1. !Mod. to interp. J, VPV on regular basis &30/11/90
              if(iint.ne.idum) then
                iint=idum
                ridum=idum
                tdum=(ridum-0.5)*rpt
                do k=1,ncon
                  j(k)=(j1(k)*(year-tdum)+j2(k)*tdum)/year
                  vpv(k)=v0(k)+(v1(k)*(year-tdum)+v2(k)*tdum)/year
                end do
              end if
              pred1(i)=tided(tim,h,g,j,vpv,1)
            end do
            do i=1,24
              if((pred1(i-1).ge.0..and.pred1(i).lt.0.).or.
     $           (pred1(i-1).le.0..and.pred1(i).gt.0.)) then
                tim=t+float(i-1)*3600.+pred1(i-1)*3600.              ! &30/11/90 
     $                /(pred1(i-1)-pred1(i))                         ! &30/11/90
                do ii=1,limit                        ! Newton's method &30/11/90
                  idum=tim/rpt+1.
                  if(iint.ne.idum) then
                    iint=idum
                    ridum=idum
                    tdum=(ridum-0.5)*rpt
                    do k=1,ncon
                      j(k)=(j1(k)*(year-tdum)+j2(k)*tdum)/year
                      vpv(k)=v0(k)+(v1(k)*(year-tdum)+v2(k)*tdum)/year
                    end do
                  endif
                  diff=tided(tim,h,g,j,vpv,1)
                  if(abs(diff).lt.eps) go to 104
                  tim=tim-diff/tided(tim,h,g,j,vpv,2)
                end do
                write(nopout,80)                 ! Failure to converge &30/11/90
                write(nlog,80)                   ! Failure to converge @30/11/90
                stop
  104           if(tim.ge.t+float(i-1)*3600..and.                    ! &30/11/90
     $             tim.lt.t+float(i)*3600.) then                     ! &30/11/90
                  height=tide(tim,h,g,j,vpv)+pmean
                  call etd(iye,mon,idy,ihr,min,sec,tim)
                  ii=ihr/10
                  ihrt(1,1,1)=ii
                  ihrt(1,1,2)=ihr-10*ii
                  ii=min/10
                  mint(1,1,1)=ii
                  mint(1,1,2)=min-10*ii
                  if(alen.eq.'l') then                                ! &3/12/90
                    write(nprout,105) ihrt(1,1,1),ihrt(1,1,2),
     $                                mint(1,1,1),mint(1,1,2),
     $                                height
  105               format(4x,4i1,f6.2)
                    nline=nline+1
                  endif
                  if(height.gt.hthat) then
                    hthat=height
                    monhat=mon
                    idyhat=idy
                    ihrhat(1)=ihrt(1,1,1)
                    ihrhat(2)=ihrt(1,1,2)
                    minhat(1)=mint(1,1,1)
                    minhat(2)=mint(1,1,2)
                  else if(height.lt.htlat) then
                    htlat=height
                    monlat=mon
                    idylat=idy
                    ihrlat(1)=ihrt(1,1,1)
                    ihrlat(2)=ihrt(1,1,2)
                    minlat(1)=mint(1,1,1)
                    minlat(2)=mint(1,1,2)
                  endif
                endif
              endif
            end do
            pred1(0)=pred1(24)
            t=t+86400.
            if(mon.ne.12.or.idy.ne.31) go to 99
            amnhat=amon(monhat)                   ! Output HAT and LAT &30/11/90
            amnlat=amon(monlat)
            write(nprout,68) hthat,ihrhat(1),ihrhat(2),              ! &30/11/90
     $                             minhat(1),minhat(2),              ! &30/11/90
     $                             idyhat,amnhat,iye,                ! &30/11/90
     $                       htlat,ihrlat(1),ihrlat(2),              ! &30/11/90
     $                             minlat(1),minlat(2),              ! &30/11/90
     $                             idylat,amnlat,iye                 ! &30/11/90
          endif       ! End of semidiurnal / (diurnal or mixed) blocks @30/11/90
   49     continue                    ! End of mods of 18/5/89 &30/11/90 &8/4/95
        endif      ! End of (hourly or diagnostic) / (HW or LW) blocks @30/11/90
      end do                                        ! End of year loop @29/11/90
      stop
      end
C
      function tided(t,h,g,j,vpv,id)
C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Mods. to array declarations           JRH:04:12:1990      *
C *                                                                            *
C *   SOURCE      :  TPD1150.FVS                                               *
C *   ROUTINE NAME:  TIDED                                                     *
C *   TYPE        :  FUNCTION                                                  *
C *                                                                            *
C *   FUNCTION    :  Calculates the IDth derivative of the tidal prediction    *
C *                  (eg. for ID=0, the result is the same as TIDE, but less   *
C *                  efficient). This is a modified version of the subroutine  *
C *                  TIDE, and is the version for multiple computations (it    *
C *                  sets up array for computation of COSINE).                 *
C *                                                                            *
C ******************************************************************************
      integer pcon
      parameter(pcon=114)      ! >/= Maximum no. of constituents used # @4/12/90
      real j
      dimension h(pcon),g(pcon),j(pcon),vpv(pcon),sigma(pcon)
      dimension sn(2001)                                ! Dimension N+1 &4/12/90
      common sigma,sn,itest,rn,n2,n4,rn2,rn4,f2,f4,ncon,n
      if(itest.eq.1) go to 2
      itest=1
      n2=2*n
      n4=4*n
      rn2=n2
      rn4=n4
      f2=rn*4./360.
      f4=3.14159265/(2.*rn)
C     CALCULATE TIDE :
    2 tided=0.
      d=id
      f3=d*rn
      tf4=t/f4
      do 1 i=1,ncon
      if(j(i).eq.0.) go to 1
      arg=abs(sigma(i)*tf4+(vpv(i)-g(i))*f2+f3)
      arg=abs(amod(arg,rn4)-rn2)-rn
      iarg=abs(arg)+0.5
      tided=tided+h(i)*j(i)*(sigma(i)**id)*sign(sn(iarg+1),arg)
    1 continue
      return
      end
C
      function tide(t,h,g,j,vpv)
C ******************************************************************************
C *                                                                            *
C *                            VAX/SUN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  TIDES                             REF:JRH:09:01:1986      *
C *                                                                            *
C *   REVISION    :  Mods. to array declarations           JRH:04:12:1990      *
C *                                                                            *
C *   SOURCE      :  TPD1150.FVS                                               *
C *   ROUTINE NAME:  TIDE                                                      *
C *   TYPE        :  FUNCTION                                                  *
C *                                                                            *
C *   FUNCTION    :  Performs tidal prediction. This is the version for        *
C *                  multiple computations (it sets up array for computation   *
C *                  of COSINE).                                               *
C *                                                                            *
C ******************************************************************************
      integer pcon
      parameter(pcon=114)      ! >/= Maximum no. of constituents used # @4/12/90
      real j
      dimension h(pcon),g(pcon),j(pcon),vpv(pcon),sigma(pcon)
      dimension sn(2001)                                ! Dimension N+1 &4/12/90
      common sigma,sn,itest,rn,n2,n4,rn2,rn4,f2,f4,ncon,n
      if (itest.eq.1) goto 2
      itest=1
      n2=2*n
      n4=4*n
      rn2=n2
      rn4=n4
      f2=rn*4./360.
      f4=3.14159265/(2.*rn)
C     CALCULATE TIDE :
    2 tide=0.
      tf4=t/f4
      do 1 i=1,ncon
      if(j(i).eq.0.) go to 1
      arg=abs(sigma(i)*tf4+(vpv(i)-g(i))*f2)
      arg=abs(amod(arg,rn4)-rn2)-rn
      iarg=abs(arg)+0.5
      tide=tide+h(i)*j(i)*sign(sn(iarg+1),arg)
    1 continue
      return
      end
C
      subroutine aread(nin,nout,nlog,anot,acc,aarray,istart,istop,form)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:04:09:1986      *
C *                                                                            *
C *   REVISION    :  Mod. for output of message only if                        *
C *                  ISTART=0 or ISTOP=0                   JRH:12:09:1986      *
C *                  Minor mods.                           JRH:31:12:1986      *
C *                  Minor mod.                            JRH:05:01:1987      *
C *                  Minor mod.                            JRH:17:05:1988      *
C *                  Mod. so that FORM is input without                        *
C *                  brackets                              JRH:21:06:1988      *
C *                  Mod. for compilation on transputer    JRH:07:09:1989      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  AREAD                                                     *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Reads in a set of CHARACTER*40 variables                  *
C *                                                                            *
C *                  NIN ..... Operator input device                           *
C *                  NOUT .... Operator output device                          *
C *                  NLOG .... Log device                                      *
C *                  ANOT .... Annotation (40 chars.)                          *
C *                  ACC ..... Carriage-control character                      *
C *                            (' ' for new line, '$' for same line)           *
C *                  AARRAY .. Resultant CHARACTER*40 array                    *
C *                  ISTART .. Start index of AARRAY                           *
C *                  ISTOP ... Stop index of AARRAY                            *
C *                  FORM .... Format for output to log device (40 chars.,     *
C *                            without brackets)                               *
C *                                                                            *
C ******************************************************************************
      character*40 aarray(*)
      character*40 anot,form
      character*1 acc
      character*42 formt                                               ! @7/9/89
      logical cont                                                   ! @31/12/86
      save cont                                                      ! @31/12/86
      character*12 formout                                             ! &5/1/87
      data cont/.false./                                              ! @11/5/88
      if(cont) then                    ! Continuation of terminal text @31/12/86
        write(formout,1) ' ',acc                                     ! &31/12/86
    1   format('(',a1,'X,A40,A1',a1,')')                               ! &5/1/87
      else                                               ! Normal text @31/12/86
        write(formout,1) '/',acc                                     ! @31/12/86
      endif                                                          ! @31/12/86
      write(nout,formout) anot,' '                                     ! &5/1/87
      write(nlog,formout) anot,' '                                     ! &5/1/87
      if(istart.ne.0..and.istop.ne.0) then                           ! &31/12/86
        read(nin,2) (aarray(i),i=istart,istop)
    2   format((a40))
        formt='('//form//')'                                           ! @7/9/89
        write(nlog,formt) (aarray(i),i=istart,istop)          ! &21/6/88 &7/9/89
        cont=.false.               ! Terminal text not to be continued @31/12/86
      else                                                           ! @31/12/86
        cont=.true.                    ! Terminal text to be continued @31/12/86
      endif                                                           ! @12/9/86
      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
      subroutine etd(iye,mon,idy,ihr,min,sec,tim)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  GENERAL                           REF:JRH:04:07:1985      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1985      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  ETD                                                       *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Converts elapsed time to daytime.                         *
C *                  Inputs are IYE and TIM.                                   *
C *                  IYE and TIM are modified if TIM is negative or            *
C *                  greater than year length.                                 *
C *                                                                            *
C *                  Daytime :                                                 *
C *                    IYE ..... Year                                          *
C *                    MON ..... Month                                         *
C *                    IDY ..... Day                                           *
C *                    IHR ..... Hour                                          *
C *                    MIN ..... Minute                                        *
C *                    SEC ..... Second                                        *
C *                                                                            *
C *                  Elapsed time :                                            *
C *                    TIM ..... Time in seconds from beginning of year        *
C *                                                                            *
C *                  NOTE that current version (with D.P. (64 bit)             *
C *                  internal computations, but S.P. (32 bit) TIM) will        *
C *                  generally be correct to +/- 2 secs. However, a time       *
C *                  in exact minutes (ie. SEC = 0.) is computed exactly.      *
C *                                                                            *
C ******************************************************************************
      double precision dtim,dyl,didy,dihr,dmin
      dimension nd(13,2)
      data nd/0,31,59,90,120,151,181,212,243,273,304,334,365,
     $        0,31,60,91,121,152,182,213,244,274,305,335,366/
      dtim=dble(tim)
      if(dtim.lt.0.d0) go to 3
    5 l=1
C     Leap year test .....
      if((iye/4)*4.eq.iye) l=2
      dyl=dble(float(nd(13,l)))*86400.d0
      if(dtim.lt.dyl) go to 4
      iye=iye+1
      dtim=dtim-dyl
      go to 5
    3 iye=iye-1
      l=1
C     Leap year test .....
      if((iye/4)*4.eq.iye) l=2
      dyl=dble(float(nd(13,l)))*86400.d0
      dtim=dtim+dyl
      if(dtim.lt.0.d0) go to 3
    4 dmin=dint(dtim/60.d0)
      sec=sngl(dtim-dmin*60.d0)
      dihr=dint(dmin/60.d0)
      min=idint(dmin-dihr*60.d0+0.5d0)
      didy=dint(dihr/24.d0)
      ihr=idint(dihr-didy*24.d0+0.5d0)
      idy=idint(didy+1.5d0)
      do 1 i=1,13
      if(nd(i,l).ge.idy) go to 2
    1 continue
    2 mon=i-1
      idy=idy-nd(mon,l)
      tim=sngl(dtim)
      return
      end
C
      subroutine iread(nin,nout,nlog,anot,acc,iarray,istart,istop,form)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:27:08:1986      *
C *                                                                            *
C *   REVISION    :  Error trap included                   JRH:04:09:1986      *
C *                  Mod. for output of message only if                        *
C *                  ISTART=0 or ISTOP=0                   JRH:12:09:1986      *
C *                  Minor mods.                           JRH:31:12:1986      *
C *                  Minor mod.                            JRH:05:01:1987      *
C *                  Minor mod.                            JRH:17:05:1988      *
C *                  Mod. so that FORM is input without                        *
C *                  brackets                              JRH:21:06:1988      *
C *                  Mod. for null read after error to cope                    *
C *                  with Sun bug                          JRH:25:07:1989      *
C *                  Mod. for compilation on transputer    JRH:07:09:1989      *
C *                                                                            *
C *   SOURCE      :  FORLIB.FOR                                                *
C *   ROUTINE NAME:  IREAD                                                     *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Reads in a set of INTEGER numbers                         *
C *                                                                            *
C *                  NIN ..... Operator input device                           *
C *                  NOUT .... Operator output device                          *
C *                  NLOG .... Log device                                      *
C *                  ANOT .... Annotation (40 chars.)                          *
C *                  ACC ..... Carriage-control character                      *
C *                            (' ' for new line, '$' for same line)           *
C *                  IARRAY .. Resultant INTEGER array                         *
C *                  ISTART .. Start index of IARRAY                           *
C *                  ISTOP ... Stop index of IARRAY                            *
C *                  FORM .... Format for output to log device (40 chars.,     *
C *                            without brackets)                               *
C *                                                                            *
C ******************************************************************************
      integer iarray(*)
      character*40 anot,form
      character*1 acc
      character*42 formt                                               ! @7/9/89
      logical cont                                                   ! @31/12/86
      save cont                                                      ! @31/12/86
      character*12 formout                                            ! &5/11/86
      data cont/.false./                                              ! @11/5/88
    2 if(cont) then           ! Continuation of terminal text @31/12/86 &25/7/89
        write(formout,1) ' ',acc                                     ! &31/12/86
    1   format('(',a1,'X,A40,A1',a1,')')                              ! &5/11/87
      else                                               ! Normal text @31/12/86
        write(formout,1) '/',acc                                     ! @31/12/86
      endif                                                          ! @31/12/86
      write(nout,formout) anot,' '                                    ! &5/11/87
      write(nlog,formout) anot,' '                                    ! &5/11/87
      if(istart.ne.0.and.istop.ne.0) then                            ! &31/12/86
        read(nin,*,err=3) (iarray(i),i=istart,istop)                   ! &4/9/86
        formt='('//form//')'                                           ! @7/9/89
        write(nlog,formt) (iarray(i),i=istart,istop)          ! &21/6/88 &7/9/89
        cont=.false.               ! Terminal text not to be continued @31/12/86
      else                                                           ! @31/12/86
        cont=.true.                    ! Terminal text to be continued @31/12/86
      endif                                                           ! @12/9/86
      return
    3 read(nin,*)                      ! Null read to cope with Sun bug @25/7/89
      go to 2                                                         ! @25/7/89
      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
