C C This program written originally August 1991 calculates probabilities of C earthquake occurrence in a rectangular grid formed by NLONG x NLATT points. C Initial data are HARVARD catalog (its first part from TSTART to TEND). C It calculates weighted average seismic moment around any point on the grid C (LOOP 1100), in a first run (ISW=1) it obtains average moment and c calculates quaternion (QUATR) corresponding to double-couple with this c moment. In the second run (ISW=2) it calculates rotation angle (ANGL) c for earthquake mechanisms surrounding the grid point and then the average c angle (ROTANGL). c c Results are written into several files: HAZH.TMP (WR1=.TRUE.) -- used by C HAZPLOT.FOR to plot risk probabilities and predicted mechanisms, C HAZTBL.TMP dataset (WR2=.TRUE.) to display probabilities and predicted focal C mechanisms in a human-readable form, and HAZMH.TMP -- to be used by C HAZNEWPROB.FOR program to calculate likelihood function and bootstrap C probabilities. C c July 26, 1992, time normalization is applied (time ~ t^TIME_D) C C 100 s CPU time for NW-Pacific 121x121 lattice, EQUAKE C C November 1998 -- program reworked to use kernel specified by Vere-Jones C (1992), slope (ALAMBDA) >= 1.0, source radius (RSOUR) 15-50 km. C C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (NLONG = 121, NLATT = 121) ! NW-Pacific region c PARAMETER (NLONG = 81, NLATT = 41) ! SW-Pacific region RLIM = 350.0D0, drct = 5.0d0 c PARAMETER (NLONG = 161, NLATT = 121) ! SW-Pacific region RLIM = 350.0D0, drct = 5.0d0 c PARAMETER (TIME_D = 0.25D0) PARAMETER (TIME_D = 0.0D0) c PARAMETER (MCUT=580) ! all regions PARAMETER (NTOT=15000) PARAMETER (NMAG = 90) PARAMETER (NDEP=100) C DIMENSION IHAZ(NLONG, NLATT), ALON(NLONG), TPLUSP(3), TMINUSP(3), 1 IDE(NDEP), ITIME(25), IDM(NMAG), quath(4, nlong, nlatt), 2 ROTANGL(nlong, nlatt), baxis(3) c REAL*8 QUAT(4), QUAT1(4), QUAT2(4), QUATR(4), QUATD(4, NTOT) REAL*8 COOR(9), MOMLS(3, 3), MOME(3, 3), EIGV(3, 3), 1 MOMD(3, 3, NTOT), RST(6), RST1(6) REAL*8 EQD(4), T(NTOT), X(NTOT), Y(NTOT) C REAL*4 HAZ(NLONG, NLATT), HAZST(NLONG, NLATT) REAL*4 EQE(8), SEC INTEGER*2 EQH(16), DE(NTOT), MW(NTOT), FPS (4, NTOT) INTEGER*2 EQH1(4) LOGICAL SWI, WR1, WR2, WR3 EQUIVALENCE (EQD(1), EQE(1), EQH(1)) EQUIVALENCE (coor(7), baxis(1)) C COMMON/POL/ PL,GA,SAL,SGA,DAL,DGA,DR1,DR2,MAG(nmag),INT,NC,NIT,NI COMMON /DATE/ SEC, IYE, IMO, IDA, IH, IMI COMMON /MOM/ RAD, PERP C WR1 = .FALSE. WR1 = .TRUE. !For plotting -- e.g., create PostScript file: HAZH.TMP c WR2 = .FALSE. WR2 = .TRUE. !Create probabilities table for display: HAZTBL.TMP c WR3 = .TRUE. WR3 = .FALSE. c timed = time_d if (timed.eq.0.0d0) TIMCD = 1.0D0 TIME_DM1 = TIMED - 1.0D0 SPC = 1.0D03/9.0D0 C HPI = DACOS(0.0D0) PI = 2.0D0*HPI c EAR = 1.0D04/HPI RAD = 90.0D0/HPI c RLIM = 1000.0D0 rsour = 15.0d0 c rsour = 25.0d0 c rsour = 35.0d0 c rsour = 50.0d0 rsour2 = rsour*rsour drct = 100.0d0 !NW_Pacif c drct = 5.0d0 !SW_Pacif hazinit = 1.0D-02 !11/13/98 c alambda = 1.0d0 c alambda = 1.10d0 c alambda = 1.25d0 c alambda = 1.5d0 c c alambda = 1.75d0 c COEF1 = pi*(1.0 - (1.0 + rlim**2/rsour2)**(1.0 - alambda)) 1 *(drct*0.5D0 + 1.0D0)/(alambda - 1.0) if(alambda.eq.1.0) COEF1 = pi*dlog(1.0 + rlim**2/rsour2)* 1 (drct*0.5D0 + 1.0D0) COEFST = 1.0d0/coef1 COEF = (1.0d0 - hazinit)/coef1 c c c ALATTUP=90.000 !WORLDWIDE c ALATTLO=-90.000 !WORLDWIDE c ALONGUP=180.00 !WORLDWIDE c ALONGLO=-180.00 !WORLDWIDE C ALATTUP=60.25 ! NW-Pacific region ALATTLO=-0.25 ! NW-Pacific region ALONGUP=170.25 ! NW-Pacific region ALONGLO=109.75 ! NW-Pacific region c c ALATTUP=0.0 ! SW-Pacific region c ALATTLO=-60.0 ! SW-Pacific region c ALONGUP=-170.0 ! SW-Pacific region c ALONGLO=110.0 ! SW-Pacific region C LONG_FLAG = 0 IF(ALONGLO.GT.ALONGUP) LONG_FLAG=1 c if (long_flag.eq.1) then hazinit = hazinit/(rad*(360.0-ALONGLO+ALONGUP)*(dsind(ALATTUP) - 1 dsind(ALATTLO))*SPC**2) !2/19/94 exact spherical else hazinit = hazinit/(rad*(ALONGUP - ALONGLO)*(dsind(ALATTUP) - 1 dsind(ALATTLO))*SPC**2) !2/19/94 exact spherical end if write (6, 185) rsour, drct, pi, coef, coefst, hazinit 185 format (' rsour, drct, pi, coef, coefst, hazinit ', 6g15.5) C ALATTUP1 = ALATTUP + RLIM/SPC ALATTLO1 = ALATTLO - RLIM/SPC ALONGUP1 = ALONGUP + RLIM/(SPC*DCOSD((ALATTUP+ALATTLO)/2.0D0)) ALONGLO1 = ALONGLO - RLIM/(SPC*DCOSD((ALATTUP+ALATTLO)/2.0D0)) write (6, 155) ALATTUP, ALATTLO, ALONGUP, ALONGLO, 1 ALATTUP1, ALATTLO1, ALONGUP1, ALONGLO1 155 format (' window:', 4f12.3, 2x, 4f12.3) C TSTART = 27395.0D0 !1/1/1977 !DZIEWONSKY'S c TEND = 31412.0D0 !1/1/1988 c TEND = 33969.0D0 !1/1/1995 c TEND = 34334.0D0 !1/1/1996 C TEND = 34700.0D0 !1/1/1997 c TEND = 35430.0D0 !1/1/1999 c TEND = 35795.0D0 !1/1/2000 TEND = 36161.0D0 !1/1/2001 c IDU = 1 !SHALLOW IDL = 70 !SHALLOW C DISTL = RLIM**2 TIME = TEND - TSTART C NRW = 100 NS = 1 ND = (IDL - IDU)/100 + 1 N = 0 C DO 7 I = 1, 100 7 IDE(I) = 0 DO 3 I = 1, NMAG 3 IDM(I) = 0 MMAX = 0 I = 0 amw = 0.0 C WRITE (6, 1) WRITE (6, 15) RLIM, RAD, TSTART, TEND, TIME WRITE (6, 30) C NINS = 0 C OPEN (UNIT=20, STATUS='OLD',DISP='KEEP', READONLY, 1 FORM='UNFORMATTED',RECORDTYPE='FIXED', 2 RECORDSIZE=8, FILE='DISK$DATA:[KAGAN]FPSHPR.BIN') !DZIEWONSKY'S C 5 READ (20, END = 20) EQD I = I + 1 IF(EQD(1).LT.TSTART) GO TO 5 IF(EQD(1).GT.TEND) GO TO 20 MAGW = DFLOAT(EQH(10)-900)*2.0D0/3.0D0 + 0.5D0 !6/17/93 IF (MAGW.LT.MCUT) GO TO 5 IDEP = EQH(9) IF(IDEP.LT.IDU.OR.IDEP.GT.IDL) GO TO 5 c IF (EQE(3).GT.ALATTUP1.OR.EQE(3).LE.ALATTLO1) go to 5 c IF (LONG_FLAG.EQ.1) then if (eqe(4).le.alongup1) go to 55 if (eqe(4).gt.alonglo1) go to 55 go to 5 else IF (EQE(4).LE.ALONGLO1.OR.EQE(4).GT.ALONGUP1) GO TO 5 !REGION EXTENDED end if 55 CONTINUE C IM = (MAGW - mcut)/10.0 + 1.0 !6/17/93 IF(IM.GT.NMAG) IM = NMAG C IYEAR = (EQD(1) - TSTART)/365.25 + 1.0 IF (LONG_FLAG.EQ.0.AND. 1 EQE(3).LT.ALATTUP.AND.EQE(3).GE.ALATTLO.AND. 2 EQE(4).GE.ALONGLO.AND.EQE(4).LT.ALONGUP) then NINS = NINS + 1 ITIME(IYEAR) = ITIME(IYEAR) + 1 IDM(IM) = IDM(IM) + 1 end if c IF (LONG_FLAG.EQ.1.AND. 1 EQE(3).LT.ALATTUP.AND.EQE(3).GE.ALATTLO.AND. 2 (EQE(4).GE.ALONGLO.OR.EQE(4).LT.ALONGUP)) then NINS = NINS + 1 ITIME(IYEAR) = ITIME(IYEAR) + 1 IDM(IM) = IDM(IM) + 1 end if c N = N + 1 IF (MAGW.GT.MMAX) MMAX = MAGW IF(N.GT.NTOT) STOP 12 T(N) = EQD(1) DE(N) = IDEP IDEP = (IDEP - IDU)/ND + 1 IF(IDEP.GT.100) IDEP = 100 X(N) = EQE(3) Y(N) = EQE(4) MW(N) = EQH(10) - 900 - MCUT*1.5 + 1.0 !6/17/93 amw = amw + mw(n)/10.0 IDE(IDEP) = IDE(IDEP) + 1 C DO 40 ISM = 1, 4 EQH1(ISM) = EQH(10 + ISM) 40 FPS(ISM, N) = EQH(10 + ISM) CALL QUATFPS (QUAT1, EQH1, 0) ICODE = 0 CALL BOXTEST (QUAT1, QUAT, QM, ICODE) DO 930 ISM = 1, 4 QUATD(ISM, N) = QUAT(ISM) 930 CONTINUE C CALL MOMAX (MOME, EQH1) DO 940 ISM1 = 1, 3 DO 940 ISM2 = 1, 3 MOMD(ISM1, ISM2, N) = MOME(ISM1, ISM2) 940 CONTINUE C IF(N.GT.25) GO TO 165 CALL CDAY (EQD(1), 2) WRITE (6, 160) I,N,IYE, IMO, IDA, IH, IMI,SEC, EQD(1), 1 (EQE(J),J=3,4),DE(N),(EQH(J),J=10,16), MAGW, MW(N) WRITE (6, 167) QUAT 167 FORMAT (' QUAT = ', 4F12.5) WRITE (6, 169) MOME 169 FORMAT (' MOMENT = ', 9F12.5) 165 CONTINUE IT = MOD(N, NRW) IF(IT.NE.0) GO TO 5 CALL CDAY (EQD(1), 2) WRITE (6, 160) I,N,IYE, IMO, IDA, IH, IMI,SEC, EQD(1), 1 (EQE(J),J=3,4),DE(N),(EQH(J),J=10,16), MAGW, MW(N) GO TO 5 20 CONTINUE CLOSE (UNIT=20) NTEMP = NLONG*NLATT*N/20 NWR = 500 IF (NWR.LT.NTEMP) NWR = NTEMP C WRITE (6, 160) WRITE (6, 162) NWR, NTEMP WRITE (6, 160) CALL CDAY (EQD(1), 2) WRITE (6, 160) I,N,IYE, IMO, IDA, IH, IMI,SEC, EQD(1), 1 (EQE(J),J=3,4),DE(N),(EQH(J),J=10,16), MAGW, MW(N) WRITE (6, 160) amws = amw/n WRITE (6, 168) amws coef = coef/amws write (6, 185) rsour, drct, pi, coef, coefst, hazinit WRITE (6, 1) ITT = I NTT = N C WRITE (6, 133) ITT, NTT, NINS, MMAX WRITE (6, 30) WRITE (6, 30) IDM WRITE (6, 30) WRITE (6, 30) IDE WRITE (6, 30) WRITE (6, 30) IDU, IDL, ND C DO 80 IM = 1, 25 MST = 1 + (IM - 1) CALL POLI (F, IDM, MST, ICODE) ANU = 10.0D0**(0.1D0*GA) GPNU = DR2*GA*DLOG(10.0D0)/10.0D0 MST = MST + MCUT/10.0 - 1 AMST = MST/10.0 IF(ICODE.GE.10) GO TO 81 80 WRITE (6, 380) IM, AMST, NC, GA, DR2, INT 380 FORMAT (' ', I5, F7.1, I7, 2F9.4, I5) 81 CONTINUE C DO 280 I = 2, NDEP IM1 = I - 1 IDE(I) = IDE(I) + IDE(IM1) 280 CONTINUE DO 50 I = 2, NMAG IT= NMAG + 1 - I IT1 = IT + 1 50 IDM(IT) = IDM(IT) + IDM(IT1) C WRITE (6, 30) WRITE (6, 30) IDM WRITE (6, 30) WRITE (6, 30) IDE C IF (WR3) OPEN (UNIT=10, STATUS='NEW', DISP='KEEP', 2 FORM = 'FORMATTED', FILE = 'HARMAG.TMP') DO 243 I = 1, NMAG AMAG = MCUT/100.0 + (I - 1)/100.0 !6/17/93 TEMP = IDM(I)*365.25/TIME IF (WR3.AND.TEMP.NE.0.0D0) WRITE (10, 247) AMAG, TEMP 243 CONTINUE 247 FORMAT (F8.2, 1PE12.4) CLOSE (UNIT=10) C K = NWR NM1 = N - 1 WRITE (6, 1) C ALATTSTEP = (ALATTUP - ALATTLO)/(NLATT-1) if (long_flag.eq.1) then ALONGSTEP = (360.0-ALONGLO+ALONGUP)/(NLONG-1) else ALONGSTEP = (ALONGUP - ALONGLO)/(NLONG-1) end if C nlah = nlatt/2 nloh = nlong/2 c DO 1100 ILATT = 1, NLATT ALATT = ALATTUP - ALATTSTEP*(ILATT - 1) XJ = ALATT DO 1100 ILONG = 1, NLONG ALONG = ALONGLO + ALONGSTEP*(ILONG - 1) IF (ALONG.GT.180.0) ALONG = ALONG - 360.0 YJ = ALONG C DO 1119 IM1 = 1, 3 DO 1119 IM2 = 1, 3 MOMLS (IM1, IM2) = 0.0D0 1119 CONTINUE ANGLALL = 0.0D0 C DO 1100 ISW = 1, 2 SWI = ISW.EQ.2 C DO 100 I = 1, N TI = T(I) TIMD = 1.0d0 if (timed.ne.0.0d0) TIMD = (TEND - TI)**TIME_DM1 !cluster rate XI = X(I) YI = Y(I) IDI = DE(I) DISTX = (XJ - XI)*SPC DISTY = YJ - YI IF (DISTY.GT.300.0) DISTY = 360.0 - DISTY IF (DISTY.LT.-300.0) DISTY = -(360.0 + DISTY) DISTY = DISTY*SPC*DCOSD((XI + XJ)/2.0D0) c DISTZ = JDI - IDI c DISTQ = DISTX*DISTX + DISTY*DISTY + DISTZ*DISTZ DISTQ = DISTX*DISTX + DISTY*DISTY C IF (DISTQ.GT.DISTL) GO TO 100 DIST = DSQRT(DISTQ) C DO 103 ISM = 1, 4 EQH1(ISM) = FPS(ISM, I) QUAT2(ISM) = QUATD(ISM, I) 103 CONTINUE C IF (.NOT.SWI) THEN DO 1033 ISM1 = 1, 3 DO 1033 ISM2 = 1, 3 MOME(ISM1, ISM2) = MOMD(ISM1, ISM2, I) 1033 CONTINUE END IF C AMWI = MW(I)/10.0 !6/17/93 MWI = AMWI + 0.5 !6/17/93 CALL COORD (COOR, TMINUSP, TPLUSP, EQH1) C NNT = NNT + 1 DISTXN = DISTX/DIST DISTYN = DISTY/DIST c TPX = baxis(1)*tminusp(3) - TMINUSP(1)*baxis(3) TPY = baxis(2)*tminusp(3) - TMINUSP(2)*baxis(3) TPNQ = TPX*TPX + TPY*TPY TPN = DSQRT(TPNQ) TPXN = TPX/TPN TPYN = TPY/TPN COSPH = TPXN*DISTXN + TPYN*DISTYN if (alambda.eq.1.0) then HAZN = TIMD*amwi*COEF*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)) HAZS = COEFST*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)) else HAZN = TIMD*amwi*COEF*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)**alambda) HAZS = COEFST*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)**alambda) end if c HAZN = TIMD*amwi*COEF*(drct*COSPH**2 + 1.0D0)/DIST !1/R,clustering c HAZS = COEFST*(drct*COSPH**2 + 1.0D0)/DIST !1/R,clustering C IF (SWI) THEN CALL ROTDET (QUATR, QUAT2, QUAT, QM, ICODE) if (dabs(quat(4)).le.1.0d0) go to 2003 iquat4 = iquat4 + 1 temp = quat(4) - 1.0 quat(4) = 1.0d0 if (iquat4.gt.10) go to 2003 WRITE (6, 166) quat(4), temp 166 format(' *** quat(4) is larger than 1.0', 2g30.20) WRITE (6, 115) NNR,I,IYE,IMO,IDA,IH,IMI,SEC,TI,XI,YI,IDI,MWI,EQH1 write (6, 168) quatr, quat2 write (6, 168) quat 2003 continue ANGL = 2.0D0*DACOSD(QUAT(4)) ANGLALL = ANGLALL + ANGL*HAZN C if (ilatt.ne.nlah.or.ilong.ne.nloh) go to 1120 CALL CDAY (TI, 2) WRITE (6, 115) NNR,I,IYE,IMO,IDA,IH,IMI,SEC,TI,XI,YI,IDI,MWI,EQH1 WRITE (6, 168) QUATR, QUAT2 WRITE (6, 168) QUAT, angl, anglall 1120 CONTINUE C go to 100 END IF C DO 1111 IM1 = 1, 3 DO 1111 IM2 = 1, 3 MOMLS(IM1, IM2) = MOMLS(IM1, IM2) + MOME(IM1, IM2)*HAZN 1111 CONTINUE C HAZ(ILONG, ILATT) = HAZ(ILONG, ILATT) + HAZN HAZST(ILONG, ILATT) = HAZST(ILONG, ILATT) + HAZS C K = K + 1 NNR = NNR + 1 if (ilatt.eq.nlah.and.ilong.eq.nloh) go to 1110 IF (K.LT.NWR) GO TO 100 K = 0 1110 CONTINUE WRITE (6, 30) CALL CDAY (TI, 2) WRITE (6, 115) NNR,I,IYE,IMO,IDA,IH,IMI,SEC,TI,XI,YI,IDI,MWI,EQH1 WRITE (6, 168) COOR c WRITE (6, 168) baxis WRITE (6, 168) MOME WRITE (6, 168) TMINUSP, TPLUSP, XJ, YJ WRITE (6, 142) DIST,DISTX,DISTY, TPN,TPX,TPY, COSPH, HAZN, HAZS 100 CONTINUE C IF (SWI) GO TO 101 RST(1) = MOMLS(1, 1) RST(2) = MOMLS(1, 2) RST(3) = MOMLS(2, 2) RST(4) = MOMLS(1, 3) RST(5) = MOMLS(2, 3) RST(6) = MOMLS(3, 3) DO 1113 IRST = 1, 6 1113 RST1(IRST) = RST(IRST) TR1 = RST(1) + RST(3) + RST(6) R1 = RST(1) - TR1/3.0D0 R3 = RST(3) - TR1/3.0D0 R6 = RST(6) - TR1/3.0D0 ANORM1 = DSQRT(-(R1*R3 + R1*R6 + R3*R6 - RST(2)**2 - 1 RST(4)**2 - RST(5)**2)) DET1 = R1*R3*R6 + 2.0*RST(2)*RST(4)*RST(5) - 2 RST(2)**2*R6 - RST(4)**2*R3 - RST(5)**2*R1 C GAMMA1 = 0.0D0 IF (ANORM1.GT.0.0D0) 1 GAMMA1 = 0.5D0*DET1/((1./3.)*ANORM1**2)**(1.5D0) C CALL DEIGEN (RST, EIGV, 3, 0) CALL QUATMOM5 (QUAT1, 2 EIGV(1,1), EIGV(2,1), EIGV(3,1), 1 EIGV(1,3), EIGV(2,3), EIGV(3,3), 3 EIGV(1,2), EIGV(2,2), EIGV(3,2), ICODE) ICODE = 0 quat1(4) = - quat1(4) ! change sign of rotation to fit QUATFPS CALL BOXTEST (QUAT1, QUATR, QM, ICODE) do 1115 iq = 1, 4 quath(iq, ILONG, ILATT) = QUATR(iq) 1115 continue C if (ilatt.eq.nlah.and.ilong.eq.nloh) then c if (ilatt.le.2.and.ilong.le.2) then WRITE (6, 168) WRITE (6, 168) momls WRITE (6, 168) rst WRITE (6, 168) rst1 WRITE (6, 168) eigv WRITE (6, 168) quatr WRITE (6, 168) ANORM1, DET1, GAMMA1 end if 101 CONTINUE C ROTANGL(ILONG, ILATT) = 0.0 IF (SWI.AND.HAZ(ILONG,ILATT).GT.0.0D0) 1 ROTANGL(ILONG, ILATT) = ANGLALL/HAZ(ILONG, ILATT) C 1100 CONTINUE C j = 0 DO 132 I = 1, N IF (x(i).GT.ALATTUP.OR.x(i).LE.ALATTLO) GO TO 132 C IF (LONG_FLAG.EQ.1) then if (Y(I).le.alongup) go to 135 if (Y(I).gt.alonglo) go to 135 else IF(y(i).LE.ALONGLO.OR.y(i).GT.ALONGUP) GO TO 132 end if 135 CONTINUE C j = j + 1 DO 1132 ILATT = 1, NLATT DO 1132 ILONG = 1, NLONG HAZ(ILONG, ILATT) = HAZ(ILONG, ILATT) + hazinit 1132 CONTINUE 132 CONTINUE write (6, 160) j C WRITE (6, 1) NLONA = 30 NLON = 30 NLONGT = NLONG IF (NLONGT.LT.30) NLON = NLONG IF (NLONGT.LT.30) NLONA = NLONG c DO 110 ILONG = 1, NLONG ALON (ILONG) = ALONGLO + ALONGSTEP*(ILONG - 1) IF (ALON(ILONG).GT.180.0) ALON(ILONG) = ALON(ILONG) - 360.0 DO 110 ILATT = 1, NLATT C IHAZ(ILONG, ILATT) = 0 IF (HAZ(ILONG, ILATT).NE.0.0D0) 1 IHAZ(ILONG, ILATT) = 100.0*ALOG10(1.0e03*HAZ(ILONG, ILATT)) + 0.5 ROTANGL(ILONG, ILATT) = ROTANGL(ILONG, ILATT)/100.0 110 CONTINUE C WRITE (6, 122) (ALON(J), J = 1, NLONA) 122 FORMAT ('0', 12X, 20F6.1, /) c ALATTUP=42.5, ALATTLO=32.5, ALONGUP=-114.00, ALONGLO=-125.00 c hazmax = 0.0d0 ihazm = 1000 DO 131 ILATT = 1, NLATT DO 131 ILONG = 1, NLONG if (hazmax.lt.HAZ(ILONG, ILATT)) hazmax = HAZ(ILONG, ILATT) if (ihazm.gt.IHAZ(ILONG, ILATT)) ihazm = IHAZ(ILONG, ILATT) 131 CONTINUE DO 1131 ILATT = 1, NLATT DO 1131 ILONG = 1, NLONG IHAZ(ILONG, ILATT) = IHAZ(ILONG, ILATT) - ihazm 1131 CONTINUE c tothaz = 0.0 hazav = 0.0 IF (WR1) OPEN (UNIT=10, STATUS='NEW', DISP='KEEP', 2 FORM = 'FORMATTED', FILE = 'HAZH.TMP') c IF (.not.WR2) go to 444 OPEN (UNIT=12, STATUS='NEW', DISP='KEEP', 2 FORM = 'FORMATTED', FILE = 'HAZTBL.TMP') c write (12, 46) 46 format (' Probability Probability ') write (12, 47) 47 format (' M>5.8 M>5.8 ', 1 ' T-axis P-axis ') write (12, 48) 48 format (' eq/day*km^2 eq/y*degr^2 ', 1 'Pl Az Pl Az ') c write (12, 146) 146 format (' Longitude | Area | ', 1 ' | | Rotation ') write (12, 147) 147 format (' | Latitude | square degree | ', 1 'Coef. | | angle ') write (12, 148) 148 format (' v v km^2 v variation v ', 1 ' v degree ',/) 444 continue c DO 111 ILATT = 1, NLATT ALATT = ALATTUP - ALATTSTEP*(ILATT - 1) c ALATT = ALATTUP - ALATTSTEP*ILATT WRITE (6, 120) ILATT, ALATT, (IHAZ(J, ILATT), J = 1, NLON) 120 FORMAT (' ', I3, F8.2, 30I4) C DO 111 ILONG = 1, NLONG ALONG = ALONGLO + ALONGSTEP*(ILONG - 1) IF (ALONG.GT.180.0) ALONG = ALONG - 360.0 IF (TIMED.EQ.0.0D0) 1 hazn = 5.0D08*HAZ(ILONG, ILATT)/time !Poiss hazard, normalized IF (TIMED.NE.0.0D0) 1 hazn = 5.0D08*HAZ(ILONG, ILATT)*TIMED/TIME**TIMED !clustering hazav = hazav + hazn IF (WR1) WRITE (10, 125) ILONG, ILATT, ALONG, ALATT, IHAZ(ILONG, 1 ILATT), HAZ(ILONG, ILATT), hazn, 2 (QUATH(IQ, ILONG, ILATT), IQ = 1, 4), ROTANGL(ILONG, ILATT) c sqarea = rad*(dsind(ALATT+0.5) - dsind(ALATT-0.5))*SPC**2 !2/19/94 exact spherical isqarea = sqarea + 0.5 hazn = 365.25d0*sqarea*HAZ(ILONG, ILATT)/(time) HAZS = -9.99 IF (HAZST(ILONG, ILATT).NE.0.0) 1 HAZS = DSQRT(1.0d0/(sqarea*HAZST(ILONG, ILATT))) if (hazs.gt.9.99) HAZS = 9.99 do 1677 iq = 1, 4 1677 quat(iq) = QUATH(IQ, ILONG, ILATT) call moment4 (quat, mome, eqh1) c IF (WR2) WRITE (12, 145) ALONG, ALATT, HAZ(ILONG, ILATT)/time, 1 isqarea, hazn, HAZS, (eqh1(iq), IQ = 1, 4), 2 1.0d02*ROTANGL(ILONG, ILATT) tothaz = tothaz + hazn c 125 FORMAT (2I5, 2F9.2, I6, 1PG12.3, 0PF9.3, 2X, 4F9.5, 2X, F8.5) 145 FORMAT (F7.1, F6.1,1PG10.2, 0Pi6, F9.5,F6.2,1X, 2(i3,i4,1x), F5.1) 111 CONTINUE CLOSE (UNIT=10) hazav = hazav/(nlatt*nlong) C C Loop 1141 added 15/2/93: in HAZPROB probabilities are needed not C densities, multiplication by COS(LATITUDE) yields proper values C DO 1141 ILATT = 1, NLATT ALATT = ALATTUP - ALATTSTEP*(ILATT - 1) DO 1141 ILONG = 1, NLONG HAZ(ILONG, ILATT) = HAZ(ILONG, ILATT)*dcosd(alatt) 1141 CONTINUE C OPEN (UNIT=10, STATUS='NEW', DISP='KEEP', 1 FORM = 'UNFORMATTED', FILE = 'HAZMH.TMP') WRITE (10) HAZ, rlim, drct, rsour, coef, hazinit, alambda, 1 ALATTUP,ALATTLO,ALONGUP,ALONGLO,MCUT CLOSE (UNIT=10) C write (6, 30) WRITE (6, 15) RLIM, RAD, TSTART, TEND, time write (6, 185) rsour, drct, pi, coef, coefst, hazinit write (6, 155) ALATTUP, ALATTLO, ALONGUP, ALONGLO, 1 ALATTUP1, ALATTLO1, ALONGUP1, ALONGLO1 hazavy = hazav*365.25 WRITE (6, 16) hazav, hazavy, tothaz WRITE (6, 1) call cday (tstart, 2) ITEMP = 0 DO 235 I = 1, 25 ITEMP = ITEMP + ITIME(I) IF (ITIME(I).GT.0) WRITE (30, 134) I, I+IYE-1, ITIME(I), ITEMP 235 CONTINUE C 1 FORMAT ('1') 15 FORMAT ('0RLIM, RAD, TSTART, TEND, TIME_TOTAL = ', 6F15.5) 16 FORMAT ('0Average hazard = ', 2F15.5, ' Total yearly rate = ', 1 4F15.5) 30 FORMAT ('0', 20I6) 115 FORMAT (' ', 3I7, 4I5, F8.2, F16.6, 2F9.3, 2I4, 4X, 2I4, 2X, 2I4) 133 FORMAT (' I, N, NINS, MMAX =', 4I10) 134 FORMAT (' ', 4I6) 140 FORMAT (' ', 10I13) 142 FORMAT (' ', 10G13.5) 160 FORMAT (' ', 3I6, 4I4, F7.2, F15.5, 2F10.3, 6I5, 2I7, 2I3, F6.3) 162 FORMAT (' NWR, NTEMP = ', 2I12) 168 FORMAT (' ', 9F14.7) 178 FORMAT (' ', 9G14.6) 220 FORMAT (' ', 10X, A10, ' CATALOG') 1550 FORMAT (' ', I2, I5, F6.1, F7.2, I7, 1X, 20I5) STOP END C SUBROUTINE ROTDET (QUAT1, QUAT2, QUAT, QM, ICODE) C IMPLICIT REAL*8 (A-H, O-Z) REAL*8 QUAT1(4), QUAT1T(4), QUAT2(4), QUAT(4), QUAT4(4, 4), 1 QUATT(4) REAL*8 QUATIJK(4, 3) /1.0D0, 0.0D0, 0.0D0, 0.0D0, 2 0.0D0, 1.0D0, 0.0D0, 0.0D0, 3 0.0D0, 0.0D0, 1.0D0, 0.0D0/ C QM = 0.0D0 DO 10 ICOD = 1, 4 DO 15 IXC = 1, 4 QUAT1T(IXC) = QUAT1(IXC) 15 CONTINUE IF (ICOD.EQ.4) GO TO 40 DO 30 IXC = 1, 4 QUATT(IXC) = QUATIJK(IXC, ICOD) 30 CONTINUE CALL QUATP (QUATT, QUAT1, QUAT1T) 40 CONTINUE CALL QUATD (QUAT1T, QUAT2, QUAT) DO 20 IXC = 1, 4 QUAT4(IXC, ICOD) = QUAT(IXC) 20 CONTINUE IF (QM.GE.DABS(QUAT(4))) GO TO 10 QM = DABS(QUAT(4)) ICODE = ICOD 10 CONTINUE C TEMP = QUAT4(4, ICODE) DO 50 IXC = 1, 4 IF (TEMP.LE.0.0D0) QUAT(IXC) = - QUAT4(IXC, ICODE) IF (TEMP.GT.0.0D0) QUAT(IXC) = QUAT4(IXC, ICODE) 50 CONTINUE C RETURN END