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