C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:15:01:1997      *
C *                                                                            *
C *   REVISION    :  Addition of time zones                JRH:17:04:1998      *
C *                  Inputs are now command line arguments JRH:17:04:1998      *
C *                  tz_t changed to real*8                JRH:10:09:1999      *
C *                  Addition of header re. software                           *
C *                  licensing                             JRH:23:12:2002      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  tpd1450                                                   *
C *   TYPE        :  MAIN                                                      *
C *                                                                            *
C *   FUNCTION    :  General tide prediction subroutine and test 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 ******************************************************************************
C
      implicit none
C
      real*8 tide_pred,tide
      real*8 tz_t
C
      real*4 sec
C
      integer*4 iye,mon,idy,ihr,min,ifail
      integer*4 n_arg,iargc
      integer*4 ios
C
      character*120 harm_const,buf
C
      n_arg=iargc()                           ! Number of command line arguments
      if(n_arg.ne.8) call error_handler(1,001)                 ! Error point 001
C
      call getarg(1,buf)
      read(buf,*,iostat=ios) iye
      call error_handler(ios,002)                              ! Error point 002
C
      call getarg(2,buf)
      read(buf,*,iostat=ios) mon
      call error_handler(ios,003)                              ! Error point 003
C
      call getarg(3,buf)
      read(buf,*,iostat=ios) idy
      call error_handler(ios,004)                              ! Error point 004
C
      call getarg(4,buf)
      read(buf,*,iostat=ios) ihr
      call error_handler(ios,005)                              ! Error point 005
C
      call getarg(5,buf)
      read(buf,*,iostat=ios) min
      call error_handler(ios,006)                              ! Error point 006
C
      call getarg(6,buf)
      read(buf,*,iostat=ios) sec
      call error_handler(ios,007)                              ! Error point 007
C
      call getarg(7,buf)
      read(buf,*,iostat=ios) tz_t
      call error_handler(ios,008)                              ! Error point 008
C
      call getarg(8,harm_const)
C
      tide=tide_pred(iye,mon,idy,ihr,min,sec,tz_t,
     $               harm_const,
     $               ifail)
C
      write(6,1) tide,ifail
    1 format(' Tidal prediction: ',d23.15,' ifail: ',i5)
C
      stop
C
      end
C     
      real*8 function tide_pred(iye,mon,idy,ihr,min,sec,tz_t,
     $                          harm_const,
     $                          ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:22:01:1997      *
C *                                                                            *
C *   REVISION    :  Mod.to allow comments (prefix "#") in                     *
C *                  input file                            JRH:17:04:1998      *
C *                  Addition of time zones                JRH:17:04:1998      *
C *                  Interpolation of j and v based on GMT JRH:17:04:1998      *
C *                  tz_t changed to real*8                JRH:10:09:1999      *
C *                  Mod. to calc. of interp. weights      JRH:10:09:1999      *
C *                  Mod. to ignore "*" in harmonic                            *
C *                  constant file                         JRH:13:09:1999      *
C *                  Addition of many "save"s              JRH:31:07:2000      *
C *                  Check Z0 occurs once and only once    JRH:31:07:2000      *
C *                  Mod. to allow any year and both "98"                      *
C *                  and "1998" type formats               JRH:06:09:2001      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  tide_pred                                                 *
C *   TYPE        :  real*8 function                                           *
C *                                                                            *
C *   FUNCTION    :  General tide prediction subroutine                        *
C *                                                                            *
C *                  iye ........ year (1970-2030)                             *
C *                  mon ........ month (1-12)                                 *
C *                  idy ........ day of month (1-31)                          *
C *                  ihr ........ hour of day (0-24)                           *
C *                  min ........ minute of hour (0-60)                        *
C *                  sec ........ second of minute (0-60)                      *
C *                  tz_t ....... time zone of time (hours, positive eastwards)*
C *                  harm_const . file of harmonic constants, in format:       *
C *                               <constituent name> <h> <g> as (a6,x,*,*)     *
C *                               NOTE THAT YOU HAVE TO LEAVE SIX COLUMNS FOR  *
C *                               THE CONSTITUENT NAME                         *
C *                  ifail ...... 0 for successful execution                   *
C *                                                                            *
C *                  The location of the directory of astronomical arguments   *
C *                  is given in the file dot_tides in the local directory.    *
C *                  The value of the variable dot_tides should be selected    *
C *                  appropriately using the data statements supplied.         *
C *                                                                            *
C *                  NOTE that this subroutine uses input device numbers 7, 9  *
C *                  and 11.                                                   *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      real*8 tz_t
      real*4 sec
      integer*4 iye,mon,idy,ihr,min,ifail
      character*120 harm_const
      character*8 dot_tides
C
C     Local variables:
C
      integer*4 pcon,ptry
      parameter(pcon=115,  ! >/= Maximum no. of constituents used (including Z0)
     $          ptry=3)                 ! No. of tries at different const. names
C
      real*8 h(pcon),g(pcon)
      real*8 j1(pcon),v1(pcon),vpv(pcon)
      real*8 j2(pcon),v2(pcon),v0(pcon)
      real*8 sigma(pcon)
C
      real*8 j_int,vpv_int
      real*8 j_day,j_day_1,j_day_2,tim
      real*8 tz_g
      real*8 d2r
      real*8 dum,year
      real*8 w1,w2
C
      integer*4 iye_ast_arg_old
      integer*4 i,i_z0,n_z0,ncon,k,nyear
      integer*4 ios,idum,itry
      integer*4 len,len_a,lenchar
C
      logical first,found,match,l_tz_g
C
      character*6 acon(pcon)
      character*6 atry(ptry,2)
C
      character*1 slash
      character*6 filnam
      character*40 buf
      character*80 formt
      character*120 harm_const_old,ast_arg
C
      save iye_ast_arg_old,harm_const_old,first
      save d2r,ast_arg,len_a,slash,tz_g,l_tz_g,ncon,acon,h,g,i_z0
      save sigma,j1,v1,vpv,j2,v2,v0
C
      data dot_tides/'.tides  '/                                          ! Unix
C     data dot_tides/'ARGS.DIR'/                                           ! DOS
      data iye_ast_arg_old/0/
      data harm_const_old( 1:30 )/'                              '/
      data harm_const_old(31:60 )/'                              '/
      data harm_const_old(61:90 )/'                              '/
      data harm_const_old(91:120)/'                              '/
C
      data atry(1,1),atry(1,2)/'sig1  ','sigma1'/
      data atry(2,1),atry(2,2)/'the1  ','theta1'/
      data atry(3,1),atry(3,2)/'lam2  ','lamda2'/
C
      data first/.true./
      data l_tz_g/.false./
C
      if(first) then
C
        d2r=datan(1.d0)*4.d0/180.d0                         ! Degrees to radians
C
C     Find directory for astronomical arguments:
C
        open(unit=9,file=dot_tides,status='old',iostat=ios)
        if(ios.ne.0) then
          ifail=1
          return
        endif
        read(9,3,iostat=ios) ast_arg      ! Directory for astronomical arguments
    3   format(a120)
        if(ios.ne.0) then
          ifail=2
          return
        endif
        len_a=lenchar(ast_arg)
        if(dot_tides.eq.'.tides  ') then                    ! Set path delimiter
          slash='/'                                                       ! Unix
        else
          slash=char(92)                                       ! (backslash) DOS
        endif
C
        first=.false.
C
      endif
C
C     Check year is in range (removed 6/9/2001):
C
C     if(iye.lt.1970.or.iye.gt.2030) then
C       ifail=3
C       return
C     endif
C
C     Read in harmonic constants, if necessary:
C
      if(harm_const.ne.harm_const_old) then
C
        harm_const_old=harm_const              ! Remember old harmonic constants
C
C       write(6,4)                                                  ! Diagnostic
C   4   format('Opening file of harmonic constants')                ! Diagnostic
        open(unit=7,file=harm_const,status='old',iostat=ios)
        if(ios.ne.0) then
          ifail=4
          return
        endif
        ncon=0
        n_z0=0
        ios=0
        do while(ios.eq.0)
          read(7,1,iostat=ios) buf
    1     format(a40)
          if(ios.gt.0) then
            ifail=5
            return
          endif
          if(ios.eq.0) then
            if(buf(1:1).eq.'#') then                                 ! A comment
              if(buf(1:5).eq.'# TZ ') then                 ! Found the time zone
                read(buf(6:40),*,iostat=ios) tz_g
                if(ios.ne.0) then
                  ifail=6
                  return
                endif
                l_tz_g=.true.
              endif
            else                                            ! Harmonic constants
              len=index(buf,' ')-1       ! Find separator between name and (h,g)
              if(len.gt.6) len=6     ! Assume format is a6 and there is no space
              ncon=ncon+1
              acon(ncon)='      '                                   ! For safety
              write(formt,'("(a",i3,")")') len
              read(buf,formt,iostat=ios) acon(ncon)(1:len)
              if(ios.ne.0) then
                ifail=7
                return
              endif
C
C     Change any "*" to space:
C
              do i=len+1,40
                if(buf(i:i).eq.'*') buf(i:i)=' '
              end do
C
              if(acon(ncon).eq.'z0    '.or.
     $           acon(ncon).eq.'Z0    ') then
                read(buf(len+1:40),*,iostat=ios) h(ncon)
                n_z0=n_z0+1
              else
                read(buf(len+1:40),*,iostat=ios) h(ncon),g(ncon)
              endif
              if(ios.ne.0) then
                ifail=8
                return
              endif
              len=lenchar(acon(ncon))
              call changecase(acon(ncon),acon(ncon),len,0) !Change to lower case
            endif
          endif
        end do
        close(unit=7)
C
        if(n_z0.ne.1) then                 ! Z0 did not occur once and once only
          ifail=9
          return
        endif
C
      endif
C
C     Read in astronomical arguments, if necessary:
C
      if(iye.ne.iye_ast_arg_old) then
C
        iye_ast_arg_old=iye                                  ! Remember old year
C
        do i=1,ncon
C
          if(acon(i).eq.'z0    ') then                           ! Trap mean, Z0
            i_z0=i
          else
            ios=0
            itry=0
            found=.false.
            do while(.not.found)
              len=lenchar(acon(i))
              filnam(1:len)=acon(i)(1:len)
C             write(6,5) filnam(1:len)                              ! Diagnostic
C   5         format('Opening ast. arg. file for ',a6)              ! Diagnostic
              open(unit=11,
     $             file=ast_arg(1:len_a)//slash//filnam(1:len)//'.con',
     $             status='old',iostat=ios)
              found=(ios.eq.0)                             ! .true. if open O.K.
              match=(ios.eq.0)               ! If open not O.K., try other names
              do while(.not.match.and.itry.lt.ptry)
                itry=itry+1
                if(acon(i).eq.atry(itry,1)) then
                  acon(i)=atry(itry,2)
                  match=.true.
                else if(acon(i).eq.atry(itry,2)) then
                  acon(i)=atry(itry,1)
                  match=.true.
                endif
              end do
              if(.not.found.and..not.match) then
                ifail=10
                return
              endif
            end do
C
            read(11,*,iostat=ios) dum
            if(ios.ne.0) then
              ifail=11
              return
            endif
            sigma(i)=dum/3600.d0                 ! Convert to degrees per second
            nyear=iye-1900             ! Note this format only works around 2000
            if(nyear.ge.100) nyear=nyear-100
            k=9999
            do while(k.ne.nyear.and.k.ne.iye)                        ! &6/9/2001
              read(11,1,iostat=ios) buf
              if(ios.ne.0) then
                ifail=12
                return
              endif
              read(buf,*,iostat=ios) k
              if(ios.ne.0) then
                ifail=13
                return
              endif
              if(k.eq.nyear.or.k.eq.iye) then                        ! &6/9/2001
                read(buf,*,iostat=ios) idum,j1(i),v1(i),vpv(i)
                if(ios.ne.0) then
                  ifail=14
                  return
                endif
                read(11,*,iostat=ios) idum,j2(i),v2(i)
                if(ios.ne.0) then
                  ifail=15
                  return
                endif
              endif
            end do
            close(unit=11)
          endif
        end do
C
        do i=1,ncon
          v0(i)=vpv(i)-v1(i)
        end do
C
      endif
C
C     Assume tz_g is the same as tz_t if not given in input file:
C
      if(.not.l_tz_g) tz_g=tz_t
C
C     Convert to time from start of year:
C
      call julday(iye,mon,idy,ihr,min,sec,j_day  ,ifail)
      if(ifail.ne.0) then
        ifail=16
        return
      endif
C
      call julday(iye,1  ,1  ,0  ,0  ,0. ,j_day_1,ifail)
      if(ifail.ne.0) then
        ifail=17
        return
      endif
C
      tim=(j_day-j_day_1)*86400.d0
C
C     Find length of current year:
C
      call julday(iye+1,1  ,1  ,0  ,0  ,0. ,j_day_2,ifail)
      if(ifail.ne.0) then
        ifail=18
        return
      endif
C
      year=(j_day_2-j_day_1)*86400.d0
C
C     Linearly interpolate nodal corrections through year,
C     and calculate tide:
C
      tide_pred=h(i_z0)                                               ! Mean, Z0
      w2=(tim-tz_t*3600.d0)/year                                   ! Time in GMT
      w1=1.d0-w2
      do i=1,ncon
        if(i.ne.i_z0) then
          j_int=j1(i)*w1+j2(i)*w2
          vpv_int=v0(i)+v1(i)*w1+v2(i)*w2
          tide_pred=tide_pred
     $      +h(i)*j_int*dcos((sigma(i)*(tim+(tz_g-tz_t)*3600.d0)
     $      +vpv_int-g(i))*d2r)
        endif
      end do
C
      ifail=0
      return
      end
C
      subroutine changecase(ain,aout,n,itype)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:03:12:1990      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1990      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  changecase                                                *
C *   TYPE        :  subroutine                                                *
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 ******************************************************************************
C
      implicit none
C
      integer*4 n,itype
      character*(*) ain,aout
C
      integer*4 idum,i
C
      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 error_handler(ifail,error_point)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  MODELLING                         REF:JRH:06:01:1995      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1995      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  error_handler                                             *
C *   TYPE        :  SUBROUTINE                                                *
C *                                                                            *
C *   FUNCTION    :  Handles error (ifail.ne.0)                                *
C *                                                                            *
C ******************************************************************************
      integer*4 ifail,error_point
C
      if(ifail.ne.0) then
        write(6,1) ifail,error_point
    1   format(/' ifail returned as ',i5,' at error point ',i5,
     $          ' ..... program terminated'/)
        stop
      endif
      return
      end
C
      integer*4 function idint2(x)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:23:04:1991      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1991      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  idint2                                                    *
C *   TYPE        :  integer*4 function                                        *
C *                                                                            *
C *   FUNCTION    :  As INT, but trunctates towards -(infinity)                *
C                    This is D.P. version of INT2                              *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      real*8 x
      idint2=idint(x)
      if(dble(idint2).eq.x) return
      if(x.lt.0.d0) idint2=idint2-1
      return
      end 
C
      subroutine julday(iye,mon,idy,ihr,min,sec,julian,ifail)
C ******************************************************************************
C *                                                                            *
C *                            FORTRAN SOURCE CODE                             *
C *                                                                            *
C *   PROGRAM SET :  UTILITY                           REF:JRH:23:05:1991      *
C *                                                                            *
C *   REVISION    :  -------------                         JRH:--:--:1991      *
C *                                                                            *
C *   SOURCE      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  julday                                                    *
C *   TYPE        :  subroutine                                                *
C *                                                                            *
C *   FUNCTION    :  Converts (year, month, day, hour, minute, second) to      *
C *                  Julian Day (based on Numerical Recipes)                   *
C *                                                                            *
C *                  IYE ........ Year                                         *
C *                  MON ........ Month                                        *
C *                  IDY ........ Day                                          *
C *                  IHR ........ Hour                                         *
C *                  MIN ........ Minute                                       *
C *                  SEC ........ Second                                       *
C *                  JULIAN ..... Julian Day (non-negative, but may be         *
C *                               non-integer)                                 *
C *                  IFAIL ...... 0 for successful execution, otherwise 1      *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      real*8 julian
      real*4 sec
      integer*4 iye,mon,idy,ihr,min,ifail
      integer*4 iyyy,jy,jm,igreg,ja,ijul
      integer*4 idint2
      parameter (igreg=15+31*(10+12*1582))
C                              ..... Gregorian Calendar was adopted 15 Oct. 1582
      if(iye.eq.0.or.                                    ! There is no year zero
     $   iye.lt.-4713) then                    ! Julian Day must be non-neagtive
        ifail=1
        return
      endif
      if(iye.lt.0) then
        iyyy=iye+1
      else
        iyyy=iye
      endif
      if(mon.gt.2) then
        jy=iyyy
        jm=mon+1
      else
        jy=iyyy-1
        jm=mon+13
      endif
      ijul=idint2(365.25d0*dble(jy))+idint2(30.6001d0*dble(jm))
     $       +idy+1720995
      if(idy+31*(mon+12*iyyy).ge.igreg) then
C                                    ..... Test for change to Gregorian Calendar
        ja=idint(0.01d0*dble(jy))
        ijul=ijul+2-ja+idint(0.25d0*dble(ja))
      endif
      julian=dble(ijul)
     $         +dble(ihr)/24.d0+dble(min)/1440.d0+dble(sec)/86400.d0
     $         -0.5d0                         ! Correction from midnight to noon
      if(julian.lt.0.d0) then                  ! Julian Day must be non-negative
        ifail=1
        return
      endif
      ifail=0
      return
      end
C
      integer*4 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      :  tpd1450.f                                                 *
C *   ROUTINE NAME:  lenchar                                                   *
C *   TYPE        :  integer*4 function                                        *
C *                                                                            *
C *   FUNCTION    :  Returns length of CHARACTER  variable (defined by         *
C *                  removing blank characters from right-hand side).          *
C *                                                                            *
C ******************************************************************************
C
      implicit none
C
      character*(*) c
C
      integer*4 itot,i
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
