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 * 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