SPEED
GET_FUNC_VALUE_SISM.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
19
34
35 real*8 function get_func_value_sism(nb_fnc, type_fnc, ind_fnc,&
36 data_fnc,nb_data_fnc, id_fnc, time, t0_delay, tau_new)
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
57 get_func_value_sism = 0.0d0
58 pi = 4.0d0 * datan(1.0d0)
59
60 select case (type_fnc(id_fnc))
61
62 case(0)
63 get_func_value_sism = 1.d0
64
65 case(1) ! - Ricker "beta" type
66 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
67 get_func_value_sism = (1.0d0 - 2.0d0*data_fnc(ind_fnc(id_fnc))*t_t0*t_t0) &
68 * dexp(-1.0d0*data_fnc(ind_fnc(id_fnc))*t_t0*t_t0)
69
70 case(2)
71 pi = 4.0d0 * datan(1.0d0)
72 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
73 get_func_value_sism = dcos(pi*data_fnc(ind_fnc(id_fnc))*t_t0) &
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)) &
84 get_func_value_sism = (v1 - v0) / (t1 - t0) * (time - t0) + v0
85 enddo
86
87 case(4) ! - First derivative of the Ricker ( d(Ricker(t))/dt )
88 t_t0 = time - data_fnc(ind_fnc(id_fnc) +1) - t0_delay
89 beta2=data_fnc(ind_fnc(id_fnc))
90 get_func_value_sism = 2.0d0*beta2*t_t0 &
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
97 get_func_value_sism = (2.d0*pi*data_fnc(ind_fnc(id_fnc))*t_t0) &
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)
100 write(33,*)time,get_func_value_sism
101
102 case(13) ! - GRENOBLE BENCHMARK
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
109 get_func_value_sism = 0.5d0 * ( 1.0d0 + erf(scaling*(t_t0)/hdur) )
110
111 case(14) ! - SCEC BENCHMARK
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 !val = amp * (1 - (1 + t_t0/TAU)*exp(-t_t0/TAU))
118 if (t_t0 .lt. 0.0d0) then
119 get_func_value_sism = 0.0d0
120 !elseif (t_t0.eq.0.0d0) then
121 ! GET_FUNC_VALUE_SISM = 0.5d0
122 else
123 get_func_value_sism = amp * (1.0d0 - (1 + t_t0/tau)*exp(-t_t0/tau))
124 endif
125
126
127
128 case(31) ! Input Source Time function using Text File
129! Supplied Input Time series file. Sufficient to just define the Unit Value Function
130! Function Value = Value_from_time_Series*Moment_tensor_component
131! Starting Time in Time series is 0sec. End Time is Rise Time
132! At the end of the Rise time, Value_from_time_Series = 1; (Unit Slip)
133 if (time.lt.t0_delay) then
134 get_func_value_sism = 0.d0;
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
142 get_func_value_sism = (v1 - v0) / (t1 - t0) * (time - t0) + v0
143 endif
144 enddo
145 elseif (time.gt.(t0_delay+tau_new)) then
146 get_func_value_sism = 1.d0;
147 endif
148
149 case(32) ! Ramp Slip-Time Function
150 if (time.lt.t0_delay) then
151 get_func_value_sism = 0.d0;
152 elseif ((time.ge.t0_delay).and.(time.le.(t0_delay+tau_new))) then
153 get_func_value_sism = (time-t0_delay)/tau_new;
154 elseif (time.gt.(t0_delay+tau_new)) then
155 get_func_value_sism = 1.d0;
156 endif
157
158
159 case(33) ! Input Source Time function using Text File
160! Supplied Input Time series file. Sufficient to just define the Unit Value Function
161! Function Value = Value_from_time_Series*Moment_tensor_component
162! Starting Time in Time series is 0sec. End Time is Rise Time
163! At the end of the Rise time, Value_from_time_Series = 1; (Unit Slip)
164 !write(*,*) time, t0_delay, tau_new
165 !read(*,*)
166
167 if (time.lt.t0_delay) then
168 get_func_value_sism = 0.d0;
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 ! write(*,*) 'qui'
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 ! write(*,*) t0,t1,v0,v1
177 ! read(*,*)
178
179 if ((time.ge.t0).and.(time.le.t1)) then
180 get_func_value_sism = (v1 - v0) / (t1 - t0) * (time - t0) + v0
181 endif
182 enddo
183 elseif (time.gt.(t0_delay+tau_new)) then
184 get_func_value_sism = 1.d0;
185 endif
186
187
188
189 case(50) ! - GRENOBLE BENCHMARK - NEW TAU
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
196 get_func_value_sism = 0.5d0 * ( 1.0d0 + erf(scaling*(t_t0)/hdur) )
197
198
199 case(55) ! - ARCHULETA FUNCTION
200 ! time = current time
201 ! t0_delay = rupture time
202 ! tau_new = rise time
203 !if ( time .lt. t0_delay + tau_new .and. time .gt. t0_delay) then
204 ! write(*,*) time, t0_delay, t0_delay + tau_new
205 !endif
206
207 if (time .le. t0_delay) then
208 get_func_value_sism = 0.d0;
209 elseif ( time .ge. t0_delay + tau_new) then
210 get_func_value_sism = 1.d0;
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 !do i = 1, np2
253 ! svf(i) = svf(i)/sum1
254 !enddo
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) = 0;
259
260 do i = 2, np2
261 int_svf(i) = 0.5*(svf(i)+svf(i-1))*dt;
262 !enddo
263 !do i = 2, np2
264 integ_svf(i) = sum(int_svf(2:i));
265 enddo
266
267 get_func_value_sism = integ_svf(np_current)
268 if (get_func_value_sism .gt. 1.05d0) then
269 write(*,*) 'Error in Archuleta slip function '
270 write(*,*) get_func_value_sism
271 stop
272 endif
273 deallocate(svf,int_svf,integ_svf)
274 !write(*,*) GET_FUNC_VALUE_SISM
275 !read(*,*)
276
277 endif
278
279
280
281 case(99) ! - CASHIMA 1
282 tau = data_fnc(ind_fnc(id_fnc))
283 amp = data_fnc(ind_fnc(id_fnc) +1)
284 t_t0 = time
285 get_func_value_sism = ( exp( - (((2.0d0 * 2.0d0*dasin(1.0d0))*1.5d0) &
286 *(t_t0 - tau)/amp)**2.0d0) * cos(((2.0d0 * 2.0d0*dasin(1.0d0))*1.5d0) &
287 *(t_t0 - tau) + 2.0d0*dasin(1.0d0)/2.0d0))
288
289 case default
290 get_func_value_sism = 0.d0
291
292 end select
293
294
295 return
296
297 end function get_func_value_sism
298
get_func_value_sism
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.
Definition GET_FUNC_VALUE_SISM.f90:37