C C HLRWM1.FOR program modified 2001/08/28 for CIT catalog C C C CALCULATION OF THE LIKELYHOOD FUNCTION, NOAA, DUDA CATALOG, C 30 NOVEMBER 1983 (USGS), MODIFICATIONS MAY 17, 1989 C MODIFIED VERSION -- ONE 'C' COEFFICIENT, C TIME INFLUENCE -- HYPERBOLIC FUNCTION C DISTANCE INFLUENCE -- RAYLEIGH DISTRIBUTION C VERTICAL INFLUENCE VERSION -- GAUSSIAN DISTRIBUTION C C calculation of weights for each earthquake, May 1989 C change for SMT catalogs (HARVARD) or CATL c c put this command so that only aftershocks are marked c [if (mi.lt.m(j)) go to 100] c C IMPLICIT REAL*8 (A-H, O-Z) C PARAMETER (NT = 15, NPAR = 7) PARAMETER (NTP1 = NT + 1) !NTP1 = NT + 1 C PARAMETER (NDATA = 10000) !NOAA c PARAMETER (NDATA = 15000) !HARVARD c PARAMETER (NDATA = 1500) !DUDA PARAMETER (NDATA = 1500) !CIT PARAMETER (NDEP = 120, NSPM = 64, NMG = 90, NMAG = 50) C PARAMETER (RMIN = 1.6D00, RMAX = 409.6D0*4.0D0) PARAMETER (RM4 = RMIN**4, RL4 = RMAX**4) PARAMETER (SCT = 3.0D0, MST = 40) PARAMETER (SIGMAS = 6.00D0, ERRORS = 1.0D01) C c REAL*8 EQD(4) !SMT HARVARD, SIPKIN c REAL*4 EQE(8) !SMT HARVARD, SIPKIN c INTEGER*2 EQH(16) !SMT HARVARD, SIPKIN C REAL*8 quat1(4) /0.0, 0.0, 0.0, 1.0/ REAL*8 EQD(3) !NOAA, DUDA REAL*4 EQE(6) !NOAA, DUDA INTEGER*2 EQH(12) !NOAA, DUDA INTEGER*2 EQH1(6) /6*0/ C DIMENSION IDD(NDEP), RD(NSPM), ISPM(NSPM), SPM(NSPM), RL(NMAG), 1 M(NDATA), SIGMH(NMAG), SIGMV(NMAG), TLT(NMAG), TDTHT(NMAG), 2 TLTHT(NMAG), TDB(NMAG), TDIFF(NMAG), TPROD(NMAG), 3 TLOGD(NMAG), TLOGL(NMAG), PIND8(NDATA), AMOM(NDATA), AMOM1(NDATA) REAL*8 T(NDATA), AFT(NDATA), AMR(NMAG), CZ(NPAR), 1 DZ(NPAR), SZ(NPAR), OCZ(NPAR), OZ(NPAR) REAL*4 AS(NDATA), BS(NDATA), CS(NDATA), XT, YT, RAD, 1 X(NDATA), Y(NDATA), DEPTH(NDATA), err(NDATA), PIND4(NDATA), SEC INTEGER*4 NCH(NPAR), MD(NMG) INTEGER*2 HX(4), HXE CHARACTER*8 RTIME CHARACTER*1 QUAL(5) /'A', 'B', 'C', 'D', 'E'/ LOGICAL FIRST, ONE, NAF C EQUIVALENCE (DTEMP, HX(1)) EQUIVALENCE (EQD(1), EQE(1), EQH(1)) EQUIVALENCE (CZ(1),AL), (CZ(2),CMU), (CZ(3),THETA), 1 (CZ(4),SIGMA), (CZ(5),ERRORH), (CZ(6),ERRORV), (CZ(7),AMU) EQUIVALENCE (DZ(1),GAL), (DZ(2),GC), (DZ(3),DTHETA), 1 (DZ(4),DSIGMA), (DZ(5),GERRH), (DZ(6),GERRV), (DZ(7),G) EQUIVALENCE (SZ(1),DALPHA), (SZ(2),DCMU), (SZ(3),STPTHE), 1 (SZ(4),STPSIG), (SZ(5),STPERH), (SZ(6),STPERV), (SZ(7),DAMU) C COMMON /POL/ PL,GA,SAL,SGA,DAL,DGA,DR1,DR2,MM(NMG),INT,N1,N2,N3 COMMON /DATE/ SEC, IYE, IMO, IDA, IH, IMI C NAMELIST /INPUT/ AMU, CMU, AL, THETA, SIGMA, ERRORH, ERRORV, 1 DAMU, DCMU, DALPHA, STPTHE, STPSIG, STPERH, STPERV, CONST, 2 EIT, TSTART, TEND, GCNGD, GCNGU, TMST, NITL, MGL, NITER C AMCOF = 1.5 NWR1 = NDATA/20 NWR2 = 5 STPERV = 5.0D-06 C ALATTUP= 37.0 !So. California ALATTLO= 32.0 !So. California ALONGUP=-114.0 !So. California ALONGLO=-122.0 !So. California C TSTRT = 3.462150D-03 c TSTRT = 3.462150D-05 c PI = 2.0D0*DACOS(0.0D0) HPI = PI/2.0D0 EAR = 1.0D04/HPI RAD = 90.0D0/HPI ERRC = 1.0D0/DSQRT(2.0D0*PI) AMST = DFLOAT(MST)/10.0D0 FIRST = .TRUE. TBASE = CTIME(1) DO 4000 I = 1, NPAR OZ(I) = 0.0D0 NCH(I) = 0 4000 CONTINUE C OPEN (UNIT = 5, STATUS = 'OLD', DISP = 'KEEP', READONLY, c 1 FILE = 'disk$user:[kagan]cit7.DAT') !CIT c 1 FILE = 'disk$user:[kagan]cit7one.DAT') !CIT 1 FILE = 'disk$user:[kagan]cit7one.tmp') !change tstart and MGL READ (5, INPUT) CLOSE (UNIT = 5) C ONE = NITL.EQ.1 TMST = TMST*0.95D0*3600.0D0 - 60.0D0*60.0D0 WRITE (6, INPUT) C OPEN (UNIT=8,STATUS='OLD',DISP='KEEP',FORM='FORMATTED',READONLY, 1 FILE = 'disk$user:[kagan.DATA]CIT32_00.DAT') ! CIT C READ (8, 155) IDU, IDL, NDEPTH READ (8, 150) IDD READ (8, 150) N READ (8, 150) (ISPM(I), I = 1, NSPM) 155 FORMAT (40X, 5I8) 150 FORMAT (8I10) CLOSE (UNIT = 8) C AN = N AMOMC = 10.0D0**(AMCOF*MGL/10.0D0) !CATL C AMCUT0 = AMCOF*MGL !SMT HARVARD C AMOMC = 10.0D0**(AMCUT0/10.0D0) !SMT HARVARD WRITE (6, 7) MGL, SCT, TSTRT, RMIN, AMOMC 7 FORMAT ('0', I5, F10.3, 3D20.6) C IF(ONE) WRITE (6, 215) IDD 215 FORMAT (' ', 20I6) IF(ONE) WRITE (6, 210) ISPM 210 FORMAT (' ', 10I10) C DO 1100 I = 1, NSPM 1100 SPM(I) = ISPM(I) IF(ONE) WRITE (6, 11) SPM 11 FORMAT (' ', 10F12.2) IF(ONE) WRITE (6, 10) DO 15 I = 1, NSPM 15 SPM (I) = SPM(I)*2.0D0/(AN*(AN - 1.0)) IF(ONE) WRITE (6, 1110) SPM 1110 FORMAT (' ', 10D12.5) IF(ONE) WRITE (6, 1110) IF(ONE) WRITE (6, 112) SPM 112 FORMAT (' ', 10F12.8) TE = DSQRT(10.D0) AT = SCT**.1D0 AR = 2.D0**0.25D0 TSTRT1 = TSTRT*(AT**(MGL - MST)) WRITE (6, 12) AT, AR 12 FORMAT ('0 AT, AR =', 2D15.5, /) C DO 20 I = 1, NMAG IT = I + MGL - 1 TDB(I) = TSTRT1*(AT**(I-1)) TLT(I) = TSTRT1*(TE**NT*AT**(I-1)) AMR(I) = 3.0D0*DSQRT(ERRORS**2 + 1 10.D0**(DFLOAT(I + MGL - 1 - MST)*0.1D0)*SIGMAS**2) RL(I) = DCOS(AMR(I)/EAR) 20 CONTINUE C DO 3090 I = 1, NSPM 3090 RD(I) = RMIN*(AR**(I - 1)) IF(ONE) WRITE (6, 10) IF(ONE) WRITE (6, 10) RD IF(ONE) WRITE (6, 10) 10 FORMAT (' ', 10F12.5) DO 3092 I = 1, NSPM-1 SPM(I) = SPM(I)/(RD(I+1) - RD(I)) 3092 CONTINUE IF(ONE) WRITE (6, 1110) SPM IF(ONE) WRITE (6, 10) IF(ONE) WRITE (6, 112) SPM IF(ONE) WRITE (6, 10) RD(1) = (2.0D0*SPM(1) + SPM(2))/3.0D0 RD(NSPM) = (2.0D0*SPM(NSPM) + SPM(NSPM - 1))/3.0D0 DO 192 I = 2, NSPM-1 RD(I) = (SPM(I-1) + 2.0D0*SPM(I) + SPM(I+1))/4.0D0 192 CONTINUE IF(ONE) WRITE (6, 1110) RD IF(ONE) WRITE (6, 10) IF(ONE) WRITE (6, 112) RD IF(ONE) WRITE (6, 10) C DO 236 J = 1, NMG 236 MD(J) = 0 TIME_INT = TEND - TSTART OLDF = - 1.0D 06 C OPEN (UNIT = 25, STATUS='OLD', DISP='KEEP', 1 READONLY, FORM = 'FORMATTED', c 2 FILE= 'DISK$USER:[KAGAN.SAVE.LIM3]cit32_00_47.DAT') !CIT c 2 FILE= 'DISK$USER:[KAGAN.SAVE.LIM3]cit32_00_47A.DAT') !CIT 2 FILE= 'DISK$USER:[KAGAN.SAVE.LIM3]cit32_00_47B.DAT') !CIT C READ (25, 3625) OCZ, CZ, OZ, SZ, CONST, OLDF READ (25, 3650) NCH, NITER 3625 FORMAT (4D20.12) 3650 FORMAT (8I10) FIRST = .FALSE. CLOSE (UNIT = 25) WRITE (6, 3777) OLDF 3777 FORMAT ('0 OLD VALUE OF LIKE = ', F15.7, /) WRITE (6, INPUT) 777 CONTINUE I = 0 N = 0 UU = 0.D0 MMAX = 0 WRITE (6, 1440) C OPEN (UNIT = 20, STATUS='OLD', DISP='KEEP', READONLY, 1 FORM='UNFORMATTED', RECORDTYPE='FIXED', 2 RECL=6,FILE='DISK$user:[kagan]cit32_00.tmp') !Run citre4.for c 2 RECL=6,FILE='DISK$DATA:[CATALOG]cit32_00.BIN') !CIT C SUMLOG1 = 0.0D0 SUMMOM1 = 0.0D0 SUMAL1 = 0.0D0 C 30 READ (20, END=40) EQD IF(EQD(1).GT.TEND) GO TO 40 N = N + 1 IF(EQD(1).LT.TSTART) GO TO 30 IF(EQH(9).LT.IDU.OR.EQH(9).GT.IDL) GO TO 30 c IEQH = DFLOAT(EQH(10)-900)*2.0D0/3.0D0 + 0.5D0 !SMT HARVARD c IEQH = (IEQH + 5)/10 !MNOAA,GREEK IEQH = eqh(10) !CIT MT = IEQH - MGL + 1 IF(MT.LE.0) GO TO 30 I = I + 1 AMOM(I) = 10.0D0**(AMCOF*(EQH(10)+0.5)/10.0D0) !CATL c AMOM(I) = 10.0D0**((EQH(10)-899.5)/100.0D0) !SMT HARVARD AMOM1(I) = amom(i) SUMAL1 = SUMAL1 + 1.0D0 SUMMOM1 = SUMMOM1 + AMOM(I)/AMOMC SUMLOG1 = SUMLOG1 + DLOG(AMOM(I)/AMOMC) UU = UU + DFLOAT(MT) C T(I) = EQD(1) X(I) = EQE(3) Y(I) = EQE(4) DEPTH(I) = EQH(9)/100.0 err(i) = EQH(11)/100.0 XT = (90.0 - EQE(3))/RAD YT = EQE(4)/RAD AS(I) = SIN(XT)*COS(YT) BS(I) = SIN(XT)*SIN(YT) CS(I) = COS(XT) M(I) = MT IF (MMAX.LT.MT) MMAX = MT MD(MT) = MD(MT) + 1 NN = MOD(I, NWR1) C IF(NN.EQ.1.AND.ONE) IF ((I.GT.100.AND.NN.EQ.1).OR.(I.LE.100)) 1 WRITE (6,1440) N,I,EQD(1),(EQE(J),J=3,4),(EQH(J),J=9,12),AMOM(I) 1440 FORMAT (' ', 2I6, F20.8, 2F12.5, 4I8, D15.6) GO TO 30 40 CONTINUE CLOSE (UNIT = 20) WRITE (6,1440) WRITE (6,1440) N,I,EQD(1),(EQE(J),J=3,4),(EQH(J),J=9,12),AMOM(I) WRITE (6,1440) BETA1 = SUMAL1/SUMLOG1 BVAL1 = BETA1*AMCOF WRITE (6, 4450) SUMAL1, SUMMOM1, SUMLOG1, BETA1, BVAL1 4450 FORMAT ('0 SUMAL1, SUMMOM1, SUMLOG1, BETA1, BVAL1 = ', 5G15.7, /) NM = I NWR = NM**2/50 WRITE (6, 1440) WRITE (6, 45) NM, N, NWR, NWR1, NWR2, MMAX WRITE (6, 45) MD 45 FORMAT (' ',10I10) C CALL POLI (CONST,MD,1,ICODE) IF (FIRST) AL = PL ANU = 10.0D0**(0.1D0*GA) GPNU = DR2*GA*DLOG(10.0D0)/10.0D0 WRITE (6,779) CONST,AL,DR1,ANU,GPNU,GA,DR2,N1,INT,N3 779 FORMAT ('0', 7F12.4, 3I5, /) IF(ICODE.GT.200) STOP 12 C DO 1000 NIT = 1, NITL C IF(ONE) WRITE (6, 779) DO 5010 I = 1, NMAG IT = I + MGL - 1 SIGMH(I) = DSQRT(ERRORH**2 + 1 10.D0**(DFLOAT(I + MGL - 1 - MST)*0.1D0)*SIGMA**2) SIGMV(I) = DSQRT(ERRORV**2 + 1 10.D0**(DFLOAT(I + MGL - 1 - MST)*0.1D0)*SIGMA**2) TDTHT(I) = TDB(I)**THETA TLTHT(I) = TLT(I)**THETA TPROD(I) = TDTHT(I)*TLTHT(I) TDIFF(I) = TLTHT(I) - TDTHT(I) TLOGD(I) = DLOG(TDB(I)) TLOGL(I) = DLOG(TLT(I)) C TLOG(I) = DLOG(TLT(I))*TDTHT(I) - DLOG(TDB(I))*TLTHT(I) IF(ONE) WRITE (6, 2065) I, IT, TDB(I), TLT(I), 1 TDTHT(I), TLTHT(I), AMR(I), SIGMH(I), SIGMV(I) 2065 FORMAT (' T = ', 2I5, 7G14.5) 5010 CONTINUE C IF(ONE) WRITE (6, 779) DO 60 I = 1, NM 60 AFT(I) = 0.0D0 U = 0.0D0 DC = 0.0D0 DMU = 0.0D0 DTHET = 0.0D0 DTHETA = 0.0D0 DSIGMA = 0.0D0 GERRH = 0.0D0 GERRV = 0.0D0 C NN = NWR RCOEF = (AN/NDEPTH)*THETA*ERRC RECOEF = AMU*RCOEF RECT = RECOEF*TIME_INT TLIMAX = TLT(MMAX) ANUM1 = ANU - 1.0D0 C NCOUNT = 0 C DO 110 JJ = 1, NM J = NM - JJ + 1 JM1 = J - 1 ASJ = AS(J) BSJ = BS(J) CSJ = CS(J) DJ = DEPTH(J) AFT(J) = AL DCJ = 0.0D0 DMUJ = 0.0D0 DTHETJ = 0.0D0 DSIGMHJ = 0.0D0 DSIGMVJ = 0.0D0 GERRHJ = 0.0D0 GERRVJ = 0.0D0 NAF = .FALSE. 95 CONTINUE IF (J.EQ.1) GO TO 111 C DO 100 II = 1, JM1 I = JM1 - II + 1 MI = M(I) c if (mi.lt.m(j)) go to 100 c MIM1 = MI - 1 TLIMIT = TLT(MI) DT = T(J) - T(I) IF (DT.GE.TLIMAX) GO TO 111 IF (DT.GE.TLIMIT) GO TO 100 RLT = RL(MI) ASI = AS(I) BSI = BS(I) CSI = CS(I) C DIST=1.0D0-((ASI-ASJ)**2+(BSI-BSJ)**2+(CSI-CSJ)**2)*.5D0 IF(DIST.LT.RLT) GO TO 100 AMI = DFLOAT(MI)/10.0D0 DI = DEPTH(I) C IDISTV = (DJ - IDU)/NDEPTH + 1.0 IF(IDISTV.LT.1) IDISTV = 1 DISTV = (DI-DJ)**2 C DIST1 = EAR*DACOS(DIST) DIST2 = DIST1*DIST1 RT1 = DIST1**4 C SIGMIH = SIGMH(MI) SIGMH2 = SIGMIH*SIGMIH SIGMIV = SIGMV(MI) SIGMV2 = SIGMIV*SIGMIV IF(RT1.GE.RL4) GO TO 100 L1 = 1 IF(RT1.LE.RM4) GO TO 101 DTEMP = RT1/RM4 HXE = IISHFT(HX(1),-7) - 127 L1 = HXE 101 CONTINUE IF (L1.GT.NSPM) STOP 15 C RE = 0.0D0 TEXPH = DIST2/(2.0D0*SIGMH2) TEXPV = DISTV/(2.0D0*SIGMV2) TEXP = TEXPH + TEXPV IF (TEXP.GT.40.0D0) GO TO 1475 C RE = DIST1*(CMU**MIM1)*TPROD(MI)/ 1 (DT*DT**(THETA)*TDIFF(MI)*RD(L1)* 2 DEXP(TEXP)*SIGMH2*SIGMIV*DFLOAT(IDD(IDISTV))) 1475 CONTINUE C IF (NAF) THEN AMOM(I) = AMOM(I) + AMOM(J)*RE*RECT/AFTJ GO TO 100 END IF C NN = NN + 1 NCOUNT = NCOUNT + 1 IF(NN.LT.NWR. OR..NOT.ONE) GO TO 1450 REC = RE*RECT DISTV1 = DSQRT(DISTV) MIT = MI + MGL - 1 MJ = M(J) + MGL - 1 WRITE (6,1460) 1 I,T(I),X(I),Y(I),DI,MIT,J,T(J),X(J),Y(J),DEPTH(J),MJ, 2 DIST1,DISTV1,REC,DT,L1,IDISTV 1460 FORMAT (' ',2(I5,F12.5,2F8.2,F5.1,I3), 4G11.4, I2, I4) NN = 0 1450 CONTINUE C AFT(J) = AFT(J) + RE*RECT DMUJ = DMUJ + RE DCJ = DCJ + DFLOAT(MIM1)*RE C TLOGT = TLOGL(MI)*TDTHT(MI) - TLOGD(MI)*TLTHT(MI) DTHET = RE*(1.0D0/THETA - DLOG(DT) - TLOGT/TDIFF(MI)) DTHETJ = DTHETJ + DTHET C DST = RE*(10.0D0**(AMI - AMST)) DSIGMNH = (DST/SIGMH2)*(DIST2/SIGMH2 - 2.0D0) DSIGMNV = (DST/SIGMV2)*(DISTV/SIGMV2 - 1.0D0) DSIGMVJ = DSIGMVJ + DSIGMNV DSIGMHJ = DSIGMHJ + DSIGMNH GERRNH = RE*(DIST2/SIGMH2 - 2.0D0)/SIGMH2 GERRNV = RE*(DISTV/SIGMV2-1.0D0)/SIGMV2 GERRVJ = GERRVJ + GERRNV GERRHJ = GERRHJ + GERRNH C IF(NN.LT.2.AND.ONE) 1 WRITE (6, 3010) I, J, MI, DTHET, DTHETA, DSIGMNH, 2 DSIGMA, GERRNH, GERRNV, AFTJ, DT 3010 FORMAT (' I, J, MI, DTHET, DTHETA, DSIGMNH, DSIGMA,' 1 ' GERRNH, GERRNV, AFTJ, DT = ', /, 3I5, 8G12.4) IF(I.GT.NWR2. OR.NCOUNT.GT.100. OR..NOT.ONE) GO TO 100 REC = RE*RECT MIT = MI + MGL - 1 MJ = M(J) + MGL - 1 WRITE (6, 3210) 1 I,T(I),X(I),Y(I),DI,MIT,J,T(J),X(J),Y(J),DEPTH(J),MJ, 2 DIST1,REC,DT,L1, AFTJ 3210 FORMAT (' ',2(I5,F12.5,2F8.2,F5.1,I3),3G11.4, I2, G11.4) 100 CONTINUE 111 CONTINUE C IF (NAF) THEN C AFT(J) = AL/AFTJ PIND8(J) = AL/AFTJ C AMOM(J) = AMOM(J)*AFT(J) AMOM(J) = AMOM(J)*PIND8(J) GO TO 110 END IF C AFTJ = AFT(J) DMU = DMU + DMUJ/AFTJ DC = DC + DCJ/AFTJ DTHETA = DTHETA + DTHETJ/AFTJ DSIGMA = DSIGMA + (DSIGMVJ + DSIGMHJ)/AFTJ GERRH = GERRH + GERRHJ/AFTJ GERRV = GERRV + GERRVJ/AFTJ U = U + DLOG(AFTJ) DALP = DALP + 1.0D0/AFTJ C DNU = DNU + AL/AFTJ C IF (.NOT.FIRST.AND.ONE) THEN NAF = .TRUE. C WRITE (6, 3456) NAF C 3456 FORMAT ('0 ## NAF = ', L4) GO TO 95 END IF C IF(.NOT.ONE) GO TO 110 IF (J.LE.25) GO TO 3350 IF ((NM-J).LE.25) GO TO 3350 IF (MOD(J, 300).NE.1) GO TO 110 3350 WRITE (6, 3440) J, T(J), X(J), Y(J), M(J), AFTJ, U, UU 3440 FORMAT (' ', I5, F20.8, 2F12.5, I8, 3G15.7) 110 CONTINUE C DC = DC*RCOEF*AMU/CMU DMU = DMU*RCOEF DTHETA = DTHETA*RECT DSIGMA = DSIGMA*RECT*SIGMA GERRH = GERRH*RECT*ERRORH GERRV = GERRV*RECT*ERRORV C YC = 0.0D0 TMU =0.0D0 TCU =0.0D0 SUM =0.0D0 SUM = SUM + AMU G = 0.0D0 C YA = 0.0D0 NIWRE = 0 C DO 3510 I = 1, NM MI = M(I) MIM1 = MI - 1 TI = T(I) C TMP = TI + TDB(MI) IF (TEND.LT.TMP) GO TO 3115 TMP = TI + TLT(MI) TIMU = CMU**MIM1 TICU = DFLOAT(MIM1)*TIMU/CMU TMU = TMU + TIMU TCU = TCU + TICU IF (TEND.GT.TMP) THEN YA = YA + SUM*TIMU YC = YC + SUM*TICU G = G - TIMU GO TO 3115 END IF C TEMP1 = TEND - TI TEMP2 = TEMP1**THETA TEMP = (TEMP2 - TDTHT(MI))*TLTHT(MI)/(TDIFF(MI)*TEMP2) IF (NIWRE.LE.10.AND.ONE) 1 WRITE (6, 3400) TEMP1, TEMP2, TEMP, TDTHT(MI), TLTHT(MI) 3400 FORMAT (' TEMP1, TEMP2, TEMP, TDTHT(MI), TLTHT(MI) = ',/, 1 5G15.5) C YA = YA + AMU*TIMU*TEMP TAUL = DLOG(TEMP1) RTE = AMU*TIMU*TEMP YA = YA + RTE DTHETA = DTHETA - 1 RTE*(TLOGL(MI) - TAUL - 2 (TLTHT(MI)*TLOGL(MI)-TDTHT(MI)*TLOGD(MI))/(TLTHT(MI)-TDTHT(MI)) + 3 (TEMP2*TAUL - TDTHT(MI)*TLOGD(MI))/(TEMP2 - TDTHT(MI))) YC = YC + AMU*TICU*TEMP G = G - TEMP*TIMU NIWRE = NIWRE + 1 IF (NIWRE.LE.10.AND.ONE) 1 WRITE (6, 3060) I, IT, MI, TI, TMP, TIMU, TICU, TEMP, 2 AMU, YA, YC 3060 FORMAT (' I,IT,MI,TI,TMP,TIMU,TICU,TEMP,AMU,YA,YC = ',/ 2 3I5, 2F15.6, 7G11.4) 3115 CONTINUE 3510 CONTINUE C SUMLOG2 = 0.0D0 SUMMOM2 = 0.0D0 SUMAL2 = 0.0D0 SUMLOG3 = 0.0D0 SUMMOM3 = 0.0D0 SUMAL3 = 0.0D0 C DO 4400 J = 1, NM c IF (AMOM(J).LE.0.0) THEN c 4350 WRITE (6, 1444) J, T(J), X(J), Y(J), M(J), AFT(J), PIND8(J), c WRITE (6, 1444) J, T(J), X(J), Y(J), M(J), AFT(J), PIND8(J), c 1 AMOM(J) c 1444 FORMAT (' ##', I6, F20.8, 2F12.5, I6, 3D15.6) c GO TO 4400 c END IF C SUMAL2 = SUMAL2 + PIND8(J) SUMMOM2 = SUMMOM2 + AMOM(J)/AMOMC SUMLOG2 = SUMLOG2 + PIND8(J)*DLOG(AMOM(J)/AMOMC) IF (AMOM(J).GE.AMOMC) THEN SUMAL3 = SUMAL3 + PIND8(J) SUMMOM3 = SUMMOM3 + AMOM(J)/AMOMC SUMLOG3 = SUMLOG3 + PIND8(J)*DLOG(AMOM(J)/AMOMC) END IF C IF (J.LE.100) GO TO 4350 C IF ((NM-J).LE.25) GO TO 4350 C IF (MOD(J, 300).NE.1) GO TO 4400 C 4350 WRITE (6, 1440) J, T(J), X(J), Y(J), M(J), AFT(J), AMOM(J) 4400 CONTINUE BETA2 = SUMAL2/SUMLOG2 BVAL2 = BETA2*AMCOF BETA3 = SUMAL3/SUMLOG3 BVAL3 = BETA3*AMCOF SUMDIF = SUMMOM1 - SUMMOM2 WRITE (6, 4455) SUMAL3, SUMLOG3, SUMMOM3, SUMAL2, SUMLOG2, 1 SUMMOM2, SUMAL1, SUMLOG1, SUMMOM1, SUMDIF 4455 FORMAT ('0 SUMAL3, SUMLOG3, SUMMOM3, SUMAL2, SUMLOG2, SUMMOM2,' 1 'SUMAL1, SUMLOG1, SUMMOM1, SUMDIF = ', 4(' ', /, ' ', 3G15.7)) WRITE (6, 4465) BVAL1, BVAL2, BVAL3, BETA1, BETA2, BETA3 4465 FORMAT ('0 BVAL1, BVAL2, BVAL3, BETA1, BETA2, BETA3 = ', /,6G15.7) C TOTALN = NM C F = U + TOTALN*DLOG(ANUM1) - UU*DLOG(ANU) - YA - AL - CONST C GAL = DALP - 1.0D0 GC = TIME_INT*DC - YC FRMRT = F/(TOTALN*DLOG(2.0D0)) C IF (ONE) THEN WRITE (6, 3300) U, UU, YA, AL, TOTALN, CONST 3300 FORMAT ('0 U, UU, YA, AL, TOTALN, CONST = ', 6G15.9, /) WRITE (6, 3007) SIGMAS, ERRORS 3007 FORMAT ('0 SIGMA(START) = ', F12.5, ' ERROR(START) = ', F12.5,/) END IF G = G + DMU*TIME_INT CVL = DLOG10(CMU)*10.0D0 BVL = DLOG10(ANU)*10.0D0 ALRAT = AL/TOTALN WRITE (6, 390) AL, GAL, F, BVL, CVL, ALRAT 390 FORMAT ('0 ALPHA =', F10.4, D12.4, 1 ' LIKE =', G20.10, ' ?', ' B-VALUE = ', F9.5, 2 ' C-VALUE = ', F9.5, ', al/al0 = ', F6.3) C BNU = SUM*CMU**10 WRITE (6, 3315) CMU,GC,TMU,TCU,DC,SUM,BNU 3315 FORMAT (' CMU,GC,TMU,TCU,DC,SUM,BNU', 5D15.6, 2F12.5) TCU = TCU*SUM WRITE (6, 4230) THETA, DTHETA, SIGMA, DSIGMA, ERRORH, GERRH, 1 ERRORV, GERRV 4230 FORMAT (' THETA, DTHETA, SIGMA, DSIGMA, ERRORH, GERRH,' 1 ' ERRORV, GERRV = ', 1 /, 4(2D15.6, 2X)) IAMU = 1.0D5*AMU + 0.5D0 WRITE (6, 410) IAMU, G 410 FORMAT (' MU, G = ', I10, 1P3D12.4) C WRITE (6, 11) C FD = F - OLDF IF(FD) 770, 780, 780 780 IF(FD.GT.EIT) GO TO 785 WRITE (6, 400) 400 FORMAT ('0 SUCCESSFUL END &&') GO TO 500 770 CONTINUE WRITE (6, 4700) CZ, OZ, DZ, SZ 4700 FORMAT ('0 *** CZ = ', 7G14.6, /, ' *** OZ = ', 7G14.6, 1 /, ' *** DZ = ', 7G14.6, /, ' *** SZ = ', 7G14.6, /) 785 CONTINUE CALL TIME (RTIME) C WRITE (6, 480) NITER, NIT, TEMP, RTIME, FRMRT WRITE (6, 480) NITER, NIT, RTIME, FRMRT 480 FORMAT (' ITERATION # ', I4, ' ,CURRENT ITERATION # ', I4, C 1 ' , COMPUTING TIME =', F15.4, ', TIME = ', A8, 1 ', TIME = ', A8, ', Inf. rate = ', F5.2, ' bit/eq') C 1000 CONTINUE NIT = NIT - 1 500 CONTINUE TEMP = CTIME(1) - TBASE CALL TIME (RTIME) WRITE (6, 480) NITER, NIT, RTIME, FRMRT C WRITE (6, INPUT) WRITE (6, 45) WRITE (6, 1580) (NCH(I), I = 1, NPAR) WRITE (6, 45) WRITE (6, 45) NCOUNT IF(.NOT.ONE) GO TO 1704 C 1580 FORMAT (' HLRBS ', 8I5, 2X, 8I5) !HARVARD C 1580 FORMAT (' NLRAD ', 8I5, 2X, 8I5) !NOAA.COR c 1580 FORMAT (' DLR7A ', 8I5, 2X, 8I5) !DUDA 1580 FORMAT (' CIT ', 8I5, 2X, 8I5) !CIT c 1580 FORMAT (' GLR7B ', 8I5, 2X, 8I5) !GREEK WRITE (6, 1580) IAMU 1704 CONTINUE C DO 1787 I = 1, NM 1787 PIND4(I) = PIND8(I) C np = 0 DO 7110 J = 1, NM TLIMIT = TDB(M(j)) RLIMIT = amr(M(j))/25.0 ASJ = AS(J) BSJ = BS(J) CSJ = CS(J) C DO 7100 I = j+1, nm MI = M(I) c if (mi.gt.m(j)) go to 7100 c DT = T(I) - T(J) ASI = AS(I) BSI = BS(I) CSI = CS(I) C DIST=1.0D0-((ASI-ASJ)**2+(BSI-BSJ)**2+(CSI-CSJ)**2)*.5D0 DIST1 = EAR*DACOS(DIST) c IF (DIST1.LT.30.0.and.dt.lt.0.001) pind8(i) = 0.0 c IF (DIST1.LT.30.0.and.dt.lt.tlimit) pind8(i) = 0.0 IF (DIST1.LT.rlimit.and.dt.lt.tlimit) then pind8(i) = 0.0 CALL CDAY(T(I), 2) write (6, 7115) j, i, IYE, IMO, IDA, IH, IMI, 1 m(j), mi, t(i), t(j), dist1, dt, tlimit, rlimit 7115 format (' @', 9i5, 2f15.6, 4f12.3) end if np = np + 1 if (mod(np,1000).eq.1) write (6, 7105) i, j, t(i), t(j), dist1, dt 7105 format (' ^', 2i5, 4f15.6) 7100 continue 7110 continue c WRITE (6, 45) WRITE (6, 45) TIND = 0.0D0 NN = 0 C NTPR = 200 NTPR = 1000 c OPEN (UNIT = 35, STATUS='NEW', DISP='KEEP', c 1 FILE= 'cit32_00_relm.tmp') 1 FILE= 'cit.tmp') OPEN (UNIT = 45, STATUS='NEW', DISP='KEEP', 1 FILE= '[kagan]cit32_00_IND.tmp') !CIT C nrelm = 0 DO 1706 I = 1, NM c ierrc = 3 magerrc = 3 amagerr = 0.3 ilocode = 4 magc = ilocode velerr = 5.0 c CALL CDAY(T(I), 2) c c write (10, 7) i, IYE, IMO, IDA, IH, IMI, SEC, c 1 DT, ALA, ALO, DEPTH, MOM, AMOMT, (EQH(LL), LL = 11, 14) c 7 FORMAT (I6, 4I3, 1X, I2, 1X, F4.1, F7.1, 1X, F6.2, 1X, c 1 F7.2, 1X, F5.1, 1X, I2, 1X, F5.2, 2(1X, I2,1X, I3)) c MT = M(I) + MGL - 1 amt = mt/10.0 mom = dlog10(amom1(i)) !old moment amomt = amom1(i)/10.0**mom !old moment c mom = dlog10(amom(i)) !new moment (redistributed) c amomt = amom(i)/10.0**mom !new moment c write (45, 77) i, IYE-1900, IMO, IDA, IH, IMI, SEC, amt, X(I), 1 Y(I), DEPTH(I), MOM+16, AMOMT, PIND8(I), amom1(i), amom(i), T(I) 77 FORMAT (I6, 4I3, 1X, I2, 1X, F4.1, F7.1, 1X, F6.2, 1X, 1 F7.2, 1X, F5.1, 1X, I2, 1X, F5.2, F10.6, 1P2G13.5, 0PF16.8) c c WRITE (45, 1441) I, IYE, IMO, IDA, IH, IMI, SEC, c 1 T(I), X(I), Y(I), DEPTH(I), MT, AFT(I), AMOM(I), PIND8(I) c TIND = TIND + PIND8(I) IF (I.GT.NTPR) NN = MOD(I,50) IF (NN.NE.1.AND.I.GT.NTPR) GO TO 1706 C if (mt.ge.50) write (19, 708) PIND8(I) !GREEK C 708 FORMAT (F8.5, ',') !GREEK WRITE (6, 1441) I, IYE, IMO, IDA, IH, IMI, SEC, 1 T(I), X(I), Y(I), DEPTH(I), MT, AFT(I), AMOM(I), PIND8(I) 1441 FORMAT (' ', I6,2X,5I4,F7.2,F17.8,2F10.3,F8.2,I4, 2G13.6,F12.8) c IF (x(I).GT.ALATTUP.OR.x(I).Lt.ALATTLO) go to 1706 IF (y(I).Lt.ALONGLO.OR.y(I).GT.ALONGUP) GO TO 1706 c ifprob = 50 ifoc = 0 igcode = 0 iref = 0 icom = 1 c if (err(i).le.0.0) ierrc = 2 if (iye.ge.1970) amagerr = 0.20 if (iye.ge.1977.and.amt.ge.6.0) then amagerr = 0.09 magc = ilocode end if c nrelm = nrelm + 1 temp = abs(err(i))*velerr mcode = PIND8(I)*100.0 + 0.5 itemp = - err(i)/0.1 + 0.5 c WRITE (35, 230) nrelm, iye, imo, ida, ih, imi, X(I), Y(I), 1 ilocode, temp, ierrc, DEPTH(I), amt, magc, amagerr, magerrc, 1 eqh1, ifprob, ifoc, mcode, igcode, iref, sec, qual(itemp) 230 FORMAT (2I5,2I3,1X,2I3,1X,F8.4,F10.4,i3,F6.1,i2,f6.1,f5.2, 1 i3, f5.2, i2, 2(i4, i3, i5), i4, i3, i4, i2, i4, f6.2, 1x, a1) c 230 FORMAT (2I5,2I3,1X,2I3,1X,F7.3,F9.3,i3,F5.1,i2,f6.1,f5.2, c 1 i3, f5.2, i2, 2(i4, i3, i5), i4, 2i3, i4, i2, 2i4) c WRITE (31, 1230) nrelm, iye, imo, ida, ih, imi, ilocode, X(I), 1 Y(I), temp, ierrc, depth(i), amt, magc, amagerr, magerrc, 2 quat1 1230 FORMAT (2I5,2I3,1X,2I3,i3,1X,F7.3,F9.3,F5.1,i2,f6.1,f5.2, 1 i3, f5.2, i2, 1x, 4f8.4) c 1706 CONTINUE WRITE (6, 1707) NM, TIND 1707 FORMAT ('0 NM, TIND = ', I10, F15.5) CLOSE (UNIT = 35) CLOSE (UNIT = 45) C STOP END