Repository navigation
Expand file tree
/
Copy pathdepGridDws_FDD09.f90
More file actions
146 lines (136 loc) · 3.55 KB
/
Copy pathdepGridDws_FDD09.f90
File metadata and controls
146 lines (136 loc) · 3.55 KB
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
142
143
144
145
146
!Output absolute velocity and relative perturbation velocity per layer
!3 Inputfiles : grid file: grid.dat; initial mod file: MOD; velocity file: Vp_model.dat
program createDepthGrid
implicit none
integer :: i,j,k
real :: bld
integer :: nx,ny,nz,nxyz2,cont
real,allocatable :: xn(:),yn(:),zn(:)
real,allocatable :: modvel(:,:,:)
real,allocatable :: grid(:,:)
real,allocatable :: absvel(:,:,:) !absolute velocity
real,allocatable :: relatvel(:,:,:) !relative velocity perturbation
real,allocatable :: dwsvalue(:,:,:)
integer :: stat=0
integer :: zdepth
character(len=2) :: zdepthstr,nxstr
character(len=15) :: absfie
character(len=15) :: perturfie
character(len=15) :: dwsfile
character(len=20) :: line
character(len=30) :: formatread
open(1,file='MOD',status='old')
read(1,*) bld,nx,ny,nz
allocate(xn(nx))
allocate(yn(ny))
allocate(zn(nz))
read(1,*) xn
read(1,*) yn
read(1,*) zn
write(nxstr,'(I2)') nx
allocate(modvel(nx,ny,nz))
formatread='('//trim(nxstr)//'f6.2)'
Do k = 1, nz
Do j = 1, ny
read (1, trim(formatread),iostat=stat) (modvel(i,j,k), i=1,nx)
if(stat /= 0) exit
End Do
End Do
close(1)
nxyz2 = (nx-2)*(ny-2)*(nz-2)
allocate(grid(nxyz2,2))
open(2,file='grid.dat',status='old')
do i=1,nxyz2
read (2, *,iostat=stat) grid(i,1),grid(i,2)
if(stat /= 0) exit
end do
close(2)
allocate(absvel(nx,ny,nz))
formatread='('//trim(nxstr)//'f7.3)'
open(3,file='Vp_model.dat',status='old')
Do k = 1, nz
Do j = 1, ny
read (3, trim(formatread),iostat=stat) (absvel(i,j,k), i=1,nx)
if(stat /= 0) exit
End Do
End Do
close(3)
allocate(dwsvalue(nx,ny,nz))
formatread='('//trim(nxstr)//'f10.2)'
open(4,file='tomofdd.vel',status='old')
do while(.true.)
10 read(4,'(a)',iostat=stat) line
if(stat /= 0) exit
if (trim(line) == ' DWS for P-wave.') then
Do k = 1, nz
Do j = 1, ny
read (4, trim(formatread)) (dwsvalue(i,j,k), i=1,nx)
End Do
End Do
else
goto 10
end if
end do
close(4)
! Output absolute velocity file
do k=2,nz-1
if(zn(k) .GE. 0 .and. zn(k) .LT. 100) then
zdepth=INT(zn(k))
if (zdepth < 10) then
write(zdepthstr,'(I1)') zdepth
else
write(zdepthstr,'(I2)') zdepth
end if
absfie='absfie'//trim(zdepthstr)//'.dat'
open(10+k,file=absfie)
do j=2,ny-1
do i=2,nx-1
cont = i+(j-2)*(nx-2)+(k-2)*(nx-2)*(ny-2)-1
write(10+k,'(f7.3,2X,f7.3,2X,f7.3)') grid(cont,1),grid(cont,2),absvel(i,j,k)
end do
end do
close(10+k)
end if
end do
! Output relative perturbation velocity file
do k=2,nz-1
if(zn(k) .GE. 0 .and. zn(k) .LT. 100) then
zdepth=INT(zn(k))
if (zdepth < 10) then
write(zdepthstr,'(I1)') zdepth
else
write(zdepthstr,'(I2)') zdepth
end if
perturfie='perturfie'//trim(zdepthstr)//'.dat'
open(20+k,file=perturfie)
do j=2,ny-1
do i=2,nx-1
cont = i+(j-2)*(nx-2)+(k-2)*(nx-2)*(ny-2)-1
write(20+k,'(f7.3,2X,f7.3,2X,f6.2)') grid(cont,1),grid(cont,2),(absvel(i,j,k)-modvel(i,j,k))*100/modvel(i,j,k)
end do
end do
close(20+k)
end if
end do
! output DWS for every node from tomoDD.vel
do k=2,nz-1
if(zn(k) .GE. 0 .and. zn(k) .LT. 100) then
zdepth=INT(zn(k))
if (zdepth < 10) then
write(zdepthstr,'(I1)') zdepth
else
write(zdepthstr,'(I2)') zdepth
end if
dwsfile='dwsfile'//trim(zdepthstr)//'.dat'
open(30+k,file=dwsfile)
do j=2,ny-1
do i=2,nx-1
cont = i+(j-2)*(nx-2)+(k-2)*(nx-2)*(ny-2)-1
write(30+k,'(f7.3,2X,f7.3,2X,f10.2)') grid(cont,1),grid(cont,2),dwsvalue(i,j,k)
end do
end do
close(30+k)
end if
end do
stop
end program