c Appendix 2 c concatenated profiles version of profile properties c Characterization method for profiles - get mean,stddev,max,min, c length and # peaks&troughs in window c gettings/Oct90; meg/Mar 93; meg/Jan95; meg May96(wavlet spectrum) c meg June 97 concatenated profiles version c meg 13 APR2000 fixed bugs, ywndow coord added;MEG 15May2000 pvals(6) c dimension x(25000),y(25000),yx(25000),zvals(6),xp(3),fld(11),pvals(6) dimension xwv(500),ywv(500), ywav(25000), ywin(25000) character ifmt*80,ofmt*80,bid*8,cid*8,oldid*8,iwfmt*80 character qu*80,qr*1,dflt1*1 logical verify real*8 rsum1,rsum2 data zvals/6*0./,xp/3*0./ ittin=5 ittout=6 io1=11 io2=12 io3=13 ianom=1 verify=.false. dflt1='n' qu='Verify keyboard responses (y/n)[n]:' call ttinaa(qu,qr,1,dflt1,ittin,ittout,verify) if (qr.eq.'y'.or.qr.eq.'Y') verify=.true. c c irec input record type c 1 binary lon,lat,z c 2 binary .post id,lon,lat,z(1-6) MGD FORMAT c 3 formatted lon,lat,z c 4 formatted .post (id,lon,lat,z(1-ianom) note variable ianom) c 5 formatted card image c id, lat deg, lat min, lon deg, lon min, z(1-ianom) c 6 binary .post MAD FORMAT c c jrec output record type where x replaces longitude and y replaces c latitude. range is jrec=1 to 4 (default is 2) c c ianom number of z anomalies to be read in a formatted input record. c default is 1. c c janom selection of one z anomaly for output in xyz record. c default is janom=2 (bouguer density 1 if gravity) c 360 write(*,'(a)')' Enter input file specs:' write(*,210) 210 format(' Input record type= 1 binary lon,lat,z'/ 1' 2 binary .post id,lon,lat,z(1-6) MGD FORMAT'/ 2' 3 formatted lon,lat,z'/ 3' 4 formatted .post (id,lon,lat,z(1-ianom) note variable ianom)'/ 4' 5 formatted card image'/ 5' id, lat deg, lat min, lon deg, lon min, z(1-ianom)'/ 6' 6 binary .post MAD FORMAT'/ 7' 7 formatted x,y'/ 8' Enter input record type:',$) read(*,'(i10)') irec if(irec.eq.4.or.irec.eq.5) then qu='Enter number of z values (<7)[1]' call ttini4(qu,ianom,1,ittin,ittout,verify) endif call opnfil(io1,ifmt,irec,0) write(*,'(a)')' Enter output file specs:' write(*,211) 211 format(' Output record type= 1 binary lon,lat,z'/ 1' 2 binary .post id,lon,lat,z(1-6) MGD FORMAT'/ 2' 3 formatted lon,lat,z'/ 3' 4 formatted .post (id,lon,lat,z(1-ianom) note variable ianom)'/ 4' 5 formatted x,yvals'/ 5' 6 binary .post id1,id2,lon,lat,z(1-6) MAD FORMAT'/ 6' Enter output record type:',$) read(*,'(i10)') jrec 370 call opnfil(io2,ofmt,jrec,1) 301 qu = ' Is this a profile dataset from program GETPRF (y/n)[n]?' dflt1 = 'n' ityp=1 call ttinaa(qu,qr,1,dflt1,ittin,ittout,verify) if(qr.eq.'y'.or.qr.eq.'Y') ityp=2 c c get field #s for x and y c call getfld(irec,ityp,nfldx,nfldy) c c loop on window width from nwstart down to nwstop points c ** nwind is the window width for each R/S estimate; nwind is c the lag tau; reducing nwind = nwind / 2 is used by Feder 1988; c 340 is the loop c which considers the data for each window width and outputs c a mean R/S for that width plus the standard deviation of mean R/S c ** nwdel is amount window of width nwind moves over on data for each c computation; for feder's case, nwdel = nwind (no overlap) c both parameters are set by : c nwind = nwind / nwden - nwoff; c the first window is always the whole dataset c nwdel = nwind * fracnw + nwdoff c For nwind = nwind / 2 set nwden = 2 and nwoff = 0 c For nwdel = nwind, set fracnw = 1.0 and nwdoff = 0; for nwdel = 1 c (maximum overlap) fracnw = 0.0 and nwdoff = 1, etc. c First get parameters nwden, nwoff, fracnw, and nwdoff c write(*,'(a,$)')'Enter largest window width: ' read(*,'(i10)') nwstart if(nwstart.eq.0) nwstart = npt write(*,'(a,$)')'Enter smallest (cutoff) window width: ' read(*,'(i10)') nwstop write(*,'(a)')'window width = width / nwd - nwoff for each loop thru data' write(*,'(a,$)')'Enter window divisor nwd and offset nwoff: ' read(*,'(2i10)') nwden, nwoff write(*,'(a)')'window movement = width * f + noff ' write(*,'(a,$)')'Enter factor f and offset noff: ' read(*,'(f10.0,i10)') fracnw, nwdoff c c Get wavelet if going to do spectral calcs; iwvlt=0 no fits; c =1 fit specified wavelet in each window c iwvlt = 0 write(*,'(a,$)')'Do you wish to fit a waveform to each window(y/n):' read(*,'(a)') qr if(qr.eq.'y'.or.qr.eq.'Y') iwvlt = 1 if(iwvlt.eq.1) then call opnfil(io3,iwfmt,3,0) i = 1 302 read(io3,iwfmt,end=300) xwv(i),ywv(i) i = i + 1 go to 302 300 nwvlt = i - 1 xwvrng = xwv(nwvlt) - xwv(1) xwvstrt = xwv(1) endif ifrst=0 ieof = 0 oldid='' c c 304 is loop over different profiles in data c 304 npt=0 if(ieof.eq.1) go to 550 nwind = (nwstart + nwoff) * nwden if(ifrst.eq.1.and.bid.ne.oldid) then npt=npt+1 x(npt)=fld(nfldx) y(npt)=fld(nfldy) oldid=bid endif c c read data c 305 call getrec(bid,fld(1),fld(2),zvals,xp,irec,ieof,ityp, 1 io1,ianom,ifmt) if(ieof.eq.1) go to 500 if(ifrst.eq.0) then oldid=bid ifrst=1 endif 306 do 50 j=1,6 fld(j+2)=zvals(j) if(j.gt.3) go to 50 fld(j+8)=xp(j) 50 continue if(bid.eq.oldid) then npt=npt+1 if (nfldx.ne.2) yx(npt) = fld(2) if (nfldx.eq.2) yx(npt) = fld(1) x(npt)=fld(nfldx) y(npt)=fld(nfldy) go to 305 endif 500 continue 340 nwind = nwind / nwden - nwoff nwindold = nwind rsum1=0.0d0 rsum2=0.0d0 numwin = 0 if (nwind.lt.nwstop) go to 304 nwdel = nwind * fracnw + nwdoff 318 npt1=1 npt2=nwind iquit=0 c c 320 is loop across data for all windows of width nwind c 320 call sstat1(y,npt1,npt2,pvals(1),pvals(2),pvals(3),pvals(4), 1 rrr,sss,rovers) deltax=x(npt2)-x(npt1) xwndow=x(npt1)+ deltax/2. deltay=yx(npt2)-yx(npt1) ywndow=yx(npt1)+ deltay/2. call npeak(x,y,npt1,npt2,npeek,ylen) pvals(5)=float(npeek)/abs(deltax) ylen = ylen / abs(deltax) if (npeek.eq.0) then c deltax = (x(npt2)-x(npt1))/2. pvals(6) = ylen / 2. else c deltax=(x(npt2)-x(npt1))/(2.*float(npeek)) pvals(6)=(ylen*abs(deltax))/(2.*float(npeek)) endif if(rovers.ne.0.0) then rsum1 = rsum1 + rovers rsum2 = rsum2 + rovers*rovers numwin = numwin + 1 endif c xp(1)=rrr xp(1)=deltax xp(2)=nwind xp(3)=rovers write(io2,226) xwndow, ywndow 226 format('xw,yw',2f10.3) call wrtrec(oldid,xwndow,ylen,pvals,xp,jrec,io2,6,ofmt) c c wavelet fit c if(iwvlt.eq.1) then do 70 i=1,nwind xx = (x(npt1 + i - 1) - x(npt1))/deltax*xwvrng + xwvstrt call interp(xwv,ywv,nwvlt,3,xx,ywav(i)) ywin(i) = y(npt1 + i - 1) 70 continue call linfit(ywav, ywin, ywin, nwind, 0, aicpt, siga, slop, sigb, rcor, sdfit) write(io2,205) nwind,x(npt2),x(npt1),aicpt, siga, slop, sigb, rcor, sdfit 205 format('LSF wvlt:',i13,1p8e13.6) endif if(npt2.eq.npt.or.iquit.eq.1) go to 330 npt1=npt1+nwdel if(npt1.gt.npt) go to 330 npt2=npt1+nwind-1 if(npt2.gt.npt) then npt2=npt nwind = npt2 - npt1 + 1 iquit=1 endif if((npt2-npt1).lt.9) go to 330 go to 320 c c end 320 loop c 330 continue nwind = nwindold c c calc mean and dispersion of rovers for this window width c roversmean = rsum1 / float(numwin) if(numwin.le.1) then roversstd = 0.0 else roversstd = rsum2*float(numwin) - (rsum1*rsum1) roversstd = sqrt( roversstd/(float(numwin)*float(numwin-1))) endif write(io2,200) nwind, numwin, roversmean, roversstd 200 format('Mean R/S:',2i13,1p2e13.6) go to 340 c c end 340 loop c c go to next profile go to 304 550 close(io1) stop end subroutine sstat1(x,npt1,npt2,xmn,xmx,xbar,xsd, r, s, rovers) c c find max,min,mean, standard deviation, range and r/s c of array x between npt1 & npt2 c rewritten Mar 93 / MEG c real*8 s1, s2 dimension x(*) xmx=-1e36 xmn=1e36 s1=0.0 s2=0.0 do 60 i=npt1,npt2 xi=x(i) s1=s1+xi c c calculate std dev in second pass to prevent loss of significance c when squaring data with large numbers e.g. mag field of 49500 gammas c c s2=s2+xi*xi if (xi.lt.xmn) xmn=xi if (xi.gt.xmx) xmx=xi 60 continue c c calc mean, std. dev, range and r/s c fnpt=float(npt2-npt1+1) xbar=s1/fnpt c s2 = s2 - fnpt*xbar*xbar do 65 i=npt1,npt2 xi = x(i) - xbar s2=s2+xi*xi 65 continue if((fnpt-1.).eq.0.0) then xsd = sqrt(s2/fnpt) else xsd=sqrt(s2/(fnpt-1)) endif r=xmx-xmn s=xsd if(xsd.eq.0.0) then rovers = 0.0 else rovers=r/xsd endif return end subroutine norm(x,npt,xmx,xmn) c c normalize an array to range xmx-xmn c dimension x(npt) f=1./(xmx-xmn) do 50 i=1,npt x(i)=(x(i)-xmn)*f 50 continue return end subroutine opnfil(iou,fmt,irtyp,k) c c Open a file and get format if formatted c iou LUN number c fmt format (<80 char) c irtyp file type: 1,2, or 6 - binary c 3-5 or 7 - formatted c k - 0 open existing file c - 1 open new file c MEG/JUL 85 c character fmt*80, fnam*50 write(*,'(a$)')' Filename: ' read(*,'(a)') fnam if(irtyp.ge.3 .and. irtyp.le.5) go to 300 if(irtyp.eq.7) go to 300 c c binary file c fmt=' ' if(k.eq.1) go to 305 open(unit=iou,file=fnam,status='old',form='unformatted') go to 900 305 open(unit=iou,file=fnam,status='unknown',form='unformatted') go to 900 c c formatted file c 300 write(*,'(a)')' Enter format including parentheses (<81c)' read(*,'(a)') fmt if(k.eq.1) go to 315 open(unit=iou,file=fnam,status='old') go to 900 315 open(unit=iou,file=fnam,status='unknown') 900 return end subroutine getrec(bid,xd,yd,b,x,irec,ieof,ityp,ior,nb,fmt) c c read a record in various geophysics formats MEG/JUL 85 c bid - 8c id c xd - lon in degrees c yd - lat in degrees c b - 6 element array of data values c x - 3 element array of x,y,s for profile data c x,y are rect. coords; s is distance down profile c irec- identifies how to read; same as prjrec program c ieof- 0 = not EOF; 1 = EOF sensed c ityp- 1 not a profile record (don't read x,y,s) c 2 profile record (read x,y,s) c ior - LUN for read c nb - number of data values (<7) to read into b c fmt - format for formatted read c character bid*8,id(2)*4,fmt*80 dimension b(6),x(3) do 50 i=1,6 50 b(i)=0.0 ieof=0 if(ityp.eq.2) go to 300 c c Non-profile record c go to (21,22,23,24,25,26,27),irec 21 read(ior,end=10)xd,yd,b(1) go to 900 22 read(ior,end=10) bid,xd,yd,b go to 900 23 read(ior,fmt,end=10) xd,yd,b(1) go to 900 24 read(ior,fmt,end=10)bid,xd,yd,(b(i),i=1,nb) go to 900 25 read(ior,fmt,end=10)bid,yd,ym,xd,xm,(b(i),i=1,nb) yd=yd+ym*1.666667e-2 xd=abs(xd)+xm*1.666667e-2 go to 900 26 read(ior,end=10) id,xd,yd,b bid=id(1)//id(2) go to 900 27 read(ior,fmt,end=10) xd, yd go to 900 c c Profile data read c 300 continue go to (121,122,123,124,125,126,127),irec 121 read(ior,end=10)xd,yd,b(1),x go to 900 122 read(ior,end=10) bid,xd,yd,b,x go to 900 123 read(ior,fmt,end=10) xd,yd,b(1),x go to 900 124 read(ior,fmt,end=10)bid,xd,yd,(b(i),i=1,nb),x go to 900 125 read(ior,fmt,end=10)bid,yd,ym,xd,xm,(b(i),i=1,nb),x yd=yd+ym*1.666667e-2 xd=abs(xd)+xm*1.666667e-2 go to 900 126 read(ior,end=10) id,xd,yd,b,x bid=id(1)//id(2) go to 900 127 read(ior,fmt,end=10) xd,yd,x 900 return 10 ieof=1 go to 900 end subroutine wrtrec(bid,xd,yd,b,x,jrec,iow,ia,fmt) c c write profile data record in various formats MEG/JUL 85 c bid - 8c id c xd - lon deg c yd - lat deg c b - 6 element data array c x - 3 element x,y,s profile coords; s= dist dn prof. c jrec- integer defining type of write, similar pgm prjrec c iow - LUN for write c ia - number of data elements (<7) to write from b c fmt - format for formatted write c character fmt*80,bid*8,cid*8 character*4 id(2) dimension b(6),x(3) equivalence (id(1),cid) cid=bid go to (3,4,5,6,5,8),jrec 3 write(iow)xd,yd,b,x go to 7 4 write(iow)bid,xd,yd,b,x go to 7 5 write(iow,fmt)xd,yd,b,x go to 7 6 write(iow,fmt) bid,xd,yd,(b(i),i=1,ia),x go to 7 8 write(iow) id(1),id(2),xd,yd,b,x 7 return end subroutine getfld(irec,ityp,nx,ny) c c Enter x and y field numbers for various input record types c if (ityp.eq.2) go to 300 c Not a profile record from GETPRF go to (10,20,10,40,60,20,50),irec 10 write(*,200) 200 format(' Fields are: 1) lon, 2) lat, 3) z.') go to 99 20 write(*,201) 201 format(' Fields are: 1) lon, 2) lat, 3-8) z(i=1-6).') go to 99 40 write(*,202) 202 format(' Fields are: 1) lon, 2) lat, 3-ianom+2) z(i=1-ianom).') go to 99 50 write(*,203) 203 format(' Fields are: 1) x, 2) y.') go to 99 60 write(*,204) 204 format(' Fields are: 1) lat, 2) lon, 3-ianom+2) z(i=1-ianom).') go to 99 c GETPRF profile record 300 go to (110,120,110,140,145,120,150),irec 110 write(*,205) 205 format(' Fields are: 1) lon, 2) lat, 3) z, 9) x, 10) y, 11) s.') go to 99 120 write(*,206) 206 format(' Fields are: 1) lon, 2) lat, 3-8) z(i=1-6), 9) x, 10)', 1 ' y, 11) s.') go to 99 140 write(*,207) 207 format(' Fields are: 1) lon, 2) lat, 3-ianom+2) z(i=1-ianom)', 1 ', 9) x, 10) y, 11) s.') go to 99 145 write(*,208) 208 format(' Fields are: 1) lat, 2) lon, 3-ianom+2) z(i=1-ianom)', 1 ', 9) x, 10) y, 11) s.') go to 99 150 write(*,209) 209 format(' Fields are: 1) xdata, 2) ydata, 9) x, 10) y, 11) s.') 99 write(*,'(a$)')' Enter field numbers for variables x and y:' read(*,'(2i10)') nx,ny return end subroutine npeak(x,y,nst,nend,npks,s) c calculate the number of peaks and troughs of y as fn of x c and the length of straight line segs connecting {x,y} c MEG/Mar 88. real*8 sum dimension x(*),y(*) npks=0 sum=sqrt((x(nst+1)-x(nst))**2+(y(nst+1)-y(nst))**2) s2=(y(nst+1)-y(nst))/(x(nst+1)-x(nst)) do 50 i=nst+1,nend-1 sum=sum+sqrt((x(i+1)-x(i))**2+(y(i+1)-y(i))**2) s1=s2 s2=(y(i+1)-y(i))/(x(i+1)-x(i)) if(sign(1.0,s1).ne.sign(1.0,s2)) npks=npks+1 50 continue s=sum return end subroutine shstat(x,npt,xmn,xmx,xbar,xsd) c c find max,min,mean and standard deviation of array x c real*8 s1 dimension x(npt) xmx=-1e36 xmn=1e36 s1=0.0 do 60 i=1,npt xi=x(i) s1=s1+xi if (xi.lt.xmn) xmn=xi if (xi.gt.xmx) xmx=xi 60 continue fnpt=float(npt) xbar=s1/fnpt s1=0.0 do 65 i=1,npt t1=x(i)-xbar s1=s1+t1*t1 65 continue xsd=sqrt(s1/(fnpt-1)) return end c c purpose c interpolate between data points to evaluate a function c c usage c call interp(x, y, npts, nterms, xin, yout) c c description of parameters c x - array of data points for independent variable c y - array of data points for dependent variable c npts - number of pairs of data points c (integer*4) c nterms - number of terms in fitting polynomial c (integer*4) c xin - input value of x c yout - interpolated value of y c c subroutines and function subprograms required c none c c modifications for fortran ii c omit double precision specifications c c comments c dimension statement valid for nterms up to 10 c value of nterms may be modified by the program c compilation with /i4 on pdp11/45 .......... m.d., usgs, dec 80 c c routine from Bevington, Data reduction and error analysis for the c physical sciences c Disclaimer: Although this program has been tested, the U.S.Geological c Survey cannot guarantee that it will give correct results nor that it c will work on all computer systems. c subroutine interp (x, y, npts, nterms, xin, yout) double precision deltax, delta, a, prod, sum dimension x(1), y(1) dimension delta(10), a(10) c c search for appropriate value of x(1) c 11 do 19 i=1, npts if (xin - x(i)) 13, 17, 19 13 i1 = i - nterms/2 if (i1) 15, 15, 21 15 i1 = 1 go to 21 17 yout = y(i) 18 go to 61 19 continue i1 = npts - nterms + 1 21 i2 = i1 + nterms - 1 if (npts - i2) 23, 31, 31 23 i2 = npts i1 = i2 - nterms + 1 25 if (i1) 26, 26, 31 26 i1 = 1 27 nterms = i2 - i1 + 1 c c evaluate deviations delta c 31 denom = x(i1+1) - x(i1) deltax = (xin - x(i1)) / denom do 35 i=1, nterms ix = i1 + i -1 35 delta(i) = (x(ix) - x(i1)) / denom c c accumulate coefficients a c 40 a(1) = y(i1) 41 do 50 k=2, nterms prod = 1. sum = 0. imax = k - 1 ixmax = i1 + imax do 49 i=1, imax j = k - i prod = prod * (delta(k) - delta(j)) 49 sum = sum - a(j)/prod 50 a(k) = sum + y(ixmax)/prod c c accumulate sum of expansion c 51 sum = a(1) do 57 j=2, nterms prod = 1. imax = j - 1 do 56 i=1, imax 56 prod = prod * (deltax - delta(i)) 57 sum = sum + a(j)*prod 60 yout = sum 61 return end subroutine linfit (x,y,sigmay,npts,mode,a,sigmaa,b,sigmab,r, 1stdv) c subroutine linfit c c purpose c make a least squares fit to data with a straight line c y=a+b*x c c usage c call linfit (x,y,sigmay,npts,mode,a,sigmaa,b,sigmab,r,stdv) c c description of parameters c x -array of data points for independent variable c y -array of data points for dependent variable c sigmay -array of standard deviations for y data points c npts -number of pairs of data points c mode -determines method of weighting least squares fit c +1 (instrumental) weight(i)=1./sigmay(i)**2 c 0 (no weight) weight(i)=1. c -1 (statistical) weight(i)=1./y(i) c a -y intercept of fitted straight line c sigmaa -standard deviation of a c b -slope of fitted stright line c sigmab -standard deviation of b c r -linear correlation coefficient c stdv -standard deviation of fitted straight line c subroutines and function subprograms required c none c c modifications for fortran ii c omit double precision specifications c change dsqrt to sqrt in statements 67, 68, and 71 real*8 sum,sumx,sumy,sumx2,sumxy,sumy2 real*8 xi,yi,weight,delta,varnce,arg1,da,db dimension x(1),y(1),sigmay(1) c c accumulate weighted sums c 11 sum=0. sumx=0. sumy=0. sumx2=0. sumxy=0. sumy2=0. 21 do 50 i=1,npts xi=x(i) yi=y(i) if (mode) 31, 36, 38 31 if (yi) 34, 36, 32 32 weight=1./yi go to 41 34 weight=1./(-yi) go to 41 36 weight=1. go to 41 38 weight=1./sigmay(i)**2 41 sum=sum+weight sumx=sumx+weight*xi sumy=sumy+weight*yi sumx2=sumx2+weight*xi*xi sumxy=sumxy+weight*xi*yi sumy2=sumy2+weight*yi*yi 50 continue c c calculate coefficients and standard deviations c 51 delta=sum*sumx2-sumx*sumx da=(sumx2*sumy-sumx*sumxy)/delta a=da 53 db=(sumxy*sum-sumx*sumy)/delta b=db 61 if (mode) 62, 64, 62 62 varnce=1. go to 67 64 c=npts-2 varnce=(sumy2+da*da*sum+db*db*sumx2 1-2.*(da*sumy+db*sumxy-da*db*sumx))/c if(varnce.lt.0.0)write(6,203) varnce 203 format(' ************** varnce lt 0:',e14.7) stdv=dsqrt(varnce) 67 arg1=varnce*sumx2/delta if(arg1.lt.0.0) write(6,200)arg1 200 format(' **************arg in sigmaa lt 0:',e14.7) sigmaa=dsqrt(arg1) 68 arg1=varnce*sum/delta if(arg1.lt.0.0) write(6,201)arg1 201 format(' ***************** arg for sigmab lt 0:',e14.7) sigmab=dsqrt(arg1) 71 arg1=delta*(sum*sumy2-sumy*sumy) if(arg1.lt.0.0)write(6,202)arg1 202 format(' **************arg lt 0 in r:',e14.7) r=(sum*sumxy-sumx*sumy)/dsqrt(arg1) return end subroutine ttinr4(qu, r4, defalt, itin,itout, verify) c Terminal entry of real*4 variable - F77 version c qu - Character string prompt que,80 C long. c r4 - variable to be input c defalt - default value if only given c itin,itout - read, write LUN's of terminal c verify - if true, asks for verification of input c if false no verification c typical call: c call ttinr4(que,a,1.67,itin,itout,verify) c MEG/ Feb 86; SUN conversion Jul 89/ MEG c logical verify character ans*10, blank*10, qu*80 character que(80)*1, qr*1 real r4, defalt blank = ' ' read(unit=qu, fmt=200) que 200 format(80a1) call lnnobl(que, nq, 80) 300 write(unit=itout, fmt=204) 204 format(1h ,$) do 40 i = 1, nq 40 write(unit=itout, fmt=201) que(i) 201 format(a,$) read(unit=itin, fmt=101, err=300) ans 101 format(a) if (ans .eq. blank) goto 910 read(unit=ans, fmt=100) r4 100 format(bn,f10.0) if (.not. verify) goto 900 write(unit=itout, fmt=202) r4 202 format(1x,6hValue ,1pe14.7,5h ok? ,$) read(unit=itin, fmt=101) qr if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300 900 return 910 r4 = defalt goto 900 end subroutine lnnobl(a, na, lngth) c determine length of string of max length lngth with no blanks on left c changed to specify string max length from 80 by MEG/10Jan91 character a(lngth), bl bl = ' ' do 50 i = lngth, 1, -1 na = i if (a(i) .ne. bl) goto 60 50 continue na = 0 60 return end subroutine ttinr8(qu, r8, defalt, itin,itout, verify) c Terminal entry of real*8 variable - F77 version c qu - Character string prompt que,80 C long. c r8 - variable to be input c def - default value if only given c itin,itout - read, write LUN's of terminal c verify - if true, asks for verification of input c if false no verification c typical call: c call ttinr8(que,a,1.67d-03,itin,itout,verify) c MEG/ Feb 86; SUN conversion Jul 89/ MEG c logical verify character ans*16, blank*16, qu*80 character que(80)*1, qr*1 double precision r8, defalt blank = ' ' read(unit=qu, fmt=200) que 200 format(80a1) call lnnobl(que, nq, 80) 300 write(unit=itout, fmt=204) 204 format(1h ,$) do 40 i = 1, nq 40 write(unit=itout, fmt=201) que(i) 201 format(a,$) read(unit=itin, fmt=101, err=300) ans 101 format(a) if (ans .eq. blank) goto 910 read(unit=ans, fmt=100) r8 100 format(bn,f16.0) if (.not. verify) goto 900 write(unit=itout, fmt=202) r8 202 format(1x,6hValue ,1pd24.17,5h ok? ,$) read(unit=itin, fmt=101) qr if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300 900 return 910 r8 = defalt goto 900 end subroutine ttini4(qu, i4, defalt, itin,itout, verify) c Terminal entry of integer*4 variable - F77 version c qu - Character string prompt que,80 C long. c i4 - variable to be input c defalt - default value if only given c itin,itout - read, write LUN's of terminal c verify - if true, asks for verification of input c if false no verification c typical call: c call ttini4(que,k,16,itin,itout,verify) c MEG/ Feb 86; SUN conversion Jul 89/ MEG c logical verify character ans*16, blank*16, qu*80 character que(80)*1, qr*1 integer*4 i4, defalt blank = ' ' read(unit=qu, fmt=200) que 200 format(80a1) call lnnobl(que, nq, 80) 300 write(unit=itout, fmt=204) 204 format(1h ,$) do 40 i = 1, nq 40 write(unit=itout, fmt=201) que(i) 201 format(a,$) read(unit=itin, fmt=101, err=300) ans 101 format(a) if (ans .eq. blank) goto 910 read(unit=ans, fmt=100) i4 100 format(bn,i16) if (.not. verify) goto 900 write(unit=itout, fmt=202) i4 202 format(1x,6hValue ,i16,5h ok? ,$) read(unit=itin, fmt=101) qr if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300 900 return 910 i4 = defalt goto 900 end subroutine ttini2(qu, i2, defalt, itin,itout, verify) c Terminal entry of integer*2 variable - F77 version c qu - Character string prompt que,80 C long. c i2 - variable to be input c defalt - default value if only given c itin,itout - read, write LUN's of terminal c verify - if true, asks for verification of input c if false no verification c typical call: c call ttini2(que,k,16,itin,itout,verify) c MEG/ Feb 86; SUN conversion Jul 89/ MEG c logical verify character ans*6, blank*6, qu*80 character que(80)*1, qr*1 integer*2 i2, defalt blank = ' ' read(unit=qu, fmt=200) que 200 format(80a1) call lnnobl(que, nq, 80) 300 write(unit=itout, fmt=204) 204 format(1h ,$) do 40 i = 1, nq 40 write(unit=itout, fmt=201) que(i) 201 format(a,$) read(unit=itin, fmt=101, err=300) ans 101 format(a) if (ans .eq. blank) goto 910 read(unit=ans, fmt=100) i2 100 format(bn,i6) if (.not. verify) goto 900 write(unit=itout, fmt=202) i2 202 format(1x,6hValue ,i6,5h ok? ,$) read(unit=itin, fmt=101) qr if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300 900 return 910 i2 = defalt goto 900 end subroutine ttinaa(qu, cs, ncs, defalt, itin,itout, verify) c Terminal entry of character string - F77 version c qu - 80C string of prompt que c cs - character string to be input-array of c*1 but can c be equivalenced to char string in calling pgm c ncs - # characters in cs c defalt - default string for cs c itin,itout - read, write LUN's of terminal c verify - if true, asks for verification of input c if false no verification c typical call: c call ttinaa(que,fnam,16,defalt,itin,itout,verify) c MEG/ Dec 85; SUN conversion Jul 89/ MEG c logical verify character css(80) character qu*80 character que(80)*1, qr*1 c character defalt(*), cs(*) replaced with following: character*(*) cs,defalt read(unit=qu, fmt=200) que 200 format(80a1) call lnnobl(que, nq, 80) 300 write(unit=itout, fmt=204) 204 format(' ',$) do 40 i = 1, nq 40 write(unit=itout, fmt=201) que(i) 201 format(a1,$) read(unit=itin, fmt=101, err=300) css 101 format(80a1) call lnnobl(css, kcss, 80) c go to default if (kcss .eq. 0) goto 910 write(cs,'(80a)')(css(i),i=1,ncs) if (.not. verify) goto 900 write(unit=itout, fmt=202) (css(i),i = 1, ncs) 202 format(1x,17hCharacter string /1x,80a:) write(unit=itout, fmt=203) 203 format(5h ok? ,$) read(unit=itin, fmt=101) qr if ((qr .eq. 'n') .or. (qr .eq. 'N')) goto 300 900 return 910 cs = defalt goto 900 end