source: anuga_core/source/anuga/shallow_water/okada.f @ 4263

Last change on this file since 4263 was 4263, checked in by sexton, 17 years ago

earthquake source function

File size: 28.1 KB
Line 
1      subroutine okada(x,y,z,n,depth,xlength,xwidth,dip,rake,slip)
2Cf2py   intent(out) z
3Cf2py   intent(in) x,y,n,depth,dip,rake,slip,xlength,xwidth
4      implicit real*8 (a-h,o-z)
5      real*8 x(n),y(n),z(n),zero
6      data alp,disl3,zero/0.5,0.,0./
7     
8      disl1 = slip*cos(rake)
9      disl2 = slip*sin(rake)
10      cd = cos(dip)
11      sd = sin(dip)
12      do 100 i=1,n
13c         dx = x(i)-x0
14c         dy = y(i)-y0
15c         xr = dx*sin(strike) + dy*cos(strike)
16c         yr = dy*sin(strike) - dx*cos(strike)
17         call srectf(alp,x(i)*.001,y(i)*.001,depth*.001,zero,xlength,
18     &        zero,xwidth,sd,cd,disl1,disl2,disl3,u1,u2,u3,u11,u12,
19     &        u21,u22,u31,u32)
20         z(i) = u3
21 100  continue
22      return
23      end
24
25
26C********************************************************************   00570000
27      SUBROUTINE  SPOINT(ALP,X,Y,DEP,SD,CD,POT1,POT2,POT3,              00530000
28     *                   U1,U2,U3,U11,U12,U21,U22,U31,U32)              00540000
29      IMPLICIT REAL*8 (A-H,O-Z)                                         00550000
30C                                                                       00560000
31C********************************************************************   00570000
32C*****                                                          *****   00580000
33C*****    SURFACE DISPLACEMENT,STRAIN,TILT                      *****   00590000
34C*****    DUE TO BURIED POINT SOURCE IN A SEMIINFINITE MEDIUM   *****   00600000
35C*****                         CODED BY  Y.OKADA ... JAN 1985   *****   00610000
36C*****                                                          *****   00620000
37C********************************************************************   00630000
38C                                                                       00640000
39C***** INPUT                                                            00650000
40C*****   ALP   : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)     00660000
41C*****   X,Y   : COORDINATE OF STATION                                  00670000
42C*****   DEP     : SOURCE DEPTH                                         00680000
43C*****   SD,CD : SIN,COS OF DIP-ANGLE                                   00690000
44C*****          (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)00700000
45C*****   POT1,POT2,POT3 : STRIKE-, DIP- AND TENSILE-POTENCY             00710000
46C*****       POTENCY=(  MOMENT OF DOUBLE-COUPLE  )/MYU    FOR POT1,2    00720000
47C*****       POTENCY=(INTENSITY OF ISOTROPIC PART)/LAMDA  FOR POT3      00730000
48C                                                                       00740000
49C***** OUTPUT                                                           00750000
50C*****   U1, U2, U3      : DISPLACEMENT ( UNIT=(UNIT OF POTENCY) /      00760000
51C*****                   :                     (UNIT OF X,Y,D)**2  )    00770000
52C*****   U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF POTENCY) /      00780000
53C*****   U31,U32         : TILT                (UNIT OF X,Y,D)**3  )    00790000
54C                                                                       00800000
55      DATA  F0,F1,F2,F3,F4,F5,F8,F9                                     00810000
56     *      /0.D0, 1.D0, 2.D0, 3.D0, 4.D0, 5.D0, 8.D0, 9.D0/            00820000
57      PI2=6.283185307179586D0                                           00830000
58C-----                                                                  00840000
59      D =DEP                                                            00850000
60      P =Y*CD + D*SD                                                    00860000
61      Q =Y*SD - D*CD                                                    00870000
62      S =P*SD + Q*CD                                                    00880000
63      X2=X*X                                                            00890000
64      Y2=Y*Y                                                            00900000
65      XY=X*Y                                                            00910000
66      D2=D*D                                                            00920000
67      R2=X2 + Y2 + D2                                                   00930000
68      R =SQRT(R2)                                                       00940000
69      R3=R *R2                                                          00950000
70      R5=R3*R2                                                          00960000
71      QR=F3*Q/R5                                                        00970000
72      XR =F5*X2/R2                                                      00980000
73      YR =F5*Y2/R2                                                      00990000
74      XYR=F5*XY/R2                                                      01000000
75      DR =F5*D /R2                                                      01010000
76      RD =R + D                                                         01020000
77      R12=F1/(R*RD*RD)                                                  01030000
78      R32=R12*(F2*R + D)/ R2                                            01040000
79      R33=R12*(F3*R + D)/(R2*RD)                                        01050000
80      R53=R12*(F8*R2 + F9*R*D + F3*D2)/(R2*R2*RD)                       01060000
81      R54=R12*(F5*R2 + F4*R*D +    D2)/R3*R12                           01070000
82C-----                                                                  01080000
83      A1= ALP*Y*(R12-X2*R33)                                            01090000
84      A2= ALP*X*(R12-Y2*R33)                                            01100000
85      A3= ALP*X/R3 - A2                                                 01110000
86      A4=-ALP*XY*R32                                                    01120000
87      A5= ALP*( F1/(R*RD) - X2*R32 )                                    01130000
88      B1= ALP*(-F3*XY*R33      + F3*X2*XY*R54)                          01140000
89      B2= ALP*( F1/R3 - F3*R12 + F3*X2*Y2*R54)                          01150000
90      B3= ALP*( F1/R3 - F3*X2/R5) - B2                                  01160000
91      B4=-ALP*F3*XY/R5 - B1                                             01170000
92      C1=-ALP*Y*(R32 - X2*R53)                                          01180000
93      C2=-ALP*X*(R32 - Y2*R53)                                          01190000
94      C3=-ALP*F3*X*D/R5 - C2                                            01200000
95C-----                                                                  01210000
96      U1 =F0                                                            01220000
97      U2 =F0                                                            01230000
98      U3 =F0                                                            01240000
99      U11=F0                                                            01250000
100      U12=F0                                                            01260000
101      U21=F0                                                            01270000
102      U22=F0                                                            01280000
103      U31=F0                                                            01290000
104      U32=F0                                                            01300000
105C======================================                                 01310000
106C=====  STRIKE-SLIP CONTRIBUTION  =====                                 01320000
107C======================================                                 01330000
108      IF(POT1.EQ.F0)   GO TO 200                                        01340000
109      UN=POT1/PI2                                                       01350000
110      QRX=QR*X                                                          01360000
111      FX=F3*X/R5*SD                                                     01370000
112      U1 =U1 - UN*( QRX*X + A1*SD )                                     01380000
113      U2 =U2 - UN*( QRX*Y + A2*SD )                                     01390000
114      U3 =U3 - UN*( QRX*D + A4*SD )                                     01400000
115      U11=U11- UN*( QRX* (F2-XR)        + B1*SD )                       01410000
116      U12=U12- UN*(-QRX*XYR      + FX*X + B2*SD )                       01420000
117      U21=U21- UN*( QR*Y*(F1-XR)        + B2*SD )                       01430000
118      U22=U22- UN*( QRX *(F1-YR) + FX*Y + B4*SD )                       01440000
119      U31=U31- UN*( QR*D*(F1-XR)        + C1*SD )                       01450000
120      U32=U32- UN*(-QRX*DR*Y     + FX*D + C2*SD )                       01460000
121C======================================                                 01470000
122C=====    DIP-SLIP CONTRIBUTION   =====                                 01480000
123C======================================                                 01490000
124  200 IF(POT2.EQ.F0)   GO TO 300                                        01500000
125      UN=POT2/PI2                                                       01510000
126      SDCD=SD*CD                                                        01520000
127      QRP=QR*P                                                          01530000
128      FS=F3*S/R5                                                        01540000
129      U1 =U1 - UN*( QRP*X - A3*SDCD )                                   01550000
130      U2 =U2 - UN*( QRP*Y - A1*SDCD )                                   01560000
131      U3 =U3 - UN*( QRP*D - A5*SDCD )                                   01570000
132      U11=U11- UN*( QRP*(F1-XR)        - B3*SDCD )                      01580000
133      U12=U12- UN*(-QRP*XYR     + FS*X - B1*SDCD )                      01590000
134      U21=U21- UN*(-QRP*XYR            - B1*SDCD )                      01600000
135      U22=U22- UN*( QRP*(F1-YR) + FS*Y - B2*SDCD )                      01610000
136      U31=U31- UN*(-QRP*DR*X           - C3*SDCD )                      01620000
137      U32=U32- UN*(-QRP*DR*Y    + FS*D - C1*SDCD )                      01630000
138C========================================                               01640000
139C=====  TENSILE-FAULT CONTRIBUTION  =====                               01650000
140C========================================                               01660000
141  300 IF(POT3.EQ.F0)   GO TO 900                                        01670000
142      UN=POT3/PI2                                                       01680000
143      SDSD=SD*SD                                                        01690000
144      QRQ=QR*Q                                                          01700000
145      FQ=F2*QR*SD                                                       01710000
146      U1 =U1 + UN*( QRQ*X - A3*SDSD )                                   01720000
147      U2 =U2 + UN*( QRQ*Y - A1*SDSD )                                   01730000
148      U3 =U3 + UN*( QRQ*D - A5*SDSD )                                   01740000
149      U11=U11+ UN*( QRQ*(F1-XR)        - B3*SDSD )                      01750000
150      U12=U12+ UN*(-QRQ*XYR     + FQ*X - B1*SDSD )                      01760000
151      U21=U21+ UN*(-QRQ*XYR            - B1*SDSD )                      01770000
152      U22=U22+ UN*( QRQ*(F1-YR) + FQ*Y - B2*SDSD )                      01780000
153      U31=U31+ UN*(-QRQ*DR*X           - C3*SDSD )                      01790000
154      U32=U32+ UN*(-QRQ*DR*Y    + FQ*D - C1*SDSD )                      01800000
155C-----                                                                  01810000
156  900 RETURN                                                            01820000
157      END                                                               01830000
158      SUBROUTINE  SRECTF(ALP,X,Y,DEP,AL1,AL2,AW1,AW2,                   01840000
159     *                   SD,CD,DISL1,DISL2,DISL3,                       01850000
160     *                   U1,U2,U3,U11,U12,U21,U22,U31,U32)              01860000
161      IMPLICIT REAL*8 (A-H,O-Z)                                         01870000
162C                                                                       01880000
163C*********************************************************              01890000
164C*****                                               *****              01900000
165C*****    SURFACE DISPLACEMENT,STRAIN,TILT           *****              01910000
166C*****    DUE TO RECTANGULAR FAULT IN A HALF-SPACE   *****              01920000
167C*****              CODED BY  Y.OKADA ... JAN 1985   *****              01930000
168C*****                                               *****              01940000
169C*********************************************************              01950000
170C                                                                       01960000
171C***** INPUT                                                            01970000
172C*****   ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)   01980000
173C*****   X,Y     : COORDINATE OF STATION                                01990000
174C*****   DEP     : SOURCE DEPTH                                         02000000
175C*****   AL1,AL2 : FAULT LENGTH RANGE                                   02010000
176C*****   AW1,AW2 : FAULT WIDTH RANGE                                    02020000
177C*****   SD,CD   : SIN,COS OF DIP-ANGLE                                 02030000
178C*****          (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)02040000
179C*****   DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION      02050000
180C                                                                       02060000
181C***** OUTPUT                                                           02070000
182C*****   U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL     )      02080000
183C*****   U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /          02090000
184C*****   U31,U32         : TILT                 UNIT OF X,Y,,,AW )      02100000
185C                                                                       02110000
186C***** SUBROUTINE USED...SRECTG                                         02120000
187C                                                                       02130000
188      DIMENSION  U(9),DU(9)                                             02140000
189      DATA  F0, F1 / 0.D0, 1.D0 /                                       02150000
190C-----                                                                  02160000
191      P = Y*CD + DEP*SD                                                 02170000
192      Q = Y*SD - DEP*CD                                                 02180000
193      DO 1111  I=1,9                                                    02190000
194 1111 U(I)=F0                                                           02200000
195C-----                                                                  02210000
196      DO 5555  K=1,2                                                    02220000
197       IF(K.EQ.1)  ET=P-AW1                                             02230000
198       IF(K.EQ.2)  ET=P-AW2                                             02240000
199       DO 4444  J=1,2                                                   02250000
200        IF(J.EQ.1)  XI=X-AL1                                            02260000
201        IF(J.EQ.2)  XI=X-AL2                                            02270000
202        JK=J+K                                                          02280000
203        IF(JK.NE.3)  SIGN= F1                                           02290000
204        IF(JK.EQ.3)  SIGN=-F1                                           02300000
205        CALL SRECTG(ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3,                02310000
206     *           DU(1),DU(2),DU(3),DU(4),DU(5),DU(6),DU(7),DU(8),DU(9)) 02320000
207        DO 3333  I=1,9                                                  02330000
208         U(I)=U(I)+SIGN*DU(I)                                           02340000
209 3333   CONTINUE                                                        02350000
210 4444  CONTINUE                                                         02360000
211 5555 CONTINUE                                                          02370000
212      U1 =U(1)                                                          02380000
213      U2 =U(2)                                                          02390000
214      U3 =U(3)                                                          02400000
215      U11=U(4)                                                          02410000
216      U12=U(5)                                                          02420000
217      U21=U(6)                                                          02430000
218      U22=U(7)                                                          02440000
219      U31=U(8)                                                          02450000
220      U32=U(9)                                                          02460000
221      RETURN                                                            02470000
222      END                                                               02480000
223      SUBROUTINE  SRECTG(ALP,XI,ET,Q,SD,CD,DISL1,DISL2,DISL3,           02490000
224     *                   U1,U2,U3,U11,U12,U21,U22,U31,U32)              02500000
225      IMPLICIT REAL*8 (A-H,O-Z)                                         02510000
226C                                                                       02520000
227C*********************************************************************  02530000
228C*****                                                           *****  02540000
229C*****  INDEFINITE INTEGRAL OF SURFACE DISPLACEMENT,STRAIN,TILT  *****  02550000
230C*****  DUE TO RECTANGULAR FAULT IN A HALF-SPACE                 *****  02560000
231C*****                          CODED BY  Y.OKADA ... JAN 1985   *****  02570000
232C*****                                                           *****  02580000
233C*********************************************************************  02590000
234C                                                                       02600000
235C***** INPUT                                                            02610000
236C*****   ALP     : MEDIUM CONSTANT  MYU/(LAMDA+MYU)=1./((VP/VS)**2-1)   02620000
237C*****   XI,ET,Q : FAULT COORDINATE                                     02630000
238C*****   SD,CD   : SIN,COS OF DIP-ANGLE                                 02640000
239C*****          (CD=0.D0, SD=+/-1.D0 SHOULD BE GIVEN FOR VERTICAL FAULT)02650000
240C*****   DISL1,DISL2,DISL3 : STRIKE-, DIP- AND TENSILE-DISLOCATION      02660000
241C                                                                       02670000
242C***** OUTPUT                                                           02680000
243C*****   U1, U2, U3      : DISPLACEMENT ( UNIT= UNIT OF DISL    )       02690000
244C*****   U11,U12,U21,U22 : STRAIN       ( UNIT= UNIT OF DISL /          02700000
245C*****   U31,U32         : TILT                 UNIT OF XI,ET,Q )       02710000
246C                                                                       02720000
247      DATA  F0,F1,F2/ 0.D0, 1.D0, 2.D0 /                                02730000
248      PI2=6.283185307179586D0                                           02740000
249C-----                                                                  02750000
250      XI2=XI*XI                                                         02760000
251      ET2=ET*ET                                                         02770000
252      Q2=Q*Q                                                            02780000
253      R2=XI2+ET2+Q2                                                     02790000
254      R =DSQRT(R2)                                                      02800000
255      R3=R*R2                                                           02810000
256      D =ET*SD-Q*CD                                                     02820000
257      Y =ET*CD+Q*SD                                                     02830000
258      RET=R+ET                                                          02840000
259      IF(RET.LT.F0)  RET=F0                                             02850000
260      RD =R+D                                                           02860000
261      RRD=F1/(R*RD)                                                     02870000
262C-----                                                                  02880000
263      IF( Q .NE.F0)  TT = DATAN( XI*ET/(Q*R) )                          02890000
264      IF( Q .EQ.F0)  TT = F0                                            02900000
265      IF(RET.NE.F0)  RE = F1/RET                                        02910000
266      IF(RET.EQ.F0)  RE = F0                                            02920000
267      IF(RET.NE.F0)  DLE= DLOG(RET)                                     02930000
268      IF(RET.EQ.F0)  DLE=-DLOG(R-ET)                                    02940000
269      RRX=F1/(R*(R+XI))                                                 02950000
270      RRE=RE/R                                                          02960000
271      AXI=(F2*R+XI)*RRX*RRX/R                                           02970000
272      AET=(F2*R+ET)*RRE*RRE/R                                           02980000
273      IF(CD.EQ.F0)  GO TO 20                                            02990000
274C==============================                                         03000000
275C=====   INCLINED FAULT   =====                                         03010000
276C==============================                                         03020000
277      TD=SD/CD                                                          03030000
278      X =DSQRT(XI2+Q2)                                                  03040000
279      IF(XI.EQ.F0)  A5=F0                                               03050000
280      IF(XI.NE.F0)                                                      03060000
281     *A5= ALP*F2/CD*DATAN( (ET*(X+Q*CD)+X*(R+X)*SD) / (XI*(R+X)*CD) )   03070000
282      A4= ALP/CD*( DLOG(RD) - SD*DLE )                                  03080000
283      A3= ALP*(Y/RD/CD - DLE) + TD*A4                                   03090000
284      A1=-ALP/CD*XI/RD        - TD*A5                                   03100000
285      C1= ALP/CD*XI*(RRD - SD*RRE)                                      03110000
286      C3= ALP/CD*(Q*RRE - Y*RRD)                                        03120000
287      B1= ALP/CD*(XI2*RRD - F1)/RD - TD*C3                              03130000
288      B2= ALP/CD*XI*Y*RRD/RD       - TD*C1                              03140000
289      GO TO 30                                                          03150000
290C==============================                                         03160000
291C=====   VERTICAL FAULT   =====                                         03170000
292C==============================                                         03180000
293   20 RD2=RD*RD                                                         03190000
294      A1=-ALP/F2*XI*Q/RD2                                               03200000
295      A3= ALP/F2*( ET/RD + Y*Q/RD2 - DLE )                              03210000
296      A4=-ALP*Q/RD                                                      03220000
297      A5=-ALP*XI*SD/RD                                                  03230000
298      B1= ALP/F2*  Q  /RD2*(F2*XI2*RRD - F1)                            03240000
299      B2= ALP/F2*XI*SD/RD2*(F2*Q2 *RRD - F1)                            03250000
300      C1= ALP*XI*Q*RRD/RD                                               03260000
301      C3= ALP*SD/RD*(XI2*RRD - F1)                                      03270000
302C-----                                                                  03280000
303   30 A2=-ALP*DLE - A3                                                  03290000
304      B3=-ALP*XI*RRE - B2                                               03300000
305      B4=-ALP*( CD/R + Q*SD*RRE ) - B1                                  03310000
306      C2= ALP*(-SD/R + Q*CD*RRE ) - C3                                  03320000
307C-----                                                                  03330000
308      U1 =F0                                                            03340000
309      U2 =F0                                                            03350000
310      U3 =F0                                                            03360000
311      U11=F0                                                            03370000
312      U12=F0                                                            03380000
313      U21=F0                                                            03390000
314      U22=F0                                                            03400000
315      U31=F0                                                            03410000
316      U32=F0                                                            03420000
317C======================================                                 03430000
318C=====  STRIKE-SLIP CONTRIBUTION  =====                                 03440000
319C======================================                                 03450000
320      IF(DISL1.EQ.F0)  GO TO 200                                        03460000
321      UN=DISL1/PI2                                                      03470000
322      REQ=RRE*Q                                                         03480000
323      U1 =U1 - UN*( REQ*XI +   TT    + A1*SD )                          03490000
324      U2 =U2 - UN*( REQ*+ Q*CD*RE + A2*SD )                          03500000
325      U3 =U3 - UN*( REQ*+ Q*SD*RE + A4*SD )                          03510000
326      U11=U11+ UN*( XI2*Q*AET - B1*SD )                                 03520000
327      U12=U12+ UN*( XI2*XI*( D/(ET2+Q2)/R3 - AET*SD ) - B2*SD )         03530000
328      U21=U21+ UN*( XI*Q/R3*CD + (XI*Q2*AET - B2)*SD )                  03540000
329      U22=U22+ UN*( Y *Q/R3*CD + (Q*SD*(Q2*AET-F2*RRE)                  03550000
330     *                            -(XI2+ET2)/R3*CD - B4)*SD )           03560000
331      U31=U31+ UN*(-XI*Q2*AET*CD + (XI*Q/R3 - C1)*SD )                  03570000
332      U32=U32+ UN*( D*Q/R3*CD + (XI2*Q*AET*CD - SD/R + Y*Q/R3 - C2)*SD )03580000
333C===================================                                    03590000
334C=====  DIP-SLIP CONTRIBUTION  =====                                    03600000
335C===================================                                    03610000
336  200 IF(DISL2.EQ.F0)  GO TO 300                                        03620000
337      UN=DISL2/PI2                                                      03630000
338      SDCD=SD*CD                                                        03640000
339      U1 =U1 - UN*( Q/R             - A3*SDCD )                         03650000
340      U2 =U2 - UN*( Y*Q*RRX + CD*TT - A1*SDCD )                         03660000
341      U3 =U3 - UN*( D*Q*RRX + SD*TT - A5*SDCD )                         03670000
342      U11=U11+ UN*( XI*Q/R3            + B3*SDCD )                      03680000
343      U12=U12+ UN*( Y *Q/R3 - SD/R     + B1*SDCD )                      03690000
344      U21=U21+ UN*( Y *Q/R3 + Q*CD*RRE + B1*SDCD )                      03700000
345      U22=U22+ UN*( Y*Y*Q*AXI - (F2*Y*RRX + XI*CD*RRE)*SD + B2*SDCD )   03710000
346      U31=U31+ UN*( D *Q/R3 + Q*SD*RRE + C3*SDCD )                      03720000
347      U32=U32+ UN*( Y*D*Q*AXI - (F2*D*RRX + XI*SD*RRE)*SD + C1*SDCD )   03730000
348C========================================                               03740000
349C=====  TENSILE-FAULT CONTRIBUTION  =====                               03750000
350C========================================                               03760000
351  300 IF(DISL3.EQ.F0)  GO TO 900                                        03770000
352      UN=DISL3/PI2                                                      03780000
353      SDSD=SD*SD                                                        03790000
354      U1 =U1 + UN*( Q2*RRE                       - A3*SDSD )            03800000
355      U2 =U2 + UN*(-D*Q*RRX - SD*(XI*Q*RRE - TT) - A1*SDSD )            03810000
356      U3 =U3 + UN*( Y*Q*RRX + CD*(XI*Q*RRE - TT) - A5*SDSD )            03820000
357      U11=U11- UN*( XI*Q2*AET             + B3*SDSD )                   03830000
358      U12=U12- UN*(-D*Q/R3 - XI2*Q*AET*SD + B1*SDSD )                   03840000
359      U21=U21- UN*( Q2*(CD/R3 + Q*AET*SD) + B1*SDSD )                   03850000
360      U22=U22- UN*((Y*CD-D*SD)*Q2*AXI - F2*Q*SD*CD*RRX                  03860000
361     *                      - (XI*Q2*AET - B2)*SDSD )                   03870000
362      U31=U31- UN*( Q2*(SD/R3 - Q*AET*CD) + C3*SDSD )                   03880000
363      U32=U32- UN*((Y*SD+D*CD)*Q2*AXI + XI*Q2*AET*SD*CD                 03890000
364     *                       - (F2*Q*RRX - C1)*SDSD )                   03900000
365C-----                                                                  03910000
366  900 RETURN                                                            03920000
367      END                                                               03930000
Note: See TracBrowser for help on using the repository browser.