forked from alistairboyce11/RJ_MCMC
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathpostprocess_binary_outputs.f90
More file actions
177 lines (119 loc) · 5 KB
/
postprocess_binary_outputs.f90
File metadata and controls
177 lines (119 loc) · 5 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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
program postprocess_binary_outputs
use mpi
implicit none
include 'params.h'
include 'data_joint.h'
real widening_prop
real widening_prop2
integer i_w,i,io_file,io
character filenamemax*300,filenamemax2*300
integer numsample,everyall_2,burn_in_2,thin_2
real d_min_2,d_max_2,width_2,xi_min_2,xi_max_2,vp_min_2,vp_max_2
DOUBLE PRECISION Ad_R_min_2,Ad_R_max_2,Ad_L_min_2,Ad_L_max_2
integer milay_2,malay_2,mk_2,ndatadmax_2
integer nptfinal,npt,npt_ani,nic,noc
DOUBLE PRECISION Ad_R,Ad_L
real, DIMENSION(mk) :: r,vsv,xi,vpv
real like_w
integer ndatad_R,ndatad_L
real, DIMENSION(ndatadmax) :: d_cR,d_cL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Read binary files containing models and converts them into ascii for python to read
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
widening_prop2=1.
do i_w=1,10
!write(*,*)widening_prop2
do i=0,200
write(filenamemax,"('/All_models_prepare_',I4.4,'_',f5.2,'.out')")i,widening_prop2
!write(*,*)dirname//filenamemax
open(100,file=dirname//filenamemax,status='old',form='unformatted',access='stream',iostat=io_file)
if (io_file/=0) goto 500
write(filenamemax2,"('/All_models_processed_prepare_',I4.4,'_',f5.2,'.out')")i,widening_prop2
write(*,*)filenamemax2
open(200,file=dirname//filenamemax2,status='replace')
read(100,IOSTAT=io)numsample
read(100,IOSTAT=io)everyall_2
write(200,*)numsample,everyall_2
!write(*,*)numsample,everyall_2
read(100,IOSTAT=io)burn_in_2
read(100,IOSTAT=io)widening_prop
read(100,IOSTAT=io)thin_2
write(200,*)burn_in_2,widening_prop,thin_2
read(100,IOSTAT=io)d_min_2
read(100,IOSTAT=io)d_max_2
write(200,*)d_min_2,d_max_2
read(100,IOSTAT=io)width_2
write(200,*)width_2
read(100,IOSTAT=io)xi_min_2
read(100,IOSTAT=io)xi_max_2
write(200,*)xi_min_2,xi_max_2
read(100,IOSTAT=io)vp_min_2
read(100,IOSTAT=io)vp_max_2
write(200,*)vp_min_2,vp_max_2
read(100,IOSTAT=io)Ad_R_min_2
read(100,IOSTAT=io)Ad_R_max_2
write(200,*)Ad_R_min_2,Ad_R_max_2
read(100,IOSTAT=io)Ad_L_min_2
read(100,IOSTAT=io)Ad_L_max_2
write(200,*)Ad_L_min_2,Ad_L_max_2
read(100,IOSTAT=io)milay_2
read(100,IOSTAT=io)malay_2
write(200,*)milay_2,malay_2
read(100,IOSTAT=io)mk_2
read(100,IOSTAT=io)ndatadmax_2
io=0
do while (io==0)
read(100,IOSTAT=io)nptfinal
!write(*,*)nptfinal
if (io/=0) goto 500
if (nptfinal>mk) CONTINUE
read(100,IOSTAT=io)nic
if (io/=0) goto 500
read(100,IOSTAT=io)noc
if (io/=0) goto 500
read(100,IOSTAT=io)npt
if (io/=0) goto 500
read(100,IOSTAT=io)npt_ani
if (io/=0) goto 500
write(200,*)nptfinal-noc,npt,npt_ani
read(100,IOSTAT=io)Ad_R
if (io/=0) goto 500
read(100,IOSTAT=io)Ad_L
if (io/=0) goto 500
write(200,*)Ad_R,Ad_L
read(100,IOSTAT=io)r
if (io/=0) goto 500
write(200,*)(rearth-r(nptfinal:noc+1:-1))/1000. ! only take what is needed for postprocessing
read(100,IOSTAT=io)vsv
if (io/=0) goto 500
write(200,*)vsv(nptfinal:noc+1:-1)/1000.
read(100,IOSTAT=io)xi
if (io/=0) goto 500
write(200,*)xi(nptfinal:noc+1:-1)
read(100,IOSTAT=io)vpv ! careful! We store vpv, so we need to check it for the python postprocessing
if (io/=0) goto 500
write(200,*)vpv(nptfinal:noc+1:-1)/1000.
read(100,IOSTAT=io)like_w
if (io/=0) goto 500
read(100,IOSTAT=io)ndatad_R
if (io/=0) goto 500
read(100,IOSTAT=io)d_cR
if (io/=0) goto 500
read(100,IOSTAT=io)ndatad_L
if (io/=0) goto 500
read(100,IOSTAT=io)d_cL
if (io/=0) goto 500
! because it wasn't put in the initial like_w
like_w=-like_w-ndatad_R*log(Ad_R)/widening_prop-ndatad_L*log(Ad_L)/widening_prop
write(200,*)like_w
write(200,*)ndatad_R
write(200,*)d_cR(:ndatad_R)
write(200,*)ndatad_L
write(200,*)d_cL(:ndatad_L)
enddo
600 close(100)
500 close(200)
enddo
widening_prop2=widening_prop2+1
enddo
end program postprocess_binary_outputs