-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathinidriver.F90
366 lines (301 loc) · 14 KB
/
inidriver.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
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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
! Code for Anisotropies in the Microwave Background
! by Antony Lewis (http://cosmologist.info/) and Anthony Challinor
! See readme.html for documentation. This is a sample driver routine that reads
! in one set of parameters and produdes the corresponding output.
program driver
use IniFile
use CAMB
use LambdaGeneral
use Lensing
use AMLUtils
use Transfer
use constants
use Bispectrum
use CAMBmain
use NonLinear
#ifdef NAGF95
use F90_UNIX
#endif
implicit none
Type(CAMBparams) P
character(LEN=Ini_max_string_len) numstr, VectorFileName, &
InputFile, ScalarFileName, TensorFileName, TotalFileName, LensedFileName,&
LensedTotFileName, LensPotentialFileName,ScalarCovFileName
integer i
character(LEN=Ini_max_string_len) TransferFileNames(max_transfer_redshifts), &
MatterPowerFileNames(max_transfer_redshifts), outroot, version_check
real(dl) output_factor, nmassive
#ifdef WRITE_FITS
character(LEN=Ini_max_string_len) FITSfilename
#endif
logical bad
InputFile = ''
if (GetParamCount() /= 0) InputFile = GetParam(1)
if (InputFile == '') error stop 'No parameter input file'
call Ini_Open(InputFile, 1, bad, .false.)
if (bad) error stop 'Error opening parameter file'
Ini_fail_on_not_found = .false.
outroot = Ini_Read_String('output_root')
if (outroot /= '') outroot = trim(outroot) // '_'
highL_unlensed_cl_template = Ini_Read_String_Default('highL_unlensed_cl_template',highL_unlensed_cl_template)
call CAMB_SetDefParams(P)
P%WantScalars = Ini_Read_Logical('get_scalar_cls')
P%WantVectors = Ini_Read_Logical('get_vector_cls',.false.)
P%WantTensors = Ini_Read_Logical('get_tensor_cls',.false.)
P%OutputNormalization=outNone
output_factor = Ini_Read_Double('CMB_outputscale',1.d0)
P%WantCls= P%WantScalars .or. P%WantTensors .or. P%WantVectors
P%PK_WantTransfer=Ini_Read_Logical('get_transfer')
AccuracyBoost = Ini_Read_Double('accuracy_boost',AccuracyBoost)
lAccuracyBoost = Ini_Read_Real('l_accuracy_boost',lAccuracyBoost)
HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',HighAccuracyDefault)
P%NonLinear = Ini_Read_Int('do_nonlinear',NonLinear_none)
P%DoLensing = .false.
if (P%WantCls) then
if (P%WantScalars .or. P%WantVectors) then
P%Max_l = Ini_Read_Int('l_max_scalar')
P%Max_eta_k = Ini_Read_Double('k_eta_max_scalar',P%Max_l*2._dl)
if (P%WantScalars) then
P%DoLensing = Ini_Read_Logical('do_lensing',.false.)
if (P%DoLensing) lensing_method = Ini_Read_Int('lensing_method',1)
end if
if (P%WantVectors) then
if (P%WantScalars .or. P%WantTensors) error stop 'Must generate vector modes on their own'
i = Ini_Read_Int('vector_mode')
if (i==0) then
vec_sig0 = 1
Magnetic = 0
else if (i==1) then
Magnetic = -1
vec_sig0 = 0
else
error stop 'vector_mode must be 0 (regular) or 1 (magnetic)'
end if
end if
end if
if (P%WantTensors) then
P%Max_l_tensor = Ini_Read_Int('l_max_tensor')
P%Max_eta_k_tensor = Ini_Read_Double('k_eta_max_tensor',Max(500._dl,P%Max_l_tensor*2._dl))
end if
endif
! Read initial parameters.
call DarkEnergy_ReadParams(DefIni)
P%h0 = Ini_Read_Double('hubble')
if (Ini_Read_Logical('use_physical',.false.)) then
P%omegab = Ini_Read_Double('ombh2')/(P%H0/100)**2
P%omegac = Ini_Read_Double('omch2')/(P%H0/100)**2
P%omegan = Ini_Read_Double('omnuh2')/(P%H0/100)**2
P%omegav = 1- Ini_Read_Double('omk') - P%omegab-P%omegac - P%omegan
else
P%omegab = Ini_Read_Double('omega_baryon')
P%omegac = Ini_Read_Double('omega_cdm')
P%omegav = Ini_Read_Double('omega_lambda')
P%omegan = Ini_Read_Double('omega_neutrino')
end if
P%tcmb = Ini_Read_Double('temp_cmb',COBE_CMBTemp)
P%yhe = Ini_Read_Double('helium_fraction',0.24_dl)
P%Num_Nu_massless = Ini_Read_Double('massless_neutrinos')
P%Nu_mass_eigenstates = Ini_Read_Int('nu_mass_eigenstates',1)
if (P%Nu_mass_eigenstates > max_nu) error stop 'too many mass eigenstates'
numstr = Ini_Read_String('massive_neutrinos')
read(numstr, *) nmassive
if (abs(nmassive-nint(nmassive))>1e-6) error stop 'massive_neutrinos should now be integer (or integer array)'
read(numstr,*, end=100, err=100) P%Nu_Mass_numbers(1:P%Nu_mass_eigenstates)
P%Num_Nu_massive = sum(P%Nu_Mass_numbers(1:P%Nu_mass_eigenstates))
if (P%Num_Nu_massive>0) then
P%share_delta_neff = Ini_Read_Logical('share_delta_neff', .true.)
numstr = Ini_Read_String('nu_mass_degeneracies')
if (P%share_delta_neff) then
if (numstr/='') write (*,*) 'WARNING: nu_mass_degeneracies ignored when share_delta_neff'
else
if (numstr=='') error stop 'must give degeneracies for each eigenstate if share_delta_neff=F'
read(numstr,*) P%Nu_mass_degeneracies(1:P%Nu_mass_eigenstates)
end if
numstr = Ini_Read_String('nu_mass_fractions')
if (numstr=='') then
if (P%Nu_mass_eigenstates >1) error stop 'must give nu_mass_fractions for the eigenstates'
P%Nu_mass_fractions(1)=1
else
read(numstr,*) P%Nu_mass_fractions(1:P%Nu_mass_eigenstates)
end if
end if
!JD 08/13 begin changes for nonlinear lensing of CMB + LSS compatibility
!P%Transfer%redshifts -> P%Transfer%PK_redshifts and P%Transfer%num_redshifts -> P%Transfer%PK_num_redshifts
!in the P%WantTransfer loop.
if (((P%NonLinear==NonLinear_lens .or. P%NonLinear==NonLinear_both) .and. P%DoLensing) &
.or. P%PK_WantTransfer) then
P%Transfer%high_precision= Ini_Read_Logical('transfer_high_precision',.false.)
else
P%transfer%high_precision = .false.
endif
if (P%NonLinear/=NonLinear_none) call NonLinear_ReadParams(DefIni)
if (P%PK_WantTransfer) then
P%WantTransfer = .true.
P%transfer%kmax = Ini_Read_Double('transfer_kmax')
P%transfer%k_per_logint = Ini_Read_Int('transfer_k_per_logint')
P%transfer%PK_num_redshifts = Ini_Read_Int('transfer_num_redshifts')
transfer_interp_matterpower = Ini_Read_Logical('transfer_interp_matterpower ', transfer_interp_matterpower)
transfer_power_var = Ini_read_int('transfer_power_var',transfer_power_var)
if (P%transfer%PK_num_redshifts > max_transfer_redshifts) error stop 'Too many redshifts'
do i=1, P%transfer%PK_num_redshifts
P%transfer%PK_redshifts(i) = Ini_Read_Double_Array('transfer_redshift',i,0._dl)
transferFileNames(i) = Ini_Read_String_Array('transfer_filename',i)
MatterPowerFilenames(i) = Ini_Read_String_Array('transfer_matterpower',i)
if (TransferFileNames(i) == '') then
TransferFileNames(i) = trim(numcat('transfer_',i))//'.dat'
end if
if (MatterPowerFilenames(i) == '') then
MatterPowerFilenames(i) = trim(numcat('matterpower_',i))//'.dat'
end if
if (TransferFileNames(i)/= '') &
TransferFileNames(i) = trim(outroot)//TransferFileNames(i)
if (MatterPowerFilenames(i) /= '') &
MatterPowerFilenames(i)=trim(outroot)//MatterPowerFilenames(i)
end do
else
P%Transfer%PK_num_redshifts = 1
P%Transfer%PK_redshifts = 0
end if
if ((P%NonLinear==NonLinear_lens .or. P%NonLinear==NonLinear_both) .and. P%DoLensing) then
P%WantTransfer = .true.
call Transfer_SetForNonlinearLensing(P%Transfer)
end if
call Transfer_SortAndIndexRedshifts(P%Transfer)
!JD 08/13 end changes
P%transfer%kmax=P%transfer%kmax*(P%h0/100._dl)
Ini_fail_on_not_found = .false.
DebugParam = Ini_Read_Double('DebugParam',DebugParam)
ALens = Ini_Read_Double('Alens',Alens)
call Reionization_ReadParams(P%Reion, DefIni)
call InitialPower_ReadParams(P%InitPower, DefIni, P%WantTensors)
call Recombination_ReadParams(P%Recomb, DefIni)
if (Ini_HasKey('recombination')) then
i = Ini_Read_Int('recombination',1)
if (i/=1) error stop 'recombination option deprecated'
end if
call Bispectrum_ReadParams(BispectrumParams, DefIni, outroot)
if (P%WantScalars .or. P%WantTransfer) then
P%Scalar_initial_condition = Ini_Read_Int('initial_condition',initial_adiabatic)
if (P%Scalar_initial_condition == initial_vector) then
P%InitialConditionVector=0
numstr = Ini_Read_String('initial_vector',.true.)
read (numstr,*) P%InitialConditionVector(1:initial_iso_neutrino_vel)
end if
if (P%Scalar_initial_condition/= initial_adiabatic) use_spline_template = .false.
end if
if (P%WantScalars) then
ScalarFileName = trim(outroot)//Ini_Read_String('scalar_output_file')
LensedFileName = trim(outroot) //Ini_Read_String('lensed_output_file')
LensPotentialFileName = Ini_Read_String('lens_potential_output_file')
if (LensPotentialFileName/='') LensPotentialFileName = concat(outroot,LensPotentialFileName)
ScalarCovFileName = Ini_Read_String_Default('scalar_covariance_output_file','scalCovCls.dat',.false.)
if (ScalarCovFileName/='') then
has_cl_2D_array = .true.
ScalarCovFileName = concat(outroot,ScalarCovFileName)
end if
end if
if (P%WantTensors) then
TensorFileName = trim(outroot) //Ini_Read_String('tensor_output_file')
if (P%WantScalars) then
TotalFileName = trim(outroot) //Ini_Read_String('total_output_file')
LensedTotFileName = Ini_Read_String('lensed_total_output_file')
if (LensedTotFileName/='') LensedTotFileName= trim(outroot) //trim(LensedTotFileName)
end if
end if
if (P%WantVectors) then
VectorFileName = trim(outroot) //Ini_Read_String('vector_output_file')
end if
#ifdef WRITE_FITS
if (P%WantCls) then
FITSfilename = trim(outroot) //Ini_Read_String('FITS_filename',.true.)
if (FITSfilename /='') then
inquire(file=FITSfilename, exist=bad)
if (bad) then
open(unit=18,file=FITSfilename,status='old')
close(18,status='delete')
end if
end if
end if
#endif
Ini_fail_on_not_found = .false.
!optional parameters controlling the computation
P%AccuratePolarization = Ini_Read_Logical('accurate_polarization',.true.)
P%AccurateReionization = Ini_Read_Logical('accurate_reionization',.false.)
P%AccurateBB = Ini_Read_Logical('accurate_BB',.false.)
P%DerivedParameters = Ini_Read_Logical('derived_parameters',.true.)
version_check = Ini_Read_String('version_check')
if (version_check == '') then
!tag the output used parameters .ini file with the version of CAMB being used now
call TNameValueList_Add(DefIni%ReadValues, 'version_check', version)
else if (version_check /= version) then
write(*,*) 'WARNING: version_check does not match this CAMB version'
end if
!Mess here to fix typo with backwards compatibility
if (Ini_HasKey('do_late_rad_trunction')) then
DoLateRadTruncation = Ini_Read_Logical('do_late_rad_trunction',.true.)
if (Ini_HasKey('do_late_rad_truncation')) error stop 'check do_late_rad_xxxx'
else
DoLateRadTruncation = Ini_Read_Logical('do_late_rad_truncation',.true.)
end if
if (HighAccuracyDefault) then
DoTensorNeutrinos = .true.
else
DoTensorNeutrinos = Ini_Read_Logical('do_tensor_neutrinos',DoTensorNeutrinos )
end if
FeedbackLevel = Ini_Read_Int('feedback_level',FeedbackLevel)
output_file_headers = Ini_Read_Logical('output_file_headers',output_file_headers)
P%MassiveNuMethod = Ini_Read_Int('massive_nu_approx',Nu_best)
ThreadNum = Ini_Read_Int('number_of_threads',ThreadNum)
use_spline_template = Ini_Read_Logical('use_spline_template',use_spline_template)
if (do_bispectrum) then
lSampleBoost = 50
else
lSampleBoost = Ini_Read_Double('l_sample_boost',lSampleBoost)
end if
if (outroot /= '') then
if (InputFile /= trim(outroot) //'params.ini') then
call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
else
write(*,*) 'Output _params.ini not created as would overwrite input'
end if
end if
call Ini_Close
if (.not. CAMB_ValidateParams(P)) error stop 'Stopped due to parameter error'
#ifdef RUNIDLE
call SetIdle
#endif
if (global_error_flag==0) call CAMB_GetResults(P)
if (global_error_flag/=0) then
write(*,*) 'Error result '//trim(global_error_message)
error stop
endif
if (P%PK_WantTransfer) then
call Transfer_SaveToFiles(MT,TransferFileNames)
call Transfer_SaveMatterPower(MT,MatterPowerFileNames)
call Transfer_output_sig8(MT)
end if
if (P%WantCls) then
call output_cl_files(ScalarFileName, ScalarCovFileName, TensorFileName, TotalFileName, &
LensedFileName, LensedTotFilename, output_factor)
call output_lens_pot_files(LensPotentialFileName, output_factor)
if (P%WantVectors) then
call output_veccl_files(VectorFileName, output_factor)
end if
#ifdef WRITE_FITS
if (FITSfilename /= '') call WriteFitsCls(FITSfilename, CP%Max_l)
#endif
end if
call CAMB_cleanup
stop
100 stop 'Must give num_massive number of integer physical neutrinos for each eigenstate'
end program driver
#ifdef RUNIDLE
!If in Windows and want to run with low priorty so can multitask
subroutine SetIdle
USE DFWIN
Integer dwPriority
Integer CheckPriority
dwPriority = 64 ! idle priority
CheckPriority = SetPriorityClass(GetCurrentProcess(), dwPriority)
end subroutine SetIdle
#endif