c****************************************************************************** C c c latest update: c c oct 17 91 by jh : increase field for ml to f4.1 in routine updatml c jan 7 92 : add routine read_stat_mod c jun 21 92 jh,cl: bug in ml calc. alog instead of alog10 c new global location changes. c jul 30 : bug in hypout c oct 9,12 92 : update_ml also calculates mc magnitude, posible c to calculate if no location if distance is there. c calculate ms c oct 13 : small changes, updateml is now called update_mag c oct 26 : check for missing stations in print.out c nov 22 : use hypocentral distance for ml in update_mag c mar 93 : test(8) neg, use log coda squared in mc c jul 19 93 : iasp to capital letters c nov : new bl program c nov 17 93 : magnitudes bug c april 5 94 : deleted hypout, added dlt and isort to hypocent window c jun 2 : no ml if no dist c oct 14 : limit mb to distances of 11000km CJAB(BGS)Jan95 : Install file & error handling. c jun 2 95 jh : only calculate Ml if balnk, S or L-phase c feb 21 96 : magnitudes on a second line c mar 6 96 : change condition for calculaiton of mb, mc and ml c mar 11 : add dist type to second mag line c mar 14 : small change to clear_old c mar 21 : bug in mb c apr 16 : --------- c sep 19, 96 : bug average corner frequency, more data on 2 mag line c mar 4 97 : increase amplitude files for amp out from 8 to 12 c mar 18 : in connection with above, move coda: right c jan 30 98 : magnitude residuals c feb 18 : spectral p velocity c mar 1 : set svelocity to 3.5 if not defined in older data c sep 18 98 jh : ---------- seisan version 7.0 check -------------- c 5 character stations c oct 28 98 bmt : linux changed, save and * c jun 10 99 jh : add AMP as a valid phase for calculating ml and mb c jul 26 : remove routine mb_att since a single file already c jul 29 : add use_eev to routine read_stat_model c sep 24 lo : put in edepth and spec_phase to upd_spec c jh : fix new spectral distance c sep 27 : replace travel time with kappa in spectral output c feb 17 2000 jh : suppres output if output false c oct 24 : clear out old H lines in clear routine c oct 26 : bug with spectral average and update: if station c with phase line missing, following spectra not used c nov 5 : small correction to write out, only write number c of spectra if more than 0 c mar 5, 2001 jh : distance as real c may 19 : more accuracy in sdrop and radius c may 27 : do not use spectra not P or S in average c may 31 : Recalculted spectra not written out !!!!! c apr 19 2002 bmt : change maxline=1000 c apr 13 2003 jh : implement AMPL, AMPb, AMPS c may 7 2003 lo : ISC amp codes changed to AML, AMb and AMS c apr 08 2003 lot : add call to iasp91_filename to get comp de c oct 12 2005 jh : only calculate mb if distance larger than test(57) c nov 28 2005 jh : put ml first, then mc (before opposite) c jul 22 2008 lo : added standard deviation to spec aver computation c oct 28 2008 jh : fix output of cornerf c oct 3 2009 jh : new amplitudes for IASP standard, ncrease distance c mb from 11 000 to 11120 km (100 deg). for mb change c upper period limit form 5 to 3 s c jan 5 2010 jh : also write out period on screen c jan 8 2010 jh : put in distance and period limits for MB c feb 28 2010 jh : pu in new iaspei amplitude phase names c nov 5 2010 lo : allow for 5 characters in SPEC line c oct 14 2011 lo : changed .35 to .37 as in mulplt for source radius calculation c jan 7 2013 jh : comment out seidim.inc since in hypparm.inc c apr 5 2013 jh : make sure no spec update of stations with no readings c apr 25 2013 jh : write out magnitudes in file c jan 30 2014 jh : fix routine spec, error with 5 letter stat, not c used in average c mar 28 2014 jh : do not read VS in spec line, replaced by q_below_1hz c dec 09 2015 jh : new test variables to limit distance in mb, Ms, c limit depth for Ms and MS. Put in far distance c limit for MB, was set wrong. different clean up, c maxline back to 20 (why 1000 ?), fix all error c messages when no magnitude calcualated even if c amplitude available c dec 22 2015 jh: add option for median magnitude, also mw c jan 13 2016 jh: allow no period for mB c feb 5 2016 jh: ensure that AMPG is not used for Ms c feb 19 2016 jh: return if no station file, do not stop c mar 12 2016 jh: add nmessage and message to hypocenter call c may 22 2016 jh: get maxline from hyp common block, do not c define locally c C################################################################# SUBROUTINE UPDATE_MAG(DATA,NHEAD,NRECORD,AGENCY,TEST,OUTPUT) C C ROUTINE CALCULATES AVERAGE MAGNITUDES FOR EVENT WRITTEN C IN STANDARD FORMAT IN TEXT ARRAY DATA. MAGNITUDES ARE c WRITTEN IN HEADER C C INPUT: DATA: EVENT IN STANDARD FORMAT IN TEXT ARRAY C NHEAD: NUMBER OF HEADER RECORDS IN DATA C NRECORD: TOTAL NUMBER OF RECORDS IN DATA C AGENCY: TYPE AND AGENCY E.G. BER C TEST : RESET TEST FROM HYPOCENT c OUTPUT : if true, write to screen C C OUTPUT: DATA: EVENT IN STANDARD FORMAT IN TEXT ARRAY C WITH AVERAGE LOCAL MAGNITUDE IN HEADER C C CHANGES C JAN 91 RAN CHANGED TO PC c aug 7 91 jh : bugs c OCT 10 91 BY J.H.: include local magnitude using the following c ml = A*log(amp) + B*log(dist) + C * dist + D c where amp is max ground motion o-p in nm, C dist is hypocentral distance in km c oct 18 : add constant d c october 92 : mc and ms calculation c nov 22 : use hypocentral distance for ml c nov 17 93 : do not overwrite 3. magnitude c apr 94 : new hypocent c apr 94 : use reset test for magnitudes c jun 94 : add mb c jul 94 : distance bug c nov 95 : error when writing out ms magntude c feb 21 96 : magnitudes on a second line c jan 28 98 jh : use hypocentral distance also for coda c jan 30 : magnitude residuals c apr 11 03 jh : implement amplitude phases AMPL, AMPB and AMPS, c implicit none save c include 'hypparm.inc' include 'rea.inc' integer nhead,nrecord c BRL 6/98: changed a,b,c,d to am,bm,cm,dm real am,bm,cm,dm logical output character*3 agency character*80 text c-- amplitude,period, distance, hypocentral distance and coda real amp,period,dist,hdist,coda c-- depth real depth c--number of stations with different magnitudes integer nml,nmc,nms,nmss,nmb,nmbb c-- help varibles integer nrec,icoda,nline,maxline real q,cmag1 real xmag(1000) ! ml real cmag(1000) ! coda magnitudes real smag(1000) ! surface wave magnitude (20 s) real ssmag(1000) ! broad band surface wave magnitude real bmag(1000) ! mn real bbmag(1000) ! mB real av_cmag ! average cmag real av_xmag ! average ml real av_smag ! average ms real av_ssmag ! average mS real av_bmag ! average of mb real av_bbmag ! average of mB logical mag_write ! true if a posiiton found for writing mag. C--data c CHARACTER*80 DATA(*) c--single local magnitude REAL MAGNITUDE c--reset test variables c real test(*) character*5 station(nstats) ! station list character*5 stat ! one station real mag_cor(5,nstats) ! magnitude residuals,1:C, 2:L, 3:B, 4:S 5:W integer kstat ! station number currently used c--counters etc integer i,k common/hyp/station,mag_cor,nline,maxline c nmc=0 nms=0 nmss=0 nmb=0 nmbb=0 nms=0.0 av_cmag=0.0 av_smag=0.0 av_ssmag=0.0 av_bmag=0.0 bmag=0.0 av_bbmag=0.0 nml=0 av_xmag=0.0 nline=0 c maxline=20 c maxline=1000 c c get ml parameters c am=test(75) bm=test(76) cm=test(77) dm=test(78) c c read depth c read(data(1)(39:43),'(f5.0)') depth c DO 1 NREC=NHEAD+1,NRECORD-1 READ(DATA(NREC),'(1x,a5,23X,i4,g7.1,1X,F4.1,25X,f5.0)' * ,err=1234) * stat,icoda,AMP,PERIOD,DIST c c if amp or coda is not zero but distance zero, try to find distance at c other line with same station if no location flag set but c magntudes calculated, test 106 c if((amp.gt.0.0.or.icoda.gt.0).and.dist.eq.0.0.and. * test(106).gt.0.0) then do i=nhead+1,nrecord-1 if(stat.eq.data(i)(2:6).and.data(i)(71:75).ne.' ') * then read(data(i)(71:75),'(f5.0)') dist goto 1235 endif enddo endif goto 1235 1234 continue write(6,*)' Something wrong with input record:' write(6,'(a)') data(nrec) amp=0.0 1235 continue c c find station number c do i=1,nstats if(station(i).eq.stat) then kstat=i goto 1245 endif enddo 1245 continue c write(6,*)kstat,mag_cor(1,kstat) c c c use hypocentral distance c hdist=sqrt(dist*dist+depth*depth) c c must be larger than 1 km c if(hdist.le.1.0) hdist=1.0 ! added april 2003 by jh c coda=icoda c c check if coda and distance available, use log coda squared if c coefficient is negative, if too far, do not calculate c if(dist.gt.0.and.icoda.gt.0) then if(dist.le.test(57)) then if(test(8).gt.0.0) then cmag1=test(7)+test(8)*alog10(coda)+hdist*test(9) else cmag1=test(7)-test(8)*alog10(coda)*alog10(coda) * +hdist*test(9) endif if(mag_cor(1,kstat).ne.99.0) then nmc=nmc+1 cmag(nmc)=cmag1+mag_cor(1,kstat) av_cmag=av_cmag+cmag(nmc) endif c c write distance, coda and mag c if(nline.ge.0.and.output) then if(mag_cor(1,kstat).ne.99.0) then if(mag_cor(1,kstat).eq.0.0) then write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,1x,f8.1,3x, * 9x,a,f5.1)') * data(nrec)(2:8),' hdist:',hdist, * 'coda: ',coda,'mc = ',cmag1 else write ! write correction * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,1x, * f8.1,3x,9x,a,f5.1,8x,a,f5.2)') * data(nrec)(2:8),' hdist:',hdist, * 'coda: ',coda,'mc = ',cmag1,' + ', * mag_cor(1,kstat) endif endif nline=nline+1 endif if(nline.gt.maxline.and.output) then nline=0 write(6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif else write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Coda available,event too far', * ' magnitude not used ' nline=nline+1 if(nline.gt.maxline.and.output) then nline=0 write(6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif endif endif c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c local magnitude cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c check if amplitude is available, assume local magnitude if c phase is blank, AMPL, L or S, AMP,IAML, period less than 5 secs c and distance less than test(57) c c IF(AMP.GT.0.0.and.dist.gt.0.0.and. * ( data(nrec)(11:11).eq.' ' * .or.data(nrec)(11:11).eq.'L' * .or.data(nrec)(11:12).eq.'S ' * .or.data(nrec)(11:12).eq.'Sg' * .or.data(nrec)(11:12).eq.'SG' * .or.data(nrec)(11:14).eq.'AMPL' * .or.data(nrec)(11:13).eq.'AML' * .or.data(nrec)(11:14).eq.'AMP ' * .or.data(nrec)(11:14).eq.'IAML')) then c c now check more c if(period.le.5.0.and. * dist.gt.0.and.dist.le.test(57))then magnitude=am*alog10(amp)+bm*alog10(hdist)+cm*hdist+dm if(mag_cor(2,kstat).ne.99.0) then nml=nml+1 xmag(nml)=magnitude + mag_cor(2,kstat) av_xmag=av_xmag+xmag(nml) endif c c write distance, amplitude and magnitude, first in a file c hypmag.out c write * (98,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' hdist:',hdist, * 'amp:',amp,' T: ',period,'ml = ',xmag(nml) c c store in rea common block c c rea_mag_type(nrec)='ML' c rea_mag(nrec)=magnitude c add res c write(data5(nrec)(62:),'(a2,f5.2)') 'ML',xmag(nml) if(nline.ge.0.and.output) then if(mag_cor(2,kstat).ne.99.0) then if(mag_cor(2,kstat).eq.0.0) then write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' hdist:',hdist, * 'amp:',amp,' T: ',period,'ml = ',magnitude else ! magnitude correction write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1,a,f5.2)') * data(nrec)(2:8),' hdist:',hdist, * 'amp:',amp,' T: ',period,'ml = ', * magnitude,' + ', * mag_cor(2,kstat) endif endif nline=nline+1 endif else if(nline.ge.0.and.output)then if(period.gt.5.0) then write(6,'(1x,a7,a,a)') data(nrec)(2:8), * ' Amp for ml available, period too large,', * ' amplitude not used' endif if(dist.gt.test(57)) then write(6,'(1x,a7,a,a)') data(nrec)(2:8), * ' Amp for ml available, distance too large,', * ' amplitude not used' endif nline=nline+1 endif endif if(nline.gt.maxline.and.output) then nline=0 write(6,'('' Return to continue, q to end listing '', * $)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif ENDIF c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c surface magnitude ms(20) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c check if surface wave magnitude, assume period gt 10 secs, c any amplitude c if(amp.gt.0.0. * and.(data(nrec)(11:11).eq.' '.or. * data(nrec)(11:14).eq.'AMPS'.or. * data(nrec)(11:14).eq.'AMS '.or. * data(nrec)(11:14).eq.'AMP '.or. * data(nrec)(11:17).eq.'IAMs_20')) then if(period.ge.10.0.and.dist.ge. * test(114).and.dist.le.160.0*111.2.and. * depth.le.test(115)) then magnitude=alog10(amp/(1000.0*period))+ * 1.66*alog10(dist/111.2)+3.3 if(mag_cor(4,kstat).ne.99.0) then nms=nms+1 smag(nms)=magnitude+mag_cor(4,kstat) av_smag=av_smag+smag(nms) endif c c write distance, amplitude and magnitude c if(nline.ge.0.and.output) then if(mag_cor(4,kstat).ne.99.0) then if(mag_cor(4,kstat).eq.0.0) then write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'Ms = ',magnitude else write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1,a,f5.2)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'Ms = ', * magnitude,' + ', * mag_cor(4,kstat) endif endif endif ! for output if(nline.gt.maxline.and.output) then nline=0 write(6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif else ! period etc not ok if(nline.ge.0.and.output) then if(dist.le.test(114).or.dist.gt.160.0*111.2) then write(6,'(1x,a7,a,a)') data(nrec)(2:8), * ' Amp for Ms available, wrong distance,', * ' amplitude not used' endif if(depth.gt.test(115)) then write(6,'(1x,a7,a,a)') data(nrec)(2:8), * ' Amp for Ms available, event too deep,', * ' amplitude not used' endif if(period.lt.10.0) then write(6,'(1x,a7,a,a)') data(nrec)(2:8), * ' Amp for Ms available, period below 10 s,', * ' amplitude not used' endif nline=nline+1 endif endif ! for period etc ok endif ! for Ms c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c surface amgnitude MS ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c check if surface wave magnitude MS c if(amp.gt.0.0. * and.data(nrec)(11:17).eq.'IVMs_BB') then c c ampllitude available c if(period.ge.3.0. * and.period.le.60.0.and.dist.ge.2.0*111.2. * and.dist.le.160.0*111.2.and.depth.le.test(115)) * then magnitude=alog10(amp/(1000.0*6.28))+ * 1.66*alog10(dist/111.2)+3.3 if(mag_cor(4,kstat).ne.99.0) then nmss=nmss+1 ssmag(nmss)=magnitude+mag_cor(4,kstat) av_ssmag=av_ssmag+ssmag(nmss) endif c c write distance, amplitude and magnitude c if(nline.ge.0.and.output) then if(mag_cor(4,kstat).ne.99.0) then if(mag_cor(4,kstat).eq.0.0) then write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'MS = ',magnitude else write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * a,f5.1,3x,a,f5.1,a,f5.2)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'MS = ', * magnitude,' + ', * mag_cor(4,kstat) endif endif nline=nline+1 endif if(nline.gt.maxline.and.output) then nline=0 write(6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif else if(nline.ge.0.and.output) then if(period.lt.3.0.or.period.gt.60.0) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for MS available, wrong period, ', * 'amplitude not used' endif if(dist.lt.2.0*111.2.or.dist.gt.160.0*111.2) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for MS available, wrong distance, ', * 'amplitude not used' endif if(depth.gt.test(115)) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for MS available, event too deep, ', * 'amplitude not used' endif nline=nline+1 if(nline.gt.maxline.and.output) then nline=0 write(6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif endif endif endif c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c mb cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c check if mb, only if period lt 3 secs and greater than 0.2 and c and distance gt test(113), and one of the phases given below c IF(AMP.GT.0.0.and. * (data(nrec)(11:11).eq.'P'.or. * data(nrec)(11:14).eq.'AMPB'.or. * data(nrec)(11:13).eq.'AMB'.or. * data(nrec)(11:14).eq.'AMPb'.or. * data(nrec)(11:13).eq.'AMb'.or. * data(nrec)(11:14).eq.'AMP '.or. * data(nrec)(11:14).eq.'IAmb'.or. * data(nrec)(11:11).eq.' ' * )) then c c check other conditions c if(period.le.3.0.and.dist.le.100.0*111.2 * .and.period.ge.0.2.and.dist.ge.test(113)) then c c get attenuation term c if(test(108).ne.1.0.and.test(108).ne.0.0) then write(6,*)' Wrong value of test(108)' write(6,*)' will set to 0.0' test(108)=0.0 endif if(test(108).eq.0.0) call mbb_att(depth,dist/111.2,q) ! richter if(test(108).eq.1.0) call mb_att(depth,dist/111.2,q) ! veith c magnitude=alog10(amp/period)+q if(mag_cor(3,kstat).ne.99.0) then nmb=nmb+1 bmag(nmb)=magnitude+mag_cor(3,kstat) av_bmag=av_bmag+bmag(nmb) endif c c write distance and depth c if(nline.ge.0.and.output) then if(mag_cor(3,kstat).ne.99.0) then if(mag_cor(3,kstat).eq.0.0) then write(6,'(1x,a7,2x,a,2x, * f8.1,3x,a,2x,f12.1,a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'mb = ',magnitude else write(6,'(1x,a7,2x,a,2x, * f8.1,3x,a,2x,f12.1,3x,a,f5.1,a,f5.2)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'mb = ', * magnitude,' + ', * mag_cor(3,kstat) endif endif nline=nline+1 endif if(nline.gt.maxline.and.output) then nline=0 write * (6,'('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif else if(nline.ge.0.and.output) then if(dist.lt.test(113))then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mb available, event too close, ', * ' amplitude not used' endif if(dist.gt.100.0*111.2) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mb available, event too far', * ' amplitude not used' endif if(period.gt.3.0.or.period.lt.0.2) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mb available, wrong period,', * ' amplitude not used' endif nline=nline+1 if(nline.gt.maxline.and.output) then nline=0 write * (6,' * ('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') * nline=-1 endif endif endif endif c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c mB only up to 100 deg, priod range 0.2-30 s ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c check if mB c IF(AMP.GT.0.0 * .and.data(nrec)(11:17).eq.'IVmB_BB') then c c check other conditions c if(period.le.30.0 * .and.(period.ge.0.2.or.period.eq.0.0) * .and.dist.ge.test(113).and.dist * .le.100.0*111.2) then c c get attenuation term c if(test(108).ne.1.0.and.test(108).ne.0.0) then write(6,*)' Wrong value of test(108)' write(6,*)' will set to 0.0' test(108)=0.0 endif if(test(108).eq.0.0) call mbb_att(depth,dist/111.2,q) ! richter if(test(108).eq.1.0) call mb_att(depth,dist/111.2,q) ! veith c magnitude=alog10(amp/6.28)+q if(mag_cor(3,kstat).ne.99.0) then nmbb=nmbb+1 bbmag(nmbb)=magnitude+mag_cor(3,kstat) av_bbmag=av_bbmag+bbmag(nmbb) endif c c write distance and depth c if(nline.ge.0.and.output) then if(mag_cor(3,kstat).ne.99.0) then if(mag_cor(3,kstat).eq.0.0) then write(6,'(1x,a7,2x,a,2x, * f8.1,3x,a,2x,f12.1,a,f5.1,3x,a,f5.1)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'mB = ',magnitude else write(6,'(1x,a7,2x,a,2x, * f8.1,3x,a,2x,f12.1,a,f5.1,3x,a,f5.1,a,f5.2)') * data(nrec)(2:8),' dist:',dist, * 'amp:',amp,' T: ',period,'mB = ', * magnitude,' + ', * mag_cor(3,kstat) endif endif nline=nline+1 endif if(nline.gt.maxline.and.output) then nline=0 write * (6,'('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif else ! mB but some conditions not met if(nline.ge.0.and.output) then if(period.lt.0.2.or.period.gt.30.0) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mB available, wrong period, ', * 'amplitude not used' endif if(dist.gt.100.0*111.2) then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mB available, event too far', * ' amplitude not used' endif if(dist.lt.test(113))then write(6,'(1x,a7,a,a)')data(nrec)(2:8), * ' Amp for mB available, event too close, ', * ' amplitude not used' endif endif nline=nline+1 c if(nline.gt.maxline) then if(nline.gt.maxline.and.output) then nline=0 write * (6,'('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif endif endif 1 CONTINUE c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c write magnitudes cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c write ml in first mag position if mag available c if(nml.ne.0) then av_xmag=av_xmag/nml c c could be median magnitude c if(test(116).eq.1.0) then call mdian(xmag,nml,av_xmag) endif if(av_xmag.lt.-10.0.or.av_xmag.gt.10.0) xmag=-9.9 write(data(1)(56:63),'(f4.1,a1,a3)')av_xmag,'L',agency c c add residual into data5 c do k=nhead+1,nrecord if (data5(k)(62:63).eq.'ML') then read(data5(k)(64:68),'(f5.2)') magnitude if (nml.gt.0.and.av_xmag.gt.-9.9) then write(data5(k)(70:74),'(f5.2)') magnitude-av_xmag endif endif enddo endif c c c check if mc magnitude is available, if mag in first pos. write c in second postion, else first. c IF(nmc.NE.0) THEN av_cmag=av_CMAG/NMC c c could be median magnitude c if(test(116).eq.1.0) then call mdian(cmag,nmc,av_cmag) endif if(av_cmag.lt.-10.0.or.av_cmag.gt.10.0) av_cmag=-9.9 if(data(1)(60:60).ne.' ') then write(data(1)(64:71),202) av_cmag,agency else write(data(1)(56:63),202) av_cmag,agency 202 format(f4.1,'C',a3) endif endif c c check if Ms is available, write in first available position c if(nms.gt.0) then av_smag=av_smag/nms c c could be median magnitude c if(test(116).eq.1.0) then call mdian(smag,nms,av_smag) endif if(av_smag.lt.-10.0.or.av_smag.gt.10.0) smag=-9.9 if(data(1)(60:60).eq.' ')then write(data(1)(56:63),203) av_smag,agency elseif(data(1)(68:68).eq.' ')then write(data(1)(64:71),203) av_smag,agency elseif(data(1)(76:76).eq.' ')then write(data(1)(72:79),203) av_smag,agency else ! all positions used on first line c c put in a new line with magnitudes c do k=nrecord,2,-1 ! move records one down data(k+1)=data(k) enddo data(2)=' ' ! new line with more magnitudes data(2)(46:48)=data(1)(46:48) ! use same agency data(2)(1:23)=data(1)(1:23) ! use same start write(data(2)(56:63),203) av_smag,agency ! put in magnitude data(2)(80:80)='1' nhead=nhead+1 nrecord=nrecord+1 203 format(f4.1,'s',a3) endif endif c c check if mS is available, write in first available position c if(nmss.gt.0) then av_ssmag=av_ssmag/nmss c c could be median magnitude c if(test(116).eq.1.0) then call mdian(ssmag,nmss,av_ssmag) endif if(av_ssmag.lt.-10.0.or.av_ssmag.gt.10.0) av_ssmag=-9.9 if(data(1)(60:60).eq.' ')then write(data(1)(56:63),703) av_ssmag,agency elseif(data(1)(68:68).eq.' ')then write(data(1)(64:71),703) av_ssmag,agency elseif(data(1)(76:76).eq.' ')then write(data(1)(72:79),703) av_ssmag,agency else ! all positions used on first line c c put in a new line with magnitudes c do k=nrecord,2,-1 ! move records one down data(k+1)=data(k) enddo data(2)=' ' ! new line with more magnitudes data(2)(46:48)=data(1)(46:48) ! use same agency data(2)(1:23)=data(1)(1:23) ! use same start write(data(2)(56:63),703) av_ssmag,agency ! put in magnitude data(2)(80:80)='1' nhead=nhead+1 nrecord=nrecord+1 703 format(f4.1,'S',a3) endif endif c c check if mb is available, write in first available position c if(nmb.gt.0) then av_bmag=av_bmag/nmb c c could be median magnitude c if(test(116).eq.1.0) then call mdian(bmag,nmb,av_bmag) endif if(av_bmag.lt.-10.0.or.av_bmag.gt.10.0) av_bmag=-9.9 if(data(1)(60:60).eq.' ')then write(data(1)(56:63),205) av_bmag,agency elseif(data(1)(68:68).eq.' ')then write(data(1)(64:71),205) av_bmag,agency elseif(data(1)(76:76).eq.' ')then write(data(1)(72:79),205) av_bmag,agency 205 format(f4.1,'b',a3) else do i=2,nhead ! check if other headers can be used if(data(i)(80:80).eq.'1'. ! must have same loc. agency * and.data(1)(46:48).eq.data(i)(46:48)) then write(data(i)(64:71),205) av_bmag,agency ! must be 2. position mag_write=.true. endif enddo c c put in a new line with magnitudes c if(.not. mag_write) then do k=nrecord,2,-1 ! move records one down data(k+1)=data(k) enddo data(2)=' ' ! new line with more magnitudes data(2)(46:48)=data(1)(46:48) ! use same agency data(2)(22:22)=data(1)(22:22) ! use same type data(2)(1:10)=data(1)(1:10) ! use same y, m, d data(2)(80:80)='1' write(data(2)(56:63),205) av_bmag,agency ! put in magnitude nhead=nhead+1 nrecord=nrecord+1 endif endif endif c c check if mB is available, write in first available position c if(nmbb.gt.0) then av_bbmag=av_bbmag/nmbb c c could be median magnitude c if(test(116).eq.1.0) then call mdian(bbmag,nmbb,av_bbmag) endif if(av_bbmag.lt.-10.0.or.av_bbmag.gt.10.0) av_bbmag=-9.9 if(data(1)(60:60).eq.' ')then write(data(1)(56:63),905) av_bbmag,agency elseif(data(1)(68:68).eq.' ')then write(data(1)(64:71),905) av_bbmag,agency elseif(data(1)(76:76).eq.' ')then write(data(1)(72:79),905) av_bbmag,agency 905 format(f4.1,'B',a3) else do i=2,nhead ! check if other headers can be used if(data(i)(80:80).eq.'1'. ! must have same loc. agency * and.data(1)(46:48).eq.data(i)(46:48)) then write(data(i)(64:71),905) av_bbmag,agency ! must be 2. position mag_write=.true. endif enddo c c put in a new line with magnitudes c if(.not. mag_write) then do k=nrecord,2,-1 ! move records one down data(k+1)=data(k) enddo data(2)=' ' ! new line with more magnitudes data(2)(46:48)=data(1)(46:48) ! use same agency data(2)(22:22)=data(1)(22:22) ! use same type data(2)(1:10)=data(1)(1:10) ! use same y, m, d data(2)(80:80)='1' write(data(2)(56:63),905) av_bbmag,agency ! put in magnitude nhead=nhead+1 nrecord=nrecord+1 endif endif endif RETURN END c--------------------------------------------------------------------------- subroutine read_stat_mod &(agency,model_code,model,minstat,minphase, &modnam,loc_indicator,iustat,iuphs,iulst,iusum,iutab,isort, &test,dlt,yter,maxline,use_eev) c c read model and stations etc, assumed from unit 1, the file name c is STATIONX.HYP where X can be any character given by model. If c model is blank, the standard model 0 is assumed. The routine will c first look for the station model in the curent directory, then c in the DAT directory. If the model file is not found, model=@ is c returned c c iustat: 0: no station list, iulst: unit for print.out file c c j havskov jan 92 c c jul 92 by jh : add new par line for event type c jul 93 : iasp to capital letters c april 94 brl : added I/O numbers - was using #'s different from main program!! c apr 94 jh : do not read addtional parameters except agency CJAB(BGS)Jan95 : Install file & error handling... CJAB(BGS)Jan95 : & remove VAX.. c implicit none save include 'libsei.inc' ! Library definitons. include 'hypparm.inc' !added 6/98 BRL external sei get file, ! Find file & open handler. & sei close, ! & closure. & sei clen, ! String length. & sei code ! Error condition handler. integer sei clen ! & function. integer code ! Condition. logical b_flag ! A problem flagged?. logical use_eev ! true if call from eev c c character*80 data(*) character*3 agency character*20 model_code character*1 yter integer minstat,minphase,init,maxline character*1 model c--- location indicator character*1 loc_indicator(3) c--- name of station-model file in current directory character*80 cur_file c--- name of stationmodel in DAT directory character*80 modnam c--- event # integer eventno c--- isort & test added c integer isort(*) logical locate c-- covarrience matrix and origin time error real covar(3,3),oterr c--- logical I/O integer iustat,iuphs,iulst,iusum,iutab integer idummy,nhead,nrecord integer nmessage ! number of messages for se character*80 message(50) ! --------------- c real test(200),dlt(narriv) model_code=' ' ! currently do no tuse c c Open the stations file... c ========================= c cur_file = 'STATION0.HYP' ! Set up filename. if(model.ne.' ') cur_file(8:8)=model ! & adjust call sei get file( open$+ignore$, ! Find file & open it. & iustat, ! On unit. & code, ! Returned condition. & 'DAT', ! Alternative directory to search & cur_file ) ! For stations file. c if( code .ne. e_ok$ ) then ! Does not exist. chr_err_msg$ = cur_file(:seiclen(cur_file)) // & ' does not exist' c call sei code( stop$, ! Halt the program with user-. call sei code( warn$, ! Halt the program with user-. & e_misf$, ! Message; dummy use of code. & 0, ! Unit (n/a). & b_flag ) ! Flag (n/a). model='@' return end if ! c c check for iaspei files c ---------------------- c call iasp91_filename(modnam) modnam = modnam(:index(modnam,' ')-1) // '.TBL' ! File name call sei get file( check$, ! Find file. & 0, ! On unit (n/a). & code, ! Returned condition. & 'DAT', ! Alternative directory to search & modnam ) ! For iaspei file. c modnam = 'IASP91.TBL' ! File name. c c File exists somewhere... c if( code .eq. e_ok$ ) then ! Exists in pathname. modnam = modnam(:index(modnam,'.TBL')-1) ! Extract root. c c File nowhere to be found... c else ! Does not exist. chr_err_msg$ = '**** WARN: ' // & modnam(:seiclen(modnam)) // & ' does not exist' write(*,*) call sei code( warn$, ! Warn the user it doesn't exist. & e_misf$, ! Message; dummy use of code. & 0, ! Unit (n/a). & b_flag ) ! Flag (n/a). write(*,*) modnam= ' ' end if ! C C READ STATION LIST AND MODEL C call hypocent(1,iustat,iuphs,iulst,iusum,iutab,init,'N', &data,modnam,eventno,dlt,isort,test,yter,maxline,idummy, &locate,covar,oterr,nhead,nrecord,use_eev,nmessage,message) c c read agency c read(iustat,'(a3)',iostat=code) agency ! Get agency. call sei code(fort$,code,iustat,b_flag) ! Process outcome. call sei close(close$,iustat,code) ! Close the file (stop on error). c c get minstat etc c minstat=test(79) minphase=test(80) loc_indicator(1)='L' loc_indicator(2)='R' loc_indicator(3)='D' if(test(81).ne.1.0) loc_indicator(1)=' ' if(test(82).ne.1.0) loc_indicator(2)=' ' if(test(83).ne.1.0) loc_indicator(3)=' ' c c make sure minphase has been set to at least 3 c changed to 2 6/98 BRL if(minphase.lt.2) minphase=2 return end c c------------------------------------------------------------------ c subroutine clear_old(data,nhead,nrecord,agency) c c clear old secondary magnitude lines, old error lines c and old id lines c implicit none save integer nrecord,nhead character*80 data(*) integer i,j,k character*3 agency c c delete secondary header lines with magnitude, only delete line c if location agency is the same as default agency c k=1 if(nhead.gt.1) then do i=2,nhead if(.not.(data(i)(80:80).eq.'1'. * and.agency.eq.data(i)(46:48)).and. * data(i)(80:80).ne.'E'.and. * data(i)(79:80).ne.'83'.and. * data(i)(79:80).ne.'93'.and. * data(i)(79:80).ne.'A3'.and. * data(i)(80:80).ne.'H' * .and.(.not. * (data(i)(80:80).eq.'3'.and.data(i)(2:7).eq.'ACTION')) * .and.(.not. * (data(i)(80:80).eq.'3'.and.data(i)(2:13) * .eq.'SPEC SD ')) * .and.(.not. * (data(i)(80:80).eq.'3'.and.data(i)(2:13) * .eq.'SPEC AVERAGE'))) then k=k+1 if(k.ne.i) data(k)=data(i) endif enddo j=k ! save new number of headers if(k.lt.nhead) then ! some header lines removed, shift rest do i=nhead+1,nrecord k=k+1 data(k)=data(i) enddo nhead=j do i=k,nrecord ! blank double lines data(i)=' ' enddo nrecord=k endif endif c c return end c-------------------------------------------------------------------- subroutine update_spec(data1,nrecord,nhead,agency,test,output) c c recalcualate spectral values so if distance has been changed, c the values are updated. Then calculate average values and write out c in nordic data array data and return new number of headers and records. c the main header is updated with average mw if there is room and agency c is written in the usual way. c implicit none save include 'hypparm.inc' c include 'seidim.inc' ! since now in hypparm.inc jan 2013 c c the following 2 otherwise defined in mulplt.inc c character*(*) data1(*) real edepth ! depth character*1 spec_phase ! nearest phase to spectrum character*3 agency ! magnitude agency character*5 station, stat5 character*4 compon ! component integer nrecord,nhead ! number of record and headers real depth ! hypocentral depth c--- spectral values, see mulplt.inc real mw(max_data),moment(max_data),sdrop(max_data), * sslope,q_below1hz, * radius(max_data),swindow, * omega0(max_data),cornerf(max_data), * spec_velocity,density, * q0,qalpha real mw_calc(1000) ! calcualted mw integer secc,travel_time,hour,min real geo_distance real kappa ! near term attenuation real edistance,xdist ! epicentral distance real avmw,avmoment,avsdrop,avsslope,avradius,avswindow,avomega0 real avcornerf real sdmw,sdmoment,sdsdrop,sdradius,sdomega0,sdcornerf logical mag_write ! true if room for magnitude on an existing line logical output ! true if screen output integer nline ! line count for magnitude output integer maxline ! maximum magnitude lines out c integer nspec,nspec1 ! number of spectra integer nmw ! number of mw, can be less than nspec if w. out real pi real surface,radpat ! surface effect and radiation pattern effect real factor character*80 text integer i,k character*5 stat(nstats) ! station list real mag_cor(5,nstats) ! magnitude residuals integer kstat ! station number currently used integer seiclen common/hyp/stat,mag_cor,nline,maxline c pi=3.14159265 c output=.true. surface=2.0 radpat=0.6 c c read depth c read(data1(1)(39:43),'(f5.1)') edepth c nspec=0 nspec1=0 nmw=0 avmw=0.0 avmoment=0.0 avsdrop=0.0 avsslope=0.0 avradius=0.0 avswindow=0.0 avcornerf=0.0 c c look for spectral values c k=1 c c back here in loop when reading spectral values c 10 continue k=k+1 if(k.gt.nhead) goto 20 ! calculate average values when all h. read if(data1(k)(80:80).eq.'3'.and.data1(k)(2:5).eq.'SPEC') then c if(data1(k)(7:14).eq.data1(k+1)(7:14).and.data1(k)(16:30). if(data1(k)(6:14).eq.data1(k+1)(6:14).and.data1(k)(16:30). * ne.data1(k+1)(16:30)) then ! check if pair of lines ok nspec1=nspec1+1 c c find current distance from station c do i=nhead+1,nrecord-1 c c get 5 letter name , jh jan 2014 c stat5=data1(k)(6:10) if(stat5(1:1).eq.' ') then stat5(1:4)=stat5(2:5) stat5(5:5)=' ' endif c if(data1(k)(7:10).eq.data1(i)(2:5)) then if(stat5.eq.data1(i)(2:6)) then if(data1(i)(71:75).eq.' ') then write(6,'(a,a)') * 'No distance for station, spec not used ',stat5 goto 2 ! no distance endif read(data1(i)(71:75),'(f5.0)') edistance goto 1 endif 2 continue enddo c c if here, no distance found, jump to next spectral values c k=k+1 goto 10 1 continue c c read old spectral values c read(data1(k),300) c * station(1:4),compon,moment(nspec1),sdrop(nspec1), * station(1:5),compon,moment(nspec1),sdrop(nspec1), * omega0(nspec1),cornerf(nspec1), * radius(nspec1),sslope,swindow,mw(nspec1) c300 format(6x,a4,a4,3x,f5.1,3x,f5.1,3x,f5.1, 300 format(5x,a5,a4,3x,f5.1,3x,f5.1,3x,f5.1, * 3x,f5.2,2x,f6.2,3x,f5.2,3x,f5.1,3x,f5.1) c station(5:5)=' ' ! only possible to use 4 char station in spec if (station(1:1).eq.' ') station=station(2:) ! allow for 5 characters read(data1(k+1),301) c * station(1:4),compon,hour,min,secc,kappa,xdist, * station(1:5),compon,hour,min,secc,kappa,xdist, * spec_velocity,density,q0,qalpha,q_below1hz spec_phase=data1(k+1)(41:41) ! type of spectrum if(data1(k+1)(24:25).eq.'TR') kappa=0.0 ! from old seisan c301 format(6x,a4,a4,2x,3i2,3x,f5.3,3x, 301 format(5x,a5,a4,2x,3i2,3x,f5.3,3x, * f5.0,3x,f5.2, * 3x,f5.2,3x,f5.1,3x,f5.2,3x,f5.2) c station(5:5)=' ' if (station(1:1).eq.' ') station=station(2:) c c correct if old data has VS, then q was frequency dependent c in whole frequency range c if(data1(k+1)(72:73).eq.'VS') q_below1hz=0.0 c c c find station number c do i=1,nstats if(stat(i).eq.station) then kstat=i goto 1245 endif enddo c c if here, station not found c write(6,'(a,a5)') * ' WARNING ************ No spectral update for station ', * station write(6,'(a)')' Station not in station list' goto 10 1245 continue c c recalculate c c c calculate geo_distance c call spec_dist * (spec_phase,edistance,edepth,factor) c geo_distance=1.0/factor factor= 4*pi*(density*1000.0)* * ((spec_velocity*1000.0)**3) moment(nspec1) = factor*(geo_distance*1000.0)* * (10.0**omega0(nspec1))/(radpat*surface*1.0e9) c radius(nspec1)=0.35*svelocity/cornerf(nspec1) radius(nspec1)=0.37*spec_velocity/cornerf(nspec1) ! lo 14/10/2011 sdrop(nspec1) = (0.44*moment(nspec1))/ * (1.0e14*radius(nspec1)**3) moment(nspec1) = alog10(moment(nspec1)) if(moment(nspec1).gt.0.0) then if(mag_cor(5,kstat).ne.99.0) then mw(nspec1)=moment(nspec1)*0.667-6.06 else if(output) * write(6,'(1x,a,a)') station, ' mw weighted out' endif endif c c write geo_distance, amplitude and magnitude c if(nline.ge.0.and.output) then if(mag_cor(5,kstat).ne.99.0) then if(mag_cor(5,kstat).eq.0.0) then write * (6,'(1x,a4,1x,2a1,2x,a,2x,f8.1, * 3x,a,2x,f12.1,3x,a,f5.1)') * data1(k)(7:10),data1(k)(11:11), * data1(k)(14:14),' gdist:',geo_distance, * 'mom:',moment(nspec1),'mw = ',mw(nspec1) else ! magnitude correction write * (6,'(1x,a7,2x,a,2x,f8.1,3x,a,2x,f12.1, * 3x,a,f5.1,a,f5.2)') * data1(k)(7:13),' gdist:',geo_distance, * 'mom:',moment(nspec1),'mw = ',mw(nspec1),' + ', * mag_cor(5,kstat) mw(nspec1)=mw(nspec1)+mag_cor(5,kstat) endif else if(output)write(6,'(1x,a7,a,a)') data1(k)(5:13), * ' Mom for mw available, magnitude not used', * ' (weighted out)' endif nline=nline+1 endif if(nline.gt.maxline.and.output) then nline=0 write(6,'('' Return to continue, q to end listing '',$)') read(5,'(a)') text if(text(1:1).eq.'q'.or.text(1:1).eq.'Q') nline=-1 endif c c check for overflow c if(moment(nspec1).gt.999.0) moment(nspec1)=999.9 if(sdrop(nspec1).gt.999.0) sdrop(nspec1)=999.9 if(geo_distance.gt.99999) geo_distance=99999 c c rewrite values c c jh 5-2001, strange seems recalculated spec not output if (seiclen(station).lt.5) then write(data1(k),400) * station(1:4),compon,moment(nspec1),sdrop(nspec1), * omega0(nspec1), * radius(nspec1),sslope,swindow,mw(nspec1) else write(data1(k),405) * station(1:5),compon,moment(nspec1),sdrop(nspec1), * omega0(nspec1), * radius(nspec1),sslope,swindow,mw(nspec1) endif 400 format(' SPEC ',a4,a4,' MO',f5.1,' ST',f5.1,' OM',f5.1, * ' f0',5x,' R',f6.2,' AL',f5.2,' WI',f5.1,' MW',f5.1) 405 format(' SPEC',a5,a4,' MO',f5.1,' ST',f5.1,' OM',f5.1, * ' f0',5x,' R',f6.2,' AL',f5.2,' WI',f5.1,' MW',f5.1) c c check size of cornerf c if(cornerf(nspec1).lt.1.0) write(data1(k)(42:46),'(f5.3)') * cornerf(nspec1) if(cornerf(nspec1).lt.10.0.and.cornerf(nspec1).ge.1.0) * write(data1(k)(42:46),'(f5.2)') cornerf(nspec1) if(cornerf(nspec1).ge.10.0.and.cornerf(nspec1).lt.100.0) * write(data1(k)(42:46),'(f5.1)') cornerf(nspec1) if(cornerf(nspec1).ge.100.0) write(data1(k)(42:46),'(f5.0)') * cornerf(nspec1) if(sdrop(nspec1).lt.1.0) * write(data1(k)(26:30),'(f5.3)') sdrop(nspec1) if(radius(nspec1).lt.10.0) * write(data1(k)(49:54),'(f6.4)') radius(nspec1) data1(k)(80:80)='3' if (seiclen(station).lt.5) then write(data1(k+1),401) * station(1:4),compon,hour,min,secc, * kappa,int(geo_distance), * spec_velocity,density,q0,qalpha,q_below1hz else write(data1(k+1),406) * station(1:5),compon,hour,min,secc, * kappa,int(geo_distance), * spec_velocity,density,q0,qalpha,q_below1hz endif 401 format(' SPEC ',a4,a4,' T',3i2,' K ',f5.3,' GD', * I5,' V ',f5.2, * ' DE',f5.2,' Q0',f5.1,' QA',f5.2,' Q1',f5.2) 406 format(' SPEC',a5,a4,' T',3i2,' K ',f5.3,' GD', * I5,' V ',f5.2, * ' DE',f5.2,' Q0',f5.1,' QA',f5.2,' Q1',f5.2) data1(k+1)(80:80)='3' data1(k+1)(41:41)=spec_phase if(geo_distance.lt.100.0) write(data1(k+1)(34:38),'(f5.1)') * geo_distance c c sum for average, do not use magnitude weighted out value c do not use unidentified spectra, must be p or s (new may 2001, jh) c if(spec_phase.eq.'P'.or.spec_phase.eq.'S') then avmoment=avmoment+moment(nspec1) avsdrop=avsdrop+sdrop(nspec1) avomega0=avomega0+omega0(nspec1) avcornerf=avcornerf+cornerf(nspec1) avradius=avradius+radius(nspec1) avsslope=avsslope+sslope avswindow=avswindow+swindow if(mag_cor(5,kstat).ne.99.0) then nmw=nmw+1 avmw=avmw+mw(nspec1) mw_calc(nmw)=mw(nspec1) ! save calculated mw else write(6,'(4a)')' Moment magnitude not used in average ', * 'since weighted out ',station(1:4),compon endif nspec=nspec+1 else write(6,'(3a)') * ' Unidentified spectrum, not used in average ', * station(1:4),compon endif endif c k=k+1 ! do not use next spectral line endif goto 10 ! get next spectral line c c enter here when all headers checked c 20 continue c c return if no spectral values c if(nspec.gt.0) write(6,'(a,2i6)') *' Number of spectra available and number used in average', *nspec1,nspec if(nspec.eq.0) return c c calculate average values c avmoment=avmoment/nspec avsdrop=avsdrop/nspec avomega0=avomega0/nspec avcornerf=avcornerf/nspec avradius=avradius/nspec avsslope=avsslope/nspec avswindow=avswindow/nspec if(nmw.gt.0) then c avmw=avmw/nspec ! lot 22/07/2008 avmw=avmw/nmw c c could be median magnitude c if(test(116).eq.1.0) then call mdian(mw_calc,nmw,avmw) endif else avmw=9.9 ! indicate no mw endif c c check for overflow c c c compute standard deviation c sdmoment=0. sdsdrop=0. sdomega0=0. sdradius=0. sdcornerf=0. sdmw=0. if(avmoment.gt.999.0) then avmoment=999.9 else call calcsd(moment,avmoment,sdmoment,nspec) endif if(avmw.gt.9.9.or.avmw.lt.-9) then avmw=9.9 else call calcsd(mw,avmw,sdmw,nmw) endif if(avsdrop.gt.999.0) then avsdrop=999.9 else call calcsd(sdrop,avsdrop,sdsdrop,nspec) endif if(avomega0.gt.999.0) then avomega0=999.9 else call calcsd(omega0,avomega0,sdomega0,nspec) endif if(avradius.gt.999.9) then avradius=999.99 else call calcsd(radius,avradius,sdradius,nspec) endif if(avcornerf.gt.999.9) then avcornerf=999.99 else call calcsd(cornerf,avcornerf,sdcornerf,nspec) endif if(abs(avsslope).gt.999.0) avsslope=999.9 if(avswindow.gt.999.9) avswindow=999.9 c c move all lines down one to make room for new sd spectral line c if not there from before c k=0 do i=2,nhead if(data1(i)(1:13).eq.' SPEC SD ') k=i enddo if(k.eq.0) then do i=nrecord,2,-1 data1(i+1)=data1(i) enddo k=2 nrecord=nrecord+1 nhead=nhead+1 endif c c check size c if(sdcornerf.gt.999.0) sdcornerf=999.0 if(sdsdrop.gt.999.0) sdsdrop=999.0 if(sdradius.gt.999.0) sdradius=999.0 c c write values c write(data1(k),500) *sdmoment,sdsdrop,sdomega0,sdradius,sdmw c if(sdcornerf.lt.1.0) write(data1(k)(42:46),'(f5.3)') * sdcornerf if(sdcornerf.lt.10.0.and.sdcornerf.ge.1.0) * write(data1(k)(42:46),'(f5.2)') sdcornerf if(sdcornerf.ge.10.0.and.sdcornerf.lt.100.0) * write(data1(k)(42:46),'(f5.1)') sdcornerf if(sdcornerf.ge.100.0) write(data1(k)(42:46),'(f5.0)') * sdcornerf if(sdsdrop.lt.1.0) *write(data1(k)(26:30),'(f5.3)') sdsdrop if(sdradius.lt.10.0) *write(data1(k)(49:54),'(f6.4)') sdradius 500 format(' SPEC ','SD ',' MO',f5.1,' ST',f5.1,' OM',f5.1, *' f0',5x,' R',f6.2,' AL ',' WI ',' MW',f5.1) data1(k)(80:80)='3' c c move all lines down one to make room for new average spectral line c if not there from before c k=0 do i=2,nhead if(data1(i)(1:13).eq.' SPEC AVERAGE') k=i enddo if(k.eq.0) then do i=nrecord,2,-1 data1(i+1)=data1(i) enddo k=2 c nrecord=nrecord+1 nhead=nhead+1 endif c c write values c write(data1(k),501) *avmoment,avsdrop,avomega0, *avradius,avsslope,avswindow,avmw c c check size c if(avcornerf.lt.1.0) write(data1(k)(42:46),'(f5.3)') avcornerf if(avcornerf.lt.10.0.and.avcornerf.ge.1.0) * write(data1(k)(42:46),'(f5.2)') avcornerf if(avcornerf.ge.10.0.and.avcornerf.lt.100.0) * write(data1(k)(42:46),'(f5.1)') avcornerf if(avcornerf.ge.100.0) write(data1(k)(42:46),'(f5.0)') * avcornerf if(avsdrop.lt.1.0) *write(data1(k)(26:30),'(f5.3)') avsdrop if(avradius.lt.10.0) *write(data1(k)(49:54),'(f6.4)') avradius 501 format(' SPEC ','AVERAGE ',' MO',f5.1,' ST',f5.1,' OM',f5.1, *' f0',5x,' R',f6.2,' AL',f5.2,' WI',f5.1,' MW',f5.1) data1(k)(80:80)='3' if(avmw.eq.9.9) data1(k)(75:78)=' ' ! no mw c c update header, if a mw,, find which header and where there is a free spot c if(avmw.ne.9.9) then if(data1(1)(60:60).eq.' ')then write(data1(1)(56:63),205) avmw,agency elseif(data1(1)(68:68).eq.' ')then write(data1(1)(64:71),205) avmw,agency elseif(data1(1)(76:76).eq.' ')then write(data1(1)(72:79),205) avmw,agency 205 format(f4.1,'W',a3) else mag_write=.false. do i=2,nhead ! check if other headers can be used if(data1(i)(80:80).eq.'1'. ! must have same loc. agency * and.data1(1)(46:48).eq.data1(i)(46:48)) then mag_write=.true. if(data1(i)(60:60).eq.' ')then write(data1(i)(56:63),205) avmw,agency elseif(data1(i)(68:68).eq.' ')then write(data1(i)(64:71),205) avmw,agency elseif(data1(i)(76:76).eq.' ')then write(data1(i)(72:79),205) avmw,agency endif endif enddo c c put in a new line with magnitudes c if(.not. mag_write) then do k=nrecord,2,-1 ! move records one down data1(k+1)=data1(k) enddo data1(2)=' ' ! new line with more magnitudes data1(2)(46:48)=data1(1)(46:48) ! use same agency data1(2)(1:23)=data1(1)(1:23) ! use same start data1(2)(80:80)='1' write(data1(2)(56:63),205) avmw,agency ! put in magnitude nhead=nhead+1 nrecord=nrecord+1 endif endif endif c return end subroutine calcsd(data,avdata,sddata,n) implicit none real data(*) real avdata,sddata integer i,n do i=1,n sddata=sddata+(data(i)-avdata)**2 enddo sddata=sqrt(sddata/n) return end