C C This program written originally November 1991 calculates probabilities of C earthquake occurrence in a rectangular grid formed by NLONG x NLATT points. c The program reads results from HAZ.FOR dataset HAZMH.TMP -- to calculate c likelihood function and bootstrap probabilities. c C Data are HARVARD catalog (its first in time part from TSTART to TEND, the c total number NM; second part from TSTART1 to TEND1, the total number NM1 C It calculates weighted average seismic moment around any earthquake in the c second part of the catalog (LOOP 1100), using the information stored in the c first part of the catalog. In a first run (ISW=1) it obtains c average moment and calculates quaternion (QUATR) corresponding to c double-couple with this moment. In the second run (ISW=2) it calculates c rotation angle (ANGL) for earthquake mechanisms surrounding the grid point c and then the average angle (ROTANGL). This part of the program is similar to c HAZ.FOR, but it calculates the probabilities for epicenters in the second c part of a catalog, not for grid points. c If TSTART1-TEND1 interval is inside TSTART-TEND, then we are essentially c doing leave-one-out (likelihood cross-validation - Silverman, 1986) analysis, c since we do not use pairs with distance=0.0 (see IDIST0). c These loops are for earthquakes, if used for ZONES calculation, change c in TSTART1-TEND1 loop several lines of code (around lines 250-5). c c In the second part of the program it calculates likelihood function for the c second part of the catalog (LOOP 2111), angles of rotation between the c predicted focal mechanisms and real mechanisms (ANGL), whereas ROTANGL1 c lists predicted average rotation angle. c c Finally in LOOP 2200 bootstrap distribution for likelihood function is c calculated: we simulate earthquake distribution with epicenters according to c HAZ(NLONG, NLATT) and calculate likelihood for it comparing the likelihood c with that obtained for real earthquakes in the second part of the catalog. c c Results are written into three files: HAZP.TMP -- for display (?), c LIKHAZ_HARV.TMP -- table of results, and HAZPR.TMP -- to be used by c FIG1.M program to plot likelihood function and bootstrap probabilities. C c CYCLOP C 13 h CPU time (b=check), nprob = 10000, NLONG = 121, NLATT = 121 NW-Pacific region, TEND = 1/1/1989; RLIM = 500, drct = 10; NISH91 C 16 m CPU time (b=check), nprob = 200, NLONG = 121, NLATT = 121 NW-Pacific region, TEND = 1/1/1989; RLIM = 500, drct = 10; NISH91 C c EQUAKE C 10 m CPU time (b=check), nprob = 10000, NLONG = 121, NLATT = 121 NW-Pacific region, TEND = 1/1/1997; RLIM = 350, drct = 100 C IMPLICIT REAL*8 (A-H,O-Z) C PARAMETER (ISEED = 921) PARAMETER (nprarr = 120) c PARAMETER (nprob = 1000) c PARAMETER (nprob = 5000) PARAMETER (nprob = 10000) c PARAMETER (nprob = 100000) PARAMETER (NLONG = 121, NLATT = 121) ! NW-Pacific region, TEND = 31778 -- 1/1/1989; RLIM = 500, drct = 10; NISH91 c PARAMETER (NTOT=15000) PARAMETER (NMAG=80) C DIMENSION IHAZ1(NTOT), TPLUSP(3), TMINUSP(3), 1 IDE(100), IDM(NMAG), quath1(4, NTOT), 2 ROTANGL1(NTOT), baxis(3), iparr(nprarr) 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), X1(NTOT), Y1(NTOT) REAL*8 RAND(128), RANDU(3), WK_ITM(128) C REAL*4 HAZ1(NTOT), HAZ(NLONG, NLATT), HAZT(NLONG, NLATT) REAL*4 EQE(8), SEC INTEGER*2 EQH(16), DE(NTOT), MW(NTOT), FPS1(4, NTOT), FPS(4, NTOT) INTEGER*2 EQH1(4), EQH2(4) LOGICAL SWI EQUIVALENCE (EQD(1), EQE(1), EQH(1)) EQUIVALENCE (coor(7), baxis(1)) C COMMON /DATE/ SEC, IYE, IMO, IDA, IH, IMI COMMON /MOM/ RAD, PERP C HPI = DACOS(0.0D0) PI = 2.0D0*HPI TPI = 2.0D0*PI RAD = 90.0D0/HPI c OPEN (UNIT=10, STATUS='OLD', DISP='KEEP', c 1 FORM = 'UNFORMATTED', FILE = 'disk$user:[kagan]HAZMH.TMP') 1 FORM = 'UNFORMATTED', FILE = 'HAZMH.TMP') read (10) HAZ, rlim, drct, rsour, coef, hazinit, alambda, 1 ALATTUP,ALATTLO,ALONGUP,ALONGLO,MCUT CLOSE (UNIT=10) rsour2 = rsour*rsour write (6, 185) mcut, rsour, drct, pi, coef, hazinit 185 format (' mcut, rsour, drct, pi, coef, hazinit ', i5, 5g15.5) c LONG_FLAG = 0 IF(ALONGLO.GT.ALONGUP) LONG_FLAG=1 IX = ISEED CALL GGU4 (IX, 128, 1, RAND, WK_ITM) SPC = 1.0D03/9.0D0 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 c TSTART1 = 31412.0D0 !1/1/1988 c TSTART1 = 34700.0D0 !1/1/1997 c TSTART1 = 35430.0D0 !1/1/1999 c TSTART1 = 35795.0D0 !1/1/2000 TSTART1 = 36161.0D0 !1/1/2001 c c TEND1 = 35399.0D0 !12/1/1998 c TEND1 = 35611.0D0 !7/1/1999 c TEND1 = 35764.0D0 !12/1/1999 c TEND1 = 35795.0D0 !12/31/1999 c TEND1 = 36161.0d0 !1/1/2001 c TEND1 = 36161.0D0 !1/1/2001 TEND1 = 36526.0D0 !1/1/2002 c TSTART = 27395.0D0 !1/1/1977 !DZIEWONSKY'S TEND = TSTART1 c TEND = 34700.0D0 !1/1/1997 c TEND = 35369.0D0 !11/1/1998 c IDU = 1 !SHALLOW IDL = 70 !SHALLOW C DISTL = RLIM**2 C ARD3 = DACOSD(-1.D0/3.D0) NRW = 100 DEGAC = 1.0D 00 C NS = 1 ND = (IDL - IDU)/100 + 1 C DO 7 I = 1, 100 7 IDE(I) = 0 DO 3 I = 1, NMAG 3 IDM(I) = 0 MMAX = 0 C WRITE (6, 30) WRITE (6, 30) WRITE (6, 15) RLIM, RAD, TSTART, TEND WRITE (6, 30) C tothaz = 0.0d0 hazmax = 0.0d0 DO 131 ILATT = 1, NLATT DO 131 ILONG = 1, NLONG tothaz = tothaz + HAZ(ILONG, ILATT) if (hazmax.lt.HAZ(ILONG, ILATT)) hazmax = HAZ(ILONG, ILATT) HAZT(ILONG, ILATT) = tothaz 131 CONTINUE write (6, 164) tothaz, hazmax 164 format ('0 tothaz, hazmax =', 2f20.5) c N = 0 I = 0 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 2005 READ (20, END = 2020) EQD I = I + 1 IF(EQD(1).LT.TSTART1) GO TO 2005 IF(EQD(1).GT.TEND1) GO TO 2020 c MAGW = DFLOAT(EQH(10)-90)*2.0D0/3.0D0 + 0.5D0 MAGW = DFLOAT(EQH(10)-900)*2.0D0/3.0D0 + 0.5D0 !6/17/93 IF (MAGW.LT.MCUT) GO TO 2005 IDEP = EQH(9) IF(IDEP.LT.IDU.OR.IDEP.GT.IDL) GO TO 2005 IF (EQE(3).GT.ALATTUP.OR.EQE(3).LE.ALATTLO) go to 2005 IF (LONG_FLAG.EQ.1) then if (eqe(4).le.alongup) go to 55 if (eqe(4).gt.alonglo) go to 55 go to 2005 else IF (EQE(4).LE.ALONGLO.OR.EQE(4).GT.ALONGUP) GO TO 2005 end if 55 CONTINUE C n = n + 1 X1(N) = EQE(3) Y1(N) = EQE(4) C DO 2040 ISM = 1, 4 EQH1(ISM) = EQH(10 + ISM) 2040 FPS1(ISM, N) = EQH(10 + ISM) IF(N.GT.100) GO TO 2005 CALL CDAY (EQD(1), 2) WRITE (6, 160) I,N,IYE, IMO, IDA, IH, IMI,SEC, EQD(1), 1 (EQE(J),J=3,4),IDEP,(EQH(J),J=10,16), MAGW go to 2005 2020 CONTINUE CLOSE (UNIT=20) 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),IDEP,(EQH(J),J=10,16), MAGW WRITE (6, 160) nm1 = n c n = 0 i = 0 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 c MAGW = DFLOAT(EQH(10)-90)*2.0D0/3.0D0 + 0.5D0 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 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 255 if (eqe(4).gt.alonglo1) go to 255 go to 5 else IF (EQE(4).LE.ALONGLO1.OR.EQE(4).GT.ALONGUP1) GO TO 5 !REGION EXTENDED end if 255 CONTINUE 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 IM = (MAGW - MCUT)/10.0 + 1.0 !6/17/93 IF(IM.GT.NMAG) IM = NMAG IDM(IM) = IDM(IM) + 1 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.50) 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) C NTEMP = N*N/2500 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) WRITE (6, 1) ITT = I nm = N C WRITE (6, 133) ITT, NM, 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 K = NWR WRITE (6, 1) c DO 1100 i1 = 1, nm1 c DO 1119 IM1 = 1, 3 DO 1119 IM2 = 1, 3 MOMLS (IM1, IM2) = 0.0D0 1119 CONTINUE ANGLALL = 0.0D0 xj = x1(i1) yj = y1(i1) C DO 1100 ISW = 1, 2 SWI = ISW.EQ.2 C DO 100 I = 1, NM TI = T(I) 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) DISTQ = DISTX*DISTX + DISTY*DISTY C IF (DISTQ.GT.DISTL) GO TO 100 DIST = DSQRT(DISTQ) IF (DIST.EQ.0.0D0) THEN if (idist0.lt.9) WRITE (6, 613) DIST 613 FORMAT ('0*** DIST.EQ.0.0D0, DIST = ', F22.16) IDIST0 = IDIST0 + 1 GO TO 100 END IF 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 c HAZN = mwi*COEF*(drct*COSPH**2 + 1.0D0)/DIST if (alambda.eq.1.0) then HAZN = amwi*COEF*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)) else HAZN = amwi*COEF*(drct*COSPH**2 + 1.0D0)/ 1 (rsour2*(1.0 + distq/rsour2)**alambda) end if 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 (i1.ne.nm1/2) 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 HAZ1(i1) = HAZ1(I1) + HAZN C K = K + 1 NNR = NNR + 1 if (i1.eq.nm1/2) 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 WRITE (6, 168) MOME WRITE (6, 168) TMINUSP, TPLUSP, XJ, YJ WRITE (6, 142) DIST,DISTX,DISTY, TPN,TPX,TPY, COSPH, HAZN 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 quath1(iq, i1) = QUATR(iq) 1115 continue C if (i1.eq.nm1/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 ROTANGL1(i1) = 0.0 IF (SWI.AND.HAZ1(I1).GT.0.0D0) 1 ROTANGL1(i1) = ANGLALL/HAZ1(I1) C 1100 CONTINUE C WRITE (6, 1) c DO 110 INM = 1, nm1 HAZ1(INM) = HAZ1(INM) + hazinit IHAZ1(INM) = 0 IF (HAZ1(INM).NE.0.0D0) 1 IHAZ1(INM) = 100.0*ALOG10(HAZ1(INM)) + 0.5 110 CONTINUE c OPEN (UNIT=10, STATUS='new', DISP='KEEP', 1 FORM = 'FORMATTED', FILE = 'HAZP.TMP') c probp = 0.0d0 angle = 0.0d0 angle1 = 0.0d0 DO 2111 INM = 1, nm1 if (haz1(INM).le.0.0) stop 35 DO 2103 ISM = 1, 4 QUATR(ISM) = QUATH1(ISM, INM) EQH2(ISM) = FPS1(ISM, INM) 2103 CONTINUE C CALL QUATFPS (QUAT1, EQH2, 0) ICODE = 0 CALL BOXTEST (QUAT1, QUAT2, QM, ICODE) CALL ROTDET (QUATR, QUAT2, QUAT, QM, ICODE) if (dabs(quat(4)).le.1.0d0) go to 3003 iquat4 = iquat4 + 1 temp = quat(4) - 1.0 quat(4) = 1.0d0 WRITE (6, 166) quat(4), temp write (6, 168) quatr, quat2 write (6, 168) quat 3003 continue ANGL = 2.0D0*DACOSD(QUAT(4)) angle = angle + angl angle1 = angle1 + rotangl1(INM) probp = probp + alog10(haz1(INM)) WRITE (6, 126) INM, x1(INM), y1(INM), IHAZ1(INM), HAZ1(INM), 2 (QUATH1(IQ, INM), IQ = 1, 4), ROTANGL1(INM), angl, probp 126 FORMAT (' ', I6, 2F8.2, i6, F10.6, 2X, 4F7.3, 2X, 2F9.4, 1 1pg16.5) WRITE (6, 128) INM, eqh2, quat1, quat2 128 FORMAT (' ', i6, 2x, 4I6, 4F7.3, 2x, 4F7.3) WRITE (10, 127) INM, x1(INM), y1(INM), IHAZ1(INM), 1 HAZ1(INM), (QUATH1(IQ, INM), IQ = 1, 4), ROTANGL1(INM), 2 angl, probp 127 FORMAT (I6, 2F9.2, i6, F9.6, 2X, 4F9.5, 2X, 2F9.4, 1pg12.4) 2111 CONTINUE CLOSE (UNIT=10) c probm = probp angle = angle/nm1 angle1 = angle1/nm1 WRITE (6, 2206) angle1, angle 2206 FORMAT ('0 average rotation angles, predicted ', f20.10, 1 ', real ', f20.10) WRITE (6, 2205) probm 2205 FORMAT ('0 probm = ', 1pg20.10) c probav = 0.0d0 probsd = 0.0d0 hazlog = nm1*dlog10(hazmax) do 2200 iprob = 1, nprob probp = 0.0d0 do 2100 i = 1, nm1 CALL GGU4 (IX, 3, 0, RANDU, WK_ITM) RANDU(1) = RANDU(1) + RANDU(2)*1.0D-03 + RANDU(3)*1.0D-06 IF(RANDU(1).GT.1.0D0) RANDU(1) = RANDU(1) - 1.0D0 RANDU(1) = RANDU(1)*tothaz c DO 132 ILATT = 1, NLATT DO 132 ILONG = 1, NLONG if (HAZT(ILONG, ILATT).gt.randu(1)) go to 134 132 CONTINUE 134 CONTINUE probp = probp + alog10(HAZ(ILONG, ILATT)) if (iprob.eq.1) write (6, 2110) iprob, i, randu(1), ilong, ilatt, 1 hazt(ilong, ilatt), haz(ilong, ilatt), probp 2110 format(' ', 2i5, f15.6, 2i5, 3f15.6) 2100 continue probp = probp - probm probav = probav + probp probsd = probsd + probp**2 ipr = probp*nprarr/100.0 + nprarr/2 + 0.5 c if (ipr.lt.1) ipr = 1 if (ipr.gt.nprarr) ipr = nprarr iparr(ipr) = iparr(ipr) + 1 if (iprob.gt.10) go to 2200 write (6, 2109) iprob, probp, probav, ipr 2109 format(' ', i5, 2f12.5, i5) 2200 continue c probav = probav/nprob probsd = dsqrt(probsd/nprob - probav**2) WRITE (6, 2201) probav, probsd, hazlog 2201 FORMAT ('0 probav, probsd, hazlog = ', 1p3g20.10) write (6, 160) WRITE (6, 2202) iparr 2202 FORMAT (' iparr = ', 20i6) do 2190 i = 1, nprarr - 1 iparr(i+1) = iparr(i+1) + iparr(i) 2190 continue WRITE (6, 160) WRITE (6, 2202) iparr c OPEN (UNIT=10, STATUS='unknown', DISP='KEEP', ACCESS='APPEND', 1 FORM = 'FORMATTED', FILE = 'likhaz_harv.TMP') write (10, 1127) rlim, drct, rsour, alambda, probm, probav, 1 probsd, angle1, angle, nm1, hazinit, probm/nm1 1127 FORMAT (2F9.2, 2F7.2, 2F9.2, 3F7.2, I5, 1pg12.3, 0PF8.3) CLOSE (UNIT=10) c OPEN (UNIT=10, STATUS='new', DISP='KEEP', 1 FORM = 'FORMATTED', FILE = 'HAZPR.TMP') do 2155 i = 1, nprarr temp = iparr(i)/float(iparr(nprarr)) write (10, 2150) temp 2150 format(f8.4) 2155 continue CLOSE (UNIT=10) C WRITE (6, 1) write (6, 185) mcut, rsour, drct, pi, coef, hazinit WRITE (6, 15) RLIM, RAD, TSTART, TEND WRITE (6, 140) idist0 C 1 FORMAT ('1') 15 FORMAT ('0RLIM, RAD, TSTART, TEND = ', 6F15.5) 30 FORMAT ('0', 20I6) 115 FORMAT (' ', 3I7, 4I5, F8.2, F16.6, 2F9.3, 2I4, 4X, 2I4, 2X, 2I4) 133 FORMAT (' I, N, MMAX =', 3I10) 140 FORMAT (' ', 10I13) 142 FORMAT (' ', 10G13.5) 160 FORMAT (' ', 3I6, 4I4, F7.2, F15.5, 2F10.3, 6I4, 2I6, 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