Computes time evolution function for seismic source.
37
38 implicit none
39
40 integer*4 :: nb_fnc,id_fnc,i,nb_data_fnc
41 integer*4 :: np0,np1,np2,np_current
42
43 integer*4, dimension(nb_fnc) :: type_fnc
44 integer*4, dimension(nb_fnc +1) :: ind_fnc
45
46 real*8 :: val,pi,t_t0,t0,t1,v0,v1
47 real*8 :: tau, scaling, hdur
48 real*8 :: amp
49 real*8 :: tau_new
50 real*8 :: time,beta2,t0_delay
51 real*8 :: dt,tr,te,tp,psv,sum1,ts,svi
52
53 real*8, dimension(nb_data_fnc) :: data_fnc
54
55 real*8, dimension(:), allocatable :: svf, int_svf, integ_svf
56
58 pi = 4.0d0 * datan(1.0d0)
59
60 select case (type_fnc(id_fnc))
61
62 case(0)
64
65 case(1)
66 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
68 * dexp(-1.0d0*data_fnc(ind_fnc(id_fnc))
69
70 case(2)
71 pi = 4.0d0 * datan(1.0d0)
72 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
74 * dexp(-0.5d0*data_fnc(ind_fnc(id_fnc))
75 * data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
76
77 case(3)
78 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
79 t0 = data_fnc(i) - t0_delay
80 t1 = data_fnc(i +2) - t0_delay
81 v0 = data_fnc(i +1)
82 v1 = data_fnc(i +3)
83 if ((time.ge.t0).and.(time.le.t1)) &
85 enddo
86
87 case(4)
88 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
89 beta2=data_fnc(ind_fnc(id_fnc))
91 * (-3.0d0 + 2.0d0*beta2*t_t0*t_t0) &
92 * dexp(-beta2*t_t0*t_t0)
93
94 case(6)
95 pi = 4.0d0 * datan(1.0d0)
96 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
98 * dexp(-0.5d0*4.d0*pi*pi*data_fnc(ind_fnc(id_fnc)) &
99 * data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
101
102 case(13)
103 tau = data_fnc(ind_fnc(id_fnc))
104 scaling = data_fnc(ind_fnc(id_fnc) +1)
105 t0 = 2.0 * tau
106 hdur = tau/2.0
107
108 t_t0 = time - t0 - t0_delay
110
111 case(14)
112 tau = data_fnc(ind_fnc(id_fnc))
113 tau = tau*0.25;
114 amp = data_fnc(ind_fnc(id_fnc) +1)
115 t_t0 = time - t0_delay
116
117
118 if (t_t0 .lt. 0.0d0) then
120
121
122 else
124 endif
125
126
127
128 case(31)
129
130
131
132
133 if (time.lt.t0_delay) then
135 elseif ((time.ge.t0_delay).and.(time.le.(t0_delay+tau_new)))then
136 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
137 t0 = data_fnc(i) + t0_delay
138 t1 = data_fnc(i +2) + t0_delay
139 v0 = data_fnc(i +1)
140 v1 = data_fnc(i +3)
141 if ((time.ge.t0).and.(time.le.t1)) then
143 endif
144 enddo
145 elseif (time.gt.(t0_delay+tau_new)) then
147 endif
148
149 case(32)
150 if (time.lt.t0_delay) then
152 elseif ((time.ge.t0_delay).and.(time.le.(t0_delay+tau_new)))then
154 elseif (time.gt.(t0_delay+tau_new)) then
156 endif
157
158
159 case(33)
160
161
162
163
164
165
166
167 if (time.lt.t0_delay) then
169 elseif ((time.ge.t0_delay).and.(time.le.(t0_delay+tau_new)))then
170 do i = ind_fnc(id_fnc),ind_fnc(id_fnc+1) -3,2
171
172 t0 = data_fnc(i)*tau_new + t0_delay
173 t1 = data_fnc(i +2)*tau_new + t0_delay
174 v0 = data_fnc(i +1)
175 v1 = data_fnc(i +3)
176
177
178
179 if ((time.ge.t0).and.(time.le.t1)) then
181 endif
182 enddo
183 elseif (time.gt.(t0_delay+tau_new)) then
185 endif
186
187
188
189 case(50)
190 tau = tau_new
191 scaling = data_fnc(ind_fnc(id_fnc) +1)
192 t0 = 2.0 * tau
193 hdur = tau/2.0
194
195 t_t0 = time - t0 - t0_delay
197
198
199 case(55)
200
201
202
203
204
205
206
207 if (time .le. t0_delay) then
209 elseif ( time .ge. t0_delay + tau_new) then
211 else
212 dt = 0.001;
213 tr = tau_new;
214 te = 0.8*tr;
215 tp = 0.2*tr;
216
217 np0 = int(tp/dt+1.0)
218 np1 = int(te/dt+1.0)
219 np2 = int(tr/dt+1.0)
220 np_current = int((time-t0_delay)/dt+1.0)
221
222 if (np_current .gt. np2) then
223 write(*,*) 'Error in Archuleta function'
224 write(*,*) np2, np_current
225 stop
226 endif
227
228 psv = sqrt(1.+100./(np0*dt))
229 sum1 =0.0
230 allocate(svf(1:np2));
231
232 do i = 1, np0
233 ts = (i-1)*dt
234 svi = ts*psv/tp*sin(0.5*pi/tp*ts)
235 svf(i) = svi
236 sum1 = sum1+svi
237 enddo
238 do i = np0 + 1, np1
239 ts = (i-1)*dt
240 svi = sqrt(1.+100./ts)
241 svf(i) = svi
242 sum1 = sum1 + svi
243 enddo
244 do i = np1 + 1, np2
245 ts = (i-1)*dt
246 svi = sqrt(1.+100./ts)*sin((np2-i)*dt*pi*0.5/(tr-te))
247 svf(i) = svi
248 sum1 = sum1 + svi
249 enddo
250 sum1 = sum1 * dt
251
252
253
254
255 svf = svf/sum1;
256
257 allocate(int_svf(1:np2)); int_svf = 0.d0; int_svf(1) = 0;
258 allocate(integ_svf(1:np2)); integ_svf = 0.d0; integ_svf(1
259
260 do i = 2, np2
261 int_svf(i) = 0.5*(svf(i)+svf(i-1))*dt;
262
263
264 integ_svf(i) = sum(int_svf(2:i));
265 enddo
266
269 write(*,*) 'Error in Archuleta slip function '
271 stop
272 endif
273 deallocate(svf,int_svf,integ_svf)
274
275
276
277 endif
278
279
280
281 case(99)
282 tau = data_fnc(ind_fnc(id_fnc))
283 amp = data_fnc(ind_fnc(id_fnc) +1)
284 t_t0 = time
286 *(t_t0 - tau)/amp)**2.0d0) * cos(((2.0d0
287 *(t_t0 - tau) + 2.0d0*dasin(1.0d0)/2.0d0
288
289 case default
291
292 end select
293
294
295 return
296
real *8 function get_func_value_sism(nb_fnc, type_fnc, ind_fnc, data_fnc, nb_data_fnc, id_fnc,
Computes time evolution function for seismic source.