SPEED
MAKE_ANELASTIC_COEFFICIENTS_NH.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
36
37 subroutine make_anelastic_coefficients_nh(nn, N_SLS, lambda_el, mu_el,&
38 rho_el, QS_nh, QP_nh, f_val, &
39 Y_lambda_nh, Y_mu_nh, frequency_range, mpi_id)
40
41 use speed_exit_codes
42
43 implicit none
44
45 integer*4 :: nn, mpi_id, N_SLS, i, j, k, &
46 PIVOT_QP(N_SLS),PIVOT_QS(N_SLS)
47
48 real*8 :: f_val, fmin, fmax, esp1, esp2, deltax, &
49 vp2, vs2, pi
50
51 real*8 :: qs_nh, qp_nh, &
52 y_lambda_nh(1,n_sls),&
53 y_mu_nh(1,n_sls), &
54 frequency_range(n_sls),&
55 frequency_range_sampling(2*n_sls-1),&
56 ys(n_sls),yp(n_sls), &
57 a_qp(2*n_sls-1,n_sls),rhs_qp(2*n_sls-1),a_qs(2*n_sls-1,n_sls),rhs_qs(2*n_sls-1),&
58 rho, rho_el(nn,nn,nn), &
59 lambda, lambda_el(nn,nn,nn), &
60 mu, mu_el(nn,nn,nn), &
61 rho_all, mu_all, lambda_all
62
63 rho_all = 0.d0; lambda_all = 0.d0; mu_all = 0.d0
64
65 do i =1, nn
66 do j =1, nn
67 do k =1, nn
68 mu_all = mu_all + mu_el(i,j,k)
69 lambda_all = lambda_all + lambda_el(i,j,k)
70 rho_all = rho_all + rho_el(i,j,k)
71 enddo
72 enddo
73 enddo
74
75 mu = mu_all/(nn**3)
76 lambda = lambda_all/(nn**3)
77 rho = rho_all/(nn**3)
78
79 pi = 4.d0*datan(1.d0);
80
81 if (f_val .le. 10) then
82 fmin = 0.05d0; !fmax = 5.d0;
83 elseif(f_val .le. 100) then
84 fmin = 5.d0; !fmax = 100.d0;
85 elseif(f_val .le. 1000) then
86 fmin = 50.d0; !fmax = 1000.d0;
87 elseif(f_val .le. 10000) then
88 fmin = 500.d0; !fmax = 10000.d0;
89 elseif(f_val .le. 100000) then
90 fmin = 5000.d0; !fmax = 100000.d0;
91 endif
92
93 f_val = 1;
94
95 fmax=100*fmin;
96
97 !if(mpi_id .eq. 0) then
98 ! write(*,*) 'Frequency range where Quality Factor is assumed to be constant is:'
99 ! write(*,*) 'FMIN =' ,fmin, ' FMAX =', fmax;
100 !endif
101
102 esp1 = log10(fmin); esp2 = log10(fmax);
103 deltax = (esp2-esp1)/(n_sls-1);
104
105
106 if (n_sls .eq. 3) then
107
108 frequency_range(1) = f_val*0.1 !2.*pi*fmin;
109 frequency_range(2) = f_val*1 !*fmin;
110 frequency_range(3) = f_val*10 !*fmax;
111 frequency_range_sampling(1) = frequency_range(1);
112 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
113 frequency_range_sampling(3) = frequency_range(2);
114 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
115 frequency_range_sampling(5) = frequency_range(3);
116
117
118 elseif (n_sls .eq. 4) then
119
120 frequency_range(1) = 2.*pi*fmin;
121 frequency_range(2) = 2.*pi*10**(esp1+deltax)
122 frequency_range(3) = 2.*pi*10**(esp1+2*deltax)
123 frequency_range(4) = 2.*pi*fmax;
124 frequency_range_sampling(1) = frequency_range(1);
125 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
126 frequency_range_sampling(3) = frequency_range(2);
127 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
128 frequency_range_sampling(5) = frequency_range(3);
129 frequency_range_sampling(6) = 0.5*(frequency_range(3)+frequency_range(4))
130 frequency_range_sampling(7) = frequency_range(4);
131
132
133 elseif (n_sls .eq. 5) then
134
135 frequency_range(1) = f_val*0.1 ! 2.*pi*fmin;
136 frequency_range(2) = f_val*0.5 !2.*pi*10**(esp1+deltax)
137 frequency_range(3) = f_val*1 !2.*pi*10**(esp1+2*deltax)
138 frequency_range(4) = f_val*5 !2.*pi*10**(esp1+3*deltax)
139 frequency_range(5) = f_val*10 !2.*pi*fmax;
140 frequency_range_sampling(1) = frequency_range(1);
141 frequency_range_sampling(2) = 0.5*(frequency_range(1)+frequency_range(2))
142 frequency_range_sampling(3) = frequency_range(2);
143 frequency_range_sampling(4) = 0.5*(frequency_range(2)+frequency_range(3))
144 frequency_range_sampling(5) = frequency_range(3);
145 frequency_range_sampling(6) = 0.5*(frequency_range(3)+frequency_range(4))
146 frequency_range_sampling(7) = frequency_range(4);
147 frequency_range_sampling(6) = 0.5*(frequency_range(4)+frequency_range(5))
148 frequency_range_sampling(7) = frequency_range(5);
149
150
151 endif
152
153
154
155
156
157 vp2 = (lambda + 2.d0*mu)/rho; vs2 = mu/rho;
158
159 y_mu_nh(1,:) = 0.d0; y_lambda_nh(1,:) = 0.d0;
160 if (qp_nh .ne. 0.d0 .and. qs_nh .ne. 0.d0) then
161
162 do i = 1, 2*n_sls-1
163 do j = 1, n_sls
164
165 a_qp(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qp_nh)) &
166 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
167
168 a_qs(i,j) = (frequency_range(j)*frequency_range_sampling(i) + frequency_range(j)**2*(1.d0/qs_nh)) &
169 / (frequency_range(j)**2 + frequency_range_sampling(i)**2);
170
171 enddo
172
173 rhs_qp(i) = 1.d0/qp_nh;
174 rhs_qs(i) = 1.d0/qs_nh;
175
176 enddo
177
178 !call FACTORIZE_MATRIX(A_QP,N_SLS,PIVOT_QP)
179 !call FACTORIZE_MATRIX(A_QS,N_SLS,PIVOT_QS)
180
181 !call DIRECT_LU_SOLVER(A_QP,RHS_QP,N_SLS,YP,PIVOT_QP)
182 !call DIRECT_LU_SOLVER(A_QS,RHS_QS,N_SLS,YS,PIVOT_QS)
183
184 call qr_solve(2*n_sls-1,n_sls, a_qp, rhs_qp, yp)
185 call qr_solve(2*n_sls-1,n_sls, a_qs, rhs_qs, ys)
186
187
188 y_mu_nh(1,:) = ys;
189 y_lambda_nh(1,:) = (vp2*yp - 2.d0*vs2*ys)/(vp2-2.d0*vs2);
190
191 endif
192
193
194! do i = 1, N_SLS
195! if (Y_lambda_nh(1,i) .lt. 0.d0 ) then
196! write(*,*) 'ERROR! ANELASTIC COEFFICIENTS HAVE TO BE POSITIVE!'
197! write(*,*) 'Y_LAMBDA = ', Y_lambda_nh(1,i)
198! write(*,*) 'SETTING Y_LAMBDA = 0'
199! Y_lambda_nh(1,i) = 0.d0
200! endif
201! if (Y_mu_nh(1,i) .lt. 0.d0) then
202! write(*,*) 'ERROR! ANELASTIC COEFFICIENTS HAVE TO BE POSITIVE!'
203! write(*,*) 'Y_MU = ', Y_mu_nh(1,i)
204! write(*,*) 'SETTING Y_MU = 0'
205! Y_mu_nh(1,i) = 0.d0
206! write(*,*) '----------------------------------------------------'
207! endif
208! enddo
209
210
211
212 return
213
214 end subroutine make_anelastic_coefficients_nh
215
216
217
make_anelastic_coefficients_nh
subroutine make_anelastic_coefficients_nh(nn, n_sls, lambda_el, mu
Compute anelastic coefficients for Standard Linear Solid model, when the not-honoring strategy is app...
Definition MAKE_ANELASTIC_COEFFICIENTS_NH.f90:38
qr_solve
subroutine qr_solve(m, n, a, b, x)
Definition QR_SOLVE.f90:2384
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25