program mag2 implicit none ! force delcaration of all variables c c This program is to establish an ML scale using c c Nev Nstat c log Ajl = -n log(Rjl/R0) - k (Rjl - R0) + sum (Mj deljl) + sum (Sk delkl) - M0 c j=1 k=1 c c where: c c Ajl - observed amplitude at station l for event j c Rjl - hypocentral distance from event j to station l c R0 - reference distance c Mj - local magnitude of event j c Nev - number of events c Sk - correction for station l c Nstat - number of stations c deljl - Kronecker delta c delkl - " " c M0 - reference magnitude at R0 c n - geometrical spreading parameter c k - anelastic attenuation parameter c c Constraint = station corrections sum to zero c c changes c 07/01/2010 lo obtained copy form slsa, some changes on variable types c 16/04/2010 lo ignore component for site term determination c 01/07/2010 lo get distance range ratio from parameter file c does not compile with gfortran, check later, lo c 05/07/2010 lo determine seperate geometrical spreading part for sg/lg c 05/07/2010 lo changed reading of mag2.par to start in col 41 c 24/08/2010 lo changed order of a in ml scale, added option to fix scale, and to add noise c 22/09/2010 lo changed to 3 part geometrical spreading and removed bootstrap c 25/02/2011 Jh changed to use all phases used for Ml, warning if distace > 2500 km c 10/04/2013 lo added fixing of a in given range (RANGE A in mag2.par) c 2013 jh+lo add name of parameter file c include 'seidim.inc' ! dimensions for rea block include 'rea.inc' ! parameter common bliock character*80 data(5000) ! s-file with data in text array character*80 infile ! input file logical all ! true: read all data, false: headers integer code ! error return code integer nevent ! number of events in file integer i ! counter integer k ! counter integer no ! counter integer flag(10000) ! flag integer maxevent,maxobs,maxsite parameter (maxevent=1500) parameter (maxobs=10*maxevent) parameter (maxsite=300) integer year(maxevent) ! array containing data for all events integer month(maxevent) integer day(maxevent) integer hour(maxevent) integer min(maxevent) real sec(maxevent) real lat(maxevent) real lon(maxevent) real mag(maxevent) real depth(maxevent) integer nobs ! total number of amplitude observations real amp(maxobs) ! arrays containing individual data points real vg(maxobs) real tempamp(maxsite),tempdis(maxsite),templdis(maxsite) real a,b,corr,rms real dist(maxobs) character*5 stat(maxobs) character*2 co(maxobs) real per(maxobs) integer event_id(maxobs) character*8 phase(maxobs) real orient ! 0. = use horiz and vert ! 1. = use horiz only ! 2. = use vert only integer j ! counter character*1 disgard(3) ! array containing comps to ignore if orient=1. c real log_dist(maxobs) ! distance logs real inv_type ! 1. = SVD ! 2. = LSQR real min_dist ! set in mag2.par real max_dist ! set in mag2.par real min_distrange ! minimum distance range real mind,maxd real x character*2 comp_array(100) ! components to ditch character*5 stat_array(100) ! stations to ditch integer n_stat_array,n_comp_array character*5 stat_list(maxsite) ! max no of stations = 300 integer stat_nobs(maxsite) ! number of observations for each site c character*2 co_list(maxsite) ! integer nstat ! no of unique stations ! not to be confused with nstat in sub character*16 evid_list(maxevent) integer nev ! number of events for inversion integer flag_stat ! flag used for compiling station list integer flag_ev ! flag used for compiling event list integer ev_index(maxobs) integer stat_index(maxobs) integer err ! error code c real magscale(5) ! a, b and c in logA0 c real magscale_err(5) ! magscale error c change lo real magscale(7) ! a, b and c in logA0 real magscale_err(7) ! magscale error real siteterm(maxsite) ! site terms real siteterm_err(maxsite) ! site term errors real sourceterm(maxevent) ! source terms real sourceterm_err(maxevent) ! source term errors real ref(3) ! 1=dist, 2=amp, 3=mag integer boot_ninv ! number of inversions for bootstrapping real ran2 ! random number function ran2 integer idum ! -ve integer for initialising ran2 real seq(maxobs) ! array containing random numbers integer int_seq ! array containing rounded randoms converted to int real amp_sim(maxobs) ! arrays containing individual data points real dist_sim(maxobs) ! in simulated datasets character*5 stat_sim(maxobs) character*2 co_sim(maxobs) real per_sim(maxobs) character*16 evid_sim(maxobs) character*16 evid_sim_list(maxevent) character*5 stat_sim_list(maxsite) character*2 co_sim_list(maxsite) character*8 phase_sim(maxobs) integer maxinv parameter (maxinv=1) character*16 evid ! event id in yyyymmddhhmmssss format character*16 evid_arr(maxobs) real siteterm_arr(maxsite,maxinv) ! arrays containing the result of the inversions real sourceterm_arr(maxevent,maxinv) ! max no. of inversions = 1 + 4999 bootstrapped real magscale_arr(4,maxinv) integer nstat_sim integer nev_sim integer stat_sim_index(maxobs) integer ev_sim_index(maxobs) integer count real sum real average real xi, xbar real stdev logical synth real noise real magsynth(9) real hdist ! hypocentral distance integer seiclen integer nstat_ev integer minobs !minimum number of obs per event to use real mindratio real gdist(2) real fixscale(4),fixsite,rangea(3) real slat,slon,stel double precision msec real var,sigma real wagain character*240 mag_par_file wagain=2080. c Open output file that contains all data points open(unit=2,file='mag2_amp_obs.out',status='unknown') open(unit=3,file='mag2_events_read.out',status='unknown') open(unit=4,file='mag2_amp_dis.out',status='unknown') open(unit=8,file='mag2_stat_list.out',status='unknown') open(unit=20,file='mag2_events_used.out',status='unknown') c slsa - write out files containing stations, events and paths open(unit=28,file='mag2_statxy.out',status='unknown') open(unit=22,file='mag2_evxy.out',status='unknown') open(unit=23,file='mag2_paths.out',status='unknown') c c get input file name, check if exist c 9 continue c slsa10 c write(6,*) 'Susanne version' write(6,*) 'Give input file' read(5,'(a)') infile if (seiclen(infile).le.0) stop open(1,file=infile,status='old',err=10) goto 11 10 continue write(6,*)' No such input file' goto 9 11 continue c c get parameter file c write(6,*) 'Give parmeter file name, mag2.par is default = enter' read(5,'(a)') mag_par_file if(mag_par_file.eq.' ') mag_par_file='mag2.par' c all=.true. ! read all parameters nevent=0 ! initialize counter c Read parameter file call get_mag2_par(inv_type,orient,min_dist,max_dist, & comp_array,n_comp_array,stat_array,n_stat_array, & ref,boot_ninv,synth,mindratio,minobs,gdist,fixscale,noise, & fixsite,rangea,mag_par_file) if (synth) write(*,*) ' synthetic mode ' c g1 scale magsynth(1)=1.1 magsynth(2)=.00189 magsynth(3)=ref(3)-alog10(ref(2)/2080.*1.E6) & -magsynth(1)*alog10(ref(1))-magsynth(2)*ref(1) c g2 scale magsynth(4)=1.1 magsynth(5)=.00189 magsynth(6)=ref(3)-alog10(ref(2)/2080.*1.E6) & -magsynth(4)*alog10(ref(1))-magsynth(5)*ref(1) c g3 scale magsynth(7)=1.1 magsynth(8)=.00189 magsynth(9)=ref(3)-alog10(ref(2)/2080.*1.E6) & -magsynth(7)*alog10(ref(1))-magsynth(8)*ref(1) idum=-1 c write(*,*) magsynth c stop c write(6,*) 'Orientations to use' c write(6,*) '0. for horiz and vert' c write(6,*) '1. for horiz only' c write(6,*) '2. for vert only' c write(6,*) c write(6,*) 'Give choice' c read(*,*) orient c Define arrays containing options for orientations c If orient=0 then don't do anything c If orient=1. then sz,bz,hz,az,ez are unacceptable write(*,*) if (orient.eq.0.) then write(*,*) 'Horiz and vert data' endif disgard(1)=' ' disgard(2)=' ' disgard(3)=' ' if (orient.eq.1.) then write(*,*) 'Horizontal data only' disgard(1)='Z' c disgard_1(2)='HZ' c disgard_1(3)='AZ' c disgard_1(4)='BZ' c disgard_1(5)='EZ' endif if (orient.eq.2) then write(*,*) 'Vertical data only' disgard(1)='E' disgard(2)='N' c disgard_2(3)='BE' c disgard_2(4)='BN' c disgard_2(5)='HE' c disgard_2(6)='HN' c disgard_2(7)='EE' c disgard_2(8)='EN' c disgard_2(9)='AE' c disgard_2(10)='AN' endif c c----------------------------------------------------------------- c Loop to read events start here c----------------------------------------------------------------- c nobs=0 nstat=0 nev=0 50 continue c c read all parameters for one event from file unit 1 c call rea_event_in(1,all,data,code) c c check if end of file (code=1), if so jump out of loop c if(code.eq.1) goto 1000 c c slsa10 c nstat_ev=0 nevent=nevent+1 ! count events year(nevent)=hyp_year(1) month(nevent)=hyp_month(1) day(nevent)=hyp_day(1) hour(nevent)=hyp_hour(1) min(nevent)=hyp_min(1) sec(nevent)=hyp_sec(1) lat(nevent)=hyp_lat(1) lon(nevent)=hyp_lon(1) mag(nevent)=hyp_mag(1,1) depth(nevent)=hyp_depth(1) call timsec(year(nevent),month(nevent),day(nevent), & hour(nevent),min(nevent),sec(nevent),msec) c c write the whole first header line c c write(6,'(a)') data(1)(1:79) write(3,101) year(nevent),month(nevent),day(nevent),hour(nevent), & min(nevent),sec(nevent),lat(nevent),lon(nevent),mag(nevent), & depth(nevent),nevent no=0 mind=9999. maxd=0. do i=1,rea_nphase !loop through phase lines flag(i)=0 if( rea_phase(i)(1:1).eq.' ' * .or.rea_phase(i)(1:1).eq.'L' * .or.rea_phase(i)(1:2).eq.'S ' * .or.rea_phase(i)(1:2).eq.'Sg' * .or.rea_phase(i)(1:2).eq.'SG' * .or.rea_phase(i)(1:4).eq.'AMPL' * .or.rea_phase(i)(1:3).eq.'AML' * .or.rea_phase(i)(1:4).eq.'AMP ' * .or.rea_phase(i)(1:4).eq.'IAML') then flag(i)=1 if (rea_amp(i).le.0.) then flag(i)=0 endif if (rea_per(i).ge.5.0) then ! don't include data with flag(i)=0 ! period ge 5.0 sec write(*,*) ' ignoring period ge 5 ' endif if (rea_dist(i).le.min_dist) then flag(i)=0 endif c c compute group velocity c vg(i)=rea_dist(i)/(rea_abs_time(i)-msec) if (vg(i).le.1.) then flag(i)=0 write(*,*) ' ignoring vg le 1 ' endif if (rea_dist(i).gt.max_dist) then flag(i)=0 endif if (rea_dist(i).gt.2500.0) then write(6,*) rea_stat(i),' ********************** warning' write(6,*) data(1)(1:79) write(6,*)' This stat/event has a distance of: ', * rea_dist(i) endif if (rea_dist(i).le.0.) then flag(i)=0 endif if (orient.eq.1..or.orient.eq.2.) then if (rea_co(i)(2:2).ne.' ') then do j=1,3 if (rea_co(i)(2:2).eq.disgard(j)) then flag(i)=0 write(*,*) ' disgard comp ',rea_co(i) endif enddo endif endif do j=1,n_stat_array if (rea_stat(i).eq.stat_array(j)) then write(*,*) rea_stat(i),stat_array(j) flag(i)=0 endif enddo do j=1,n_comp_array if (rea_co(i).eq.comp_array(j)) then flag(i)=0 endif enddo c c remove amplitudes that are too large for distance c do j=1,rea_nphase !loop through phase lines if (rea_phase(i).eq.'AML'.and.rea_phase(j).eq.'AML')then if (rea_amp(i).gt.0..and.rea_amp(j).gt.0) then if (rea_dist(i).gt.rea_dist(j)) then if (rea_amp(i).gt.rea_amp(j)*10.) then flag(i)=0 c write(*,*) ' amplitude too large for distance ' write(6,*) ' amplitude too large for distance ', data(1)(1:10),' ' & ,rea_stat(i), rea_co(i),rea_dist(i),' ', & rea_amp(i),' ',rea_stat(j),rea_co(j),' ',rea_dist(j), & ' ',rea_amp(j) endif endif endif endif enddo do j=1,i-1 if (rea_stat(i).eq.rea_stat(j).and. & rea_co(i)(2:2).eq.rea_co(j)(2:2).and. & flag(j).eq.1) then flag(i)=0 write(*,*) ' duplicate component ' endif enddo if (flag(i).eq.1) then c count number of aml used no=no+1 c get min and max distance if (rea_dist(i).gt.maxd) maxd=rea_dist(i) if (rea_dist(i).lt.mind) mind=rea_dist(i) tempamp(no)=alog10(rea_amp(i)) tempdis(no)=sqrt(rea_dist(i)**2+hyp_depth(1)**2) templdis(no)=alog10(tempdis(no)) c slsa10 - count number of stations used c nstat_ev=nstat_ev+1 endif endif enddo if (no.lt.minobs) then write(6,*) ' ignoring event, too few amplitudes ', no, & ' ', data(1)(1:10) goto 50 endif c min_distrange=(max_dist-min_dist)/10. c min_distrange=(max_dist-min_dist)/5. min_distrange=(max_dist-min_dist)*mindratio if (maxd-mind.lt.min_distrange) then write(*,*) ' ignoring event, distance range too short ', & (maxd-mind) goto 50 endif c slsa10 c write(6,*) 'nstat_ev =',nstat_ev c c check that lsq fit shows decay of log(amplitudes) c call lsqlin(no,tempdis,tempamp,a,b,corr,rms) if (b.gt.0.) then goto 50 write(*,*) ' amplitudes not decaying with distance ', & a,b,corr endif call lsqlin(no,templdis,tempamp,a,b,corr,rms) if (b.gt.0.) then goto 50 write(*,*) ' amplitudes not decaying with distance ', & a,b,corr endif call rea_event_out(20,.true.,data,code) c slsa - 19/7/10 write out lat and lons of events used in inversion write(22,*) lat(nevent), lon(nevent), depth(nevent) do i=1,rea_nphase !loop through phase lines if (flag(i).eq.1) then nobs=nobs+1 ! increment nobs for each point dist(nobs)=sqrt(rea_dist(i)**2+depth(nevent)**2) ! hypocentral distance c if (dist(nobs).eq.ref(1)) c & dist(nobs)=dist(nobs)+.1 c log_dist(nobs)=log10(rea_dist(i)) c log_dist(nobs)=alog10(dist(nobs)) stat(nobs)=rea_stat(i) co(nobs)=rea_co(i) if (synth) then if (dist(nobs).le.gdist(1)) then c g1 amp(nobs)=mag(nevent) & -magsynth(1)*alog10(dist(nobs)) & -magsynth(2)*dist(nobs)-magsynth(3) elseif (dist(nobs).le.gdist(2)) then c g2 amp(nobs)=mag(nevent) & -magsynth(4)*alog10(dist(nobs)) & -magsynth(5)*dist(nobs)-magsynth(6) else c g3 amp(nobs)=mag(nevent) & -magsynth(7)*alog10(dist(nobs)) & -magsynth(8)*dist(nobs)-magsynth(9) endif amp(nobs)=10.**amp(nobs) c c add noise c if (noise.gt.0.) then c write(77,*) amp(nobs) write(77,*) amp(nobs) amp(nobs)=amp(nobs)+amp(nobs)*(ran2(idum)-.5)*2.*noise c write(77,*) amp(nobs) write(77,*) amp(nobs) write(77,*) endif else amp(nobs)=rea_amp(i) endif c event_id is assigned for each observation and corresponds to nevent. It is not c connected to a date and time event_id(nobs)=nevent call make_ev_id(year(nevent),month(nevent), & day(nevent),hour(nevent),min(nevent),sec(nevent),evid) evid_arr(nobs)=evid c Compile station list (this is a list of all the stat/comps for which corrections c will be determined in the inversion) c Also compile list of event ids (these are the events for which source terms will c be determined in the inversion) c if (nobs.eq.1) then c nstat=1 c nev=1 c stat_list and co_list are lists of all the stations and components in the dataset c stat_list(nstat)=stat(nobs) c co_list(nstat)=co(nobs) c event_list is a list of integers which represent events c I'm not sure if I need to keep event_list(nev)=event_id(nobs) c I think I can just use evid(nev) instead. c Make an event ID from the year,month,day,hour,minute and second. This is the c index that will hopefully glue everything together.... I'm writing the ID directly c into an array called evid, which has 5000 elements. c c evid_list(nev)=evid_arr(nobs) c endif c if (nobs.gt.1) then flag_stat=1 do j=1,nstat if (stat(nobs).eq.stat_list(j)) then c lo 16/4/2010 if (co(nobs).eq.co_list(j)) then flag_stat=0 c endif endif enddo if (flag_stat.eq.1) then nstat=nstat+1 stat_list(nstat)=stat(nobs) c lo 16/4/2010 co_list(nstat)=co(nobs) endif flag_ev=1 do j=1,nev if (evid_arr(nobs).eq.evid_list(j)) then flag_ev=0 c if (event_id(nobs).eq.event_list(j)) then c flag_ev=0 endif enddo if (flag_ev.eq.1) then nev=nev+1 evid_list(nev)=evid_arr(nobs) endif per(nobs)=rea_per(i) phase(nobs)=rea_phase(i) vg(nobs)=vg(i) if (per(nobs).lt.5.) ! lo 10 oct 2012 & write(2,100) amp(nobs),dist(nobs),stat(nobs),co(nobs), & per(nobs),phase(nobs),event_id(nobs),vg(nobs) write(4,*) amp(nobs),dist(nobs) c slsa - 19/7/10 write(23,*) lat(nevent),lon(nevent) call stat_loc(stat(nobs),'0',slat,slon,stel) write(23,*) slat,slon write(23,*) '>' endif enddo c c get next event c write(4,'(a)') '>' goto 50 c c end of file c 1000 continue c c Note the indices in event_list, stat_list and co_list for each observation c This is the master index to which all subsequent simulated datasets c should be linked do i=1,nobs do j=1,nev if (evid_arr(i).eq.evid_list(j)) then ev_index(i)=j endif enddo enddo do i=1,nobs do j=1,nstat if (stat(i).eq.stat_list(j)) then c lo 16/4/2010 if (co(i).eq.co_list(j)) then stat_index(i)=j c count number of observations per station stat_nobs(j)=stat_nobs(j)+1 c endif endif enddo enddo do i=1,nstat c write(8,*) stat_list(i),co_list(i) write(8,*) stat_list(i),stat_nobs(i) call stat_loc(stat_list(i),'0',slat,slon,stel) write(28,*) slat,slon enddo write(6,*) ! blank line close(2) ! close output file close(3) close(4) close(8) close(11) close(20) c slsa close(28) close(22) close(23) write(6,*) 'Number of events in input file', nevent write(6,*) 'Number of events used ', nev write(6,*) 'Number of observations is ', nobs write(6,*) 'Number of stations is ',nstat if (rangea(1).ne.0..and.rangea(2).ne.0.and. & rangea(3).ne.0.) then open(15,file='mag2_range_a.out',status='unknown') fixscale(1)=rangea(1) do while(fixscale(1).le.rangea(2)) write(6,*) ' SVD inversion a= ',fixscale(1) call svd_ml(nobs,amp,dist,stat_index,ev_index,ref, & magscale,magscale_err,siteterm,siteterm_err,sourceterm, & sourceterm_err,gdist,fixscale,fixsite,err) c c evaluate result c magscale(6)=ref(3)-alog10(ref(2)/wagain*1E6) & -magscale(1)*alog10(ref(1))-magscale(4)*ref(1) call mag2var(nobs,amp,dist,nev,ev_index, & nstat,stat_index,siteterm,magscale, & var,sigma) c c output c write(15,*) ' values ',fixscale(1),var,sigma do k=1,7 write(15,'(a,f4.2,a,i1,a,f8.5)') & ' a=',fixscale(1),' magscale(',k,')=',magscale(k) enddo fixscale(1)=fixscale(1)+rangea(3) enddo close(15) else if (inv_type.eq.1.) then write(6,*) 'SVD inversion' call svd_ml(nobs,amp,dist,stat_index,ev_index,ref, & magscale,magscale_err,siteterm,siteterm_err,sourceterm, & sourceterm_err,gdist,fixscale,fixsite,err) endif c if (inv_type.eq.2) then c write(*,*) 'LSQR inversion' c call lsqr_inv(dist,amp,ev_index,stat_index,nobs,nev, c & nstat,ref,magscale,siteterm,sourceterm) c endif do k=1,7 write(*,'(a,i1,a,f8.5)') ' magscale(',k,')=',magscale(k) enddo endif c c evaluate result c magscale(6)=ref(3)-alog10(ref(2)/wagain*1E6) & -magscale(1)*alog10(ref(1))-magscale(4)*ref(1) call mag2var(nobs,amp,dist,nev,ev_index, & nstat,stat_index,siteterm,magscale, & var,sigma) write(*,*) ' variance and standard deviation: ', & var,sigma c Write out results to mag2.out 9999 continue call write_mag2(inv_type,nev,nstat,nobs,magscale, & magscale_err,siteterm,siteterm_err,sourceterm,sourceterm_err, & ref,stat_list,evid_list,gdist) c & ref,stat_list,co_list,evid_list) 100 format(f8.1,x,f8.1,x,a5,x,a2,f6.2,x,a8,x,i4,x,f6.2) 101 format(i4,x,i2,x,i2,x,i2,x,i2,x,f5.2,x,f6.2,x,f7.2,x, & f3.1,x,f4.1,x,i5) 102 format(a2) 103 format(a1,a1) 104 format(a9,f8.6,a5,f8.6) 105 format(a13,i3,x,f6.3,a5,f6.4) 106 format(a10,i5,x,a5,x,f4.2,a5,f5.3) stop end c------------------------------------------------------------- c Subroutines subroutine get_mag2_par(inv_type,orient,min_dist,max_dist, & comp_array,ncomp,stat_array,nstat,ref,boot_ninv,synth, & mindratio,minobs,gdist,fixscale,noise,fixsite,rangea, & mag_par_file) implicit none include 'libsei.inc' ! Seisan definitions include 'seidim.inc' ! ------------------ include 'seisan.inc' ! ------------------ integer read1 ! input unit integer code character*80 line real var real inv_type real min_dist real max_dist real orient real ref(*) real noise integer ncomp integer nstat integer seiclen character*2 comp_array(*) character*5 stat_array(*) character*2 var_char real fixscale(*),fixsite,rangea(*) integer boot_ninv logical synth integer minobs real mindratio real gdist(*) character*(*) mag_par_file c Up to 10 components can be set to be ignored (2 letter code) c Up to 20 stations can be set to be ignored (up to 4 letter code) call sei get file( open$+ignore$, ! Open waveform file. & read1, ! On unit. & code, ! Returned condition. & 'DAT', ! Alternative search director & mag_par_file) ! For this filename. if(code.ne.e_ok$) then write(*,*) ' definition file does not exist: mag2.par' stop endif ncomp=0 synth=.false. mindratio=.5 gdist(1)=0. gdist(2)=0. fixscale(1)=0. fixscale(2)=0. fixscale(3)=0. fixscale(4)=0. noise=0. 10 continue read(read1,'(a80)',end=999) line c write(*,*) line if (line(1:14).eq.'INVERSION TYPE') then read(line(41:50),*) var if (var.ne.0.) then inv_type=var endif elseif (line(1:9).eq.'BOOTSTRAP') then read(line(41:50),*) var if (var.ne.0) then boot_ninv=var endif elseif (line(1:9).eq.'SYNTHETIC') then read(line(41:50),*) var if (var.ne.0) then synth=.true. endif elseif (line(1:9).eq.'DISTANCES') then read(line(41:50),*) var if (var.ne.0.) then min_dist=var endif read(line(51:60),*) var if (var.ne.0.) then if (var.gt.min_dist) then max_dist=var endif endif elseif (line(1:11).eq.'ORIENTATION') then read(line(41:50),*) var orient=var elseif (line(1:5).eq.'NOISE') then read(line(41:50),*) noise elseif (line(1:18).eq.'REFERENCE DISTANCE') then read(line(41:50),*) var ref(1)=var elseif (line(1:14).eq.'SCALE DISTANCE') then read(line(41:50),*) var gdist(1)=var read(line(51:60),*) var gdist(2)=var elseif (line(1:11).eq.'FIX SCALE A') then read(line(41:50),*) fixscale(1) read(line(51:60),*) fixscale(2) read(line(61:70),*) fixscale(3) elseif (line(1:7).eq.'RANGE A') then read(line(41:50),*) rangea(1) read(line(51:60),*) rangea(2) read(line(61:70),*) rangea(3) elseif (line(1:11).eq.'FIX SCALE B') then read(line(41:50),*) fixscale(4) elseif (line(1:8).eq.'FIX SITE') then read(line(41:50),*) fixsite elseif (line(1:19).eq.'REFERENCE AMPLITUDE') then read(line(41:50),*) var ref(2)=var elseif (line(1:19).eq.'REFERENCE MAGNITUDE') then read(line(41:50),*) var ref(3)=var elseif (line(1:11).eq.'IGNORE COMP') then if (seiclen(line(41:50)).ne.0) then ncomp=ncomp+1 comp_array(ncomp)=line(41:50) endif elseif (line(1:11).eq.'IGNORE STAT') then if (seiclen(line(41:50)).ne.0) then nstat=nstat+1 stat_array(nstat)=line(41:50) endif c slsa10 elseif (line(1:27).eq.'MINIMUM NUMBER OF OBS/EVENT') then read(line(41:50),*) var minobs=var elseif (line(1:23).eq.'MIN DISTANCERANGE RATIO') then read(line(41:50),*) var mindratio=var endif goto 10 999 continue call sei close(close$,read1,code) if (minobs.lt.2) minobs=2 ! lo 16 May 2012 return end c-------------------------------------------------------------------- subroutine write_mag2(inv_type,nev,nstat,nobs,magscale, & magscale_err,siteterm,siteterm_err,sourceterm,sourceterm_err, & ref,stat_list,evid_list,gdist) c & ref,stat_list,co_list,evid_list) implicit none c Subroutine to write output file mag2.out c Declare variables integer i real inv_type integer nev integer nstat integer nobs real magscale(*), magscale_err(*) real siteterm(*), siteterm_err(*) real sourceterm(*), sourceterm_err(*) real c(3) character*16 evid_list(*) character*5 stat_list(*) c character*2 co_list(*) real ref(*) real slat,slon,latm,lonm real elev integer latd,lond character*1 latc,lonc integer seiclen real gdist(*) real siteav real wagain wagain=2080. open(13,file='mag2.out',status='unknown') open(14,file='mag2_station_hyp.out',status='unknown') c write(*,*) inv_type write(13,*) 'ML inversion output' write(13,*) if (inv_type.eq.1) then write(13,*) 'SVD inversion' endif write(13,*) write(13,*) 'Total number of events: ',nev write(13,*) 'Total number of stations: ',nstat write(13,*) 'Total number of observations: ',nobs write(13,*) write(13,*) 'Reference distance = ', ref(1) write(13,*) 'Reference amplitude = ', ref(2) write(13,*) 'Reference magnitude = ', ref(3) write(13,*) c write(13,'(a)') ' Ml = log A + a log(dist/refdist) '// c & '+ b (dist-refdist) + c + S ' c write(13,104) 'a1= ', magscale(1), ' +/- ',magscale_err(1) c write(13,104) 'a2= ', magscale(2), ' +/- ',magscale_err(2) c write(13,104) 'a3= ', magscale(3), ' +/- ',magscale_err(3) c write(13,104) 'b = ', magscale(4), ' +/- ',magscale_err(4) c write(13,104) 'c1= ', magscale(5), ' +/- ',magscale_err(5) c write(13,104) 'c2= ', magscale(6), ' +/- ',magscale_err(6) c write(13,104) 'c3= ', magscale(7), ' +/- ',magscale_err(7) c c(1)=magscale(5)-magscale(1)*alog10(ref(1))-magscale(4)*ref(1) c c(2)=magscale(6)-magscale(2)*alog10(ref(1))-magscale(4)*ref(1) c c(3)=magscale(7)-magscale(3)*alog10(ref(1))-magscale(4)*ref(1) c(1)=ref(3)-alog10(ref(2)/wagain*1E6) & -magscale(1)*alog10(ref(1))-magscale(4)*ref(1) c(2)=ref(3)-alog10(ref(2)/wagain*1E6) & -magscale(2)*alog10(ref(1))-magscale(4)*ref(1) c(3)=ref(3)-alog10(ref(2)/wagain*1E6) & -magscale(7)-magscale(3)*alog10(ref(1))-magscale(4)*ref(1) write(13,*) write(13,'(a)') ' Ml = log A + a log(dist) '// & '+ b (dist) + c + S ' write(13,104) 'a1 = ', magscale(1), ' +/- ',magscale_err(1) write(13,104) 'a2 = ', magscale(2), ' +/- ',magscale_err(2) write(13,104) 'a3 = ', magscale(3), ' +/- ',magscale_err(3) write(13,104) 'b = ', magscale(4), ' +/- ',magscale_err(4) write(13,104) 'c1 = ',c(1) write(13,104) 'c2 = ',c(2) write(13,104) 'c3 = ',c(3) write(13,*) write(6,*) write(6,'(a)') ' Ml = log A + a log(dist) '// & '+ b (dist) + c + S ' write(6,104) 'a1 = ', magscale(1), ' +/- ',magscale_err(1) write(6,104) 'a2 = ', magscale(2), ' +/- ',magscale_err(2) write(6,104) 'a3 = ', magscale(3), ' +/- ',magscale_err(3) write(6,104) 'b = ', magscale(4), ' +/- ',magscale_err(4) write(6,104) 'c1 = ',c(1) write(6,104) 'c2 = ',c(2) write(6,104) 'c3 = ',c(3) c c write mag scale to station file c write(14,'(a,f5.3)') 'RESET TEST(76)=',magscale(1) write(14,'(a,f7.5)') 'RESET TEST(77)=',magscale(4) write(14,'(a,f5.2)') 'RESET TEST(78)=',c(1) write(14,'(a,f5.3)') 'RESET TEST(176)=',magscale(2) write(14,'(a,f7.5)') 'RESET TEST(177)=',magscale(4) write(14,'(a,f5.2)') 'RESET TEST(178)=',c(2) write(14,'(a,f6.1)') 'RESET TEST(179)=',gdist(1) write(14,'(a,f5.3)') 'RESET TEST(180)=',magscale(3) write(14,'(a,f7.5)') 'RESET TEST(181)=',magscale(4) write(14,'(a,f5.2)') 'RESET TEST(182)=',c(3) write(14,'(a,f6.1)') 'RESET TEST(183)=',gdist(2) write(14,'(a)') ' ' siteav=0. do i=1,nstat c c write out residuals c call stat_loc(stat_list(i),' ',slat,slon,elev) write(13,105) 'Station # ',i, stat_list(i), c & co_list(i),siteterm(i),' +/- ', siteterm_err(i), & siteterm(i),' +/- ', siteterm_err(i), & slat,slon siteav=siteav+siteterm(i)/nstat c write(*,*) ' siteav ',siteav c c write stations with location and residual to seisan station file format c latc='N' lonc='E' if (slat.lt.0.) latc='S' if (slon.lt.0.) lonc='W' slat=abs(slat) slon=abs(slon) latd=int(slat) lond=int(slon) latm=(slat-float(latd))*60. lonm=(slon-float(lond))*60. if (seiclen(stat_list(i)).lt.5) then write(14,'(2x,a4,i2,f5.2,a1,i3,f5.2,a1,i4,6x,f5.2)') & stat_list(i),latd,latm,latc,lond,lonm,lonc,int(elev), & siteterm(i) else write(14,'(1x,a5,i2,f5.2,a1,i3,f5.2,a1,i4,6x,f5.2)') & stat_list(i),latd,latm,latc,lond,lonm,lonc,int(elev), & siteterm(i) endif enddo write(13,'(a,1x,f5.2)') ' Average site term: ',siteav do i=1,nev write(13,106) 'Event # ', i,evid_list(i), ' ML = ', c put in proper magnitude, xxx & sourceterm(i),' +/- ', sourceterm_err(i) enddo close(13) close(14) c Format statements 104 format(a9,f8.5,a5,f8.5) c105 format(a10,x,i3,x,a5,x,a2,x,f6.3,a5,x,f6.4,x,f8.3,x,f8.3) 105 format(a10,x,i3,x,a5,x,f6.3,a5,x,f6.4,x,f8.3,x,f8.3) 106 format(a10,i5,x,a16,x,a5,x,f5.2,a5,f5.3) return end c----------------------------------------------------------------------- function ran2(idum) implicit none c------------------------------------------------------------ c Long period (> 2 x 10^18) random number generator of c L'Ecuyer with Bays-Durham shuffle and added safeguards. c Returns a uniform random deviate between 0.0 and 1.0 c (exclusive of the endpoint values). Call with idum a -ve c number to initialise; thereafter, do not alter idum c between successive deviates in a sequence. RNMZ should c approximate the largest floating value that is less than c 1 c c Numerical Recipes in Fortran 77: The Art of Scientific c Computing (1986-1992) c------------------------------------------------------------ integer idum integer im1,im2 integer imm1 integer ia1,ia2 integer iq1,iq2 integer ir1,ir2 integer ntab integer ndiv real ran2 real am real eps real rnmx parameter (im1=2147483563,im2=2147483399,am=1./im1, & imm1=im1-1,ia1=40014,ia2=40692,iq1=53668,iq2=52774, & ir1=12211,ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7, & rnmx=1.-eps) integer idum2 integer j,k integer iv(ntab) integer iy save iv,iy,idum2 data idum2/123456789/, iv/ntab*0/, iy/0/ c Initialise if (idum.le.0) then c Be sure to prevent idum=0 idum=max(-idum,1) idum2=idum c Load the shuffle table (after 8 warm ups) do j=ntab+8,1,-1 k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 if (idum.lt.0) idum=idum+im1 if (j.le.ntab) iv(j)=idum enddo iy=iv(1) endif c Start here when not initialising k=idum/iq1 c Compute idum=mod(ia1*idum,im1) without overflows by Schrage's method idum=ia1*(idum-k*iq1)-k*ir1 if (idum.lt.0) idum=idum+im1 k=idum2/iq2 c Compute idum2=mod(ia2*idum2,im2) likewise idum2=ia2*(idum2-k*iq2)-k*ir2 if (idum2.lt.0) idum2=idum2+im2 j=1+iy/ndiv ! will be in the range 1:ntab c Here idum is shuffled, idum and idum2 are combined to generate c output iy=iv(j)-idum2 iv(j)=idum if (iy.lt.1) iy=iy+imm1 c Because users don't expect endpoint values ran2=min(am*iy,rnmx) return end subroutine mag2var(nobs,amp,dist,nev,ev_index, & nstat,stat_index,siteterm,magscale, & var,sigma) implicit none integer nobs,nev,ev_index(*),stat_index(*) real amp(*),dist(*),magscale(*),siteterm(*) integer i,j,k,evobs,nstat,n real evav,ml(100),var,sigma,stc,evvar var=0. n=0 do i=1,nev evobs=0 evav=0. evvar=0. do j=1,nobs if (ev_index(j).eq.i) then c find station correction stc=siteterm(stat_index(j)) c magnitude evobs=evobs+1 n=n+1 ml(evobs)=a & log10(amp(j))+magscale(1)*alog10(dist(j))+ & magscale(4)*dist(j)+magscale(6)+stc evav=evav+ml(evobs) write(55,*) i,j,stc,ml(evobs),stat_index(j) endif enddo evav=evav/float(evobs) do j=1,evobs var=var+(ml(j)-evav)**2 evvar=evvar+(ml(j)-evav)**2 enddo write(55,*) i,evav,evobs,var,sqrt(evvar/float(evobs)) enddo var=var/n sigma=sqrt(var) return end c-------------------------------------------------------------------- subroutine make_ev_id(year,month,day,hour,min,sec,evid) implicit none c Routine to construct a 16 character ID for each event by c concatenating the year, month, day, hour, minute and second character*4 string1 character*2 string2 character*2 string3 character*2 string4 character*2 string5 character*4 string6 character*16 evid integer year integer month integer day integer hour integer min real sec integer sec_mul write(unit=string1,fmt='(i4.4)') year write(unit=string2,fmt='(i2.2)') month write(unit=string3,fmt='(i2.2)') day write(unit=string4,fmt='(i2.2)') hour write(unit=string5,fmt='(i2.2)') min sec_mul=sec*100 write(unit=string6,fmt='(i4.4)') sec_mul evid=string1//string2//string3//string4//string5//string6 c return end c------------------------------------------------------------------ subroutine standev(data,n,stdev,boot_ninv,average) implicit none c Declare variables real sum integer count integer n integer i integer j integer jstop real xi real xbar real stdev(*) integer maxinv parameter (maxinv=1000) real data(maxinv,maxinv) real average(n) integer boot_ninv jstop=boot_ninv+1 do i=1,n sum=0 count=0 do j=1,jstop if (data(i,j).ne.0.) then count=count+1 xi=data(i,j) xbar=average(i) sum=sum+((xi-xbar)**2) endif enddo stdev(i)=((1/count)*sum)**0.5 enddo return end