parameter(nndim=5000,ipow=4096,icol=20) data io/3/,jo/4/ itype=-1 call gfs_opena(io,'input file : ',itype,iret) if(iret.gt.0.and.itype.eq.1) goto 5 print *,'file does not exist or is of wrong type' goto 50 5 call gfs_opena(jo,'output file : ',itype,iret) c** initialise for first file *** call gfs_phdef(io) call gfs_sdef(io) call sel3(io,jo) call gfs_close(jo) 50 call gfs_close(io) stop end c c subroutine sel3(io,jo) c The main purpose of select is to allow data sequences to be extracted from c time series. There is a zoom feature to aid selection-the original version c stored all zoom levels and scrolled back through a nested zoom sequence. c this was usually found to be irritating so this version saves only the c original screen and the last zoom level. parameter(nndim=5000,ipow=4096,icol=20,no=7) integer*4 spoint,epoint,xlo,xhi,indata,xl character*256 string logical autosc,nxtf,log1,log2 common/scayle/zmin,zmax common/sdata/comp(nndim,3) common/slect$/autosc common/head$/nsta,nchn,iy,id,ih,im,ss,si,nscan,title common/flag$$/flag,badun common/env$$$/spoint,epoint,nser,npt,ymax,ymin,screen(1022,3) + ,il,ix common/scrtch/scrbuf(nndim) common/myhed/jscan,jsta,jchn,jtyp,jy,jd,jh,jm,ssj,dtj,spare(20) common/trav/trminr(no),trmaxr(no),trminl(no),trmaxl(no) dimension hbuf(30),fnm(20) character*32 text(12) character*80 title,ffiler,ffilet,fint character*120 line character*1 edc common/hcopy/ncops character*9 namhcopy equivalence (jscan,hbuf) data dflag,badf/5.e15,0.7/ ffiler='disp.sph' ffilet='disp.tor' c... set numbers for polar pi=4.*atan(1.d0) ndim=nndim print*,' enter # windows and time b-w product (4 3)' read(*,*)nwin,fw print*,' enter unwrapping parameter for hor. angles' print*,' 0:no, 1:rayl (-90,90),2:love (0,180), 3:(-360,360)' read(*,*)iwave print*,' flip sign of comp-3?' read(*,*)iflip iflip=1-2*iflip print*,' enter frequency range (in Hz)' read(*,*)fmin,fmax print*,' enter # points for 2D matrix' read(*,*)np2d text(1)=' choose range to view ' text(2)=' select ranges ' text(3)=' time a point ' text(4)=' A) set scale or B) self scale ' text(5)=' A) next+plot B) next C) skip ' text(6)=' quit ' text(7)=' A) advance or B) re-plot ' text(8)=' set flag ' text(9)=' hardcopy ' text(10)=' change model for arrivals ' text(11)=' calculate 2D matrix of d1 ' text(12)=' horiz. part. motion plot ' flag=dflag badun=badf*flag autosc=.true. ymin=0. ymax=0. ncops=0 lenv=0 kfl=0 edc='A' c... set initial menu c... (print title of seismogram) c... (get source, calculate distance and begin of record c... rel. to origin time, get dispersion model) c... (labels 35,351,352,15,152) call scr_init(text,12) c... now read the first file goto 35 c... get data, plot them, plot expected arrivals and wait for mouse input c... (label 152,1000,200,subroutine getran,label 10) 1000 if(autosc) ymax=ymin call scr_p3(0,0,screen,npt,spoint,nser,ymin,ymax,flag,title) write(line,'(a,f6.2,a,f8.2)') & ' dist ',distg,' azi ',azi call scr_pline(line) xmina=(spoint-1)*dtj/60. xmaxa=(spoint+nser-1)*dtj/60. dsam=60./dtj call axis(xmina,xmaxa,zmin,zmax,dsam,1) call plt_trav(spoint) c... you need to activate and deactivate segment to get c... everything on hardcopy goto 200 c*** get option **** 5 call scr_xhair(xlo,xval,yval,edc,iopt) if(iopt.gt.0) goto 10 call scr_pline('please choose an option below') goto 5 10 jopt=iopt goto (15,20,25,30,35,40,45,50,55,65,75,85),iopt c... option 1: choose range to view 15 if(edc.eq.'C') then spoint=1 nser=nscan elseif(edc.eq.'A')then 151 call scr_rdnums('input range to view :$',fnm,num) if(num.eq.0) goto 152 if(num.ne.2) goto 151 if(fnm(1)+1.ge.fnm(2)) goto 151 i2=nint(fnm(2)) i1=nint(fnm(1)) nser=i2-i1+1 spoint=i1 endif c... read gfs data into buffer 152 npts=indata(io,kfl,spoint,nser) if(npts.gt.1) goto 1000 if(spoint.eq.1.and.nser.gt.1) goto 40 call scr_pline('invalid range-try again') goto 15 c... option 2: select range 20 if(edc.ne.'C') goto 200 call scr_pline(' mouse key A or B = extract ranges') goto 5 200 call scr_pline('mouse key options with the x-hair are :') call scr_pline('A=select range ; B=zoom ; C=unzoom or quit') call getran(xlo,xhi,edc,iopt) if(iopt.ne.0) goto 10 if(edc.eq.'C') then * ... unzoom if(lenv.ne.0)goto 206 if(nser.le.0)goto 200 lenv=0 goto 15 elseif(edc.eq.'B')then c*** zoom in region of plot *** len=xhi-xlo+1 call pshenv(lenv) npts=indata(io,kfl,xlo,len) goto 1000 endif if(xhi-xlo.gt.0.and.xhi-xlo.le.5000) goto 202 call scr_pline('selected range invalid-try again') goto 200 c**************************************************************** c... MTPA c... calling POLAR and plotting results with PLTRES c... input of parameters have been changed to interactive c... I/O c... call of POLAR fits to updated subroutines in c... /net/guppy/gabi/polaris/polsub?.o (?=1,2) c... subroutine EPTIS evaluates ellipticity and unwrapps angles c... if option IWAVE is set 202 kscan=xhi-xlo+1 do 666 jser=1,3 mfl=kfl+jser-1 666 call gfs_rwdata(io,comp(1,jser),mfl,xlo,kscan,'r',ierr) fmax=amin1(fmax,1./(2.*si)) do 667 js=1,kscan 667 comp(js,3)=comp(js,3)*float(iflip) call polar(ndim,kscan,si,nwin,fw,fmin,fmax,nout) call eptis(iwave,nout) call pltres(nout,si,fmin,fmax,iwave,xlo,nwin,fw,kscan,0, & azi,distg,jo,title) call scr_init(text,12) goto 1000 c************************************************************************ c*** unzoom *** 206 if(lenv.eq.0) goto 5 call popenv(lenv) goto 1000 c... option 3: time a point 25 call scr_xhair(xl,xval,yval,edc,iopt) if(iopt.ne.0) goto 10 if(edc.eq.'C') goto 206 call updat(xl,iy1,id1,ih1,im1,ss1) write(line,250) xl,iy1,id1,ih1,im1,ss1,yval 250 format('time of point ',i10,' is : ',i5,i4,2i3,f6.2, + ' y value is :',g14.6,'$') call scr_pline(line) goto 25 c... option 4: choose type of scaling 30 if(edc.eq.'B') goto 302 if(edc.eq.'A') goto 301 call scr_pline(' mouse key A = input scale from screen') call scr_pline(' mouse key B = default to self-scale') goto 5 301 call scr_rdnums('min and max (return for self-scale) :$',fnm,num) if(num.eq.0) goto 302 if(num.ne.2) goto 301 if(fnm(1).ge.fnm(2)) goto 301 ymin=fnm(1) ymax=fnm(2) autosc=.false. goto 1000 302 ymin=0. ymax=0. autosc=.true. goto 1000 c... option 5: choose time series 35 if(edc.eq.'A'.or.edc.eq.'B') goto 351 350 call scr_rdnums('enter time series number :',fnm,num) if(num.ne.1) goto 350 kfl=fnm(1)-1 kfl=max0(0,kfl) 351 if(nxtf(io,kfl)) goto 352 call scr_pline('unable to move to requested file') goto 5 352 call scr_pline(title) call dist(dista,azi,tstart,fmin,fmax,ffiler,ffilet) distg=dista/6371./pi*180. do 701 i=1,no trminr(i)=trminr(i)/dtj trmaxr(i)=trmaxr(i)/dtj trminl(i)=trminl(i)/dtj trmaxl(i)=trmaxl(i)/dtj c print 1001,trminr(i),trmaxr(i),trminl(i),trmaxl(i) c1001 format(4f20.5) 701 continue lenv=0 if (edc.eq.'A') goto 152 goto 5 c... option 6: quit interactive graphics 40 call scr_end return c... option 7: advance in time series or re-plot 45 if(edc.eq.'A') goto 452 if(edc.eq.'B') goto 451 call scr_pline(' mouse key A = plot next segment of data') call scr_pline(' mouse key B = re-plot this segment of data') goto 5 451 npts=indata(io,kfl,spoint,nser) goto 453 452 if(spoint+nser.ge.nscan) goto 454 npts=indata(io,kfl,spoint+nser,nser) 453 if(npts.gt.1) goto 1000 454 call scr_pline('end of file reached') goto 5 c... option 8: set flag 50 call scr_rdnums('enter flag value (return for default) :',fnm,num) flag=dflag badun=badf*flag if(num.eq.0) goto 501 if(num.ne.1) goto 50 flag=fnm(1) if(flag.le.0.) goto 50 badun=badf*flag 501 write(line,502) flag 502 format('flag value is ',g12.4,'$') call scr_pline(line) goto 5 c... option 9: hardcopy 55 ncops=ncops+1 do 56 jjj=1,9 56 namhcopy(jjj:jjj)=' ' if(ncops.le.9) then write(namhcopy,57) ncops elseif(ncops.ge.10.and.ncops.le.99) then write(namhcopy,58) ncops else write(namhcopy,59) ncops end if 57 format('hrdcpy',i1) 58 format('hrdcpy',i2) 59 format('hrdcpy',i3) call scr_rdtext('enter output name:',string,nstrin) call scr_dump(7,string) goto 5 c... option 10: change model for arrival times 65 continue log1=.true. fint=' ' call scr_rdtext('enter file of Rw disp. c./dummy$',fint,nfr) if(fint.eq.'dummy')then do 702 i=1,no trminr(i)=1 trmaxr(i)=1 trminl(i)=1 trmaxl(i)=1 702 continue goto 1000 endif ffiler=fint inquire(file=ffiler,exist=log1) if(.not.log1)then call scr_pline('file not existent, enter again') goto 65 endif 66 continue log2=.true. ffilet=' ' call scr_rdtext('enter file of Lw disp. c.$',ffilet,nft) inquire(file=ffilet,exist=log2) if(.not.log2)then call scr_pline('file not existent, enter again') goto 66 endif call dist(dista,azi,tstart,fmin,fmax,ffiler,ffilet) distg=dista/6371./pi*180. do 700 i=1,no trminr(i)=trminr(i)/dtj trmaxr(i)=trmaxr(i)/dtj trminl(i)=trminl(i)/dtj trmaxl(i)=trmaxl(i)/dtj 700 continue goto 1000 c... option 11: show 2D plot of largest singular value 75 continue call scr_pline(' ') call scr_pline('mouse key opitions with the x-hair are:') call scr_pline('A=select range; C=unzoom or quit') call getran(xlo,xhi,edc,iopt) if(edc.eq.'C')goto 206 if(edc.eq.'B')goto 75 if(xhi-xlo.gt.0.and.xhi-xlo.le.nndim)goto 601 call scr_pline('selected range invalid-try again') goto 75 601 kscan=np2d c... decrease istep, if you want to have a denser plot on c... y-axis istep=int(float(xhi-xlo)/20.) if(istep.eq.0)then call scr_pline('chosen window too small, choose again') goto 75 endif ires=0 call plt_mat(io,kfl,xlo,xhi,kscan,nwin,fw, & fmin,fmax,istep,ires,nout,iwave,iflip) c... plot results, if travel times were picked in 2D matrix c... of d1 if(ires.eq.1)then call eptis(iwave,nout) call pltres(nout,si,fmin,fmax,iwave,xlo,nwin,fw,kscan,1, & azi,distg,jo,title) endif call scr_init(text,12) goto 1000 c... option 12: particle motion plot 85 continue call scr_pline(' ') call scr_pline('mouse key opitions with the x-hair are:') call scr_pline('A=select range; C=unzoom or quit') call getran(xlo,xhi,edc,iopt) if(edc.eq.'C')goto 206 if(edc.eq.'B')goto 85 if(xhi-xlo.gt.0.and.xhi-xlo.le.nndim)goto 86 call scr_pline('selected range invalid-try again') goto 85 86 continue kscan=xhi-xlo+1 do 87 jser=1,3 mfl=kfl+jser-1 87 call gfs_rwdata(io,comp(1,jser),mfl,xlo,kscan,'r',ierr) xpmit=0. ypmit=0. xpmax=0. do 88 jjj=1,kscan xpmit=xpmit+comp(jjj,3) ypmit=ypmit+comp(jjj,2) 88 continue xpmit=xpmit/kscan ypmit=ypmit/kscan do 89 jjj=1,kscan comp(jjj,3)=(comp(jjj,3)-xpmit) comp(jjj,2)=(comp(jjj,2)-ypmit) if(abs(comp(jjj,3)).gt.xpmax)xpmax=abs(comp(jjj,3)) if(abs(comp(jjj,2)).gt.xpmax)xpmax=abs(comp(jjj,2)) 89 continue do 91 jjj=1,kscan comp(jjj,3)=comp(jjj,3)/xpmax comp(jjj,2)=comp(jjj,2)/xpmax 91 continue c... x/y on screen is about 1.84 xpmax=1.84 ypmax=1.0 write(fint,'(10x,44a1,2i5)')(title(j:j),j=1,44),xlo,xhi call scr_pl(0,0,comp(1,2),dum,comp(1,3),dum,kscan,-ypmax,ypmax, + -xpmax,xpmax,fint) call scr_dline(-xpmax,0.,xpmax,0.) call scr_dline(0.,-ypmax,0.,ypmax) 90 continue call scr_pline(' press "C" to quit or hardcopy option') call scr_xhair(xlo,xval,yval,edc,iopt) if(iopt.eq.9)goto 55 if(edc.ne.'C')goto 90 line=' ' goto 1000 end **************************************** subroutine plt_mat(io,kfl,xlo,xhi,kscan,nwin,fw, & fmin,fmax,istep,ires,nout,iwave,iflip) **************************************** c... two dimensional plot of largest singular value d1 c... (time-frequency) c... after plotting d1, amplitudes of VERT or TRA can c... be added c... - the length of the moving window is given by c... interactive I/O but can be changed by option #3 c... - option #1 provides picking of travel times, followed c... by an moving window polarization analysis for the c... interpolated frequencies along the travel time c... curve c... parameter(nndim=5000,ipow=4096,icol=20) integer*4 spoint,epoint,xlo,xhi common/result/ans(ipow,icol) common/resule/epshh(ipow),epsvh(ipow) common/resuli/itmat(ipow) common/head$/nsta,nchn,iy,id,ih,im,ss,si,nscan,title common/env$$$/spoint,epoint,nser,npt,ymax,ymin,screen(1022,3) & ,il,ix common/scayle/zmin,zmax common/sdata/comp(nndim,3) common/slect$/autosc common/flag$$/flag,badun common/scrtch/scrbuf(nndim) common/myhed/jscan,jsta,jchn,jtyp,jy,jd,jh,jm,ssj,dtj,spare(20) dimension hbuf(30),fnm(30) dimension fmat(50),tmat(50) dimension ans1(ipow,icol) character*32 text(12) character*80 title,line character*1 edc equivalence (jscan,hbuf) data dflag,badf/5.e15,0.7/ ndim=nndim text(1)=' pick travel times ' text(2)=' quit ' text(3)=' change window length ' text(4)=' add VERT amplitude ' text(5)=' add TRA amplitude ' open(44,file='ampv',status='scratch') open(45,file='ampt',status='scratch') call scr_init(text,5) ampvm=-1. amptm=-1. c... calculate first 2D matrix of largest singular value c... and wait for mouse input goto 15 5 if(ires.eq.1)return line='x-hair options: A=command;C=quit' call scr_pline(line) call scr_xhair(lo,xval,yval,edc,iopt) if(iopt.ne.0)goto 10 if(edc.eq.'C')return 10 continue goto (25,35,45,65,75),iopt c... option 1: pick travel times and do MTPA after pressing c... mouse button "C" (goto 55) 25 continue npf=0 line='x-hair options: A=get value;B=cancel last;C=quit' call scr_pline(line) 24 call scr_pline(' ') call scr_pline(' B=cancel last; C=stop') call scr_xhair(lo,xval,yval,edc,iopt) if(iopt.ne.0)goto 10 if(edc.eq.'C')goto 55 if(npf.ge.50)goto 55 if(edc.eq.'B'.and.npf.gt.0)then write(line,102)xval,yval call scr_pline(line) npf=npf-1 goto 24 endif npf=npf+1 fmat(npf)=xval tmat(npf)=yval write(line,100)xval,yval call scr_pline(line) goto 24 c... calculate 2D matrix of d1 c... ymin, ymax are seismogram index limits for the matrix c... ind1, ind2 are the actual index limits for the seismogram 15 ind1=max0(1,xlo-kscan/2) ind2=min0(npt,xhi+kscan/2) ymin=xlo ymax=xhi fmax=amin1(fmax,1./(2.*si)) xmin=fmin*1000. xmax=fmax*1000. nstep=(xhi-xlo)/float(istep) c... loop over moving window do 602 i=1,nstep ind=ind1+(i-1)*istep yval=ymin+(i-1)*istep c... get data do 603 jser=1,3 mfl=kfl+jser-1 603 call gfs_rwdata(io,comp(1,jser),mfl,ind,kscan,'r',ierr) do 607 js=1,kscan 607 comp(js,3)=comp(js,3)*float(iflip) c... do MTPA call polar(ndim,kscan,si,nwin,fw,fmin,fmax,nout) if(i.eq.1)then do 604 jj=1,nout 604 scrbuf(jj)=ans(jj,1)*1000. call scr_pl2(0,ans(1,7),scrbuf,nout,ymin,ymax,yval, & xmin,xmax,title,0) else call scr_pl2(1,ans(1,7),scrbuf,nout,ymin,ymax,yval, & xmin,xmax,title,0) endif c... write out amplitudes on temporary file do 605 jj=1,nout if(ampvm.lt.ans(jj,2))ampvm=ans(jj,2) if(amptm.lt.ans(jj,4))amptm=ans(jj,4) 605 continue write(44,*)(ans(jj,2),jj=1,nout) write(45,*)(ans(jj,4),jj=1,nout) 602 continue call axis(xmin,xmax,ymin,ymax,1.,0) ampvm=max(ampvm,amptm) amptm=ampvm rewind(44) rewind(45) goto 5 c... option 3: choose another window length for the moving c... window 45 continue write(line,'(a,i5)')'old window length: ',kscan call scr_pline(line) call scr_rdnums('enter # samples for new window:$',fnm,num) if(num.ne.1)goto 45 kscan=fnm(1) ampvm=-1. amptm=-1. rewind(44) rewind(45) if(kscan.eq.0)goto 5 goto 15 c... do MTPA along interpolated travel time curve c... and set option for plotting results (ires=1) 55 continue do 200 ifr=1,nout itmat(ifr)=nint(ynterp(fmat,tmat,scrbuf(ifr),npf)) 200 continue do 201 i=1,nout ind=itmat(i)-kscan/2 do 202 jser=1,3 mfl=kfl+jser-1 202 call gfs_rwdata(io,comp(1,jser),mfl,ind,kscan,'r',ierr) do 608 js=1,kscan 608 comp(js,3)=comp(js,3)*float(iflip) call polar(ndim,kscan,si,nwin,fw,fmin,fmax,nout) do 203 j=1,icol 203 ans1(i,j)=ans(i,j) 201 continue do 204 i=1,nout do 204 j=1,icol 204 ans(i,j)=ans1(i,j) ires=1 goto 5 c... option 4: read VERT amplitudes from temporary file c... and add them to 2D plot with an offset in c... frequency 65 continue dscr=(scrbuf(2)-scrbuf(1))/4. do 303 i=1,nout 303 scrbuf(i)=scrbuf(i)+dscr do 301 i=1,nstep read(44,*)(ans(jj,2),jj=1,nout) do 302 jj=1,nout 302 ans(jj,2)=ans(jj,2)/ampvm ind=ind1+(i-1)*istep yval=ymin+(i-1)*istep call scr_pl2(1,ans(1,2),scrbuf,nout,ymin,ymax,yval, & xmin,xmax,title,1) 301 continue call axis(xmin,xmax,ymin,tmax,1.,0) do 304 i=1,nout 304 scrbuf(i)=scrbuf(i)-dscr goto 5 c... option 5: read TRA amplitudes from temporary file c... and add them to 2D plot with an offset in c... frequency 75 continue dscr=(scrbuf(2)-scrbuf(1))/2. do 403 i=1,nout 403 scrbuf(i)=scrbuf(i)+dscr do 401 i=1,nstep read(45,*)(ans(jj,4),jj=1,nout) do 402 jj=1,nout 402 ans(jj,4)=ans(jj,4)/amptm ind=ind1+(i-1)*istep yval=ymin+(i-1)*istep call scr_pl2(1,ans(1,4),scrbuf,nout,ymin,ymax,yval, & xmin,xmax,title,2) 401 continue call axis(xmin,xmax,ymin,tmax,1.,0) do 404 i=1,nout 404 scrbuf(i)=scrbuf(i)-dscr goto 5 c... option 2: close temporary files and quit 35 continue close(44,status='delete') close(45,status='delete') return 100 format('crshair coordinates: ',2f8.3) 102 format('cancelled coordinates: ',2f8.3) end ********************************************************* subroutine pltres(nout,si,fmin,fmax,iwave,xlo, & nwin,fw,kscan,ires,dazi,ddist,jo,title1) ********************************************************* parameter(nndim=5000,ipow=4096,icol=20) integer*4 xlo dimension smth(ipow) common/result/ans(ipow,icol) common/resule/epshh(ipow),epsvh(ipow) common/resuli/itmat(ipow) common/scrtch/freq(nndim) ********** c... MTPA results and headers in gfs format common/myhed/jscan,jsta,jchn,jtyp,jy,jd,dtj,nfr, & ffmin,ffmax,iwa,ixlo,iwin,ffw,iwinl, & azi,dist,iawa,nwt,xmin,xmax,scola,slong,jih,spare(6) common/srce/ky,kd,kh,km,ssec,st0,sp0,sd0 c... description of gfs header: c jscan: nout*21 (ires=0), nout*22 (ires=1) c jsta: station name c jchn: channel (VHT, VHE, TRA or E-W or BE) c jtyp: 4char station type c jy: year c jd: julian day c dtj: sampling interval c nfr: # of frequencies (nout) c ffmin: minimum frequency c ffmax: maximum frequency c iwa: wrapping (0,1,2,3) c ixlo: 0 for (ires=1), first index of window (ires=0) c ires=1 for moving window; ires=0 for sing. wind. c iwin: # of used multiple tapers c ffw: b-w product c iwinl: window length in samples c azi: azimuth of station c dist: epicentral distance c iawa: wave type (1 rayl, 2 love) c nwt: order # of wave train c xmin,xmax: last zoomed in window c scola: geograph. source colatitude c slong: longitude c jih: hour of event c... results are stored: each value for all frequencies c ans(1-8),ans(17),ans(9),ans(12),ans(13),ans(16), c epshh,ans(10),ans(11),epsvh,ans(14),ans(15),ans(19), c ans(20),itmat (for ires=1, only) ********** dimension hbuf(30),hbuf1(30),fnm(30) character*32 text(12) character*80 title,title1,namhcopy*9 character*120 line,line1,string*256 character*1 edc equivalence (jscan,hbuf) data dflag,badf/5.e15,0.7/ text(1)=' plot sing vals ' text(2)=' plot ad. spectra ' text(3)=' plot theta-2H/-3H ' text(4)=' plot theta-3V/-3T ' text(5)=' 3D-ellipticity ' text(6)=' eps-hh, eps-vh ' text(7)=' quit ' text(8)=' zoom in theta-2H ' text(9)=' hardcopy ' text(10)=' write out results (rayl) ' text(11)=' write out results (love) ' if(ires.eq.1)xlo=0 fac=1000. if(si.lt.0.999) fac=1. do 1 i=1,nout 1 freq(i)=ans(i,1)*fac jopt=1 xmin=fmin*fac xmax=fmax*fac xmn=xmin plen=xmax-xmin ymin=0. ymax=0. winl=1./(ans(2,1)-ans(1,1))/60. winb=xlo*si/60. scola=st0 slong=sp0 jih=kh line=' ' call scr_init(text,11) goto 15 5 continue line1='x-hair options: A=get val; B=zoom; C=unzoom or quit' call scr_pline(line1) CALL scr_xhair(lo,xval,yval,edc,iopt) if(iopt.ne.0) goto 10 if(edc.eq.'C') goto 201 if(edc.eq.'B') goto 202 if(edc.eq.'A')then call scr_pline(line) write(line1,506) xval,yval,winb,winl call scr_pline(line1) goto 5 endif c*** unzoom *** 201 if(xmax-xmin.eq.plen) goto 5 xmin=xmn xmax=xmn+plen goto (15,20,25,35,55,65,45,85,95,105,115),jopt c*** zoom *** 202 xsav=xval call scr_xhair(lo,xval,yval,edc,iopt) if(iopt.ne.0) goto 10 if(xval.le.xsav) goto 5 xmin=xsav xmax=xval goto (15,20,25,35,55,65,45,85,95,105,115),jopt c*** option selection 10 jopt=iopt goto (15,20,25,35,55,65,45,85,95,105,115),iopt c... plotting singular values 15 continue icolon=0 do 16 i=1,80 if(title1(i:i).eq.':')then icolon=i goto 17 endif 16 continue 17 continue write(title,'(a9,a12,a)')title1(1:9),title1(icolon-11:icolon), &' singular values$' ymin=0. ymax=1. call scr_pl(0,-2,ans(1,5),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call scr_pl(1,-3,ans(1,6),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call scr_pl(1,-7,ans(1,7),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call axis(xmin,xmax,ymin,ymax,1.,0) line=' ' goto 5 c... plotting spectra 20 continue title(23:40)='adaptive spectra$' call limit(ymin,ymax,ans(1,2),ans(1,3),ans(1,4),nout) call scr_pl(0,0,ans(1,2),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call scr_pl(1,20,ans(1,3),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call scr_pl(1,1,ans(1,4),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call axis(xmin,xmax,ymin,ymax,1.,0) line=' ' goto 5 c... plotting theta-2h and theta-3h 25 title(23:54)='theta-2H (cir), theta-3H (crs)$' c... calculate limits of angles ymin=1.e+30 ymax=-1.e30 call limit1(ymin,ymax,xmin,xmax,freq,ans(1,8),nout) call limit1(ymin,ymax,xmin,xmax,freq,ans(1,17),nout) c... plot theta-2h and evaluate mean value call scr_pl(0,-17,ans(1,8),ans(1,9),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,ans(1,8),ans(1,9),xmin,xmax,ava,era,nout) call scr_line(xmin,ava,xmax,ava) c... plot theta-3h and evaluate mean value c... also plot smoothed curve call scr_pl(1,-12,ans(1,17),ans(1,9),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call smooth(ans(1,17),smth,nout) call scr_pl(1,1,smth,dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,ans(1,17),ans(1,9),xmin,xmax,avp,erp,nout) write(line,502) ava,era,avp,erp call scr_pline(line) call scr_line(xmin,avp,xmax,avp) call scr_dline(xmin,0.,xmax,0.) call axis(xmin,xmax,ymin,ymax,1.,0) goto 5 c... plotting theta-3v and theta-3t 35 title(23:53)='theta-3V (cir)/theta-3T (crs)$' ymin=-90. ymax=90. call scr_pl(0,-17,ans(1,12),ans(1,13),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,ans(1,12),ans(1,13),xmin,xmax,ava,era,nout) call scr_line(xmin,ava,xmax,ava) call scr_pl(1,-2,ans(1,16),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avevar(freq,ans(1,16),xmin,xmax,ave,var,nout) call scr_line(xmin,ave,xmax,ave) write(line,503) ava,era,ave,var call scr_pline(line) call scr_dline(xmin,0.,xmax,0.) call axis(xmin,xmax,ymin,ymax,1.,0) goto 5 c... plotting 3D-ellipticity 55 title(23:38)='3D-ellipticity$' ymin=0. ymax=1. call scr_pl(0,-7,ans(1,19),dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avevar(freq,ans(1,19),xmin,xmax,ave,var,nout) call scr_line(xmin,ave,xmax,ave) call smooth(ans(1,19),smth,nout) call scr_pl(1,1,smth,dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) write(line,504) ave,var call scr_pline(line) call axis(xmin,xmax,ymin,ymax,1.,0) goto 5 c... plotting eps-hh and eps-vh 65 title(23:50)='eps-hh (cir), eps-vh (crs)$' ymin=0. ymax=1. call scr_pl(0,-17,epshh,ans(1,11),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,epshh,ans(1,11),xmin,xmax,av1,er1,nout) call scr_line(xmin,av1,xmax,av1) call scr_pl(1,-12,epsvh,ans(1,15),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,epsvh,ans(1,15),xmin,xmax,av2,er2,nout) call scr_line(xmin,av2,xmax,av2) write(line,505) av1,er1,av2,er2 call scr_pline(line) call axis(xmin,xmax,ymin,ymax,1.,0) goto 5 c... zoom in theta-2h 85 title(23:54)='theta-2H (cir), theta-3H (crs)$' c... calculate limits of angles ymin=1.e+30 ymax=-1.e30 call limit2(ymin,ymax,xmin,xmax,freq,ans(1,8),ans(1,9),nout) c... plot theta-2h and evaluate mean value call scr_pl(0,-17,ans(1,8),ans(1,9),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,ans(1,8),ans(1,9),xmin,xmax,ava,era,nout) call scr_line(xmin,ava,xmax,ava) c... plot theta-3h and evaluate mean value c... also plot smoothed curve call scr_pl(1,-12,ans(1,17),ans(1,9),freq,dum,nout,ymin,ymax, + xmin,xmax,title) call smooth(ans(1,17),smth,nout) call scr_pl(1,1,smth,dum,freq,dum,nout,ymin,ymax, + xmin,xmax,title) call avbnd(freq,ans(1,17),ans(1,9),xmin,xmax,avp,erp,nout) write(line,502) ava,era,avp,erp call scr_pline(line) call scr_line(xmin,avp,xmax,avp) call scr_dline(xmin,ava+era,xmax,ava+era) call scr_dline(xmin,ava-era,xmax,ava-era) call scr_dline(xmin,0.,xmax,0.) call axis(xmin,xmax,ymin,ymax,1.,0) goto 5 c... option 9: hardcopy 95 continue npcops=npcops+1 namhcopy=' ' if(npcops.le.9)then write(namhcopy,57)npcops elseif(ncops.ge.10.and.npcops.le.99)then write(namhcopy,58) npcops else write(namhcopy,59) npcops endif 57 format('hrdcpp',i1) 58 format('hrdcpp',i2) 59 format('hrdcpp',i3) call scr_rdtext('enter output name:',string,nstrin) call scr_dump(7,string) goto 5 c... options 10 and 11: write out results to gfs file 105 iawa=1 goto 116 115 iawa=2 116 continue call scr_rdnums('enter order # of wave train or 99 to quit:$' +,fnm,num) if(num.ne.1)goto 5 nwt=fnm(1) if(nwt.eq.99)goto 5 iopt=1 jopt=1 call scr_pline(' ') call scr_pline(' writing out results') do 106 i=1,30 106 hbuf1(i)=hbuf(i) jscan=nout*21 if(ires.eq.1)jscan=nout*22 if(jscan.gt.ipow)then line1=' no output, # of results too large' call scr_pline(line1) goto 108 endif nfr=nout ffmin=fmin ffmax=fmax iwa=iwave ixlo=xlo dtj=si iwin=nwin ffw=fw iwinl=kscan azi=dazi dist=ddist do 109 j=1,8 do 109 i=1,nout 109 smth(nout*(j-1)+i)=ans(i,j) do 110 i=1,nout smth(8*nout+i)=ans(i,17) smth(9*nout+i)=ans(i,9) smth(10*nout+i)=ans(i,12) smth(11*nout+i)=ans(i,13) smth(12*nout+i)=ans(i,16) smth(13*nout+i)=epshh(i) smth(14*nout+i)=ans(i,10) smth(15*nout+i)=ans(i,11) smth(16*nout+i)=epsvh(i) smth(17*nout+i)=ans(i,14) smth(18*nout+i)=ans(i,15) smth(19*nout+i)=ans(i,19) smth(20*nout+i)=ans(i,20) 110 continue if(ires.eq.1)then do 111 i=1,nout 111 smth(21*nout+i)=itmat(i) endif call gfs_rwentry(jo,hbuf,smth,0,'w') 108 continue do 107 i=1,30 107 hbuf(i)=hbuf1(i) goto 5 c... finishing of plotting POLAR results 45 return 502 format('average value of theta-2H is ',f8.3,'+-',f7.3,5x, + 'average value of theta-3H is ',f8.3,'+-',f7.3,'$') 503 format('average value of theta-3V is ',f8.3,'+-',f7.3,5x, + 'average value of theta-3T is ',f8.3,'+-',f8.3,'$') 504 format('average value of 3D-ellip is ',f8.3,'+-',f7.3) 505 format('average value of eps-hh is ',f8.3,'+-',f7.3,5x, + 'average value of eps-vh is ',f8.3,'+-',f7.3,'$') 506 format('crshair coordinates: ',2f8.3, +' window beg. length:',2f5.0) end ********************************* function ynterp(x,y,xi,n) ********************************* c x(i) must be all different c HIROO KANAMORI AUG.4 1978 c 13.07.92 extracted from ITT1SUB.FOR dimension x(1), y(1) i=1 j=n k=1 if ( x(1) .gt. x(n) ) goto 11 1 if( xi-x(k)) 2,2,3 2 j=k goto 4 3 i=k 4 k=(i+j)/2 if(k.ne.i) goto 1 if( xi.gt.x(n)) j=n+1 goto 21 11 if( x(k)-xi) 12,12,13 12 j=k goto 14 13 i=k 14 k=(i+j)/2 if(k.ne.i) goto 11 if( xi.lt. x(n)) j=n+1 21 if ( j.eq. 1) goto 23 if(j.eq.n+1) goto 24 ynterp=y(j-1)+(y(j)-y(j-1))*(xi-x(j-1))/(x(j)-x(j-1)) return 23 ynterp=y(1) return 24 ynterp=y(n) return end