SPEED
CHECK_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
45
46 subroutine check_sism(cs_nnz,cs,&
47 nm,tm,sd,&
48 nl_sism,&
49 num_ns,max_num_ns,&
50 sour_ns,dist_sour_ns,&
51 pos_sour_nx, pos_sour_ny, pos_sour_nz, &
52 check_ns,check_dist_ns,&
53 length_cns,&
54 check_pos_ns, &
55 fun_sism,nf,tag_func,val_sism,&
56 nn_loc, local_n_num,srcmodflag,szsism)
57
58
59
60! use speed_par, only: slip_type
61
62 implicit none
63
64 integer*4 :: cs_nnz,nm,ne,nl_sism,nf, nn_loc
65 integer*4 :: max_num_ns
66 integer*4 :: im,ie,isism,nn
67 integer*4 :: is,in,ip
68 integer*4 :: i,j,k,h
69 integer*4 :: fn
70 integer*4 :: length_cns,srcmodflag,szsism,val_st
71
72 integer*4, dimension(0:cs_nnz) :: cs
73 integer*4, dimension(nm) :: tm,sd
74 integer*4, dimension(nl_sism) :: num_ns
75 integer*4, dimension(nl_sism) :: fun_sism
76 integer*4, dimension(nn_loc) :: local_n_num
77 integer*4, dimension(nf) :: tag_func
78
79 integer*4, dimension(max_num_ns,nl_sism) :: sour_ns
80 integer*4, dimension(length_cns,4) :: check_ns
81
82 real*8 :: vel_prop, dist_b, xb, yb, zb, trup_b, dist_p, trup_p
83 real*8 :: dumvr, delay_tr
84
85 real*8, dimension(:), allocatable :: ct,ww
86
87 real*8, dimension(:,:), allocatable :: dd
88 real*8, dimension(max_num_ns,nl_sism) :: dist_sour_ns
89 real*8, dimension(max_num_ns,nl_sism) :: pos_sour_nx, pos_sour_ny, pos_sour_nz
90 real*8, dimension(length_cns,1) :: check_dist_ns
91 real*8, dimension(length_cns,3) :: check_pos_ns
92 real*8, dimension(nl_sism,szsism) :: val_sism
93
94
95
96
97 h = 0
98 nn = 2
99 ne = cs(0) -1
100
101 if (srcmodflag.eq.0) then
102 val_st = 12
103 elseif (srcmodflag.eq.1) then
104 val_st = 6
105 endif
106
107 allocate(ct(nn),ww(nn),dd(nn,nn))
108 call make_lgl_nw(nn,ct,ww,dd)
109
110 do im = 1,nm
111 if ((sd(im) +1).ne.nn) then
112 deallocate(ct,ww,dd)
113 nn = sd(im) +1
114 allocate(ct(nn),ww(nn),dd(nn,nn))
115 call make_lgl_nw(nn,ct,ww,dd)
116 endif
117
118 do ie = 1,ne
119 if (cs(cs(ie -1) +0).eq.tm(im)) then
120 do k = 1,nn
121 do j = 1,nn
122 do i = 1,nn
123 is = nn*nn*(k -1) +nn*(j -1) +i
124 in = cs(cs(ie -1) +is)
125
126 do isism = 1,nl_sism
127 do ip = 1,num_ns(isism)
128
129 if (local_n_num(in) .eq. sour_ns(ip,isism)) then
130 h = h + 1
131 check_ns(h,1) = sour_ns(ip,isism) !source node
132 check_ns(h,2) = fun_sism(isism) !fun type
133 check_ns(h,3) = isism !fault number
134 check_ns(h,4) = ie !local element
135 ! new --- added position of nodes
136 check_pos_ns(h,1) = pos_sour_nx(ip, isism)
137 check_pos_ns(h,2) = pos_sour_ny(ip, isism)
138 check_pos_ns(h,3) = pos_sour_nz(ip, isism)
139
140 !!distance from hypo / rupture velocity = rupture time --> std
141 !!rupture time ---> Archuleta (line 120 ~ rupt. velocity)
142 !!write(*,*) slip_type
143 !!read(*,*)
144 !!if (slip_type .eq. 'STD') then
145 ! !check_dist_ns(h,1) = dist_sour_ns(ip,isism) / val_sism(isism,19)
146 ! !old version val_sism(isism,19) = vrup
147 ! !new_version val_sism(isism,19) = trup
148 check_dist_ns(h,1) = val_sism(isism,val_st+7) ! Trupt for the baricenter of the triangle
149 !!elseif (slip_type .eq. 'ARC') then
150 !! check_dist_ns(h,1) = val_sism(isism,19)
151 !!elseif (slip_type .eq. 'GAL') then
152 !! check_dist_ns(h,1) = val_sism(isism,19)
153 !!endif
154 !! New implementation for trup (node by node)
155 !trup_b : dist_b = trup_p : dist_p --> trup_p = (trup_b * dist_p)/dist_b
156 if (srcmodflag.eq.0) then
157 xb = (val_sism(isism,4) + val_sism(isism,7) + val_sism(isism,10))/3.d0;
158 yb = (val_sism(isism,5) + val_sism(isism,8) + val_sism(isism,11))/3.d0;
159 zb = (val_sism(isism,6) + val_sism(isism,9) + val_sism(isism,12))/3.d0;
160
161 trup_b = val_sism(isism,19)
162 dist_b = dsqrt((val_sism(isism,1)-xb)**2.d0+(val_sism(isism,2)-yb)**2.d0+(val_sism(isism,3)-zb)**2.d0);
163 dist_p = dsqrt((val_sism(isism,1)-check_pos_ns(h,1))**2.d0 &
164 +(val_sism(isism,2)-check_pos_ns(h,2))**2.d0 &
165 +(val_sism(isism,3)-check_pos_ns(h,3))**2.d0);
166
167 trup_p = trup_b*dist_p/dist_b;
168
169 !write(*,*) h, dist_b, trup_b, dist_p, trup_p
170
171 if (dabs(dist_b) .le. 100) then
172 check_dist_ns(h,1) = val_sism(isism,19)
173 else
174 check_dist_ns(h,1) = (val_sism(isism,19) * dist_sour_ns(ip,isism)) / dist_b;
175 endif
176
177 elseif (srcmodflag.eq.1) then
178 xb = val_sism(isism,4)
179 yb = val_sism(isism,5)
180 zb = val_sism(isism,6)
181 dist_b = dsqrt((val_sism(isism,1)-xb)**2.d0+(val_sism(isism,2)-yb)**2.d0+(val_sism(isism,3)-zb)**2.d0)
182 dist_p = dsqrt((val_sism(isism,1)-check_pos_ns(h,1))**2.d0 &
183 +(val_sism(isism,2)-check_pos_ns(h,2))**2.d0 &
184 +(val_sism(isism,3)-check_pos_ns(h,3))**2.d0);
185
186 if (dabs(dist_b) .le. 100) then
187 check_dist_ns(h,1) = val_sism(isism,val_st+7);
188 else
189 ! delay_tr = Rupture time delay for respective segment + time window delay
190 !dumvr = 2500;
191 !delay_tr = val_sism(isism,val_st+7) - (dist_b/dumvr);
192 !check_dist_ns(h,1) = delay_tr + (dist_p/dumvr);
193 check_dist_ns(h,1) = val_sism(isism,val_st+7)*dist_p/dist_b;
194 endif
195 endif
196
197
198
199
200 endif
201
202 enddo !ip
203 enddo !isism
204
205 enddo !i
206 enddo !j
207 enddo !k
208
209
210 endif !if (cs....)
211
212 enddo !ie = 1,ne
213
214 enddo !im = 1,nm
215
216 return
217
218 end subroutine check_sism
219
220
check_sism
subroutine check_sism(cs_nnz, cs, nm, tm, sd, nl_sism, num_ns, max_num_ns, sour_ns, dist_sour_ns, pos_sour_nx, pos_sour_ny, pos_sour_nz, check_ns, check_dist_ns, length_cns, check_pos_ns, fun_sism, nf, tag_func, val_sism, nn_loc, local_n_num, srcmodflag, szsism)
Fills array check_ns for seismic force.
Definition CHECK_SISM.f90:57
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29