C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : TIDES REF:JRH:13:09:1999 * C * * C * REVISION : Removal of non-standard "<>" construct * C * in formats JRH:17:05:2000 * C * Minor mods. to format JRH:17:05:2000 * C * Mod. to allow any year JRH:07:09:2001 * C * All necessary subroutines added to * C * source code JRH:14:09:2001 * C * Change of format for output of Z0 JRH:18:09:2001 * C * Modified formats to cope with "bug" in * C * latest g77 -- on reaching the end of * C * a format statement, it does return to * C * second outermost bracket BUT * C * (incorrectly, I think) inserts a * C * carriage return. The old statements * C * have been commented out for retrieval * C * later. JRH:01:11:2002 * C * Addition of header re. software * C * licensing JRH:23:12:2002 * C * * C * SOURCE : tda1537.f * C * ROUTINE NAME: tda1537 * C * TYPE : MAIN * C * * C * FUNCTION : Latest tidal analysis program (based on tda1154 and * C * tda1527). Major updates are: * C * * C * 1. It uses an input parameter file. * C * 2. It uses a "dot_tides" file to point to the directory * C * of the astronomical arguments, constituent lists and * C * harmonic constants for reference port. * C * 3. It has had potential Y2K problems fixed. * C * 4. It has been changed to real*8. * C * 5. It uses SVD * C * 6. It handles time zones properly * C * * C * The location of the directory of astronomical arguments, * C * constituent lists and harmonic constants for the reference* C * port 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 ****************************************************************************** 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 integer*4 pcon,pequ,pset,phead C parameter(pcon=114, ! >/= Maximum no. of constituents used $ pequ=pcon*2+1, ! >/= (maximum no. of constituents)*2+1 $ pset=20, ! >/= Maximum number of input data sets $ phead=10) ! >/= Maximum number of header lines C real*8 j1(pcon),j2(pcon) real*8 v1(pcon),v2(pcon) real*8 v0(pcon),vpv(pcon) real*8 j(pcon),sigma(pcon) real*8 hrel(pcon),grel(pcon) C real*8 a(pequ,pequ),acov(pequ,pequ),x(pequ) real*8 usvd(pequ,pequ),vsvd(pequ,pequ),wsvd(pequ) real*8 e(pequ),tmpsvd(pequ) real*8 rmean(pset),dat(pset),res(pset),svar(pset) real*8 var(pset) real*8 sres(pset),r1(pset),r2(pset) real*8 h(pcon,pset),g(pcon,pset) real*8 b(pequ,pset),bcov(pequ,pset),c(pequ,pset) real*8 b1(pequ),c1(pequ) C real*8 resiset real*8 jul_start,jul_year,year real*8 tz_t,tz_g,tz_out,del_tz real*8 t,tlast real*8 arg real*8 r,delta real*8 rc,rs real*8 w1,w2 real*8 d2r C real*4 sec C integer*4 aisp(pcon),index_ref(pcon) C integer*4 aispi integer*4 nhead,ihead integer*4 ncon,icon,nref,iref,isp integer*4 ndum integer*4 nequ,nequm1,iequ1,iequ2,ndat,idat,nset,iset integer*4 len,lenchar,lentid,lhschar integer*4 ilhs integer*4 indexi integer*4 iye,iye_start,iye_last,mon,idy,ihr,min integer*4 ic,ics,icr,is,iss,isr integer*4 ichar integer*4 job,info integer*4 ios,ifail C logical logs,logr,lres,lyearchange,ltrep logical end_of_cons,end_of_file logical l_tz_t,l_tz_g C character*120 head(phead) character*6 asplit,aref,atest,acon(pcon),arel(pcon) C character*120 parameter_file,conslist_file,input_file character*120 hgout_file,resout_file,temp_file,port_ref character*120 conbuf,time_type,datbuf,tidal_dir character*80 formt character*40 refbuf character*8 dot_tides character*1 astar,tf character*1 slash C data dot_tides/'.tides '/ ! Unix C data dot_tides/'ARGS.DIR'/ ! DOS C d2r=4.d0*datan(1.d0)/180.d0 ! Degrees to radians C C Find directory for astronomical arguments, constituent lists and C harmonic constants for the reference port: C open(unit=7,file=dot_tides,status='old',iostat=ios) call error_handler(ios,001) ! Error point 001 C read(7,8,iostat=ios) tidal_dir ! Required directory 8 format(a120) call error_handler(ios,002) ! Error point 002 C close(unit=7) C if(dot_tides.eq.'.tides ') then ! Set path delimiter slash='/' ! Unix else slash=char(92) ! (backslash) DOS endif C lentid=lenchar(tidal_dir) tidal_dir=tidal_dir(1:lentid)//slash ! Add appropriate slash lentid=lenchar(tidal_dir) C C Read name of parameter file and open it: C read (5,1,iostat=ios) parameter_file 1 format(a120) call error_handler(ios,003) ! Error point 003 C open(unit=7,file=parameter_file,status='old', $ iostat=ios) call error_handler(ios,004) ! Error point 004 C C Output file of residuals? C call areadfil3(7,'OPRESID',7,tf,ifail) call error_handler(ifail,005) ! Error point 005 C if(tf.ne.'t'.and.tf.ne.'T'.and.tf.ne.'f'.and.tf.ne.'F') $ call error_handler(1,006) ! Error point 006 lres=(tf.eq.'t'.or.tf.eq.'T') C C Filename for constituent lists: C call areadfil3(7,'CONSLIST_FILE',13,conslist_file,ifail) call error_handler(ifail,007) ! Error point 007 C open(unit=11,file=tidal_dir(1:lentid)//conslist_file, $ status='old',iostat=ios) call error_handler(ios,008) ! Error point 008 C C Filename for input data: C call areadfil3(7,'INPUT_FILE',10,input_file,ifail) call error_handler(ifail,009) ! Error point 009 C open(unit=13,file=input_file,status='old', $ iostat=ios) call error_handler(ios,010) ! Error point 010 C C Filename for (h,g) output data: C call areadfil3(7,'HGOUT_FILE',10,hgout_file,ifail) call error_handler(ifail,011) ! Error point 011 C open(unit=8,file=hgout_file,status='unknown', $ iostat=ios) call error_handler(ios,012) ! Error point 012 C if(lres) then C C Filename for residual output data: C call areadfil3(7,'RESOUT_FILE',11,resout_file,ifail) call error_handler(ifail,013) ! Error point 013 C open(unit=10,file=resout_file,status='unknown', $ iostat=ios) call error_handler(ios,014) ! Error point 014 C C Filename for temporary storage (sequential access, binary): C call areadfil3(7,'TEMP_FILE',9,temp_file,ifail) call error_handler(ifail,015) ! Error point 015 C open(unit=15,file=temp_file,status='unknown', $ access='sequential',form='unformatted', $ iostat=ios) call error_handler(ios,016) ! Error point 016 C endif C C Header: C call ireadfil2(7,'NHEAD',5,nhead,ifail) ! Number of header lines call error_handler(ifail,017) ! Error point 017 C read(7,2) (head(ihead),ihead=1,nhead) 2 format(a120) C C Time zone of time series (if present): C call dreadfil2(7,'TZ_T',4,tz_t,ifail) l_tz_t=(ifail.eq.0) C C Type of time coordinate ("REP": time repeats every year; C "NONREP": time is continuous from year to year): C call areadfil3(7,'TIME_TYPE',9,time_type,ifail) call error_handler(ifail,018) ! Error point 018 C if(time_type.ne.'REP'.and.time_type.ne.'rep'.and. $ time_type.ne.'NONREP'.and.time_type.ne.'nonrep') $ call error_handler(1,019) ! Error point 019 ltrep=(time_type.eq.'REP'.or.time_type.eq.'rep') C C Input constituent lists: C ncon=0 nref=0 end_of_cons=.false. end_of_file=.false. C do while(ios.eq.0.and..not.end_of_cons.and..not.end_of_file) C read(11,10,iostat=ios) conbuf 10 format(a120) C if(ios.eq.0) then ! Read is O.K. C if(conbuf(1:1).ne.'#') then ! Ignore comments C len=lenchar(conbuf) if(len.le.6) then ! Only one constituent name C ncon=ncon+1 acon(ncon)=conbuf(1:6) call changecase(acon(ncon),acon(ncon),len,0) !Change to lower case C else !Assume more than one constituent name, so related constituents C end_of_cons=.true. C endif C endif C else if(ios.gt.0) then ! Error C call error_handler(1,020) ! Error point 020 C else ! Assume end of file C end_of_file=.true. C endif C end do C write(6,17) ncon 17 format('No. of constituents = ',i4) C nequ=2*ncon+1 nequm1=nequ-1 C do icon=1,ncon ! Set index array for relating constutuents index_ref(icon)=0 end do C if(.not.end_of_file) then ! There are some related constituents to read C nref=0 ios=0 C do while(ios.eq.0) C if(nref.ne.0) read(11,10,iostat=ios) conbuf !Have already read 1st one C if(ios.eq.0) then ! Read is O.K. C if(conbuf(1:1).ne.'#') then ! Ignore comments C asplit=conbuf(1:6) C ilhs=lhschar(conbuf(7:120)) aref=conbuf(ilhs+6:ilhs+11) C len=lenchar(asplit) ! Change to lower case call changecase(asplit,asplit,len,0) len=lenchar(aref) ! Change to lower case call changecase(aref,aref,len,0) C do icon=1,ncon if(acon(icon).eq.asplit) isp=icon if(acon(icon).eq.aref) iref=icon end do C index_ref(isp)=iref nref=nref+1 aisp(nref)=isp C endif C else if(ios.gt.0) then ! Error C call error_handler(1,021) ! Error point 021 C endif C end do C else ! There are no related constituents C nref=0 C endif C write(6,16) nref 16 format('No. of related constituents = ',i4) C l_tz_g=.false. C if(nref.ne.0) then C C Input reference data from port C call areadfil3(7,'PORT_REF',8,port_ref,ifail) call error_handler(ifail,022) ! Error point 022 C len=lenchar(port_ref) call changecase(port_ref,port_ref,len,0) ! Convert to lower case C open(unit=7,file=tidal_dir(1:lentid)//port_ref,status='old', $ iostat=ios) call error_handler(ios,023) ! Error point 023 C C Note that constituent list may be in any order .... C ndum=0 ios=0 do while(ios.eq.0) read(7,6,iostat=ios) refbuf 6 format(a40) if(ios.gt.0) call error_handler(ios,024) ! Error point 024 if(ios.eq.0) then if(refbuf(1:1).eq.'#') then ! A comment if(refbuf(1:5).eq.'# TZ ') then ! Found the time zone read(refbuf(6:40),*,iostat=ios) tz_g call error_handler(ios,025) ! Error point 025 l_tz_g=.true. endif else ! Harmonic constants if(refbuf(1:3).ne.'z0 '.and.refbuf(1:3).ne.'Z0 ') then ! Not Z0 len=index(refbuf,' ')-1 ! Find separator between name and (h,g) if(len.gt.6) len=6 ! Assume format is a6 and there is no space ndum=ndum+1 arel(ndum)=' ' ! For safety write(formt,'("(a",i3,")")') len read(refbuf,formt,iostat=ios) arel(ndum)(1:len) call error_handler(ios,026) ! Error point 026 C do ichar=len+1,40 ! Change any "*" to space if(refbuf(ichar:ichar).eq.'*') $ refbuf(ichar:ichar)=' ' end do C read(refbuf(len+1:40),*,iostat=ios) h(ndum,1),g(ndum,1) call error_handler(ios,027) ! Error point 027 len=lenchar(arel(ndum)) call changecase(arel(ndum),arel(ndum),len,0) ! Change to L.C endif endif endif end do C C Set up array for related constituents: C do icon=1,pcon hrel(icon)=0.d0 grel(icon)=0.d0 end do C do iref=1,nref C aispi=aisp(iref) asplit=acon(aispi) indexi=index_ref(aispi) aref=acon(indexi) logs=.false. logr=.false. icon=1 C do while((.not.logs.or..not.logr) $ .and.icon.le.ndum) C atest=arel(icon) C if(asplit.eq.atest) then hrel(aispi)=h(icon,1) ! Found split const. - convert to rad. grel(aispi)=g(icon,1) logs=.true. endif C if(aref.eq.atest) then hrel(indexi)=h(icon,1) ! Found ref. const. - convert to rad. grel(indexi)=g(icon,1) logr=.true. endif C icon=icon+1 C end do C if(.not.logs.or..not.logr) $ call error_handler(ios,029) ! Error point 029 C end do C endif C C Sort out time zone: C if(.not.l_tz_t.and..not.l_tz_g) then write(6,11) 11 format('No time zone has been defined') call error_handler(ios,028) ! Error point 028 else if(.not.l_tz_t.and.l_tz_g) then ! Assume time zone of related cons. tz_out=tz_g write(6,12) tz_g 12 format('Time zone of input and output (',f6.2, $ ') assumed to be that of related constituents') else if(l_tz_t.and..not.l_tz_g) then ! Assume time zone in param. file tz_out=tz_t write(6,13) tz_t 13 format('Time zone of input and output (',f6.2, $ ') assumed to be that of parameter file') else tz_out=tz_t if(tz_t.ne.tz_g) then write(6,14) 14 format('Different time zone between related constituents ', $ 'and parameter file.') endif write(6,15) tz_t 15 format('Time zone of output (',f6.2, $ ') selected to be that of parameter file') endif C C Read number of data columns and starting year: C read(13,*) nset,iye iye_start=iye ! Remember year at start iye_last=iye call julday(iye,1,1,0,0,0.,jul_start,ifail) call error_handler(ifail,030) ! Error point 030 jul_year=jul_start C C Clear covariance matrices: C do iset=1,pset var(iset)=0.d0 end do C do iequ1=1,pequ C do iset=1,pset bcov(iequ1,iset)=0.d0 b(iequ1,iset)=0.d0 end do C do iequ2=1,pequ acov(iequ1,iequ2)=0.d0 a(iequ1,iequ2)=0.d0 end do C end do C x(1)=1.d0 C ndat=0 ios=0 C do while(ios.eq.0) C read(13,4,iostat=ios) datbuf 4 format(a120) C if(ios.gt.0) call error_handler(ios,031) ! Error point 031 C if(ios.eq.0) then C read(datbuf,*,iostat=ios) t,(dat(iset),iset=1,nset) C ndat=ndat+1 if((ndat/500)*500.eq.ndat) then write(6,29) ndat 29 format('Processing data point ',i5) endif C if(ltrep) then C if(ndat.eq.1) then lyearchange=.false. else lyearchange=(t.lt.tlast) endif C if(lyearchange) then iye=iye+1 ! Update year endif C else C C Find year: C call caldat(iye,mon,idy,ihr,min,sec, $ dble(t)/86400.d0+jul_start, $ ifail) call error_handler(ifail,032) ! Error point 032 C lyearchange=(iye.ne.iye_last) C endif C if(ndat.eq.1.or.lyearchange) then ! Load astronomical data C C Removed 7/9/2001: C C if(iye.lt.1970.or.iye.gt.2030) C $ call error_handler(1,033) ! Error point 033 C iye_last=iye C if(.not.ltrep) then call julday(iye,1,1,0,0,0.,jul_year, $ ifail) call error_handler(ifail,034) ! Error point 034 C endif C call ldast(acon,iye,ncon,sigma,j1,j2,v1,v2,v0,year, $ tidal_dir,lentid, $ ifail) call error_handler(ifail,035) ! Error point 035 C endif C C Change time so that it is relative to start of current year, if necessary: C if(.not.ltrep) t=t-sngl((jul_year-jul_start)*86400.d0) C C Interpolate j and vpv: C w2=(t-tz_out*3600.d0)/year ! Time in GMT C w2=t/year ! Local time w1=1.d0-w2 do icon=1,ncon j(icon)=j1(icon)*w1+j2(icon)*w2 vpv(icon)=v0(icon)+v1(icon)*w1+v2(icon)*w2 end do C do icon=1,ncon ic=2*icon is=ic+1 arg=vpv(icon)+sigma(icon)*t x(ic)=j(icon)*dcos(arg*d2r) x(is)=j(icon)*dsin(arg*d2r) end do C do iset=1,nset var(iset)=var(iset)+dat(iset)**2 end do C do iequ1=1,nequ C do iequ2=iequ1,nequ acov(iequ1,iequ2)=acov(iequ1,iequ2)+x(iequ1)*x(iequ2) end do C do iset=1,nset bcov(iequ1,iset)=bcov(iequ1,iset)+x(iequ1)*dat(iset) end do end do C C Output X for temp. storage: C if(lres) write(15)(x(iequ1),iequ1=1,nequ) C tlast=t C endif C end do C write(6,5) ndat 5 format(/'No. of data points = ',i5/) C do iequ1=1,nequm1 ! Generate other half of ACOV do iequ2=iequ1+1,nequ acov(iequ2,iequ1)=acov(iequ1,iequ2) end do end do C C Divide VAR, ACOV and BCOV by number of data points: C do iset=1,nset var(iset)=var(iset)/dble(ndat) end do C do iequ1=1,nequ C do iequ2=1,nequ acov(iequ1,iequ2)=acov(iequ1,iequ2)/dble(ndat) end do C do iset=1,nset bcov(iequ1,iset)=bcov(iequ1,iset)/dble(ndat) end do C end do C C Generate A and B: C do iequ1=1,nequ C logs=.false. if(iequ1.ne.1.and.nref.ne.0) then icon=iequ1/2 ! Test for split constituent if(index_ref(icon).gt.0) logs=.true. endif C do iequ2=1,nequ if(.not.logs) a(iequ1,iequ2)=acov(iequ1,iequ2) end do C do iset=1,nset if(.not.logs) b(iequ1,iset)=bcov(iequ1,iset) end do C end do C if(nref.ne.0.) then C C Deal with split constutuents: C do icon=1,ncon C if(index_ref(icon).ne.0) then C r=hrel(icon)/hrel(index_ref(icon)) delta=grel(icon)-grel(index_ref(icon)) C C Correct reference phases to final time zone: C if(l_tz_g) then del_tz=(tz_out-tz_g)*3600.d0 ! Change in time zone (secs.) delta=delta+(sigma(icon)-sigma(index_ref(icon)))*del_tz endif C ics=2*icon iss=ics+1 icr=2*index_ref(icon) isr=icr+1 C rc=r*cos(delta*d2r) rs=r*sin(delta*d2r) C a(ics,ics)=1.d0 a(ics,icr)=-rc a(ics,isr)=rs a(iss,iss)=1.d0 a(iss,icr)=-rs a(iss,isr)=-rc C endif C end do C endif C C Solve A*C=B, using SVD: C job=11 C call dsvdc(a,pequ,nequ,nequ, $ wsvd,e,usvd,pequ,vsvd,pequ, $ tmpsvd,job,info) C write(6,7) info 7 format('Info from LINPACK dsvdc: ',i5) C call error_handler(info,036) ! Error point 036 C write(6,3) (wsvd(iequ1),iequ1=1,nequ) 3 format('Singular values : '/ $ (5x,5e12.4)) C C Back substitute: C do iset=1,nset C do iequ1=1,nequ b1(iequ1)=b(iequ1,iset) end do C call dsvbksb2(usvd,wsvd,vsvd,nequ,nequ,nequ, $ pequ,pequ,b1,c1,tmpsvd) C do iequ1=1,nequ c(iequ1,iset)=c1(iequ1) end do C end do C do iset=1,nset C rmean(iset)=c(1,iset) C do icon=1,ncon C ic=2*icon is=ic+1 h(icon,iset)=sqrt(c(ic,iset)**2+c(is,iset)**2) C if(c(is,iset).ne.0.d0.or.c(ic,iset).ne.0.d0) then arg=atan2(c(is,iset),c(ic,iset))/d2r if(arg.lt.0.d0) arg=arg+360.d0 else arg=0.d0 endif C g(icon,iset)=arg C end do C end do C C Output results: C do ihead=1,nhead len=lenchar(head(ihead)) write(formt,'("(a2,a",i3,")")') len write(8,formt) '# ',head(ihead)(1:len) end do write(8,30) '# TZ ',tz_out, $ '# ',((' H G '),iset=1,nset) C 30 format(a6,f6.2/ C $ a9,(a14)) 30 format(a6,f6.2/ $ a9,100(a14)) C write(8,19)(rmean(iset),iset=1,nset) C 19 format('Z0 ',(f7.4,7x)) 19 format('Z0 ',100(f7.4,7x)) C do icon=1,ncon astar=' ' if(index_ref(icon).gt.0) astar='*' len=lenchar(acon(icon)) ! Convert to upper case call changecase(acon(icon),acon(icon),len,1) write(8,20) acon(icon),astar,(h(icon,iset),g(icon,iset), $ iset=1,nset) C 20 format(a6,1x,a1,(1x,f5.3,2x,f5.1,1x)) 20 format(a6,1x,a1,100(1x,f5.3,2x,f5.1,1x)) end do C C Compute residual sum of squares; C do iset=1,nset C resiset=var(iset) C do iequ1=1,nequ resiset=resiset-2.*c(iequ1,iset)*bcov(iequ1,iset) do iequ2=1,nequ resiset=resiset $ +c(iequ1,iset)*c(iequ2,iset)*acov(iequ1,iequ2) end do end do C res(iset)=resiset C end do C C Compute statistics: C do iset=1,nset C var(iset)=var(iset)-bcov(1,iset)**2 svar(iset)=sqrt(var(iset)) sres(iset)=sqrt(abs(res(iset))) r1(iset)=res(iset)/var(iset) r2(iset)=sres(iset)/svar(iset) C end do C write(6,21)(iset,iset=1,nset) C 21 format(//'STATISTICS'//'DATA SET NO :',19x,(i3,9x)/) 21 format(//'STATISTICS'//'DATA SET NO :',19x,100(i3,9x)/) C write(6,22)(var(iset),iset=1,nset) C 22 format(/'VARIANCE ORIGINAL (VO) ',(f12.3)) 22 format(/'VARIANCE ORIGINAL (VO) ',100(f12.3)) C write(6,23)(res(iset),iset=1,nset) C 23 format('VARIANCE RESIDUAL (VR) ',(f12.5)) 23 format('VARIANCE RESIDUAL (VR) ',100(f12.5)) C write(6,24)(r1(iset),iset=1,nset) C 24 format('RATIO VR/VO',14x,(f12.4)) 24 format('RATIO VR/VO',14x,100(f12.4)) C write(6,25)(svar(iset),iset=1,nset) C 25 format('SQRT(VO)',17x,(f12.3)) 25 format('SQRT(VO)',17x,100(f12.3)) C write(6,26)(sres(iset),iset=1,nset) C 26 format('SQRT(VR)',17x,(f12.5)) 26 format('SQRT(VR)',17x,100(f12.5)) C write(6,27)(r2(iset),iset=1,nset) C 27 format('RATIO SQRT(VR)/SQRT(VO) ',(f12.4)) 27 format('RATIO SQRT(VR)/SQRT(VO) ',100(f12.4)) C write(6,28) 28 format(//) C C Compute residuals: C if(lres) then C rewind 13 rewind 15 read(13,*) write(10,32) nset,iye_start 32 format(1x,2i5) C do idat=1,ndat C read(13,*) t,(dat(iset),iset=1,nset) read(15) (x(iequ1),iequ1=1,nequ) C do iset=1,nset do iequ1=1,nequ dat(iset)=dat(iset)-c(iequ1,iset)*x(iequ1) end do end do C write(10,31) t,(dat(iset),iset=1,nset) C 31 format(f10.0,(f10.3)) 31 format(f10.0,100(f10.3)) C end do C endif C stop C end C subroutine areadfil3(nin,anot,nchar,avar,ifail) C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : UTILITY REF:JRH:08:09:1999 * C * * C * REVISION : ------------- JRH:--:--:1999 * C * * C * SOURCE : tda1537.f * C * ROUTINE NAME: areadfil3 * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Reads CHARACTER*120 variable from file based on keyword. * C * (Modified version of areadfil, without output to operator * C * or log file) * C * * C * nin ..... File input device * C * anot .... Keyword in file * C * nchar ... Number of characters in ANOT (max. 40) * C * avar .... Resultant CHARACTER*120 variable * C * ifail ... 0 for successful execution, otherwise 1 * C * * C * (A modified version of areadfil2 for 120 characters.) * C * * C ****************************************************************************** C implicit none C integer*4 nin,nchar,ifail character*(*) anot character*120 avar C integer*4 ios logical found character*161 buff C rewind(nin) ios=0 found=.false. do while(ios.eq.0.and..not.found) read(nin,1,iostat=ios) buff 1 format(a160) if(buff(1:nchar).eq.anot.and. $ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found avar(1:120)=buff(nchar+2:nchar+121) found=.true. endif end do if(ios.eq.0) then ifail=0 ! Match found and no error else ifail=1 ! No match found endif end C subroutine ldast(acon,iye,ncon,sigma,j1,j2,v1,v2,v0,year, $ tidal_dir,lentid, $ ifail) C ****************************************************************************** C * * C * VAX/SUN SOURCE CODE * C * * C * PROGRAM SET : TIDES REF:JRH:10:09:1999 * C * * C * REVISION : Mod. to allow any year and both "98" * C * and "1998" type formats JRH:07:09:2001 * C * * C * SOURCE : tda1537.f * C * ROUTINE NAME: ldast * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Loads astronomical data from disc (converts to radians * C * and radians/sec). * C * * C * NOTE this uses unit 9 for input. * C * * C ****************************************************************************** C implicit none C integer pcon,ptry C parameter(pcon=114, ! >/= Max. no. of constituents used $ ptry=3) ! No. of tries at different const. names C real*8 sigma(pcon) real*8 j1(pcon),j2(pcon),v1(pcon),v2(pcon),v0(pcon) real*8 year C integer*4 iye,ncon,lentid,ifail C character*120 tidal_dir character*80 formt character*6 acon(pcon) C real*8 jul_this,jul_next real*8 d2r real*8 vpv real*8 dum C integer*4 icon,itry,nyear,iyear integer*4 len,lenchar integer*4 idum integer*4 ios C logical first logical match,found C character*6 atry(ptry,2) C character*40 buf character*6 filnam C data first/.true./ 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 save first C if(first) then C d2r=datan(1.d0)*4.d0/180.d0 ! Degrees to radians C endif C call julday(iye,1,1,0,0,0.,jul_this,ifail) C if(ifail.ne.0) then ifail=1 return endif C call julday(iye+1,1,1,0,0,0.,jul_next,ifail) C if(ifail.ne.0) then ifail=2 return endif C year=(jul_next-jul_this)*86400.d0 ! Length of current year in seconds C write(6,4) 4 format(/) C do icon=1,ncon C ios=0 itry=0 found=.false. do while(.not.found) len=lenchar(acon(icon)) filnam(1:len)=acon(icon)(1:len) write(formt,'("(a27,a",i3,")")') len write(6,formt) 'Opening ast. arg. file for ',filnam(1:len)! Diagnostic C open(unit=9, $ file=tidal_dir(1:lentid)//filnam(1:len)//'.con', $ status='old',iostat=ios) C found=(ios.eq.0) ! .true. if open O.K. match=(ios.eq.0) ! If open not O.K., try other names C do while(.not.match.and.itry.lt.ptry) itry=itry+1 if(acon(icon).eq.atry(itry,1)) then acon(icon)=atry(itry,2) match=.true. else if(acon(icon).eq.atry(itry,2)) then acon(icon)=atry(itry,1) match=.true. endif end do C if(.not.found.and..not.match) then ifail=3 return endif end do C read(9,*,iostat=ios) dum C if(ios.ne.0) then ifail=4 return endif C if(first) sigma(icon)=dum/3600.d0 ! Convert to degrees per second C nyear=iye-1900 ! Note this format only works around 2000 if(nyear.ge.100) nyear=nyear-100 iyear=9999 do while(iyear.ne.nyear.and.iyear.ne.iye) ! &7/9/2001 C read(9,1,iostat=ios) buf 1 format(a40) C if(ios.ne.0) then ifail=5 return endif C read(buf,*,iostat=ios) iyear if(ios.ne.0) then ifail=6 return endif C if(iyear.eq.nyear.or.iyear.eq.iye) then ! &7/9/2001 C read(buf,*,iostat=ios) idum,j1(icon),v1(icon),vpv if(ios.ne.0) then ifail=7 return endif C v0(icon)=vpv-v1(icon) C read(9,*,iostat=ios) idum,j2(icon),v2(icon) if(ios.ne.0) then ifail=8 return endif C endif C end do C close(unit=9) C end do C first=.false. C return C end C C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : MODELLING REF:JRH:10:01:1994 * C * * C * REVISION : All "1" array declarations changed to * C * "*" (so that subscript check works) JRH:10:01:1994 * C * * C * SOURCE : tad1537.f * C * ROUTINE NAME: -- * C * TYPE : SUBROUTINES * C * * C * FUNCTION : Following are : * C * * C * LINPACK routines dsvdc and drot * C * BLAS routines daxpy,ddot,dscal,dswap,dnrm2,drotg * C * * C ****************************************************************************** C subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info) integer ldx,n,p,ldu,ldv,job,info double precision x(ldx,*),s(*),e(*),u(ldu,*),v(ldv,*),work(*) c c c dsvdc is a subroutine to reduce a double precision nxp matrix x c by orthogonal transformations u and v to diagonal form. the c diagonal elements s(i) are the singular values of x. the c columns of u are the corresponding left singular vectors, c and the columns of v the right singular vectors. c c on entry c c x double precision(ldx,p), where ldx.ge.n. c x contains the matrix whose singular value c decomposition is to be computed. x is c destroyed by dsvdc. c c ldx integer. c ldx is the leading dimension of the array x. c c n integer. c n is the number of rows of the matrix x. c c p integer. c p is the number of columns of the matrix x. c c ldu integer. c ldu is the leading dimension of the array u. c (see below). c c ldv integer. c ldv is the leading dimension of the array v. c (see below). c c work double precision(n). c work is a scratch array. c c job integer. c job controls the computation of the singular c vectors. it has the decimal expansion ab c with the following meaning c c a.eq.0 do not compute the left singular c vectors. c a.eq.1 return the n left singular vectors c in u. c a.ge.2 return the first min(n,p) singular c vectors in u. c b.eq.0 do not compute the right singular c vectors. c b.eq.1 return the right singular vectors c in v. c c on return c c s double precision(mm), where mm=min(n+1,p). c the first min(n,p) entries of s contain the c singular values of x arranged in descending c order of magnitude. c c e double precision(p), c e ordinarily contains zeros. however see the c discussion of info for exceptions. c c u double precision(ldu,k), where ldu.ge.n. if c joba.eq.1 then k.eq.n, if joba.ge.2 c then k.eq.min(n,p). c u contains the matrix of left singular vectors. c u is not referenced if joba.eq.0. if n.le.p c or if joba.eq.2, then u may be identified with x c in the subroutine call. c c v double precision(ldv,p), where ldv.ge.p. c v contains the matrix of right singular vectors. c v is not referenced if job.eq.0. if p.le.n, c then v may be identified with x in the c subroutine call. c c info integer. c the singular values (and their corresponding c singular vectors) s(info+1),s(info+2),...,s(m) c are correct (here m=min(n,p)). thus if c info.eq.0, all the singular values and their c vectors are correct. in any event, the matrix c b = trans(u)*x*v is the bidiagonal matrix c with the elements of s on its diagonal and the c elements of e on its super-diagonal (trans(u) c is the transpose of u). thus the singular c values of x and b are the same. c c linpack. this version dated 08/14/78 . c correction made to shift 2/84. c g.w. stewart, university of maryland, argonne national lab. c c dsvdc uses the following functions and subprograms. c c external drot c blas daxpy,ddot,dscal,dswap,dnrm2,drotg c fortran dabs,dmax1,max0,min0,mod,dsqrt c c internal variables c integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit, * mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1 double precision ddot,t double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn, * smm1,t1,test,ztest logical wantu,wantv c c c set the maximum number of iterations. c maxit = 30 c c determine what is to be computed. c wantu = .false. wantv = .false. jobu = mod(job,100)/10 ncu = n if (jobu .gt. 1) ncu = min0(n,p) if (jobu .ne. 0) wantu = .true. if (mod(job,10) .ne. 0) wantv = .true. c c reduce x to bidiagonal form, storing the diagonal elements c in s and the super-diagonal elements in e. c info = 0 nct = min0(n-1,p) nrt = max0(0,min0(p-2,n)) lu = max0(nct,nrt) if (lu .lt. 1) go to 170 do 160 l = 1, lu lp1 = l + 1 if (l .gt. nct) go to 20 c c compute the transformation for the l-th column and c place the l-th diagonal in s(l). c s(l) = dnrm2(n-l+1,x(l,l),1) if (s(l) .eq. 0.0d0) go to 10 if (x(l,l) .ne. 0.0d0) s(l) = dsign(s(l),x(l,l)) call dscal(n-l+1,1.0d0/s(l),x(l,l),1) x(l,l) = 1.0d0 + x(l,l) 10 continue s(l) = -s(l) 20 continue if (p .lt. lp1) go to 50 do 40 j = lp1, p if (l .gt. nct) go to 30 if (s(l) .eq. 0.0d0) go to 30 c c apply the transformation. c t = -ddot(n-l+1,x(l,l),1,x(l,j),1)/x(l,l) call daxpy(n-l+1,t,x(l,l),1,x(l,j),1) 30 continue c c place the l-th row of x into e for the c subsequent calculation of the row transformation. c e(j) = x(l,j) 40 continue 50 continue if (.not.wantu .or. l .gt. nct) go to 70 c c place the transformation in u for subsequent back c multiplication. c do 60 i = l, n u(i,l) = x(i,l) 60 continue 70 continue if (l .gt. nrt) go to 150 c c compute the l-th row transformation and place the c l-th super-diagonal in e(l). c e(l) = dnrm2(p-l,e(lp1),1) if (e(l) .eq. 0.0d0) go to 80 if (e(lp1) .ne. 0.0d0) e(l) = dsign(e(l),e(lp1)) call dscal(p-l,1.0d0/e(l),e(lp1),1) e(lp1) = 1.0d0 + e(lp1) 80 continue e(l) = -e(l) if (lp1 .gt. n .or. e(l) .eq. 0.0d0) go to 120 c c apply the transformation. c do 90 i = lp1, n work(i) = 0.0d0 90 continue do 100 j = lp1, p call daxpy(n-l,e(j),x(lp1,j),1,work(lp1),1) 100 continue do 110 j = lp1, p call daxpy(n-l,-e(j)/e(lp1),work(lp1),1,x(lp1,j),1) 110 continue 120 continue if (.not.wantv) go to 140 c c place the transformation in v for subsequent c back multiplication. c do 130 i = lp1, p v(i,l) = e(i) 130 continue 140 continue 150 continue 160 continue 170 continue c c set up the final bidiagonal matrix or order m. c m = min0(p,n+1) nctp1 = nct + 1 nrtp1 = nrt + 1 if (nct .lt. p) s(nctp1) = x(nctp1,nctp1) if (n .lt. m) s(m) = 0.0d0 if (nrtp1 .lt. m) e(nrtp1) = x(nrtp1,m) e(m) = 0.0d0 c c if required, generate u. c if (.not.wantu) go to 300 if (ncu .lt. nctp1) go to 200 do 190 j = nctp1, ncu do 180 i = 1, n u(i,j) = 0.0d0 180 continue u(j,j) = 1.0d0 190 continue 200 continue if (nct .lt. 1) go to 290 do 280 ll = 1, nct l = nct - ll + 1 if (s(l) .eq. 0.0d0) go to 250 lp1 = l + 1 if (ncu .lt. lp1) go to 220 do 210 j = lp1, ncu t = -ddot(n-l+1,u(l,l),1,u(l,j),1)/u(l,l) call daxpy(n-l+1,t,u(l,l),1,u(l,j),1) 210 continue 220 continue call dscal(n-l+1,-1.0d0,u(l,l),1) u(l,l) = 1.0d0 + u(l,l) lm1 = l - 1 if (lm1 .lt. 1) go to 240 do 230 i = 1, lm1 u(i,l) = 0.0d0 230 continue 240 continue go to 270 250 continue do 260 i = 1, n u(i,l) = 0.0d0 260 continue u(l,l) = 1.0d0 270 continue 280 continue 290 continue 300 continue c c if it is required, generate v. c if (.not.wantv) go to 350 do 340 ll = 1, p l = p - ll + 1 lp1 = l + 1 if (l .gt. nrt) go to 320 if (e(l) .eq. 0.0d0) go to 320 do 310 j = lp1, p t = -ddot(p-l,v(lp1,l),1,v(lp1,j),1)/v(lp1,l) call daxpy(p-l,t,v(lp1,l),1,v(lp1,j),1) 310 continue 320 continue do 330 i = 1, p v(i,l) = 0.0d0 330 continue v(l,l) = 1.0d0 340 continue 350 continue c c main iteration loop for the singular values. c mm = m iter = 0 360 continue c c quit if all the singular values have been found. c c ...exit if (m .eq. 0) go to 620 c c if too many iterations have been performed, set c flag and return. c if (iter .lt. maxit) go to 370 info = m c ......exit go to 620 370 continue c c this section of the program inspects for c negligible elements in the s and e arrays. on c completion the variables kase and l are set as follows. c c kase = 1 if s(m) and e(l-1) are negligible and l.lt.m c kase = 2 if s(l) is negligible and l.lt.m c kase = 3 if e(l-1) is negligible, l.lt.m, and c s(l), ..., s(m) are not negligible (qr step). c kase = 4 if e(m-1) is negligible (convergence). c do 390 ll = 1, m l = m - ll c ...exit if (l .eq. 0) go to 400 test = dabs(s(l)) + dabs(s(l+1)) ztest = test + dabs(e(l)) if (ztest .ne. test) go to 380 e(l) = 0.0d0 c ......exit go to 400 380 continue 390 continue 400 continue if (l .ne. m - 1) go to 410 kase = 4 go to 480 410 continue lp1 = l + 1 mp1 = m + 1 do 430 lls = lp1, mp1 ls = m - lls + lp1 c ...exit if (ls .eq. l) go to 440 test = 0.0d0 if (ls .ne. m) test = test + dabs(e(ls)) if (ls .ne. l + 1) test = test + dabs(e(ls-1)) ztest = test + dabs(s(ls)) if (ztest .ne. test) go to 420 s(ls) = 0.0d0 c ......exit go to 440 420 continue 430 continue 440 continue if (ls .ne. l) go to 450 kase = 3 go to 470 450 continue if (ls .ne. m) go to 460 kase = 1 go to 470 460 continue kase = 2 l = ls 470 continue 480 continue l = l + 1 c c perform the task indicated by kase. c go to (490,520,540,570), kase c c deflate negligible s(m). c 490 continue mm1 = m - 1 f = e(m-1) e(m-1) = 0.0d0 do 510 kk = l, mm1 k = mm1 - kk + l t1 = s(k) call drotg(t1,f,cs,sn) s(k) = t1 if (k .eq. l) go to 500 f = -sn*e(k-1) e(k-1) = cs*e(k-1) 500 continue if (wantv) call drot(p,v(1,k),1,v(1,m),1,cs,sn) 510 continue go to 610 c c split at negligible s(l). c 520 continue f = e(l-1) e(l-1) = 0.0d0 do 530 k = l, m t1 = s(k) call drotg(t1,f,cs,sn) s(k) = t1 f = -sn*e(k) e(k) = cs*e(k) if (wantu) call drot(n,u(1,k),1,u(1,l-1),1,cs,sn) 530 continue go to 610 c c perform one qr step. c 540 continue c c calculate the shift. c scale = dmax1(dabs(s(m)),dabs(s(m-1)),dabs(e(m-1)), * dabs(s(l)),dabs(e(l))) sm = s(m)/scale smm1 = s(m-1)/scale emm1 = e(m-1)/scale sl = s(l)/scale el = e(l)/scale b = ((smm1 + sm)*(smm1 - sm) + emm1**2)/2.0d0 c = (sm*emm1)**2 shift = 0.0d0 if (b .eq. 0.0d0 .and. c .eq. 0.0d0) go to 550 shift = dsqrt(b**2+c) if (b .lt. 0.0d0) shift = -shift shift = c/(b + shift) 550 continue f = (sl + sm)*(sl - sm) + shift g = sl*el c c chase zeros. c mm1 = m - 1 do 560 k = l, mm1 call drotg(f,g,cs,sn) if (k .ne. l) e(k-1) = f f = cs*s(k) + sn*e(k) e(k) = cs*e(k) - sn*s(k) g = sn*s(k+1) s(k+1) = cs*s(k+1) if (wantv) call drot(p,v(1,k),1,v(1,k+1),1,cs,sn) call drotg(f,g,cs,sn) s(k) = f f = cs*e(k) + sn*s(k+1) s(k+1) = -sn*e(k) + cs*s(k+1) g = sn*e(k+1) e(k+1) = cs*e(k+1) if (wantu .and. k .lt. n) * call drot(n,u(1,k),1,u(1,k+1),1,cs,sn) 560 continue e(m-1) = f iter = iter + 1 go to 610 c c convergence. c 570 continue c c make the singular value positive. c if (s(l) .ge. 0.0d0) go to 580 s(l) = -s(l) if (wantv) call dscal(p,-1.0d0,v(1,l),1) 580 continue c c order the singular value. c 590 if (l .eq. mm) go to 600 c ...exit if (s(l) .ge. s(l+1)) go to 600 t = s(l) s(l) = s(l+1) s(l+1) = t if (wantv .and. l .lt. p) * call dswap(p,v(1,l),1,v(1,l+1),1) if (wantu .and. l .lt. n) * call dswap(n,u(1,l),1,u(1,l+1),1) l = l + 1 go to 590 600 continue iter = 0 m = m - 1 610 continue go to 360 620 continue return end C subroutine drot (n,dx,incx,dy,incy,c,s) c c applies a plane rotation. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp,c,s integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = c*dx(ix) + s*dy(iy) dy(iy) = c*dy(iy) - s*dx(ix) dx(ix) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n dtemp = c*dx(i) + s*dy(i) dy(i) = c*dy(i) - s*dx(i) dx(i) = dtemp 30 continue return end C subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end C double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end C subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end C subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end C double precision function dnrm2 ( n, dx, incx) integer i, incx, ix, j, n, next double precision dx(*), cutlo, cuthi, hitest, sum, xmax,zero,one data zero, one /0.0d0, 1.0d0/ c c euclidean norm of the n-vector stored in dx() with storage c increment incx . c if n .le. 0 return with result = 0. c if n .ge. 1 then incx must be .ge. 1 c c c.l.lawson, 1978 jan 08 c modified to correct failure to update ix, 1/25/92. c modified 3/93 to return if incx .le. 0. c c four phase method using two built-in constants that are c hopefully applicable to all machines. c cutlo = maximum of dsqrt(u/eps) over all known machines. c cuthi = minimum of dsqrt(v) over all known machines. c where c eps = smallest no. such that eps + 1. .gt. 1. c u = smallest positive no. (underflow limit) c v = largest no. (overflow limit) c c brief outline of algorithm.. c c phase 1 scans zero components. c move to phase 2 when a component is nonzero and .le. cutlo c move to phase 3 when a component is .gt. cutlo c move to phase 4 when a component is .ge. cuthi/m c where m = n for x() real and m = 2*n for complex. c c values for cutlo and cuthi.. c from the environmental parameters listed in the imsl converter c document the limiting values are as follows.. c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are c univac and dec at 2**(-103) c thus cutlo = 2**(-51) = 4.44089e-16 c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. c thus cuthi = 2**(63.5) = 1.30438e19 c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. c thus cutlo = 2**(-33.5) = 8.23181d-11 c cuthi, d.p. same as s.p. cuthi = 1.30438d19 c data cutlo, cuthi / 8.232d-11, 1.304d19 / c data cutlo, cuthi / 4.441e-16, 1.304e19 / data cutlo, cuthi / 8.232d-11, 1.304d19 / c if(n .gt. 0 .and. incx.gt.0) go to 10 dnrm2 = zero go to 300 c 10 assign 30 to next sum = zero i = 1 ix = 1 c begin main loop 20 go to next,(30, 50, 70, 110) 30 if( dabs(dx(i)) .gt. cutlo) go to 85 assign 50 to next xmax = zero c c phase 1. sum is zero c 50 if( dx(i) .eq. zero) go to 200 if( dabs(dx(i)) .gt. cutlo) go to 85 c c prepare for phase 2. assign 70 to next go to 105 c c prepare for phase 4. c 100 continue ix = j assign 110 to next sum = (sum / dx(i)) / dx(i) 105 xmax = dabs(dx(i)) go to 115 c c phase 2. sum is small. c scale to avoid destructive underflow. c 70 if( dabs(dx(i)) .gt. cutlo ) go to 75 c c common code for phases 2 and 4. c in phase 4 sum is large. scale to avoid overflow. c 110 if( dabs(dx(i)) .le. xmax ) go to 115 sum = one + sum * (xmax / dx(i))**2 xmax = dabs(dx(i)) go to 200 c 115 sum = sum + (dx(i)/xmax)**2 go to 200 c c c prepare for phase 3. c 75 sum = (sum * xmax) * xmax c c c for real or d.p. set hitest = cuthi/n c for complex set hitest = cuthi/(2*n) c 85 hitest = cuthi/float( n ) c c phase 3. sum is mid-range. no scaling. c do 95 j = ix,n if(dabs(dx(i)) .ge. hitest) go to 100 sum = sum + dx(i)**2 i = i + incx 95 continue dnrm2 = dsqrt( sum ) go to 300 c 200 continue ix = ix + 1 i = i + incx if( ix .le. n ) go to 20 c c end of main loop. c c compute square root and adjust for scaling. c dnrm2 = xmax * dsqrt(sum) 300 continue return end C subroutine drotg(da,db,c,s) c c construct givens plane rotation. c jack dongarra, linpack, 3/11/78. c modified 9/27/86. c double precision da,db,c,s,roe,scale,r,z c roe = db if( dabs(da) .gt. dabs(db) ) roe = da scale = dabs(da) + dabs(db) if( scale .ne. 0.0d0 ) go to 10 c = 1.0d0 s = 0.0d0 r = 0.0d0 go to 20 10 r = scale*dsqrt((da/scale)**2 + (db/scale)**2) r = dsign(1.0d0,roe)*r c = da/r s = db/r 20 z = s if( dabs(c) .gt. 0.0d0 .and. dabs(c) .le. s ) z = 1.0d0/c da = r db = z return end C subroutine dsvbksb2(u,w,v,n,p,minnp,nmax,pmax,b,x,tmp) C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : MODELLING REF:JRH:11:01:1994 * C * * C * REVISION : Re-ordered declarations JRH:17:05:2000 * C * * C * SOURCE : tda1537 * C * ROUTINE NAME: DSVBKSB2 * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Singular value decomposition solution * C * (from Numerical Recipes - D.P. version of SVBKSB, * C * modified for use with LINPACK) * C * * C * U ..... Input N by MINNP matrix * C * W ..... Input MINNP element weight vector * C * V ..... Input P by MINNP matrix * C * ( NOT the transpose of "V") * C * N ..... Number of rows of U * C * P ..... Number of rows of V * C * MINNP . MIN(N,P) * C * NMAX .. Physical number of rows of U * C * PMAX .. Physical number of rows of V * C * B ..... Input N element vector (right hand side) * C * X ..... Output P element solution vector * C * TMP ... Workspace MINNP element vector * C * * C ****************************************************************************** C implicit none C integer*4 n,p,minnp,nmax,pmax real*8 u(nmax,*),w(*),v(pmax,*),b(*),x(*),tmp(*) C real*8 s integer*4 i,j,jj do j=1,minnp s=0. if(w(j).ne.0.)then do i=1,n s=s+u(i,j)*b(i) end do s=s/w(j) endif tmp(j)=s end do do j=1,p s=0. do jj=1,minnp s=s+v(j,jj)*tmp(jj) end do x(j)=s end do return end C subroutine caldat(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 : FORLIB.FVS * C * ROUTINE NAME: CALDAT * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Converts Julian Day to (year, month, day, hour, minute, * C * second) (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 ****************************************************************************** real*8 julian real*4 sec integer*4 iye,mon,idy,ihr,min,ifail integer*4 igreg,jalpha,ja,jb,jc,jd,je,ijul parameter (igreg=2299161) ! Cross-over to Gregorian Calendar if(julian.lt.0.d0) then ! Negative Julian Day not allowed ifail=1 return endif ijul=idint(julian) ! Integral Julian Day sec=sngl((julian-dble(ijul))*86400.d0)! Seconds from beginning of Jul. Day if(sec.ge.43200.) then ! In next calendar day ijul=ijul+1 sec=sec-43200. ! Adjust from noon to midnight else ! In same calendar day sec=sec+43200. ! Adjust from noon to midnight endif if(sec.ge.86400.) then ! Final check to prevent time 24:00:00 ijul=ijul+1 sec=sec-86400. endif min=int(sec/60.) ! Integral minutes from beginning of day sec=sec-float(min*60) ! Seconds from beginning of minute ihr=min/60 ! Integral hours from beginning of day min=min-ihr*60 ! Integral minutes from beginning of hour if(ijul.ge.igreg)then ! Correction for Gregorian Calendar jalpha=idint((dble(ijul-1867216)-0.25d0)/36524.25d0) ja=ijul+1+jalpha-idint(0.25d0*dble(jalpha)) else ! No correction ja=ijul endif jb=ja+1524 jc=idint(6680.d0+(dble(jb-2439870)-122.1d0)/365.25d0) jd=365*jc+idint(0.25d0*dble(jc)) je=idint(dble(jb-jd)/30.6001d0) idy=jb-jd-idint(30.6001d0*dble(je)) mon=je-1 if(mon.gt.12)mon=mon-12 iye=jc-4715 if(mon.gt.2)iye=iye-1 if(iye.le.0)iye=iye-1 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 : 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 dreadfil2(nin,anot,nchar,var,ifail) C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : UTILITY REF:JRH:27:01:1995 * C * * C * REVISION : ------------- JRH:--:--:1995 * C * * C * SOURCE : forlib.f * C * ROUTINE NAME: dreadfil2 * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Reads REAL*8 variable from file based on keyword. * C * (Modified version of areadfil2) * C * * C * nin ..... File input device * C * anot .... Keyword in file * C * nchar ... Number of characters in ANOT (max. 40) * C * var ..... Resultant REAL*8 variable * C * ifail ... 0 for successful execution, otherwise 1 * C * * C ****************************************************************************** real*8 var integer*4 nin,nchar,ifail character*(*) anot C integer*4 ios logical found character*80 buff C rewind(nin) ios=0 found=.false. do while(ios.eq.0.and..not.found) read(nin,1,iostat=ios) buff 1 format(a80) if(buff(1:nchar).eq.anot.and. $ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found read(buff(nchar+2:80),*,iostat=ios) var found=.true. endif end do if(ios.eq.0) then ifail=0 ! Match found and no error else ifail=1 ! No match found 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 : forlib.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 : FORLIB.FVS * C * ROUTINE NAME: IDINT2 * C * TYPE : FUNCTION * C * * C * FUNCTION : As INT, but trunctates towards -(infinity) * C This is D.P. version of INT2 * C * * 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 ireadfil2(nin,anot,nchar,ivar,ifail) C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : UTILITY REF:JRH:06:01:1995 * C * * C * REVISION : ------------- JRH:--:--:1995 * C * * C * SOURCE : forlib.f * C * ROUTINE NAME: ireadfil2 * C * TYPE : SUBROUTINE * C * * C * FUNCTION : Reads INTEGER*4 variable from file based on keyword. * C * (Modified version of ireadfil, without output to operator * C * or log file) * C * * C * nin ..... File input device * C * anot .... Keyword in file * C * nchar ... Number of characters in ANOT (max. 40) * C * ivar .... Resultant INTEGER*4 variable * C * ifail ... 0 for successful execution, otherwise 1 * C * * C ****************************************************************************** integer*4 nin,nchar,ivar,ifail character*(*) anot C integer*4 ios logical found character*80 buff C rewind(nin) ios=0 found=.false. do while(ios.eq.0.and..not.found) read(nin,1,iostat=ios) buff 1 format(a80) if(buff(1:nchar).eq.anot.and. $ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found read(buff(nchar+2:80),*,iostat=ios) ivar found=.true. endif end do if(ios.eq.0) then ifail=0 ! Match found and no error else ifail=1 ! No match found endif 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 : FORLIB.FVS * 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 ****************************************************************************** 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 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 integer function lhschar(c) C ****************************************************************************** C * * C * FORTRAN SOURCE CODE * C * * C * PROGRAM SET : UTILITY REF:JRH:07:09:1989 * C * * C * REVISION : ------------- JRH:--:--:1989 * C * * C * SOURCE : FORLIB.FOR * C * ROUTINE NAME: LHSCHAR * C * TYPE : INTEGER FUNCTION * C * * C * FUNCTION : Returns index of first non-blank character of CHARACTER * C * variable (defined by removing blank characters from * C * left-hand side). * C * * C ****************************************************************************** character*(*) c itot=len(c) do i=1,itot if(c(i:i).ne.' ') go to 1 end do lhschar=0 ! String is all blanks return 1 lhschar=i return end C