C
C FORTRAN LIBRARY SOURCE CODE
C
C
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITIES REF:JRH:14:01:2001 *
C * *
C * SOURCE : forlib.f *
C * TYPE : LIBRARY *
C * *
C * FUNCTION : Library of Fortran utility routines *
C * *
C * CONTENTS LATEST VERSION *
C * ******** ************** *
C * *
C * aconv 12:05:89 *
C * aread 07:09:89 *
C * areadfil 17:06:88 *
C * areadfil2 06:01:95 *
C * banner 25:05:93 *
C * banner2 04:04:95 *
C * bisec 10:10:85 *
C * caldat 23:05:91 *
C * changecase 03:12:90 *
C * clip 23:09:86 *
C * clop 07:07:89 *
C * datin 19:03:91 *
C * degmin 22:08:89 *
C * deg2degmin 24:03:91 *
C * dgasdev3 31:03:95 *
C * dheapsort 26:08:92 *
C * div1d 01:04:87 *
C * div2d 01:08:86 *
C * dopen 31:12:86 *
C * dran3 24:03:95 *
C * dreadfil2 27:01:95 *
C * dskopen 02:05:90 *
C * dsvbksb 26:08:92 *
C * dsvdcmp 26:08:92 *
C * dte 04:07:85 *
C * dtu 23:05:91 *
C * d2rcalc 25:07:91 *
C * error_handler 06:01:95 *
C * etd 04:07:85 *
C * fft 02:09:86 *
C * fft2 25:11:87 *
C * fit2d 07:07:89 *
C * gasdev3 31:03:95 *
C * geoginp 24:04:91 *
C * geoginpbuff 21:05:90 *
C * gsnumber 14:07:88 *
C * gssymbol 07:11:90 *
C * heapsort 29:08:91 *
C * idint2 23:04:91 *
C * idplace 18:11:91 *
C * inbuf 15:09:93 *
C * inilink 11:10:91 *
C * inslink 23:03:91 *
C * inslink2 05:08:93 *
C * intpoint 17:08:94 *
C * int2 23:04:91 *
C * invmerc 14:10:88 *
C * invrob 12:11:98 *
C * iosproc 28:08:91 *
C * ipdig 03:12:85 *
C * iread 07:09:89 *
C * ireadfil 07:09:89 *
C * ireadfil2 06:01:95 *
C * isearch 28:03:95 *
C * julday 23:05:91 *
C * lenchar 30:09:96 *
C * lhschar 07:09:89 *
C * lpclr 03:12:85 *
C * lpload 30:12:86 *
C * lpplot 30:12:86 *
C * merc 14:10:88 *
C * mxcopy 14:01:86 *
C * mxgaus 14:01:86 *
C * nchar 31:12:86 *
C * newpen 13:07:89 *
C * newton 23:07:86 *
C * nnonspace 08:08:89 *
C * number 17:07:89 *
C * opdig 03:12:85 *
C * overlayplot 05:06:90 *
C * picalc 25:07:91 *
C * pinsd 23:08:89 *
C * plot 16:11:89 *
C * plota 11:08:89 *
C * plots 17:07:89 *
C * ran3 24:03:95 *
C * ra2xy 20:08:85 *
C * remlink 11:10:91 *
C * rostp 03:07:86 *
C * rread 07:09:89 *
C * rreadfil 07:09:89 *
C * rreadfil2 06:01:95 *
C * r2dcalc 25:07:91 *
C * sim 23:07:86 *
C * sima 15:05:86 *
C * simb 15:05:86 *
C * solvmx 14:01:86 (minor mod. 18:12:89) *
C * srect 14:05:90 *
C * svbksb 09:01:92 *
C * svdcmp 08:01:92 *
C * symbol 13:06:91 *
C * timin 14:01:2001 *
C * utd 23:05:91 *
C * writearray 07:07:89 *
C * writearray2 06:01:92 *
C * xy2ra 02:07:85 *
C * ydif 31:07:86 *
C * ydif2d 29:09:86 *
C * yesno 19:02:87 *
C * yint 30:07:86 *
C * yint2d 29:09:86 *
C * yrect 22:09:86 *
C * yrect2d 29:09:86 *
C * *
C ******************************************************************************
CBack up
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
CBack up
subroutine aread(nin,nout,nlog,anot,acc,aarray,istart,istop,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:04:09:1986 *
C * *
C * REVISION : Mod. for output of message only if *
C * ISTART=0 or ISTOP=0 JRH:12:09:1986 *
C * Minor mods. JRH:31:12:1986 *
C * Minor mod. JRH:05:01:1987 *
C * Minor mod. JRH:17:05:1988 *
C * Mod. so that FORM is input without *
C * brackets JRH:21:06:1988 *
C * Mod. for compilation on transputer JRH:07:09:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: AREAD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads in a set of CHARACTER*40 variables *
C * *
C * NIN ..... Operator input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Annotation (40 chars.) *
C * ACC ..... Carriage-control character *
C * (' ' for new line, '$' for same line) *
C * AARRAY .. Resultant CHARACTER*40 array *
C * ISTART .. Start index of AARRAY *
C * ISTOP ... Stop index of AARRAY *
C * FORM .... Format for output to log device (40 chars., *
C * without brackets) *
C * *
C ******************************************************************************
character*40 aarray(*)
character*40 anot,form
character*1 acc
character*42 formt ! @7/9/89
logical cont ! @31/12/86
save cont ! @31/12/86
character*12 formout ! &5/1/87
data cont/.false./ ! @11/5/88
if(cont) then ! Continuation of terminal text @31/12/86
write(formout,1) ' ',acc ! &31/12/86
1 format('(',a1,'X,A40,A1',a1,')') ! &5/1/87
else ! Normal text @31/12/86
write(formout,1) '/',acc ! @31/12/86
endif ! @31/12/86
write(nout,formout) anot,' ' ! &5/1/87
write(nlog,formout) anot,' ' ! &5/1/87
if(istart.ne.0..and.istop.ne.0) then ! &31/12/86
read(nin,2) (aarray(i),i=istart,istop)
2 format((a40))
formt='('//form//')' ! @7/9/89
write(nlog,formt) (aarray(i),i=istart,istop) ! &21/6/88 &7/9/89
cont=.false. ! Terminal text not to be continued @31/12/86
else ! @31/12/86
cont=.true. ! Terminal text to be continued @31/12/86
endif ! @12/9/86
return
end
CBack up
subroutine areadfil(nin,nout,nlog,anot,nchar,avar,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:06:01:1987 *
C * *
C * REVISION : Mod. for keyword checking JRH:19:01:1987 *
C * Simplification JRH:17:06:1988 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: AREADFIL *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads CHARACTER*40 variable from file based on keyword. *
C * *
C * NIN ..... File input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Keyword in file *
C * NCHAR ... Number of characters in ANOT (max. 40) *
C * AVAR .... Resultant CHARACTER*40 variable *
C * FORM .... Format for output (40 chars., without brackets) *
C * *
C ******************************************************************************
character*(*) anot
character*40 avar
character*40 form
character*80 buff
character*40 blank ! &17/6/88
character*51 formt2
character*14 formt3
data blank/' '/
write(formt2,6) nchar,form
6 format('(1X,A',i2,',A3,',a40,')')
write(formt3,5) nchar
5 format('(/A9,A',i2,',A34/)')
rewind(nin)
4 read(nin,2,end=3) buff
2 format(a80)
if(buff(1:nchar).eq.anot.and. ! &17/6/88
$ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found &19/1/87
avar(1:40)=buff(nchar+2:nchar+41)
write(nout,formt2) anot,' = ',avar
write(nlog,formt2) anot,' = ',avar
return
else
go to 4 ! Read more data
endif
3 write(nout,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
write(nlog,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
stop
end
CBack up
subroutine areadfil2(nin,anot,nchar,avar,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: areadfil2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads CHARACTER*40 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*40 variable *
C * ifail ... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 nin,nchar,ifail
character*(*) anot
character*40 avar
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
avar(1:40)=buff(nchar+2:nchar+41)
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
CBack up
subroutine banner(nout,nlog,nline,lines)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:22:05:1987 *
C * *
C * REVISION : Changed JMG to CSIRO JRH:06:07:1989 *
C * Mods. to licence conditions JRH:15:07:1992 *
C * Mods. to make it more general JRH:25:05:1993 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: BANNER *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Writes program banner to terminal and to log file. *
C * *
C * NOUT ..... Operator output LU *
C * NLOG ..... Log output LU *
C * NLINE .... No. of input lines *
C * LINES .... NLINE lines of text (each of 40 chars.) *
C * *
C ******************************************************************************
character*40 lines(*)
write(nout,1)
write(nlog,1)
1 format(////11x,
$ '************************************************************'
$ /11x,
$ '* *')
write(nout,2) (lines(i),i=1,nline)
write(nlog,2) (lines(i),i=1,nline)
2 format((11x,'* ',a40,' *'/11x,
$ '* *'))
nline=0
write(nout,3)
write(nlog,3)
3 format(11x,
$ '* *'
$ /11x,
$ '* This software is supplied subject to the terms of the *'! &15/7/92
$ /11x,
$ '* licence agreement supplied with it. The terms may be *'! &15/7/92
$ /11x,
$ '* summarised as : *'! &15/7/92
$ /11x,
$ '* *'
$ /11x,
$ '* 1. The software and associated manual(s) are *'! &15/7/92
$ /11x,
$ '* supplied with no warranty, expressed or implied. *'! &15/7/92
$ /11x,
$ '* *'
$ /11x,
C Following 5 lines removed 25/5/93 :
C $ '* 2. The software is supplied for non-profit *'! &15/7/92
C $ /11x,
C $ '* teaching and research purposes only. Sufficient *'! &15/7/92
C $ /11x,
C $ '* copies may be made for these purposes only. *'! &15/7/92
$ '* 2. The software is supplied for the Licensee''s *'!@25/5/93
$ /11x, ! @25/5/93
$ '* computing activities only. Sufficient copies *'! @25/5/93
$ /11x, ! @25/5/93
$ '* may be made for these purposes only. *'! @25/5/93
$ /11x,
$ '* *'
$ /11x,
C Following 5 lines removed 25/5/93 :
C $ '* 3. Any modifications to the software source code *'! &15/7/92
C $ /11x,
C $ '* should carry a prominent notice stating the *'! &15/7/92
C $ /11x,
C $ '* change and the date of that change. *'! &15/7/92
$ '* 3. Any modifications to the software source code *'! @25/5/93
$ /11x, ! @25/5/93
$ '* should be approved by CSIRO and should carry a *'! @25/5/93
$ /11x, ! @25/5/93
$ '* prominent notice stating the change and the date *'! @25/5/93
$ /11x, ! @25/5/93
$ '* of that change. *'! @25/5/93
$ /11x,
$ '* *'
$ /11x,
$ '* 4. This banner is not to be modified and its *'! &15/7/92
$ /11x,
$ '* display is not to be inactivated. *'! &15/7/92
$ /11x,
$ '* *'
$ /11x,
$ '* OCEAN MODELLING GROUP, *' ! &6/7/89
$ /11x,
$ '* CSIRO DIVISION OF OCEANOGRAPHY, *' ! &6/7/89
$ /11x,
$ '* HOBART, TASMANIA, 7000. *' ! &6/7/89
$ /11x,
$ '* *'
$ /11x,
$ '************************************************************')
return
end
CBack up
subroutine banner2(nout,nline,lines)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:04:04:1995 *
C * *
C * REVISION : ---------- JRH:--:--:1995 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: banner2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Writes program banner to terminal. *
C * *
C * nout ..... Operator output LU *
C * nline .... No. of input lines *
C * lines .... NLINE lines of text (each of 40 chars.) *
C * *
C ******************************************************************************
integer*4 nout,nline
character*40 lines(*)
C
write(nout,1)
1 format(////11x,
$ '************************************************************'
$ /11x,
$ '* *')
write(nout,2) (lines(i),i=1,nline)
2 format((11x,'* ',a40,' *'/11x,
$ '* *'))
nline=0
write(nout,3)
3 format(11x,
$ '* *'
$ /11x,
$ '* This software is supplied subject to the terms of the *'
$ /11x,
$ '* licence agreement supplied with it. The terms may be *'
$ /11x,
$ '* summarised as : *'
$ /11x,
$ '* *'
$ /11x,
$ '* 1. The software and associated manual(s) are *'
$ /11x,
$ '* supplied with no warranty, expressed or implied. *'
$ /11x,
$ '* *'
$ /11x,
$ '* 2. The software is supplied for the Licensee''s *'
$ /11x,
$ '* computing activities only. Sufficient copies *'
$ /11x,
$ '* may be made for these purposes only. *'
$ /11x,
$ '* *'
$ /11x,
$ '* 3. Any modifications to the software source code *'
$ /11x,
$ '* should be approved by CSIRO and should carry a *'
$ /11x,
$ '* prominent notice stating the change and the date *'
$ /11x,
$ '* of that change. *'
$ /11x,
$ '* *'
$ /11x,
$ '* 4. This banner is not to be modified and its *'
$ /11x,
$ '* display is not to be inactivated. *'
$ /11x,
$ '* *'
$ /11x,
$ '* OCEAN MODELLING GROUP, *'
$ /11x,
$ '* CSIRO DIVISION OF OCEANOGRAPHY, *'
$ /11x,
$ '* HOBART, TASMANIA, 7000. *'
$ /11x,
$ '* *'
$ /11x,
$ '************************************************************')
return
end
CBack up
subroutine bisec(f,y,xmin,xmax,nit1,nit2,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:10:10:1986 *
C * *
C * REVISION : Speeded up JRH:29:06:1997 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: BISEC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Finds approximate solution of Y = F(X) using the ,method *
C * of bisection. Initially range XMIN to XMAX is expanded *
C * by a factor of 3, up to NIT1 times, until the range *
C * brackets the required solution. The range XMIN to XMAX *
C * is then reduced by a factor of 2, NIT2 times, ensuring *
C * the range still brackets the required solution. *
C * *
C * F ..... FUNCTION of X *
C * Y ..... Value of F(X) for solution *
C * XMIN .. Min. value of initial range *
C * XMAX .. Max. value of initial range *
C * NIT1 .. Max. number of initial expansions *
C * NIT2 .. Number of bisections *
C * IFAIL . 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
do i=1,nit1
fmin=f(xmin)
fmax=f(xmax)
if((fmin.le.y.and.fmax.gt.y).or.
$ (fmin.gt.y.and.fmax.le.y)) go to 1
xdum=2.*xmin-xmax
xmax=2.*xmax-xmin
xmin=xdum
end do
ifail=1 ! Could not find suitable range
return
1 do i=1,nit2
xdum=(xmin+xmax)/2.
fdum=f(xdum)
if((fmin.le.y.and.fdum.gt.y).or.
$ (fmin.gt.y.and.fdum.le.y)) then
xmax=xdum
C if(i.ne.nit2) fmax=f(xmax) ! &29/6/97
if(i.ne.nit2) fmax=fdum ! @29/6/97
else
xmin=xdum
C if(i.ne.nit2) fmin=f(xmin) ! &29/6/97
if(i.ne.nit2) fmin=fdum ! @29/6/97
endif
end do
ifail=0
return
end
CBack up
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
CBack up
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
CBack up
subroutine clip(x,y,xlc,xrc,ybc,ytc,lc,lplot)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:23:09:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : PLT1199.FOR *
C * ROUTINE NAME: CLIP *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Clips line (X(1),Y(1)) to (X(2),Y(2)) to rectangular *
C * window defined by X = XLC to XRC and Y = YBC to YTC. *
C * If point clipped, LC() returned as .TRUE.. *
C * If line completely off plot, LPLOT returned as .FALSE.. *
C * *
C ******************************************************************************
dimension x(2),y(2)
logical lc(2),lplot
logical lout(4,2)
do i=1,2
lout(1,i)=(x(i).lt.xlc)
lout(2,i)=(x(i).gt.xrc)
lout(3,i)=(y(i).lt.ybc)
lout(4,i)=(y(i).gt.ytc)
lc(i)=.false.
end do
do i=1,2
3 if(lout(1,i).or.lout(2,i).or.lout(3,i).or.lout(4,i)) then
if((lout(1,1).and.lout(1,2)).or.
$ (lout(2,1).and.lout(2,2)).or.
$ (lout(3,1).and.lout(3,2)).or.
$ (lout(4,1).and.lout(4,2))) then
lplot=.false. ! Off plot
return
endif
if(lout(1,i)) then ! Crosses LH edge
y(i)=y(1)+(y(2)-y(1))*(xlc-x(1))/(x(2)-x(1))
x(i)=xlc
lc(i)=.true.
else
if(lout(2,i)) then ! Crosses RH edge
y(i)=y(1)+(y(2)-y(1))*(xrc-x(1))/(x(2)-x(1))
x(i)=xrc
lc(i)=.true.
else
if(lout(3,i)) then ! Crosses bottom edge
x(i)=x(1)+(x(2)-x(1))*(ybc-y(1))/(y(2)-y(1))
y(i)=ybc
lc(i)=.true.
else
if(lout(4,i)) then ! Crosses top edge
x(i)=x(1)+(x(2)-x(1))*(ytc-y(1))/(y(2)-y(1))
y(i)=ytc
lc(i)=.true.
endif
endif
endif
endif
lout(1,i)=(x(i).lt.xlc)
lout(2,i)=(x(i).gt.xrc)
lout(3,i)=(y(i).lt.ybc)
lout(4,i)=(y(i).gt.ytc)
go to 3
endif
end do
lplot=.true.
return ! Visible line is (X(1),Y(1)) to (X(2),Y(2))
end
CBack up
subroutine clop(ndev)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:09:11:1988 *
C * *
C * REVISION : Removal of "shared" option in "open" JRH:06:07:1989 *
C * "name" changed to "file" in open *
C * statement JRH:07:07:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: CLOP *
C * TYPE : MAIN *
C * *
C * FUNCTION : Closes and re-opens file on device NDEV, so as to empty *
C * buffer when sharing file. THIS SHOULD ONLY BE USED WITH *
C * SEQUENTIAL FILES. *
C * *
C ******************************************************************************
character*10 formt
character*80 filnam
inquire(unit=ndev,form=formt,name=filnam,recl=nr)
close(unit=ndev)
open(unit=ndev,access='APPEND',form=formt,file=filnam,recl=nr, ! &7/7/89
$ status='UNKNOWN') ! &6/7/89
return
end
CBack up
subroutine datin(nin,nout,nlog,anot,iye,mon,idy)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:19:03:1991 *
C * *
C * REVISION : Minor mod. to formats JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: DATIN *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads in a date, in various formats *
C * *
C * NIN ..... Operator input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Annotation (40 chars.) *
C * IYE ..... Year *
C * MON ..... Month *
C * IDY ..... Day *
C * *
C * NOTE : Main constraint is that day is input BEFORE year *
C * *
C ******************************************************************************
integer*4 nin,nout,nlog,iye,mon,idy
character*40 anot
integer*4 idmon(12)
integer*4 i,idum,ind
logical alpha
character*80 buff
character*3 amon(12)
data amon/'JAN','FEB','MAR','APR','MAY','JUN',
$ 'JUL','AUG','SEP','OCT','NOV','DEC'/
data idmon/ 31 , 28 , 31 , 30 , 31 , 30 ,
$ 31 , 31 , 30 , 31 , 30 , 31 /
4 write(nout,1) anot
write(nlog,1) anot
1 format(1x,a40)
read(nin,2,err=4) buff(1:80)
2 format(a80)
write(nlog,3) buff(1:80)
3 format(1x,a80)
call changecase(buff(1:80),buff(1:80),80,1) ! Change to upper case
alpha=.false.
do i=1,12
ind=index(buff(1:80),amon(i))
if(ind.ne.0) then ! Found a matching month
mon=i
buff(ind:ind+2)=' ' ! Delete month
alpha=.true. ! Alphabetic month
go to 6
endif
end do
6 do i=1,80 ! First remove all non-numeric characters
idum=ichar(buff(i:i))
if(idum.lt.48.or.idum.gt.57) buff(i:i)=' '
end do
if(alpha) then ! Alphabetic month
read(buff(1:80),*,err=4) idy,iye
else
read(buff(1:80),*,err=4) idy,mon,iye
endif
if(iye.le.99) iye=iye+1900 ! Cope with two digit year (assume 20th C.)
if(iye.le.1900.or.iye.ge.2100) go to 4 ! Set range of years
if(mon.lt.1.or.mon.gt.12) go to 4 ! Check range of months
if(mon.eq.2) then ! Special check for Feb.
if((iye/4)*4.eq.iye) then ! Leap year
if(idy.lt.1.or.idy.gt.29) go to 4
else ! Non-leap year
if(idy.lt.1.or.idy.gt.28) go to 4
endif
else
if(idy.lt.1.or.idy.gt.idmon(mon)) go to 4 ! Error ..... re-input
endif
return
end
CBack up
character*29 function degmin(ang)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:27:10:1988 *
C * *
C * REVISION : Calculation of ANG/D2R outside *
C * intrinsic functions (otherwise Sun *
C * sometimes gives different results for *
C * IDINT(ANG/D2R) and DINT(ANG/D2R)) JRH:22:08:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: DEGMIN *
C * TYPE : CHARACTER*29 FUNCTION *
C * *
C * FUNCTION : Converts ANG (in radians) to DEGMIN (CHARACTER*29) as *
C * degrees and decimal minutes. *
C * *
C ******************************************************************************
real*8 ang,d2r,ddum ! &22/8/89
character*1 asign
d2r=datan(1.d0)/45.d0 ! Degrees to radians
if(ang.ge.0.) then
asign=' '
else
asign='-'
endif
ddum=ang/d2r ! @22/8/89
write(degmin,1) asign,iabs(idint(ddum)), ! @22/8/89
$ (dabs(ddum)-dabs(dint(ddum)))*60.d0 ! @22/8/89
1 format(a1,i3,' degrees, ',f7.4,' minutes')
return
end
CBack up
subroutine deg2degmin(ang,atype,ndec,ideg,min1,min2,adir)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:24:03:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: DEG2DEGMIN *
C * TYPE : SUBOUTINE *
C * *
C * FUNCTION : Converts lat. or long. (degrees) to degrees and decimal *
C * minutes, and alphabetic direction *
C * *
C * ANG ..... input lat. or long. (degrees) *
C * ATYPE ... "LONG" for longitude or "LAT " for latitude *
C * NDEC .... required no. of decimal places for minutes *
C * IDEG .... output degrees *
C * MIN1 .... output integral minutes *
C * MIN2 .... output fractional minutes *
C * ADIR .... 'E' (+ve) or 'W' (-ve) for longitude or *
C * 'N' (+ve) or 'S' (-ve) for latitude *
C * *
C ******************************************************************************
real*4 ang,rmin
integer*4 ndec,ideg,min1,min2
character*4 atype
character*1 adir
ideg=int(abs(ang))
rmin=(abs(ang)-float(ideg))*60.
min1=int(rmin)
min2=nint((rmin-float(min1))*(10**ndec))
if(min2.eq.10**ndec) then
min2=0
min1=min1+1
endif
if(min1.eq.60) then
min1=0
ideg=ideg+1
endif
if(atype.eq.'LONG'.or.atype.eq.'long') then
if(ang.gt.0.) then
adir='E'
else
adir='W'
endif
else if(atype.eq.'LAT '.or.atype.eq.'lat ') then
if(ang.gt.0.) then
adir='N'
else
adir='S'
endif
else
adir=' ' ! Leave as blank
endif
return
end
CBack up
real*8 function dgasdev3(idum)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:31:03:1995 *
C * *
C * REVISION : ------------- JRH:--:--:1995 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: dgasdev3 *
C * TYPE : real*8 function *
C * *
C * FUNCTION : Random number generator yielding normal distribution *
C * with zero mean and unit variance (real*8 version based on *
C * Numerical Recipes). *
C * *
C ******************************************************************************
integer*4 idum
C
real*8 v1,v2,r,fac,gset
real*8 dran3
integer*4 iset
data iset/0/
save iset,gset
if (iset.eq.0) then
1 v1=2.d0*dran3(idum)-1.d0
v2=2.d0*dran3(idum)-1.d0
r=v1**2+v2**2
if(r.ge.1.d0)go to 1
fac=dsqrt(-2.d0*dlog(r)/r)
gset=v1*fac
dgasdev3=v2*fac
iset=1
else
dgasdev3=gset
iset=0
endif
return
end
CBack up
subroutine dheapsort(ain,aout,ind,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:26:08:1992 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB .FVS *
C * ROUTINE NAME: DHEAPSORT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Heapsort based on Numerical Recipes subroutine SORT2 *
C * - D.P. version *
C * *
C * AIN ..... Array to be sorted (unchanged on output) *
C * AOUT .... Resultant sorted array (increasing) *
C * IND ..... Array of sorted indices *
C * N ....... Size of arrays AIN and AOUT *
C * *
C ******************************************************************************
real*8 ain(*),aout(*)
integer*4 ind(*)
integer*4 n
real*8 rra
integer*4 i,j,l,ir,idum
l=n/2+1
ir=n
do i=1,n
aout(i)=ain(i) ! Copy input array to output array
ind(i)=i ! Generate initial idum array
end do
if(n.eq.1) return ! Special for only one record
10 continue
if(l.gt.1)then
l=l-1
rra=aout(l)
idum=ind(l)
else
rra=aout(ir)
idum=ind(ir)
aout(ir)=aout(1)
ind(ir)=ind(1)
ir=ir-1
if(ir.eq.1)then
aout(1)=rra
ind(1)=idum
return
endif
endif
i=l
j=l+l
20 if(j.le.ir)then
if(j.lt.ir)then
if(aout(j).lt.aout(j+1))j=j+1
endif
if(rra.lt.aout(j))then
aout(i)=aout(j)
ind(i)=ind(j)
i=j
j=j+j
else
j=ir+1
endif
go to 20
endif
aout(i)=rra
ind(i)=idum
go to 10
end
CBack up
subroutine div1d(ain,imaxin,aout,imaxout,ifac,zero,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:01:08:1986 *
C * *
C * REVISION : Minor mod. JRH:01:04:1987 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: DIV1D *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Sub-divides 1-D array of data using parabolic *
C * interpolation. *
C * *
C * AIN(I) ..... input array *
C * IMAXIN ..... max. index of AIN for calculations *
C * AOUT(I) .... output array *
C * IMAXOUT .... max. index of AOUT for calculations *
C * IFAC ....... division factor in I *
C * ZERO ....... "base" level of AIN *
C * IFAIL ...... 0 on exit if successful *
C * 1 on exit if unsuccessful *
C * *
C ******************************************************************************
dimension ain(*),aout(*)
do ii=1,imaxout
aout(ii)=zero ! Clear AOUT to base level
end do
imaxd=(imaxin-1)*ifac+1 ! Required space in AOUT
if(imaxd.gt.imaxout) then
ifail=1
return
endif
do i=1,imaxin ! Step through I
dum=ain(i)-zero ! Values relative to base
if(dum.lt.0.) dum=0. ! Set values below base to base level
ain(i)=dum
end do
aout(1)=ain(1)+zero
do i=1,imaxin-1 ! Step through I in AIN
d1l=d1r ! Use "old" RH derivative
d2l=d2r ! " " " "
if(i.ne.imaxin-1) then ! Not at RH end
d1r=ain(i+2)-ain(i) ! 1st derivative, RH curve
d2r=ain(i+2)+ain(i)-2*ain(i+1) ! 2nd " " "
endif
ioff=(i-1)*ifac+1
do ii=1,ifac ! Step within an interpolation interval
xl=float(ii)/float(ifac)
xr=xl-1.
c1l=xl/2.
c2l=xl*xl/2.
c1r=xr/2.
c2r=xr*xr/2.
C Note checks for interpolation in vicinity of boundary of "real data"
C and "no data" (the latter to be output at the level of the bottom of
C the base) - following scheme preserves sharp step here .....
if(ain(i).eq.0..or.ain(i+1).eq.0.) then ! One end at bottom of base
if(ii.ne.ifac) then
dum=0. ! Set all to bottom of base
else
dum=ain(i+1) ! except for RH end
endif
else
if(i.eq.1) then ! LH end
dum=ain(i+1)+d1r*c1r+d2r*c2r ! RH parabola
else
if(i.eq.imaxin-1) then ! RH end
dum=ain(i)+d1l*c1l+d2l*c2l ! LH parabola
else ! Other
if(ain(i-1).eq.0.) then ! Bottom of base is to left
if(ain(i+2).eq.0.) then ! Bottom of base is also to right
dum=-ain(i)*xr+ain(i+1)*xl ! Linear interp.
else ! Bottom of base is not to right
dum=ain(i+1)+d1r*c1r+d2r*c2r ! RH parabola
endif
else
if(ain(i+2).eq.0.) then ! Bottom of base is to right
dum=ain(i)+d1l*c1l+d2l*c2l ! LH parablola
else
dum=(ain(i)+d1l*c1l ! LH parabola
$ +d2l*c2l)*(-xr) ! LH parabola
$ +(ain(i+1)+d1r*c1r ! RH parabola
$ +d2r*c2r)*xl ! RH parabola
endif
endif
endif
endif
endif
if(dum.lt.0.) dum=0.
aout(ii+ioff)=dum+zero
end do
end do
return
end
CBack up
subroutine div2d(ain,imaxin,imaxina,jmaxin,
$ aout,imaxout,imaxouta,jmaxout,
$ ifac,jfac,zero,ifail,
$ work1,work2,work3,work4,work5)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:01:08:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: DIV2D *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Sub-divides 2-D array of data using parabolic *
C * interpolation. *
C * *
C * AIN(I,J) ... input array *
C * IMAXIN ..... max. index of AIN for calculations *
C * IMAXINA .... first dimension of AIN *
C * JMAXIN ..... max. index of AIN for calculations *
C * AOUT(I,J) .. output array *
C * IMAXOUT .... max. index of AOUT for calculations *
C * IMAXOUTA ... first dimension of AOUT *
C * JMAXOUT .... max. index of AOUT for calculations *
C * IFAC ....... division factor in I *
C * JFAC ....... division factor in J *
C * ZERO ....... "base" level of AIN *
C * IFAIL ...... 0 on exit if successful *
C * 1 on exit if unsuccessful *
C * WORK1 ...... workspace array (dim. at least IMAXIN) *
C * WORK2 ...... workspace array (dim. at least IMAXOUT) *
C * WORK3....... workspace array (dim. at least JMAXIN) *
C * WORK4....... workspace array (dim. at least JMAXOUT) *
C * WORK5 ...... workspace array (first dim. = IMAXOUTA, *
C * second dim. at least JMAXIN) *
C * *
C ******************************************************************************
dimension ain(imaxina,*),aout(imaxouta,*)
dimension work1(*),work2(*),work3(*)
dimension work4(*),work5(imaxouta,*)
do j=1,jmaxin
do i=1,imaxin
work1(i)=ain(i,j)
end do
call div1d(work1,imaxin,work2,imaxout,ifac,zero,ifail)
if(ifail.eq.1) return
do i=1,imaxout
work5(i,j)=work2(i)
end do
end do
do i=1,imaxout
do j=1,jmaxin ! &4/9/86
work3(j)=work5(i,j)
end do
call div1d(work3,jmaxin,work4,jmaxout,jfac,zero,ifail)
if(ifail.eq.1) return
do j=1,jmaxout
aout(i,j)=work4(j)
end do
end do
return
end
CBack up
subroutine dopen(nin,nout,ndev,anot,stat)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:02:07:1985 *
C * *
C * REVISION : Minor mod. JRH:31:12:1986 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: DOPEN *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Interactive routine for opening disc file. *
C * *
C ******************************************************************************
character*40 anot
character*10 stat
character*20 filnam
write(nout,1) anot
1 format(/1x,a40)
read(nin,2) filnam
2 format(a20)
open(unit=ndev,file=filnam,status=stat)
return
end
CBack up
real*8 function dran3(idum)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:24:03:1995 *
C * *
C * REVISION : Change in order of save and data *
C * statements for safety JRH:29:03:2001 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: dran3 *
C * TYPE : real*8 function *
C * *
C * FUNCTION : Random number generator (real*8 version based on *
C * Numerical Recipes). *
C * *
C ******************************************************************************
integer*4 idum
C
real*8 fac
integer*4 mbig,mseed,mz
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.d-9)
integer*4 ma(55)
integer*4 iff,inext,inextp,mj,mk,i,ii,k
save inext,inextp,ma,iff ! Mod. to Numerical Recipes code
data iff /0/
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do 11 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
11 continue
do 13 k=1,4
do 12 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12 continue
13 continue
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
dran3=mj*fac
return
end
CBack up
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
CBack up
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
CBack up
subroutine dsvbksb(u,w,v,m,n,mp,np,b,x,tmp)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:26:08:1992 *
C * *
C * REVISION : ------------- JRH:--:--:1992 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: DSVBKSB *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Singular value decomposition solution *
C * (from Numerical Recipes - D.P. version of SVBKSB) *
C * *
C * U ..... Input M by N matrix *
C * W ..... Input N element weight vector *
C * V ..... Input N by N matrix ( NOT the transpose of "V") *
C * M ..... Number of rows of U (M.ge.N) *
C * N ..... Number of columns of U (M.ge.N) *
C * MP .... Physical number of rows of U *
C * NP .... Physical number of columns of U *
C * B ..... Input M element vector (right hand side) *
C * X ..... Output N element solution vector *
C * TMP ... Workspace N element vector *
C * *
C ******************************************************************************
real*8 u(mp,np),w(np),v(np,np),b(mp),x(np),tmp(*)
integer*4 m,n,mp,np
real*8 s
integer*4 i,j,jj
do 12 j=1,n
s=0.
if(w(j).ne.0.)then
do 11 i=1,m
s=s+u(i,j)*b(i)
11 continue
s=s/w(j)
endif
tmp(j)=s
12 continue
do 14 j=1,n
s=0.
do 13 jj=1,n
s=s+v(j,jj)*tmp(jj)
13 continue
x(j)=s
14 continue
return
end
CBack up
subroutine dsvdcmp(a,m,n,mp,np,u,w,v,rv1,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:26:08:1992 *
C * *
C * REVISION : ------------- JRH:--:--:1992 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: DSVDCMP *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Singular value decomposition solution *
C * (from Numerical Recipes - D.P. version of SVDCMP) *
C * *
C * A ..... Input M by N matrix *
C * M ..... Number of rows of A (M.ge.N) *
C * N ..... Number of columns of A (M.ge.N) *
C * MP .... Physical number of rows of A *
C * NP .... Physical number of columns of A *
C * U ..... Output M by N matrix *
C * W ..... Output N element weight vector *
C * V ..... Output N by N matrix (NOT the transpose of "V") *
C * RV1 ... Workspace N element vector *
C * IFAIL . 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
real*8 a(mp,np),u(mp,np),w(np),v(np,np),rv1(*)
integer*4 m,n,mp,np,ifail
real*8 f,g,h,scale,anorm,s,c,x,y,z
integer*4 i,j,k,l,its,nm
if(m.lt.n) then
ifail=1
return
endif
g=0.0
scale=0.0
anorm=0.0
do i=1,m ! Copy A to U
do j=1,n
u(i,j)=a(i,j)
end do
end do
do 25 i=1,n
l=i+1
rv1(i)=scale*g
g=0.0
s=0.0
scale=0.0
if (i.le.m) then
do 11 k=i,m
scale=scale+abs(u(k,i))
11 continue
if (scale.ne.0.0) then
do 12 k=i,m
u(k,i)=u(k,i)/scale
s=s+u(k,i)*u(k,i)
12 continue
f=u(i,i)
g=-sign(sqrt(s),f)
h=f*g-s
u(i,i)=f-g
if (i.ne.n) then
do 15 j=l,n
s=0.0
do 13 k=i,m
s=s+u(k,i)*u(k,j)
13 continue
f=s/h
do 14 k=i,m
u(k,j)=u(k,j)+f*u(k,i)
14 continue
15 continue
endif
do 16 k= i,m
u(k,i)=scale*u(k,i)
16 continue
endif
endif
w(i)=scale *g
g=0.0
s=0.0
scale=0.0
if ((i.le.m).and.(i.ne.n)) then
do 17 k=l,n
scale=scale+abs(u(i,k))
17 continue
if (scale.ne.0.0) then
do 18 k=l,n
u(i,k)=u(i,k)/scale
s=s+u(i,k)*u(i,k)
18 continue
f=u(i,l)
g=-sign(sqrt(s),f)
h=f*g-s
u(i,l)=f-g
do 19 k=l,n
rv1(k)=u(i,k)/h
19 continue
if (i.ne.m) then
do 23 j=l,m
s=0.0
do 21 k=l,n
s=s+u(j,k)*u(i,k)
21 continue
do 22 k=l,n
u(j,k)=u(j,k)+s*rv1(k)
22 continue
23 continue
endif
do 24 k=l,n
u(i,k)=scale*u(i,k)
24 continue
endif
endif
anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25 continue
do 32 i=n,1,-1
if (i.lt.n) then
if (g.ne.0.0) then
do 26 j=l,n
v(j,i)=(u(i,j)/u(i,l))/g
26 continue
do 29 j=l,n
s=0.0
do 27 k=l,n
s=s+u(i,k)*v(k,j)
27 continue
do 28 k=l,n
v(k,j)=v(k,j)+s*v(k,i)
28 continue
29 continue
endif
do 31 j=l,n
v(i,j)=0.0
v(j,i)=0.0
31 continue
endif
v(i,i)=1.0
g=rv1(i)
l=i
32 continue
do 39 i=n,1,-1
l=i+1
g=w(i)
if (i.lt.n) then
do 33 j=l,n
u(i,j)=0.0
33 continue
endif
if (g.ne.0.0) then
g=1.0/g
if (i.ne.n) then
do 36 j=l,n
s=0.0
do 34 k=l,m
s=s+u(k,i)*u(k,j)
34 continue
f=(s/u(i,i))*g
do 35 k=i,m
u(k,j)=u(k,j)+f*u(k,i)
35 continue
36 continue
endif
do 37 j=i,m
u(j,i)=u(j,i)*g
37 continue
else
do 38 j= i,m
u(j,i)=0.0
38 continue
endif
u(i,i)=u(i,i)+1.0
39 continue
do 49 k=n,1,-1
do 48 its=1,30
do 41 l=k,1,-1
nm=l-1
if ((abs(rv1(l))+anorm).eq.anorm) go to 2
if ((abs(w(nm))+anorm).eq.anorm) go to 1
41 continue
1 c=0.0
s=1.0
do 43 i=l,k
f=s*rv1(i)
if ((abs(f)+anorm).ne.anorm) then
g=w(i)
h=sqrt(f*f+g*g)
w(i)=h
h=1.0/h
c= (g*h)
s=-(f*h)
do 42 j=1,m
y=u(j,nm)
z=u(j,i)
u(j,nm)=(y*c)+(z*s)
u(j,i)=-(y*s)+(z*c)
42 continue
endif
43 continue
2 z=w(k)
if (l.eq.k) then
if (z.lt.0.0) then
w(k)=-z
do 44 j=1,n
v(j,k)=-v(j,k)
44 continue
endif
go to 3
endif
if (its.eq.30) then
ifail=1
return
endif
x=w(l)
nm=k-1
y=w(nm)
g=rv1(nm)
h=rv1(k)
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
g=sqrt(f*f+1.0)
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
c=1.0
s=1.0
do 47 j=l,nm
i=j+1
g=rv1(i)
y=w(i)
h=s*g
g=c*g
z=sqrt(f*f+h*h)
rv1(j)=z
c=f/z
s=h/z
f= (x*c)+(g*s)
g=-(x*s)+(g*c)
h=y*s
y=y*c
do 45 nm=1,n
x=v(nm,j)
z=v(nm,i)
v(nm,j)= (x*c)+(z*s)
v(nm,i)=-(x*s)+(z*c)
45 continue
z=sqrt(f*f+h*h)
w(j)=z
if (z.ne.0.0) then
z=1.0/z
c=f*z
s=h*z
endif
f= (c*g)+(s*y)
x=-(s*g)+(c*y)
do 46 nm=1,m
y=u(nm,j)
z=u(nm,i)
u(nm,j)= (y*c)+(z*s)
u(nm,i)=-(y*s)+(z*c)
46 continue
47 continue
rv1(l)=0.0
rv1(k)=f
w(k)=x
48 continue
3 continue
49 continue
ifail=0
return
end
CBack up
subroutine dte(iye,mon,idy,ihr,min,sec,tim)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : GENERAL REF:JRH:04:07:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: DTE *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts daytime to elapsed time *
C * *
C * Daytime : *
C * IYE ..... Year *
C * MON ..... Month *
C * IDY ..... Day *
C * IHR ..... Hour *
C * MIN ..... Minute *
C * SEC ..... Second *
C * *
C * Elapsed time : *
C * TIM ..... Time in seconds from beginning of year *
C * *
C * NOTE that current version (with D.P. (64 bit) *
C * internal computations, but S.P. (32 bit) TIM) will *
C * generally be correct to +/- 2 secs. However, a time *
C * in exact minutes (ie. SEC = 0.) is computed exactly. *
C * *
C ******************************************************************************
double precision dtim
dimension nd(12)
data nd/0,31,59,90,120,151,181,212,243,273,304,334/
dtim=((dble(float(nd(mon)+idy-1))*24.d0
$ +dble(float(ihr)))*60.d0
$ +dble(float(min)))*60.d0
$ +dble(sec)
C Leap year correction .....
if((iye/4)*4.eq.iye.and.mon.gt.2) dtim=dtim+86400.d0
tim=sngl(dtim)
return
end
CBack up
subroutine dtu(iye,mon,idy,ihr,min,sec,utim,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: DTU *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts daytime to Unix Time *
C * *
C * (Unix Time starts at 0000 on 1 Jan. 1970) *
C * *
C * IYE ........ Year *
C * MON ........ Month *
C * IDY ........ Day *
C * IHR ........ Hour *
C * MIN ........ Minute *
C * SEC ........ Second *
C * UTIM ....... Unix time (seconds) *
C * IFAIL ...... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
real*8 utim
real*4 sec
integer*4 iye,mon,idy,ihr,min,ifail
real*8 jd0,julian
logical first
save jd0,first
data first/.true./
if(first) then ! Compute zero of Unix Time in Julian days
call julday(1970,1,1,0,0,0.,jd0,ifail)
if(ifail.ne.0) return ! Error
first=.false.
endif
call julday(iye,mon,idy,ihr,min,sec,julian,ifail)
if(ifail.ne.0) return ! Error
utim=(julian-jd0)*86400.d0
return
end
CBack up
subroutine d2rcalc(d2r)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:25:07:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: D2RCALC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Returns the value of D2R (degrees to radians constant) *
C * *
C ******************************************************************************
real*4 d2r
d2r=atan(1.)*4./180.
return
end
CBack up
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
CBack up
subroutine etd(iye,mon,idy,ihr,min,sec,tim)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : GENERAL REF:JRH:04:07:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: ETD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts elapsed time to daytime. *
C * Inputs are IYE and TIM. *
C * IYE and TIM are modified if TIM is negative or *
C * greater than year length. *
C * *
C * Daytime : *
C * IYE ..... Year *
C * MON ..... Month *
C * IDY ..... Day *
C * IHR ..... Hour *
C * MIN ..... Minute *
C * SEC ..... Second *
C * *
C * Elapsed time : *
C * TIM ..... Time in seconds from beginning of year *
C * *
C * NOTE that current version (with D.P. (64 bit) *
C * internal computations, but S.P. (32 bit) TIM) will *
C * generally be correct to +/- 2 secs. However, a time *
C * in exact minutes (ie. SEC = 0.) is computed exactly. *
C * *
C ******************************************************************************
double precision dtim,dyl,didy,dihr,dmin
dimension nd(13,2)
data nd/0,31,59,90,120,151,181,212,243,273,304,334,365,
$ 0,31,60,91,121,152,182,213,244,274,305,335,366/
dtim=dble(tim)
if(dtim.lt.0.d0) go to 3
5 l=1
C Leap year test .....
if((iye/4)*4.eq.iye) l=2
dyl=dble(float(nd(13,l)))*86400.d0
if(dtim.lt.dyl) go to 4
iye=iye+1
dtim=dtim-dyl
go to 5
3 iye=iye-1
l=1
C Leap year test .....
if((iye/4)*4.eq.iye) l=2
dyl=dble(float(nd(13,l)))*86400.d0
dtim=dtim+dyl
if(dtim.lt.0.d0) go to 3
4 dmin=dint(dtim/60.d0)
sec=sngl(dtim-dmin*60.d0)
dihr=dint(dmin/60.d0)
min=idint(dmin-dihr*60.d0+0.5d0)
didy=dint(dihr/24.d0)
ihr=idint(dihr-didy*24.d0+0.5d0)
idy=idint(didy+1.5d0)
do 1 i=1,13
if(nd(i,l).ge.idy) go to 2
1 continue
2 mon=i-1
idy=idy-nd(mon,l)
tim=sngl(dtim)
return
end
CBack up
subroutine fft(fin,fout,ln,itype)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:02:09:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: FFT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Fast Fourier Transform *
C * *
C * FIN .... Array to be transformed *
C * FOUT ... Transformed array *
C * LN ..... LOG2 of no. of data points *
C * ITYPE .. +1 for forward transform *
C * -1 for reverse transform *
C * *
C * Adapted from : *
C * *
C * Gonzalez, R.C. and Wintz, P., 1977. *
C * Digital Image Procssing, Addison-Wesley Publishing Co. *
C * *
C ******************************************************************************
complex fin(*),fout(*),u,w,t
pi=3.14159265
n=2**ln
nv2=n/2
nm1=n-1
do i=1,n
fout(i)=fin(i)
end do
j=1
do i=1,nm1
if(i.ge.j) go to 1
t=fout(j)
fout(j)=fout(i)
fout(i)=t
1 k=nv2
2 if(k.ge.j) go to 3
j=j-k
k=k/2
go to 2
3 j=j+k
end do
do l=1,ln
le=2**l
le1=le/2
u=(1.0,0.0)
w=cmplx(cos(pi/le1),-sin(pi/le1)*float(itype))
do j=1,le1
do i=j,n,le
ip=i+le1
t=fout(ip)*u
fout(ip)=fout(i)-t
fout(i)=fout(i)+t
end do
u=u*w
end do
end do
if(itype.eq.1) then ! Forward transform
do i=1,n
fout(i)=fout(i)/float(n)
end do
endif
return
end
CBack up
subroutine fft2(fin2,fout2,lnx,lny,nxmax,itype,fin1,fout1)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:02:09:1986 *
C * *
C * REVISION : Special code for square arrays *
C * removed JRH:25:11:1987 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: FFT2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : 2-D fast Fourier transform *
C * *
C * FIN2 .... input 2-D COMPLEX array *
C * FOUT2 ... output 2-D COMPLEX array *
C * LNX ..... LOG2(no. values in X to be transformed) *
C * LNY ..... LOG2(no. values in Y to be transformed) *
C * NXMAX ... first dimension of FIN2, FOUT2 *
C * ITYPE ... 1 .... forward transform *
C * -1 ... reverse transform *
C * FIN1 .... COMPLEX 1-D work array of dimension larger *
C * than : 2**LNX AND 2**LNY *
C * FOUT1 ... " " " " " " " *
C * " " " " *
C * *
C ******************************************************************************
complex fin2(nxmax,*),fout2(nxmax,*)
complex fin1(*),fout1(*)
nx=2**lnx
ny=2**lny
do i=1,nx
do j=1,ny
fin1(j)=fin2(i,j)
end do
call fft(fin1,fout1,lny,itype)
do j=1,ny
fout2(i,j)=fout1(j)
end do
end do
do j=1,ny
do i=1,nx
fin1(i)=fout2(i,j)
end do
call fft(fin1,fout1,lnx,itype)
do i=1,nx
fout2(i,j)=fout1(i)
end do
end do
C Following code removed 25/11/87 :
C IF(LNX.EQ.LNY) THEN ! Special for square arrays
C DO I=1,NX
C DO J=1,NY
C IF(ITYPE.EQ.1) THEN ! Forward transform
C FOUT2(I,J)=FOUT2(I,J)*CMPLX(FLOAT(NX))
C ELSE ! Reverse transform
C FOUT2(I,J)=FOUT2(I,J)/CMPLX(FLOAT(NX))
C ENDIF
C END DO
C END DO
C ENDIF
return
end
CBack up
subroutine fit2d(xa,ya,za,ca,ipow,ndata,nterm,error,x,y,z,
$ imode,ixint,iyint,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILTY REF:JRH:15:04:1987 *
C * *
C * REVISION : Mod. to parameter JRH:07:07:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: FIT2D *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : 2-D polynomial fit subroutine. *
C * *
C * XA ..... Input array of X-coordinates *
C * YA ..... Input array of Y-coordinates *
C * ZA ..... Input array of Z-values *
C * CA ..... Intermediate array of polynomial constants *
C * IPOW ... Array of powers of X and Y (returned by FIT2D) *
C * NDATA .. Number of input data points *
C * NTERM .. Number of terms in polynomial *
C * ERROR .. RMS deviation between data points and polynomial *
C * X ...... Input X-coordinate for evaluation using *
C * polynomial *
C * Y ...... Input Y-coordinate for evaluation using *
C * polynomial *
C * Z ...... Output Z-value from evaluation using polynomial *
C * IMODE .. Mode of operation : *
C * 0 ..... Evaluate constants in polynomial *
C * 1 ..... Evaluate Z(X,Y) *
C * IXINT .. Level of integration to be applied in *
C * X-direction (0 : no integration, *
C * 1 : integrate, *
C * -1 : differentiate) *
C * IYINT .. Level of integration to be applied in *
C * Y-direction (as for IXINT) *
C * IFAIL .. 0 for successful execution, otherwise 1 *
C * *
C * Uses NAG. *
C * *
C ******************************************************************************
integer pndatmax ! @7/7/89
parameter (pndatmax=100) ! Max. no. of data points
real xa(*),ya(*),za(*),ca(*),x,y,z
integer ipow(2,*)
real a(pndatmax,pndatmax),term(pndatmax,pndatmax),zz(pndatmax)
dimension aa(pndatmax,pndatmax),wks1(pndatmax),wks2(pndatmax) ! For NAG
data ndatmax/pndatmax/
if(n.gt.ndatmax) then
ifail=1
return
endif
if(imode.eq.0) then
ip=0
i=0
2 do j=0,i
ip=ip+1
ipow(1,ip)=i-j
ipow(2,ip)=j
if(ip.eq.nterm) go to 1
end do
i=i+1
go to 2
1 continue
do i=1,ndata ! Set up array of dependant variable
do j=1,nterm
ipx=ipow(1,j)
ipy=ipow(2,j)
if(ipx.eq.0.) then
facx=1.
else
facx=xa(i)**ipx
endif
if(ipy.eq.0.) then
facy=1.
else
facy=ya(i)**ipy
endif
term(i,j)=facx*facy
end do
end do
if(ndata.eq.nterm) then ! Exact fit
call f04ate(term,ndatmax,za,nterm,ca,aa,ndatmax,wks1,wks2,
$ ifail)
if(ifail.ne.0) then
ifail=1
return
endif
error=0.
else ! Least squares fit
do j=1,nterm
zz(j)=0.
do i=1,ndata ! Form right-hand side
zz(j)=zz(j)+term(i,j)*za(i)
end do
do k=j,nterm ! Form left-hand square matrix
a(j,k)=0.
do i=1,ndata
a(j,k)=a(j,k)+term(i,j)*term(i,k)
if(j.ne.k) a(k,j)=a(j,k) ! Fill in other half of A
end do
end do
end do
call f04ase(a,ndatmax,zz,nterm,ca,wks1,wks2,ifail)
if(ifail.ne.0) then
ifail=1
return
endif
error=0.
do i=1,ndata
error=error+za(i)**2
end do
do j=1,nterm
error=error-ca(j)*zz(j)
end do
error=sqrt(error/float(ndata))
endif
else
if(imode.eq.1) then
z=0.
do j=1,nterm
ipx=ipow(1,j)+ixint
ipy=ipow(2,j)+iyint
if(ixint.eq.1) then ! Integrate
facx=1./float(ipx)
else
if(ixint.eq.0) then
facx=1.
else
if(ixint.eq.-1) then ! Differentiate
facx=float(ipx+1)
else
ifail=1
return
endif
endif
endif
if(iyint.eq.1) then ! Integrate
facy=1./float(ipy)
else
if(iyint.eq.0) then
facy=1.
else
if(iyint.eq.-1) then ! Differentiate
facy=float(ipy+1)
else
ifail=1
return
endif
endif
endif
if(facx.ne.0..and.ipx.ne.0.) then
facx=facx*(x**ipx)
endif
if(facy.ne.0..and.ipy.ne.0.) then
facy=facy*(y**ipy)
endif
z=z+ca(j)*facx*facy
end do
else
ifail=1
return
endif
endif
ifail=0
return
end
CBack up
real*4 function gasdev3(idum)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:31:03:1995 *
C * *
C * REVISION : ------------- JRH:--:--:1995 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: gasdev3 *
C * TYPE : real*4 function *
C * *
C * FUNCTION : Random number generator yielding normal distribution *
C * with zero mean and unit variance (based on Numerical *
C * Recipes). *
C * *
C ******************************************************************************
integer*4 idum
C
real*4 v1,v2,r,fac,gset
real*4 ran3
integer*4 iset
data iset/0/
save iset,gset
if (iset.eq.0) then
1 v1=2.*ran3(idum)-1.
v2=2.*ran3(idum)-1.
r=v1**2+v2**2
if(r.ge.1.)go to 1
fac=sqrt(-2.*log(r)/r)
gset=v1*fac
gasdev3=v2*fac
iset=1
else
gasdev3=gset
iset=0
endif
return
end
CBack up
subroutine geoginp(nin,nout,nlog,annot,ncoord,long,lat,x,y,atype)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:20:10:1988 *
C * *
C * REVISION : Mod. for null read after error to cope *
C * with Sun bug JRH:25:07:1989 *
C * Minor mods. to this header JRH:09:10:1989 *
C * Mod. to accept lower case "e" and "n" JRH:24:04:1991 *
C * Minor mod. to formats JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: GEOGINP *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Inputs geographical coordinates (grid (X,Y) or longitude *
C * and latitude (LONG,LAT)) and indicates which has been *
C * input. Each grid coordinate must be preceded by "E" or *
C * "N" (they may be in any order). LONG must be input before *
C * LAT. West longitude or South latitude is indicated by a *
C * negative sign prior any value in longitude or latitude. *
C * NCOORD is the number of coordinates to be entered (ie. 1 *
C * or 2). *
C * *
C * NIN ...... Operator input device *
C * NOUT ..... Operator output device *
C * NLOG ..... Log file device *
C * ANNOT .... Annotation (40 chars.) *
C * NCOORD ... Number of coordinates to be input : *
C * 1 : only LONG or LAT, or X or Y to be input - *
C * if LONG or LAT, result is output to both *
C * LONG and LAT, *
C * if X or Y, result is output to X or Y as *
C * appropriate. *
C * 2 : both LONG and LAT, or X and Y to be input. *
C * LONG ..... Longitude (radians internally, *
C * deg. *
C * or (deg., min.) *
C * or (deg., min., sec.) at input) *
C * LAT ..... Latitude (radians internally, *
C * deg. *
C * or (deg., min.) *
C * or (deg., min., sec.) at input) *
C * X ........ East grid coordinates *
C * Y ........ North grid coordinates *
C * ATYPE .... 'L' for input of longitude and latitude *
C * 'G' for input of grid coordinates, *
C * *
C ******************************************************************************
real*8 long,lat,x,y
character*40 annot
real*8 longd,longm,longs,latd,latm,lats,d2r,sig1,sig2
dimension ind(6)
character*1 atype
character*81 buff
d2r=datan(1.d0)/45.d0 ! Degrees to radians
buff(1:1)=' '
7 write(nout,1) annot
write(nlog,1) annot
1 format(1x,a40)
read(nin,2) buff(2:81)
2 format(a80)
write(nlog,3) buff(2:81)
3 format(1x,a80)
do i=2,81 ! @24/4/91
if(buff(i:i).eq.'e') buff(i:i)='E' ! @24/4/91
if(buff(i:i).eq.'n') buff(i:i)='N' ! @24/4/91
end do ! @24/4/91
ie=index(buff,'E')
in=index(buff,'N')
if(ie.eq.0.and.in.eq.0) then ! Longitude and latitude
nval=0
do i=1,80
ip1=i+1
if((buff(i:i).eq.' '.or.buff(i:i).eq.',').and. ! Data separator
$ (buff(ip1:ip1).ne.' '.and.buff(ip1:ip1).ne.',')) then ! Data
nval=nval+1 ! Number of values in buffer
if((ncoord.eq.1.and.nval.gt.3).or.
$ (ncoord.eq.2.and.nval.gt.6)) go to 7 ! Too many input values
ind(nval)=ip1 ! Pointer to first character of value
endif
end do
if(nval.eq.0) go to 7 ! No values input
if(ncoord.eq.1) then
if(index(buff,'-').eq.0) then ! No minus sign
sig1=1.d0
else ! Minus sign
sig1=-1.d0
endif
if(nval.eq.3) then
read(buff,*,err=4) longd,longm,longs ! NVAL = 3
C WRITE(NLOG,8) LONGD,LONGM,LONGS
long=dsign(dabs(longd)+dabs(longm)/60.d0
$ +dabs(longs)/3600.d0,
$ sig1)*d2r
else if(nval.eq.2) then ! NVAL = 2
read(buff,*,err=4) longd,longm
C WRITE(NLOG,9) LONGD,LONGM
long=dsign(dabs(longd)+dabs(longm)/60.d0,sig1)*d2r
else ! NVAL = 1
read(buff,*,err=4) longd
C WRITE(NLOG,10) LONGD
long=dsign(dabs(longd),sig1)*d2r
endif
lat=long
else ! Assume NCOORD = 2
if((nval/2)*2.ne.nval) go to 7 ! Not even number of values
idum=ind(nval/2+1) ! Pointer to first char. of first value of latitude
if(index(buff(2:idum-1),'-').eq.0) then ! No minus sign
sig1=1.d0
else ! Minus sign
sig1=-1.d0
endif
if(index(buff(idum:81),'-').eq.0) then ! No minus sign
sig2=1.d0
else ! Minus sign
sig2=-1.d0
endif
if(nval.eq.6) then ! NVAL = 6
read(buff,*,err=4) longd,longm,longs,latd,latm,lats
C WRITE(NLOG,8) LONGD,LONGM,LONGS,LATD,LATM,LATS
C 8 FORMAT(1X,F5.0,X,F3.0,X,F7.4,X,F5.0,X,F3.0,X,F7.4)
long=dsign(dabs(longd)+dabs(longm)/60.d0
$ +dabs(longs)/3600.d0,
$ sig1)*d2r
lat= dsign(dabs(latd) +dabs(latm)/60.d0
$ +dabs(lats)/3600.d0 ,
$ sig2)*d2r
else if(nval.eq.4) then ! NVAL = 4
read(buff,*,err=4) longd,longm,latd,latm
C WRITE(NLOG,9) LONGD,LONGM,LATD,LATM
C 9 FORMAT(X,F5.0,X,F7.4,X,F5.0,X,F7.4)
long=dsign(dabs(longd)+dabs(longm)/60.d0,sig1)*d2r
lat= dsign(dabs(latd) +dabs(latm)/60.d0 ,sig2)*d2r
else ! NVAL = 2
read(buff,*,err=4) longd,latd
C WRITE(NLOG,10) LONGD,LATD
C 10 FORMAT(1X,F11.6,X,F11.6)
long=dsign(dabs(longd),sig1)*d2r
lat= dsign(dabs(latd) ,sig2)*d2r
endif
endif
atype='L'
else if(ie.ne.0.or.in.ne.0) then ! Grid coordinates
if(ncoord.eq.1) then
if(ie.ne.0..and.in.ne.0) then ! Both grid coordinates - failure
go to 7
else
if(ie.ne.0) then ! Eastings
read(buff(ie+1:80),*,err=4) x
C WRITE(NLOG,15) ' E',X
C 15 FORMAT(A2,F12.2)
else ! Northings
read(buff(in+1:80),*,err=4) y
C WRITE(NLOG,15) ' N',Y
endif
endif
else ! Assume NCOORD = 2
if(ie.ne.0..and.in.ne.0) then ! Both grid coordinates
read(buff(ie+1:80),*,err=4) x
read(buff(in+1:80),*,err=4) y
C WRITE(NLOG,11) X,Y
C 11 FORMAT(' E',F12.2,' N',F12.2)
atype='G'
else ! Failure
go to 7
endif
endif
atype='G'
endif
return
4 read(buff,*) ! Null read to cope with Sun bug @25/7/89
go to 7 ! @25/7/89
end
CBack up
subroutine geoginpbuff(buffin,ncoord,long,lat,x,y,atype,
$ ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:10:10:1989 *
C * *
C * REVISION : Correction to name in header JRH:21:05:1990 *
C * Correction so that ifail is set to *
C * zero for successful execution JRH:28:09:1995 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: GEOGINPBUFF *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Inputs geographical coordinates (grid (X,Y) or longitude *
C * and latitude (LONG,LAT)) and indicates which has been *
C * input. Each grid coordinate must be preceded by "E" or *
C * "N" (they may be in any order). LONG must be input before *
C * LAT. West longitude or South latitude is indicated by a *
C * negative sign prior to any value in longitude or latitude.*
C * NCOORD is the number of coordinates to be entered (ie. 1 *
C * or 2). *
C * *
C * BUFFIN ... Input buffer *
C * NCOORD ... Number of coordinates to be input : *
C * 1 : only LONG or LAT, or X or Y to be input - *
C * if LONG or LAT, result is output to both *
C * LONG and LAT, *
C * if X or Y, result is output to X or Y as *
C * appropriate. *
C * 2 : both LONG and LAT, or X and Y to be input. *
C * LONG ..... Longitude (radians internally, *
C * deg. *
C * or (deg., min.) *
C * or (deg., min., sec.) at input) *
C * LAT ..... Latitude (radians internally, *
C * deg. *
C * or (deg., min.) *
C * or (deg., min., sec.) at input) *
C * X ........ East grid coordinates *
C * Y ........ North grid coordinates *
C * ATYPE .... 'L' for input of longitude and latitude *
C * 'G' for input of grid coordinates, *
C * IFAIL .... returned 0 if execution satisfactory, *
C * otherwise 1 *
C * *
C * NOTE : This is a special version of subroutine GEOGINP, *
C * for input from a buffer *
C * *
C ******************************************************************************
real*8 long,lat,x,y
character*80 buffin
real*8 longd,longm,longs,latd,latm,lats,d2r,sig1,sig2
dimension ind(6)
character*1 atype
character*81 buff
d2r=datan(1.d0)/45.d0 ! Degrees to radians
buff(1:1)=' '
buff(2:81)=buffin(1:80)
ie=index(buff,'E')
in=index(buff,'N')
if(ie.eq.0.and.in.eq.0) then ! Longitude and latitude
nval=0
do i=1,80
ip1=i+1
if((buff(i:i).eq.' '.or.buff(i:i).eq.',').and. ! Data separator
$ (buff(ip1:ip1).ne.' '.and.buff(ip1:ip1).ne.',')) then ! Data
nval=nval+1 ! Number of values in buffer
if((ncoord.eq.1.and.nval.gt.3).or.
$ (ncoord.eq.2.and.nval.gt.6)) then ! Too many input values
ifail=1
return
endif
ind(nval)=ip1 ! Pointer to first character of value
endif
end do
if(nval.eq.0) then ! No values input
ifail=1
return
endif
if(ncoord.eq.1) then
if(index(buff,'-').eq.0) then ! No minus sign
sig1=1.d0
else ! Minus sign
sig1=-1.d0
endif
if(nval.eq.3) then
read(buff,*,err=4) longd,longm,longs ! NVAL = 3
long=dsign(dabs(longd)+dabs(longm)/60.d0
$ +dabs(longs)/3600.d0,
$ sig1)*d2r
else if(nval.eq.2) then ! NVAL = 2
read(buff,*,err=4) longd,longm
long=dsign(dabs(longd)+dabs(longm)/60.d0,sig1)*d2r
else ! NVAL = 1
read(buff,*,err=4) longd
long=dsign(dabs(longd),sig1)*d2r
endif
lat=long
else ! Assume NCOORD = 2
if((nval/2)*2.ne.nval) then ! Not even number of values
ifail=1
return
endif
idum=ind(nval/2+1) ! Pointer to first char. of first value of latitude
if(index(buff(2:idum-1),'-').eq.0) then ! No minus sign
sig1=1.d0
else ! Minus sign
sig1=-1.d0
endif
if(index(buff(idum:81),'-').eq.0) then ! No minus sign
sig2=1.d0
else ! Minus sign
sig2=-1.d0
endif
if(nval.eq.6) then ! NVAL = 6
read(buff,*,err=4) longd,longm,longs,latd,latm,lats
long=dsign(dabs(longd)+dabs(longm)/60.d0
$ +dabs(longs)/3600.d0,
$ sig1)*d2r
lat= dsign(dabs(latd) +dabs(latm)/60.d0
$ +dabs(lats)/3600.d0 ,
$ sig2)*d2r
else if(nval.eq.4) then ! NVAL = 4
read(buff,*,err=4) longd,longm,latd,latm
long=dsign(dabs(longd)+dabs(longm)/60.d0,sig1)*d2r
lat= dsign(dabs(latd) +dabs(latm)/60.d0 ,sig2)*d2r
else ! NVAL = 2
read(buff,*,err=4) longd,latd
long=dsign(dabs(longd),sig1)*d2r
lat= dsign(dabs(latd) ,sig2)*d2r
endif
endif
atype='L'
else if(ie.ne.0.or.in.ne.0) then ! Grid coordinates
if(ncoord.eq.1) then
if(ie.ne.0..and.in.ne.0) then ! Both grid coordinates - failure
ifail=1
return
else
if(ie.ne.0) then ! Eastings
read(buff(ie+1:80),*,err=4) x
else ! Northings
read(buff(in+1:80),*,err=4) y
endif
endif
else ! Assume NCOORD = 2
if(ie.ne.0..and.in.ne.0) then ! Both grid coordinates
read(buff(ie+1:80),*,err=4) x
read(buff(in+1:80),*,err=4) y
atype='G'
else ! Failure
ifail=1
return
endif
endif
atype='G'
endif
ifail=0 ! @28/9/95
return
4 ifail=1
return
end
CBack up
subroutine gsnumber(x,y,wid,r,ang,idplace)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:14:07:1988 *
C * *
C * REVISION : Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: GSNUMBER *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Version of NUMBER for GS software. *
C * *
C * Labelled COMMON GSCOM must be included in main program, *
C * in which LUPLT, LULOG, NSET, NCHARP and NDPLP must be set.*
C * *
C * X,Y ...... coordinates of bottom LH corner of text *
C * WID ...... character width *
C * R ........ REAL number to be written *
C * ANG ...... angle (deg.) of text anticlockwise from X-axis *
C * IDPLACE ... no. of decimal places to be plotted *
C * *
C ******************************************************************************
character*10 form
character*80 a
common/gscom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ nset, ! Symbol set no.
$ nopen, ! Current pen number
$ ncharp, ! No. of chars. in coordinate output
$ ndplp ! No. of decimal places in coordinate output
if(idplace.ge.0.) then ! Output as REAL no.
call nchar(r,idplace,n,nasc) ! Tot. chars. (incl. sign and dec. point)
write(form,4) nasc,idplace
4 format('(F',i2,'.',i2,')')
write(a,form) r
else ! Output as INTEGER
call nchar(r,0,n,ndum) ! Tot. chars. (incl. sign and dec. point)
nasc=n+1 ! Total characters (incl. sign)
write(form,5) nasc
5 format('(1X,I',i2,')')
write(a,form) nint(r)
endif
call gssymbol(x,y,wid,a,ang,nasc)
return
end
CBack up
subroutine gssymbol(x,y,wid,a,ang,nasc)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:14:07:1988 *
C * *
C * REVISION : Reads font on 99 rather than 7 JRH:27:11:1988 *
C * Reads font on 101 rather than 99 JRH:16:03:1989 *
C * Mod. to parameter JRH:07:07:1989 *
C * Mods. where byte variable used as *
C * array index JRH:07:07:1989 *
C * Mod. for use on Sun JRH:07:11:1990 *
C * Minor mod. to format JRH:16:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: GSSYMBOL *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Version of SYMBOL for GS software. *
C * *
C * Labelled COMMON GSCOM must be included in main program, *
C * in which LUPLT, LULOG, NSET, NCHARP and NDPLP must be set.*
C * *
C * X,Y ...... coordinates of bottom LH corner of text *
C * WID ...... character width *
C * A ........ text to be written *
C * ANG ...... angle (deg.) of text anticlockwise from X-axis *
C * NASC ..... no. of ASCII characters *
C * *
C ******************************************************************************
integer pbytemax ! @7/7/89
parameter(pbytemax=20000) ! Max. no. of bytes in file
character*(*) a
dimension xa(255),ya(255)
dimension idir(0:255)
byte byt(pbytemax)
character*1 abyt(pbytemax),asc
character*1 aset ! &7/11/90
character*512 adum ! &7/11/90
character*18 form
common/gscom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ nset, ! Symbol set no.
$ nopen, ! Current pen number
$ ncharp, ! No. of chars. in coordinate output
$ ndplp ! No. of decimal places in coordinate output
equivalence (byt,abyt)
save nsetlst,idir,byt
data nsetlst/0/
dtr=3.14159265/180.
write(form,8) ncharp,ndplp,ncharp,ndplp
8 format('(I3/(F',i1,'.',i1,',X,F',i1,'.',i1,'))')
if(nsetlst.ne.nset) then ! Load symbol set NSET
write(aset,7) nset
7 format(i1) ! &7/11/90
open(unit=101,file='set'//aset//'.bin',status='UNKNOWN', ! &7/11/90
$ form='UNFORMATTED',recl=512,access='DIRECT') ! &7/11/90
do i=0,255
idir(i)=0 ! Clear directory
end do
irec=0
2 irec=irec+1
read(101,rec=irec,iostat=ios) adum(1:512) ! &7/11/90
if(ios.ne.0) go to 1 ! End of file &7/11/90
do i=1,510
ind=(irec-1)*510+i
if(ind.gt.pbytemax) then
write(6,6)
write(lulog,6)
6 format(' Array DIMENSIONs too small',
$ ' ..... program terminated')
stop
endif
abyt(ind)=adum(i:i)
end do
go to 2
1 ipoint=4 ! Pointer to BYT array
5 if(byt(ipoint+1).ne.32) then ! Space signifies end of data
index=byt(ipoint) ! @7/7/89
idir(index)=ipoint ! Set up directory &7/7/89
idir(index+128)=ipoint ! Set up for either parity &7/7/89
ipoint=ipoint+3
nppoint=byt(ipoint) ! No. of plotting points
ipoint=ipoint+nppoint*3+1
if(ipoint.lt.(irec-1)*510) go to 5 ! Not yet at end of data
endif
close(unit=101) ! &16/3/89
endif
offset=0.
do iasc=1,nasc ! Step through characters
asc=a(iasc:iasc)
ibyt=ichar(asc) ! Byte value
ipoint=idir(ibyt) ! Pointer to character in BYT
if(ipoint.eq.0) then ! Character not in font
ibyt=32 ! Force a space
ipoint=idir(ibyt)
endif
iwid=byt(ipoint+2)
npt=byt(ipoint+3)
j=0
do i=1,npt ! Step through plotting points
xx= byt(ipoint+(i-1)*3+4)
yy= byt(ipoint+(i-1)*3+5)
ipen=byt(ipoint+(i-1)*3+6)
if(ipen.eq.3) then ! IPEN = 3
if(i.ne.1) then ! Output last segment
do k=1,j
call xy2ra(xa(k)+offset,ya(k),r,az,
$ .true.,.true.) ! Convert to polars
az=az-ang*dtr ! Rotate by ANG anticlockwise
r=r*wid/20. ! Scale (nominal width of 20 units)
call ra2xy(r,az,xa(k),ya(k)) ! Convert to rectangulars
xa(k)=xa(k)+x
ya(k)=ya(k)+y
end do
if(j.eq.2) then ! GS software will not allow 2 point only
j=3
xa(3)=xa(2)
ya(3)=ya(2)
endif
write(luplt,form) j,(xa(k),ya(k),k=1,j) ! Output last segment
endif
j=1
xa(1)=xx
ya(1)=yy
else ! IPEN = 2
j=j+1
xa(j)=xx
ya(j)=yy
endif
end do
if(j.ne.0.) then
do k=1,j
call xy2ra(xa(k)+offset,ya(k),r,az,
$ .true.,.true.) ! Convert to polars
az=az-ang*dtr ! Rotate by ANG anticlockwise
r=r*wid/20. ! Scale (nominal width of 20 units)
call ra2xy(r,az,xa(k),ya(k)) ! Convert to rectangulars
xa(k)=xa(k)+x
ya(k)=ya(k)+y
end do
if(j.eq.2) then ! GS software will not allow 2 point only
j=3
xa(3)=xa(2)
ya(3)=ya(2)
endif
write(luplt,form) j,(xa(k),ya(k),k=1,j) ! Output last segment
endif
offset=offset+iwid
end do
return
end
CBack up
subroutine heapsort(ain,aout,ind,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:29:08:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB .FVS *
C * ROUTINE NAME: HEAPSORT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Heapsort based on Numerical Recipes subroutine SORT2 *
C * *
C * AIN ..... Array to be sorted (unchanged on output) *
C * AOUT .... Resultant sorted array (increasing) *
C * IND ..... Array of sorted indices *
C * N ....... Size of arrays AIN and AOUT *
C * *
C ******************************************************************************
real*4 ain(*),aout(*)
integer*4 ind(*)
integer*4 n
real*4 rra
integer*4 i,j,l,ir,idum
l=n/2+1
ir=n
do i=1,n
aout(i)=ain(i) ! Copy input array to output array
ind(i)=i ! Generate initial idum array
end do
if(n.eq.1) return ! Special for only one record
10 continue
if(l.gt.1)then
l=l-1
rra=aout(l)
idum=ind(l)
else
rra=aout(ir)
idum=ind(ir)
aout(ir)=aout(1)
ind(ir)=ind(1)
ir=ir-1
if(ir.eq.1)then
aout(1)=rra
ind(1)=idum
return
endif
endif
i=l
j=l+l
20 if(j.le.ir)then
if(j.lt.ir)then
if(aout(j).lt.aout(j+1))j=j+1
endif
if(rra.lt.aout(j))then
aout(i)=aout(j)
ind(i)=ind(j)
i=j
j=j+j
else
j=ir+1
endif
go to 20
endif
aout(i)=rra
ind(i)=idum
go to 10
end
CBack up
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
CBack up
integer*4 function idplace(r,isig)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:10:08:1989 *
C * *
C * REVISION : Mod. using INT2 so that values like *
C * 0.1 are handled correctly JRH:18:11:1991 *
C * Full declarations added JRH:18:11:1991 *
C * Addition of external JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: IDPLACE *
C * TYPE : INTEGER*4 FUNCTION *
C * *
C * FUNCTION : Returns required number of decimal places to output a *
C * REAL number, R, to ISIG significant figures *
C * *
C ******************************************************************************
real*4 r ! @18/11/91
integer*4 isig ! @18/11/91
integer*4 nleft,int2 ! @18/11/91
external int2
if(r.ne.0.) then
nleft=int2(alog10(abs(r))+1) ! No. of figs. to L of dec. pt. @18/11/91
else ! r = 0
nleft=1
endif
idplace=isig-nleft ! ISIG significant figs
if(idplace.lt.0) idplace=0
return
end
CBack up
subroutine inbuf(type,key,buf,nbuf,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:15:09:1993 *
C * *
C * REVISION : ------------- JRH:--:--:1993 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: INBUF *
C * TYPE : MAIN *
C * *
C * FUNCTION : Inputs buffer from device 5, based on keyword. *
C * *
C * type ..... '+' for inclusion of matched record, *
C * '-' for exclusion of matched record *
C * key ...... Keyword *
C * buf ...... Buffer selected *
C * nbuf ..... Number of characters in buf *
C * ifail .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 nbuf,ifail
character*1 type
character*(*) key
character*(*) buf
C
integer*4 nkey,nbufin,lenchar,ios
character*7 fmt
C
nkey=len(key)
nbufin=len(buf)
write(fmt,1) nbufin
1 format('(a',i4,')')
2 read(5,fmt,iostat=ios) buf(1:nbufin)
if(ios.ne.0) then ! Error or end of file
ifail=1
return
endif
if((type.eq.'+'.and. ! Inclusion required
$ (buf(1:nkey).eq.key(1:nkey).and. ! Found a match
$ buf(nkey+1:nkey+1).eq.' ')).or. ! followed by a space
$ (type.eq.'-'.and..not. ! Exclusion required
$ (buf(1:nkey).eq.key(1:nkey).and. ! Found a match
$ buf(nkey+1:nkey+1).eq.' '))) then ! followed by a space
nbuf=lenchar(buf)
ifail=0
return
endif
go to 2 ! Try another record
end
CBack up
subroutine inilink(prev,next,stack,ifirst,ilast,topstack,maxlink)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:11:10:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: INILINK *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Initialises link list (to empty) *
C * *
C * PREV ..... array of pointers to previous item in list *
C * NEXT ..... array of pointers to next item in list *
C * STACK .... stack of unused list pointers *
C * IFIRST ... pointer to first item in list *
C * ILAST .... pointer to last item in list *
C * TOPSTACK . pointer to top of stack *
C * MAXLINK .. maximum number of items in list *
C * *
C ******************************************************************************
integer*4 prev(maxlink),next(maxlink),stack(maxlink)
integer*4 ifirst,ilast,topstack,maxlink
integer*4 i
do i=1,maxlink ! Clear arrays
prev(i)=0
next(i)=0
end do
ifirst=0 ! Empty list
ilast=0
topstack=maxlink
do i=1,maxlink
stack(i)=maxlink+1-i ! "1" is on top of stack
end do
return
end
CBack up
subroutine inslink(prev,next,stack,ifirst,ilast,topstack,maxlink,
$ ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:11:10:1991 *
C * *
C * REVISION : Correction to error return JRH:23:03:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: INSLINK *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Inserts item in link list (at position ILAST) *
C * *
C * PREV ..... array of pointers to previous item in list *
C * NEXT ..... array of pointers to next item in list *
C * STACK .... stack of unused list pointers *
C * IFIRST ... pointer to first item in list *
C * ILAST .... pointer to last item in list *
C * TOPSTACK . pointer to top of stack *
C * MAXLINK .. maximum number of items in list *
C * IFAIL .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 prev(maxlink),next(maxlink),stack(maxlink)
integer*4 ifirst,ilast,topstack,maxlink,ifail
integer*4 new
if(topstack.eq.0) then ! No more room
ifail=1
return
endif
new=stack(topstack) ! Take new pointer off stack
if(new.eq.0) then ! Invalid pointer
ifail=1 ! &23/3/93
return
endif
stack(topstack)=0 ! Clear pointer
topstack=topstack-1
if(ilast.eq.0) then ! Empty list
ifirst=new
else
next(ilast)=new
endif
prev(new)=ilast
next(new)=0
ilast=new
ifail=0
return
end
CBack up
subroutine inslink2(iins,new,
$ prev,next,stack,ifirst,ilast,topstack,maxlink,
$ ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:05:08:1993 *
C * *
C * REVISION : -------------------------- JRH:--:--:1993 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: INSLINK2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Inserts item in link list (at position IINS) *
C * *
C * IINS ..... pointer after which item to be inserted *
C * (if IINS = 0, item is to be inserted at start *
C * of link list) *
C * NEW ...... pointer for new item *
C * PREV ..... array of pointers to previous item in list *
C * NEXT ..... array of pointers to next item in list *
C * STACK .... stack of unused list pointers *
C * IFIRST ... pointer to first item in list *
C * ILAST .... pointer to last item in list *
C * TOPSTACK . pointer to top of stack *
C * MAXLINK .. maximum number of items in list *
C * IFAIL .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 prev(maxlink),next(maxlink),stack(maxlink)
integer*4 iins,new,ifirst,ilast,topstack,maxlink,ifail
integer*4 idum
if(topstack.eq.0) then ! No more room
ifail=1
return
endif
new=stack(topstack) ! Take new pointer off stack
if(new.eq.0) then ! Invalid pointer
ifail=1
return
endif
stack(topstack)=0 ! Clear pointer
topstack=topstack-1
if((iins.ne.0.and.ilast.eq.0)) then ! IINS is invalid pointer
ifail=1
return
endif
if(ilast.eq.0) then ! Empty list
ifirst=new
ilast=new
next(new)=0
prev(new)=0
else if(iins.eq.0) then ! Insert item at start of link list
next(new)=ifirst
prev(new)=0
prev(ifirst)=new
ifirst=new
else ! IINS cannot be zero
idum=next(iins)
if(idum.eq.0.and.iins.ne.ilast) then ! IINS is invalid pointer
ifail=1
return
endif
next(iins)=new
next(new)=idum
prev(new)=iins
if(idum.eq.0) then ! End of link list
ilast=new
else
prev(idum)=new
endif
endif
ifail=0
return
end
CBack up
subroutine intpoint(xin,xout,n,m)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:29:12:1993 *
C * *
C * REVISION : Mod. for one- and two-point segments JRH:17:08:1994 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: INTPOINT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Subroutine to interpolate point data *
C * (simplified real*8 version of div1d in forlib). *
C * *
C * xin ..... array of n input values *
C * xout .... array of nm output interpolated values *
C * n ....... number of input values *
C * m ....... interpolation factor *
C * *
C ******************************************************************************
real*8 xin(*),xout(*)
integer*4 n,m
C
real*8 d1l,d2l,d1r,d2r,xl,xr,c1l,c2l,c1r,c2r
integer*4 i,ioff,ii
C
xout(1)=xin(1)
if(n.eq.1) then ! Only one point @17/8/94
return ! @17/8/94
else if(n.eq.2) then ! Two points @17/8/94
d1r=(xin(2)-xin(1))*2.d0 ! @17/8/94
d2r=0.d0 ! @17/8/94
endif ! @17/8/94
do i=1,n-1 ! Do "parabolic" interpolation
if(i.ne.1) then
d1l=d1r ! Use "old" RH derivative
d2l=d2r ! " " " "
endif
if(i.ne.n-1) then ! Not at RH end
d1r=xin(i+2)-xin(i) ! 1st derivative, RH curve
d2r=xin(i+2)+xin(i)-2*xin(i+1) ! 2nd " " "
endif
ioff=(i-1)*m+1
do ii=1,m ! Step within an interpolation interval
xl=dble(ii)/dble(m)
xr=xl-1.d0
c1l=xl/2.d0
c2l=xl*xl/2.d0
c1r=xr/2.d0
c2r=xr*xr/2.d0
if(i.eq.1) then ! LH end
xout(ii+ioff)=xin(i+1)+d1r*c1r+d2r*c2r ! RH parabola
else
if(i.eq.n-1) then ! RH end
xout(ii+ioff)=xin(i)+d1l*c1l+d2l*c2l ! LH parabola
else ! Other
xout(ii+ioff)=(xin(i) +d1l*c1l+d2l*c2l)*(-xr) ! LH parabola
$ +(xin(i+1)+d1r*c1r+d2r*c2r)*xl ! RH parabola
endif
endif
end do
end do
C
return
end
CBack up
function int2(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: INT2 *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : As INT, but trunctates towards -(infinity) *
C * *
C ******************************************************************************
real*4 x
int2=int(x)
if(float(int2).eq.x) return
if(x.lt.0.) int2=int2-1
return
end
CBack up
subroutine invmerc(x,y,long,lat,ae,ecs,lat0,s0,longor,lator,
$ eps,maxit,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:14:10:1988 *
C * *
C * REVISION : eps declared as real*8 JRH:09:06:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: INVMERC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts Mercator grid coordinates (X,Y) to longitude and *
C * latitude (LONG,LAT). LAT is computed by Newton's method. *
C * *
C * X ........ East grid coordinates *
C * Y ........ North grid coordinates *
C * LONG ..... Longitude (radians) *
C * LAT ..... Latitude (radians) *
C * AE ....... Semi-major axis of Earth *
C * ECS ...... Eccentricity squared of figure of Earth *
C * LAT0 ..... Latitude at which scale of projection is *
C * defined *
C * S0 ....... Scale of projection at LAT0 *
C * LONGOR ... Longitude of origin of projection *
C * LATOR .... Latitude of origin of projection *
C * EPS ...... Tolerance in Y *
C * MAXIT .... Max. no. of iterations *
C * IFAIL .... returned 0 if convergence satisfactory, *
C * otherwise 1 *
C * *
C ******************************************************************************
real*8 long,lat,x,y,ae,ecs,lat0,s0,longor,lator
real*8 alpha,e,slor,eslor,yor,yy,ee,sl,esl,diff,eps
alpha=ae*s0/dsqrt(1.d0+(1.d0-ecs)*(dtan(lat0)**2))
e=dsqrt(ecs)
slor=dsin(lator)
eslor=e*slor
yor=alpha*(dlog((1.d0+slor)
$ /(1.d0-slor))
$ -e*dlog((1.d0+eslor)
$ /(1.d0-eslor)))/2.d0
long=x/alpha+longor
yy=y+yor ! Compute initial approximation to LAT
ee=dexp(2.d0*yy/alpha)
lat=dasin((ee-1.d0)/(ee+1.d0))
do i=1,maxit
sl=dsin(lat)
esl=e*sl
diff=alpha*(dlog((1.d0+sl)
$ /(1.d0-sl))
$ -e*dlog((1.d0+esl)
$ /(1.d0-esl)))/2.-yy
if(dabs(diff).lt.eps) then
ifail=0
return
endif
lat=lat-diff*(1.d0-ecs*(dsin(lat)**2))*cos(lat)
$ /(alpha*(1.d0-ecs))
end do
ifail=1
return
end
CBack up
double precision function invrob(
$ x1,y1,z1,x2,y2,z2,saz,caz,a,ecs,eps)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : GEODESY JRH:18:02:1986 *
C * *
C * REVISION : Mod. to replace trig. functions for *
C * greater accuracy for small ranges JRH:18:10:1988 *
C * Removal of redundant variables JRH:27:11:1996 *
C * Quickie for horizontal coincidence JRH:12:11:1998 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: INVROB *
C * TYPE : DOUBLE PRECISION FUNCTION *
C * *
C * FUNCTION : Inverse Robbins computation for distance, and azimuth on *
C * spheroid. *
C * *
C * (X1,Y1,Z1) ..... Long., lat., height of point (1) *
C * (X2,Y2,Z2) ..... Long., lat., height of point (2) *
C * SAZ ............ Sin of azimuth (1) to (2) *
C * CAZ ............ Cos of azimuth (1) to (2) *
C * A .............. Semi-major axis of Earth *
C * ECS ............ First eccentricity squared of Earth *
C * EPS ............ Second eccentricity squared of Earth *
C * *
C * No. of calls of various operators : *
C * *
C * Operator No. of calls *
C * *
C * + or - 25 *
C * * or / 49 *
C * ** 1 *
C * SQRT 5 *
C * TRIGS 7 *
C * *
C * Checked against example in "The Austalian Map Grid" *
C * (National Mapping Council of Australia, Special *
C * Publication 7), page 10, 18 Feb. 1986 - accurate to *
C * +/- .001 metre, +/- .01 sec.. *
C * *
C ******************************************************************************
real z1,z2
double precision x1,y1,x2,y2,saz,caz,a,ecs,eps
double precision eabs,eprop,sp1,sp2,cp1,cdl,sdl,t,tau1
double precision chi,sig,f1,f2,dum,dum1,dum2,dum3,dums,h2,gh,g2
double precision sigp ! &27/11/96
C
data eabs,eprop/1.d-4,1.d-9/ ! Absolute and proportional errors
C
if(x1.eq.x2.and.y1.eq.y2) then ! Horizontal coincidence @12/11/98
saz=0.d0 ! Default to zero azimuth @12/11/98
caz=1.d0 ! Default to zero azimuth @12/11/98
invrob=dble(abs(z1-z2)) ! @12/11/98
return ! @12/11/98
endif ! @12/11/98
C
sp1=dsin(y1)
sp2=dsin(y2)
f1=a/dsqrt(1.d0-ecs*sp1*sp1)
f2=a/dsqrt(1.d0-ecs*sp2*sp2)
C CP1=DSQRT(1.D0-SP1*SP1) ! Faster than DCOS &18/10/88
cp1=dcos(y1) ! @18/10/88
cdl=dcos(x2-x1)
C SDL=DSIGN(DSQRT(1.D0-CDL*CDL),X2-X1) ! Faster than DSIN &18/10/88
sdl=dsin(x2-x1) ! @18/10/88
C T=((1.D0-ECS)*SP2+ECS*SP1*F1/F2)/DSQRT(1.D0-SP2*SP2) ! Faster &18/10/88
t=((1.d0-ecs)*sp2+ecs*sp1*f1/f2)/dcos(y2) ! @18/10/88
tau1=cp1*t-sp1*cdl
chi=dsqrt(tau1*tau1+sdl*sdl)
sig=dasin(chi/dsqrt(1.d0+t*t))
saz=sdl/chi
caz=tau1/chi
dum=cp1*caz
h2=dum*dum*eps
dum1=f1*sig ! Multiplier for brackets
dums=dum1 ! Sum of terms
sigp=sig*sig ! Sigma to required power
dum3=-sigp*h2*(1.d0-h2)*0.16666666666666667 ! One term in brackets
dum2=dum1*dum3 ! One term
dums=dums+dum2 ! Sum of terms
if(dabs(dum2).gt.eabs.and.dabs(dum3).gt.eprop) then ! Error too large
sigp=sigp*sig ! Increase power of sigma by one
gh=sp1*dum*eps ! We now need G*H
dum3=sigp*gh*(1.d0-2.d0*h2)*0.125d0 ! One term in brackets
dum2=dum1*dum3 ! One term
dums=dums+dum2 ! Sum of terms
if(dabs(dum2).gt.eabs.and.dabs(dum3).gt.eprop) then !Error too large
sigp=sigp*sig ! Increase power of sigma by one
g2=sp1*sp1*eps ! We now need G squared
dum3=sigp*(h2*(4.d0-7.d0*h2)-3.d0*g2*(1.d0-7.d0*h2))
$ *8.3333333333333333d-3
dum2=dum1*dum3 ! One term
dums=dums+dum2 ! Sum of terms
if(dabs(dum2).gt.eabs.and.dabs(dum3).gt.eprop) then ! Error too large
sigp=sigp*sig ! Increase power of sigma by one
dum3=-sigp*gh*2.0833333333333333d-2 ! One term in brackets
dum2=dum1*dum3 ! One term
dums=dums+dum2 ! Sum of terms
endif
endif
endif
C Correction for heights above sea level .....
dums=dums*(1.d0+dble(z1+z2)/(f1+f2)) ! Mean height
invrob=dsqrt(dums*dums+dble(z1-z2)**2) ! Slant
return
end
CBack up
subroutine iosproc(nout,nlog,ios,lfatal,lerr,leof)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:28:08:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: IOSPROC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Processes IOSTAT from READ *
C * *
C * NOUT ..... Operator output device number (6) *
C * NLOG ..... Log file device number *
C * IOS ...... Value returned by "IOSTAT=" in READ *
C * LFATAL ... .TRUE. if termination required on error or end *
C * of file, otherwise .FALSE. *
C * LERR ..... (If LFATAL=.FALSE.) .TRUE. if read error, *
C * otherwise .FALSE. *
C * LEOF ..... (If LFATAL=.FALSE.) .TRUE. if end of file, *
C * otherwise .FALSE. *
C * *
C ******************************************************************************
integer*4 nout,nlog,ios
logical lfatal,lerr,leof
if(ios.gt.0) then ! Read error
if(lfatal) then
write(nout,1)
write(nlog,1)
1 format(/' Read error ..... program terminated'/)
stop
else
lerr=.true.
leof=.false.
endif
else if(ios.eq.-1) then ! End of file
if(lfatal) then
write(nout,2)
write(nlog,2)
2 format(/' End of file ..... program terminated'/)
stop
else
lerr=.false.
leof=.true.
endif
else
lerr=.false.
leof=.false.
endif
return
end
CBack up
subroutine ipdig(ilp,i,j,ic)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:03:12:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: IPDIG *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Inputs digit, IC, at (I,J) to array ILP *
C * *
C ******************************************************************************
integer ilp(33,100)
i1=(i-1)/4+1
i2=5-i+(i1-1)*4
i3=ilp(i1,j)
ip=i2-1
i4=i3/(10**ip)
i5=i3/(10**i2)
ii=i4-i5*10
i3=i3+(ic-ii)*(10**ip)
ilp(i1,j)=i3
return
end
CBack up
subroutine iread(nin,nout,nlog,anot,acc,iarray,istart,istop,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:27:08:1986 *
C * *
C * REVISION : Error trap included JRH:04:09:1986 *
C * Mod. for output of message only if *
C * ISTART=0 or ISTOP=0 JRH:12:09:1986 *
C * Minor mods. JRH:31:12:1986 *
C * Minor mod. JRH:05:01:1987 *
C * Minor mod. JRH:17:05:1988 *
C * Mod. so that FORM is input without *
C * brackets JRH:21:06:1988 *
C * Mod. for null read after error to cope *
C * with Sun bug JRH:25:07:1989 *
C * Mod. for compilation on transputer JRH:07:09:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: IREAD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads in a set of INTEGER numbers *
C * *
C * NIN ..... Operator input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Annotation (40 chars.) *
C * ACC ..... Carriage-control character *
C * (' ' for new line, '$' for same line) *
C * IARRAY .. Resultant INTEGER array *
C * ISTART .. Start index of IARRAY *
C * ISTOP ... Stop index of IARRAY *
C * FORM .... Format for output to log device (40 chars., *
C * without brackets) *
C * *
C ******************************************************************************
integer iarray(*)
character*40 anot,form
character*1 acc
character*42 formt ! @7/9/89
logical cont ! @31/12/86
save cont ! @31/12/86
character*12 formout ! &5/11/86
data cont/.false./ ! @11/5/88
2 if(cont) then ! Continuation of terminal text @31/12/86 &25/7/89
write(formout,1) ' ',acc ! &31/12/86
1 format('(',a1,'X,A40,A1',a1,')') ! &5/11/87
else ! Normal text @31/12/86
write(formout,1) '/',acc ! @31/12/86
endif ! @31/12/86
write(nout,formout) anot,' ' ! &5/11/87
write(nlog,formout) anot,' ' ! &5/11/87
if(istart.ne.0.and.istop.ne.0) then ! &31/12/86
read(nin,*,err=3) (iarray(i),i=istart,istop) ! &4/9/86
formt='('//form//')' ! @7/9/89
write(nlog,formt) (iarray(i),i=istart,istop) ! &21/6/88 &7/9/89
cont=.false. ! Terminal text not to be continued @31/12/86
else ! @31/12/86
cont=.true. ! Terminal text to be continued @31/12/86
endif ! @12/9/86
return
3 read(nin,*) ! Null read to cope with Sun bug @25/7/89
go to 2 ! @25/7/89
end
CBack up
subroutine ireadfil(nin,nout,nlog,anot,nchar,ivar,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:06:01:1987 *
C * *
C * REVISION : Mod. for keyword checking JRH:19:01:1987 *
C * Simplification JRH:17:06:1988 *
C * Mod. for compilation on transputer JRH:07:09:1989 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: IREADFIL *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads INTEGER variable from file based on keyword. *
C * *
C * NIN ..... File input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Keyword in file *
C * NCHAR ... Number of characters in ANOT (max. 40) *
C * IVAR .... Resultant INTEGER variable *
C * FORM .... Format for output (40 chars., without brackets) *
C * *
C ******************************************************************************
character*(*) anot
character*40 form
character*80 buff
character*40 blank ! &17/6/88
character*40 formt ! @7/9/89
character*51 formt2
character*14 formt3
data blank/' '/
write(formt2,6) nchar,form
6 format('(1X,A',i2,',A3,',a40,')')
write(formt3,5) nchar
5 format('(/A9,A',i2,',A34/)')
rewind(nin)
4 read(nin,2,end=3) buff
2 format(a80)
if(buff(1:nchar).eq.anot.and. ! &17/6/88
$ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found &19/1/87
buff(1:nchar)=blank(1:nchar) ! Delete keyword
idum1=lhschar(buff) ! Index of first character of value @7/9/89
idum2=lenchar(buff) ! Index of last character of value @7/9/89
write(formt,7) idum2-idum1+1 ! @7/9/89
7 format('(i',i2,')') ! @7/9/89
read(buff(idum1:idum2),formt) ivar ! Read variable &7/9/89
write(nout,formt2) anot,' = ',ivar
write(nlog,formt2) anot,' = ',ivar
return
else
go to 4 ! Read more data
endif
3 write(nout,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
write(nlog,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
stop
end
CBack up
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
CBack up
integer*4 function isearch(x,x0,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:28:03:1995 *
C * *
C * REVISION : ------------- JRH:--:--:1995 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: isearch *
C * TYPE : integer*4 function *
C * *
C * FUNCTION : Inputs array x(i) of n monotonically increasing real*8 *
C * values and scalar x0, and returns lower index of bounding *
C * pair. *
C * *
C * If x0 < x(1), returned value = 0 *
C * If x0 >= x(n), returned value = n *
C * *
C * x ....... input array of monotonic values *
C * x0 ...... required value of x *
C * n ....... number of values in x *
C * *
C ******************************************************************************
real*8 x(*)
real*8 x0
integer*4 n
C
C Local declarations:
C
integer*4 ilow,ihigh,imid
C
C Quickies:
C
if(x0.lt.x(1)) then
isearch=0
return
else if(x0.ge.x(n)) then
isearch=n
return
endif
C
C Now we know x0 lies within x:
C
ilow=1
ihigh=n
C
do while(ihigh.ne.ilow+1)
imid=(ilow+ihigh)/2
if(x0.lt.x(imid)) then
ihigh=imid
else
ilow=imid
endif
end do
isearch=ilow
C
return
end
CBack up
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
CBack up
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
CBack up
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
CBack up
subroutine lpclr(ilp)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:03:12:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: LPCLR *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Clears array ILP *
C * *
C * NOTE : (1) ILP(33,100) must be declared externally *
C * - this is a 16-bit integer variable *
C * (2) LPCLR(ILP) must initially be called to clear *
C * ILP *
C * (3) LPLOAD loads ILP with superimposed data *
C * (4) LPPLOT plots data *
C * *
C ******************************************************************************
integer ilp(33,100)
do i=1,33
do j=1,100
ilp(i,j)=0
end do
end do
return
end
CBack up
subroutine lpload(ilp,x,y,nip,xmin,xmax,ymin,ymax,iwid,ilen,ich,!&19/12/85
$ nout) ! @19/12/85
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:03:12:1985 *
C * *
C * REVISION : Addition of NOUT to parameters JRH:19:12:1985 *
C * Replacement of IROUND with NINT JRH:30:12:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: LPLOAD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Loads ILP with superimposed data *
C * *
C ******************************************************************************
integer ilp(33,100)
dimension x(100),y(100)
if(iwid.gt.131.or.ilen.gt.100) then
write(nout,2)
2 format(' Output too large')
return
endif
if(nip.gt.100) then
write(nout,1)
1 format(' Too many input values')
return
endif
rwid=iwid-1 ! Compute XSCL, YSCL
rlen=ilen-1
nipm1=nip-1
xscl=rwid/(xmax-xmin)
yscl=rlen/(ymax-ymin)
do m=1,nipm1
ia=nint((x(m)-xmin)*xscl+1.) ! &30/12/86
ib=nint((x(m+1)-xmin)*xscl+1.) ! &30/12/86
ja=nint((y(m)-ymin)*yscl+1.) ! &30/12/86
jb=nint((y(m+1)-ymin)*yscl+1.) ! &30/12/86
idel=iabs(ia-ib)+1
jdel=iabs(ja-jb)+1
nlp=max0(idel,jdel)+1
rnlp=nlp
do iilp=1,nlp
riilp=iilp
w1=(rnlp-riilp)/(rnlp-1.)
w2=1.-w1
xi=x(m)*w1+x(m+1)*w2
yi=y(m)*w1+y(m+1)*w2
i=nint((xi-xmin)*xscl+1.) ! &30/12/86
j=nint((yi-ymin)*yscl+1.) ! &30/12/86
if(i.ge.1..and.i.le.iwid.and.j.ge.1..and.j.le.ilen)
$ call ipdig(ilp,i,j,ich)
end do
end do
return
end
CBack up
subroutine lpplot(ilp,xmin,xmax,ymin,ymax,iwid,ilen,nout,
$ ichar,index)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:03:12:1985 *
C * *
C * REVISION : Minor mod. to FORMAT JRH:11:12:1985 *
C * Mod. so that ICHAR defined externally JRH:30:12:1986 *
C * Minor mod. to FORMAT JRH:30:12:1986 *
C * Minor mod. to formats JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: LPPLOT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Plots data to file in form suitable for lineprinter *
C * *
C * Normal values for ICHAR : *
C * /' ','1','2','3','4','5','6','7','8','9','+'/ *
C * *
C ******************************************************************************
integer ilp(33,100),ichar(11,*) ! &30/12/86
dimension idum(131) ! &30/12/86
character*39 fmt ! &11/12/85
if(iwid.gt.131.or.ilen.gt.100) then
write(nout,2)
2 format(' Output too large')
return
endif
write(fmt,11) iwid-14,iwid-14
11 format('(//',i3,'X,4HX = ,E11.4/',i3,'X,4HY = ,E11.4)') ! &11/12/85
write(nout,fmt) xmax,ymax
do jj=1,ilen
j=ilen-jj+1
do i=1,iwid
call opdig(ilp,i,j,ic)
ii=ic+1
idum(i)=ichar(ii,index) ! &30/12/86
if(i.eq.1.or.i.eq.iwid.or.j.eq.1.or.j.eq.ilen) then
if(ii.eq.1) idum(i)=ichar(11,index) ! &30/12/86
endif
end do
write(nout,7) (idum(i),i=1,iwid)
7 format(1x,131a1)
5 end do
write(nout,8) xmin,ymin
8 format(1x,4hx = ,e11.4/1x,4hy = ,e11.4) ! &30/12/
return
end
CBack up
subroutine merc(long,lat,x,y,ae,ecs,lat0,s0,longor,lator)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:14:10:1988 *
C * *
C * REVISION : ------------- JRH:--:--:1988 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: MERC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts longitude and latitude (LONG,LAT) to Mercator *
C * grid coordinates (X,Y). *
C * *
C * LONG ..... Longitude (radians) *
C * LAT ..... Latitude (radians) *
C * X ........ East grid coordinates *
C * Y ........ North grid coordinates *
C * AE ....... Semi-major axis of Earth *
C * ECS ...... Eccentricity squared of figure of Earth *
C * LAT0 ..... Latitude at which scale of projection is *
C * defined *
C * S0 ....... Scale of projection at LAT0 *
C * LONGOR ... Longitude of origin of projection *
C * LATOR .... Latitude of origin of projection *
C * *
C ******************************************************************************
real*8 long,lat,x,y,ae,ecs,lat0,s0,longor,lator
real*8 alpha,e,sl,slor,esl,eslor
alpha=ae*s0/dsqrt(1.d0+(1.d0-ecs)*(dtan(lat0)**2))
e=dsqrt(ecs)
x=alpha*(long-longor)
sl=dsin(lat)
slor=dsin(lator)
esl=e*sl
eslor=e*slor
y=alpha*(dlog(((1.d0+sl)*(1.d0-slor))
$ /((1.d0-sl)*(1.d0+slor)))
$ -e*dlog(((1.d0+esl)*(1.d0-eslor))
$ /((1.d0-esl)*(1.d0+eslor))))/2.d0
return
end
CBack up
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
CBack up
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 : ------------- JRH:--:--:1986 *
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(140),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
CBack up
subroutine nchar(r,idplace,n,ntot)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:31:12:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: NCHAR *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns number of figs. to left of decimal point, N, and *
C * and total no. of characters (including sign and decimal *
C * point), NTOT, in real number, R, defined to IDPLACE *
C * decimal places. *
C * *
C ******************************************************************************
if(r.ne.0.) then ! First find N, no. of figs. to left of decimal point
n=int(alog10(abs(r)))+1
if(n.le.0) n=1
else
n=1
endif
ntot=n+idplace+2 ! Total no. of characters (including sign and dec. point)
return
end
CBack up
subroutine newpen(n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:12:1986 *
C * *
C * REVISION : Mods. for Sun JRH:13:07:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: NEWPEN *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simulates Calcomp plotting routine, directing data to *
C * disc or plotter. *
C * *
C * Labelled COMMON PCOM must be included in main program, *
C * in which LUPLT and LULOG must be set. *
C * *
C ******************************************************************************
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
nopen=mod(n-1,9)+1 ! Pen number in range 1 to 9
call jnewpen(nopen) ! @13/7/89
return
end
CBack up
subroutine newton(f,fprime,y,xguess,x,epsa,epsp,maxit,nit,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:23:07:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: NEWTON *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Function F(X), and it's derivative FPRIME(X) are defined *
C * externally. *
C * NEWTON uses Newton's method to return value of X, based *
C * on initial guess, XGUESS, and input value, Y in : *
C * *
C * Y=F(X) *
C * *
C * EPSA ..... absolute tolerance in X *
C * EPSP ..... proportional tolerance in X *
C * MAXIT .... max. number of iterations *
C * NIT ...... number of iterations at end *
C * IFAIL .... returned 0 if convergence satisfactory, *
C * otherwise 1 *
C * *
C ******************************************************************************
C WRITE(6,2) Y,XGUESS,EPSA,EPSP,MAXIT !!!!!!!!!!
C 2 FORMAT(' Entering NEWTON with Y = ',F10.5/
C $ ' XGUESS = ',F10.5/
C $ ' EPSA = ',F10.8/
C $ ' EPSP = ',F10.8/
C $ ' MAXIT = ',I5) !!!!!!!!!!
xold=xguess
do i=1,maxit ! Step through iterations
xnew=xold+(y-f(xold))/fprime(xold)
if(abs(xnew-xold).lt.epsa) go to 1
if(xnew.ne.0.) then
if(abs((xnew-xold)/xnew).lt.epsp) go to 1
endif
xold=xnew
end do
x=0. ! Failure to converge
nit=i-1 ! Failure to converge
ifail=1 ! Failure to converge
return
1 x=xnew ! Successful convergence
nit=i ! Successful convergence
ifail=0 ! Successful convergence
return
end
CBack up
function nnonspace(a,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:08:08:1989 *
C * *
C * REVISION : Minor mod. to test JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: NNONSPACE *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns the number of groups of "non-space" characters in *
C * the character variable A from A(1:1) to A(N:N) *
C * *
C ******************************************************************************
character*(*) a
logical lspace
icount=0
lspace=.true.
do i=1,n
C if(a(i:i).ne.' '.and.lspace.eq..true.) icount=icount+1
if(a(i:i).ne.' '.and.lspace) icount=icount+1
lspace=(a(i:i).eq.' ')
end do
nnonspace=icount
return
end
CBack up
subroutine number(x,y,wid,r,ang,idplace)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:12:1986 *
C * *
C * REVISION : Minor correction JRH:31:01:1987 *
C * Mods. for direct plotting JRH:03:04:1987 *
C * Mods. for rotation JRH:07:04:1987 *
C * Minor mods. JRH:07:04:1987 *
C * Mod. to output format JRH:28:05:1987 *
C * Mods. for Sun JRH:13:07:1989 *
C * Mods. for Sun JRH:17:07:1989 *
C * Mod. to clipping logic JRH:17:07:1989 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: NUMBER *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simulates Calcomp plotting routine, directing data to *
C * disc or plotter. *
C * *
C * Labelled COMMON PCOM must be included in main program, *
C * in which LUPLT and LULOG must be set. *
C * *
C * X,Y ...... coordinates of bottom LH corner of text *
C * WID ...... character width *
C * R ........ REAL number to be written *
C * ANG ...... angle (deg.) of text anticlockwise from X-axis *
C * IDPLACE .. no. of decimal places to be plotted *
C * *
C ******************************************************************************
character*10 form
character*80 a ! @3/4/87
dimension xa(2),ya(2) ! @3/4/87
logical lc(2),lplot ! @3/4/87
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
logical lrot ! @3/4/87
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units @3/4/87
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units @3/4/87
$ xoff,yoff, ! Plotter offsets @3/4/87
$ pfac, ! Plotter scale factor @3/4/87
$ cw,ch, ! Character width and height @3/4/87
$ lrot ! Rotation logical @3/4/87
if(luplt.eq.-1) then ! Direct plotting @3/4/87
if(idplace.ge.0.) then ! Output as REAL no. @3/4/87
call nchar(r,idplace,n,nasc)! Tot. chs. (incl. sign and dec. point) "
write(form,4) nasc,idplace ! @3/4/87
4 format('(F',i2,'.',i2,')') ! @3/4/87
write(a,form) r ! @3/4/87
else ! Output as INTEGER @3/4/87
call nchar(r,0,n,ndum) ! Tot. chs. (incl. sign and dec. point) @3/4/87
nasc=n+1 ! Total characters (incl. sign) @3/4/87
write(form,5) nasc ! @3/4/87
5 format('(1X,I',i2,')') ! @3/4/87
write(a,form) nint(r) ! @3/4/87
endif ! @3/4/87
dtr=3.14159265/180. ! @3/4/87
widdum=wid*pfac ! @7/4/87
if(lrot) then ! Rotate plot @3/4/87
xdum=-y ! @7/4/87
ydum=x ! @7/4/87
iang=90+nint(ang) ! @3/4/87
else ! @3/4/87
xdum=x ! @7/4/87
ydum=y ! @7/4/87
iang=nint(ang) ! @3/4/87
endif ! @3/4/87
xa(1)=(xdum-xmin)*pfac+xoff ! &7/4/87
ya(1)=(ydum-ymin)*pfac+yoff ! &7/4/87
clen=float(nasc)*widdum !Length of string on plot (in.) @7/4/87 &13/7/89
xa(2)=xa(1)+clen*cos(float(iang)*dtr) ! @3/4/87
ya(2)=ya(1)+clen*sin(float(iang)*dtr) ! @3/4/87
call clip(xa,ya,xlc,xrc,ybc,ytc,lc,lplot) ! @3/4/87
if(.not.lc(1).and..not.lc(2).and.lplot) then
C ..... All on plot @3/4/87 &17/7/89
C call move(xa(1),ya(1)) ! @3/4/87 &13/7/89
C call charot(iang) ! @3/4/87 &13/7/89
C call chasiz(widdum,widdum*1.5) ! @7/4/87 &13/7/89
C call chahol(a,nasc) ! @3/4/87 &13/7/89
call jsymbol(xa(1),ya(1),widdum*1.0,a,float(iang),nasc) ! @17/7/89
endif ! @3/4/87
else ! @3/4/87
az=90.-ang ! Azimuth of text
if(az.lt.0.) az=az+360.
if(idplace.ge.0.) then ! Output as REAL no.
call nchar(r,idplace,n,nasc)! Total chars. (incl. sign and dec. point)
write(form,1) nasc,idplace
1 format('(1X,F',i2,'.',i2,')')
write(luplt,2) x,y,wid,az,nasc,nopen
2 format(' -3 ',2e12.4,e11.4,f7.1,i3,i2) ! &28/5/87
write(luplt,form) r
else ! Output as INTEGER
call nchar(r,0,n,ndum) ! Total chars. (incl. sign and dec. point)
nasc=n+1 ! Total characters (incl. sign)
write(form,3) nasc
3 format('(1X,I',i2,')')
write(luplt,2) x,y,wid,az,nasc,nopen ! &31/1/87
idum=nint(r)
write(luplt,form) idum
endif
endif ! @3/4/87
if(luplt.ne.6) ndata=ndata+1 ! Do not update if to user's terminal
return
end
CBack up
subroutine opdig(ilp,i,j,ic)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:03:12:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: OPDIG *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Outputs digit, IC, to (I,J) from array ILP *
C * *
C ******************************************************************************
integer ilp(33,100)
i1=(i-1)/4+1
i2=5-i+(i1-1)*4
i3=ilp(i1,j)
ip=i2-1
i4=i3/(10**ip)
i5=i3/(10**i2)
ic=i4-i5*10
return
end
CBack up
subroutine overlayplot(nin,ngeo,nout,nlog,lgeo,scale,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:12:10:1989 *
C * *
C * REVISION : Correction to input format JRH:17:11:1989 *
C * Mod. to this header JRH:05:06:1990 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: OVERLAYPLOT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Overlays a plot, using data in input file *
C * *
C * NIN ..... Overlay file input device *
C * NGEO .... Geodesy file input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * LGEO .... .TRUE. if geographical conversion required, *
C * otherwise .FALSE. *
C * SCALE ... Plot scale ((true)/(plot)) *
C * IFAIL ... returned 0 if execution satisfactory, *
C * otherwise 1 *
C * *
C * Codes in overlay file : *
C * *
C * First column : *
C * S or s ..... End of file *
C * T or t ..... Textual data *
C * P or p ..... Polygon data *
C * N or n ..... Newpen data *
C * W or w ..... New window *
C * R or r ..... Reset original window *
C * *
C * Second column (for textual, polygon or new window data):*
C * G or g ..... Grid coordinate data *
C * L or l ..... Longitude and latitude data *
C * *
C * Above character(s) are then followed (in the same *
C * record) by : *
C * *
C * (1) For textual data, WID,AZ,NASC *
C * *
C * (2) For polygon data, NPOINT, the number of points to *
C * be plotted *
C * *
C * (3) For newpen data, the required pen number *
C * *
C * (1) For textual data, this is followed by a record of *
C * the text (80 characters), followed by a record of *
C * the coordinate pair for the text position *
C * *
C * (2) For polygon data, this is followed by a NPOINT *
C * records of the point coordinate pairs *
C * *
C * (3) For new window data, this is followed by two *
C * records of the bottom LH and top RH corner *
C * coordinate pairs *
C * *
C ******************************************************************************
logical lgeo
real*8 long,lat,xd,yd,xdumd
real*8 ae,iflat,be,ecs,eps,lat0,s0,longor,lator,ang,sang,cang
character*1 atype
character*40 sname
character*80 buffer,text
logical lfirst,lorigwind
logical lrot
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units
$ xoff,yoff, ! Plotter offsets
$ pfac, ! Plotter scale factor
$ cw,ch, ! Character width and height
$ lrot ! Rotation logical
save lfirst
data lfirst,lorigwind/.true.,.true./
if(lfirst.and.lgeo) then ! First time through and geog. conv. required
lfirst=.false.
rewind(ngeo) ! For safety
read(ngeo,iostat=ios) sname,ae,iflat,be,ecs,eps,
$ lat0,s0,longor,lator,ang
if(ios.ne.0) then ! Read error or end of file
write(nout,9)
write(nlog,9)
9 format(/' Error in reading geodesy')
ifail=0
return
endif
sang=dsin(ang)
cang=dcos(ang)
endif
rewind(nin) ! For safety
4 read(nin,1,iostat=ios) buffer
1 format(a80)
if(ios.eq.-1) then ! End of file
rewind(nin)
ifail=0
return
else if(ios.gt.0) then ! Read error
write(nout,5)
write(nlog,5)
5 format(/' Read error in subroutine OVERLAYPLOT')
ifail=1
return
endif
if(buffer(1:1).eq.'S'.or.buffer(1:1).eq.'s') then ! End of overlay
rewind(nin)
ifail=0
if(.not.lorigwind) then ! Not currently original window
xlc=xlcold ! Reset original window
ybc=ybcold ! " " "
xrc=xrcold ! " " "
ytc=ytcold ! " " "
endif
return
else if(buffer(1:1).eq.'T'.or.buffer(1:1).eq.'t') then ! Textual data
if(buffer(2:2).eq.'G'.or.buffer(2:2).eq.'g') then ! Grid coordinates
buffer(1:2)=' '
read(buffer,*,iostat=ios) wid,az,nasc
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
7 format(/
$ ' Error reading textual data in subroutine OVERLAYPLOT')
ifail=1
return
endif
read(nin,2,iostat=ios) text
2 format(a80)
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
ifail=1
return
endif
read(nin,*,iostat=ios) x,y ! &17/11/89
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
ifail=1
return
endif
x=x/scale
y=y/scale
wid=wid/scale
angle=90.-az
C Note that WID is char. width in input units,
C AZ is azimuth (deg.) of line of text :
call symbol(x,y,wid,text,angle,nasc)
else if(buffer(2:2).eq.'L'.or.buffer(2:2).eq.'l') then ! Long. and lat.
if(.not.lgeo) then
write(nout,10)
write(nlog,10)
10 format(/
$ ' Geographical input without geodesy in OVERLAYPLOT')
ifail=1
return
endif
buffer(1:2)=' '
read(buffer,*,iostat=ios) wid,az,nasc
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
ifail=1
return
endif
read(nin,2,iostat=ios) text
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
ifail=1
return
endif
read(nin,1,iostat=ios) buffer
if(ios.ne.0) then
write(nout,7)
write(nlog,7)
ifail=1
return
endif
call geoginpbuff(buffer,2,long,lat,xd,yd,atype,ifail)
if(ifail.ne.0.or.atype.ne.'L') then
write(nout,8)
write(nlog,8)
8 format(/' GEOGINPBUFF error in subroutine OVERLAYPLOT')
return
endif
call merc(long,lat,xd,yd,ae,ecs,lat0,s0,longor,lator)
xdumd=xd*cang+yd*sang
yd=-xd*sang+yd*cang
xd=xdumd
x=sngl(xd)/scale
y=sngl(yd)/scale
wid=wid/scale
angle=90.-az
C Note that WID is char. width in input units,
C AZ is azimuth (deg.) of line of text :
call symbol(x,y,wid,text,angle,nasc)
else ! Invalid code
write(nout,3)
write(nlog,3)
3 format(/' Invalid code in subroutine OVERLAYPLOT')
ifail=1
endif
else if(buffer(1:1).eq.'P'.or.buffer(1:1).eq.'p') then ! Polygon data
if(buffer(2:2).eq.'G'.or.buffer(2:2).eq.'g') then ! Grid coordinates
buffer(1:2)=' '
read(buffer,*,iostat=ios) npoint
if(ios.ne.0) then
write(nout,6)
write(nlog,6)
6 format(/
$ ' Error reading polygon data in subroutine OVERLAYPLOT')
ifail=1
return
endif
do i=1,npoint
read(nin,*,iostat=ios) x,y
if(ios.ne.0) then
write(nout,6)
write(nlog,6)
ifail=1
return
endif
x=x/scale
y=y/scale
if(i.eq.1) then
call plot(x,y,3)
else
call plot(x,y,2)
endif
end do
else if(buffer(2:2).eq.'L'.or.buffer(2:2).eq.'l') then ! Long. and lat.
if(.not.lgeo) then
write(nout,10)
write(nlog,10)
ifail=1
return
endif
buffer(1:2)=' '
read(buffer,*,iostat=ios) npoint
if(ios.ne.0) then
write(nout,6)
write(nlog,6)
ifail=1
return
endif
do i=1,npoint
read(nin,1,iostat=ios) buffer
if(ios.ne.0) then
write(nout,6)
write(nlog,6)
ifail=1
return
endif
call geoginpbuff(buffer,2,long,lat,xd,yd,atype,ifail)
if(ifail.ne.0.or.atype.ne.'L') then
write(nout,8)
write(nlog,8)
return
endif
call merc(long,lat,xd,yd,ae,ecs,lat0,s0,longor,lator)
xdumd=xd*cang+yd*sang
yd=-xd*sang+yd*cang
xd=xdumd
x=sngl(xd)/scale
y=sngl(yd)/scale
if(i.eq.1) then
call plot(x,y,3)
else
call plot(x,y,2)
endif
end do
else ! Invalid code
write(nout,3)
write(nlog,3)
ifail=1
endif
else if(buffer(1:1).eq.'N'.or.buffer(1:1).eq.'n') then ! Newpen data
buffer(1:1)=' '
read(buffer,*,iostat=ios) nopen
if(ios.ne.0) then
write(nout,11)
write(nlog,11)
11 format(/
$ ' Error reading NEWPEN data in subroutine OVERLAYPLOT')
ifail=1
return
endif
call newpen(nopen)
else if(buffer(1:1).eq.'W'.or.buffer(1:1).eq.'w') then ! New window data
if(lorigwind) then ! Currently uses original window
xlcold=xlc ! Save original window
ybcold=ybc ! " " "
xrcold=xrc ! " " "
ytcold=ytc ! " " "
lorigwind=.false.
endif
if(buffer(2:2).eq.'G'.or.buffer(2:2).eq.'g') then ! Grid coordinates
do i=1,2
read(nin,*,iostat=ios) x,y
if(ios.ne.0) then
write(nout,12)
write(nlog,12)
12 format(/
$ ' Error reading window data in subroutine OVERLAYPLOT')
ifail=1
return
endif
if(i.eq.1) then
xlc=x
ybc=y
else
xrc=x
ytc=y
endif
end do
else if(buffer(2:2).eq.'L'.or.buffer(2:2).eq.'l') then ! Long. and lat.
do i=1,2
read(nin,1,iostat=ios) buffer
if(ios.ne.0) then
write(nout,12)
write(nlog,12)
ifail=1
return
endif
call geoginpbuff(buffer,2,long,lat,xd,yd,atype,ifail)
if(ifail.ne.0.or.atype.ne.'L') then
write(nout,8)
write(nlog,8)
return
endif
call merc(long,lat,xd,yd,ae,ecs,lat0,s0,longor,lator)
xdumd=xd*cang+yd*sang
yd=-xd*sang+yd*cang
xd=xdumd
if(i.eq.1) then
xlc=sngl(xd)
ybc=sngl(yd)
else
xrc=sngl(xd)
ytc=sngl(yd)
endif
end do
else ! Invalid code
write(nout,3)
write(nlog,3)
ifail=1
endif
if(lrot) then ! Rotate coordinates by 90 deg.
xlcd=xlc
xrcd=xrc
xlc=-ytc
xrc=-ybc
ybc=xlcd
ytc=xrcd
endif
xlc=(xlc/scale-xmin)*pfac+xoff
xrc=(xrc/scale-xmin)*pfac+xoff
ybc=(ybc/scale-ymin)*pfac+yoff
ytc=(ytc/scale-ymin)*pfac+yoff
else if(buffer(1:1).eq.'R'.or.buffer(1:1).eq.'r') then ! Original window
if(.not.lorigwind) then ! Not currently original window
xlc=xlcold ! Reset original window
ybc=ybcold ! " " "
xrc=xrcold ! " " "
ytc=ytcold ! " " "
lorigwind=.true.
endif
else
write(nout,3)
write(nlog,3)
ifail=1
return
endif
go to 4 ! Read more data
end
CBack up
subroutine picalc(pi)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:25:07:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: PICALC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Returns the value of PI *
C * *
C ******************************************************************************
real*4 pi
pi=atan(1.)*4.
return
end
CBack up
logical function pinsd(xp,yp,n,x0,y0)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:23:08:1989 *
C * *
C * REVISION : Minor mod. to comment JRH:29:03:1993 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: PINSD *
C * TYPE : LOGICAL FUNCTION *
C * *
C * FUNCTION : Returns .TRUE. if (X0,Y0) inside polygon, otherwise *
C * returns .FALSE. *
C * No. sides of polygon = N-1. *
C * Total boundary of polygon is defined by N values in *
C * arrays X and Y (XP(N)=XP(1), YP(N)=YP(1)) *
C * *
C ******************************************************************************
real*4 xp(1),yp(1)
logical lxpos,lxneg ! Logicals for determining intersec. dir. on X=X0
m=0
do 2 i=1,n-1
ip1=i+1
C Do quickies first ......
if(xp(i).eq.x0.and.yp(i).eq.y0) go to 1 ! Point at start of segment
if(xp(i).eq.x0.and.xp(ip1).eq.x0) then ! Segment N-S, alligned with pt.
if((yp(i).ge.y0.and.yp(ip1).le.y0).or.
$ (yp(i).le.y0.and.yp(ip1).ge.y0)) go to 1 ! Point on segment
endif
C Now do proper intersection test ......
lxpos=(xp(i).lt.x0.and.xp(ip1).ge.x0)
lxneg=(xp(i).ge.x0.and.xp(ip1).lt.x0)
if(lxpos.or.lxneg) then ! A valid intersection on X=X0
dum=yp(i)-y0-(xp(i)-x0)*(yp(ip1)-yp(i))/
$ (xp(ip1)-xp(i)) ! Intersection on X=X0
if(dum.eq.0.) go to 1 ! Point on segment
if(dum.gt.0.) then ! Valid intersection N of point
if(lxpos) then
m=m+1 ! Positive intersection N of point
else
m=m-1 ! Negative intersection N of point
endif
endif
endif
2 continue
if(m.eq.0) then
pinsd=.false. ! Equal no. of +ve and -ve intersects. on 0 deg. az.
else
pinsd=.true. ! Unequal no. of +ve and -ve intersects. on 0 deg. az.
endif
return
C Exit point for point on boundary.
C Put PINSD =.TRUE. if boundary point considered inside,
C Put PINSD =.FALSE. if boundary point considered outside :
1 pinsd=.true. ! Boundary point considered inside
return
end
CBack up
subroutine plot(x,y,ipen)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:12:1986 *
C * *
C * REVISION : Terminal input of PTIT removed JRH:01:01:1987 *
C * Mods. for direct plotting JRH:03:04:1987 *
C * Minor mods. JRH:06:04:1987 *
C * Mods. for rotation JRH:07:04:1987 *
C * Minor mod. JRH:15:04:1987 *
C * Mods. for Sun JRH:13:07:1989 *
C * Mod. to pause at end of SUNCGI plot JRH:17:07:1989 *
C * Mod. so that jplot(0.,0.,999) may *
C * follow jplot(0.,0.,-999) JRH:18:07:1989 *
C * Mod. for "sleep" option with CGI at *
C * end of plot (for use with Unix script *
C * files) JRH:15:11:1989 *
C * Continuation of above mod. JRH:21:12:1989 *
C * Mod. to "sleep" option JRH:08:11:1990 *
C * Minor mod. to formats JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: PLOT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simulates Calcomp plotting routine, directing data to *
C * disc, or plotter. *
C * *
C * Labelled COMMON PCOM must be included in main program, *
C * in which LUPLT and LULOG must be set. *
C * *
C * Values of IPEN : *
C * *
C * -3 ..... Terminate old plot, start new plot, set *
C * default pen number, NOPEN, reset origin, move *
C * pen to origin *
C * *
C * 2 ...... Draw line from current pen position to point *
C * (X,Y) *
C * *
C * 3 ...... Move pen from current pen position to point *
C * (X,Y) *
C * *
C * 99 ..... Terminate old plot only *
C * *
C * 999 .... Terminate plotting completely *
C * *
C ******************************************************************************
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
logical lrot ! @3/4/87
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units @3/4/87
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units @3/4/87
$ xoff,yoff, ! Plotter offsets @3/4/87
$ pfac, ! Plotter scale factor @3/4/87
$ cw,ch, ! Character width and height @3/4/87
$ lrot ! Rotation logical @3/4/87
common/device/idevice ! For Andrewartha's plot routines @17/7/89
dimension rarray(2) ! @3/4/87
character*20 adum ! @15/11/89 &8/11/90
if(ipen.eq.999) then ! Terminate plotting completely
if(luplt.eq.-1.) then ! Direct plotting @3/4/87
C if(ndata.ne.-2) call endgrf ! &13/7/89
if(ndata.ne.-2) then ! &17/7/89
if(idevice.eq.0) then ! SUNCGI @17/7/89
write(6,4) ! @17/7/89
4 format(/' Press RETURN to continue, ', ! &15/11/89
$ '''s'' RETURN to sleep') ! @15/11/89
read(5,5) adum ! @15/11/89
5 format(a20) ! @15/11/89 &8/11/90
if(adum(1:1).eq.'s'.or. ! @15/11/89 &8/11/90
$ adum(1:1).eq.'S') then ! @15/11/89 &8/11/90
read(adum(2:20),*) idum ! @8/11/90
call sleep(idum) ! @8/11/90
endif ! @8/11/90
endif ! @17/7/89
C Following mods. (next three lines) to allow call of jplot(0.,0.,-999)
C followed by call of jplot(0.,0.,999), since jplots requires a previous
C call of jplot(0.,0.,999) (rather than simply jplot(0.,0.,-999) &18/7/89
endif ! @18/7/89
call jplot(0.,0.,999) ! @13/7/89
C endif ! @17/7/89 &18/7/89
C ..... Terminate old plot (check that plot not previously terminated) @6/4/87
else ! @3/4/87
if(ndata.ne.-2) write(luplt,3)
C ..... Terminate old plot (check that plot not previously terminated) @6/4/87
3 format(' -1 0 0 0')
close(luplt)
endif
else
if(ipen.eq.99) then ! Terminate plot only @3/4/87
if(luplt.eq.-1) then ! Direct plotting @6/4/87
if(idevice.eq.0) then ! SUNCGI @17/7/89
write(6,4) ! @17/7/89
read(5,5) adum ! @15/11/89
if(adum(1:1).eq.'s'.or. ! @15/11/89 &8/11/90
$ adum(1:1).eq.'S') then ! @15/11/89 &8/11/90
read(adum(2:20),*) idum ! @8/11/90
call sleep(idum) ! @8/11/90
endif ! @8/11/90
endif ! @17/7/89
C call endgrf ! &13/7/89
call jplot(0.,0.,-999) ! @13/7/89
else ! @6/4/87
write(luplt,3) ! @6/4/87
endif ! @6/4/87
ndata=-2 ! Flag that plot has been terminated @6/4/87
else
if(ipen.eq.-3) then ! New plot
nopen=1 ! @3/4/87
if(luplt.eq.-1) then ! Direct plotting @3/4/87
C if(ndata.ne.-1.and.ndata.ne.-2) call endgrf ! &13/7/89
if(ndata.ne.-1.and.ndata.ne.-2) then ! &17/7/89
if(idevice.eq.0) then ! SUNCGI @17/7/89
write(6,4) ! @17/7/89
read(5,5) adum ! @16/11/89
if(adum(1:1).eq.'s'.or. ! @15/11/89 &8/11/90
$ adum(1:1).eq.'S') then ! @15/11/89 &8/11/90
read(adum(2:20),*) idum ! @8/11/90
call sleep(idum) ! @8/11/90
endif ! @8/11/90
endif ! @17/7/89
call jplot(0.,0.,-3) ! @13/7/89
endif ! @17/7/89
C ..... Terminate old plot (check that this is not first plot and that
C ..... plot not previously terminated) @6/4/87
call rread(5,6,lulog, ! @3/4/87
$ 'X and Y offsets of plot (metres) ? ','$',!@3/4/87
$ rarray,1,2, ! @3/4/87
$ '((1X,2F10.4)) ') ! &15/4/87
xoff=rarray(1) ! @3/4/87
yoff=rarray(2) ! @3/4/87
xoff=xoff*1000./25.4 ! Convert to in. @3/4/87 &13/7/89
yoff=yoff*1000./25.4 ! Convert to in. @3/4/87 &13/7/89
C call init(idum) ! @3/4/87 &13/7/89
xlc=xoff ! LH side of window @3/4/87
xrc=(xmax-xmin)*pfac+xoff ! RH side of window @3/4/87
ybc=yoff ! Bottom of window @3/4/87
ytc=(ymax-ymin)*pfac+yoff ! Top of window @3/4/87
C call move(xlc,ybc) ! @3/4/87 &13/7/89
call jplot(xlc,ybc,3) ! @13/7/89
C call draw(xrc,ybc) ! @3/4/87 &13/7/89
call jplot(xrc,ybc,2) ! @13/7/89
C call draw(xrc,ytc) ! @3/4/87 &13/7/89
call jplot(xrc,ytc,2) ! @13/7/89
C call draw(xlc,ytc) ! @3/4/87 &13/7/89
call jplot(xlc,ytc,2) ! @13/7/89
C call draw(xlc,ybc) ! @3/4/87 &13/7/89
call jplot(xlc,ybc,2) ! @13/7/89
C call move(10.+xlc,ytc+10.) ! @3/4/87 &13/7/89
C call charot(0) ! @3/4/87 &13/7/89
C call chasiz(cw,ch) ! Set character size @3/4/87 &13/7/89
C call chahol(ptit,40) ! @3/4/87 &13/7/89
call jsymbol(xlc+10./25.4,ytc+10./25.4,ch,ptit,0.,40) ! @13/7/89
C call move(xlc,ybc) ! Go to origin @3/4/87 &13/7/89
call jplot(xlc,ybc,3) ! Go to origin @13/7/89
else ! @3/4/87
if(ndata.ne.-1.and.ndata.ne.-2) write(luplt,3)
C ..... Terminate old plot (check that this is not first plot and that
C ..... plot not previously terminated) @6/4/87
write(luplt,1) ptit
1 format(1x,a40)
write(luplt,2) 0.,0.,3,nopen ! Go to origin
endif ! @3/4/87
ndata=0 ! @3/4/87
else ! Normal plot instruction
if(luplt.eq.-1) then ! Direct plotting @3/4/87
if(lrot) then ! Rotation by 90 deg. required @3/4/87
xdum=-y ! @7/4/87
ydum=x ! @7/4/87
else ! @7/4/87
xdum=x ! @7/4/87
ydum=y ! @7/4/87
endif ! @3/4/87
if(ipen.eq.2.or.ipen.eq.3) then ! @3/4/87
call plota((xdum-xmin)*pfac+xoff, ! &7/4/87
$ (ydum-ymin)*pfac+yoff,ipen, ! &7/4/87
$ xlc,xrc,ybc,ytc) ! @3/4/87
endif
else ! @3/4/87
write(luplt,2) x,y,ipen,nopen
2 format(1x,2e11.4,i4,i2)
endif ! @3/4/87
if(luplt.ne.6) ndata=ndata+1 ! Do not update if to user's terminal
endif
endif ! @3/4/87
endif
return
end
CBack up
subroutine plota(xx,yy,ipen,xlc,xrc,ybc,ytc)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:23:09:1986 *
C * *
C * REVISION : Modification of name PLOT to PLOTA JRH:29:12:1986 *
C * Mods. for Sun JRH:13:07:1989 *
C * Mod. so that NEWPEN is called *
C * regularly when IPEN = 2 (the purpose *
C * of this is for Postcript plotting) JRH:18:07:1989 *
C * Inclusion of additional call to *
C * NEWPEN JRH:11:08:1989 *
C * *
C * SOURCE : PLT1199.FOR *
C * ROUTINE NAME: PLOTA *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Calcomp-type plotting routine, with clipping to window *
C * (XLC,YBC) to (XRC,YTC) *
C * *
C ******************************************************************************
C COMMON/PCOM/ added for regular calls of NEWPEN : @18/7/89
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
dimension x(2),y(2)
logical l3,lc(2),lplot
save xlst,ylst,l3
if(ipen.eq.3) then ! IPEN = 3
xlst=xx
ylst=yy
l3=.true.
if((ndata/100)*100.eq.ndata) call newpen(nopen) ! @11/8/89
else ! IPEN = 2
x(1)=xlst
y(1)=ylst
x(2)=xx
y(2)=yy
xlst=x(2)
ylst=y(2)
call clip(x,y,xlc,xrc,ybc,ytc,lc,lplot)
if(lplot) then
C if(l3.or.lc(1)) call move(x(1),y(1)) ! &13/7/89
if(l3.or.lc(1)) call jplot(x(1),y(1),3) ! @13/7/89
C call draw(x(2),y(2)) ! &13/7/89
call jplot(x(2),y(2),2) ! &13/7/89
if((ndata/100)*100.eq.ndata) then ! @18/7/89
call newpen(nopen) ! @18/7/89
call jplot(x(2),y(2),3)! Move back to point is now required @18/7/89
endif ! @18/7/89
endif
l3=.false. ! Clear "move" flag
endif
return
end
CBack up
subroutine plots(idum)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:12:1986 *
C * *
C * REVISION : DOPEN changed to DSKOPEN JRH:13:01:1987 *
C * Mod. to DSKOPEN call JRH:22:01:1987 *
C * Mods. for direct plotting JRH:03:04:1987 *
C * Mods. for Sun JRH:13:07:1989 *
C * Mods. for Sun JRH:17:07:1989 *
C * Minor mod. to terminal O/P JRH:17:07:1989 *
C * Minor mod. to formats JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: PLOTS *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simulates Calcomp plotting routine, directing data to *
C * disc or plotter. *
C * *
C * Labelled COMMON PCOM must be included in main program, *
C * in which LUPLT and LULOG must be set. *
C * *
C ******************************************************************************
character*20 filnam ! @3/4/87
dimension rarray(4) ! @3/4/87
logical yesno ! @3/4/87
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
logical lrot ! @3/4/87
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units @3/4/87
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units @3/4/87
$ xoff,yoff, ! Plotter offsets @3/4/87
$ pfac, ! Plotter scale factor @3/4/87
$ cw,ch, ! Character width and height @3/4/87
$ lrot ! Rotation logical @3/4/87
if(luplt.ne.6) then ! Only do following if unit no. not user's terminal
filnam=' ' ! @13/7/89
if(luplt.ne.-1) then ! Check if already set for direct plotting @3/4/87
write(6,1) ! @3/4/87
1 format(/' Filename for plot file ' ! @3/4/87
$ /' ("DIRECT" for direct plotting) ? ') ! @13/7/89
C $ /' ("PLOTTER" for plotter) ? '$)!@3/4/87&13/7/89
read(5,2) filnam ! @3/4/87
2 format(a20) ! @3/4/87
write(lulog,3) filnam ! @3/4/87
3 format(/' Filename for plot file ' ! @3/4/87
$ /' ("DIRECT" for direct plotting) ? ',a20)
C ..... @3/4/87 @13/7/89
C $ /' ("PLOTTER" for plotter) ? ',a20)
C ..... @3/4/87 &13/7/89
endif ! @3/4/87
C if(filnam.eq.'PLOTTER'.or.filnam.eq.'plotter'.or. ! @3/4/87 &13/7/89
if(filnam.eq.'DIRECT'.or.filnam.eq.'direct'.or. ! @3/4/87 @13/7/89
$ luplt.eq.-1) then ! Direct plotting @3/4/87
luplt=-1 ! Flag direct plotting @3/4/87
iorient=0 ! Landscape @13/7/89
jfil=50 ! This unit number is generally unused @13/7/89
kdev=0 ! SUNCGI @13/7/89
if(yesno(5,6,lulog, ! @13/7/89
$ 'Postscript (Y) or SUNCGI (N) ? ')) ! @13/7/89
$ kdev=1 ! Postscript @13/7/89
call jplots(iorient,jfil,kdev) ! @13/7/89
call rread(5,6,lulog, ! @3/4/87
$ 'Units of input (eg. ".001" for mm.) ? ','$', ! @3/4/87
$ rarray,1,1, ! @3/4/87
$ '((1X,F10.4)) ') ! @3/4/87
units=rarray(1) ! @3/4/87
call rread(5,6,lulog, ! @17/7/89
$ 'Clipping boundaries : ',' ', ! @17/7/89
$ rarray,0,0, ! @17/7/89
$ '((1X,4F10.4)) ') ! @17/7/89
call rread(5,6,lulog, ! @3/4/87
$ 'XMIN,XMAX,YMIN,YMAX (in units of I/P) ? ','$', ! @3/4/87
$ rarray,1,4, ! @3/4/87
$ '((1X,4F10.4)) ') ! @3/4/87
xmin=rarray(1) ! @3/4/87
xmax=rarray(2) ! @3/4/87
ymin=rarray(3) ! @3/4/87
ymax=rarray(4) ! @3/4/87
lrot=yesno(5,6,lulog, ! @3/4/87
$ 'Rotate plot anticlockwise by 90 deg. ? ') ! @3/4/87
if(lrot) then ! Swop round XMIN,XMAX,YMIN,YMAX @3/4/87
xmind=xmin ! @3/4/87
xmaxd=xmax ! @3/4/87
xmin=-ymax ! @3/4/87
xmax=-ymin ! @3/4/87
ymin=xmind ! @3/4/87
ymax=xmaxd ! @3/4/87
endif ! @3/4/87
call rread(5,6,lulog, ! @3/4/87
$ 'Scale for plotter ',' ', ! &17/7/89
$ rarray,0,0, ! @3/4/87
$ '((1X,F10.0)) ') ! @3/4/87
call rread(5,6,lulog, ! @3/4/87
$ '(input size)/(plot size)) ? ','$', ! @3/4/87
$ rarray,1,1, ! @3/4/87
$ '((1X,F10.0)) ') ! @3/4/87
pscale=rarray(1) ! @3/4/87
pfac=units*1000./(pscale*25.4) ! Convert to in. @3/4/87 &13/7/89
cw=(xmax-xmin)*1000./(pscale*50.*25.4) ! Set character width
C ..... @3/4/87 &13/7/89
C ch=cw*1.5 ! Set character height @3/4/87 &17/7/89
ch=cw*1.0 ! Set character height @17/7/89
else ! @3/4/87
open(unit=luplt,file=filnam,status='UNKNOWN') ! @3/4/87
endif ! @3/4/87
endif
idum=0 ! For time being, do not check if disc OPEN O.K.
ndata=-1 ! Flag call of PLOTS but not yet PLOT(*,*,-3)
return
end
CBack up
real*4 function ran3(idum)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:24:03:1995 *
C * *
C * REVISION : ------------- JRH:--:--:1995 *
C * *
C * SOURCE : forlib.f *
C * ROUTINE NAME: ran3 *
C * TYPE : real*4 function *
C * *
C * FUNCTION : Random number generator (from Numerical Recipes). *
C * *
C ******************************************************************************
integer*4 idum
C
real*4 fac
integer*4 mbig,mseed,mz
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
integer*4 ma(55)
integer*4 iff,inext,inextp,mj,mk,i,ii,k
data iff /0/
save inext,inextp,ma,iff ! Mod. to Numerical Recipes code
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do 11 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
11 continue
do 13 k=1,4
do 12 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12 continue
13 continue
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
return
end
CBack up
subroutine ra2xy(r,a,x,y)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:20:08:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: RA2XY *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts (R,A(azimuth) to (X,Y). *
C * *
C ******************************************************************************
x=r*sin(a)
y=r*cos(a)
return
end
CBack up
subroutine remlink(irem,
$ prev,next,stack,ifirst,ilast,topstack,maxlink,
$ ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:11:10:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: REMLINK *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Removes item from link list *
C * *
C * IREM ..... pointer to item to be removed *
C * PREV ..... array of pointers to previous item in list *
C * NEXT ..... array of pointers to next item in list *
C * STACK .... stack of unused list pointers *
C * IFIRST ... pointer to first item in list *
C * ILAST .... pointer to last item in list *
C * TOPSTACK . pointer to top of stack *
C * MAXLINK .. maximum number of items in list *
C * IFAIL .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 prev(maxlink),next(maxlink),stack(maxlink)
integer*4 irem,ifirst,ilast,topstack,maxlink,ifail
integer*4 prevdum,nextdum
if(irem.lt.1.or.irem.gt.maxlink) then ! Invalid pointer
ifail=1
return
endif
if(topstack.eq.maxlink) then ! List must be empty
ifail=1
return
endif
topstack=topstack+1
if(stack(topstack).ne.0) then ! Overwriting existing pointer
ifail=1
return
endif
stack(topstack)=irem ! Put old pointer back on stack
if(irem.eq.ifirst.and.irem.eq.ilast) then ! Do this easy one first
ifirst=0
ilast=0
else
prevdum=prev(irem)
nextdum=next(irem)
if(irem.eq.ifirst) then ! Remove first item in list
if(nextdum.eq.0) then ! Invalid pointer
ifail=1
return
endif
prev(nextdum)=0
ifirst=nextdum
else if(irem.eq.ilast) then ! Remove last item in list
if(prevdum.eq.0) then ! Invalid pointer
ifail=1
return
endif
next(prevdum)=0
ilast=prevdum
else
if(nextdum.eq.0.or.prevdum.eq.0) then ! Invalid pointer
ifail=1
return
endif
prev(nextdum)=prevdum
next(prevdum)=nextdum
endif
endif
prev(irem)=0 ! Clear link
next(irem)=0 ! Clear link
ifail=0
return
end
CBack up
double precision function rostp(s,t,p)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:03:07:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: ROSTP *
C * TYPE : DOUBLE PRECISION FUNCTION *
C * *
C * FUNCTION : Returns density of seawater as function of salinity *
C * (Practical), temperature (degrees C) and pressure *
C * (bars, above one standard atmosphere). *
C * Reference : *
C * UNESCO, 1981. Tenth report of the joint panel on *
C * oceanographic tables and standards, UNESCO Technical *
C * Papers in Marine Science No. 36, UNESCO, Paris. *
C * *
C ******************************************************************************
double precision s,t,p
double precision b(0:4),c(0:2),d0,a(0:5),f(0:3),g(0:2),i(0:2)
double precision j0,m(0:2),e(0:4),h(0:3),k(0:2)
double precision row,rost0,kst0,kstp,kw,aw,bw,aa,bb,dum
data b
$/ 8.24493d-1, -4.0899d-3, 7.6438d-5, -8.2467d-7, 5.3875d-9/
data c
$/ -5.72466d-3, 1.0227d-4, -1.6546d-6/
data d0
$/ 4.8314d-4/
data a
$/ 999.842594, 6.793952d-2,-9.095290d-3, 1.001685d-4,-1.120083d-6,
$ 6.536332d-9/
data f
$/ 54.6746, -0.603459, 1.09987d-2, -6.1670d-5/
data g
$/ 7.944d-2, 1.6483d-2, -5.3009d-4/
data i
$/ 2.2838d-3, -1.0981d-5, -1.6078d-6/
data j0
$/ 1.91075d-4/
data m
$/ -9.9348d-7, 2.0816d-8, 9.1697d-10/
data e
$/ 19652.21, 148.4206, -2.327105, 1.360477d-2,-5.155288d-5/
data h
$/ 3.239908, 1.43713d-3, 1.16092d-4, -5.77905d-7/
data k
$/ 8.50935d-5, -6.12293d-6, 5.2787d-8/
C Calculate density of reference pure water :
row=a(0)
do ii=1,5
row=row+a(ii)*t**ii ! Density of reference pure water
end do
C Calculate density at one standard atmosphere :
dum=b(0)
do ii=1,4
dum=dum+b(ii)*t**ii
end do
rost0=row+dum*s
dum=c(0)
do ii=1,2
dum=dum+c(ii)*t**ii
end do
rost0=rost0+dum*dsqrt(s**3)+d0*s*s ! Density at one standard atmosphere
C Calculate pure water terms of the secant bulk modulus :
kw=e(0)
do ii=1,4
kw=kw+e(ii)*t**ii
end do
aw=h(0)
do ii=1,3
aw=aw+h(ii)*t**ii
end do
bw=k(0)
do ii=1,2
bw=bw+k(ii)*t**ii
end do
C Calculate secant bulk modulus at one standard atmosphere :
dum=f(0)
do ii=1,3
dum=dum+f(ii)*t**ii
end do
kst0=kw+dum*s
dum=g(0)
do ii=1,2
dum=dum+g(ii)*t**ii
end do
kst0=kst0+dum*dsqrt(s**3) ! Secant bulk modulus at one standard atmosphere
C Calculate secant bulk modulus :
dum=i(0)
do ii=1,2
dum=dum+i(ii)*t**ii
end do
aa=aw+dum*s+j0*dsqrt(s**3)
dum=m(0)
do ii=1,2
dum=dum+m(ii)*t**ii
end do
bb=bw+dum*s
kstp=kst0+aa*p+bb*p*p ! Secant bulk modulus
C Calculate density of seawater :
rostp=rost0/(1.d0-p/kstp) ! Density of seawater
return
end
CBack up
subroutine rread(nin,nout,nlog,anot,acc,rarray,istart,istop,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:27:08:1986 *
C * *
C * REVISION : Error trap included JRH:04:09:1986 *
C * Mod. for output of message only if *
C * ISTART=0 or ISTOP=0 JRH:12:09:1986 *
C * Minor mods. JRH:31:12:1986 *
C * Minor mod. JRH:05:01:1987 *
C * Minor mod. JRH:17:05:1988 *
C * Mod. so that FORM is input without *
C * brackets JRH:21:06:1988 *
C * Mod. for null read after error to cope *
C * with Sun bug JRH:25:07:1989 *
C * Mod. for compilation on transputer JRH:07:09:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: RREAD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads in a set of REAL numbers *
C * *
C * NIN ..... Operator input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Annotation (40 chars.) *
C * ACC ..... Carriage-control character *
C * (' ' for new line, '$' for same line) *
C * RARRAY .. Resultant REAL array *
C * ISTART .. Start index of RARRAY *
C * ISTOP ... Stop index of RARRAY *
C * FORM .... Format for output to log device (40 chars., *
C * without brackets) *
C * *
C ******************************************************************************
real rarray(*)
character*40 anot,form
character*1 acc
character*42 formt ! @7/9/89
logical cont ! @31/12/86
save cont ! @31/12/86
character*12 formout ! &5/1/87
data cont/.false./ ! @11/5/88
2 if(cont) then ! Continuation of terminal text @31/12/86 &25/7/89
write(formout,1) ' ',acc ! &31/12/86
1 format('(',a1,'X,A40,A1',a1,')') ! &5/1/87
else ! Normal text @31/12/86
write(formout,1) '/',acc ! @31/12/86
endif ! @31/12/86
write(nout,formout) anot,' ' ! &5/1/87
write(nlog,formout) anot,' ' ! @5/1/87
if(istart.ne.0.and.istop.ne.0) then ! &31/12/86
read(nin,*,err=3) (rarray(i),i=istart,istop) ! &4/9/86
formt='('//form//')' ! @7/9/89
write(nlog,formt) (rarray(i),i=istart,istop) ! &21/6/88 &7/9/89
cont=.false. ! Terminal text not to be continued @31/12/86
else ! @31/12/86
cont=.true. ! Terminal text to be continued @31/12/86
endif ! @12/9/86
return
3 read(nin,*) ! Null read to cope with Sun bug @25/7/89
go to 2 ! @25/7/89
end
CBack up
subroutine rreadfil(nin,nout,nlog,anot,nchar,rvar,form)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:06:01:1987 *
C * *
C * REVISION : Mod. for keyword checking JRH:19:01:1987 *
C * Simplification JRH:17:06:1988 *
C * Minor correction JRH:07:07:1988 *
C * Mod. for compilation on transputer JRH:07:09:1989 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: RREADFIL *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads REAL variable from file based on keyword. *
C * *
C * NIN ..... File input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Keyword in file *
C * NCHAR ... Number of characters in ANOT (max. 40) *
C * RVAR .... Resultant REAL variable *
C * FORM .... Format for output (40 chars., without brackets) *
C * *
C ******************************************************************************
character*(*) anot
character*40 form
character*80 buff
character*40 blank ! &17/6/88
character*40 formt ! @7/9/89
character*51 formt2
character*14 formt3
data blank/' '/
write(formt2,6) nchar,form
6 format('(1X,A',i2,',A3,',a40,')')
write(formt3,5) nchar
5 format('(/A9,A',i2,',A34/)')
rewind(nin)
4 read(nin,2,end=3) buff
2 format(a80)
if(buff(1:nchar).eq.anot.and. ! &17/6/88
$ buff(nchar+1:nchar+1).eq.' ') then ! Match has been found &19/1/87
buff(1:nchar)=blank(1:nchar) ! Delete keyword
idum1=lhschar(buff) ! Index of first character of value @7/9/89
idum2=lenchar(buff) ! Index of last character of value @7/9/89
write(formt,7) idum2-idum1+1 ! @7/9/89
7 format('(f',i2,'.0)') ! @7/9/89
read(buff(idum1:idum2),formt) rvar ! Read variable &7/9/89
write(nout,formt2) anot,' = ',rvar
write(nlog,formt2) anot,' = ',rvar
return
else
go to 4 ! Read more data
endif
3 write(nout,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
write(nlog,formt3) ' Keyword ',anot,
$ ' not found .... program terminated'
stop
end
CBack up
subroutine rreadfil2(nin,anot,nchar,var,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: rreadfil2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads REAL*4 variable from file based on keyword. *
C * (Modified version of rreadfil, 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 * var ..... Resultant REAL*4 variable *
C * ifail ... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
real*4 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
CBack up
subroutine r2dcalc(r2d)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:25:07:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1991 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: R2DCALC *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Returns the value of R2D (radians to degrees constant) *
C * *
C ******************************************************************************
real*4 r2d
r2d=180./(atan(1.)*4.)
return
end
CBack up
subroutine sim(f,xl,xh,epsa,epsp,s,m4)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MATHEMATICS REF:JRH:23:07:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: SIM *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Performs Simpson integration of function F between limits *
C * limits XL and XH. *
C * EPSA ..... required absolute difference between *
C * successive approximations *
C * EPSP ..... required proportional difference between *
C * successive approximations *
C * S ........ solution *
C * M4 ....... final number of integration steps *
C * *
C ******************************************************************************
C WRITE(6,6) !!!!!!!!!!
C 6 FORMAT(/) !!!!!!!!!!
C WRITE(6,8) XL,XH,EPSA,EPSP !!!!!!!!!!
C 8 FORMAT(' Entering SIM with XL = ',F10.5/
C $ ' XH = ',F10.5/
C $ ' EPSA = ',F10.8/
C $ ' EPSP = ',F10.8) !!!!!!!!!!!
sold=0.
fl=f(xl)
fh=f(xh)
m4=2
4 m2=m4-1
C WRITE(6,5) !!!!!!!!!!
C 5 FORMAT('+.'$) !!!!!!!!!!
rm=m4
del=(xh-xl)/rm
s4=0.
s2=0.
do 1 i=1,m4
ri=i
1 s4=s4+f(xl+(ri-0.5)*del)
do 2 i=1,m2
ri=i
2 s2=s2+f(xl+ri*del)
s=(fl+fh+4.*s4+2.*s2)*del/6.
if(abs(s-sold).lt.epsa) return
if(s.ne.0.) then
if(abs((s-sold)/s).lt.epsp) return
endif
sold=s
m4=2*m4
go to 4
end
CBack up
subroutine sima(aa,s,m4,range)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:15:05:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: SIMA *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simpson integration - definite - routine A *
C * *
C * AA .... Array to be integrated *
C * S ..... Definite integral *
C * M4 .... Number of "double" increments *
C * RANGE . Integration range *
C * *
C ******************************************************************************
dimension aa(1)
C WRITE(6,1000)
C1000 FORMAT(' A'$)
m2=m4-1
del=range/float(m4)
s4=0.
s2=0.
do i=1,m4
ii=i*2
C WRITE(6,4000) S4
C4000 FORMAT(' S4 = ',E10.3)
s4=s4+aa(ii)
end do
do i=1,m2
ii=i*2+1
C WRITE(6,5000) S2
C5000 FORMAT(' S2 = ',E10.3)
s2=s2+aa(ii)
end do
idum=m4*2+1
s=(aa(1)+aa(idum)+4.*s4+2.*s2)*del/6.
C WRITE(6,2000) S,M4,DEL
C2000 FORMAT(' S,M4,DEL = ',F10.5,I5,F10.5)
C WRITE(6,3000) AA(1),AA(IDUM),S4,S2
C3000 FORMAT(' AA(1),AA(IDUM),S4,S2 = ',4E10.3)
return
end
CBack up
subroutine simb(aa,bb,m4,range)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:15:05:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: SIMB *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simpson integration - indefinite - routine B *
C * *
C * AA .... Input array *
C * BB .... Indefinite integral array *
C * M4 .... Numberof "double" increments *
C * RANGE . Integration range *
C * *
C ******************************************************************************
dimension aa(1),bb(1)
C WRITE(6,1000)
C1000 FORMAT(' B'$)
m2=m4-1
del=range/float(m4) ! Width of "double" increment
f1=del*5./24.
f2=del*8./24.
f3=-del/24.
g1=del/6.
g2=del*4./6.
g3=del/6.
bb(1)=0.
do i=1,m4
i2=i*2
i1=i2-1
i3=i2+1
bb(i2)=bb(i1)+aa(i1)*f1+aa(i2)*f2+aa(i3)*f3
bb(i3)=bb(i1)+aa(i1)*g1+aa(i2)*g2+aa(i3)*g3
end do
return
end
CBack up
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 * *
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(140)
double precision d(140),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
CBack up
subroutine srect(x,y,shade)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:07:1989 *
C * *
C * REVISION : Correction to declaration of PAT1 etc.JRH:16:08:1989 *
C * Temporary mod. to cope with Sun 4 bug JRH:14:05:1990 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: SRECT *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Shades a rectangle. *
C * *
C * (x(1),y(1)) ..... coordinates of bottom LH corner *
C * (x(2),y(2)) ..... coordinates of top RH corner *
C * shade ........... shade (0 = white, 1 = black) *
C * (NOTE : the reverse of Postcript) *
C * *
C * NOTE : SRECT should be called with SHADE = 0 imediately *
C * after PLOT(*,*,-3), in order to initialise system. *
C * *
C * NOTE : The last statement in this routine (NDATA=NDATA+1) *
C * is not executed using the f68881 optimiser. Hence from *
C * now on, I will not use f68881 optimiser. *
C * JRH 16/8/89 *
C * *
C ******************************************************************************
dimension x(2),y(2)
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
logical lrot
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units
$ xoff,yoff, ! Plotter offsets
$ pfac, ! Plotter scale factor
$ cw,ch, ! Character width and height
$ lrot ! Rotation logical
common/device/idevice ! (JRA common)
common/scale/iscale ! (JRA common)
common/plotpos/ifile,ipath,iorient,ihlast,iplast ! (JRA common)
common/plotcgi/ix(2),iy(2),ixoff,iyoff ! (JRA common)
dimension ipx(4),ipy(4),xdum(2),ydum(2)
dimension ipxold(2),ipyold(2)
save ipxold,ipyold
dimension irxoff(2),iryoff(2),ixfig(2),iyfig(2)
integer pat1(16), pat2(16), pat3(16), pat4(16), ! &16/8/89
$ pat5(16), pat6(16), pat7(16), pat8(16),
$ pat9(16),pat10(16),pat11(16),pat12(16),
$ pat13(16),pat14(16),pat15(16),pat16(16),
$ pat17(16)
character*29 formt ! @14/5/90
data pat1/0,0,0,0,
$ 0,0,0,0,
$ 0,0,0,0,
$ 0,0,0,0/
data pat2/0,0,0,0,
$ 0,0,1,0,
$ 0,0,0,0,
$ 0,0,0,0/
data pat3/0,0,0,0,
$ 0,0,1,0,
$ 0,1,0,0,
$ 0,0,0,0/
data pat4/0,0,0,0,
$ 0,0,1,0,
$ 0,1,1,0,
$ 0,0,0,0/
data pat5/0,0,0,0,
$ 0,1,1,0,
$ 0,1,1,0,
$ 0,0,0,0/
data pat6/0,0,0,1,
$ 0,1,1,0,
$ 0,1,1,0,
$ 0,0,0,0/
data pat7/0,0,0,1,
$ 0,1,1,0,
$ 0,1,1,0,
$ 1,0,0,0/
data pat8/0,0,1,1,
$ 0,1,1,0,
$ 0,1,1,0,
$ 1,0,0,0/
data pat9/0,0,1,1,
$ 0,1,1,0,
$ 0,1,1,0,
$ 1,1,0,0/
data pat10/0,0,1,1,
$ 0,1,1,0,
$ 1,1,1,0,
$ 1,1,0,0/
data pat11/0,0,1,1,
$ 0,1,1,1,
$ 1,1,1,0,
$ 1,1,0,0/
data pat12/0,1,1,1,
$ 0,1,1,1,
$ 1,1,1,0,
$ 1,1,0,0/
data pat13/0,1,1,1,
$ 0,1,1,1,
$ 1,1,1,0,
$ 1,1,1,0/
data pat14/0,1,1,1,
$ 1,1,1,1,
$ 1,1,1,0,
$ 1,1,1,0/
data pat15/0,1,1,1,
$ 1,1,1,1,
$ 1,1,1,1,
$ 1,1,1,0/
data pat16/1,1,1,1,
$ 1,1,1,1,
$ 1,1,1,1,
$ 1,1,1,0/
data pat17/1,1,1,1,
$ 1,1,1,1,
$ 1,1,1,1,
$ 1,1,1,1/
if(ndata.eq.0) then ! Initialisation
if(idevice.eq.1) then ! Postcript
write(ifile,1) ! Write procedure "o" ("offset" rectangle)
1 format(' /o '/
$ ' {/sh exch 100 div def '/
$ ' /oy2 exch def '/
$ ' /ox2 exch def '/
$ ' /oy1 exch def '/
$ ' /ox1 exch def '/
$ ' /x1 x1 ox1 add def '/
$ ' /y1 y1 oy1 add def '/
$ ' /x2 x2 ox2 add def '/
$ ' /y2 y2 oy2 add def '/
$ ' newpath '/
$ ' x1 y1 moveto '/
$ ' x1 y2 lineto '/
$ ' x2 y2 lineto '/
$ ' x2 y1 lineto '/
$ ' closepath sh setgray fill '/
$ ' 0 setgray} def ')
write(ifile,2) ! Set x1, y1, x2, y2 to zero
2 format(' /x1 0 def '/
$ ' /y1 0 def '/
$ ' /x2 0 def '/
$ ' /y2 0 def ')
do i=1,2
ipxold(i)=0
ipyold(i)=0
end do
else ! SunCGI
call cfpattable(1,4,4,pat1)
call cfpattable(2,4,4,pat2)
call cfpattable(3,4,4,pat3)
call cfpattable(4,4,4,pat4)
call cfpattable(5,4,4,pat5)
call cfpattable(6,4,4,pat6)
call cfpattable(7,4,4,pat7)
call cfpattable(8,4,4,pat8)
call cfpattable(9,4,4,pat9)
call cfpattable(10,4,4,pat10)
call cfpattable(11,4,4,pat11)
call cfpattable(12,4,4,pat12)
call cfpattable(13,4,4,pat13)
call cfpattable(14,4,4,pat14)
call cfpattable(15,4,4,pat15)
call cfpattable(16,4,4,pat16)
call cfpattable(17,4,4,pat17)
endif
endif
do i=1,2
if(lrot) then ! Rotation by 90 deg. required
xdum(i)=-y(i)
ydum(i)=x(i)
else
xdum(i)=x(i)
ydum(i)=y(i)
endif
xdum(i)=(xdum(i)-xmin)*pfac+xoff
ydum(i)=(ydum(i)-ymin)*pfac+yoff
end do
do i=1,2
if(xdum(i).lt.xlc) xdum(i)=xlc ! Clip onto plot
if(xdum(i).gt.xrc) xdum(i)=xrc ! Clip onto plot
if(ydum(i).lt.ybc) ydum(i)=ybc ! Clip onto plot
if(ydum(i).gt.ytc) ydum(i)=ytc ! Clip onto plot
end do
if(xdum(1).ne.xdum(2).and.ydum(1).ne.ydum(2)) then
if(idevice.eq.1) then ! Postcript
do i=1,2
ipx(i)=xdum(i)*float(iscale)
ipy(i)=ydum(i)*float(iscale)
end do
if(ipath.eq.1) then
write(ifile,3)
3 format('stroke')
ipath=0
endif
grey=aint((1.-shade)*17.)/16. ! 17 shades
if(grey.lt.0.) grey=0.
if(grey.gt.1.) grey=1.
igrey=nint(grey*100.)
do i=1,2
irxoff(i)=ipx(i)-ipxold(i)
iryoff(i)=ipy(i)-ipyold(i)
if(irxoff(i).eq.0) then
ixfig(i)=1
else if(irxoff(i).gt.0) then
ixfig(i)=int(alog10(float(irxoff(i))))+1
else
ixfig(i)=int(alog10(float(-irxoff(i))))+2 !To allow for minus sign
endif
if(iryoff(i).eq.0) then
iyfig(i)=1
else if(iryoff(i).gt.0) then
iyfig(i)=int(alog10(float(iryoff(i))))+1
else
iyfig(i)=int(alog10(float(-iryoff(i))))+2 !To allow for minus sign
endif
end do
if(igrey.eq.0) then
igfig=1
else ! Note that IGREY must be positive
igfig=int(alog10(float(igrey)))+1
endif
write(formt,5) ixfig(1),iyfig(1),ixfig(2),iyfig(2),igfig ! @14/5/90
5 format('(i',i1,',x,i',i1,',x,i',i1,',x,i',i1,',x,i',i1, ! @14/4/90
$ ','' o'')') ! @14/5/90
C write(ifile,4) (irxoff(i),iryoff(i),i=1,2),igrey ! &14/5/90
write(ifile,formt) (irxoff(i),iryoff(i),i=1,2),igrey ! @14/5/90
C 4 format(i,x,i,x, ! &14/5/90
C $ i,x,i,x, ! &14/5/90
C $ i,' o') ! &14/5/90
do i=1,2
ipxold(i)=ipx(i)
ipyold(i)=ipy(i)
end do
else ! SunCGI
ind=aint(shade*17.+1.)
if(ind.lt.1) ind=1
if(ind.gt.17) ind=17
call cfhatchix(ind)
ipx(1)=xdum(1)*float(iscale)+float(ixoff)
ipx(2)=ipx(1)
ipx(3)=xdum(2)*float(iscale)+float(ixoff)
ipx(4)=ipx(3)
ipy(1)=ydum(1)*float(iscale)+float(iyoff)
ipy(2)=ydum(2)*float(iscale)+float(iyoff)
ipy(3)=ipy(2)
ipy(4)=ipy(1)
call cfpolygon(ipx,ipy,4)
endif
endif
ndata=ndata+1
return
end
CBack up
subroutine svbksb(u,w,v,m,n,mp,np,b,x,tmp)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:08:01:1992 *
C * *
C * REVISION : ------------- JRH:--:--:1992 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: SVBKSB *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Singular value decomposition solution *
C * (from Numerical Recipes) *
C * *
C * U ..... Input M by N matrix *
C * W ..... Input N element weight vector *
C * V ..... Input N by N matrix ( NOT the transpose of "V") *
C * M ..... Number of rows of U (M.ge.N) *
C * N ..... Number of columns of U (M.ge.N) *
C * MP .... Physical number of rows of U *
C * NP .... Physical number of columns of U *
C * B ..... Input M element vector (right hand side) *
C * X ..... Output N element solution vector *
C * TMP ... Workspace N element vector *
C * *
C ******************************************************************************
real*4 u(mp,np),w(np),v(np,np),b(mp),x(np),tmp(*)
integer*4 m,n,mp,np
real*4 s
integer*4 i,j,jj
do 12 j=1,n
s=0.
if(w(j).ne.0.)then
do 11 i=1,m
s=s+u(i,j)*b(i)
11 continue
s=s/w(j)
endif
tmp(j)=s
12 continue
do 14 j=1,n
s=0.
do 13 jj=1,n
s=s+v(j,jj)*tmp(jj)
13 continue
x(j)=s
14 continue
return
end
CBack up
subroutine svdcmp(a,m,n,mp,np,u,w,v,rv1,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:08:01:1992 *
C * *
C * REVISION : ------------- JRH:--:--:1992 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: SVDCMP *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Singular value decomposition solution *
C * (from Numerical Recipes) *
C * *
C * A ..... Input M by N matrix *
C * M ..... Number of rows of A (M.ge.N) *
C * N ..... Number of columns of A (M.ge.N) *
C * MP .... Physical number of rows of A *
C * NP .... Physical number of columns of A *
C * U ..... Output M by N matrix *
C * W ..... Output N element weight vector *
C * V ..... Output N by N matrix (NOT the transpose of "V") *
C * RV1 ... Workspace N element vector *
C * IFAIL . 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
real*4 a(mp,np),u(mp,np),w(np),v(np,np),rv1(*)
integer*4 m,n,mp,np,ifail
real*4 f,g,h,scale,anorm,s,c,x,y,z
integer*4 i,j,k,l,its,nm
if(m.lt.n) then
ifail=1
return
endif
g=0.0
scale=0.0
anorm=0.0
do i=1,m ! Copy A to U
do j=1,n
u(i,j)=a(i,j)
end do
end do
do 25 i=1,n
l=i+1
rv1(i)=scale*g
g=0.0
s=0.0
scale=0.0
if (i.le.m) then
do 11 k=i,m
scale=scale+abs(u(k,i))
11 continue
if (scale.ne.0.0) then
do 12 k=i,m
u(k,i)=u(k,i)/scale
s=s+u(k,i)*u(k,i)
12 continue
f=u(i,i)
g=-sign(sqrt(s),f)
h=f*g-s
u(i,i)=f-g
if (i.ne.n) then
do 15 j=l,n
s=0.0
do 13 k=i,m
s=s+u(k,i)*u(k,j)
13 continue
f=s/h
do 14 k=i,m
u(k,j)=u(k,j)+f*u(k,i)
14 continue
15 continue
endif
do 16 k= i,m
u(k,i)=scale*u(k,i)
16 continue
endif
endif
w(i)=scale *g
g=0.0
s=0.0
scale=0.0
if ((i.le.m).and.(i.ne.n)) then
do 17 k=l,n
scale=scale+abs(u(i,k))
17 continue
if (scale.ne.0.0) then
do 18 k=l,n
u(i,k)=u(i,k)/scale
s=s+u(i,k)*u(i,k)
18 continue
f=u(i,l)
g=-sign(sqrt(s),f)
h=f*g-s
u(i,l)=f-g
do 19 k=l,n
rv1(k)=u(i,k)/h
19 continue
if (i.ne.m) then
do 23 j=l,m
s=0.0
do 21 k=l,n
s=s+u(j,k)*u(i,k)
21 continue
do 22 k=l,n
u(j,k)=u(j,k)+s*rv1(k)
22 continue
23 continue
endif
do 24 k=l,n
u(i,k)=scale*u(i,k)
24 continue
endif
endif
anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25 continue
do 32 i=n,1,-1
if (i.lt.n) then
if (g.ne.0.0) then
do 26 j=l,n
v(j,i)=(u(i,j)/u(i,l))/g
26 continue
do 29 j=l,n
s=0.0
do 27 k=l,n
s=s+u(i,k)*v(k,j)
27 continue
do 28 k=l,n
v(k,j)=v(k,j)+s*v(k,i)
28 continue
29 continue
endif
do 31 j=l,n
v(i,j)=0.0
v(j,i)=0.0
31 continue
endif
v(i,i)=1.0
g=rv1(i)
l=i
32 continue
do 39 i=n,1,-1
l=i+1
g=w(i)
if (i.lt.n) then
do 33 j=l,n
u(i,j)=0.0
33 continue
endif
if (g.ne.0.0) then
g=1.0/g
if (i.ne.n) then
do 36 j=l,n
s=0.0
do 34 k=l,m
s=s+u(k,i)*u(k,j)
34 continue
f=(s/u(i,i))*g
do 35 k=i,m
u(k,j)=u(k,j)+f*u(k,i)
35 continue
36 continue
endif
do 37 j=i,m
u(j,i)=u(j,i)*g
37 continue
else
do 38 j= i,m
u(j,i)=0.0
38 continue
endif
u(i,i)=u(i,i)+1.0
39 continue
do 49 k=n,1,-1
do 48 its=1,30
do 41 l=k,1,-1
nm=l-1
if ((abs(rv1(l))+anorm).eq.anorm) go to 2
if ((abs(w(nm))+anorm).eq.anorm) go to 1
41 continue
1 c=0.0
s=1.0
do 43 i=l,k
f=s*rv1(i)
if ((abs(f)+anorm).ne.anorm) then
g=w(i)
h=sqrt(f*f+g*g)
w(i)=h
h=1.0/h
c= (g*h)
s=-(f*h)
do 42 j=1,m
y=u(j,nm)
z=u(j,i)
u(j,nm)=(y*c)+(z*s)
u(j,i)=-(y*s)+(z*c)
42 continue
endif
43 continue
2 z=w(k)
if (l.eq.k) then
if (z.lt.0.0) then
w(k)=-z
do 44 j=1,n
v(j,k)=-v(j,k)
44 continue
endif
go to 3
endif
if (its.eq.30) then
ifail=1
return
endif
x=w(l)
nm=k-1
y=w(nm)
g=rv1(nm)
h=rv1(k)
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
g=sqrt(f*f+1.0)
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
c=1.0
s=1.0
do 47 j=l,nm
i=j+1
g=rv1(i)
y=w(i)
h=s*g
g=c*g
z=sqrt(f*f+h*h)
rv1(j)=z
c=f/z
s=h/z
f= (x*c)+(g*s)
g=-(x*s)+(g*c)
h=y*s
y=y*c
do 45 nm=1,n
x=v(nm,j)
z=v(nm,i)
v(nm,j)= (x*c)+(z*s)
v(nm,i)=-(x*s)+(z*c)
45 continue
z=sqrt(f*f+h*h)
w(j)=z
if (z.ne.0.0) then
z=1.0/z
c=f*z
s=h*z
endif
f= (c*g)+(s*y)
x=-(s*g)+(c*y)
do 46 nm=1,m
y=u(nm,j)
z=u(nm,i)
u(nm,j)= (y*c)+(z*s)
u(nm,i)=-(y*s)+(z*c)
46 continue
47 continue
rv1(l)=0.0
rv1(k)=f
w(k)=x
48 continue
3 continue
49 continue
ifail=0
return
end
CBack up
subroutine symbol(x,y,wid,a,ang,nasc)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : PLOTTING REF:JRH:31:12:1986 *
C * *
C * REVISION : Mods. for direct plotting JRH:03:04:1987 *
C * Mods. for rotation JRH:07:04:1987 *
C * Minor mods. JRH:07:04:1987 *
C * Mod. to output format JRH:28:05:1987 *
C * Mods. for Sun JRH:13:07:1989 *
C * Mods. for Sun JRH:17:07:1989 *
C * Mod. to clipping logic JRH:17:07:1989 *
C * Mod. to cope with -ve NASC *
C * (i.e. "symbol" plotting) JRH:13:03:1991 *
C * Minor mod. to format JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: SYMBOL *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Simulates Calcomp plotting routine, directing data to *
C * disc or plotter. *
C * *
C * Labelled COMMON PCOM must be included in main program, *
C * in which LUPLT and LULOG must be set. *
C * *
C * X,Y ...... coordinates of bottom LH corner of text *
C * WID ...... character width *
C * A ........ text to be written *
C * ANG ...... angle (deg.) of text anticlockwise from X-axis *
C * NASC ..... no. of ASCII characters *
C * *
C ******************************************************************************
character*(*) a
character*7 form
dimension xa(2),ya(2) ! @3/4/87
logical lc(2),lplot ! @3/4/87
character*40 ptit(1) ! Plot title
common/pcom/luplt, ! Plot file unit no.
$ lulog, ! Log file unit no.
$ ndata, ! No. instructions in plot
$ ptit, ! Current plot title
$ nopen ! Current pen number
logical lrot ! @3/4/87
common/wcom/xmin,xmax,ymin,ymax, ! Plot bounds in input units @3/4/87
$ xlc, xrc, ybc, ytc, ! Plot bounds in plotter units @3/4/87
$ xoff,yoff, ! Plotter offsets @3/4/87
$ pfac, ! Plotter scale factor @3/4/87
$ cw,ch, ! Character width and height @3/4/87
$ lrot ! Rotation logical @3/4/87
if(luplt.eq.-1) then ! Direct plotting @3/4/87
dtr=3.14159265/180. ! @3/4/87
widdum=wid*pfac ! @7/4/87
if(lrot) then ! Rotate plot @3/4/87
xdum=-y ! @7/4/87
ydum=x ! @7/4/87
iang=90+nint(ang) ! @3/4/87
else ! @3/4/87
xdum=x ! @7/4/87
ydum=y ! @7/4/87
iang=nint(ang) ! @3/4/87
endif ! @3/4/87
xa(1)=(xdum-xmin)*pfac+xoff ! @7/4/87
ya(1)=(ydum-ymin)*pfac+yoff ! @7/4/87
if(nasc.ge.1) then ! String plotting &13/6/91
clen=float(nasc)*widdum ! Length of string on plot (mm.) @7/4/87
else ! &13/6/91
clen=widdum ! &13/6/91
endif ! &13/6/91
xa(2)=xa(1)+clen*cos(float(iang)*dtr) ! @3/4/87
ya(2)=ya(1)+clen*sin(float(iang)*dtr) ! @3/4/87
call clip(xa,ya,xlc,xrc,ybc,ytc,lc,lplot) ! @3/4/87
if(.not.lc(1).and..not.lc(2).and.lplot) then
C ..... All on plot @3/4/87 &17/7/89
C call move(xa(1),ya(1)) ! @3/4/87 &13/7/89
C call charot(iang) ! @3/4/87 &13/7/89
C call chasiz(widdum,widdum*1.5) ! @7/4/87 &13/7/89
C call chahol(a,nasc) ! @3/4/87 &13/7/89
call jsymbol(xa(1),ya(1),widdum*1.0,a,float(iang),nasc) ! @17/7/89
endif ! @3/4/87
else ! @3/4/87
az=90.-ang ! Azimuth of text
if(az.lt.0.) az=az+360.
if(nasc.ge.1) then ! String plotting &13/6/91
write(form,1) nasc
else ! &13/6/91
write(form,1) 1 ! &13/6/91
endif ! &13/6/91
1 format('(1X,A',i2,')')
write(luplt,2) x,y,wid,az,nasc,nopen
2 format(' -3 ',2e12.4,e11.4,f7.1,i3,i2) ! &28/5/87
write(luplt,form) a
endif ! @3/4/87
if(luplt.ne.6) ndata=ndata+1 ! Do not update if to user's terminal
return
end
CBack up
subroutine timin(nin,nout,nlog,anot,ihr,min,sec)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:20:03:1991 *
C * *
C * REVISION : Minor mod. to formats JRH:15:05:2000 *
C * "end=" added to read JRH:14:01:2001 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: TIMIN *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Reads in a time, in various formats *
C * *
C * NIN ..... Operator input device *
C * NOUT .... Operator output device *
C * NLOG .... Log device *
C * ANOT .... Annotation (40 chars.) *
C * IHR ..... Hour *
C * MIN ..... Minute *
C * SEC ..... Second *
C * *
C * NOTE : Main constraint is that input order is hour, minute, second *
C * *
C ******************************************************************************
real*4 sec
integer*4 nin,nout,nlog,ihr,min
character*40 anot
integer*4 idum,icount,isec
character*80 buff
4 write(nout,1) anot
write(nlog,1) anot
1 format(1x,a40)
read(nin,2,err=4) buff(1:80)
2 format(a80)
write(nlog,3) buff(1:80)
3 format(1x,a80)
do i=1,80 ! First remove all non-numeric characters except '.'
idum=ichar(buff(i:i))
if((idum.lt.48.or.idum.gt.57).and.
$ buff(i:i).ne.'.') buff(i:i)=' '
end do
icount=0 ! First look for 4-or 6-digit time
do i=1,80
idum=ichar(buff(i:i))
if(idum.ge.48.and.idum.le.57) then ! A digit
icount=icount+1
else
icount=0
endif
if(icount.eq.4.and. ! 4 digits
$ (buff(i+1:i+1).eq.' '.or. ! Followed by space
$ buff(i+1:i+1).eq.'.')) then ! or '.'
read(buff(i-3:i),5,err=4) ihr,min
5 format(2i2)
sec=0.
go to 7
endif
if(icount.eq.6.and. ! 6 digits
$ (buff(i+1:i+1).eq.' '.or. ! Followed by space
$ buff(i+1:i+1).eq.'.')) then ! or '.'
read(buff(i-5:i),6,err=4) ihr,min,isec
6 format(3i2)
sec=float(isec)
go to 7
endif
end do
read(buff(1:80),*,err=8,end=8) ihr,min,sec
C ..... Read with delimiters - first try HMS
go to 7
8 read(buff(1:80),*,err=4) ihr,min ! Then try HM
sec=0.
7 if(ihr.lt.0.or.ihr.gt.24) go to 4
if(min.lt.0.or.min.gt.59) go to 4
if(sec.lt.0..or.sec.gt.60.) go to 4
return
end
CBack up
subroutine utd(iye,mon,idy,ihr,min,sec,utim,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:23:05:1991 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: UTD *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts Unix Time to daytime *
C * *
C * (Unix Time starts at 0000 on 1 Jan. 1970) *
C * *
C * IYE ........ Year *
C * MON ........ Month *
C * IDY ........ Day *
C * IHR ........ Hour *
C * MIN ........ Minute *
C * SEC ........ Second *
C * UTIM ....... Unix time (seconds) *
C * IFAIL ...... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
real*8 utim
real*4 sec
integer*4 iye,mon,idy,ihr,min,ifail
real*8 jd0,julian
logical first
save jd0,first
data first/.true./
if(first) then ! Compute zero of Unix Time in Julian days
call julday(1970,1,1,0,0,0.,jd0,ifail)
if(ifail.ne.0) return ! Error
first=.false.
endif
julian=utim/86400.d0+jd0
call caldat(iye,mon,idy,ihr,min,sec,julian,ifail)
return
end
CBack up
subroutine writearray(nout,a,athresh,imin,imax,
$ form1,form2,nchar1,nchar2,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:27:11:1988 *
C * *
C * REVISION : Mod. to parameter JRH:07:07:1989 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: WRITEARRAY *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Writes array A from indices IMIN to IMAX in two possible *
C * formats, depending on whether the value is smaller or *
C * larger than a threshold, ATHRESH. 10 values are output *
C * per record. *
C * *
C * NOUT ..... Output device *
C * A ........ Array (NOTE that lower bound is 0) *
C * ATHRESH .. Threshold value *
C * IMIN ..... Minimum index of array, A *
C * IMAX ..... Maximum index of array, A *
C * FORM1 .... Format if (A.LT.ATHRESH) *
C * FORM2 .... Format if (A.GE.ATHRESH) *
C * NCHAR1 ... Number of characters in FORM1 *
C * NCHAR2 ... Number of characters in FORM2 *
C * IFAIL .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer pmax ! @7/7/89
parameter(pmax=160)
dimension a(0:*)
character*40 form1,form2
character*(pmax) formt
data max/pmax/
nrec=(imax-imin)/10+1
do j=1,nrec
i1=(j-1)*10+imin ! Index of first variable of record
i2=i1+9 ! Index of last variable of record
if(i2.gt.imax) i2=imax ! " " " " "
formt(1:1)='('
index=2
do i=i1,i2
if(abs(a(i)).lt.athresh) then
formt(index:index+nchar1-1)=form1(1:nchar1)
index=index+nchar1
else
formt(index:index+nchar2-1)=form2(1:nchar2)
index=index+nchar2
endif
if(i.ne.i2) then
formt(index:index)=','
index=index+1
endif
end do
formt(index:index)=')'
index=index+1
C WRITE(6,9999) INDEX !!!
C9999 FORMAT(' INDEX = ',I3) !!!
write(nout,formt(1:index)) (a(i),i=i1,i2)
if(index.le.max) then
ifail=0
else
ifail=1
endif
end do
return
end
CBack up
subroutine writearray2(nout,a,imin,imax,nwid,ifail)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:18:11:1991 *
C * *
C * REVISION : Minor mod. to header JRH:06:01:1992 *
C * Addition of external JRH:15:05:2000 *
C * *
C * SOURCE : FORLIB.FVS *
C * ROUTINE NAME: WRITEARRAY2 *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Writes array A from indices IMIN to IMAX in "F" format, *
C * adjusting the format according the size of each value. *
C * 10 values are output per record. *
C * *
C * NOUT ..... Output device *
C * A ........ Array (NOTE that lower bound is 1) *
C * IMIN ..... Minimum index of array, A *
C * IMAX ..... Maximum index of array, A *
C * NWID ..... No. of characters to be occupied by each value *
C * (including one seperator space and one sign *
C * space) *
C * IFAIL .... 0 for successful execution, otherwise 1 *
C * *
C ******************************************************************************
integer*4 pmax
parameter(pmax=160)
real*4 a(*)
integer*4 nout,imin,imax,nwid,ifail
integer*4 nrec,j,i1,i2,index,i,idp
integer*4 int2,nleft
external int2
character*(pmax) formt
nrec=(imax-imin)/10+1
do j=1,nrec
i1=(j-1)*10+imin ! Index of first variable of record
i2=i1+9 ! Index of last variable of record
if(i2.gt.imax) i2=imax ! " " " " "
formt(1:1)='('
index=2
do i=i1,i2
if(a(i).ne.0.) then
nleft=int2(alog10(abs(a(i)))+1) ! No. of figs. to L of dec. pt.
else ! r = 0
nleft=1
endif
if(nleft.lt.0) nleft=0
idp=nwid-nleft-3 ! No. of decimal places
write(formt(index:index+6),1) nwid-1,idp
1 format('f',i2,'.',i2,'x')
index=index+7
if(i.ne.i2) then
formt(index:index)=','
index=index+1
endif
end do
formt(index:index)=')'
index=index+1
C WRITE(6,9999) INDEX !!!
C9999 FORMAT(' INDEX = ',I3) !!!
write(nout,formt(1:index)) (a(i),i=i1,i2)
if(index.le.pmax) then
ifail=0
else
ifail=1
endif
end do
return
end
CBack up
subroutine xy2ra(x,y,r,a,lr,la)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : UTILITY REF:JRH:02:07:1985 *
C * *
C * REVISION : ------------- JRH:--:--:1985 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: XY2RA *
C * TYPE : SUBROUTINE *
C * *
C * FUNCTION : Converts (X,Y) to (R,A(azimuth)).
C * *
C * LR = .TRUE. ..... Returns R *
C * LA = .TRUE. ..... Returns A *
C * *
C ******************************************************************************
logical lr,la
if(lr) then
r=sqrt(x*x+y*y)
else
r=0.
endif
if(la) then
if(x.ne.0..or.y.ne.0) then
a=atan2(x,y)
if(a.lt.0.) a=a+6.2831853
else
a=0.
endif
endif
return
end
CBack up
function ydif(x,xa,ya,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:31:07:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YDIF *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns first derivative of function Y(X) defined *
C * piecewise by arrays XA and YA, each containing N values. *
C * *
C * X ..... input X-value *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values *
C * N ..... number of values in each array *
C * *
C ******************************************************************************
dimension xa(*),ya(*)
if(x.lt.xa(1)) then ! X before start of data - use zero
ydif=0.
return
else
if(x.ge.xa(n)) then ! X after end of data - use zero
ydif=0.
return
else ! Linearly interpolate
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
ydif=(ya(i)-ya(i-1))/(xa(i)-xa(i-1))
return
endif
end do
endif
endif
end
CBack up
function ydif2d(x,iset,xa,ya,n,nmax)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:29:09:1986 *
C * *
C * REVISION : --------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YDIF2D *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns first derivative of function Y(X) defined *
C * piecewise by arrays XA and YA, each containing N values. *
C * *
C * X ..... input X-value *
C * ISET .. 2nd index of array YA *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values *
C * N ..... number of values in each array *
C * NMAX .. max. first index of array YA *
C * *
C * Special version for 2-D input array, YA *
C * *
C ******************************************************************************
dimension xa(*),ya(nmax,*)
if(x.lt.xa(1)) then ! X before start of data - use zero
ydif2d=0.
return
else
if(x.ge.xa(n)) then ! X after end of data - use zero
ydif2d=0.
return
else ! Linearly interpolate
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
ydif2d=(ya(i,iset)-ya(i-1,iset))/(xa(i)-xa(i-1))
return
endif
end do
endif
endif
end
CBack up
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
CBack up
function yint(x,xa,ya,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:30:07:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YINT *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns value of function Y(X) defined piecewise by *
C * arrays XA and YA, each containing N values, using linear *
C * interpolation where possible. *
C * *
C * X ..... input X-value *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values *
C * N ..... number of values in each array *
C * *
C ******************************************************************************
dimension xa(*),ya(*)
if(x.lt.xa(1)) then ! X before start of data - use first value
yint=ya(1)
return
else
if(x.ge.xa(n)) then ! X after end of data - use last value
yint=ya(n)
return
else ! Linearly interpolate
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
w1=(xa(i)-x)/(xa(i)-xa(i-1))
yint=w1*ya(i-1)+(1.-w1)*ya(i)
return
endif
end do
endif
endif
end
CBack up
function yint2d(x,iset,xa,ya,n,nmax)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:27:08:1986 *
C * *
C * REVISION : NMAX added JRH:29:09:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YINT2D *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns value of function Y(X) defined piecewise by *
C * arrays XA and YA, each containing N values, using linear *
C * interpolation where possible. *
C * *
C * X ..... input X-value *
C * ISET .. 2nd index of array YA *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values (first dimension = N) *
C * N ..... number of values in each array *
C * NMAX .. max. first index of array YA *
C * *
C * Special version for 2-D input array, YA *
C * *
C ******************************************************************************
dimension xa(*),ya(nmax,*) ! &29/9/86
if(x.lt.xa(1)) then ! X before start of data - use first value
yint2d=ya(1,iset)
return
else
if(x.ge.xa(n)) then ! X after end of data - use last value
yint2d=ya(n,iset)
return
else ! Linearly interpolate
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
w1=(xa(i)-x)/(xa(i)-xa(i-1))
yint2d=w1*ya(i-1,iset)+(1.-w1)*ya(i,iset)
return
endif
end do
endif
endif
end
CBack up
function yrect(x,xa,ya,n)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:22:09:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YRECT *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns value of function Y(X) defined as rectangular *
C * elements by arrays XA and YA, XA containing N values, *
C * YA containing N-1 values. *
C * *
C * X ..... input X-value *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values *
C * N ..... number of values in each array *
C * *
C ******************************************************************************
dimension xa(*),ya(*)
if(x.lt.xa(1)) then ! X before start of data - use first value
yrect=ya(1)
return
else
if(x.ge.xa(n)) then ! X after end of data - use last value
yrect=ya(n-1)
return
else ! Select rectangular element
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
yrect=ya(i-1)
return
endif
end do
endif
endif
end
CBack up
function yrect2d(x,iset,xa,ya,n,nmax)
C ******************************************************************************
C * *
C * FORTRAN SOURCE CODE *
C * *
C * PROGRAM SET : MODELLING REF:JRH:29:09:1986 *
C * *
C * REVISION : ------------- JRH:--:--:1986 *
C * *
C * SOURCE : FORLIB.FOR *
C * ROUTINE NAME: YRECT2D *
C * TYPE : FUNCTION *
C * *
C * FUNCTION : Returns value of function Y(X) defined as rectangular *
C * elements by arrays XA and YA, XA containing N values, *
C * YA containing N-1 values. *
C * *
C * X ..... input X-value *
C * ISET .. 2nd index of array YA *
C * XA .... array of X-values (these must monotonically *
C * increase with index) *
C * YA .... array of Y-values *
C * N ..... number of values in each array *
C * NMAX .. max. first index of array YA *
C * *
C * Special version for 2-D input array, YA *
C * *
C ******************************************************************************
dimension xa(*),ya(nmax,*)
if(x.lt.xa(1)) then ! X before start of data - use first value
yrect2d=ya(1,iset)
return
else
if(x.ge.xa(n)) then ! X after end of data - use last value
yrect2d=ya(n-1,iset)
return
else ! Select rectangular element
do i=2,n
if(x.lt.xa(i)) then ! X must be in interval I-1 to I
yrect2d=ya(i-1,iset)
return
endif
end do
endif
endif
end
C#
C
C
C