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
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
|
! Copyright (C) 2005 Barbara Ercolano
!
! Version 2.02
module set_input_mod
use common_mod
implicit none
character(len=50) :: shapeDiffuse ! diffuse spectrum shape
contains
subroutine readInput()
implicit none
! local variables
integer :: err ! allocation error status
integer :: i, j ! counters
integer :: ios ! I/O error status
integer :: maxLim = 1000 ! max limit for reading loop
character(len=50) :: cValue ! character value corresponding to keyword
character(len=50) :: in_file ! input file
character(len=50) :: keyword ! input parameter keyword
character(len=50) :: multiPhotoSources ! stars file
logical :: lgOutputExists ! to check for existence of output directory
real :: nphotonsin ! to allow scientific notation in input
! set default values and set non oprional values to 0 or 0. or "zero"
lgIsotropic = .false.
lgRecombination = .false.
lgAutoPackets = .false.
lgMultiDustChemistry = .false.
lgMultiChemistry = .false.
lgTalk = .false.
lgHdenConstant= .false.
lgDfile = .false.
lgDebug = .false.
lgDlaw = .false.
lgDust = .false.
lgDustConstant= .false.
lgGas = .true.
lgMdMg = .false.
lgMdMh = .false.
lgNeInput = .false.
lgNeutral = .true.
lgOutput = .false.
lgPlaneIonization = .false.
lg1D = .false.
lgDustScattering = .true.
lgSymmetricXYZ= .false.
lgEcho = .false.
lgNosource = .false.
lgforceTDust = .false.
lgOutputExists= .false.
nPhotonsDiffuse = 0
nStars = 0
nxIn(:) = 0
nyIn(:) = 0
nzIn(:) = 0
nxIn(1) = 30
nyIn(1) = 30
nzIn(1) = 30
maxIterateMC = 30
maxPhotons = 0
nbins = 600
MdMgValue = 0.
MdMgFile = "none"
nAngleBins = 0
nGrids = 1
NdustValue = 0.
NdustFile = "none"
nstages = 7
! qheat
Tmax =700.
nTbins =300
lgWritePss =.false.
minaQHeat = 1.e-3
minConvQHeat = 99.
fillingFactor = 1.
contCube = -1.
convIncPercent= 0.
convWriteGrid = 0.
nPhotIncrease = 0.
densityLaw = 0.
multiPhotoSources = "none"
densityFile = "none"
dustFile = "none"
gridList = "none"
nu0 = 0.
Ldiffuse = 0.
Tdiffuse = 0.
shapeDiffuse = "none"
Hdensity = 0.
H0Start = 3.e-5
LPhot = 0.
minConvergence= 95.
NeStart = 0.
nuMax = 15.
nuMin = 1.001e-5
nuStepSize = 0.075
resLinesTransfer = 101.
Rnx = -1.
Rny = -1.
Rnz = -1.
R_in = -1.
R_out = 0.
TeStart = 10000.
XHILimit = 0.05
! ask the user to input file name and check for errors
in_file = 'input/input.in'
close(10)
open (unit = 10, file = in_file, status = "old", position="rewind", iostat = ios, action="read")
if (ios /= 0) then
print*, "!out_in: can't open file input/input.in - terminating"
stop
end if
inquire(file='./output/.', exist=lgOutputExists)
if ( .not. lgOutputExists ) then
print *,"!output directory output/ doesn't exist - terminating"
stop
end if
print*, taskid, "input file: ", in_file
do i = 1, maxLim
! check if the end of file has been reached
read(unit=10, fmt=*, iostat=ios) keyword
if (ios < 0) exit ! end of file reached
select case (keyword)
case ("TDust") ! BEKS2010
backspace 10
read(unit=10,fmt=*, iostat=ios) keyword, forceTDust
lgforceTDust = .true.
!print*,keyword,lgforceTDust,forceTDust
case ("echo")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, echot1, echot2, echoTemp
!print "(1X,A4,1X,3es12.4)",keyword,echot1,echot2,echoTemp
print*, '! readInput [Reminder]: Light echo option - is the inclination correctly specified? [see Readme]'
print*, '! readInput [Warning]: Light echo option - '
print*, '! The actual volume of material illuminated by the echo is probably'
print*, '! much smaller than that illuminated in this grid. Be sure to perform'
print*, '! volume correction. See echo.out file for volume illuminated in this grid.'
! echot1 = echot1 * 2.59020684e15 !ct
! echot2 = echot2 * 2.59020684e15
lgEcho = .true.
case("NoSourceSED")
lgNosource = .true.
!print*, keyword, lgNosource
case("isotropicScattering")
lgIsotropic = .true.
!print*, keyword, lgIsotropic
case("2D")
lg2D = .true.
!print*, keyword, lg2D
case("nstages")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nstages
if (nstages>10) then
print*, " !out_in: Max 10 ionisation stages allowed. &
&If more are required (and you have the relevant atomic data), &
& please update the data/ directory (and let B. Ercolano know ;-) &
& and edit the data/fileNames.dat file"
stop
end if
!print*, keyword, nstages
case ("getEquivalentTau")
lgEquivalentTau = .true.
!print*, keyword, lgEquivalentTau
case ("diffuseSource")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, Ldiffuse, Tdiffuse, shapeDiffuse, nPhotonsDiffuse, emittingGrid
!print*, keyword, Ldiffuse, Tdiffuse, shapeDiffuse, nPhotonsDiffuse, emittingGrid
case ("traceHeating")
lgTraceHeating = .true.
!print*, keyword, lgTraceHeating
case("quantumHeatGrain")
lgQheat = .true.
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword,minaQHeat,minConvQHeat
!print*, keyword, minaQHeat,minConvQHeat
case("quantumHeatGrainParameters")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword,Tmax,nTbins,lgWritePss
!print*, keyword,Tmax,nTbins,lgWritePss
case ("noPhotoelectric")
!print*, keyword
lgPhotoelectric = .false.
case ("multiPhotoSources")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword,multiPhotoSources
!print*, keyword,multiPhotoSources
lgMultiStars=.true.
case("continuumCube")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, contCube(1), contCube(2)
!print*, keyword, contCube
case("noScattering")
!print*, keyword
lgDustScattering = .false.
case("resLinesTransfer")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, resLinesTransfer
!print*, keyword, resLinesTransfer
case("multiGrids")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nGrids, gridList
!print*, keyword, nGrids, gridList
if (nGrids > maxGrids) then
print*, "! readInput: please increse the maxGrids parameter in constants_mod.f90", &
& nGrids, maxGrids
stop
end if
call readGridList(gridList)
case("recombinationLines")
lgRecombination=.true.
!print*, keyword
case ("fillingFactor")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, fillingFactor
!print*, keyword, fillingFactor
case ("slit")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, dxSlit, dySlit
!print*, keyword, dxSlit, dySlit
case ("inputNe")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lgNeInput = .true.
!print*, keyword
case ("edges")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, Rnx, Rny, Rnz
!print*, keyword, Rnx, Rny, Rnz
case ("output")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lgOutput = .true.
!print*, keyword
case ("multiChemistry")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nAbComponents
lgMultiChemistry = .true.
backspace(10)
allocate(abundanceFile(1:nAbComponents))
read(unit=10, fmt=*, iostat=ios) keyword, nAbComponents, &
& (abundanceFile(j), j=1,nAbComponents)
!print*, keyword, nAbComponents, (abundanceFile(j), j=1,nAbComponents)
case ("multiDustChemistry")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nDustComponents
lgMultiDustChemistry = .true.
lgDust = .true.
backspace(10)
allocate(dustSpeciesFile(1:nDustComponents))
allocate(nSpeciesPart(1:nDustComponents))
read(unit=10, fmt=*, iostat=ios) keyword, nDustComponents, (dustSpeciesFile(j), j=1,nDustComponents), dustFile(2)
!print*, keyword, nDustComponents, (dustSpeciesFile(j), j=1,nDustComponents)
print*, dustFile(2)
case ("planeIonization")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, meanField, nu0
lgPlaneIonization = .true.
!print*, keyword, meanField, nu0
allocate(Lstar(1))
Lstar=0.
case ("autoPackets")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, convIncPercent, nPhotIncrease, nphotonsin
maxPhotons=int(nphotonsin,int64)
lgAutoPackets = .true.
!print*, keyword, convIncPercent, nPhotIncrease, maxPhotons
if (maxPhotons .eq. 0) then
print*,'! readInput: autoPackets input invalid - maximum number of photons is zero'
stop
endif
case ("oneD")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lg1D = .true.
!print*, keyword
print*, 'ERROR: oneD option is not currently avalable!!!!'
stop
case ("densityFile")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, densityFile
lgDfile = .true.
!print*, keyword, densityFile
open(unit=11,file=densityFile,status='old')
read(unit=11, fmt=*, iostat=ios) keyword
if (keyword=="#") then
backspace 11
read(unit=11,fmt=*, iostat=ios) keyword,nxin(1),nyin(1),nzin(1)
print*,'grid axes read in from file:',trim(densityFile),nxin(1),nyin(1),nzin(1)
endif
close(11)
case ("dustFile")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, dustFile(1), dustFile(2)
!print*, keyword, trim(dustFile(1))," | ",trim(dustFile(2))," |"
allocate(dustSpeciesFile(1:1))
allocate(nspeciesPart(1:1))
nDustComponents = 1
dustSpeciesFile(1) = dustFile(1)
case ("debug")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lgDebug = .true.
!print*, keyword
case ("nbins")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nbins
!print*, keyword, nbins
case ("talk")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lgTalk = .true.
!print*, keyword
case ("symmetricXYZ")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
lgSymmetricXYZ = .true.
!print*, keyword
case ("contShape")
backspace 10
nstars=1
allocate(contShape(0:1))
allocate(contShapeIn(0:1))
allocate(spID(0:1))
allocate(tstep(0:1))
contShape = 'none'
read(unit=10, fmt=*, iostat=ios) keyword, contShape(1)
if (contShape(1) == 'powerlaw') then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, contShape(1), pwlIndex
!print*, keyword, contShape(1), pwlIndex
else
!print*, keyword, contShape(1)
end if
if (contShape(1) == 'blackbody' .or. contShape(1) == 'powerlaw') then
spID(1) = contShape(1)
else
spID(1) = 'rauch'
end if
tstep = 0.
case ("nebComposition")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, cValue
if (cValue=="noGas") lgGas=.false.
!print*, keyword, cValue
allocate(abundanceFile(1:nAbComponents))
abundanceFile = cValue
case ("Hdensity")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, Hdensity
lgHdenConstant = .true.
!print*, keyword, Hdensity
case ("densityLaw")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, (densityLaw(j),j=1,3)
lgDlaw = .true.
!print*, keyword, densityLaw
case ("maxIterateMC")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, maxIterateMC, minConvergence
!print*, keyword, maxIterateMC, minConvergence
case ("nPhotons")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nphotonsin
nPhotonsTot=int(nphotonsin,int64)
!print*, keyword, nPhotonsTot
case ("nx")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nxin(1)
case ("ny")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nyin(1)
case ("nz")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nzin(1)
case ("starPosition")
backspace 10
allocate(starPosition(1))
read(unit=10, fmt=*, iostat=ios) keyword, starPosition(1)%x, starPosition(1)%y, starPosition(1)%z
print*, keyword, starPosition(1)%x, starPosition(1)%y, starPosition(1)%z
case ("TeStart")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, TeStart
!print*, keyword, TeStart
case ("NeStart")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, NeStart
!print*, keyword, NeStart
case ("TStellar")
backspace 10
allocate(TStellar(0:1))
read(unit=10, fmt=*, iostat=ios) keyword, TStellar(1)
!print*, keyword, TStellar(1)
case ("convLimit")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, XHILimit
!print*, keyword, XHILimit
case ("H0Start")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, H0Start
!print*, keyword, H0Start
case ("LStar")
backspace 10
allocate(Lstar(0:1))
Lstar = 0.
read(unit=10, fmt=*, iostat=ios) keyword, LStar(1)
!print*, keyword, LStar(1)
case ("LPhot")
backspace 10
if (.not. allocated(LStar)) allocate(Lstar(0:1))
Lstar = 0.
read(unit=10, fmt=*, iostat=ios) keyword, LPhot
!print*, keyword, LPhot
case ("nuMax")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nuMax
!print*, keyword, nuMax
case ("nuMin")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nuMin
!print*, keyword, nuMin
case ("Rin")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, R_in
!print*, keyword, R_in
case ("Rout")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, R_out
!print*, keyword, R_out
case ("MdMg")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword
if (keyword=="constant") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, MdMgValue
print*, 'MdMg constant ', MdMgValue
lgDustConstant = .true.
lgDust = .true.
else if (keyword=="file") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, MdMgFile
print*, 'MdMg file ', MdMgFile
lgDust = .true.
open(unit=11,file=MdMgFile,status='old')
read(unit=11, fmt=*, iostat=ios) keyword
if (keyword=="#") then
backspace 11
read(unit=11,fmt=*, iostat=ios) keyword,nxin(1),nyin(1),nzin(1)
print*,'grid axes read in from file:',trim(MdMgFile),nxin(1),nyin(1),nzin(1)
endif
close(11)
else
print*, "! readInput: invalid keyword in MdMg field ", keyword
stop
end if
lgMdMg = .true.
case ("MdMh")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword
if (keyword=="constant") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, MdMgValue
print*, 'MdMh constant ', MdMgValue
lgDustConstant = .true.
lgDust = .true.
else if (keyword=="file") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, MdMgFile
print*, 'MdMh file ', MdMgFile
lgDust = .true.
open(unit=11,file=MdMgFile,status='old')
read(unit=11, fmt=*, iostat=ios) keyword
if (keyword=="#") then
backspace 11
read(unit=11,fmt=*, iostat=ios) keyword,nxin(1),nyin(1),nzin(1)
print*,'grid axes read in from file:',trim(MdMgFile),nxin(1),nyin(1),nzin(1)
endif
close(11)
else
print*, "! readInput: invalid keyword in MdMh field ", keyword
stop
end if
lgMdMh = .true.
case ("Ndust")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword
if (keyword=="constant") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, NdustValue
print*, 'Ndust constant ', NdustValue
lgDustConstant = .true.
lgDust = .true.
else if (keyword=="file") then
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, keyword, NdustFile
print*, 'Ndust file ', NdustFile
open(unit=11,file=NdustFile,status='old')
read(unit=11, fmt=*, iostat=ios) keyword
if (keyword=="#") then
backspace 11
read(unit=11,fmt=*, iostat=ios) keyword,nxin(1),nyin(1),nzin(1)
print*,'grid axes read in from file:',trim(NdustFile),nxin(1),nyin(1),nzin(1)
endif
close(11)
lgDust = .true.
else
print*, "! readInput: invalid keyword in Ndust field ", keyword
stop
end if
case ("writeGrid")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, convWriteGrid
!print*, keyword, convWriteGrid
case ("inclination")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nAngleBins
if (nAngleBins>2) then
print*, '! readInput: only two inclination anges are allowed per simulation. &
& If more inclinations are required please run mocassinWarm for the other&
& viewing angles.'
stop
end if
allocate(viewPointTheta(0:nAngleBins), stat=err)
if (err /= 0) then
print*, '! readInput: allocation error for viewPointTheta pointer'
stop
end if
viewPointTheta=0.
allocate(viewPointPhi(0:nAngleBins), stat=err)
if (err /= 0) then
print*, '! readInput: allocation error for viewPointPhi pointer'
stop
end if
viewPointPhi=0.
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, nAngleBins, (viewPointTheta(j), viewPointPhi(j), j=1,nAngleBins)
!print*, keyword, nAngleBins, (viewPointTheta(j), viewPointPhi(j), j = 1, nAngleBins)
case ("dustMass")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, inputDustMass
lginputDustMass = .true.
case ("gasMass")
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword, inputGasMass
lginputGasMass = .true.
print*,"Using gasMass keyword. Be sure to specify gas density file or Hdensity 1.0"
case default
backspace 10
read(unit=10, fmt=*, iostat=ios) keyword
if (keyword(1:1).ne."#") then
print*, "! readInput: unrecognised keyword in model parameter input file", &
& in_file, keyword
endif
end select
end do
close(10)
if (nPhotonsTot .eq. 0 .and. nPhotonsDiffuse .eq. 0) then
print*,'! readInput: no photons in model.'
stop
endif
allocate(lgDataAvailable(3:nElements, 1:nstages), stat=err)
if (err /= 0) then
print*, '! readInput: allocation error for lgDataAvailable pointer'
stop
end if
if (lginputDustMass) then
inputDustMass = inputDustMass / 5.028e11
if (lgSymmetricXYZ) then
inputDustMass = inputDustMass / 8
end if
end if
if (lginputGasMass) then
inputGasMass = inputGasMass / 5.028e11
if (lgSymmetricXYZ) then
inputGasMass = inputGasMass / 8
end if
end if
if (lginputDustMass .and. .not.lgDust) then
print*,"! readInput: you have asked to scale dustMass but there is no dust"
stop
endif
if (lginputGasMass .and. .not.lgGas) then
print*,"! readInput: you have asked to scale gasMass but there is no gas"
stop
endif
if (lginputDustMass .and. lginputGasMass) then
if (lgDust .and. lgGas) then
if (lgMdMg .or. lgMdMh) then
print*,'! readInput: dust mass is defined by your gas-to-dust ratio'
print*,'! readInput: please remove dustMass keyword from input.in'
stop
else
print*,'! readInput: Simulation will scale gas to gasMass and dust to dustMass.'
print*,'! readInput: This will change your gas-to-dust ratio so BE SURE THIS IS WHAT YOU WANT...'
endif
endif
endif
if (lgMultiStars) then
call setMultiPhotoSources(multiPhotoSources)
else if (lgMultiStars .and. nStars==1) then
print*, '! readInput: multiPhotoSources keyword and Lstar, Tstellar, ContShape are cannot be specified together'
stop
else
allocate(nPhotons(1))
allocate(deltaE(0:1))
deltaE = 0.
nPhotons(1)=nPhotonsTot
if (Ldiffuse>0. .and. nPhotonsDiffuse>0) then
if (nStars<1) then
if(.not.allocated(contShape)) allocate(contShape(0:0))
if(.not.allocated(spID)) allocate(spID(0:0))
if(.not.allocated(LStar)) allocate(LStar(0:0))
end if
deltaE(0) = Ldiffuse/nPhotonsDiffuse
contShape(0) = shapeDiffuse
if (contShape(0) == 'blackbody' .or. contShape(0) == 'powerlaw') then
spID(0) = contShape(0)
else
spID(0) = 'rauch'
end if
end if
if (.not.allocated(starPosition)) then
allocate (starPosition(1))
starPosition(1)%x=0.
starPosition(1)%y=0.
starPosition(1)%z=0.
print* , '! readInput: [talk] ionising source positioned at 0.,0.,0.'
end if
end if
if (.not.allocated(contShape)) then
print*, "! readInput: no information about the continuum shape - contShape -"
stop
end if
if (.not.allocated(TStellar)) then
print*, "! readInput: no information about the stellar effective temperature - Tstellar -"
stop
end if
if (.not.allocated(nPhotons)) then
print*, "! readInput: no information about the number of energy packets to be used - nPhotons -"
stop
end if
if (.not.allocated(Lstar) .and. Lphot<=0.) then
print*, "! readInput: no information about the stellar luminosity - Lstar or Lphot -"
stop
end if
contShapeIn = contShape
print*, "mother nx, ny, nz" , nxin(1), nyin(1), nzin(1)
! check for missing or invalid values in the model parameters input file
if (lg1D) then
nyIn = 1
nzIn = 1
lgSymmetricXYZ = .true.
else if (lgPlaneIonization .and. lg2d) then
print*, "! readInput: planeIonization and 2D options are not compatible"
stop
else if (resLinesTransfer >= minConvergence .and. resLinesTransfer /= 101. .and. lgDust) then
print*, "! readInput: the min convergence level assigned to the calculation of &
&the resonant lines transfer is higher or equal to that assigned by maxIterateMC"
stop
else if (resLinesTransfer >= 100.5 .and. lgDust .and. lgGas) then
print*, "! readInput: [warning] both dust and gas processes have been activated, but the &
&resonanceLinesTransfer keyword is absent. The dust temperatures may be underestimated!!!!"
else if (resLinesTransfer <= 100 .and. (.not.lgDust .or. .not.lgGas)) then
print*, "! readInput: [warning] both dust and gas processes must be activated, to use the&
&resonanceLinesTransfer keyword. Reset to 101."
resLinesTransfer = 101.
else if (.not.allocated(abundanceFile) .and. lgGas) then
print*, "! readInput: no elemental abundance file has been specified"
stop
else if ((lgMdMh .and. lgMdMg) .and. (NdustFile /= "none" .or. NdustValue /=0) )then
print*, "! readInput: MdMg and Ndust cannot be specified together"
stop
else if ( (.not.lgDust) .and. (.not.lgGas) )then
print*, "! readInput: the grid is completely empty. no gas or dust present.", &
lgGas,lgDust
stop
else if ( lgDebug .and. .not.(lgGas) ) then
print*, "! readInput: debugging mode can only be used when a gas component is present. &
& debugging mode will be turned off."
lgDebug = .false.
else if (fillingFactor > 1.) then
print*, "! readInput: filling factor, epsilon, specified is greater than unity.", &
& fillingFactor
stop
else if ((.not.lgDfile) .and. (Rnx<0. .or. Rny<0. .or. Rnz<0.) .and. (NdustFile == "none")) then
print*, "! readInput: Grid edges unspecified or non-valid grid &
&edges", Rnx,Rny,Rnz
stop
else if (lgDfile .and. (Rnx>=0. .or. Rny>=0. .or. Rnz>=0.) ) then
print*, "! readInput: [warning] when the density distribution is&
& specified via an external density file, no grid edges need &
& to be specified. The specified values will be ignored."
else if (lgPlaneIonization .and. .not.lgDfile) then
print*, "! readInput: plane ionizing field specified - must use user defined &
&density distribution via the densityFile option"
stop
else if (lgPlaneIonization .and. lgSymmetricXYZ) then
print*, "! readInput [warning]: plane ionizing field specified - cannot use&
& symmetricXYZ - removed."
lgSymmetricXYZ = .false.
else if (TStellar(1) == 0. .and. Tdiffuse==0. .and. size(contShape).ge.1) then
if (contShape(1) /= 'powerlaw') then
print*, "! readInput: TStellar missing from model parameter input file"
stop
endif
else if (R_in < 0.) then
print*, "! readInput: Invalid Rin parameter in the input file", R_in
stop
else if (LPhot == 0. .and. .not.allocated(LStar) .and. .not.lgPlaneIonization &
& .and. Ldiffuse<=0.) then
print*, "! readInput: LPhot and LStar missing from model parameter input file"
stop
else if (LPhot /= 0. .and. allocated(LStar)) then
print*, "! readInput: [warning] LPhot and LStar both entered. "
else if (maxIterateMC < 1) then
print*, "! readInput: invalid maxIterateMC in model parameter input &
& file", in_file, maxIterateMC
stop
else if (minConvergence < 0. .or. minConvergence > 100.) then
print*, "! readInput: invalid minConvergence input (must be between 0. and 100.) ", &
& in_file, minConvergence
stop
else if (nPhotons(1) < 0) then
print*, "! readInput: invalid nPhotons in model parameter input &
& file", in_file, nPhotons
stop
else if (lgDust) then
if ((dustFile(1) == "none" .or. dustFile(2)=="none") .and. &
&.not.lgMultiDustChemistry) then
print*, "! readInput: dust files were not specified; please include dustFile field &
& in the input files and specify name and location of dust files"
stop
end if
else if (nxin(1) < 1 .or. nyin(1) < 1 .or. nzin(1) < 1) then
print*, "! readInput: invalid grid boundaries in model parameter input &
& file", in_file, nxin(1), nyin(1), nzin(1)
stop
else if (TStellar(1) < 0.) then
print*, "! readInput: invalid TStellar in model parameter input &
& file", in_file, TStellar
stop
else if (XHILimit < 0.) then
print*, "! readInput: invalid XHILimit in model parameter input &
& file", in_file, XHILimit
stop
else if (LPhot < 0.) then
print*, "! readInput: invalid LPhot in model parameter input &
& file", in_file, LPhot
stop
else if (LStar(1) < 0.) then
print*, "! readInput: invalid LStar in model parameter input &
& file", in_file, LStar
stop
else if (nuMax < 0.) then
print*, "! readInput: invalid nuMax in model parameter input &
& file", in_file, nuMax
stop
else if (nuMin < 0.) then
print*, "! readInput: invalid nuMin in model parameter input &
& file", in_file, nuMin
stop
else if (nuMin > 1.4e-6) then
print*, "! readInput [waring]: nuMin in model parameter input &
& file is too large to include 5GHz flux calculations", in_file,nuMin
else if (R_out < 0.) then
print*, "! readInput: invalid R_out in model parameter input &
& file", in_file, R_out
stop
else if (R_in <0.) then
print*, "! readInput: invalid R_in in model parameter input &
& file", in_file, R_in
stop
else if ((R_in >= R_out) .and. (R_out > 0.) ) then
print*, "! readInput: the inner radius must be smaller than &
&the outer radius ", in_file, R_in, R_out
stop
else if (lgMultiChemistry .and. .not.lgDfile) then
print*, "! readInput: the MultiChemistry option has been specified &
&therefore an user density file of the type x,y,z,den,abIndex must also &
&be specified."
stop
end if
TStellar(0) = Tdiffuse
contains
subroutine readGridList(filename)
implicit none
real :: skipR
character(len=50), intent(in) :: filename
character(len=50) :: skipC
integer :: i
integer :: ios
integer :: skip
close(17)
open (unit = 17, file = filename, status = "old", position="rewind", iostat = ios, action="read")
if (ios/=0) then
print*, "! readGridList: can't open file for reading ", filename
stop
end if
do i = 2, nGrids
read(17, *) skip, nxIn(i), nyIn(i), nzIn(i), skipC, skipR
read(17, *) skipR, skipR, skipR, skipR, skipR, skipR
end do
close(17)
end subroutine readGridList
end subroutine readInput
subroutine setMultiPhotoSources(infile)
implicit none
character(len=50), intent(in) :: infile
integer :: ios, i
close(13)
open(file=infile, unit=13, iostat=ios, action="read")
if (ios /= 0) then
print*, "! setMultiPhotoSources: can't open ionising sources file", infile
stop
end if
read(13,*) nStars
print*, "Multiple Ionising sources", nStars
print*, "(i, Tstellar, Lstar, ContShape, Nphotons, position)"
allocate(deltaE(0:nStars))
allocate(spID(0:nStars))
allocate(tStep(nStars))
allocate(TStellar(0:nStars))
allocate(Lstar(nStars))
allocate(ContShape(0:nStars))
allocate(ContShapeIn(0:nStars))
allocate(nPhotons(nStars))
allocate(starPosition(nStars))
TStellar=0.
Lstar=0.
ContShape='none'
ContShapeIn='none'
spID='none'
tStep=0.
nPhotons=0
do i = 1, nStars
read(13,*) TStellar(i), Lstar(i), ContShape(i), starPosition(i)%x, starPosition(i)%y, starPosition(i)%z, &
& spID(i), tStep(i)
contShapeIn(i) = contShape(i)
end do
print*, 'i, Tstellar(i), Lstar(i), contShape(i), Nphotons(i), starPosition(i), deltaE(i)'
do i = 1, nStars
nPhotons(i) = NphotonsTot/nStars
deltaE(i) = Lstar(i)/nPhotons(i)
print*, i, Tstellar(i), Lstar(i), contShape(i), Nphotons(i), starPosition(i), deltaE(i)
end do
if (Ldiffuse>0. .and. nPhotonsDiffuse>0) then
deltaE(0) = Ldiffuse/nPhotonsDiffuse
contShape(0) = shapeDiffuse
if (contShape(0) == 'blackbody' .or. contShape(0) == 'powerlaw') then
spID(0) = contShape(0)
else
spID(0) = 'rauch'
end if
end if
end subroutine setMultiPhotoSources
end module set_input_mod
|
