-
Notifications
You must be signed in to change notification settings - Fork 1
/
set_velocity_model.F90
144 lines (125 loc) · 7.2 KB
/
set_velocity_model.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
module set_velocity_model
!!gives velocity and attenuation structure for AmplitudeSourceLocation_PulseWidth.F90
!!Author: Masashi Ogiso (masashi.ogiso@gmail.com)
!!Copyright: (c) Masashi Ogiso 2020
!!License : MIT License https://opensource.org/licenses/MIT
private
public :: set_velocity
contains
subroutine set_velocity(z_min, dz, velocity, qinv)
use nrtype, only : fp
use constants, only : pi
implicit none
real(kind = fp), intent(in) :: z_min, dz !!size of depth direction
real(kind = fp), intent(out) :: velocity(:, :, :, :), qinv(:, :, :, :)
real(kind = fp), parameter :: VpVs = sqrt(3.0_fp)
integer :: i, nlon, nlat, nz
real(kind = fp) :: depth
nlon = ubound(velocity, 1)
nlat = ubound(velocity, 2)
nz = ubound(velocity, 3)
qinv(1 : nlon, 1 : nlat, 1 : nz, 1 : 2) = 0.0_fp
velocity(1 : nlon, 1 : nlat, 1 : nz, 1 : 2) = 0.0_fp
do i = 1, nz
depth = z_min + dz * real(i - 1, kind = fp)
#if defined (V_MEA1D)
if(i .eq. 1) write(0, '(a)') "velocity/Qinv = meakan 1d structure"
if(depth .lt. 0.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (3.2_fp - 3.0_fp) / (0.5_fp - 0.0_fp) * depth + 3.0_fp
elseif(depth .ge. 0.5_fp .and. depth .lt. 1.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (4.5_fp - 3.2_fp) / (1.0_fp - 0.5_fp) * (depth - 0.5_fp) + 3.2_fp
elseif(depth .ge. 1.0_fp .and. depth .lt. 1.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.0_fp - 4.5_fp) / (1.5_fp - 1.0_fp) * (depth - 1.0_fp) + 4.5_fp
elseif(depth .ge. 1.5_fp .and. depth .lt. 2.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.133_fp - 5.0_fp) / (2.0_fp - 1.5_fp) * (depth - 1.5_fp) + 5.0_fp
elseif(depth .ge. 2.0_fp .and. depth .lt. 2.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.262_fp - 5.133_fp) / (2.5_fp - 2.0_fp) * (depth - 2.0_fp) + 5.133_fp
elseif(depth .ge. 2.5_fp .and. depth .lt. 3.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.383_fp - 5.262_fp) / (3.0_fp - 2.5_fp) * (depth - 2.5_fp) + 5.262_fp
elseif(depth .ge. 3.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.494_fp - 5.383_fp) / (3.5_fp - 3.0_fp) * (depth - 3.0_fp) + 5.383_fp
endif
!!Vp to Vs
velocity(1 : nlon, 1 : nlat, i, 2) = velocity(1 : nlon, 1 : nlat, i, 1) / VpVs
!!the values of Qinv are taken from Kumagai et al (2019): case of Nevado
!del Ruiz volcano (layer thickness is different)
if(depth .lt. 1.0_fp) then
qinv(1 : nlon, 1 : nlat, i, 2) = 1.0_fp / 40.0_fp
else
qinv(1 : nlon, 1 : nlat, i, 2) = 1.0_fp / 180.0_fp
endif
qinv(1 : nlon, 1 : nlat, i, 1) = qinv(1 : nlon, 1 : nlat, i, 2) / (9.0_fp / 4.0_fp)
#elif defined (V_CONST)
if(i .eq. 1) write(0, '(a)') "Velocity/Qinv = constant"
velocity(1 : nlon, 1 : nlat, i, 1) = 2.5_fp
!!Vp to Vs
velocity(1 : nlon, 1 : nlat, i, 2) = velocity(1 : nlon, 1 : nlat, i, 1) / VpVs
qinv(1 : nlon, 1 : nlat, i, 2) = 1.0_fp / 500.0_fp
qinv(1 : nlon, 1 : nlat, i, 1) = qinv(1 : nlon, 1 : nlat, i, 2) / (9.0_fp / 4.0_fp)
#elif defined (V_OFFKII)
!!velocity structure read visually from Nakanishi et al., 2002, JGR, doi: 10.1029/2001JB000424
if(i .eq. 1) write(0, '(a)') "velocity = off Kii peninsula"
if(depth .le. 7.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (4.0_fp - 2.0_fp) / (7.0_fp - 0.0_fp) * (depth - 0.0_fp) + 2.0_fp
elseif(depth .gt. 7.0_fp .and. depth .le. 9.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (5.6_fp - 4.0_fp) / (9.0_fp - 7.0_fp) * (depth - 7.0_fp) + 4.0_fp
elseif(depth .gt. 9.0_fp .and. depth .le. 13.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (6.8_fp - 5.6_fp) / (13.0_fp - 9.0_fp) * (depth - 9.0_fp) + 5.6_fp
elseif(depth .gt. 13.0_fp .and. depth .le. 14.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (8.0_fp - 6.8_fp) / (14.0_fp - 13.0_fp) * (depth - 13.0_fp) + 6.8_fp
elseif(depth .gt. 14.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = (8.2_fp - 8.0_fp) / (33.0_fp - 14.0_fp) * (depth - 14.0_fp) + 8.0_fp
!elseif(depth .gt. 7.0_fp .and. depth .le. 11.0_fp) then
! velocity(1 : nlon, 1 : nlat, i, 1) = (8.0_fp - 2.0_fp) / (11.0_fp - 7.0_fp) * (depth - 7.0_fp) + 2.0_fp
!elseif(depth .gt. 11.0_fp) then
! velocity(1 : nlon, 1 : nlat, i, 1) = (8.2_fp - 8.0_fp) / (45.0_fp - 11.0_fp) * (depth - 11.0_fp) + 8.0_fp
endif
!!Vp to Vs
velocity(1 : nlon, 1 : nlat, i, 2) = velocity(1 : nlon, 1 : nlat, i, 1) / VpVs
!!freq = 5.0 (Hz, 2-8 Hz BPF), velocity = 3.5 km/s, C = 0.02 (1/km) from
!!Yabe et al., 2020, Techtonophysics, doi: 10.1016/j.tecto.2020.228714
!!modify for depth > 14 km
qinv(1 : nlon, 1 : nlat, i, 2) = (0.02_fp * 3.5_fp) / (pi * 5.0_fp)
qinv(1 : nlon, 1 : nlat, i, 1) = qinv(1 : nlon, 1 : nlat, i, 2) / (9.0_fp / 4.0_fp)
if(depth .gt. 14.0_fp) qinv(1 : nlon, 1 : nlat, i, 1 : 2) = qinv(1 : nlon, 1 : nlat, i, 1 : 2) * 0.5_fp
#elif defined (JMA2001)
if(i .eq. 1) write(0, '(a)') "velocity = (approx.) JMA2001"
if(depth .le. 4.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = ((5.780_fp - 4.800_fp) / 4.5_fp) * depth + 4.800_fp
velocity(1 : nlon, 1 : nlat, i, 2) = ((3.409_fp - 2.844_fp) / 4.5_fp) * depth + 2.844_fp
elseif(depth .gt. 4.5_fp .and. depth .le. 21.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = &
& ((6.502_fp - 5.780_fp) / (21.0_fp - 4.5_fp)) * (depth - 4.5_fp) + 5.780_fp
velocity(1 : nlon, 1 : nlat, i, 2) = &
& ((3.782_fp - 3.409_fp) / (21.0_fp - 4.5_fp)) * (depth - 4.5_fp) + 3.409_fp
elseif(depth .gt. 21.0_fp .and. depth .le. 45.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = &
& ((7.638_fp - 6.502_fp) / (45.5_fp - 21.0_fp)) * (depth - 21.0_fp) + 6.502_fp
velocity(1 : nlon, 1 : nlat, i, 2) = &
& ((4.362_fp - 3.782_fp) / (45.5_fp - 21.0_fp)) * (depth - 21.0_fp) + 3.782_fp
elseif(depth .gt. 45.5_fp .and. depth .le. 83.0_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = &
& ((7.900_fp - 7.638_fp) / (83.0_fp - 45.5_fp)) * (depth - 45.5_fp) + 7.638_fp
velocity(1 : nlon, 1 : nlat, i, 2) = &
& ((4.424_fp - 4.362_fp) / (83.0_fp - 45.5_fp)) * (depth - 45.5_fp) + 4.362_fp
elseif(depth .gt. 83.0_fp .and. depth .le. 175.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = &
& ((8.184_fp - 7.900_fp) / (175.5_fp - 83.0_fp)) * (depth - 83.0_fp) + 7.900_fp
velocity(1 : nlon, 1 : nlat, i, 2) = &
& ((4.558_fp - 4.424_fp) / (175.5_fp - 83.0_fp)) * (depth - 83.0_fp) + 4.424_fp
elseif(depth .gt. 175.5_fp) then
velocity(1 : nlon, 1 : nlat, i, 1) = &
& ((8.269_fp - 8.184_fp) / (200.0_fp - 175.5_fp)) * (depth - 175.5_fp) + 8.184_fp
velocity(1 : nlon, 1 : nlat, i, 2) = &
& ((4.602_fp - 4.558_fp) / (200.0_fp - 175.5_fp)) * (depth - 175.5_fp) + 4.558_fp
endif
qinv(1 : nlon, 1 : nlat, i, 2) = 1.0_fp / 200.0_fp
qinv(1 : nlon, 1 : nlat, i, 1) = qinv(1 : nlon, 1 : nlat, i, 2) / (9.0_fp / 4.0_fp)
#else
write(0, *) "please set -DV_CONST or appropriate definition"
stop
#endif
enddo
return
end subroutine set_velocity
end module set_velocity_model