program he_relp c c ...... This version of the Regional Earthquake Location Program (RELP) c is constructed for the historical earthquake relocation project. c c Last Modified - 23 Feb 2000 - RAU c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) include "config.h" character*72 id character*80 tscfile,astrng character*1 polar character*4 stn character*5 chpid real*8 t1p_time integer*4 ncr,nsp,ns1,ns2,ns3,ns4,nsreq,nstns,ipid,getpid integer*4 m,n,nmax,nit,nitmax,system,nr_switch common/files/ ncr,nsp,ns1,ns2,ns3,ns4 common/comp/a(1000,4),b(1000),w(1000),q(4),qi(4),u(1000,4), 1 v(4,4),cov(4,4),x(4),bmax(1000),se common/obs/pt(1000),st(1000),nsreq common/obs1/stn(1000),polar(1000) common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns common/epid/id common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et common/dim/m,n,nmax,nit common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd common/nm/nitmax common/qf/qf common/mmax/mmax integer*4 swhdf,swhdx,swnit,swnmx,swcon common/swtch/swhdf,swhdx,swnit,swnmx,swcon common/step/step common /mags/ amla,seamla,amlb,seamlb,nmla,nmlb common /eqelp/ del_ns,hor_se,iaz_gap,np_loc,ns_loc c character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam, 1 epi_dt,me_stn_adj_file,stn_adj_file, 2 vel_mod_file character*9 epi_id common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam, 1 epi_id,epi_dt character*80 geo_info common /ginfo/geo_info integer*4 ipar,iplot,imagn,istna,ivelm,imevt,iniloc character*1 mcode c real*8 meglat,melat,melon,medep,meot,mesec integer*4 mehr,memn,meflag,mwflag common /mehypo/ meglat,melat,melon,medep,meot,mehr,memn, 1 mesec,meflag real*4 il_lat,il_lon,fl_lat,fl_lon common /tester/ il_lat,il_lon,fl_lat,fl_lon data n_retry,n_zdse /0,0/ c c ...... debugging aid c c call abrupt_underflow() c i = ieee_handler("set","common",SIGFPE_ABORT ) c c ...... name I/O files c call fnamer(iplot,imagn,istna,ivelm,imevt,mwflag, 1 iniloc,nr_switch,me_stn_adj_file,stn_adj_file, 2 vel_mod_file) c c ...... debugging check on file names c c print*, 'Model directory: ',model_dir(1:lnblnk(model_dir)) c print*, 'Model types: ',models(1:lnblnk(models)) c print*, 'Instr file: ',instr_coord(1:lnblnk(instr_coord)) c print*, 'Location file: ',landmark_file(1:lnblnk(landmark_file)) c print*, 'BDSN stn adj: ',mladj_bdsn(1:lnblnk(mladj_bdsn)) c print*, 'BRK stn adj: ',mladj_brk(1:lnblnk(mladj_brk)) c c ...... input seistool pickfile data c call st_pfd(t1p_time) c c do 432 ijkl=1,nstns c if(pt(ijkl).gt.0.d0) write(*,*) 'pt(',ijkl,') = ',pt(ijkl) c if(st(ijkl).gt.0.d0) write(*,*) 'st(',ijkl,') = ',st(ijkl) c 432 continue c call parset(n,nit,nmax,mmax,nitmax,ncr,nsp) c ...... input initial velocity model call mdl(ivelm,vel_mod_file) neq=neq+1 c c ...... if master event location supplied, skip to final solution c if(meflag.eq.1) then ex=melon ey=meglat ez=medep edepth=ez et=meot c if(et.gt.t1p_time) et=et-86400.d0 go to 25 endif c c ...... calculate or input initial epicenter c if(iniloc.eq.0) then call st_ilc c write(*,*) 'ilc - et = ',et else ex=melon ey=meglat ez=medep edepth=ez et=meot go to 25 c if(et.gt.t1p_time) et=et-86400.d0 endif il_lat=ey il_lon=ex c c ...... extract velocity model and station adjustments c (the program uses geocentric coordinates internally) c elatd=datand(dtand(ey)/0.993277d0) elongd=ex edepth=ez ipar=1 c ...... extract velocity model call st_msa(ipar,istna,stn_adj_file,ivelm,vel_mod_file, 1 elatd,elongd,stn(1),mcode,vpa,vpb,vpc, 1 vsa,vsb,vsc,vh,psaf(1),ssaf(1)) c write(*,'(a,2f10.4,3x,a1)') 'lat, long & mcode (1) = ', c 1 elatd,elongd,mcode ipar=2 c ...... extract station adjustments do 1 i=1,nstns call st_msa(ipar,istna,stn_adj_file,ivelm,vel_mod_file, 1 elatd,elongd,stn(i),mcode,vpa,vpb,vpc, 1 vsa,vsb,vsc,vh,psaf(i),ssaf(i)) c write(*,'(a,a4,2f10.3)') 'stn, padj & sadj = ', c 1 stn(i),psaf(i),ssaf(i) 1 continue c c ...... determine location with velocity model and station c adjustments inferred from initial location c nit=0 2 nit=nit+1 c write(*,*) '01 - nit & et = ',nit,' ',et call nfparam(meflag) call abcomp call weight(meflag,mwflag) call solve(meflag) call pertur(n_zdse) c write(*,*) '02 - nit & et = ',nit,' ',et c write(*,*) 'st_relp: call conver - 2' call conver if(nit.le.nitmax) then if(nit.le.swnit) then go to 2 elseif(swcon.eq.0) then go to 2 endif endif c c ...... extract velocity model and station adjustments again c based on refined location c edepth=8.d0 ez=edepth 25 continue elatd=datand(dtand(ey)/0.993277d0) elongd=ex ipar=1 c ...... extract velocity model call st_msa(ipar,istna,stn_adj_file,ivelm,vel_mod_file, 1 elatd,elongd,stn(1),mcode,vpa,vpb,vpc, 1 vsa,vsb,vsc,vh,psaf(1),ssaf(1)) c write(*,'(a,2f10.4,3x,a1)') 'lat, long & mcode (2) = ', c 1 elatd,elongd,mcode ipar=2 c ...... extract station adjustments do 3 i=1,nstns call st_msa(ipar,istna,stn_adj_file,ivelm,vel_mod_file, 1 elatd,elongd,stn(i),mcode,vpa,vpb,vpc, 1 vsa,vsb,vsc,vh,psaf(i),ssaf(i)) c write(*,'(a,a4,2f10.3)') 'stn, padj & sadj = ', c 1 stn(i),psaf(i),ssaf(i) 3 continue c c write(*,*) 'V model = ',vpa,vpb,vpc,vsa,vsb,vsc,vd c do 65 i=1,nstns c write(*,*) 'Stn adj = ',i,stn(i),psaf(i),ssaf(i) c 65 continue c c c ...... compute final earthquake location c call parset(n,nit,nmax,mmax,nitmax,ncr,nsp) if(n_retry.eq.1) n=3 nit=0 4 nit=nit+1 c write(*,*) '03 - nit & et = ',nit,' ',et c write(*,*) 'st_relp: nit = ',nit c write(*,*) 'st_relp: ex, ey, ez & et = ',ex,ey,ez,et call nfparam(meflag) if(n_zdse.eq.1) n=3 if(meflag.eq.1) n=4 call abcomp c write(*,*) 'st_relp: n & m = ',n,m call weight(meflag,mwflag) call solve(meflag) c write(*,*) '04 - nit & et = ',nit,' ',et if(meflag.lt.1) then call pertur(n_zdse) c write(*,*) '05 - nit & et = ',nit,' ',et c write(*,*) 'st_relp: call conver - 1' call conver else go to 45 endif if(nit.le.nitmax) then if((swcon.ne.0).and.(n_retry.eq.1)) go to 45 if(nit.le.swnit) then go to 4 elseif(swcon.eq.0) then go to 4 endif endif 45 continue c write(*,*) '06 - nit & et = ',nit,' ',et c c if(meflag.lt.1) then c if((nr_switch.eq.1).and.(n_retry.eq.0)) then c c ...... check for suspect solution c - del_ns > 250 km c - nit > 40 c - se > 1 sec c - cov(4,4) > 225 c c if(((del_ns.gt.250.d0).or.(nit.gt.40)).or. c 1 ((se.gt.1.d0).or.(cov(4,4).gt.225.d0))) then c n_retry=1 c mwflag=1 c call rinloc(ex,ey,pt_min) c et=pt_min-2.d0 c edepth=8.d0 c ez=edepth c if(cov(4,4).gt.225.d0) then c n_zdse=1 c n=3 c endif c go to 25 c endif c endif c endif c c ...... calculate ML c elatd=ey elongd=ex fl_lat=elatd fl_lon=elongd call eif_ldif if(imagn.lt.1) then edepx=edepth c elat_1=datand(dtand(elatd)/0.993277d0) elat_1=elatd call st_dml(elat_1,elongd,edepx,amla,seamla,nmla,amlb,seamlb, 1 nmlb,fnamam,fnamml,epi_id) endif c c ...... determine geographical landmark info c elongd=ex elatd=datand(dtand(ey)/0.993277d0) call st_glc(elatd,elongd,geo_info) c c ...... output solution c c write(*,'(a,2f10.4)') 'call st_sln - elatd & elongd = ', c 1 elatd,elongd call st_sln(imagn,imevt,mcode,me_stn_adj_file) c c ...... delete output files and abort program with return code = 1 c if the number of stations used in the solution is < 3 c if(ns_loc.lt.3) then write(*,*) 'Number of stations < 3, program STOPPED' if(imagn.gt.0) then open(33,file=fnamml(1:lnblnk(fnamml)),status='unknown') close(33,status='delete') endif open(33,file=fnamsl(1:lnblnk(fnamsl)),status='unknown') close(33,status='delete') open(33,file=fnamel(1:lnblnk(fnamel)),status='unknown') close(33,status='delete') ipid=getpid() write(chpid,'(i5.5)') ipid write(tscfile,'(a15,a5)') '/tmp/stn_coord.',chpid open(39,file=tscfile(1:lnblnk(tscfile)),status='unknown') close(39,status='delete') stop 1 endif c c ...... plot solution c c if(iplot.lt.1) then c if(imagn.lt.1) then c call st_plot(imagn,elatd,elongd,edepth,amla,epi_dt,geo_info, c 1 fnamps) c else c aml=-1. c call st_plot(imagn,elatd,elongd,edepth,aml,epi_dt,geo_info, c 1 fnamps) c endif c endif c c ...... normal procedure termination c ipid=getpid() write(chpid,'(i5.5)') ipid write(tscfile,'(a15,a5)') '/tmp/stn_coord.',chpid open(39,file=tscfile(1:lnblnk(tscfile)),status='unknown') close(39,status='delete') c write(astrng,'(a4,a,a1)') 'cat ',fnamsl(1:lnblnk(fnamsl)),' ' i=system(astrng) c stop end subroutine eif_ldif c c ...... determine difference between initial & final epicenter c implicit real*4 (a-h,o-z) implicit integer*4 (i-n) real*4 il_lat,il_lon,fl_lat,fl_lon common /tester/ il_lat,il_lon,fl_lat,fl_lon character*80 aline integer*4 getpid,pid character*5 chpid c pid=getpid() write(chpid,'(i5.5)') pid aline='./st_relp_epi.'//chpid dlat=111.195*(fl_lat-il_lat) dlon=111.195*cosd((fl_alt+il_lat)/2.)*(fl_lon-il_lon) c open(44,file=aline(1:lnblnk(aline)),status='unknown') c write(44,'(a,2f10.3)') 'lat & lon diff (km) = ',dlat,dlon c close(44) c aline='chmod 777 ./st_relp_epi.'//chpid c i=system(aline(1:lnblnk(aline))) c return end subroutine rinloc(ex,ey,pt_min) c c ...... try center of observing network for initial location c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns common/obs/pt(1000),st(1000),nsreq elatd=0.d0 elond=0.d0 do 1 i=1,nstns elatd=elatd+syf(i) elond=elond+sxf(i) 1 continue ey=elatd/dble(nstns) ex=elond/dble(nstns) pt_min=100000.d0 do 2 i=1,nsreq if(pt(i).gt.0.d0) then if(pt_min.gt.pt(i)) pt_min=pt(i) endif 2 continue c return end subroutine st_info c c ...... print out information about st_relp c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c write(*,'(1h ,a)') 1 'st_relp - Regional Earthquake Location Program' c write(*,'(1h ,a)') 'Syntax:' write(*,'(1h ,7x,a)') 'st_relp [-h] [-p] [-b] [-f sc_file] ' write(*,'(1h ,14x,a)') '[-I lat lon dep hr mn sec] ' write(*,'(1h ,14x,a)') '[-L lat lon dep hr mn sec] ' write(*,'(1h ,14x,a)') 1 '[-g gl_file] [-q ep_vt ep_dt ep_tt ep_ct] ' write(*,'(1h ,14x,a)') '[-m] [-M mesa_file] [-s sa_file] ' write(*,'(1h ,14x,a)') '[-v vm_file] [-w] pickfile' c write(*,'(1h ,a)') 'where:' write(*,'(1h ,7x,a)') 1 '-h Help - prints this help message' write(*,'(1h ,7x,a)') 1 '-p Plot - used to inhibit plot generation' write(*,'(1h ,7x,a)') 1 ' (for use on non-plotting terminals)' write(*,'(1h ,7x,a)') 1 '-b BDSN - adds "b" to beginning of file names' write(*,'(1h ,7x,a)') 1 ' (for distinguishing solutions obtained' write(*,'(1h ,7x,a)') 1 ' using only Berkeley data)' write(*,'(1h ,7x,a)') 1 ' (for signifing broadband data solution)' c write(*,'(1h ,7x,a)') 1 '-f Ascf - input alternate station coordinate file' write(*,'(1h ,7x,a)') 1 ' (REDI needs a different station ' write(*,'(1h ,7x,a)') 1 ' coordinate file to avoid ambiguity)' write(*,'(1h ,7x,a)') 1 ' (assumes same format as the instr.coord' write(*,'(1h ,7x,a)') 1 ' ie: stn code, lat & long (a4,1x,2f9.4))' c write(*,'(1h ,7x,a)') 1 '-g Glfl - input alternate geographical location file' write(*,'(1h ,7x,a)') 1 ' (lat, long & place name (2f10.4,3x,a)' write(*,'(1h ,7x,a)') 1 ' (lat & lon are positive for N & E)' c write(*,'(1h ,7x,a)') 1 '-I Iloc - initial event location ' write(*,'(1h ,7x,a)') 1 ' input "lat lon dep hr mn sec" ' write(*,'(1h ,7x,a)') 1 ' where lat, lon, dep & sec are floats ' write(*,'(1h ,7x,a)') 1 ' and hr & mn are integers ' write(*,'(1h ,7x,a)') 1 ' (lat & lon are positive for N & E)' c write(*,'(1h ,7x,a)') 1 '-L Lmev - master event location ' write(*,'(1h ,7x,a)') 1 ' input "lat lon dep hr mn sec" ' write(*,'(1h ,7x,a)') 1 ' where lat, lon, dep & sec are floats ' write(*,'(1h ,7x,a)') 1 ' and hr & mn are integers ' write(*,'(1h ,7x,a)') 1 ' (lat & lon are positive for N & E)' c write(*,'(1h ,7x,a)') 1 '-m Magn - used to inhibit magnitude calculation' write(*,'(1h ,7x,a)') 1 ' (for compatability with REDI)' c write(*,'(1h ,7x,a)') 1 '-M Mevt - master event location proceedure' write(*,'(1h ,7x,a)') 1 ' write station adjustments to mesa_file' c write(*,'(1h ,7x,a)') 1 '-q Qpar - event type threshold parameters ' write(*,'(1h ,7x,a)') 1 ' input "ep_vt ep_dt ep_tt ep_ct" ' write(*,'(1h ,7x,a)') 1 ' where ep_vt, ep_dt, ep_tt & ep_ct are floats' write(*,'(1h ,7x,a)') 1 ' ep_vt = phase velocity threshold (km/sec), ' write(*,'(1h ,7x,a)') 1 ' ep_dt = cluster distance threshold ' write(*,'(1h ,7x,a)') 1 ' (fraction of inter-quartile distance) & ' write(*,'(1h ,7x,a)') 1 ' ep_tt = cluster time threshold (sec).' write(*,'(1h ,7x,a)') 1 ' ep_ct = curvature threshold (rad/km), ' c write(*,'(1h ,7x,a)') 1 '-r Rilo - retry initial location if suspect solution' write(*,'(1h ,7x,a)') 1 ' is well outside network, ie:' write(*,'(1h ,21x,a)') '- dist to nearest stn > 250 km' write(*,'(1h ,21x,a)') '- no. of iterations > 40' write(*,'(1h ,21x,a)') '- standard error > 1 sec' write(*,'(1h ,21x,a)') '- depth standard error > 15 km' c write(*,'(1h ,7x,a)') 1 '-R REDI - set switches for REDI use, ie:' write(*,'(1h ,21x,a)') 'set "-p" to inhibit plottting' write(*,'(1h ,21x,a)') 'set "-m" to inhibit ML calculation' write(*,'(1h ,21x,a)') 'set "-r" to retry initial location' write(*,'(1h ,21x,a)') ' for suspect solutions' c write(*,'(1h ,7x,a)') 1 '-s Sadj - input sa_file station adjustment file' write(*,'(1h ,7x,a)') 1 ' which is used in place of default station' write(*,'(1h ,7x,a)') 1 ' adjustment file (if sa_file .eq. "0", all' write(*,'(1h ,7x,a)') 1 ' stations adjustments are set to 0.d0)' c write(*,'(1h ,7x,a)') 1 '-v Vmod - input vm_file velocity model file' write(*,'(1h ,7x,a)') 1 ' which is used in place of default velocity' write(*,'(1h ,7x,a)') 1 ' model file' c write(*,'(1h ,7x,a)') 1 '-w Wght - force all weights to unity' c write(*,'(1h ,7x,a)') 'pickfile' write(*,'(1h ,14x,a)') 'seistool generated pickfile' c write(*,'(1h ,a)') 'Notes:' write(*,'(1h ,7x,a)') 'The pickfile must end in ".pf" and be' write(*,'(1h ,14x,a)') 'located in the current working directory' write(*,'(1h ,7x,a)') 'The associated "s*.af" file must also be' write(*,'(1h ,14x,a)') 'in the current working directory if ML' write(*,'(1h ,14x,a)') 'is to be calculated' write(*,'(1h ,7x,a)') 'Do not specify both "-I" & "-L" at the' write(*,'(1h ,14x,a)') 'same time' write(*,'(1h ,7x,a)') 'No *.el, *.sl, *.ml or *.ps file is' write(*,'(1h ,14x,a)') 'generated when the are insufficient' write(*,'(1h ,14x,a)') 'observing stations to obtain a solution' c write(*,'(1h ,7x,a)') 'The format for the station adjustment file' write(*,'(1h ,14x,a)') '(for the mesa_file and sa_file files)' write(*,'(1h ,14x,a)') 'is stn code, P adj & S adj (1x,a4,2f6.2)' write(*,'(1h ,14x,a)') 'where:' write(*,'(1h ,21x,a)') 'stn code = 4 character station code' write(*,'(1h ,21x,a)') 'P adj = P station adjustment (sec)' write(*,'(1h ,21x,a)') 'S adj = S station adjustment (sec)' c write(*,'(1h ,7x,a)') 'The format for the velocity model file is' write(*,'(1h ,14x,a)') 1 'vpa, vpb, vpc, vsa, vsb, vsc & vd (7f10.2)' write(*,'(1h ,14x,a)') 'where:' write(*,'(1h ,21x,a)') 'vpa = surface P velocity (km/sec)' write(*,'(1h ,21x,a)') 1 'vpb = layer P velocity gradient (km/sec)/km' write(*,'(1h ,21x,a)') 'vpc = halfspace P velocity (km/sec)' write(*,'(1h ,21x,a)') 'vsa = surface S velocity (km/sec)' write(*,'(1h ,21x,a)') 1 'vsb = layer S velocity gradient (km/sec)/km' write(*,'(1h ,21x,a)') 'vsc = halfspace S velocity (km/sec)' write(*,'(1h ,21x,a)') 'vd = thickness of layer (km)' write(*,'(1h ,7x,a)') 1 'If "-v" or "-I" are used "-s" must also be specified' c write(*,'(1h ,7x,a)') 'Once a mesa_file is generated for a master' write(*,'(1h ,14x,a)') 1 'event, use that same file as the station adjustment' write(*,'(1h ,14x,a)') 1 'file (sa_file) to locate the other events relative' write(*,'(1h ,14x,a)') 1 'to the master event' write(*,'(1h ,7x,a)') 'The event type threshold parameters: ' write(*,'(1h ,14x,a)') 'ep_vt, ep_dt & ep_tt are used to flag ' write(*,'(1h ,14x,a)') 'an event as a telemetry glitch, as ' write(*,'(1h ,14x,a)') 'having unclustered stations or as a ' write(*,'(1h ,14x,a)') 'teleseism. The default values are' write(*,'(1h ,14x,a)') '10.0, 2.5, 0.015 & 0.005, respectively,' write(*,'(1h ,14x,a)') 'for ep_vt, ep_dt, ep_tt & ep_ct.' c write(*,'(1h ,a)') 'Example:' write(*,'(1h ,7x,a)') 'st_relp ARC.HHE.D.93.016.0628.pf' c write(*,'(1h ,a)') 'Output:' write(*,'(1h ,7x,a)') 'r930160628.el' write(*,'(1h ,14x,a)') 'earthquake location file' write(*,'(1h ,7x,a)') 'r930160628.sl' write(*,'(1h ,14x,a)') 'solution listing file' write(*,'(1h ,7x,a)') 'r930160628.ml' write(*,'(1h ,14x,a)') 'magnitude information file' write(*,'(1h ,7x,a)') 'r930160628.ps' write(*,'(1h ,14x,a)') 'plot postscript file' c return end subroutine fnamer(iplot,imagn,istna,ivelm,imevt,mwflag, 1 iniloc,nr_switch,me_stn_adj_file,stn_adj_file, 2 vel_mod_file) c c ...... name I/O files c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c cccc lsg mod - add config files include "config.h" cccc lsg mod - add config files c integer*4 iplot,ibdsn,imagn,istna,ivelm,imevt,mwflag,nr_switch, 1 iniloc character*80 ficarg(30) character*12 template,ftype character*80 fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam, 1 epi_dt,me_stn_adj_file,stn_adj_file, 2 vel_mod_file character*9 epi_id common /fnamr/ fnampf,fnamsl,fnamel,fnamml,fnamps,fnamam, 1 epi_id,epi_dt c real*8 meglat,melat,melon,medep,meot,mesec integer*4 mehr,memn,meflag common /mehypo/ meglat,melat,melon,medep,meot,mehr,memn, 1 mesec,meflag real*8 ep_dt,ep_tt,ep_st,ep_vt,ep_ct common /eqtyp5/ ep_dt,ep_tt,ep_st,ep_ct c data template /'Pickfile V.2'/ c iplot=0 ibdsn=0 imagn=0 istna=0 ivelm=0 imevt=0 meflag=0 mwflag=0 nr_switch=0 iniloc=0 c c ...... initialize default event parameter thresholds c ep_dt=4.0d0 ep_tt=0.015d0 ep_vt=10.d0 ep_st=1.d0/ep_vt ep_ct=0.005d0 c nfic=iargc() if(nfic.gt.29) then write(*,'(a)') 'Illegal number of arguments' go to 1234 endif i=0 1 i=i+1 if((iniloc.eq.1).and.(meflag.eq.1)) then write(*,*) 'Illegal to specify both "-I" & "-L", STOPPED' call exit(1) endif c if(i.eq.nfic) go to 2 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) if(ficarg(i)(1:lfa).eq.'-h') then c c ...... print st_relp info c call st_info stop elseif(ficarg(i)(1:lfa).eq.'-p') then c c ...... inhibit plot (for non=plotting VDT's) c iplot=1 elseif(ficarg(i)(1:lfa).eq.'-b') then c c ...... adds "b" to beginning of output file names c to indicate broadband data solution c ibdsn=1 elseif(ficarg(i)(1:lfa).eq.'-f') then c c ...... specify alternate station coordinate file c i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(a)') instr_coord c write(*,*) '-f option: ',instr_coord(1:lfa) elseif(ficarg(i)(1:lfa).eq.'-g') then c c ...... specify alternate geographical location file c i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(a)') landmark_file c write(*,*) '-g option: ',landmark_file(1:lfa) elseif(ficarg(i)(1:lfa).eq.'-I') then c c ...... input initial location c iniloc=1 i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') melat meglat=datand(0.993277d0*dtand(melat)) i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') melon i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') medep i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(i)') mehr i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(i)') memn i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') mesec meot=3600.d0*dble(mehr)+60.d0*dble(memn)+mesec c c write(*,'(a,3f10.4,2i5,f10.4)') c 1 'iniloc = ',melat,melon,medep,mehr,memn,mesec c elseif(ficarg(i)(1:lfa).eq.'-L') then c c ...... input master event location c meflag=1 i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') melat meglat=datand(0.993277d0*dtand(melat)) i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') melon i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') medep i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(i)') mehr i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(i)') memn i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') mesec meot=3600.d0*dble(mehr)+60.d0*dble(memn)+mesec c elseif(ficarg(i)(1:lfa).eq.'-m') then c c ...... inhibit calculation of magnitudes c imagn=1 elseif(ficarg(i)(1:lfa).eq.'-M') then c c ...... calculate master event location and output corresponding c stations adjustments c imevt=1 if(i.eq.nfic) then go to 1234 endif i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) me_stn_adj_file = ficarg(i)(1:lfa) c elseif(ficarg(i)(1:lfa).eq.'-q') then c c ...... input threshold parameters to identify event type c i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') ep_vt ep_st=1.d0/ep_vt i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') ep_dt i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') ep_tt i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) read(ficarg(i)(1:lfa),'(f)') ep_ct c elseif(ficarg(i)(1:lfa).eq.'-r') then c c ...... retry initial location c nr_switch=1 elseif(ficarg(i)(1:lfa).eq.'-R') then c c ...... set switches for REDI use c nr_switch=1 imagn=1 iplot=1 c istna=2 c write(*,*) '-R option selected' elseif(ficarg(i)(1:lfa).eq.'-s') then c c ...... input station adjustment file (overrides default station c adjustment file) c istna=1 if(i.eq.nfic) then go to 1234 endif i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) stn_adj_file=ficarg(i)(1:lfa) c c ...... if stn_adj_file .eq. '0' set all station adjustments to zero c if((lfa.eq.1).and.(stn_adj_file(1:1).eq.'0')) then istna=2 endif elseif(ficarg(i)(1:lfa).eq.'-v') then c c ...... input gradient layer over halfspace velocity model c (overrides default velocity models) c ivelm=1 if(i.eq.nfic) then go to 1234 endif i=i+1 call getarg(i,ficarg(i)) lfa=lnblnk(ficarg(i)) vel_mod_file=ficarg(i)(1:lfa) elseif(ficarg(i)(1:lfa).eq.'-w') then c c ...... force all weights to unity c mwflag=1 elseif(i.eq.nfic) then c c ...... Determine if input file is in fact a seistool c format pickfile. Validity checked by looking c for "Pickfile V.2" on the first line of the file. c lfa=lnblnk(ficarg(i)) if(lfa.lt.4) then go to 1234 else if(ficarg(i)(lfa-2:lfa).ne.'.pf') then write(*,*) 'NOTE - input file must end with ".pf"' go to 1234 else open(19,file=ficarg(i)(1:lfa),status='old') read(19,'(a12)') ftype do 2 j=1,12 if(ichar(ftype(j:j)).ne.ichar(template(j:j))) then write(*,'(3a)') 'Input file: ', 1 ficarg(i)(1:lnblnk(ficarg(i))), 2 ' is not a valid pickfile' close(19) go to 1234 endif 2 continue close(19) endif go to 3 endif else go to 1234 endif go to 1 3 continue c c ...... extract filename c call fnextr(ficarg(nfic),epi_id) c c ...... convert date/time format c call dtconv(epi_id,epi_dt) c lfa=lnblnk(ficarg(nfic)) fnampf=ficarg(nfic)(1:lfa) if(ibdsn.lt.1) then fnamsl=fnampf(1:lfa-3)//'.sl' fnamel=fnampf(1:lfa-3)//'.el' fnamml=fnampf(1:lfa-3)//'.ml' fnamps=fnampf(1:lfa-3)//'.ps' fnamam=fnampf(1:lfa-3)//'.af' else fnamsl=fnampf(1:lfa-3)//'.sl' fnamel=fnampf(1:lfa-3)//'.el' fnamml=fnampf(1:lfa-3)//'.ml' fnamps=fnampf(1:lfa-3)//'.ps' fnamam=fnampf(1:lfa-3)//'.af' endif c c write(*,*) 'fnamer - ficarg = ',ficarg c write(*,*) 'fnamer - fnampf = ',fnampf c write(*,*) 'fnamer - fnamsl = ',fnamsl c write(*,*) 'fnamer - fnamel = ',fnamel c write(*,*) 'fnamer - fnamml = ',fnamml c write(*,*) 'fnamer - fnamps = ',fnamps c write(*,*) 'fnamer - fnamam = ',fnamam c return c c ...... invalid command line info if this statement is reached c 1234 write(*,*) 'Invalid command line arguments:' nfic=iargc() do 1235 ijkl=1,nfic call getarg(ijkl,ficarg(ijkl)) write(*,'(a,i3,2a)') ' ficarg(',ijkl,') = ', 1 ficarg(ijkl)(1:lnblnk(ficarg(ijkl))) 1235 continue stop 1 end subroutine dtconv(epi_id,epi_dt) c c ...... system call to caldate to convert format c implicit none integer year, yr, doy, month, day, hr, min character*80 epi_dt character*9 epi_id c read(epi_id(1:2),'(i)') yr read(epi_id(3:5),'(i)') doy read(epi_id(6:7),'(i)') hr read(epi_id(8:9),'(i)') min year=yr+1900 if(yr.lt.100) year=year+100 call dy_to_mdy(doy,year,month,day) write(epi_dt,20) yr,month,day,hr,min 20 format(i2.2,'/',i2.2,'/',i2.2,',',i2.2,':',i2.2) c return end subroutine fnextr(fninpt,epi_id) c c ...... extract date and time from first phase arrival c in pickfile c c NOTE: it is assumed that the first phase in the file is c also the earliest arriving phase! c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*80 fninpt,aline character*9 epi_id character*1 dummy character chyr*2,chdoy*3,chhr*2,chmn*2 c c write(*,*) 'fnextr - fninpt = ',fninpt(1:lnblnk(fninpt)) do 1 i=1,9 epi_id(i:i)=' ' 1 continue open(19,file=fninpt(1:lnblnk(fninpt)),status='old') read(19,'(a)') dummy read(19,'(a)') dummy read(19,'(a)') aline close(19) c c ...... decode pickfile phase time c call pfpt_dec(aline,chyr,chdoy,chhr,chmn) c epi_id=chyr//chdoy//chhr//chmn c epi_id=aline(iclb(2)-2:iclb(2)-1)// c 1 aline(iclb(2)+1:iclb(1)-1)// c 2 aline(iclb(1)+1:iclc(2)-1)// c 3 aline(iclc(2)+1:iclc(1)-1) epi_id=epi_id(1:lnblnk(epi_id)) c c j=10 c lfa=lnblnk(fninpt) c i=lfa+1 c 2 i=i-1 c if(i.lt.1) then c write(*,*) 'ILLEGAL file name - ',fninpt(1:lfa) c stop c endif c if(j.eq.1) go to 3 c achf=fninpt(i:i) c ichn=ichar(achf) c if((ichn.ge.48).and.(ichn.le.57)) then c j=j-1 c epi_id(j:j)=achf c endif c go to 2 c 3 continue c write(*,*) 'fnextr - epi_id = ',epi_id(1:lnblnk(epi_id)) c return end subroutine pfpt_dec(aline,chyr,chdoy,chhr,chmn) c c ...... decode phase onset time c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*80 aline integer*4 ich(4,2) character chyr*2,chdoy*3,chhr*2,chmn*2 c lnb=lnblnk(aline) i=lnb 1 i=i-1 if(aline(i:i).ne.':') go to 1 ich(4,2)=i-1 2 i=i-1 if(aline(i:i).ne.':') go to 2 ich(4,1)=i+1 ich(3,2)=i-1 3 i=i-1 if(aline(i:i).ne.' ') go to 3 ich(3,1)=i+1 ich(2,2)=i-1 4 i=i-1 if(aline(i:i).ne.' ') go to 4 ich(2,1)=i+1 5 continue if(aline(i:i).eq.' ') then i=i-1 go to 5 endif ich(1,2)=i ich(1,1)=i-1 c chyr=aline(ich(1,1):ich(1,2)) jd=ich(2,2)-ich(2,1)+1 if(jd.eq.3) then chdoy=aline(ich(2,1):ich(2,2)) elseif(jd.eq.2) then chdoy='0'//aline(ich(2,1):ich(2,2)) else chdoy='00'//aline(ich(2,1):ich(2,2)) endif jd=ich(3,2)-ich(3,1)+1 if(jd.eq.2) then chhr=aline(ich(3,1):ich(3,2)) else chhr='0'//aline(ich(3,1):ich(3,2)) endif jd=ich(4,2)-ich(4,1)+1 if(jd.eq.2) then chmn=aline(ich(4,1):ich(4,2)) else chmn='0'//aline(ich(4,1):ich(4,2)) endif c return end subroutine nfparam(meflag) c c ...... determine number of free parameters c c (3 for x, y & t and 4 for x, y, z & t) c implicit real*8 (a-h,o-z) integer*4 swhdf,swhdx,swnit,swnmx,swcon,meflag common/swtch/swhdf,swhdx,swnit,swnmx,swcon common/dim/m,n,nmax,nit c if(nit.eq.1) then swhdf=0 swhdx=0 swnit=100 swnmx=4 swcon=0 if(meflag.eq.1) then n=4 else n=3 endif iq=0 if(n.gt.m) n=m return endif if(swhdx.eq.1) then n=swnmx if(n.gt.m) n=m if(iq.eq.0) then iq=1 swnit=nit+2 endif endif c return end subroutine parset(n,nit,nmax,mmax,nitmax,ncr,nsp) c c ...... initialize parameters c implicit real*8 (a-h,o-z) n=3 nit=0 nmax = 4 mmax=1000 nitmax=50 ncr=1 nsp=2 c return end subroutine mdl(ivelm,vel_mod_file) c c ...... input velocity model c implicit real*8 (a-h,o-z) integer*4 ivelm character*80 vel_mod_file common/files/ncr,nsp,ns1,ns2,ns3,ns4 common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd c c ...... gradient layer over a half-space velocity model c c P-wave velocity = vpa + vpb * z in layer (km/sec) c = vpc in half-space (km/sec) c c S-wave velocity = vsa + vsb * z in layer (km/sec) c = vsc in half-space (km/sec) c c layer thickness = vd (in km) c if(ivelm.eq.0) then vpa=5.28 vpb=0.075 vpc=7.70 vsa=2.98 vsb=0.043 vsc=4.36 vd=25.0 else open(27,file=vel_mod_file(1:lnblnk(vel_mod_file)), 1 status='old') read(27,'(7f10.0)') vpa,vpb,vpc,vsa,vsb,vsc,vd close(27) endif c return end real*8 function dsqrtf(n,arg) c c ...... local dsqrtf function to track down domain errorr c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c if(arg.le.0.d0) then write(*,*) 'WARNING - sqrt DOMAIN ERROR -> n, arg = ', 1 n,' ',arg stop 1 else dsqrtf=dsqrt(arg) endif c return end