19 double precision ::
t_aia(1:101)
22 double precision ::
f_335(1:101)
35 data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
36 4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
37 4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
38 5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
39 5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
40 6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
41 6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
42 6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
43 7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
44 7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
45 8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
46 8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
47 8.8, 8.85, 8.9, 8.95, 9. /
49 data f_94 / 4.25022959
d-37, 4.35880298
d-36, 3.57054296
d-35, 2.18175426
d-34, &
50 8.97592571
d-34, 2.68512961
d-33, 7.49559346
d-33, 2.11603751
d-32, &
51 5.39752853
d-32, 1.02935904
d-31, 1.33822307
d-31, 1.40884290
d-31, &
52 1.54933156
d-31, 2.07543102
d-31, 3.42026227
d-31, 6.31171444
d-31, &
53 1.16559416
d-30, 1.95360497
d-30, 2.77818735
d-30, 3.43552578
d-30, &
54 4.04061803
d-30, 4.75470982
d-30, 5.65553769
d-30, 6.70595782
d-30, &
55 7.80680354
d-30, 8.93247715
d-30, 1.02618156
d-29, 1.25979030
d-29, &
56 1.88526483
d-29, 3.62448572
d-29, 7.50553279
d-29, 1.42337571
d-28, &
57 2.37912813
d-28, 3.55232305
d-28, 4.84985757
d-28, 6.20662827
d-28, &
58 7.66193687
d-28, 9.30403645
d-28, 1.10519802
d-27, 1.25786927
d-27, &
59 1.34362634
d-27, 1.33185242
d-27, 1.22302081
d-27, 1.05677973
d-27, &
60 9.23064720
d-28, 8.78570994
d-28, 8.02397416
d-28, 5.87681142
d-28, &
61 3.82272695
d-28, 3.11492649
d-28, 3.85736090
d-28, 5.98893519
d-28, &
62 9.57553548
d-28, 1.46650267
d-27, 2.10365847
d-27, 2.79406671
d-27, &
63 3.39420087
d-27, 3.71077520
d-27, 3.57296767
d-27, 2.95114380
d-27, &
64 2.02913103
d-27, 1.13361825
d-27, 5.13405629
d-28, 2.01305089
d-28, &
65 8.15781482
d-29, 4.28366817
d-29, 3.08701543
d-29, 2.68693906
d-29, &
66 2.51764203
d-29, 2.41773103
d-29, 2.33996083
d-29, 2.26997246
d-29, &
67 2.20316143
d-29, 2.13810001
d-29, 2.07424438
d-29, 2.01149189
d-29, &
68 1.94980213
d-29, 1.88917920
d-29, 1.82963583
d-29, 1.77116920
d-29, &
69 1.71374392
d-29, 1.65740593
d-29, 1.60214447
d-29, 1.54803205
d-29, &
70 1.49510777
d-29, 1.44346818
d-29, 1.39322305
d-29, 1.34441897
d-29, &
71 1.29713709
d-29, 1.25132618
d-29, 1.20686068
d-29, 1.14226584
d-29, &
72 1.09866413
d-29, 1.05635524
d-29, 1.01532444
d-29, 9.75577134
d-30, &
73 9.37102736
d-30, 8.99873335
d-30, 8.63860172
d-30, 8.29051944
d-30, &
76 data f_131 / 3.18403601
d-37, 3.22254703
d-36, 2.61657920
d-35, &
77 1.59575286
d-34, 6.65779556
d-34, 2.07015132
d-33, &
78 6.05768615
d-33, 1.76074833
d-32, 4.52633001
d-32, &
79 8.57121883
d-32, 1.09184271
d-31, 1.10207963
d-31, &
80 1.11371658
d-31, 1.29105226
d-31, 1.80385897
d-31, &
81 3.27295431
d-31, 8.92002136
d-31, 3.15214579
d-30, &
82 9.73440787
d-30, 2.22709702
d-29, 4.01788984
d-29, &
83 6.27471832
d-29, 8.91764995
d-29, 1.18725647
d-28, &
84 1.52888040
d-28, 2.05082946
d-28, 3.47651873
d-28, &
85 8.80482184
d-28, 2.66533063
d-27, 7.05805149
d-27, &
86 1.46072515
d-26, 2.45282476
d-26, 3.55303726
d-26, &
87 4.59075911
d-26, 5.36503515
d-26, 5.68444094
d-26, &
88 5.47222296
d-26, 4.81119761
d-26, 3.85959059
d-26, &
89 2.80383406
d-26, 1.83977650
d-26, 1.11182849
d-26, &
90 6.50748885
d-27, 3.96843481
d-27, 2.61876319
d-27, &
91 1.85525324
d-27, 1.39717024
d-27, 1.11504283
d-27, &
92 9.38169611
d-28, 8.24801234
d-28, 7.43331919
d-28, &
93 6.74537063
d-28, 6.14495760
d-28, 5.70805277
d-28, &
94 5.61219786
d-28, 6.31981777
d-28, 9.19747307
d-28, &
95 1.76795732
d-27, 3.77985446
d-27, 7.43166191
d-27, &
96 1.19785603
d-26, 1.48234676
d-26, 1.36673114
d-26, &
97 9.61047146
d-27, 5.61209353
d-27, 3.04779780
d-27, &
98 1.69378976
d-27, 1.02113491
d-27, 6.82223774
d-28, &
99 5.02099099
d-28, 3.99377760
d-28, 3.36279037
d-28, &
100 2.94767378
d-28, 2.65740865
d-28, 2.44396277
d-28, &
101 2.28003967
d-28, 2.14941419
d-28, 2.04178995
d-28, &
102 1.95031045
d-28, 1.87011994
d-28, 1.79777869
d-28, &
103 1.73093957
d-28, 1.66795789
d-28, 1.60785455
d-28, &
104 1.55002399
d-28, 1.49418229
d-28, 1.44022426
d-28, &
105 1.38807103
d-28, 1.33772767
d-28, 1.28908404
d-28, &
106 1.24196208
d-28, 1.17437501
d-28, 1.12854330
d-28, &
107 1.08410498
d-28, 1.04112003
d-28, 9.99529904
d-29, &
108 9.59358806
d-29, 9.20512291
d-29, 8.83009123
d-29, &
109 8.46817043
d-29, 8.11921928
d-29 /
111 data f_171 / 2.98015581
d-42, 1.24696230
d-40, 3.37614652
d-39, 5.64103034
d-38, &
112 5.20550266
d-37, 2.77785939
d-36, 1.16283616
d-35, 6.50007689
d-35, &
113 9.96177399
d-34, 1.89586076
d-32, 2.10982799
d-31, 1.36946479
d-30, &
114 6.27396553
d-30, 2.29955134
d-29, 7.13430211
d-29, 1.91024282
d-28, &
115 4.35358848
d-28, 7.94807808
d-28, 1.07431875
d-27, 1.08399488
d-27, &
116 9.16212938
d-28, 7.34715770
d-28, 6.59246382
d-28, 9.13541375
d-28, &
117 2.05939035
d-27, 5.08206555
d-27, 1.10148083
d-26, 2.01884662
d-26, &
118 3.13578384
d-26, 4.14367719
d-26, 5.36067711
d-26, 8.74170213
d-26, &
119 1.64161233
d-25, 2.94587860
d-25, 4.76298332
d-25, 6.91765639
d-25, &
120 9.08825111
d-25, 1.08496183
d-24, 1.17440114
d-24, 1.13943939
d-24, &
121 9.71696981
d-25, 7.09593688
d-25, 4.31376399
d-25, 2.12708486
d-25, &
122 8.47429567
d-26, 3.17608104
d-26, 1.95898842
d-26, 1.98064242
d-26, &
123 1.67706555
d-26, 8.99126003
d-27, 3.29773878
d-27, 1.28896127
d-27, &
124 8.51169698
d-28, 7.53520167
d-28, 6.18268143
d-28, 4.30034650
d-28, &
125 2.78152409
d-28, 1.95437088
d-28, 1.65896278
d-28, 1.68740181
d-28, &
126 1.76054383
d-28, 1.63978419
d-28, 1.32880591
d-28, 1.00833205
d-28, &
127 7.82252806
d-29, 6.36181741
d-29, 5.34633869
d-29, 4.58013864
d-29, &
128 3.97833422
d-29, 3.49414760
d-29, 3.09790940
d-29, 2.76786227
d-29, &
129 2.48806269
d-29, 2.24823367
d-29, 2.04016653
d-29, 1.85977413
d-29, &
130 1.70367499
d-29, 1.56966125
d-29, 1.45570643
d-29, 1.35964565
d-29, &
131 1.27879263
d-29, 1.21016980
d-29, 1.15132499
d-29, 1.09959628
d-29, &
132 1.05307482
d-29, 1.01040261
d-29, 9.70657096
d-30, 9.33214234
d-30, &
133 8.97689427
d-30, 8.63761192
d-30, 8.31149879
d-30, 7.85162401
d-30, &
134 7.53828281
d-30, 7.23559452
d-30, 6.94341530
d-30, 6.66137038
d-30, &
135 6.38929156
d-30, 6.12669083
d-30, 5.87346434
d-30, 5.62943622
d-30, &
138 data f_193 / 6.40066486
d-32, 4.92737300
d-31, 2.95342934
d-30, 1.28061594
d-29, &
139 3.47747667
d-29, 5.88554792
d-29, 7.72171179
d-29, 9.75609282
d-29, &
140 1.34318963
d-28, 1.96252638
d-28, 2.70163878
d-28, 3.63192965
d-28, &
141 5.28087341
d-28, 8.37821446
d-28, 1.39089159
d-27, 2.31749718
d-27, &
142 3.77510689
d-27, 5.85198594
d-27, 8.26021568
d-27, 1.04870405
d-26, &
143 1.25209374
d-26, 1.47406787
d-26, 1.77174067
d-26, 2.24098537
d-26, &
144 3.05926105
d-26, 4.50018853
d-26, 6.84720216
d-26, 1.00595861
d-25, &
145 1.30759222
d-25, 1.36481773
d-25, 1.15943558
d-25, 1.01467304
d-25, &
146 1.04092532
d-25, 1.15071251
d-25, 1.27416033
d-25, 1.38463476
d-25, &
147 1.47882726
d-25, 1.57041238
d-25, 1.69786224
d-25, 1.94970397
d-25, &
148 2.50332918
d-25, 3.58321431
d-25, 5.18061550
d-25, 6.60405549
d-25, &
149 6.64085365
d-25, 4.83825816
d-25, 2.40545020
d-25, 8.59534098
d-26, &
150 2.90920638
d-26, 1.33204845
d-26, 9.03933926
d-27, 7.78910836
d-27, &
151 7.29342321
d-27, 7.40267022
d-27, 8.05279981
d-27, 8.13829291
d-27, &
152 6.92634262
d-27, 5.12521880
d-27, 3.59527615
d-27, 2.69617560
d-27, &
153 2.84432713
d-27, 5.06697306
d-27, 1.01281903
d-26, 1.63526978
d-26, &
154 2.06759342
d-26, 2.19482312
d-26, 2.10050611
d-26, 1.89837248
d-26, &
155 1.66347131
d-26, 1.43071097
d-26, 1.21518419
d-26, 1.02078343
d-26, &
156 8.46936184
d-27, 6.93015742
d-27, 5.56973237
d-27, 4.38951754
d-27, &
157 3.38456457
d-27, 2.55309556
d-27, 1.88904224
d-27, 1.38057546
d-27, &
158 1.00718330
d-27, 7.43581116
d-28, 5.63562931
d-28, 4.43359435
d-28, &
159 3.63923535
d-28, 3.11248143
d-28, 2.75586846
d-28, 2.50672237
d-28, &
160 2.32419348
d-28, 2.18325682
d-28, 2.06834486
d-28, 1.93497044
d-28, &
161 1.84540751
d-28, 1.76356504
d-28, 1.68741425
d-28, 1.61566157
d-28, &
162 1.54754523
d-28, 1.48249410
d-28, 1.42020176
d-28, 1.36045230
d-28, &
165 data f_211 / 4.74439912
d-42, 1.95251522
d-40, 5.19700194
d-39, 8.53120166
d-38, &
166 7.72745727
d-37, 4.04158559
d-36, 1.64853511
d-35, 8.56295439
d-35, &
167 1.17529722
d-33, 2.16867729
d-32, 2.40472264
d-31, 1.56418133
d-30, &
168 7.20032889
d-30, 2.65838271
d-29, 8.33196904
d-29, 2.26128236
d-28, &
169 5.24295811
d-28, 9.77791121
d-28, 1.35913489
d-27, 1.43957785
d-27, &
170 1.37591544
d-27, 1.49029886
d-27, 2.06183401
d-27, 3.31440622
d-27, &
171 5.42497318
d-27, 8.41100374
d-27, 1.17941366
d-26, 1.49269794
d-26, &
172 1.71506074
d-26, 1.71266353
d-26, 1.51434781
d-26, 1.36766622
d-26, &
173 1.33483562
d-26, 1.36834518
d-26, 1.45829002
d-26, 1.62575306
d-26, &
174 1.88773347
d-26, 2.22026986
d-26, 2.54930499
d-26, 2.80758138
d-26, &
175 3.06176409
d-26, 3.62799792
d-26, 5.13226109
d-26, 8.46260744
d-26, &
176 1.38486586
d-25, 1.86192535
d-25, 1.78007934
d-25, 1.16548409
d-25, &
177 5.89293257
d-26, 2.69952884
d-26, 1.24891081
d-26, 6.41273176
d-27, &
178 4.08282914
d-27, 3.26463328
d-27, 2.76230280
d-27, 2.08986882
d-27, &
179 1.37658470
d-27, 8.48489381
d-28, 5.19304217
d-28, 3.19312514
d-28, &
180 2.02968197
d-28, 1.50171666
d-28, 1.39164218
d-28, 1.42448821
d-28, &
181 1.41714519
d-28, 1.33341059
d-28, 1.20759270
d-28, 1.07259692
d-28, &
182 9.44895400
d-29, 8.29030041
d-29, 7.25440631
d-29, 6.33479483
d-29, &
183 5.51563757
d-29, 4.79002469
d-29, 4.14990482
d-29, 3.59384972
d-29, &
184 3.12010860
d-29, 2.72624742
d-29, 2.40734791
d-29, 2.15543565
d-29, &
185 1.95921688
d-29, 1.80682882
d-29, 1.68695662
d-29, 1.59020936
d-29, &
186 1.50940886
d-29, 1.43956179
d-29, 1.37731622
d-29, 1.32049043
d-29, &
187 1.26771875
d-29, 1.21803879
d-29, 1.17074716
d-29, 1.10507836
d-29, &
188 1.06022834
d-29, 1.01703080
d-29, 9.75436986
d-30, 9.35349257
d-30, &
189 8.96744546
d-30, 8.59527489
d-30, 8.23678940
d-30, 7.89144480
d-30, &
192 data f_304 / 3.62695850
d-32, 2.79969087
d-31, 1.68340584
d-30, 7.32681440
d-30, &
193 1.99967770
d-29, 3.41296785
d-29, 4.55409104
d-29, 5.94994635
d-29, &
194 8.59864963
d-29, 1.39787633
d-28, 3.17701965
d-28, 1.14474920
d-27, &
195 4.44845958
d-27, 1.54785841
d-26, 4.70265345
d-26, 1.24524365
d-25, &
196 2.81535352
d-25, 5.10093666
d-25, 6.83545307
d-25, 6.82110329
d-25, &
197 5.66886188
d-25, 4.36205513
d-25, 3.29265688
d-25, 2.49802368
d-25, &
198 1.92527113
d-25, 1.51058572
d-25, 1.20596047
d-25, 9.76884267
d-26, &
199 7.89979266
d-26, 6.18224289
d-26, 4.67298332
d-26, 3.57934505
d-26, &
200 2.84535785
d-26, 2.32853022
d-26, 1.95228514
d-26, 1.67880071
d-26, &
201 1.47608785
d-26, 1.32199691
d-26, 1.20070960
d-26, 1.09378177
d-26, &
202 1.00031730
d-26, 9.62434001
d-27, 1.05063954
d-26, 1.27267143
d-26, &
203 1.45923057
d-26, 1.36746707
d-26, 1.03466970
d-26, 6.97647829
d-27, &
204 4.63141039
d-27, 3.19031994
d-27, 2.33373613
d-27, 1.81589079
d-27, &
205 1.48446917
d-27, 1.26611478
d-27, 1.12617468
d-27, 1.03625148
d-27, &
206 9.61400595
d-28, 8.79016231
d-28, 7.82612130
d-28, 6.73762960
d-28, &
207 5.59717956
d-28, 4.53010243
d-28, 3.65712196
d-28, 3.00958686
d-28, &
208 2.54011502
d-28, 2.18102277
d-28, 1.88736437
d-28, 1.63817539
d-28, &
209 1.42283147
d-28, 1.23631916
d-28, 1.07526003
d-28, 9.36797928
d-29, &
210 8.18565660
d-29, 7.18152734
d-29, 6.32523238
d-29, 5.59513985
d-29, &
211 4.96614048
d-29, 4.42518826
d-29, 3.95487628
d-29, 3.54690294
d-29, &
212 3.18953930
d-29, 2.87720933
d-29, 2.60186750
d-29, 2.36011522
d-29, &
213 2.14717806
d-29, 1.95905217
d-29, 1.79287981
d-29, 1.64562262
d-29, &
214 1.51489425
d-29, 1.39876064
d-29, 1.29496850
d-29, 1.18665438
d-29, &
215 1.10240474
d-29, 1.02643099
d-29, 9.57780996
d-30, 8.95465151
d-30, &
216 8.38950190
d-30, 7.87283711
d-30, 7.40136507
d-30, 6.96804279
d-30, &
219 data f_335 / 2.46882661
d-32, 1.89476632
d-31, 1.13216502
d-30, 4.89532008
d-30, &
220 1.32745970
d-29, 2.25390335
d-29, 3.00511672
d-29, 3.96035934
d-29, &
221 5.77977656
d-29, 8.58600736
d-29, 1.14083000
d-28, 1.48644411
d-28, &
222 2.15788823
d-28, 3.51628877
d-28, 6.12200698
d-28, 1.08184987
d-27, &
223 1.85590697
d-27, 2.91679107
d-27, 3.94405396
d-27, 4.63610680
d-27, &
224 5.13824456
d-27, 5.66602209
d-27, 6.30009232
d-27, 7.03422868
d-27, &
225 7.77973918
d-27, 8.32371831
d-27, 8.56724316
d-27, 8.62601374
d-27, &
226 8.13308844
d-27, 6.53188216
d-27, 4.55197029
d-27, 3.57590087
d-27, &
227 3.59571707
d-27, 4.03502770
d-27, 4.54366411
d-27, 4.96914990
d-27, &
228 5.24601170
d-27, 5.39979250
d-27, 5.43023669
d-27, 5.26235042
d-27, &
229 4.91585495
d-27, 4.52628362
d-27, 4.13385020
d-27, 3.67538967
d-27, &
230 3.39939742
d-27, 3.81284533
d-27, 5.02332701
d-27, 6.19438602
d-27, &
231 6.49613071
d-27, 6.04010475
d-27, 5.24664275
d-27, 4.37225997
d-27, &
232 3.52957182
d-27, 2.76212276
d-27, 2.08473158
d-27, 1.50850518
d-27, &
233 1.04602472
d-27, 7.13091243
d-28, 5.34289645
d-28, 5.21079581
d-28, &
234 6.22246365
d-28, 6.99555864
d-28, 6.29665489
d-28, 4.45077026
d-28, &
235 2.67046793
d-28, 1.52774686
d-28, 9.18061770
d-29, 6.09116074
d-29, &
236 4.48562572
d-29, 3.59463696
d-29, 3.05820218
d-29, 2.70766652
d-29, &
237 2.46144034
d-29, 2.27758450
d-29, 2.13331183
d-29, 2.01537836
d-29, &
238 1.91566180
d-29, 1.82893912
d-29, 1.75167748
d-29, 1.68136168
d-29, &
239 1.61615595
d-29, 1.55481846
d-29, 1.49643236
d-29, 1.44046656
d-29, &
240 1.38657085
d-29, 1.33459068
d-29, 1.28447380
d-29, 1.23615682
d-29, &
241 1.18963296
d-29, 1.14478976
d-29, 1.10146637
d-29, 1.04039479
d-29, &
242 9.98611410
d-30, 9.58205147
d-30, 9.19202009
d-30, 8.81551313
d-30, &
243 8.45252127
d-30, 8.10224764
d-30, 7.76469090
d-30, 7.43954323
d-30, &
249 data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
250 4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
251 5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
252 5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
253 6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
254 6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
255 7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
256 7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
259 data f_1354 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
260 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
261 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
262 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
263 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
264 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.09503647
d-39, &
265 5.47214550
d-36, 2.42433983
d-33, 2.75295034
d-31, 1.21929718
d-29, &
266 2.48392125
d-28, 2.33268145
d-27, 8.68623633
d-27, 1.00166284
d-26, &
267 3.63126633
d-27, 7.45174807
d-28, 1.38224064
d-28, 2.69270994
d-29, &
268 5.53314977
d-30, 1.15313092
d-30, 2.34195788
d-31, 4.48242942
d-32, &
274 data t_eis1 / 1.99526231
d+05, 2.23872114
d+05, 2.51188643
d+05, 2.81838293
d+05, &
275 3.16227766
d+05, 3.54813389
d+05, 3.98107171
d+05, 4.46683592
d+05, &
276 5.01187234
d+05, 5.62341325
d+05, 6.30957344
d+05, 7.07945784
d+05, &
277 7.94328235
d+05, 8.91250938
d+05, 1.00000000
d+06, 1.12201845
d+06, &
278 1.25892541
d+06, 1.41253754
d+06, 1.58489319
d+06, 1.77827941
d+06, &
279 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
280 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
281 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
282 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
283 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
284 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
285 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
286 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
287 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
288 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08 /
290 data t_eis2 / 1.99526231
d+06, 2.23872114
d+06, 2.51188643
d+06, 2.81838293
d+06, &
291 3.16227766
d+06, 3.54813389
d+06, 3.98107171
d+06, 4.46683592
d+06, &
292 5.01187234
d+06, 5.62341325
d+06, 6.30957344
d+06, 7.07945784
d+06, &
293 7.94328235
d+06, 8.91250938
d+06, 1.00000000
d+07, 1.12201845
d+07, &
294 1.25892541
d+07, 1.41253754
d+07, 1.58489319
d+07, 1.77827941
d+07, &
295 1.99526231
d+07, 2.23872114
d+07, 2.51188643
d+07, 2.81838293
d+07, &
296 3.16227766
d+07, 3.54813389
d+07, 3.98107171
d+07, 4.46683592
d+07, &
297 5.01187234
d+07, 5.62341325
d+07, 6.30957344
d+07, 7.07945784
d+07, &
298 7.94328235
d+07, 8.91250938
d+07, 1.00000000
d+08, 1.12201845
d+08, &
299 1.25892541
d+08, 1.41253754
d+08, 1.58489319
d+08, 1.77827941
d+08, &
300 1.99526231
d+08, 2.23872114
d+08, 2.51188643
d+08, 2.81838293
d+08, &
301 3.16227766
d+08, 3.54813389
d+08, 3.98107171
d+08, 4.46683592
d+08, &
302 5.01187234
d+08, 5.62341325
d+08, 6.30957344
d+08, 7.07945784
d+08, &
303 7.94328235
d+08, 8.91250938
d+08, 1.00000000
d+09, 1.12201845
d+09, &
304 1.25892541
d+09, 1.41253754
d+09, 1.58489319
d+09, 1.77827941
d+09 /
306 data f_263 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
307 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, &
308 0.00000000
d+00, 4.46454917
d-45, 3.26774829
d-42, 1.25292566
d-39, &
309 2.66922338
d-37, 3.28497742
d-35, 2.38677554
d-33, 1.03937729
d-31, &
310 2.75075687
d-30, 4.47961733
d-29, 4.46729177
d-28, 2.64862689
d-27, &
311 8.90863800
d-27, 1.72437548
d-26, 2.22217752
d-26, 2.27999477
d-26, &
312 2.08264363
d-26, 1.78226687
d-26, 1.45821699
d-26, 1.14675379
d-26, &
313 8.63082492
d-27, 6.15925429
d-27, 4.11252514
d-27, 2.51530564
d-27, &
314 1.37090986
d-27, 6.42443134
d-28, 2.48392636
d-28, 7.59187874
d-29, &
315 1.77852938
d-29, 3.23945221
d-30, 4.90533903
d-31, 6.75458158
d-32, &
316 9.06878868
d-33, 1.23927474
d-33, 1.75769395
d-34, 2.60710914
d-35, &
317 4.04318030
d-36, 6.53500581
d-37, 1.09365022
d-37, 1.88383322
d-38, &
318 3.31425233
d-39, 5.90964084
d-40, 1.06147549
d-40, 1.90706170
d-41, &
319 3.41331584
d-42, 6.07310718
d-43, 1.07364738
d-43, 1.89085498
d-44, &
320 3.32598922
d-45, 5.87125640
d-46, 0.00000000
d+00, 0.00000000
d+00 /
322 data f_264 / 0.00000000
d+00, 2.81670057
d-46, 1.28007268
d-43, 2.54586603
d-41, &
323 2.67887256
d-39, 1.68413285
d-37, 6.85702304
d-36, 1.91797284
d-34, &
324 3.84675839
d-33, 5.69939170
d-32, 6.36224608
d-31, 5.39176489
d-30, &
325 3.45478458
d-29, 1.64848693
d-28, 5.71476364
d-28, 1.39909997
d-27, &
326 2.37743056
d-27, 2.86712530
d-27, 2.65206348
d-27, 2.07175767
d-27, &
327 1.47866767
d-27, 1.01087374
d-27, 6.79605811
d-28, 4.54746770
d-28, &
328 3.04351751
d-28, 2.03639149
d-28, 1.35940991
d-28, 9.01451939
d-29, &
329 5.91289972
d-29, 3.81821178
d-29, 2.41434696
d-29, 1.48871220
d-29, &
330 8.93362094
d-30, 5.21097445
d-30, 2.95964719
d-30, 1.64278748
d-30, &
331 8.95571660
d-31, 4.82096011
d-31, 2.57390991
d-31, 1.36821781
d-31, &
332 7.27136350
d-32, 3.87019426
d-32, 2.06883430
d-32, 1.11228884
d-32, &
333 6.01883313
d-33, 3.27790676
d-33, 1.79805012
d-33, 9.93085346
d-34, &
334 5.52139556
d-34, 3.08881387
d-34, 1.73890315
d-34, 9.84434964
d-35, &
335 5.60603378
d-35, 3.20626492
d-35, 1.84111068
d-35, 0.00000000
d+00, &
336 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
338 data f_192 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 4.35772105
d-44, &
339 1.26162319
d-41, 1.97471205
d-39, 1.83409019
d-37, 1.08206288
d-35, &
340 4.27914363
d-34, 1.17943846
d-32, 2.32565755
d-31, 3.33087991
d-30, &
341 3.47013260
d-29, 2.60375866
d-28, 1.37737127
d-27, 5.01053913
d-27, &
342 1.23479810
d-26, 2.11310542
d-26, 2.71831513
d-26, 2.89851163
d-26, &
343 2.77312376
d-26, 2.50025229
d-26, 2.18323661
d-26, 1.86980322
d-26, &
344 1.58035034
d-26, 1.31985651
d-26, 1.08733133
d-26, 8.81804906
d-27, &
345 7.00417973
d-27, 5.43356567
d-27, 4.09857884
d-27, 2.99651764
d-27, &
346 2.11902962
d-27, 1.45014127
d-27, 9.62291023
d-28, 6.21548647
d-28, &
347 3.92807578
d-28, 2.44230375
d-28, 1.50167782
d-28, 9.17611405
d-29, &
348 5.58707641
d-29, 3.40570915
d-29, 2.08030862
d-29, 1.27588676
d-29, &
349 7.86535588
d-30, 4.87646151
d-30, 3.03888897
d-30, 1.90578649
d-30, &
350 1.20195947
d-30, 7.61955060
d-31, 4.85602199
d-31, 3.11049969
d-31, &
351 2.00087065
d-31, 1.29223740
d-31, 8.37422008
d-32, 0.00000000
d+00, &
352 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
354 data f_255 / 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 1.76014287
d-44, &
355 5.07057938
d-42, 7.90473970
d-40, 7.31852999
d-38, 4.30709255
d-36, &
356 1.70009061
d-34, 4.67925160
d-33, 9.21703546
d-32, 1.31918676
d-30, &
357 1.37393161
d-29, 1.03102379
d-28, 5.45694018
d-28, 1.98699648
d-27, &
358 4.90346776
d-27, 8.40524725
d-27, 1.08321456
d-26, 1.15714525
d-26, &
359 1.10905152
d-26, 1.00155023
d-26, 8.75799161
d-27, 7.50935839
d-27, &
360 6.35253533
d-27, 5.30919268
d-27, 4.37669455
d-27, 3.55185164
d-27, &
361 2.82347055
d-27, 2.19257595
d-27, 1.65589541
d-27, 1.21224987
d-27, &
362 8.58395132
d-28, 5.88163935
d-28, 3.90721447
d-28, 2.52593407
d-28, &
363 1.59739995
d-28, 9.93802874
d-29, 6.11343388
d-29, 3.73711135
d-29, &
364 2.27618743
d-29, 1.38793199
d-29, 8.48060787
d-30, 5.20305940
d-30, &
365 3.20867365
d-30, 1.99011277
d-30, 1.24064551
d-30, 7.78310544
d-31, &
366 4.91013681
d-31, 3.11338381
d-31, 1.98451675
d-31, 1.27135460
d-31, &
367 8.17917486
d-32, 5.28280497
d-32, 3.42357159
d-32, 0.00000000
d+00, &
368 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00, 0.00000000
d+00 /
373 integer,
intent(in) :: ixI^L, ixO^L
374 double precision,
intent(in) :: w(ixI^S,nw)
375 double precision,
intent(in) :: x(ixI^S,1:ndim)
376 double precision,
intent(out):: res(ixI^S)
384 integer,
intent(in) :: ixI^L, ixO^L
385 double precision,
intent(in) :: w(ixI^S, nw)
386 double precision,
intent(out):: val1(ixI^S), val2(ixI^S)
392 procedure(
get_subr1),
pointer,
nopass :: get_rho => null()
393 procedure(
get_subr1),
pointer,
nopass :: get_pthermal => null()
394 procedure(
get_subr1),
pointer,
nopass :: get_var_rfactor => null()
405 double precision,
allocatable :: source(:^d&)
406 double precision,
allocatable :: opacity(:^d&)
407 double precision,
allocatable :: sourcev(:^d&)
408 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
409 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
410 double precision,
allocatable :: rface2(:),theta_cos(:),phi_sin(:),phi_cos(:)
411 double precision :: box_min(1:3)=0.d0
412 double precision :: box_max(1:3)=0.d0
417 logical :: has_pixels=.false.
429 character(len=*),
intent(in) :: datatype
432 call mpistop(
"bad radiation_transfer")
448 case(
'cart',
'cart_dda')
450 case(
'spherical',
'sph_intersection')
452 case(
'sph_dda',
'spherical_dda')
462 call mpistop(
"bad emission_model")
466 call mpistop(
"tau and absorption-fraction output need thick transfer")
470 call mpistop(
"radsyn_pixel_batch must be positive")
473 call mpistop(
"radsyn_segment_batch_factor must be non-negative")
476 call mpistop(
"radsyn_segment_memory_mb must be positive")
479 call mpistop(
"radsyn_segment_comm_factor must be positive")
484 call mpistop(
"instrument_postprocess currently needs dat-resolution EUV images")
487 call mpistop(
"instrument_postprocess is not yet supported for spherical rays")
490 call mpistop(
"instrument_postprocess currently supports only EUV AIA or radio_ff images")
493 call mpistop(
"radio_ff instrument_postprocess needs radio_beam_fwhm > 0 arcsec")
501 if (datatype /=
'image_euv' .and. datatype /=
'spectrum_euv')
then
502 call mpistop(
"emission_model=euv_aia is only valid for EUV synthesis")
505 if (datatype /=
'image_whitelight')
then
506 call mpistop(
"emission_model=white_light is only valid for white-light synthesis")
509 if (datatype /=
'image_euv')
then
510 call mpistop(
"emission_model=radio_ff currently reuses EUV-image convert types")
513 call mpistop(
"emission_model=radio_ff needs radio_frequency > 0")
515 case(
'pseudo_current')
516 if (datatype /=
'image_euv')
then
517 call mpistop(
"emission_model=pseudo_current is only valid for EUV-image convert types")
520 call mpistop(
"emission_model=pseudo_current currently supports only thin transfer")
526 call mpistop(
"ray_method=cart needs Cartesian EUV slab images")
531 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
534 call mpistop(
"ray_method=spherical currently needs 3D spherical grids")
539 call mpistop(
"bad ray_method=spherical mode")
542 call mpistop(
"ray_method=spherical currently supports only EUV AIA emission")
544 if (xprobmin2<=1.
d-10 .or. xprobmax2>=dpi-1.
d-10)
then
545 call mpistop(
"ray_method=spherical does not support polar-axis crossing domains")
547 if (xprobmax3<=xprobmin3 .or. xprobmax3-xprobmin3>=2.d0*dpi-1.
d-10)
then
548 call mpistop(
"ray_method=spherical does not support phi-wrapping domains")
553 if (datatype /=
'image_euv')
then
554 call mpistop(
"thick transfer is only defined for EUV images")
559 if (.not.
slab)
call mpistop(
"cartesian thick EUV currently needs slab output")
561 call mpistop(
"thick EUV currently needs Cartesian dat_resolution output")
567 call mpistop(
"thick EUV currently needs x/y/z-aligned LOS")
574 double precision,
intent(in) :: emissivity,opacity,path_length
575 double precision,
intent(inout) :: intensity,tau
577 double precision :: dtau
579 if (path_length<=zero)
return
581 dtau=max(zero,opacity)*path_length
587 trim(emission_model)/=
'radio_ff' .and. &
592 logical,
intent(in) :: has_doppler,has_thick
595 if (has_doppler) num_outputs=num_outputs+1
596 if (has_thick .and. output_tau) num_outputs=num_outputs+1
597 if (has_thick .and. output_absorption_fraction) num_outputs=num_outputs+1
601 integer,
intent(in) :: nI1,nI2
602 double precision,
intent(in) :: EUV(nI1,nI2),unitv
603 double precision,
intent(inout) :: Dpl(nI1,nI2)
609 if (euv(ix1,ix2)/=zero)
then
610 dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
614 if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
620 integer,
intent(in) :: nI1,nI2
621 double precision,
intent(in) :: EUV(nI1,nI2),EUVthin(nI1,nI2),smallflux
622 double precision,
intent(out) :: Absorption(nI1,nI2)
623 logical,
intent(in),
optional :: cap_to_one
626 logical :: cap_absorption
629 cap_absorption=.false.
630 if (
present(cap_to_one)) cap_absorption=cap_to_one
633 if (euvthin(ix1,ix2)>smallflux)
then
634 absorption(ix1,ix2)=max(zero,(euvthin(ix1,ix2)-euv(ix1,ix2))/euvthin(ix1,ix2))
635 if (cap_absorption) absorption(ix1,ix2)=min(one,absorption(ix1,ix2))
641 subroutine pack_euv_image_outputs(nI1,nI2,EUV,wI,smallflux,has_doppler,has_thick,Dpl,Tau,EUVthin,&
643 integer,
intent(in) :: nI1,nI2
644 double precision,
intent(in) :: EUV(nI1,nI2),smallflux
645 double precision,
intent(inout) :: wI(:,:,:)
646 logical,
intent(in) :: has_doppler,has_thick
647 double precision,
intent(in),
optional :: Dpl(nI1,nI2),Tau(nI1,nI2),EUVthin(nI1,nI2)
648 logical,
intent(in),
optional :: cap_absorption
651 double precision,
allocatable :: Absorption(:,:)
656 if (has_doppler)
then
657 if (.not.
present(dpl))
call mpistop(
"Doppler output requested without Doppler image")
661 if (has_thick .and. output_tau)
then
662 if (.not.
present(tau))
call mpistop(
"tau output requested without tau image")
666 if (has_thick .and. output_absorption_fraction)
then
667 if (.not.
present(euvthin))
call mpistop(
"absorption output requested without thin image")
668 allocate(absorption(ni1,ni2))
671 wi(:,:,iw)=absorption(:,:)
672 deallocate(absorption)
677 integer,
intent(out) :: pixel_batch_target,segment_batch_target,segment_comm_target
679 pixel_batch_target=max(1,radsyn_pixel_batch)
680 if (radsyn_segment_batch_factor>0)
then
681 segment_batch_target=max(128,radsyn_segment_batch_factor*pixel_batch_target)
683 segment_batch_target=max(128,int(min(dble(huge(segment_batch_target)),&
684 max(128.d0,radsyn_segment_memory_mb*1048576.d0/256.d0))))
686 segment_comm_target=max(128,radsyn_segment_comm_factor*pixel_batch_target)
690 double precision,
intent(in) :: tau
700 double precision,
intent(in) :: argument
702 if (argument<-700.d0)
then
704 else if (argument>700.d0)
then
712 double precision,
intent(in) :: exponent
714 if (exponent>300.d0)
then
716 else if (exponent<-300.d0)
then
724 double precision,
intent(in) :: temperature
725 integer,
intent(in) :: n_table
726 double precision,
intent(in) :: t_table(n_table),f_table(n_table)
727 logical,
intent(in) :: log_temperature,log_response
729 integer :: ilo,ihi,imid
730 double precision :: temp_lookup,response_lookup,flo,fhi
733 if (temperature<=zero)
return
734 if (log_temperature)
then
735 temp_lookup=log10(temperature)
737 temp_lookup=temperature
739 if (temp_lookup<t_table(1) .or. temp_lookup>t_table(n_table))
return
740 if (temp_lookup==t_table(n_table))
then
741 if (log_response)
then
742 response_lookup=log10(max(f_table(n_table),1.d-99))
744 response_lookup=f_table(n_table)
751 if (temp_lookup>=t_table(imid))
then
757 if (log_response)
then
758 flo=log10(max(f_table(ilo),1.d-99))
759 fhi=log10(max(f_table(ilo+1),1.d-99))
764 response_lookup=flo*(temp_lookup-t_table(ilo+1))/(t_table(ilo)-t_table(ilo+1))+&
765 fhi*(temp_lookup-t_table(ilo))/(t_table(ilo+1)-t_table(ilo))
768 if (log_response)
then
777 integer,
intent(in) :: ixI^L, ixO^L, n_table
778 double precision,
intent(in) :: Te(ixI^S),t_table(n_table),f_table(n_table)
779 double precision,
intent(inout) :: flux(ixI^S)
780 logical,
intent(in) :: log_temperature,log_response
783 double precision :: GT
785 {
do ix^db=ixomin^db,ixomax^db\}
787 flux(ix^d)=flux(ix^d)*gt
788 if (flux(ix^d)<zero) flux(ix^d)=zero
797 integer,
intent(in) :: wl
798 integer,
intent(in) :: ixI^L, ixO^L
799 double precision,
intent(in) :: x(ixI^S,1:ndim)
800 double precision,
intent(in) :: w(ixI^S,1:nw)
802 double precision,
intent(out) :: kappa(ixI^S)
805 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
806 double precision :: wave_ratio,s_H1,s_He1,s_He2,Pe
807 double precision :: log_H21,log_He21,log_He32,log_He321,logScaleHe,w_H21
808 double precision :: term0,term1,term2,denHe,i0,j1,j2,be
809 double precision :: N_H,N_H1,N_He1,N_He2
810 double precision,
parameter :: Xe_H21=13.6d0, xe_he21=24.587d0, xe_he32=54.416d0
811 double precision,
parameter :: rHe=0.1d0
812 double precision,
parameter :: sigma_H1=5.16d-20, sigma_he1=9.25d-19, sigma_he2=7.17d-19
814 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
815 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
816 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
817 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
819 double precision :: nH_dummy(ixI^S)
820 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
823 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
825 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
828 wave_ratio=dble(wl)/171.d0
832 if (wl<=912) s_h1=wave_ratio**3*sigma_h1
833 if (wl<=504) s_he1=wave_ratio**2*sigma_he1
834 if (wl<=228) s_he2=wave_ratio**2.75d0*sigma_he2
837 {
do ix^db=ixomin^db,ixomax^db\}
838 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
839 pe=ne(ix^d)*kb_cgs*te(ix^d)
841 log_h21=2.5d0*log10(te(ix^d))-5040.d0*xe_h21/te(ix^d)-log10(pe)-0.48d0
842 log_he21=log10(4.d0)+2.5d0*log10(te(ix^d))-5040.d0*xe_he21/te(ix^d)-log10(pe)-0.48d0
843 log_he32=2.5d0*log10(te(ix^d))-5040.d0*xe_he32/te(ix^d)-log10(pe)-0.48d0
845 i0=w_h21/(1.d0+w_h21)
846 log_he321=log_he21+log_he32
847 logscalehe=max(zero,log_he21,log_he321)
851 denhe=term0+term1+term2
859 be=i0+rhe*(j1+2.d0*j2)
860 if (be>smalldouble)
then
863 n_he1=(1.d0-j1-j2)*rhe*n_h
865 kappa(ix^d)=max(zero,n_h1*s_h1+n_he1*s_he1+n_he2*s_he2)
873 integer,
intent(in) :: igrid
874 integer,
intent(in) :: ixI^L, ixO^L
875 double precision,
intent(in) :: w(ixI^S,1:nw)
876 double precision,
intent(out) :: source(ixI^S)
878 integer :: ix^D,idir,idirmin,idirmin0
879 double precision :: current(ixI^S,7-2*ndir:3)
881 if (.not.
allocated(iw_mag))
then
882 call mpistop(
"emission_model=pseudo_current needs magnetic-field variables")
887 call curlvector(w(ixi^s,iw_mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
889 current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+ps(igrid)%J0(ixo^s,idirmin0:3)
893 {
do ix^db=ixomin^db,ixomax^db\}
895 source(ix^d)=source(ix^d)+current(ix^d,idir)**2
903 integer,
intent(in) :: ixI^L, ixO^L
904 double precision,
intent(in) :: x(ixI^S,1:ndim)
905 double precision,
intent(in) :: w(ixI^S,1:nw)
907 double precision,
intent(out) :: source(ixI^S),kappa(ixI^S)
910 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
911 double precision :: nH_dummy(ixI^S),gff
913 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
914 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
915 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
916 te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
917 call eos%get_ne_nH(ixi^l,ixo^l,w,ne,nh_dummy)
919 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
921 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
926 {
do ix^db=ixomin^db,ixomax^db\}
927 if (te(ix^d)>zero .and. ne(ix^d)>zero)
then
928 if (te(ix^d)<2.d5)
then
929 gff=18.2d0+1.5d0*log(te(ix^d))-log(radio_frequency)
931 gff=24.5d0+log(te(ix^d))-log(radio_frequency)
934 kappa(ix^d)=9.78d-3*ne(ix^d)**2*gff/(radio_frequency**2*te(ix^d)**1.5d0)
935 source(ix^d)=te(ix^d)*kappa(ix^d)
940 subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
952 integer,
intent(in) :: wl
953 integer,
intent(out) :: mass
954 character(len=30),
intent(out) :: ion
955 double precision,
intent(out) :: logTe,line_center,spatial_px,spectral_px
956 double precision,
intent(out) :: sigma_PSF,width_slit
1026 line_center=1354.1d0
1028 spectral_px=12.98
d-3
1035 line_center=262.976d0
1044 line_center=263.765d0
1053 line_center=192.028d0
1062 line_center=255.113d0
1068 call mpistop(
"No information about this line")
1082 integer,
intent(in) :: wl
1083 integer,
intent(in) :: ixI^L, ixO^L
1084 double precision,
intent(in) :: x(ixI^S,1:ndim)
1085 double precision,
intent(in) :: w(ixI^S,1:nw)
1087 double precision,
intent(out) :: flux(ixI^S)
1090 double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
1092 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1093 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1094 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1098 double precision :: nH_dummy(ixI^S)
1099 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
1102 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1103 flux(ixo^s)=ne(ixo^s)**2
1105 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1106 flux(ixo^s)=ne(ixo^s)**2
1135 call mpistop(
"Unknown wavelength")
1147 integer,
intent(in) :: ixI^L,ixO^L
1148 integer,
intent(in) :: El,Eu
1149 double precision,
intent(in) :: x(ixI^S,1:ndim)
1150 double precision,
intent(in) :: w(ixI^S,nw)
1152 double precision,
intent(out) :: flux(ixI^S)
1154 integer :: ix^D,ixO^D
1156 double precision :: I0,kb,keV,dE,Ei
1157 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1158 double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
1159 double precision :: EM(ixI^S)
1165 nume=floor((eu-el)/de)
1166 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1167 call fl%get_rho(w,x,ixi^l,ixo^l,ne)
1168 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
1172 double precision :: nH_dummy(ixI^S)
1173 call eos%get_ne_nH(ixi^l, ixo^l, w, ne, nh_dummy)
1176 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1177 em(ixo^s)=(ne(ixo^s))**2*1.d6
1179 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1180 em(ixo^s)=(ne(ixo^s))**2
1182 kbt(ixo^s)=kb*te(ixo^s)/kev
1187 {
do ix^db=ixomin^db,ixomax^db\}
1188 if (kbt(ix^d)>0.01*ei)
then
1189 if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
1190 fi(ix^d)=(em(ix^d)*gff(ix^d))*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1195 flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
1197 flux(ixo^s)=flux(ixo^s)*i0
1204 double precision,
intent(in) :: xbox^L
1206 double precision,
intent(out) :: eflux
1208 double precision :: dxb^D,xb^L
1209 integer :: iigrid,igrid,j
1210 integer :: ixO^L,ixI^L,ix^D
1211 double precision :: eflux_grid,eflux_pe
1213 ^d&iximin^d=
ixglo^d;
1214 ^d&iximax^d=
ixghi^d;
1215 ^d&ixomin^d=ixmlo^d;
1216 ^d&ixomax^d=ixmhi^d;
1218 do iigrid=1,igridstail; igrid=igrids(iigrid);
1222 call get_goes_flux_grid(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,ps(igrid)%dvolume(ixi^s),xbox^l,xb^l,fl,eflux_grid)
1223 eflux_pe=eflux_pe+eflux_grid
1225 call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
1232 integer,
intent(in) :: ixI^L,ixO^L
1233 double precision,
intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
1234 double precision,
intent(in) :: w(ixI^S,nw)
1235 double precision,
intent(in) :: xbox^L,xb^L
1237 double precision,
intent(out) :: eflux_grid
1239 integer :: ix^D,ixO^D,ixb^L
1240 integer :: iE,numE,j,inbox
1241 double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
1242 double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
1243 double precision :: Ne(ixI^S),EM(ixI^S)
1244 double precision :: gff,fi,erg_SI
1248 {
if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
1250 if (inbox==ndim)
then
1252 ^d&ixbmin^d=ixomin^d;
1253 ^d&ixbmax^d=ixomax^d;
1254 {
if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/
dxlevel(^d))\}
1255 {
if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/
dxlevel(^d))+ixomin^d\}
1262 el=const_h*const_c/(8.d0*a_cgs)/kev
1263 eu=const_h*const_c/(1.d0*a_cgs)/kev
1265 nume=floor((eu-el)/de)
1266 call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
1267 call fl%get_rho(w,x,ixi^l,ixb^l,ne)
1268 call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
1272 double precision :: nH_dummy(ixI^S)
1273 call eos%get_ne_nH(ixi^l, ixb^l, w, ne, nh_dummy)
1276 ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6
1277 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3
1279 ne(ixo^s)=ne(ixo^s)*unit_numberdensity
1280 em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
1282 kbt(ixb^s)=kb*te(ixb^s)/kev
1287 {
do ix^db=ixbmin^db,ixbmax^db\}
1288 if (kbt(ix^d)>1.d-2*ei)
then
1289 if(kbt(ix^d)<ei)
then
1290 gff=(kbt(ix^d)/ei)**0.4
1294 fi=(em(ix^d)*gff)*
exp_clamped(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
1295 eflux_grid=eflux_grid+fi*de*ei
1299 eflux_grid=eflux_grid*kev*erg_si
1308 integer,
intent(in) :: qunit
1310 character(20) :: datatype
1313 character (30) :: ion
1314 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1315 double precision :: xslit,arcsec
1317 datatype=
'spectrum_euv'
1322 if (
mype==0) print *,
'###################################################'
1325 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by IRIS).'
1326 case (263,264,192,255)
1327 if (
mype==0) print *,
'Systhesizing EUV spectrum (observed by Hinode/EIS).'
1329 call mpistop(
'Wrong wavelength!')
1333 call mpistop(
'Wrong spectrum window!')
1336 if (
mype==0)
write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1337 if (
mype==0) print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1341 write(*,
'(a,f5.3,a,f5.1,a)')
' Supposed pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1342 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1344 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1346 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1348 write(*,
'(a,f8.1,a)')
' Supposed width of slit: ',wslit*725.0,
' km'
1353 print *,
'Unit of wavelength: Angstrom (0.1 nm) '
1355 write(*,
'(a,f5.3,a,f5.1,a)')
' Pixel: ',wlrsl,
' Angstrom x ',spacersl*725.0,
' km'
1356 print *,
'Unit of length: arcsec (~725 km)'
1357 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' arcsec'
1358 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit,
' arcsec'
1361 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1363 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
1365 write(*,
'(a,f8.1,a)')
' Location of slit: xI1 = ',
location_slit,
' Unit_length'
1366 write(*,
'(a,f8.1,a)')
' Width of slit: ',wslit*725.d0,
' km'
1369 if (
mype==0) print *,
'Direction of the slit: parallel to xI2 vector'
1370 if (coordinate==cartesian .or. coordinate==spherical)
then
1373 call mpistop(
"EUV spectrum synthesis: support for sperical coordinates is to be added!")
1377 if (
mype==0) print *,
'###################################################'
1383 integer,
intent(in) :: qunit
1384 character(20),
intent(in) :: datatype
1387 integer :: numWL,numXS,iwL,ixS,numWI,numS
1388 double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1389 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1390 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1391 integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1392 double precision :: qs,dxfirst,dxmid,lenstr
1394 integer :: iigrid,igrid,j,dir_loc
1395 double precision :: xbmin(1:ndim),xbmax(1:ndim)
1398 numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1399 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1400 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1401 allocate(wl(numwl),dwl(numwl))
1404 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1407 select case(direction_slit)
1409 numxs=domain_nx1*2**(refine_max_level-1)
1414 strtype=stretch_type(1)
1415 nstrb=nstretchedblocks_baselevel(1)
1416 qs=qstretch_baselevel(1)
1417 if (mype==0) print *,
'Direction of the slit: x'
1419 numxs=domain_nx2*2**(refine_max_level-1)
1424 strtype=stretch_type(2)
1425 nstrb=nstretchedblocks_baselevel(2)
1426 qs=qstretch_baselevel(2)
1427 if (mype==0) print *,
'Direction of the slit: y'
1429 numxs=domain_nx3*2**(refine_max_level-1)
1434 strtype=stretch_type(3)
1435 nstrb=nstretchedblocks_baselevel(3)
1436 qs=qstretch_baselevel(3)
1437 if (mype==0) print *,
'Direction of the slit: z'
1439 call mpistop(
'Wrong direction_slit')
1442 allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1444 allocate(wi(numwl,numxs,numwi))
1446 select case(strtype)
1448 dxs(:)=(xsmax-xsmin)/numxs
1450 xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1453 qs=qs**(one/2**(refine_max_level-1))
1454 dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1457 dxs(ixs)=dxfirst*qs**(ixs-1)
1458 xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1464 lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1465 dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1468 nstr=nstr*2**(refine_max_level-1)
1469 nuni=nuni*2**(refine_max_level-1)
1470 qs=qs**(one/2**(refine_max_level-1))
1471 dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1472 dxmid=dxmid/2**(refine_max_level-1)
1474 if(nuni .gt. 0)
then
1475 do ixs=nstr+1,nstr+nuni
1477 xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1482 dxs(ixs)=dxfirst*qs**(nstr-ixs)
1483 xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1486 do ixs=nstr+nuni+1,numxs
1487 dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1488 xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1491 call mpistop(
"unknown stretch type")
1494 if (los_phi==0 .and. los_theta==90 .and. direction_slit==2)
then
1497 else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3)
then
1500 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1)
then
1503 else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3)
then
1506 else if (los_theta==0 .and. direction_slit==1)
then
1509 else if (los_theta==0 .and. direction_slit==2)
then
1513 call mpistop(
'Wrong combination of LOS and slit direction!')
1516 if (dir_loc==1)
then
1517 if (location_slit>xprobmax1 .or. location_slit<xprobmin1)
then
1518 call mpistop(
'Wrong value for location_slit!')
1520 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: x = ',location_slit,
' Unit_length'
1521 else if (dir_loc==2)
then
1522 if (location_slit>xprobmax2 .or. location_slit<xprobmin2)
then
1523 call mpistop(
'Wrong value for location_slit!')
1525 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: y = ',location_slit,
' Unit_length'
1527 if (location_slit>xprobmax3 .or. location_slit<xprobmin3)
then
1528 call mpistop(
'Wrong value for location_slit!')
1530 if(mype==0)
write(*,
'(a,f8.1,a)')
' Location of slit: z = ',location_slit,
' Unit_length'
1535 do iigrid=1,igridstail; igrid=igrids(iigrid);
1536 ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1537 ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1538 if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc))
then
1544 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1545 mpi_sum,icomm,ierrmpi)
1548 if (spectra_rc(iwl,ixs)>smalldouble)
then
1549 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1556 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1558 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1565 integer,
intent(in) :: igrid,numWL,numXS,dir_loc
1567 double precision,
intent(in) :: wL(numWL),dwL(numWL)
1568 double precision,
intent(inout) :: spectra(numWL,numXS)
1570 integer :: direction_LOS
1571 integer :: ixO^L,ixI^L,ix^D,ixOnew
1572 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1573 double precision :: wlc,wlwd
1576 double precision :: logTe,lineCent
1577 character (30) :: ion
1578 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1580 integer :: levelg,rft,ixSmin,ixSmax,iwL
1581 double precision :: flux_pix,dL
1583 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1585 if (los_phi==0 .and. los_theta==90)
then
1587 else if (los_phi==90 .and. los_theta==90)
then
1593 ^d&ixomin^d=ixmlo^d\
1594 ^d&ixomax^d=ixmhi^d\
1595 ^d&iximin^d=ixglo^d\
1596 ^d&iximax^d=ixghi^d\
1597 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1600 if (dir_loc==1)
then
1601 do ix1=ixomin1,ixomax1
1602 if (location_slit>=(ps(igrid)%x(ix^d,1)-
half*ps(igrid)%dx(ix^d,1)) .and. &
1603 location_slit<(ps(igrid)%x(ix^d,1)+
half*ps(igrid)%dx(ix^d,1)))
then
1609 else if (dir_loc==2)
then
1610 do ix2=ixomin2,ixomax2
1611 if (location_slit>=(ps(igrid)%x(ix^d,2)-
half*ps(igrid)%dx(ix^d,2)) .and. &
1612 location_slit<(ps(igrid)%x(ix^d,2)+
half*ps(igrid)%dx(ix^d,2)))
then
1619 do ix3=ixomin3,ixomax3
1620 if (location_slit>=(ps(igrid)%x(ix^d,3)-
half*ps(igrid)%dx(ix^d,3)) .and. &
1621 location_slit<(ps(igrid)%x(ix^d,3)+
half*ps(igrid)%dx(ix^d,3)))
then
1629 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1630 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1631 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1632 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1633 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1634 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1635 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1638 levelg=ps(igrid)%level
1639 rft=2**(refine_max_level-levelg)
1641 {
do ix^d=ixomin^d,ixomax^d\}
1644 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/
const_c)
1646 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/
const_c)
1648 wlwd=sqrt(
kb_cgs*te(ix^d)*unit_temperature/(mass*
mp_cgs))
1651 select case(direction_slit)
1653 ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1654 ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1656 ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1657 ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1659 ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1660 ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1663 select case(direction_los)
1665 dl=ps(igrid)%dx(ix^d,1)*unit_length
1667 dl=ps(igrid)%dx(ix^d,2)*unit_length
1669 dl=ps(igrid)%dx(ix^d,3)*unit_length
1671 if (si_unit) dl=dl*1.d2
1674 flux_pix=flux(ix^d)*wlrsl*dl*
exp_clamped(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*
dpi)*wlwd)
1676 flux_pix=flux_pix*wslit/spacersl
1677 spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1683 deallocate(flux,v,pth,te,rho)
1689 integer,
intent(in) :: qunit
1690 character(20),
intent(in) :: datatype
1693 integer :: numWL,numXS,iwL,ixS,numWI,ix^D
1694 double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
1695 double precision,
allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1696 double precision,
allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1697 double precision :: vec_cor(1:3),xI_cor(1:2)
1698 double precision :: res,r_loc,r_max
1701 character (30) :: ion
1702 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1703 double precision :: unitv,arcsec,RHESSI_rsl,pixel
1704 integer :: iigrid,igrid,i,j,numS
1705 double precision :: xLmin,xLmax,xslit
1707 if (coordinate==spherical)
then
1715 if (coordinate==spherical)
then
1716 xsmin=-abs(xprobmax1)
1717 xsmax=abs(xprobmax1)
1720 if (ix1==1) vec_cor(1)=xprobmin1
1721 if (ix1==2) vec_cor(1)=xprobmax1
1723 if (ix2==1) vec_cor(2)=xprobmin2
1724 if (ix2==2) vec_cor(2)=xprobmax2
1726 if (ix3==1) vec_cor(3)=xprobmin3
1727 if (ix3==2) vec_cor(3)=xprobmax3
1729 r_loc=(vec_cor(1)-x_origin(1))**2
1730 r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1731 r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1733 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1736 r_max=max(r_max,r_loc)
1740 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1744 xsmin=min(xsmin,xi_cor(2))
1745 xsmax=max(xsmax,xi_cor(2))
1756 xscent=(xsmin+xsmax)/2.d0
1760 arcsec=7.25d5/unit_length
1762 arcsec=7.25d7/unit_length
1764 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1765 dxsg=spacersl*arcsec
1766 numxs=ceiling((xsmax-xscent)/dxsg)
1767 xsmin=xscent-numxs*dxsg
1768 xsmax=xscent+numxs*dxsg
1771 numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
1772 wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1773 wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1774 allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
1776 allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1778 wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1782 xs(ixs)=xsmin+dxsg*(ixs-half)
1788 do iigrid=1,igridstail; igrid=igrids(iigrid);
1790 if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1791 if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1793 if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1794 if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1796 if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1797 if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1799 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
1803 xlmin=min(xlmin,xi_cor(1))
1804 xlmax=max(xlmax,xi_cor(1))
1810 if (activate_unit_arcsec)
then
1811 xslit=location_slit*arcsec
1815 if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec)
then
1821 call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1822 mpi_sum,icomm,ierrmpi)
1825 if (spectra_rc(iwl,ixs)>smalldouble)
then
1826 wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1833 if (activate_unit_arcsec)
then
1838 call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1840 deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1846 integer,
intent(in) :: igrid,numWL,numXS
1847 double precision,
intent(in) :: wL(numWL),xS(numXS)
1848 double precision,
intent(in) :: dwLg,dxSg
1849 double precision,
intent(inout) :: spectra(numWL,numXS)
1852 integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1853 double precision,
allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1854 double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1855 double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1856 integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1857 double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1858 double precision :: xSubC(1:3),xCent(1:2)
1861 double precision :: logTe,lineCent
1862 character (30) :: ion
1863 double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1864 double precision :: sigma_wl,sigma_xs,factor
1867 arcsec=7.25d5/unit_length
1869 arcsec=7.25d7/unit_length
1871 if (activate_unit_arcsec)
then
1872 xslit=location_slit*arcsec
1877 call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1879 ^d&ixomin^d=ixmlo^d\
1880 ^d&ixomax^d=ixmhi^d\
1881 ^d&iximin^d=ixglo^d\
1882 ^d&iximax^d=ixghi^d\
1883 allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1885 call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1886 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
1887 call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1888 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1889 call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1890 te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1891 {
do ix^d=ixomin^d,ixomax^d\}
1893 vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1901 slit_width=wslit*arcsec
1902 sigma_wl=sigma_psf*dwlg
1903 sigma_xs=sigma_psf*dxsg
1904 {
do ix^d=ixomin^d,ixomax^d\}
1905 if (flux(ix^d)>smalldouble)
then
1906 xloc(1:3)=ps(igrid)%x(ix^d,1:3)
1907 dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
1911 if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
1912 xiloc(1)<=xslit+half*(slit_width+dxiloc(1)))
then
1914 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(slit_width/16.d0)));
1915 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxsg/4.d0)));
1916 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1919 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg
1920 wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1922 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg
1923 wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1925 wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1926 wlwd=wlwd*linecent/const_c
1928 {
do isubc^d=1,nsubc^d\}
1929 ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
1931 dst_slit=abs(xcent(1)-xslit)
1932 if (dst_slit<=half*slit_width)
then
1933 ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
1935 ixsmax=min(ixs+3,numxs)
1936 iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
1937 nwl=3*ceiling(wlwd/dwlg+1)
1938 iwlmin=max(1,iwl-nwl)
1939 iwlmax=min(iwl+nwl,numwl)
1941 do iwl=iwlmin,iwlmax
1942 do ixs=ixsmin,ixsmax
1943 xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1944 xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1945 xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1946 xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1947 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1948 spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
1958 deallocate(flux,v,pth,te)
1966 integer,
intent(in) :: qunit
1968 character(20) :: datatype
1971 character (30) :: ion
1972 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1973 double precision :: t0,t1
1976 datatype=
'image_euv'
1981 print *,
'###################################################'
1982 print *,
'Systhesizing EUV image'
1983 write(*,
'(a,f8.3,a)')
' Wavelength: ',linecent,
' Angstrom'
1984 print *,
'Unit of EUV flux: DN s^-1 pixel^-1'
1989 call mpistop(
'EUV dat-resolution needs Cartesian or spherical native rays')
1991 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Supposed Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0, &
1992 ' km (', spacersl,
' arcsec x ', spacersl,
' arcsec)'
1994 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
1996 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2010 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2014 write(*,
'(a,f7.1,a,f7.1,a,f5.1,a,f5.1,a)')
' Pixel: ',spacersl*725.0,
' km x ',spacersl*725.0,
' km (', &
2015 spacersl,
' arcsec x ', spacersl,
' arcsec)'
2017 print *,
'Unit of length: arcsec (~725 km)'
2020 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2022 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2026 if (coordinate==cartesian)
then
2028 '] of the simulation box is located at [X=0,Y=0] of the image'
2030 else if (coordinate==spherical)
then
2031 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2034 call mpistop(
"EUV synthesis: this coordinate is not supported!")
2039 if (
mype==0) print *,
'time comsuming: ',t1-t0,
' s'
2040 if (
mype==0) print *,
'###################################################'
2047 integer,
intent(in) :: qunit
2049 character(20) :: datatype
2050 double precision :: RHESSI_rsl
2051 double precision :: t0,t1
2054 datatype=
'image_sxr'
2059 print *,
'###################################################'
2060 print *,
'Systhesizing SXR image (observed at 1 AU).'
2065 if (coordinate/=cartesian)
call mpistop(
'SXR synthesis: only cartesian is supported for .dat resolution!')
2067 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2068 write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Supposed Pixel: ',rhessi_rsl*0.725,
' Mm x ',rhessi_rsl*0.725, &
2069 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2071 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2073 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2083 call mpistop(
'ERROR: Wrong LOS for synthesizing emission!')
2087 print *,
'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
2088 write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Pixel: ',rhessi_rsl*0.725,
' Mm x ',rhessi_rsl*0.725, &
2089 ' Mm (', rhessi_rsl,
' arcsec x ', rhessi_rsl,
' arcsec)'
2091 print *,
'Unit of length: arcsec (~725 km)'
2094 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2096 write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2100 if (coordinate==cartesian)
then
2102 '] of the simulation box is located at [X=0,Y=0] of the image'
2104 else if (coordinate==spherical)
then
2105 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2108 call mpistop(
"SXR synthesis: this coordinate is not supported!")
2113 if (
mype==0) print *,
'time comsuming:',t1-t0
2114 if (
mype==0) print *,
'###################################################'
2121 integer,
intent(in) :: qunit
2123 character(20) :: datatype
2124 double precision :: LASCO_rsl
2126 if (
mype==0) print *,
'###################################################'
2130 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C1).'
2133 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C2).'
2136 if (
mype==0) print *,
'Systhesizing white light image (observed by LASCO/C3).'
2138 call mpistop(
'Whitelight synthesis: instrument is not supported!')
2141 if (
mype==0)
write(*,
'(a,f5.1,a,f5.1,a,f5.1,a,f5.1,a)')
' Pixel: ',lasco_rsl*0.725,
' Mm x ',lasco_rsl*0.725,
' Mm (', &
2142 lasco_rsl,
' arcsec x ', lasco_rsl,
' arcsec) '
2143 if (
mype==0) print *,
'Unit of white light flux: average Sun brightness'
2145 datatype=
'image_whitelight'
2150 print *,
'Unit of length: arcsec (~725 km)'
2153 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d6,
' Mm'
2155 if (
mype==0)
write(*,
'(a,f8.1,a)')
' Unit of length: ',
unit_length/1.d8,
' Mm'
2160 if (coordinate==spherical)
then
2161 if (
mype==0)
write(*,
'(a,f6.3,f8.3,f8.3,a)')
' Mapping: R=0 of the simulation box is located at [X=0,Y=0] of the image'
2164 call mpistop(
"Whitelight synthesis: this coordinate is not supported!")
2167 if (
mype==0) print *,
'###################################################'
2172 EUV,Dpl,nOut1,nOut2,xOut1,xOut2,&
2173 dxOut1,dxOut2,wOut,numWOut,Tau,EUVthin)
2177 integer,
intent(in) :: nSrc1,nSrc2
2178 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2179 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2180 double precision,
intent(in) :: EUV(nSrc1,nSrc2),Dpl(nSrc1,nSrc2)
2181 integer,
intent(out) :: nOut1,nOut2,numWOut
2182 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2183 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2184 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),EUVthin(nSrc1,nSrc2)
2186 integer :: mass,ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw
2187 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2188 character(30) :: ion
2189 double precision :: logTe,lineCent,spaceRsl,wlRsl,sigma_PSF,wslit
2190 double precision :: arcsec,dxInst,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2191 double precision :: sigma0,xerfmin1,xerfmax1,xerfmin2,xerfmax2
2192 double precision :: factor,weightSum,weightNorm,thinVal,tauVal
2193 double precision,
allocatable :: dplNum(:,:),thinOut(:,:),tauOut(:,:),tauWeight(:,:)
2201 dxinst=spacersl*arcsec
2202 if (dxinst<=
zero)
call mpistop(
"instrument_postprocess has non-positive pixel size")
2204 xmin1=minval(xsrc1-
half*dxsrc1)
2205 xmax1=maxval(xsrc1+
half*dxsrc1)
2206 xmin2=minval(xsrc2-
half*dxsrc2)
2207 xmax2=maxval(xsrc2+
half*dxsrc2)
2208 xcent1=
half*(xmin1+xmax1)
2209 xcent2=
half*(xmin2+xmax2)
2210 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*dxinst)))
2211 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*dxinst)))
2212 xmin1=xcent1-
half*dble(nout1)*dxinst
2213 xmin2=xcent2-
half*dble(nout2)*dxinst
2215 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2217 xout1(ixp1)=xmin1+dxinst*(dble(ixp1)-
half)
2221 xout2(ixp2)=xmin2+dxinst*(dble(ixp2)-
half)
2226 if (
present(tau) .and.
output_tau) numwout=numwout+1
2228 allocate(wout(nout1,nout2,numwout),dplnum(nout1,nout2))
2231 if (
present(euvthin))
then
2232 allocate(thinout(nout1,nout2))
2236 allocate(tauout(nout1,nout2),tauweight(nout1,nout2))
2241 sigma0=sigma_psf*dxinst
2246 if (
present(euvthin)) thinval=euvthin(ixs1,ixs2)
2247 if (
present(tau)) tauval=tau(ixs1,ixs2)
2251 ixc1=floor((xsrc1(ixs1)-(xout1(1)-
half*dxinst))/dxinst)+1
2252 ixc2=floor((xsrc2(ixs2)-(xout2(1)-
half*dxinst))/dxinst)+1
2253 ixpmin1=max(1,ixc1-3)
2254 ixpmax1=min(nout1,ixc1+3)
2255 ixpmin2=max(1,ixc2-3)
2256 ixpmax2=min(nout2,ixc2+3)
2259 do ixp1=ixpmin1,ixpmax1
2260 do ixp2=ixpmin2,ixpmax2
2261 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2262 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2263 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2264 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2265 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2266 weightsum=weightsum+factor
2269 if (weightsum<=
zero) cycle
2271 do ixp1=ixpmin1,ixpmax1
2272 do ixp2=ixpmin2,ixpmax2
2273 xerfmin1=((xout1(ixp1)-
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2274 xerfmax1=((xout1(ixp1)+
half*dxinst)-xsrc1(ixs1))/(sqrt(2.d0)*sigma0)
2275 xerfmin2=((xout2(ixp2)-
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2276 xerfmax2=((xout2(ixp2)+
half*dxinst)-xsrc2(ixs2))/(sqrt(2.d0)*sigma0)
2277 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
2278 weightnorm=factor/weightsum
2279 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+euv(ixs1,ixs2)*weightnorm
2280 dplnum(ixp1,ixp2)=dplnum(ixp1,ixp2)+euv(ixs1,ixs2)*dpl(ixs1,ixs2)*weightnorm
2281 if (
present(euvthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weightnorm
2283 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weightnorm
2284 tauweight(ixp1,ixp2)=tauweight(ixp1,ixp2)+weightnorm
2294 wout(ixp1,ixp2,2)=dplnum(ixp1,ixp2)/wout(ixp1,ixp2,1)
2296 wout(ixp1,ixp2,2)=
zero
2305 if (tauweight(ixp1,ixp2)>
zero)
then
2306 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/tauweight(ixp1,ixp2)
2308 wout(ixp1,ixp2,iw)=
zero
2318 wout(ixp1,ixp2,iw)=min(
one,max(
zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2320 wout(ixp1,ixp2,iw)=
zero
2327 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,1pe12.5)') &
2328 ' instrument_postprocess EUV grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,
' dx=',dxinst
2332 if (
allocated(thinout))
deallocate(thinout)
2333 if (
allocated(tauout))
deallocate(tauout,tauweight)
2337 Bright,nOut1,nOut2,xOut1,xOut2,&
2338 dxOut1,dxOut2,wOut,numWOut,Tau,BrightThin)
2341 integer,
intent(in) :: nSrc1,nSrc2
2342 double precision,
intent(in) :: xSrc1(nSrc1),xSrc2(nSrc2)
2343 double precision,
intent(in) :: dxSrc1(nSrc1),dxSrc2(nSrc2)
2344 double precision,
intent(in) :: Bright(nSrc1,nSrc2)
2345 integer,
intent(out) :: nOut1,nOut2,numWOut
2346 double precision,
allocatable,
intent(out) :: xOut1(:),xOut2(:),dxOut1(:),dxOut2(:)
2347 double precision,
allocatable,
intent(out) :: wOut(:,:,:)
2348 double precision,
intent(in),
optional :: Tau(nSrc1,nSrc2),BrightThin(nSrc1,nSrc2)
2350 integer :: ixS1,ixS2,ixP1,ixP2,ixC1,ixC2,iw,nStencil
2351 integer :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
2352 double precision :: arcsec,beamPixel,beamSigma,xMin1,xMax1,xMin2,xMax2,xCent1,xCent2
2353 double precision :: distance1,distance2,weight,cellArea,thinVal,tauVal
2354 double precision,
allocatable :: norm(:,:),thinOut(:,:),tauOut(:,:),tauNorm(:,:)
2367 if (beamsigma<=zero .or. beampixel<=zero)
then
2368 call mpistop(
"radio beam postprocess has non-positive beam or pixel size")
2371 xmin1=minval(xsrc1-half*dxsrc1)
2372 xmax1=maxval(xsrc1+half*dxsrc1)
2373 xmin2=minval(xsrc2-half*dxsrc2)
2374 xmax2=maxval(xsrc2+half*dxsrc2)
2375 xcent1=half*(xmin1+xmax1)
2376 xcent2=half*(xmin2+xmax2)
2377 nout1=16*max(1,ceiling((xmax1-xmin1)/(16.d0*beampixel)))
2378 nout2=16*max(1,ceiling((xmax2-xmin2)/(16.d0*beampixel)))
2379 xmin1=xcent1-half*dble(nout1)*beampixel
2380 xmin2=xcent2-half*dble(nout2)*beampixel
2382 allocate(xout1(nout1),xout2(nout2),dxout1(nout1),dxout2(nout2))
2384 xout1(ixp1)=xmin1+beampixel*(dble(ixp1)-half)
2385 dxout1(ixp1)=beampixel
2388 xout2(ixp2)=xmin2+beampixel*(dble(ixp2)-half)
2389 dxout2(ixp2)=beampixel
2393 if (
present(tau) .and.
output_tau) numwout=numwout+1
2395 allocate(wout(nout1,nout2,numwout),norm(nout1,nout2))
2398 if (
present(brightthin))
then
2399 allocate(thinout(nout1,nout2))
2403 allocate(tauout(nout1,nout2),taunorm(nout1,nout2))
2408 nstencil=max(3,ceiling(4.d0*beamsigma/beampixel)+1)
2413 if (
present(brightthin)) thinval=brightthin(ixs1,ixs2)
2414 if (
present(tau)) tauval=tau(ixs1,ixs2)
2415 if (abs(bright(ixs1,ixs2))<=smalldouble .and. abs(thinval)<=smalldouble .and. &
2416 abs(tauval)<=smalldouble) cycle
2418 ixc1=floor((xsrc1(ixs1)-(xout1(1)-half*beampixel))/beampixel)+1
2419 ixc2=floor((xsrc2(ixs2)-(xout2(1)-half*beampixel))/beampixel)+1
2420 ixpmin1=max(1,ixc1-nstencil)
2421 ixpmax1=min(nout1,ixc1+nstencil)
2422 ixpmin2=max(1,ixc2-nstencil)
2423 ixpmax2=min(nout2,ixc2+nstencil)
2424 cellarea=max(smalldouble,dxsrc1(ixs1)*dxsrc2(ixs2))
2426 do ixp1=ixpmin1,ixpmax1
2427 distance1=xout1(ixp1)-xsrc1(ixs1)
2428 do ixp2=ixpmin2,ixpmax2
2429 distance2=xout2(ixp2)-xsrc2(ixs2)
2430 weight=
exp_clamped(-half*(distance1**2+distance2**2)/beamsigma**2)*cellarea
2431 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)+bright(ixs1,ixs2)*weight
2432 norm(ixp1,ixp2)=norm(ixp1,ixp2)+weight
2433 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)+thinval*weight
2435 tauout(ixp1,ixp2)=tauout(ixp1,ixp2)+tauval*weight
2436 taunorm(ixp1,ixp2)=taunorm(ixp1,ixp2)+weight
2445 if (norm(ixp1,ixp2)>zero)
then
2446 wout(ixp1,ixp2,1)=wout(ixp1,ixp2,1)/norm(ixp1,ixp2)
2447 if (
present(brightthin)) thinout(ixp1,ixp2)=thinout(ixp1,ixp2)/norm(ixp1,ixp2)
2449 wout(ixp1,ixp2,1)=zero
2450 if (
present(brightthin)) thinout(ixp1,ixp2)=zero
2460 if (taunorm(ixp1,ixp2)>zero)
then
2461 wout(ixp1,ixp2,iw)=tauout(ixp1,ixp2)/taunorm(ixp1,ixp2)
2463 wout(ixp1,ixp2,iw)=zero
2472 if (thinout(ixp1,ixp2)>smalldouble)
then
2473 wout(ixp1,ixp2,iw)=min(one,max(zero,(thinout(ixp1,ixp2)-wout(ixp1,ixp2,1))/thinout(ixp1,ixp2)))
2475 wout(ixp1,ixp2,iw)=zero
2482 write(*,
'(a,2(i8,1x),a,2(i8,1x),a,2(1pe12.5,1x))') &
2483 ' radio_beam_postprocess grid src/out: ',nsrc1,nsrc2,
' -> ',nout1,nout2,&
2488 if (
allocated(thinout))
deallocate(thinout)
2489 if (
allocated(tauout))
deallocate(tauout,taunorm)
2498 integer,
intent(in) :: qunit
2499 character(20),
intent(in) :: datatype
2502 double precision :: dx^D
2503 integer :: numX^D,ix^D
2504 double precision,
allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
2505 double precision,
allocatable :: EUVthin(:,:),Tau(:,:)
2506 double precision,
allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
2507 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
2508 integer :: numXI1,numXI2,numSI,numWI,iw
2509 double precision :: xI^L
2510 integer :: iigrid,igrid,i,j
2511 double precision,
allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
2512 double precision,
allocatable :: xIP1(:),xIP2(:),dxIP1(:),dxIP2(:),wIP(:,:,:)
2513 integer :: nXIF1,nXIF2
2514 integer :: nXIP1,nXIP2,numWIP
2515 double precision :: xIF^L
2516 double precision :: vec_cor(1:3),xI_cor(1:2),dxDDA,xIcent1,xIcent2
2518 double precision :: unitv,arcsec,RHESSI_rsl
2519 integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
2520 double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
2521 logical :: has_doppler_output,has_thick_output
2534 xicent1=
half*(xifmin1+xifmax1)
2535 xicent2=
half*(xifmin2+xifmax2)
2536 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2537 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2538 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2539 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2540 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2541 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2552 if (
mype==0)
write(*,
'(a,1pe12.5,a,2(i8,1x))') &
2553 ' spherical native dat-resolution image-plane dx: ',dxdda,
' n=',nxif1,nxif2
2558 if (ix1==1) vec_cor(1)=xprobmin1
2559 if (ix1==2) vec_cor(1)=xprobmax1
2561 if (ix2==1) vec_cor(2)=xprobmin2
2562 if (ix2==2) vec_cor(2)=xprobmax2
2564 if (ix3==1) vec_cor(3)=xprobmin3
2565 if (ix3==2) vec_cor(3)=xprobmax3
2567 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
2573 xifmin1=min(xifmin1,xi_cor(1))
2574 xifmax1=max(xifmax1,xi_cor(1))
2575 xifmin2=min(xifmin2,xi_cor(2))
2576 xifmax2=max(xifmax2,xi_cor(2))
2581 dxdda=min((xprobmax1-xprobmin1)/dble(numx1),&
2582 (xprobmax2-xprobmin2)/dble(numx2),&
2583 (xprobmax3-xprobmin3)/dble(numx3))
2584 xicent1=
half*(xifmin1+xifmax1)
2585 xicent2=
half*(xifmin2+xifmax2)
2586 nxif1=max(1,ceiling((xifmax1-xifmin1)/dxdda))
2587 nxif2=max(1,ceiling((xifmax2-xifmin2)/dxdda))
2588 xifmin1=xicent1-
half*dble(nxif1)*dxdda
2589 xifmax1=xicent1+
half*dble(nxif1)*dxdda
2590 xifmin2=xicent2-
half*dble(nxif2)*dxdda
2591 xifmax2=xicent2+
half*dble(nxif2)*dxdda
2619 if (
mype==0)
write(*,
'(a)')
' LOS vector: [-1.00 0.00 0.00]'
2620 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 0.00 1.00 0.00]'
2621 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2639 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 -1.00 0.00]'
2640 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [-1.00 0.00 0.00]'
2641 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 0.00 1.00]'
2659 if (
mype==0)
write(*,
'(a)')
' LOS vector: [ 0.00 0.00 -1.00]'
2660 if (
mype==0)
write(*,
'(a)')
' xI1 vector: [ 1.00 0.00 0.00]'
2661 if (
mype==0)
write(*,
'(a)')
' xI2 vector: [ 0.00 1.00 0.00]'
2663 allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
2666 select case(strtype1)
2668 dxif1(:)=(xifmax1-xifmin1)/nxif1
2670 xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-
half)
2674 dxfirst1=(xifmax1-xifmin1)*(
one-qs1)/(
one-qs1**nxif1)
2677 dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2678 xif1(ix1)=dxif1(1)/(
one-qs1)*(
one-qs1**(ix1-1))+
half*dxif1(ix1)
2683 nuni1=nbb1-nstrb1*bnx1
2684 lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(
one-qs1)/(
one-qs1**nstr1))
2685 dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(
one-qs1)*(
one-qs1**nstr1))
2691 dxfirst1=lenstr1*(
one-qs1)/(
one-qs1**nstr1)
2694 if(nuni1 .gt. 0)
then
2695 do ix1=nstr1+1,nstr1+nuni1
2697 xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2702 dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2703 xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*
half-dxfirst1*(
one-qs1**(nstr1-ix1))/(
one-qs1)
2706 do ix1=nstr1+nuni1+1,nxif1
2707 dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2708 xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*
half+dxfirst1*(
one-qs1**(ix1-nstr1-nuni1-1))/(
one-qs1)
2711 call mpistop(
"unknown stretch type")
2714 select case(strtype2)
2716 dxif2(:)=(xifmax2-xifmin2)/nxif2
2718 xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-
half)
2722 dxfirst2=(xifmax2-xifmin2)*(
one-qs2)/(
one-qs2**nxif2)
2725 dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2726 xif2(ix2)=dxif2(1)/(
one-qs2)*(
one-qs2**(ix2-1))+
half*dxif2(ix2)
2731 nuni2=nbb2-nstrb2*bnx2
2732 lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(
one-qs2)/(
one-qs2**nstr2))
2733 dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(
one-qs2)*(
one-qs2**nstr2))
2739 dxfirst2=lenstr2*(
one-qs2)/(
one-qs2**nstr2)
2742 if(nuni2 .gt. 0)
then
2743 do ix2=nstr2+1,nstr2+nuni2
2745 xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2750 dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2751 xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*
half-dxfirst2*(
one-qs2**(nstr2-ix2))/(
one-qs2)
2754 do ix2=nstr2+nuni2+1,nxif2
2755 dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2756 xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*
half+dxfirst2*(
one-qs2**(ix2-nstr2-nuni2-1))/(
one-qs2)
2759 call mpistop(
"unknown stretch type")
2763 if (datatype==
'image_euv')
then
2772 allocate(wi(nxif1,nxif2,numwi))
2773 allocate(euv(nxif1,nxif2),dpl(nxif1,nxif2))
2775 allocate(euvthin(nxif1,nxif2),tau(nxif1,nxif2))
2784 if (has_doppler_output)
then
2789 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
2798 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2803 do iigrid=1,igridstail; igrid=igrids(iigrid);
2807 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2809 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2812 if (has_doppler_output)
then
2816 deallocate(euvs,dpls)
2818 if (has_thick_output)
then
2819 if (has_doppler_output)
then
2821 has_thick_output,dpl=dpl,tau=tau,euvthin=euvthin)
2824 has_thick_output,tau=tau,euvthin=euvthin)
2826 else if (has_doppler_output)
then
2828 has_thick_output,dpl=dpl)
2838 euv,nxip1,nxip2,xip1,xip2,&
2839 dxip1,dxip2,wip,numwip,tau=tau,brightthin=euvthin)
2842 euv,nxip1,nxip2,xip1,xip2,&
2843 dxip1,dxip2,wip,numwip)
2847 euv,dpl,nxip1,nxip2,xip1,xip2,&
2848 dxip1,dxip2,wip,numwip,tau=tau,euvthin=euvthin)
2851 euv,dpl,nxip1,nxip2,xip1,xip2,&
2852 dxip1,dxip2,wip,numwip)
2854 call output_data(qunit,xip1,xip2,dxip1,dxip2,wip,nxip1,nxip2,numwip,datatype)
2855 deallocate(xip1,xip2,dxip1,dxip2,wip)
2857 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2860 deallocate(wi,euv,dpl,euvthin,tau)
2862 deallocate(wi,euv,dpl)
2867 if (datatype==
'image_sxr')
then
2875 allocate(wi(nxif1,nxif2,numwi))
2876 allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
2879 do iigrid=1,igridstail; igrid=igrids(iigrid);
2883 call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
2886 sxr=sxr*(rhessi_rsl*arcsec)**2
2894 call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2895 deallocate(wi,sxr,sxrs)
2898 deallocate(xif1,xif2,dxif1,dxif2)
2905 integer,
intent(in) :: igrid,nXIF1,nXIF2
2906 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2907 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2909 double precision,
intent(out) :: SXR(nXIF1,nXIF2)
2911 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2912 double precision :: xb^L,xd^D
2913 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
2914 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2915 double precision,
allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2916 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2917 double precision :: SXRt,xc^L,xg^L,r2,area_1AU
2918 integer :: ixP^L,ixP^D
2919 integer :: direction_LOS
2929 ^d&ixomin^d=ixmlo^d\
2930 ^d&ixomax^d=ixmhi^d\
2931 ^d&iximin^d=
ixglo^d\
2932 ^d&iximax^d=
ixghi^d\
2936 allocate(flux(ixi^s))
2937 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2938 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2939 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2940 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2945 levelg=ps(igrid)%level
2949 select case(direction_los)
2960 allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2966 select case(direction_los)
2968 do ix2=ixomin2,ixomax2
2969 ixgmin1=(ix2-1)*rft+1
2971 do ix3=ixomin3,ixomax3
2972 ixgmin2=(ix3-1)*rft+1
2975 do ix1=ixomin1,ixomax1
2978 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2982 do ix3=ixomin3,ixomax3
2983 ixgmin1=(ix3-1)*rft+1
2985 do ix1=ixomin1,ixomax1
2986 ixgmin2=(ix1-1)*rft+1
2989 do ix2=ixomin2,ixomax2
2992 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2996 do ix1=ixomin1,ixomax1
2997 ixgmin1=(ix1-1)*rft+1
2999 do ix2=ixomin2,ixomax2
3000 ixgmin2=(ix2-1)*rft+1
3003 do ix3=ixomin3,ixomax3
3006 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
3016 select case(direction_los)
3018 ixgmin1=(ixomin2-1)*rft+1
3020 ixgmin2=(ixomin3-1)*rft+1
3023 ixgmin1=(ixomin3-1)*rft+1
3025 ixgmin2=(ixomin1-1)*rft+1
3028 ixgmin1=(ixomin1-1)*rft+1
3030 ixgmin2=(ixomin2-1)*rft+1
3034 select case(direction_los)
3036 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3037 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3038 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3039 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3041 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3042 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3043 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3044 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3046 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3047 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3048 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3049 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3051 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3052 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3053 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3054 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3055 sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3056 sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3058 deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
3065 integer,
intent(in) :: igrid,nXIF1,nXIF2
3066 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3067 double precision,
intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
3069 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3071 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
3072 double precision :: xb^L,xd^D
3073 double precision,
allocatable :: flux(:^D&),v(:^D&),rho(:^D&),opacity(:^D&)
3074 double precision,
allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
3075 double precision,
allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
3076 integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
3077 double precision :: EUVt,Fvt,xc^L,xg^L,r2
3078 integer :: ixP^L,ixP^D
3079 integer :: direction_LOS
3089 ^d&ixomin^d=ixmlo^d\
3090 ^d&ixomax^d=ixmhi^d\
3091 ^d&iximin^d=
ixglo^d\
3092 ^d&iximax^d=
ixghi^d\
3096 allocate(flux(ixi^s),v(ixi^s),rho(ixi^s),opacity(ixi^s))
3097 allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
3098 dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
3099 dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
3100 dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
3113 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3114 v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
3119 levelg=ps(igrid)%level
3123 select case(direction_los)
3134 allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
3141 select case(direction_los)
3143 do ix2=ixomin2,ixomax2
3144 ixgmin1=(ix2-1)*rft+1
3146 do ix3=ixomin3,ixomax3
3147 ixgmin2=(ix3-1)*rft+1
3151 do ix1=ixomin1,ixomax1
3155 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3156 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3160 do ix3=ixomin3,ixomax3
3161 ixgmin1=(ix3-1)*rft+1
3163 do ix1=ixomin1,ixomax1
3164 ixgmin2=(ix1-1)*rft+1
3168 do ix2=ixomin2,ixomax2
3172 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3173 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3177 do ix1=ixomin1,ixomax1
3178 ixgmin1=(ix1-1)*rft+1
3180 do ix2=ixomin2,ixomax2
3181 ixgmin2=(ix2-1)*rft+1
3185 do ix3=ixomin3,ixomax3
3189 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
3190 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
3201 select case(direction_los)
3203 ixgmin1=(ixomin2-1)*rft+1
3205 ixgmin2=(ixomin3-1)*rft+1
3208 ixgmin1=(ixomin3-1)*rft+1
3210 ixgmin2=(ixomin1-1)*rft+1
3213 ixgmin1=(ixomin1-1)*rft+1
3215 ixgmin2=(ixomin2-1)*rft+1
3219 select case(direction_los)
3221 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3222 ixpmax1=
node(pig2_,igrid)*rft*block_nx2
3223 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3224 ixpmax2=
node(pig3_,igrid)*rft*block_nx3
3226 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+1
3227 ixpmax1=
node(pig3_,igrid)*rft*block_nx3
3228 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3229 ixpmax2=
node(pig1_,igrid)*rft*block_nx1
3231 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+1
3232 ixpmax1=
node(pig1_,igrid)*rft*block_nx1
3233 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+1
3234 ixpmax2=
node(pig2_,igrid)*rft*block_nx2
3236 xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
3237 xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
3238 dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
3239 dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
3240 euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3241 euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3242 dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
3243 fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
3245 deallocate(flux,v,opacity,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
3254 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),box_min(1:3),box_max(1:3)
3255 logical,
intent(out) :: hit
3256 double precision,
intent(out) :: t_enter,t_exit
3259 double precision :: t1,t2,td
3265 if (abs(ray_dir(idir))<=smalldouble)
then
3266 if (ray_origin(idir)<box_min(idir) .or. ray_origin(idir)>box_max(idir))
then
3271 t1=(box_min(idir)-ray_origin(idir))/ray_dir(idir)
3272 t2=(box_max(idir)-ray_origin(idir))/ray_dir(idir)
3278 t_enter=max(t_enter,t1)
3279 t_exit=min(t_exit,t2)
3280 if (t_enter>=t_exit)
then
3289 integer,
intent(in) :: ixI^L, ixO^L
3290 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
3291 double precision,
allocatable,
intent(out) :: xface1(:),xface2(:),xface3(:)
3295 allocate(xface1(ixomin1:ixomax1+1),xface2(ixomin2:ixomax2+1),xface3(ixomin3:ixomax3+1))
3299 do ix1=ixomin1,ixomax1
3300 xface1(ix1)=x(ix^d,1)-half*dx(ix^d,1)
3302 xface1(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
3306 do ix2=ixomin2,ixomax2
3307 xface2(ix2)=x(ix^d,2)-half*dx(ix^d,2)
3309 xface2(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
3313 do ix3=ixomin3,ixomax3
3314 xface3(ix3)=x(ix^d,3)-half*dx(ix^d,3)
3316 xface3(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
3320 integer,
intent(in) :: imin,imax
3321 double precision,
intent(in) :: pos,faces(imin:imax+1)
3323 integer :: ilo,ihi,imid
3325 if (pos<=faces(imin))
then
3329 if (pos>=faces(imax+1))
then
3336 do while (ihi-ilo>1)
3338 if (pos>=faces(imid))
then
3344 idx=min(imax,max(imin,ilo))
3348 integer,
intent(in) :: imin,imax,idx
3349 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3350 integer,
intent(out) :: step
3351 double precision,
intent(out) :: tMax
3353 if (ray_dir_axis>zero)
then
3355 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3356 else if (ray_dir_axis<zero)
then
3358 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3366 integer,
intent(in) :: imin,imax,step
3367 double precision,
intent(in) :: ray_origin_axis,ray_dir_axis,faces(imin:imax+1)
3368 integer,
intent(inout) :: idx
3369 double precision,
intent(inout) :: tMax
3370 logical,
intent(out) :: done
3374 if (idx<imin .or. idx>imax)
then
3379 tmax=(faces(idx+1)-ray_origin_axis)/ray_dir_axis
3380 else if (step<0)
then
3381 tmax=(faces(idx)-ray_origin_axis)/ray_dir_axis
3388 ray_origin,xface1,xface2,xface3,t_enter,t_exit,EUVp,Dplp)
3389 integer,
intent(in) :: ixI^L, ixO^L
3390 double precision,
intent(in) :: source(ixI^S),sourcev(ixI^S)
3391 double precision,
intent(in) :: ray_origin(1:3)
3392 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3393 xface3(ixOmin3:ixOmax3+1)
3394 double precision,
intent(in) :: t_enter,t_exit
3395 double precision,
intent(inout) :: EUVp,Dplp
3397 integer :: ix^D,step(1:3)
3398 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay
3401 if (t_exit<=t_enter)
return
3402 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3403 pos=ray_origin+(t_enter+epsray)*
vec_los
3414 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3415 if (tnext>tnow)
then
3416 ds_cm=(tnext-tnow)*unit_length
3417 if (si_unit) ds_cm=ds_cm*1.d2
3418 euvp=euvp+source(ix^d)*ds_cm
3419 dplp=dplp+sourcev(ix^d)*ds_cm
3422 if (tnow>=t_exit-epsray)
exit
3424 if (tmax(1)<=tnow+epsray)
then
3428 if (tmax(2)<=tnow+epsray)
then
3432 if (tmax(3)<=tnow+epsray)
then
3440 double precision,
allocatable,
intent(inout) :: segments(:,:)
3441 integer,
intent(inout) :: nseg,capacity
3442 integer,
intent(in) :: pixel_id
3443 double precision,
intent(in) :: tseg,jds,kds,jvds
3445 double precision,
allocatable :: tmp(:,:)
3446 integer :: new_capacity
3448 if (capacity<=0)
then
3450 allocate(segments(5,capacity))
3451 else if (nseg>=capacity)
then
3452 new_capacity=2*capacity
3453 allocate(tmp(5,new_capacity))
3454 tmp(:,1:capacity)=segments(:,1:capacity)
3455 call move_alloc(tmp,segments)
3456 capacity=new_capacity
3460 segments(1,nseg)=dble(pixel_id)
3461 segments(2,nseg)=tseg
3462 segments(3,nseg)=jds
3463 segments(4,nseg)=kds
3464 segments(5,nseg)=jvds
3468 pixel_id,ray_origin,xface1,xface2,xface3,t_enter,t_exit,&
3469 segments,nseg,capacity)
3470 integer,
intent(in) :: ixI^L, ixO^L
3471 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S),sourcev(ixI^S)
3472 integer,
intent(in) :: pixel_id
3473 double precision,
intent(in) :: ray_origin(1:3)
3474 double precision,
intent(in) :: xface1(ixOmin1:ixOmax1+1),xface2(ixOmin2:ixOmax2+1),&
3475 xface3(ixOmin3:ixOmax3+1)
3476 double precision,
intent(in) :: t_enter,t_exit
3477 double precision,
allocatable,
intent(inout) :: segments(:,:)
3478 integer,
intent(inout) :: nseg,capacity
3480 integer :: ix^D,step(1:3)
3481 double precision :: pos(1:3),tMax(1:3),tNow,tNext,ds_cm,epsRay,tseg
3482 double precision :: jds,kds,jvds
3485 if (t_exit<=t_enter)
return
3486 epsray=max(1.d-12,1.d-10*abs(t_exit-t_enter))
3487 pos=ray_origin+(t_enter+epsray)*
vec_los
3498 tnext=min(t_exit,tmax(1),tmax(2),tmax(3))
3499 if (tnext>tnow)
then
3500 ds_cm=(tnext-tnow)*unit_length
3501 if (si_unit) ds_cm=ds_cm*1.d2
3502 jds=source(ix^d)*ds_cm
3503 kds=opacity(ix^d)*ds_cm
3504 jvds=sourcev(ix^d)*ds_cm
3505 if (jds/=zero .or. kds/=zero .or. jvds/=zero)
then
3506 tseg=half*(tnow+tnext)
3511 if (tnow>=t_exit-epsray)
exit
3513 if (tmax(1)<=tnow+epsray)
then
3517 if (tmax(2)<=tnow+epsray)
then
3521 if (tmax(3)<=tnow+epsray)
then
3529 double precision,
intent(in) :: segments(:,:)
3530 integer,
intent(inout) :: idx(:)
3531 integer,
intent(in) :: nidx
3542 double precision,
intent(in) :: segments(:,:)
3543 integer,
intent(inout) :: idx(:)
3544 integer,
intent(in) :: ilo,ihi
3551 do while (j>=ilo .and. segments(2,idx(j))>segments(2,key))
3560 double precision,
intent(in) :: segments(:,:)
3561 integer,
intent(inout) :: idx(:)
3562 integer,
intent(in) :: ilo,ihi
3565 double precision :: pivot
3567 if (ihi-ilo<=32)
then
3574 pivot=segments(2,idx((ilo+ihi)/2))
3576 do while (segments(2,idx(i))<pivot)
3579 do while (segments(2,idx(j))>pivot)
3597 integer,
intent(in) :: pixel_id
3599 owner=mod(pixel_id-1,npe)
3603 double precision,
intent(in) :: segments(:,:)
3604 integer,
intent(in) :: is,nvars
3610 if (segments(iv,is)/=segments(iv,is) .or. abs(segments(iv,is))>=1.d90)
then
3617 subroutine cart_dda_block_pixel_range(box_min,box_max,nXIF1,nXIF2,xIF1,xIF2,ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
3618 double precision,
intent(in) :: box_min(1:3),box_max(1:3)
3619 integer,
intent(in) :: nXIF1,nXIF2
3620 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3621 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
3622 logical,
intent(out) :: has_pixels
3625 double precision :: vec_cor(1:3),xI_cor(1:2)
3626 double precision :: xmin1,xmax1,xmin2,xmax2,dx1,dx2
3629 if (i1==1) vec_cor(1)=box_min(1)
3630 if (i1==2) vec_cor(1)=box_max(1)
3632 if (i2==1) vec_cor(2)=box_min(2)
3633 if (i2==2) vec_cor(2)=box_max(2)
3635 if (i3==1) vec_cor(3)=box_min(3)
3636 if (i3==2) vec_cor(3)=box_max(3)
3638 if (i1==1 .and. i2==1 .and. i3==1)
then
3644 xmin1=min(xmin1,xi_cor(1))
3645 xmax1=max(xmax1,xi_cor(1))
3646 xmin2=min(xmin2,xi_cor(2))
3647 xmax2=max(xmax2,xi_cor(2))
3654 dx1=abs(xif1(2)-xif1(1))
3656 dx1=max(one,abs(xmax1-xmin1))
3659 dx2=abs(xif2(2)-xif2(1))
3661 dx2=max(one,abs(xmax2-xmin2))
3664 ixpmin1=max(1,floor((xmin1-xif1(1))/dx1)+1-1)
3665 ixpmax1=min(nxif1,ceiling((xmax1-xif1(1))/dx1)+1+1)
3666 ixpmin2=max(1,floor((xmin2-xif2(1))/dx2)+1-1)
3667 ixpmax2=min(nxif2,ceiling((xmax2-xif2(1))/dx2)+1+1)
3668 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
3674 integer,
intent(in) :: nXIF1,nXIF2
3675 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3677 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3679 integer :: ixO^L,ixO^D,ixI^L,ix^D
3680 integer :: iigrid,igrid,ixP1,ixP2,numSI,ixPmin1,ixPmax1,ixPmin2,ixPmax2
3681 double precision :: box_min(1:3),box_max(1:3),ray_origin(1:3)
3682 double precision :: t_enter,t_exit,vlos
3683 double precision :: profile_local(2),profile_global(2)
3684 logical :: hit,has_pixels
3685 double precision,
allocatable :: source(:^D&),sourcev(:^D&),rho(:^D&),opacity(:^D&)
3686 double precision,
allocatable :: xface1(:),xface2(:),xface3(:)
3687 double precision,
allocatable :: EUVs(:,:),Dpls(:,:)
3689 allocate(euvs(nxif1,nxif2),dpls(nxif1,nxif2))
3694 do iigrid=1,igridstail; igrid=igrids(iigrid);
3695 ^d&ixomin^d=ixmlo^d\
3696 ^d&ixomax^d=ixmhi^d\
3697 ^d&iximin^d=
ixglo^d\
3698 ^d&iximax^d=
ixghi^d\
3700 box_min(1)=
rnode(rpxmin1_,igrid)
3701 box_min(2)=
rnode(rpxmin2_,igrid)
3702 box_min(3)=
rnode(rpxmin3_,igrid)
3703 box_max(1)=
rnode(rpxmax1_,igrid)
3704 box_max(2)=
rnode(rpxmax2_,igrid)
3705 box_max(3)=
rnode(rpxmax3_,igrid)
3708 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
3709 if (.not. has_pixels)
then
3710 deallocate(xface1,xface2,xface3)
3714 allocate(source(ixi^s),sourcev(ixi^s),rho(ixi^s),opacity(ixi^s))
3724 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3725 do ix1=ixomin1,ixomax1
3726 do ix2=ixomin2,ixomax2
3727 do ix3=ixomin3,ixomax3
3728 if (rho(ix^d)>smalldouble)
then
3729 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3730 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3731 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3732 sourcev(ix^d)=source(ix^d)*vlos
3738 deallocate(rho,opacity)
3740 do ixp1=ixpmin1,ixpmax1
3741 do ixp2=ixpmin2,ixpmax2
3743 profile_local(1)=profile_local(1)+one
3746 profile_local(2)=profile_local(2)+one
3748 ray_origin,xface1,xface2,xface3,t_enter,t_exit,euvs(ixp1,ixp2),dpls(ixp1,ixp2))
3753 deallocate(source,sourcev,xface1,xface2,xface3)
3757 call mpi_allreduce(euvs,euv,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3758 call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3759 call mpi_allreduce(profile_local,profile_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
3761 write(*,
'(a,2(es12.5,1x))')
' cart_dda thin profile ray_tests ray_hits: ',profile_global
3763 deallocate(euvs,dpls)
3769 integer,
intent(in) :: nXIF1,nXIF2
3770 double precision,
intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
3772 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
3773 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
3775 integer,
parameter :: nSegVars=5
3776 integer :: ixO^L,ixO^D,ixI^L,ix^D
3777 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
3778 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
3779 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
3780 integer :: nPixBatchTarget
3781 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
3782 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
3783 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
3784 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
3785 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
3786 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
3787 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
3788 double precision :: ray_origin(1:3)
3789 double precision :: t_enter,t_exit,vlos,atten
3790 double precision :: profile_local(5),profile_global(5),profile_batch(5)
3791 logical :: hit,has_pixels,batchAccepted,batchReduced
3792 double precision,
allocatable :: rho(:^D&)
3793 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
3794 double precision,
allocatable :: segments_recv_round(:,:)
3795 double precision,
allocatable :: image_reduce(:,:)
3803 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
3804 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
3805 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
3806 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
3807 allocate(cache(igridstail))
3809 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
3810 bucketfill(npixbatchtarget))
3812 do iigrid=1,igridstail; igrid=igrids(iigrid);
3813 ^d&ixomin^d=ixmlo^d\
3814 ^d&ixomax^d=ixmhi^d\
3815 ^d&iximin^d=
ixglo^d\
3816 ^d&iximax^d=
ixghi^d\
3818 cache(iigrid)%igrid=igrid
3819 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
3820 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
3821 cache(iigrid)%source=zero
3822 cache(iigrid)%opacity=zero
3823 cache(iigrid)%sourcev=zero
3826 cache(iigrid)%source,cache(iigrid)%opacity)
3828 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
3831 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
3832 do ix1=ixomin1,ixomax1
3833 do ix2=ixomin2,ixomax2
3834 do ix3=ixomin3,ixomax3
3835 if (rho(ix^d)>smalldouble)
then
3836 vlos=(ps(igrid)%w(ix^d,iw_mom(1))*
vec_los(1)+&
3837 ps(igrid)%w(ix^d,iw_mom(2))*
vec_los(2)+&
3838 ps(igrid)%w(ix^d,iw_mom(3))*
vec_los(3))/rho(ix^d)
3839 cache(iigrid)%sourcev(ix^d)=cache(iigrid)%source(ix^d)*vlos
3846 cache(iigrid)%box_min(1)=
rnode(rpxmin1_,igrid)
3847 cache(iigrid)%box_min(2)=
rnode(rpxmin2_,igrid)
3848 cache(iigrid)%box_min(3)=
rnode(rpxmin3_,igrid)
3849 cache(iigrid)%box_max(1)=
rnode(rpxmax1_,igrid)
3850 cache(iigrid)%box_max(2)=
rnode(rpxmax2_,igrid)
3851 cache(iigrid)%box_max(3)=
rnode(rpxmax3_,igrid)
3853 cache(iigrid)%xface1,cache(iigrid)%xface2,&
3854 cache(iigrid)%xface3)
3856 nxif1,nxif2,xif1,xif2,cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
3857 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
3860 npixtotal=nxif1*nxif2
3862 do while (ipixstart<=npixtotal)
3863 ipixend=min(nxif1*nxif2,ipixstart+npixbatchtarget-1)
3864 npixbatch=ipixend-ipixstart+1
3865 batchaccepted=.false.
3866 batchreduced=.false.
3868 do while (.not. batchaccepted)
3873 do iigrid=1,igridstail; igrid=igrids(iigrid);
3874 ^d&ixomin^d=ixmlo^d\
3875 ^d&ixomax^d=ixmhi^d\
3876 ^d&iximin^d=
ixglo^d\
3877 ^d&iximax^d=
ixghi^d\
3879 ixpmin1=cache(iigrid)%ixPmin1
3880 ixpmax1=cache(iigrid)%ixPmax1
3881 ixpmin2=cache(iigrid)%ixPmin2
3882 ixpmax2=cache(iigrid)%ixPmax2
3883 has_pixels=cache(iigrid)%has_pixels
3884 if (.not. has_pixels) cycle
3886 do ixp2=ixpmin2,ixpmax2
3887 ifirst=max(ipixstart,(ixp2-1)*nxif1+ixpmin1)
3888 ilast=min(ipixend,(ixp2-1)*nxif1+ixpmax1)
3889 if (ifirst>ilast) cycle
3890 do ipix=ifirst,ilast
3891 ixp1=1+mod(ipix-1,nxif1)
3893 profile_batch(1)=profile_batch(1)+one
3895 cache(iigrid)%box_max,hit,t_enter,t_exit)
3897 profile_batch(2)=profile_batch(2)+one
3900 cache(iigrid)%opacity,cache(iigrid)%sourcev,&
3901 ipix,ray_origin,cache(iigrid)%xface1,&
3902 cache(iigrid)%xface2,cache(iigrid)%xface3,&
3904 segments,nseg,capacity)
3905 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
3911 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
3912 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
3913 npixbatch=max(1,npixbatch/2)
3914 ipixend=ipixstart+npixbatch-1
3915 if (
allocated(segments))
deallocate(segments)
3918 batchaccepted=.true.
3922 profile_local=profile_local+profile_batch
3924 write(*,
'(a,3(i0,1x))')
' cart_dda thick adaptive batch: ',&
3925 ipixstart,ipixend,maxnsegbatch
3928 if (.not.
allocated(segments))
then
3930 allocate(segments(nsegvars,capacity))
3935 ownersegcounts(owner)=ownersegcounts(owner)+1
3937 sendcounts=nsegvars*ownersegcounts
3940 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
3943 allocate(segments_send(nsegvars,max(1,nseg)))
3947 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
3948 segments_send(:,isegdest)=segments(:,is)
3949 owneroffsets(owner)=owneroffsets(owner)+1
3952 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
3955 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
3957 totalcount=sum(recvcounts)
3958 totalseg=totalcount/nsegvars
3959 profile_local(4)=profile_local(4)+dble(totalcount)
3960 allocate(segments_recv(nsegvars,max(1,totalseg)))
3963 maxownersegcountlocal=maxval(ownersegcounts)
3964 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
3965 do segoffset=0,maxownersegcount-1,maxsegcommtarget
3967 roundsenddispls=senddispls
3969 if (ownersegcounts(ipe)>segoffset)
then
3970 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
3971 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
3975 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
3976 roundrecvdispls(0)=0
3978 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
3980 totalroundcount=sum(roundrecvcounts)
3981 totalroundseg=totalroundcount/nsegvars
3982 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
3984 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
3985 segments_recv_round,roundrecvcounts,roundrecvdispls,&
3988 if (totalroundseg>0)
then
3989 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
3990 recvfill=recvfill+totalroundseg
3992 deallocate(segments_recv_round)
3995 if (recvfill/=totalseg)
call mpistop(
"cart_dda thick segmented receive mismatch")
3997 if (totalseg>0)
then
3998 allocate(idx(totalseg))
3999 bucketcounts(1:npixbatch)=0
4002 ipix=nint(segments_recv(1,is))
4004 ilocal=ipix-ipixstart+1
4005 bucketcounts(ilocal)=bucketcounts(ilocal)+1
4011 do ilocal=1,npixbatch
4012 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
4014 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
4017 ipix=nint(segments_recv(1,is))
4019 ilocal=ipix-ipixstart+1
4020 idx(bucketfill(ilocal))=is
4021 bucketfill(ilocal)=bucketfill(ilocal)+1
4026 do ipix=ipixstart,ipixend
4028 ilocal=ipix-ipixstart+1
4029 nidx=bucketcounts(ilocal)
4031 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
4033 ixglobal=1+mod(ipix-1,nxif1)
4034 iyglobal=1+(ipix-1)/nxif1
4035 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
4037 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
4039 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
4040 dpl(ixglobal,iyglobal)=dpl(ixglobal,iyglobal)+atten*segments_recv(5,is)
4041 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
4048 deallocate(segments_send,segments_recv)
4049 if (
allocated(segments))
deallocate(segments)
4053 do iigrid=1,igridstail
4054 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4055 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4056 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4057 if (
allocated(cache(iigrid)%xface1))
deallocate(cache(iigrid)%xface1)
4058 if (
allocated(cache(iigrid)%xface2))
deallocate(cache(iigrid)%xface2)
4059 if (
allocated(cache(iigrid)%xface3))
deallocate(cache(iigrid)%xface3)
4062 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
4063 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
4064 bucketoffsets,bucketfill)
4065 allocate(image_reduce(nxif1,nxif2))
4066 call mpi_allreduce(euv,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4068 call mpi_allreduce(dpl,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4070 call mpi_allreduce(tau,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4072 call mpi_allreduce(euvthin,image_reduce,nxif1*nxif2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4073 euvthin=image_reduce
4074 deallocate(image_reduce)
4075 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4077 write(*,
'(a,5(es12.5,1x))') &
4078 ' cart_dda thick profile: ',profile_global
4085 integer,
intent(in) :: nXIF1,nXIF2
4087 double precision,
intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
4088 double precision,
intent(out) :: Tau(nXIF1,nXIF2),EUVthin(nXIF1,nXIF2)
4090 integer :: ixO^L,ixO^D,ixI^L,ix^D
4091 integer :: iigrid,igrid,levelg,rft,direction_LOS,nLOS,numSeg,nLayerVars,nLayerSeg
4092 integer :: ixP1,ixP2,ixL,iSub1,iSub2,relL
4093 integer :: nLosBatch,nBatch,iBatch,ixLstart,ixLend,ixLgridStart,ixLgridEnd
4094 integer :: ixPmin1,ixPmin2
4095 double precision :: ds_cm,jds,kds,jvds,atten,layerBytes,targetBytes
4096 double precision,
allocatable :: rho(:^D&)
4097 double precision,
allocatable :: layer_ds(:,:,:,:),layer_all(:,:,:,:)
4112 if (nxif1>huge(numseg)/max(1,nxif2) .or. nxif1*nxif2>huge(numseg)/nlayervars)
then
4113 call mpistop(
"thick EUV layer buffer is too large for one MPI reduction")
4115 nlayerseg=nxif1*nxif2*nlayervars
4116 targetbytes=256.d0*1024.d0*1024.d0
4117 layerbytes=dble(nlayerseg)*8.d0*2.d0
4118 nlosbatch=max(1,min(16,int(targetbytes/max(one,layerbytes))))
4119 if (nlayerseg>huge(numseg)/nlosbatch)
then
4120 call mpistop(
"thick EUV batched layer buffer is too large for one MPI reduction")
4123 allocate(cache(igridstail))
4124 do iigrid=1,igridstail; igrid=igrids(iigrid);
4125 ^d&ixomin^d=ixmlo^d\
4126 ^d&ixomax^d=ixmhi^d\
4127 ^d&iximin^d=
ixglo^d\
4128 ^d&iximax^d=
ixghi^d\
4130 cache(iigrid)%igrid=igrid
4131 levelg=ps(igrid)%level
4133 cache(iigrid)%level=levelg
4134 cache(iigrid)%rft=rft
4136 select case(direction_los)
4138 cache(iigrid)%los_min=(
node(pig1_,igrid)-1)*rft*block_nx1+1
4139 cache(iigrid)%los_max=
node(pig1_,igrid)*rft*block_nx1
4141 cache(iigrid)%los_min=(
node(pig2_,igrid)-1)*rft*block_nx2+1
4142 cache(iigrid)%los_max=
node(pig2_,igrid)*rft*block_nx2
4144 cache(iigrid)%los_min=(
node(pig3_,igrid)-1)*rft*block_nx3+1
4145 cache(iigrid)%los_max=
node(pig3_,igrid)*rft*block_nx3
4148 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s),&
4149 cache(iigrid)%sourcev(ixi^s),rho(ixi^s))
4150 cache(iigrid)%source=zero
4151 cache(iigrid)%opacity=zero
4152 cache(iigrid)%sourcev=zero
4155 cache(iigrid)%source,cache(iigrid)%opacity)
4157 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
4160 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
4161 cache(iigrid)%sourcev(ixo^s)=cache(iigrid)%source(ixo^s)*&
4162 (-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s))
4168 allocate(layer_ds(nxif1,nxif2,nlayervars,nlosbatch),layer_all(nxif1,nxif2,nlayervars,nlosbatch))
4174 do ixlstart=1,nlos,nlosbatch
4175 ixlend=min(nlos,ixlstart+nlosbatch-1)
4176 nbatch=ixlend-ixlstart+1
4177 layer_ds(:,:,:,1:nbatch)=zero
4179 do iigrid=1,igridstail
4180 ixlgridstart=max(ixlstart,cache(iigrid)%los_min)
4181 ixlgridend=min(ixlend,cache(iigrid)%los_max)
4182 if (ixlgridstart>ixlgridend) cycle
4183 igrid=cache(iigrid)%igrid
4184 rft=cache(iigrid)%rft
4185 ^d&ixomin^d=ixmlo^d\
4186 ^d&ixomax^d=ixmhi^d\
4187 ^d&iximin^d=
ixglo^d\
4188 ^d&iximax^d=
ixghi^d\
4190 do ixl=ixlgridstart,ixlgridend
4191 ibatch=ixl-ixlstart+1
4192 rell=ixl-cache(iigrid)%los_min
4194 select case(direction_los)
4196 ix1=ixomin1+rell/rft
4197 do ix2=ixomin2,ixomax2
4198 ixpmin1=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4199 do ix3=ixomin3,ixomax3
4200 ixpmin2=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4203 jds=cache(iigrid)%source(ix^d)*ds_cm
4204 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4205 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4210 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4211 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4212 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4218 ix2=ixomin2+rell/rft
4219 do ix3=ixomin3,ixomax3
4220 ixpmin1=(
node(pig3_,igrid)-1)*rft*block_nx3+(ix3-ixomin3)*rft+1
4221 do ix1=ixomin1,ixomax1
4222 ixpmin2=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4225 jds=cache(iigrid)%source(ix^d)*ds_cm
4226 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4227 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4232 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4233 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4234 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4240 ix3=ixomin3+rell/rft
4241 do ix1=ixomin1,ixomax1
4242 ixpmin1=(
node(pig1_,igrid)-1)*rft*block_nx1+(ix1-ixomin1)*rft+1
4243 do ix2=ixomin2,ixomax2
4244 ixpmin2=(
node(pig2_,igrid)-1)*rft*block_nx2+(ix2-ixomin2)*rft+1
4247 jds=cache(iigrid)%source(ix^d)*ds_cm
4248 kds=cache(iigrid)%opacity(ix^d)*ds_cm
4249 jvds=cache(iigrid)%sourcev(ix^d)*ds_cm
4254 layer_ds(ixp1,ixp2,1,ibatch)=layer_ds(ixp1,ixp2,1,ibatch)+jds
4255 layer_ds(ixp1,ixp2,2,ibatch)=layer_ds(ixp1,ixp2,2,ibatch)+kds
4256 layer_ds(ixp1,ixp2,3,ibatch)=layer_ds(ixp1,ixp2,3,ibatch)+jvds
4265 numseg=nlayerseg*nbatch
4266 call mpi_allreduce(layer_ds,layer_all,numseg,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
4272 euvthin(ixp1,ixp2)=euvthin(ixp1,ixp2)+layer_all(ixp1,ixp2,1,ibatch)
4274 euv(ixp1,ixp2)=euv(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,1,ibatch)
4275 dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+atten*layer_all(ixp1,ixp2,3,ibatch)
4276 tau(ixp1,ixp2)=tau(ixp1,ixp2)+layer_all(ixp1,ixp2,2,ibatch)
4283 do iigrid=1,igridstail
4284 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
4285 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
4286 if (
allocated(cache(iigrid)%sourcev))
deallocate(cache(iigrid)%sourcev)
4288 deallocate(cache,layer_ds,layer_all)
4297 integer,
intent(in) :: ixI^L, ixO^L
4298 double precision,
intent(in) :: x(ixI^S,1:ndim),dx(ixI^S,1:ndim)
4299 double precision,
allocatable,
intent(out) :: rface(:),thetaface(:),phiface(:)
4300 integer :: ix1,ix2,ix3
4302 allocate(rface(ixomin1:ixomax1+1),thetaface(ixomin2:ixomax2+1),phiface(ixomin3:ixomax3+1))
4303 do ix1=ixomin1,ixomax1
4304 rface(ix1)=x(ix1,ixomin2,ixomin3,1)-half*dx(ix1,ixomin2,ixomin3,1)
4306 rface(ixomax1+1)=x(ixomax1,ixomin2,ixomin3,1)+half*dx(ixomax1,ixomin2,ixomin3,1)
4307 do ix2=ixomin2,ixomax2
4308 thetaface(ix2)=x(ixomin1,ix2,ixomin3,2)-half*dx(ixomin1,ix2,ixomin3,2)
4310 thetaface(ixomax2+1)=x(ixomin1,ixomax2,ixomin3,2)+half*dx(ixomin1,ixomax2,ixomin3,2)
4311 do ix3=ixomin3,ixomax3
4312 phiface(ix3)=x(ixomin1,ixomin2,ix3,3)-half*dx(ixomin1,ixomin2,ix3,3)
4314 phiface(ixomax3+1)=x(ixomin1,ixomin2,ixomax3,3)+half*dx(ixomin1,ixomin2,ixomax3,3)
4318 double precision,
allocatable,
intent(inout) :: tvals(:)
4319 integer,
intent(inout) :: nt,capacity
4320 double precision,
intent(in) :: t
4322 double precision,
allocatable :: tmp(:)
4324 if (t /= t .or. abs(t)>1.d90)
return
4325 if (.not.
allocated(tvals))
then
4327 allocate(tvals(capacity))
4328 else if (nt>=capacity)
then
4329 allocate(tmp(capacity))
4332 allocate(tvals(2*capacity))
4333 tvals(1:capacity)=tmp
4342 double precision,
intent(inout) :: tvals(:)
4343 integer,
intent(inout) :: nt
4344 double precision,
intent(in) :: t
4346 if (t /= t .or. abs(t)>1.d90)
return
4347 if (nt>=
size(tvals))
return
4353 double precision,
intent(inout) :: tvals(:)
4354 integer,
intent(inout) :: nt
4357 double precision :: key,epsT
4363 do while (j>=1 .and. tvals(j)>key)
4369 epst=max(1.d-12,1.d-10*max(one,abs(tvals(nt)-tvals(1))))
4372 if (abs(tvals(i)-tvals(nout))>epst)
then
4374 tvals(nout)=tvals(i)
4381 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),rface
4382 double precision,
allocatable,
intent(inout) :: tvals(:)
4383 integer,
intent(inout) :: nt,capacity
4385 double precision :: aa,bb,cc,disc,root
4388 bb=2.d0*sum(ray_origin*ray_dir)
4389 cc=sum(ray_origin**2)-rface**2
4390 disc=bb**2-4.d0*aa*cc
4391 if (disc<zero)
return
4392 root=sqrt(max(zero,disc))
4393 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4394 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4398 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),thetaface
4399 double precision,
allocatable,
intent(inout) :: tvals(:)
4400 integer,
intent(inout) :: nt,capacity
4402 double precision :: cth,aa,bb,cc,disc,root
4405 aa=ray_dir(3)**2-cth**2*sum(ray_dir**2)
4406 bb=2.d0*(ray_origin(3)*ray_dir(3)-cth**2*sum(ray_origin*ray_dir))
4407 cc=ray_origin(3)**2-cth**2*sum(ray_origin**2)
4408 if (abs(aa)<1.d-14)
then
4409 if (abs(bb)>1.d-14)
call sph_add_t(tvals,nt,capacity,-cc/bb)
4412 disc=bb**2-4.d0*aa*cc
4413 if (disc<zero)
return
4414 root=sqrt(max(zero,disc))
4415 call sph_add_t(tvals,nt,capacity,(-bb-root)/(2.d0*aa))
4416 call sph_add_t(tvals,nt,capacity,(-bb+root)/(2.d0*aa))
4420 double precision,
intent(in) :: ray_origin(1:3),ray_dir(1:3),phiface
4421 double precision,
allocatable,
intent(inout) :: tvals(:)
4422 integer,
intent(inout) :: nt,capacity
4424 double precision :: normal(1:3),denom,numer
4426 normal(1)=-sin(phiface)
4427 normal(2)=cos(phiface)
4429 denom=sum(normal*ray_dir)
4430 if (abs(denom)<1.d-14)
return
4431 numer=sum(normal*ray_origin)
4432 call sph_add_t(tvals,nt,capacity,-numer/denom)
4436 double precision,
intent(in) :: pos(1:3)
4437 double precision,
intent(out) :: sph(1:3)
4439 sph(1)=sqrt(sum(pos**2))
4440 if (sph(1)>zero)
then
4441 sph(2)=acos(max(-one,min(one,pos(3)/sph(1))))
4445 sph(3)=atan2(pos(2),pos(1))
4449 integer,
intent(in) :: imin,imax
4450 double precision,
intent(in) ::
value,faces(imin:imax+1)
4452 integer :: ilo,ihi,imid
4455 if (
value<faces(imin)-1.d-12 .or.
value>faces(imax+1)+1.d-12)
return
4456 if (
value<=faces(imin))
then
4460 if (
value>=faces(imax+1))
then
4467 do while (ihi-ilo>1)
4469 if (
value>=faces(imid))
then
4475 idx=min(imax,max(imin,ilo))
4479 integer,
intent(in) :: imin,imax
4480 double precision,
intent(in) ::
value,faces(imin:imax+1)
4482 integer :: ilo,ihi,imid
4485 if (
value>faces(imin)+1.d-12 .or.
value<faces(imax+1)-1.d-12)
return
4486 if (
value>=faces(imin))
then
4490 if (
value<=faces(imax+1))
then
4497 do while (ihi-ilo>1)
4499 if (
value<=faces(imid))
then
4505 idx=min(imax,max(imin,ilo))
4509 double precision,
intent(in) :: pos(1:3)
4510 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4511 phiface(ixOmin3:ixOmax3+1)
4512 integer,
intent(in) :: ixO^L
4513 integer,
intent(out) :: ix1,ix2,ix3
4514 logical,
intent(out) :: inside
4516 double precision :: sph(1:3),phi
4520 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4521 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4525 inside=ix1>0 .and. ix2>0 .and. ix3>0
4530 double precision,
intent(in) :: pos(1:3),ximg1,ximg2
4532 double precision :: dotp,rc,rthick,rloc
4535 rc=sqrt(ximg1**2+ximg2**2)
4536 rloc=sqrt(sum(pos**2))
4539 if (dotp>=
zero)
then
4540 if (rc<=rthick) visible=.false.
4542 if (rloc<=rthick) visible=.false.
4547 ixPmin1,ixPmax1,ixPmin2,ixPmax2,has_pixels)
4548 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4549 phiface(ixOmin3:ixOmax3+1)
4550 integer,
intent(in) :: ixO^L,nXI1,nXI2
4551 double precision,
intent(in) :: xI1(nXI1),xI2(nXI2),dxI
4552 integer,
intent(out) :: ixPmin1,ixPmax1,ixPmin2,ixPmax2
4553 logical,
intent(out) :: has_pixels
4555 integer,
parameter :: nsample=5
4557 double precision :: sph(1:3),xcent(1:2)
4558 double precision :: xmin1,xmax1,xmin2,xmax2
4559 double precision :: wr,wt,wp,pad
4567 wr=dble(ir)/dble(nsample-1)
4568 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
4570 wt=dble(it)/dble(nsample-1)
4571 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
4573 wp=dble(ip)/dble(nsample-1)
4574 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
4575 ip/=0 .and. ip/=nsample-1) cycle
4576 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
4578 xmin1=min(xmin1,xcent(1))
4579 xmax1=max(xmax1,xcent(1))
4580 xmin2=min(xmin2,xcent(2))
4581 xmax2=max(xmax2,xcent(2))
4590 ixpmin1=max(1,floor((xmin1-(xi1(1)-half*dxi))/dxi)+1)
4591 ixpmax1=min(nxi1,ceiling((xmax1-(xi1(1)-half*dxi))/dxi))
4592 ixpmin2=max(1,floor((xmin2-(xi2(1)-half*dxi))/dxi)+1)
4593 ixpmax2=min(nxi2,ceiling((xmax2-(xi2(1)-half*dxi))/dxi))
4594 has_pixels=ixpmin1<=ixpmax1 .and. ixpmin2<=ixpmax2
4598 rface,thetaface,phiface,EUVp)
4599 integer,
intent(in) :: ixI^L,ixO^L
4600 double precision,
intent(in) :: source(ixI^S),ray_origin(1:3),ximg1,ximg2
4601 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4602 phiface(ixOmin3:ixOmax3+1)
4603 double precision,
intent(inout) :: EUVp
4605 integer :: nt,capacity,i,ix^D
4606 double precision,
allocatable :: tvals(:)
4607 double precision :: posMid(1:3),ds_cm,tMid,t0,t1
4612 do ix1=ixomin1,ixomax1+1
4615 do ix2=ixomin2,ixomax2+1
4618 do ix3=ixomin3,ixomax3+1
4622 if (
allocated(tvals))
deallocate(tvals)
4632 posmid=ray_origin+tmid*
vec_los
4634 call sph_locate_cell(posmid,rface,thetaface,phiface,ixo^l,ix1,ix2,ix3,inside)
4635 if (.not. inside) cycle
4636 ds_cm=(t1-t0)*unit_length
4637 if (si_unit) ds_cm=ds_cm*1.d2
4638 euvp=euvp+source(ix^d)*ds_cm
4644 pixel_id,ray_origin,ximg1,ximg2,&
4645 rface,thetaface,phiface,rface2,&
4646 theta_cos,phi_sin,phi_cos,&
4647 segments,nseg,capacity)
4650 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4651 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4652 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4653 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4654 phiface(ixOmin3:ixOmax3+1)
4655 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4656 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4657 double precision,
allocatable,
intent(inout) :: segments(:,:)
4658 integer,
intent(inout) :: nseg,capacity
4660 integer :: nt,i,ix^D
4661 double precision :: tvals(2*(ixOmax1-ixOmin1+2)+2*(ixOmax2-ixOmin2+2)+&
4662 (ixOmax3-ixOmin3+2))
4663 double precision :: posMid(1:3),ds_cm,tMid,t0,t1,jds,kds
4664 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2
4665 double precision :: denom,numer,r2,mu,phi,dotp,rthick2,rc2
4670 origin2=sum(ray_origin**2)
4672 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4673 rc2=ximg1**2+ximg2**2
4675 do ix1=ixomin1,ixomax1+1
4678 cc=origin2-rface2(ix1)
4679 disc=bb**2-4.d0*aa*cc
4680 if (disc>=
zero)
then
4681 root=sqrt(max(
zero,disc))
4686 do ix2=ixomin2,ixomax2+1
4687 cth2=theta_cos(ix2)**2
4689 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4690 cc=ray_origin(3)**2-cth2*origin2
4691 if (abs(aa)<1.d-14)
then
4694 disc=bb**2-4.d0*aa*cc
4695 if (disc>=
zero)
then
4696 root=sqrt(max(
zero,disc))
4702 do ix3=ixomin3,ixomax3+1
4704 if (abs(denom)>=1.d-14)
then
4705 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
4717 posmid=ray_origin+tmid*
vec_los
4720 if (dotp>=
zero)
then
4721 if (rc2<=rthick2) cycle
4723 if (r2<=rthick2) cycle
4728 mu=posmid(3)/sqrt(r2)
4733 phi=atan2(posmid(2),posmid(1))
4734 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*
dpi
4735 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*
dpi
4737 inside=ix1>0 .and. ix2>0 .and. ix3>0
4738 if (.not. inside) cycle
4739 ds_cm=(t1-t0)*unit_length
4740 if (si_unit) ds_cm=ds_cm*1.d2
4741 jds=max(
zero,source(ix^d))*ds_cm
4742 kds=max(
zero,opacity(ix^d))*ds_cm
4748 double precision,
intent(in) :: pos(1:3)
4749 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4750 phiface(ixOmin3:ixOmax3+1)
4751 integer,
intent(in) :: ixO^L
4752 integer,
intent(out) :: ix1,ix2,ix3
4753 logical,
intent(out) :: inside
4755 double precision :: r2,mu,phi
4765 phi=atan2(pos(2),pos(1))
4766 if (phi<phiface(ixomin3)-1.d-12) phi=phi+2.d0*dpi
4767 if (phi>phiface(ixomax3+1)+1.d-12) phi=phi-2.d0*dpi
4769 inside=ix1>0 .and. ix2>0 .and. ix3>0
4773 double precision,
intent(in) :: t,tNow,tExit,epsRay
4774 double precision,
intent(inout) :: tNext
4775 logical,
intent(inout) :: found
4777 if (t>tnow+epsray .and. t<=texit+epsray .and. t<tnext)
then
4784 double precision,
intent(in) :: t,theta_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4785 double precision,
intent(inout) :: tNext
4786 logical,
intent(inout) :: found
4788 double precision :: pos(1:3),r2
4790 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4793 if (r2<=zero)
return
4794 if (theta_face_cos>1.d-12 .and. pos(3)<-1.d-10)
return
4795 if (theta_face_cos<-1.d-12 .and. pos(3)>1.d-10)
return
4796 if (abs(pos(3)**2-theta_face_cos**2*r2)>1.d-6*max(one,r2))
return
4802 double precision,
intent(in) :: t,phi_face_sin,phi_face_cos,ray_origin(1:3),tNow,tExit,epsRay
4803 double precision,
intent(inout) :: tNext
4804 logical,
intent(inout) :: found
4806 double precision :: pos(1:3),radialDot
4808 if (t<=tnow+epsray .or. t>texit+epsray .or. t>=tnext)
return
4810 radialdot=phi_face_cos*pos(1)+phi_face_sin*pos(2)
4811 if (radialdot<-1.d-10)
return
4817 ix1,ix2,ix3,tNow,tExit,epsRay,tNext,found)
4818 double precision,
intent(in) :: ray_origin(1:3)
4819 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4820 phiface(ixOmin3:ixOmax3+1)
4821 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4822 integer,
intent(in) :: ixO^L,ix1,ix2,ix3
4823 double precision,
intent(in) :: tNow,tExit,epsRay
4824 double precision,
intent(out) :: tNext
4825 logical,
intent(out) :: found
4828 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4833 origin2=sum(ray_origin**2)
4839 cc=origin2-rface2(iface)
4840 disc=bb**2-4.d0*aa*cc
4841 if (disc>=zero)
then
4842 root=sqrt(max(zero,disc))
4849 cth2=theta_cos(iface)**2
4851 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4852 cc=ray_origin(3)**2-cth2*origin2
4853 if (abs(aa)<1.d-14)
then
4855 ray_origin,tnow,texit,epsray,tnext,found)
4857 disc=bb**2-4.d0*aa*cc
4858 if (disc>=zero)
then
4859 root=sqrt(max(zero,disc))
4861 ray_origin,tnow,texit,epsray,tnext,found)
4863 ray_origin,tnow,texit,epsray,tnext,found)
4870 if (abs(denom)>=1.d-14)
then
4871 numer=-phi_sin(iface)*ray_origin(1)+phi_cos(iface)*ray_origin(2)
4873 tnow,texit,epsray,tnext,found)
4879 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,&
4881 t_enter,t_exit,segments,nseg,capacity,ok)
4884 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4885 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4886 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4887 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4888 phiface(ixOmin3:ixOmax3+1)
4889 double precision,
intent(in) :: phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4890 double precision,
intent(in) :: t_enter,t_exit
4891 double precision,
allocatable,
intent(inout) :: segments(:,:)
4892 integer,
intent(inout) :: nseg,capacity
4893 logical,
intent(out) :: ok
4895 integer :: ix^D,nstep,maxSteps
4896 double precision :: tNow,tNext,tEnd,tMid,epsRay,ds_cm,jds,kds
4897 double precision :: pos(1:3),r2,dotp,rthick2,rc2
4898 logical :: inside,found
4901 if (t_exit<=t_enter)
return
4902 epsray=max(1.d-12,1.d-10*max(
one,abs(t_exit-t_enter)))
4903 pos=ray_origin+(t_enter+epsray)*
vec_los
4905 if (.not. inside)
then
4910 rthick2=(r_opt_thick*
const_rsun/unit_length)**2
4911 rc2=ximg1**2+ximg2**2
4914 maxsteps=8*((ixomax1-ixomin1+1)+(ixomax2-ixomin2+1)+(ixomax3-ixomin3+1)+3)
4916 do while (tnow<t_exit-epsray)
4918 ix1,ix2,ix3,tnow,t_exit,epsray,tnext,found)
4919 if (.not. found)
then
4923 tend=min(tnext,t_exit)
4925 tmid=
half*(tnow+tend)
4929 if (.not. ((dotp>=
zero .and. rc2<=rthick2) .or. &
4930 (dotp<
zero .and. r2<=rthick2)))
then
4931 ds_cm=(tend-tnow)*unit_length
4932 if (si_unit) ds_cm=ds_cm*1.d2
4933 jds=max(
zero,source(ix^d))*ds_cm
4934 kds=max(
zero,opacity(ix^d))*ds_cm
4939 if (tnow>=t_exit-epsray)
exit
4940 pos=ray_origin+(tnow+epsray)*
vec_los
4942 if (.not. inside)
then
4947 if (nstep>maxsteps)
then
4955 pixel_id,ray_origin,ximg1,ximg2,&
4956 rface,thetaface,phiface,rface2,&
4957 theta_cos,phi_sin,phi_cos,&
4958 segments,nseg,capacity,fallback)
4959 integer,
intent(in) :: ixI^L,ixO^L,pixel_id
4960 double precision,
intent(in) :: source(ixI^S),opacity(ixI^S)
4961 double precision,
intent(in) :: ray_origin(1:3),ximg1,ximg2
4962 double precision,
intent(in) :: rface(ixOmin1:ixOmax1+1),thetaface(ixOmin2:ixOmax2+1),&
4963 phiface(ixOmin3:ixOmax3+1)
4964 double precision,
intent(in) :: rface2(ixOmin1:ixOmax1+1),theta_cos(ixOmin2:ixOmax2+1),&
4965 phi_sin(ixOmin3:ixOmax3+1),phi_cos(ixOmin3:ixOmax3+1)
4966 double precision,
allocatable,
intent(inout) :: segments(:,:)
4967 integer,
intent(inout) :: nseg,capacity
4968 logical,
intent(out) :: fallback
4970 integer :: nt,i,nsegStart,ix^D
4971 double precision :: tvals(12),posMid(1:3),t0,t1,tMid
4972 double precision :: dir2,origin2,odotd,aa,bb,cc,disc,root,cth2,denom,numer
4973 logical :: inside,ok
4979 origin2=sum(ray_origin**2)
4982 do ix1=ixomin1,ixomax1+1,ixomax1-ixomin1+1
4985 cc=origin2-rface2(ix1)
4986 disc=bb**2-4.d0*aa*cc
4987 if (disc>=zero)
then
4988 root=sqrt(max(zero,disc))
4993 do ix2=ixomin2,ixomax2+1,ixomax2-ixomin2+1
4994 cth2=theta_cos(ix2)**2
4996 bb=2.d0*(ray_origin(3)*
vec_los(3)-cth2*odotd)
4997 cc=ray_origin(3)**2-cth2*origin2
4998 if (abs(aa)<1.d-14)
then
5001 disc=bb**2-4.d0*aa*cc
5002 if (disc>=zero)
then
5003 root=sqrt(max(zero,disc))
5009 do ix3=ixomin3,ixomax3+1,ixomax3-ixomin3+1
5011 if (abs(denom)>=1.d-14)
then
5012 numer=-phi_sin(ix3)*ray_origin(1)+phi_cos(ix3)*ray_origin(2)
5025 posmid=ray_origin+tmid*
vec_los
5027 if (.not. inside) cycle
5029 ray_origin,ximg1,ximg2,rface2,theta_cos,phiface,phi_sin,phi_cos,&
5030 t0,t1,segments,nseg,capacity,ok)
5035 pixel_id,ray_origin,ximg1,ximg2,rface,thetaface,phiface,rface2,&
5036 theta_cos,phi_sin,phi_cos,segments,nseg,capacity)
5046 integer,
intent(in) :: numXI1,numXI2
5047 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5049 double precision,
intent(inout) :: EM(numXI1,numXI2)
5051 integer :: ixO^L,ixI^L,ix^D
5052 integer :: iigrid,igrid,ixP1,ixP2,ixPmin1,ixPmax1,ixPmin2,ixPmax2
5053 double precision :: ray_origin(1:3),profile_local(3),profile_global(3)
5054 double precision,
allocatable :: source(:^D&)
5055 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5056 logical :: has_pixels
5059 do iigrid=1,igridstail; igrid=igrids(iigrid);
5060 ^d&ixomin^d=ixmlo^d\
5061 ^d&ixomax^d=ixmhi^d\
5062 ^d&iximin^d=
ixglo^d\
5063 ^d&iximax^d=
ixghi^d\
5065 allocate(source(ixi^s))
5070 ixpmin1,ixpmax1,ixpmin2,ixpmax2,has_pixels)
5071 if (has_pixels)
then
5072 do ixp1=ixpmin1,ixpmax1
5073 do ixp2=ixpmin2,ixpmax2
5075 profile_local(1)=profile_local(1)+one
5077 rface,thetaface,phiface,em(ixp1,ixp2))
5080 profile_local(2)=profile_local(2)+dble((ixpmax1-ixpmin1+1)*(ixpmax2-ixpmin2+1))
5082 profile_local(3)=profile_local(3)+one
5083 deallocate(source,rface,thetaface,phiface)
5085 call mpi_allreduce(profile_local,profile_global,3,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5087 write(*,
'(a,3(es12.5,1x))')
' sph_intersection thin profile rays pixels blocks: ',profile_global
5094 integer,
intent(in) :: numXI1,numXI2
5095 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2),dxI
5097 double precision,
intent(out) :: EUV(numXI1,numXI2),Tau(numXI1,numXI2),EUVthin(numXI1,numXI2)
5099 integer,
parameter :: nSegVars=5
5100 integer :: ixO^L,ixI^L,ix^D
5101 integer :: iigrid,igrid,ixP1,ixP2,ipix,ipixStart,ipixEnd,nPixBatch,pixel_id
5102 integer :: nseg,capacity,totalCount,totalSeg,ipe,is,iseg,nidx,owner,isegDest,nsegBefore
5103 integer :: ixGlobal,iyGlobal,ixPmin1,ixPmax1,ixPmin2,ixPmax2,iFirst,iLast,iLocal
5104 integer :: nPixBatchTarget
5105 integer :: maxSegBatchTarget,maxSegCommTarget,maxNsegBatch,nPixTotal
5106 integer :: maxOwnerSegCount,maxOwnerSegCountLocal,segOffset,recvFill,totalRoundCount,totalRoundSeg
5107 integer :: sphDdaFallbackLocal,sphDdaFallbackGlobal
5108 integer,
allocatable :: sendCounts(:),recvCounts(:),sendDispls(:),recvDispls(:)
5109 integer,
allocatable :: roundSendCounts(:),roundRecvCounts(:)
5110 integer,
allocatable :: roundSendDispls(:),roundRecvDispls(:)
5111 integer,
allocatable :: ownerSegCounts(:),ownerOffsets(:),idx(:)
5112 integer,
allocatable :: bucketCounts(:),bucketOffsets(:),bucketFill(:)
5113 double precision :: ray_origin(1:3),atten
5114 double precision :: profile_local(5),profile_global(5),profile_batch(5)
5115 double precision :: phys_max_local(2),phys_max_global(2)
5116 double precision :: phys_sum_local(2),phys_sum_global(2),phys_sum_batch(2)
5117 double precision,
allocatable :: segments(:,:),segments_send(:,:),segments_recv(:,:)
5118 double precision,
allocatable :: segments_recv_round(:,:)
5119 double precision,
allocatable :: image_reduce(:,:)
5120 logical :: has_pixels,batchAccepted,batchReduced,ddaFallback
5129 sphddafallbacklocal=0
5130 allocate(sendcounts(0:
npe-1),recvcounts(0:
npe-1),senddispls(0:
npe-1),recvdispls(0:
npe-1))
5131 allocate(roundsendcounts(0:
npe-1),roundrecvcounts(0:
npe-1))
5132 allocate(roundsenddispls(0:
npe-1),roundrecvdispls(0:
npe-1))
5133 allocate(ownersegcounts(0:
npe-1),owneroffsets(0:
npe-1))
5134 allocate(cache(igridstail))
5136 allocate(bucketcounts(npixbatchtarget),bucketoffsets(npixbatchtarget+1),&
5137 bucketfill(npixbatchtarget))
5139 do iigrid=1,igridstail; igrid=igrids(iigrid);
5140 ^d&ixomin^d=ixmlo^d\
5141 ^d&ixomax^d=ixmhi^d\
5142 ^d&iximin^d=
ixglo^d\
5143 ^d&iximax^d=
ixghi^d\
5145 cache(iigrid)%igrid=igrid
5146 allocate(cache(iigrid)%source(ixi^s),cache(iigrid)%opacity(ixi^s))
5147 cache(iigrid)%source=zero
5148 cache(iigrid)%opacity=zero
5149 call get_euv(
wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,cache(iigrid)%source)
5152 phys_max_local(1)=max(phys_max_local(1),maxval(cache(iigrid)%source(ixo^s)))
5153 phys_max_local(2)=max(phys_max_local(2),maxval(cache(iigrid)%opacity(ixo^s)))
5155 cache(iigrid)%rface,cache(iigrid)%thetaface,&
5156 cache(iigrid)%phiface)
5157 allocate(cache(iigrid)%rface2(ixomin1:ixomax1+1),&
5158 cache(iigrid)%theta_cos(ixomin2:ixomax2+1),&
5159 cache(iigrid)%phi_sin(ixomin3:ixomax3+1),&
5160 cache(iigrid)%phi_cos(ixomin3:ixomax3+1))
5161 cache(iigrid)%rface2=cache(iigrid)%rface**2
5162 cache(iigrid)%theta_cos=cos(cache(iigrid)%thetaface)
5163 cache(iigrid)%phi_sin=sin(cache(iigrid)%phiface)
5164 cache(iigrid)%phi_cos=cos(cache(iigrid)%phiface)
5166 cache(iigrid)%phiface,ixo^l,numxi1,numxi2,xi1,xi2,dxi,&
5167 cache(iigrid)%ixPmin1,cache(iigrid)%ixPmax1,&
5168 cache(iigrid)%ixPmin2,cache(iigrid)%ixPmax2,cache(iigrid)%has_pixels)
5171 npixtotal=numxi1*numxi2
5173 do while (ipixstart<=npixtotal)
5174 ipixend=min(numxi1*numxi2,ipixstart+npixbatchtarget-1)
5175 npixbatch=ipixend-ipixstart+1
5176 batchaccepted=.false.
5177 batchreduced=.false.
5179 do while (.not. batchaccepted)
5185 do iigrid=1,igridstail; igrid=igrids(iigrid);
5186 ^d&ixomin^d=ixmlo^d\
5187 ^d&ixomax^d=ixmhi^d\
5188 ^d&iximin^d=
ixglo^d\
5189 ^d&iximax^d=
ixghi^d\
5191 ixpmin1=cache(iigrid)%ixPmin1
5192 ixpmax1=cache(iigrid)%ixPmax1
5193 ixpmin2=cache(iigrid)%ixPmin2
5194 ixpmax2=cache(iigrid)%ixPmax2
5195 has_pixels=cache(iigrid)%has_pixels
5196 if (.not. has_pixels) cycle
5198 do ixp2=ixpmin2,ixpmax2
5199 ifirst=max(ipixstart,(ixp2-1)*numxi1+ixpmin1)
5200 ilast=min(ipixend,(ixp2-1)*numxi1+ixpmax1)
5201 if (ifirst>ilast) cycle
5202 do ipix=ifirst,ilast
5203 ixp1=1+mod(ipix-1,numxi1)
5206 profile_batch(1)=profile_batch(1)+one
5210 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5211 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5212 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5213 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5214 segments,nseg,capacity,ddafallback)
5215 if (ddafallback) sphddafallbacklocal=sphddafallbacklocal+1
5218 cache(iigrid)%opacity,pixel_id,ray_origin,xi1(ixp1),xi2(ixp2),&
5219 cache(iigrid)%rface,cache(iigrid)%thetaface,cache(iigrid)%phiface,&
5220 cache(iigrid)%rface2,cache(iigrid)%theta_cos,&
5221 cache(iigrid)%phi_sin,cache(iigrid)%phi_cos,&
5222 segments,nseg,capacity)
5224 if (nseg>nsegbefore) profile_batch(2)=profile_batch(2)+one
5225 profile_batch(3)=profile_batch(3)+dble(nseg-nsegbefore)
5226 do iseg=nsegbefore+1,nseg
5227 phys_sum_batch(1)=phys_sum_batch(1)+segments(3,iseg)
5228 phys_sum_batch(2)=phys_sum_batch(2)+segments(4,iseg)
5234 call mpi_allreduce(nseg,maxnsegbatch,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5235 if (maxnsegbatch>maxsegbatchtarget .and. npixbatch>1)
then
5236 npixbatch=max(1,npixbatch/2)
5237 ipixend=ipixstart+npixbatch-1
5238 if (
allocated(segments))
deallocate(segments)
5241 batchaccepted=.true.
5245 profile_local=profile_local+profile_batch
5246 phys_sum_local=phys_sum_local+phys_sum_batch
5249 write(*,
'(a,3(i0,1x))')
' sph_dda thick adaptive batch: ',&
5250 ipixstart,ipixend,maxnsegbatch
5252 write(*,
'(a,3(i0,1x))')
' sph_intersection thick adaptive batch: ',&
5253 ipixstart,ipixend,maxnsegbatch
5257 if (.not.
allocated(segments))
then
5259 allocate(segments(nsegvars,capacity))
5264 ownersegcounts(owner)=ownersegcounts(owner)+1
5266 sendcounts=nsegvars*ownersegcounts
5269 senddispls(ipe)=senddispls(ipe-1)+sendcounts(ipe-1)
5272 allocate(segments_send(nsegvars,max(1,nseg)))
5276 isegdest=senddispls(owner)/nsegvars+owneroffsets(owner)+1
5277 segments_send(:,isegdest)=segments(:,is)
5278 owneroffsets(owner)=owneroffsets(owner)+1
5281 call mpi_alltoall(sendcounts,1,mpi_integer,recvcounts,1,mpi_integer,
icomm,
ierrmpi)
5284 recvdispls(ipe)=recvdispls(ipe-1)+recvcounts(ipe-1)
5286 totalcount=sum(recvcounts)
5287 totalseg=totalcount/nsegvars
5288 profile_local(4)=profile_local(4)+dble(totalcount)
5289 allocate(segments_recv(nsegvars,max(1,totalseg)))
5292 maxownersegcountlocal=maxval(ownersegcounts)
5293 call mpi_allreduce(maxownersegcountlocal,maxownersegcount,1,mpi_integer,mpi_max,
icomm,
ierrmpi)
5294 do segoffset=0,maxownersegcount-1,maxsegcommtarget
5296 roundsenddispls=senddispls
5298 if (ownersegcounts(ipe)>segoffset)
then
5299 roundsendcounts(ipe)=nsegvars*min(maxsegcommtarget,ownersegcounts(ipe)-segoffset)
5300 roundsenddispls(ipe)=senddispls(ipe)+nsegvars*segoffset
5304 call mpi_alltoall(roundsendcounts,1,mpi_integer,roundrecvcounts,1,mpi_integer,
icomm,
ierrmpi)
5305 roundrecvdispls(0)=0
5307 roundrecvdispls(ipe)=roundrecvdispls(ipe-1)+roundrecvcounts(ipe-1)
5309 totalroundcount=sum(roundrecvcounts)
5310 totalroundseg=totalroundcount/nsegvars
5311 allocate(segments_recv_round(nsegvars,max(1,totalroundseg)))
5313 call mpi_alltoallv(segments_send,roundsendcounts,roundsenddispls,mpi_double_precision,&
5314 segments_recv_round,roundrecvcounts,roundrecvdispls,&
5317 if (totalroundseg>0)
then
5318 segments_recv(:,recvfill+1:recvfill+totalroundseg)=segments_recv_round(:,1:totalroundseg)
5319 recvfill=recvfill+totalroundseg
5321 deallocate(segments_recv_round)
5324 if (recvfill/=totalseg)
call mpistop(
"ray-segment receive mismatch")
5326 if (totalseg>0)
then
5327 allocate(idx(totalseg))
5328 bucketcounts(1:npixbatch)=0
5331 ipix=nint(segments_recv(1,is))
5333 ilocal=ipix-ipixstart+1
5334 bucketcounts(ilocal)=bucketcounts(ilocal)+1
5340 do ilocal=1,npixbatch
5341 bucketoffsets(ilocal+1)=bucketoffsets(ilocal)+bucketcounts(ilocal)
5343 bucketfill(1:npixbatch)=bucketoffsets(1:npixbatch)
5346 ipix=nint(segments_recv(1,is))
5348 ilocal=ipix-ipixstart+1
5349 idx(bucketfill(ilocal))=is
5350 bucketfill(ilocal)=bucketfill(ilocal)+1
5355 do ipix=ipixstart,ipixend
5357 ilocal=ipix-ipixstart+1
5358 nidx=bucketcounts(ilocal)
5360 profile_local(5)=profile_local(5)+dble(nidx)*dble(nidx)
5362 idx(bucketoffsets(ilocal):bucketoffsets(ilocal+1)-1),nidx)
5363 ixglobal=1+mod(ipix-1,numxi1)
5364 iyglobal=1+(ipix-1)/numxi1
5365 do iseg=bucketoffsets(ilocal),bucketoffsets(ilocal+1)-1
5367 euvthin(ixglobal,iyglobal)=euvthin(ixglobal,iyglobal)+segments_recv(3,is)
5369 euv(ixglobal,iyglobal)=euv(ixglobal,iyglobal)+atten*segments_recv(3,is)
5370 tau(ixglobal,iyglobal)=tau(ixglobal,iyglobal)+max(zero,segments_recv(4,is))
5377 deallocate(segments_send,segments_recv)
5378 if (
allocated(segments))
deallocate(segments)
5382 do iigrid=1,igridstail
5383 if (
allocated(cache(iigrid)%source))
deallocate(cache(iigrid)%source)
5384 if (
allocated(cache(iigrid)%opacity))
deallocate(cache(iigrid)%opacity)
5385 if (
allocated(cache(iigrid)%rface))
deallocate(cache(iigrid)%rface)
5386 if (
allocated(cache(iigrid)%thetaface))
deallocate(cache(iigrid)%thetaface)
5387 if (
allocated(cache(iigrid)%phiface))
deallocate(cache(iigrid)%phiface)
5388 if (
allocated(cache(iigrid)%rface2))
deallocate(cache(iigrid)%rface2)
5389 if (
allocated(cache(iigrid)%theta_cos))
deallocate(cache(iigrid)%theta_cos)
5390 if (
allocated(cache(iigrid)%phi_sin))
deallocate(cache(iigrid)%phi_sin)
5391 if (
allocated(cache(iigrid)%phi_cos))
deallocate(cache(iigrid)%phi_cos)
5394 deallocate(sendcounts,recvcounts,senddispls,recvdispls,roundsendcounts,roundrecvcounts,&
5395 roundsenddispls,roundrecvdispls,ownersegcounts,owneroffsets,bucketcounts,&
5396 bucketoffsets,bucketfill)
5397 allocate(image_reduce(numxi1,numxi2))
5398 call mpi_allreduce(euv,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5400 call mpi_allreduce(tau,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5402 call mpi_allreduce(euvthin,image_reduce,numxi1*numxi2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5403 euvthin=image_reduce
5404 deallocate(image_reduce)
5405 call mpi_allreduce(profile_local,profile_global,5,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5406 call mpi_allreduce(phys_max_local,phys_max_global,2,mpi_double_precision,mpi_max,
icomm,
ierrmpi)
5407 call mpi_allreduce(phys_sum_local,phys_sum_global,2,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5408 call mpi_allreduce(sphddafallbacklocal,sphddafallbackglobal,1,mpi_integer,mpi_sum,
icomm,
ierrmpi)
5411 write(*,
'(a,5(es12.5,1x))')
' sph_dda thick profile: ',profile_global
5412 write(*,
'(a,i0)')
' sph_dda thick fallback rays: ',sphddafallbackglobal
5413 write(*,
'(a,4(es12.5,1x))')
' sph_dda thick physics maxj maxk sumjds sumkds: ',&
5414 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5416 write(*,
'(a,5(es12.5,1x))')
' sph_intersection thick profile: ',profile_global
5417 write(*,
'(a,4(es12.5,1x))')
' sph_intersection thick physics maxj maxk sumjds sumkds: ',&
5418 phys_max_global(1),phys_max_global(2),phys_sum_global(1),phys_sum_global(2)
5428 double precision,
intent(out) :: xImin1,xImax1,xImin2,xImax2
5430 integer,
parameter :: nsample=5
5431 integer :: iigrid,igrid,ir,it,ip
5432 integer :: ixI^L,ixO^L
5433 double precision,
allocatable :: rface(:),thetaface(:),phiface(:)
5434 double precision :: local_min1,local_max1,local_min2,local_max2
5435 double precision :: sph(1:3),xcent(1:2),wr,wt,wp
5437 local_min1=huge(one)
5438 local_max1=-huge(one)
5439 local_min2=huge(one)
5440 local_max2=-huge(one)
5442 do iigrid=1,igridstail
5443 igrid=igrids(iigrid)
5444 ^d&ixomin^d=ixmlo^d\
5445 ^d&ixomax^d=ixmhi^d\
5446 ^d&iximin^d=ixglo^d\
5447 ^d&iximax^d=ixghi^d\
5450 rface,thetaface,phiface)
5452 wr=dble(ir)/dble(nsample-1)
5453 sph(1)=(one-wr)*rface(ixomin1)+wr*rface(ixomax1+1)
5455 wt=dble(it)/dble(nsample-1)
5456 sph(2)=(one-wt)*thetaface(ixomin2)+wt*thetaface(ixomax2+1)
5458 wp=dble(ip)/dble(nsample-1)
5459 if (ir/=0 .and. ir/=nsample-1 .and. it/=0 .and. it/=nsample-1 .and. &
5460 ip/=0 .and. ip/=nsample-1) cycle
5461 sph(3)=(one-wp)*phiface(ixomin3)+wp*phiface(ixomax3+1)
5463 local_min1=min(local_min1,xcent(1))
5464 local_max1=max(local_max1,xcent(1))
5465 local_min2=min(local_min2,xcent(2))
5466 local_max2=max(local_max2,xcent(2))
5470 deallocate(rface,thetaface,phiface)
5473 call mpi_allreduce(local_min1,ximin1,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5474 call mpi_allreduce(local_max1,ximax1,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5475 call mpi_allreduce(local_min2,ximin2,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5476 call mpi_allreduce(local_max2,ximax2,1,mpi_double_precision,mpi_max,icomm,ierrmpi)
5477 if (ximin1>0.5d0*huge(one) .or. ximax1<-0.5d0*huge(one) .or. &
5478 ximin2>0.5d0*huge(one) .or. ximax2<-0.5d0*huge(one))
then
5479 call mpistop(
"sph_intersection could not determine image bounds")
5484 double precision,
intent(out) :: dxI
5486 integer :: iigrid,igrid,ixI^L,ixO^L,ix^D
5487 double precision :: local_min,global_min,dr,ds_theta,ds_phi,rval,theta
5490 do iigrid=1,igridstail
5491 igrid=igrids(iigrid)
5492 ^d&ixomin^d=ixmlo^d\
5493 ^d&ixomax^d=ixmhi^d\
5494 ^d&iximin^d=ixglo^d\
5495 ^d&iximax^d=ixghi^d\
5497 do ix1=ixomin1,ixomax1
5498 do ix2=ixomin2,ixomax2
5499 do ix3=ixomin3,ixomax3
5500 rval=max(smalldouble,ps(igrid)%x(ix^d,1))
5501 theta=ps(igrid)%x(ix^d,2)
5502 dr=ps(igrid)%dx(ix^d,1)
5503 ds_theta=rval*ps(igrid)%dx(ix^d,2)
5504 ds_phi=rval*max(smalldouble,sin(theta))*ps(igrid)%dx(ix^d,3)
5505 local_min=min(local_min,dr,ds_theta,ds_phi)
5511 call mpi_allreduce(local_min,global_min,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
5512 if (global_min<=zero .or. global_min>half*huge(one))
then
5513 call mpistop(
"sph_intersection could not determine dat-resolution image spacing")
5524 integer,
intent(in) :: qunit
5526 character(20),
intent(in) :: datatype
5528 integer :: ix^D,numXI1,numXI2,numWI
5529 double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
5530 double precision,
allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
5531 double precision,
allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),Dpl(:,:),Tau(:,:),EMthin(:,:),WLB(:,:,:)
5532 double precision :: vec_temp1(1:3),vec_temp2(1:3)
5533 double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
5534 double precision :: res,LOS_psi,r_max,r_loc
5537 character (30) :: ion
5538 double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
5539 double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
5540 integer :: iigrid,igrid,i,j,numSI,iw
5541 logical :: emit,ray_image_global,has_thick_output
5543 if (coordinate==spherical)
then
5551 if (coordinate==spherical)
then
5555 ximin1=-abs(xprobmax1)
5556 ximin2=-abs(xprobmax1)
5557 ximax1=abs(xprobmax1)
5558 ximax2=abs(xprobmax1)
5563 if (ix1==1) vec_cor(1)=xprobmin1
5564 if (ix1==2) vec_cor(1)=xprobmax1
5566 if (ix2==1) vec_cor(2)=xprobmin2
5567 if (ix2==2) vec_cor(2)=xprobmax2
5569 if (ix3==1) vec_cor(3)=xprobmin3
5570 if (ix3==2) vec_cor(3)=xprobmax3
5573 r_loc=r_loc+(vec_cor(2)-
x_origin(2))**2
5574 r_loc=r_loc+(vec_cor(3)-
x_origin(3))**2
5576 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5579 r_max=max(r_max,r_loc)
5583 if (ix1==1 .and. ix2==1 .and. ix3==1)
then
5589 ximin1=min(ximin1,xi_cor(1))
5590 ximax1=max(ximax1,xi_cor(1))
5591 ximin2=min(ximin2,xi_cor(2))
5592 ximax2=max(ximax2,xi_cor(2))
5605 xicent1=(ximin1+ximax1)/2.d0
5606 xicent2=(ximin2+ximax2)/2.d0
5614 if (datatype==
'image_euv')
then
5618 else if (datatype==
'image_sxr')
then
5620 dxi=rhessi_rsl*arcsec
5622 else if (datatype==
'image_whitelight')
then
5633 call mpistop(
'Whitelight synthesis: instrument is not supported!')
5635 dxi=lasco_rsl*arcsec
5640 numxi1=8*ceiling((ximax1-xicent1)/dxi/8.d0)
5641 ximin1=xicent1-numxi1*dxi
5642 ximax1=xicent1+numxi1*dxi
5644 numxi2=8*ceiling((ximax2-xicent2)/dxi/8.d0)
5645 ximin2=xicent2-numxi2*dxi
5646 ximax2=xicent2+numxi2*dxi
5648 allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
5650 xi1(ix1)=ximin1+dxi*(ix1-
half)
5654 xi2(ix2)=ximin2+dxi*(ix2-
half)
5659 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
5660 has_thick_output=datatype==
'image_euv' .and. trim(
radiation_transfer)==
'thick' .and. &
5663 if (datatype==
'image_euv')
then
5668 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),em(numxi1,numxi2))
5672 ray_image_global=.false.
5673 if (has_thick_output)
then
5674 allocate(tau(numxi1,numxi2),emthin(numxi1,numxi2))
5678 if (coordinate==cartesian .and. datatype==
'image_euv' .and. &
5680 ray_image_global=.true.
5681 allocate(dpl(numxi1,numxi2))
5689 else if (coordinate==cartesian)
then
5690 do iigrid=1,igridstail; igrid=igrids(iigrid);
5693 else if (trim(
ray_method_active) ==
'spherical' .and. datatype ==
'image_euv')
then
5695 ray_image_global=.true.
5701 do iigrid=1,igridstail; igrid=igrids(iigrid);
5705 if (ray_image_global)
then
5706 if (has_thick_output)
then
5708 has_thick_output,tau=tau,euvthin=emthin,&
5709 cap_absorption=.true.)
5716 if (em(ix1,ix2)>smallflux) wis(ix1,ix2,1)=em(ix1,ix2)
5720 if (.not. ray_image_global)
then
5721 numsi=numxi1*numxi2*numwi
5722 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5730 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5731 if (
allocated(tau))
deallocate(tau)
5732 if (
allocated(emthin))
deallocate(emthin)
5733 deallocate(wi,wis,em)
5734 else if (datatype==
'image_whitelight')
then
5736 allocate(wi(numxi1,numxi2,numwi),wis(numxi1,numxi2,numwi),wlb(numxi1,numxi2,numwi))
5740 if (coordinate==spherical)
then
5741 do iigrid=1,igridstail; igrid=igrids(iigrid);
5747 if (wlb(ix1,ix2,1)>smallflux)
then
5748 wis(ix1,ix2,1)=wlb(ix1,ix2,1)
5749 wis(ix1,ix2,2)=wlb(ix1,ix2,2)
5753 numsi=numxi1*numxi2*numwi
5754 call mpi_allreduce(wis,wi,numsi,mpi_double_precision,mpi_sum,
icomm,
ierrmpi)
5761 call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
5762 deallocate(wi,wis,wlb)
5765 deallocate(xi1,xi2,dxi1,dxi2)
5770 integer,
intent(in) :: igrid,numXI1,numXI2
5771 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
5772 double precision,
intent(in) :: dxI
5774 character(20),
intent(in) :: datatype
5775 double precision,
intent(inout) :: EM(numXI1,numXI2)
5777 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
5778 double precision :: xb^L,xd^D
5779 double precision,
allocatable :: flux(:^D&),opacity(:^D&)
5780 double precision :: res
5781 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
5782 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
5783 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
5786 double precision :: logTe
5787 character (30) :: ion
5788 double precision :: lineCent
5789 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
5790 double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
5791 double precision :: aa,bb
5793 ^d&ixomin^d=ixmlo^d\
5794 ^d&ixomax^d=ixmhi^d\
5795 ^d&iximin^d=ixglo^d\
5796 ^d&iximax^d=ixghi^d\
5797 ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
5798 ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
5801 arcsec=7.25d5/unit_length
5803 arcsec=7.25d7/unit_length
5806 allocate(flux(ixi^s),opacity(ixi^s))
5807 if (datatype==
'image_euv')
then
5808 if (trim(emission_model)==
'pseudo_current')
then
5810 else if (trim(emission_model)==
'radio_ff')
then
5814 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
5815 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
5817 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
5818 pixel=spacersl*arcsec
5819 sigma0=sigma_psf*pixel
5820 else if (datatype==
'image_sxr')
then
5822 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
5823 rhessi_rsl=2.3d0/instrument_resolution_factor
5825 pixel=rhessi_rsl*arcsec
5826 sigma0=sigma_psf*pixel
5831 {
do ix^d=ixomin^d,ixomax^d\}
5833 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi1(^d))/(dxi/2.d0)));
5834 ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(
vec_xi2(^d))/(dxi/2.d0)));
5835 ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
5836 if (datatype==
'image_euv')
then
5838 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi
5840 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi
5842 else if (datatype==
'image_sxr')
then
5844 fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
5846 if (fluxsubc>smalldouble)
then
5848 {
do isubc^d=1,nsubc^d\}
5849 ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
5853 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
5854 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
5855 ixpmin1=max(1,ixp1-3)
5856 ixpmax1=min(ixp1+3,numxi1)
5857 ixpmin2=max(1,ixp2-3)
5858 ixpmax2=min(ixp2+3,numxi2)
5859 do ixp1=ixpmin1,ixpmax1
5860 do ixp2=ixpmin2,ixpmax2
5861 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
5862 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
5863 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
5864 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
5865 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
5866 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
5873 deallocate(flux,opacity)
5877 integer,
intent(in) :: igrid,numXI1,numXI2
5878 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
5879 double precision,
intent(in) :: dxI
5881 character(20),
intent(in) :: datatype
5882 double precision,
intent(inout) :: EM(numXI1,numXI2)
5884 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
5885 double precision,
allocatable :: flux(:^D&),Ne(:^D&),opacity(:^D&)
5886 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
5887 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
5888 double precision :: TBsubC,PBsubC
5889 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
5890 double precision :: R_thick,dotp,dvolume,R_occult,Rc
5891 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
5892 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
5893 logical :: sun_back_side,emit
5896 double precision :: logTe
5897 character (30) :: ion
5898 double precision :: lineCent
5899 double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
5900 double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
5902 ^d&ixomin^d=ixmlo^d;
5903 ^d&ixomax^d=ixmhi^d;
5904 ^d&iximin^d=ixglo^d;
5905 ^d&iximax^d=ixghi^d;
5908 arcsec=7.25d5/unit_length
5910 arcsec=7.25d7/unit_length
5913 allocate(flux(ixi^s),opacity(ixi^s))
5914 if (datatype==
'image_euv')
then
5915 if (trim(emission_model)==
'pseudo_current')
then
5917 else if (trim(emission_model)==
'radio_ff')
then
5921 call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
5922 flux(ixo^s)=flux(ixo^s)/instrument_resolution_factor**2
5924 call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
5925 pixel=spacersl*arcsec
5926 sigma0=sigma_psf*pixel
5927 else if (datatype==
'image_sxr')
then
5929 call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
5930 rhessi_rsl=2.3d0/instrument_resolution_factor
5932 pixel=rhessi_rsl*arcsec
5933 sigma0=sigma_psf*pixel
5938 r_thick=r_opt_thick*const_rsun/unit_length
5939 {
do ix^d=ixomin^d,ixomax^d\}
5940 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
5941 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
5943 dxl(2)=x_sph(1)*dx_sph(2)
5944 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
5949 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
5951 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
5953 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
5955 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
5957 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
5959 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
5964 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
5966 dxsubc1=dx_sph(1)/nsubc1
5969 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
5970 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
5971 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
5972 dvolume=dxsubc1*dxsubc2*dxsubc3
5973 if (datatype==
'image_euv')
then
5975 fluxsubc=flux(ix^d)*dvolume*unit_length*1.d2/dxi/dxi
5977 fluxsubc=flux(ix^d)*dvolume*unit_length/dxi/dxi
5979 else if (datatype==
'image_sxr')
then
5981 fluxsubc=flux(ix^d)*dvolume*unit_length**3/area_1au
5984 if (fluxsubc>smalldouble)
then
5987 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
5989 rc=dsqrt(xcent(1)**2+xcent(2)**2)
5994 sun_back_side=.true.
5995 if (dotp<0.d0) sun_back_side=.false.
5997 if (sun_back_side)
then
5999 if (rc>r_thick) emit=.true.
6002 if (xsubc(1)<=r_thick) emit=.false.
6008 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6009 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6010 ixpmin1=max(1,ixp1-3)
6011 ixpmax1=min(ixp1+3,numxi1)
6012 ixpmin2=max(1,ixp2-3)
6013 ixpmax2=min(ixp2+3,numxi2)
6014 do ixp1=ixpmin1,ixpmax1
6015 do ixp2=ixpmin2,ixpmax2
6016 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6017 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6018 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6019 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6020 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6021 em(ixp1,ixp2)=em(ixp1,ixp2)+fluxsubc*factor
6031 deallocate(flux,opacity)
6038 integer,
intent(in) :: igrid,numXI1,numXI2,numWI
6039 double precision,
intent(in) :: xI1(numXI1),xI2(numXI2)
6040 double precision,
intent(in) :: dxI
6042 character(20),
intent(in) :: datatype
6043 double precision,
intent(inout) :: WLB(numXI1,numXI2,numWI)
6045 integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
6046 double precision,
allocatable :: flux(:^D&),Ne(:^D&)
6047 integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
6048 double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC,RsubC
6049 double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
6050 double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
6051 double precision :: R_thick,dotp,dvolume,R_occult
6052 double precision :: xSubC(1:3),dxSubC^D,xCent(1:2),xSubC_car(1:3)
6053 double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
6054 double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6057 ^d&ixomin^d=ixmlo^d;
6058 ^d&ixomax^d=ixmhi^d;
6059 ^d&iximin^d=ixglo^d;
6060 ^d&iximax^d=ixghi^d;
6063 arcsec=7.25d5/unit_length
6065 arcsec=7.25d7/unit_length
6069 if (whitelight_instrument==
'LASCO/C1')
then
6070 lasco_rsl=5.6d0/instrument_resolution_factor
6072 else if (whitelight_instrument==
'LASCO/C2')
then
6073 lasco_rsl=11.4d0/instrument_resolution_factor
6075 else if (whitelight_instrument==
'LASCO/C3')
then
6076 lasco_rsl=56.d0/instrument_resolution_factor
6079 if (r_occultor>1.d0) r_occult=r_occultor
6080 r_occult=r_occult*const_rsun/unit_length
6081 call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,ne)
6084 double precision :: nH_dummy(ixI^S)
6085 call eos%get_ne_nH(ixi^l, ixo^l, ps(igrid)%w, ne, nh_dummy)
6088 pixel=lasco_rsl*arcsec
6089 sigma0=sigma_psf*pixel
6092 r_thick=r_opt_thick*const_rsun/unit_length
6093 {
do ix^d=ixomin^d,ixomax^d\}
6094 x_sph(1:3)=ps(igrid)%x(ix^d,1:3)
6095 dx_sph(1:3)=ps(igrid)%dx(ix^d,1:3)
6097 dxl(2)=x_sph(1)*dx_sph(2)
6098 dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3)
6099 ne0=ne(ix^d)*unit_numberdensity
6104 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6106 nsubc1=max(nsubc1,ceiling(dxl(1)*abs(dotp)/(dxi/2.d0)))
6108 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6110 nsubc2=max(nsubc2,ceiling(dxl(2)*abs(dotp)/(dxi/2.d0)))
6112 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6114 nsubc3=max(nsubc3,ceiling(dxl(3)*abs(dotp)/(dxi/2.d0)))
6119 xsubc(1)=x_sph(1)-half*dx_sph(1)+(isubc1-half)*dx_sph(1)/nsubc1
6121 dxsubc1=dx_sph(1)/nsubc1
6125 xsubc(2)=x_sph(2)-half*dx_sph(2)+(isubc2-half)*dx_sph(2)/nsubc2
6126 dxsubc2=xsubc(1)*dx_sph(2)/nsubc2
6127 dxsubc3=xsubc(1)*dsin(xsubc(2))*dx_sph(3)/nsubc3
6128 dvolume=dxsubc1*dxsubc2*dxsubc3
6131 xsubc(3)=x_sph(3)-half*dx_sph(3)+(isubc3-half)*dx_sph(3)/nsubc3
6133 rc=dsqrt(xcent(1)**2+xcent(2)**2)
6136 if (rc>r_occult)
then
6140 tbsubc=tbsubc*dvolume*unit_length/dxi/dxi
6141 pbsubc=pbsubc*dvolume*unit_length/dxi/dxi
6142 if (tbsubc<1.d-20) emit=.false.
6147 ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
6148 ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
6149 ixpmin1=max(1,ixp1-3)
6150 ixpmax1=min(ixp1+3,numxi1)
6151 ixpmin2=max(1,ixp2-3)
6152 ixpmax2=min(ixp2+3,numxi2)
6153 do ixp1=ixpmin1,ixpmax1
6154 do ixp2=ixpmin2,ixpmax2
6155 xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6156 xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
6157 xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6158 xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
6159 factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
6160 wlb(ixp1,ixp2,1)=wlb(ixp1,ixp2,1)+tbsubc*factor
6161 wlb(ixp1,ixp2,2)=wlb(ixp1,ixp2,2)+pbsubc*factor
6177 double precision,
intent(in) :: Rl
6178 double precision,
intent(inout) :: A,B,C,D
6180 double precision :: sinO,cosO,sinO2,cosO2,tmp
6185 coso=abs(dsqrt(coso2))
6186 tmp=log((1.d0+sino)/coso)
6188 b=-(1.d0-3.d0*sino2-(coso2/sino)*(1.d0+3.d0*sino2)*tmp)/8.d0
6189 c=4.d0/3.d0-coso-coso*coso2/3.d0
6190 d=(5.d0+sino2-(coso2/sino)*(5.d0-sino2)*tmp)/8.d0
6196 double precision,
intent(in) :: Rl,Rin,Ne,A,B,C,D
6197 double precision,
intent(inout) :: fluxTB,fluxPB
6199 double precision :: const,u,Bt,Br,PB,TB,sinchi2
6202 const=1.24878d-25/(1.d0-u/3.d0)
6204 bt=const*(c+u*(d-c))
6205 pb=const*sinchi2*((a+u*(b-a)))
6214 double precision,
intent(in) :: x_sph(1:3)
6215 double precision,
intent(inout) :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
6217 unitv_r(1)=dsin(x_sph(2))*dcos(x_sph(3))
6218 unitv_r(2)=dsin(x_sph(2))*dsin(x_sph(3))
6219 unitv_r(3)=dcos(x_sph(2))
6220 unitv_theta(1)=dcos(x_sph(2))*dcos(x_sph(3))
6221 unitv_theta(2)=dcos(x_sph(2))*dsin(x_sph(3))
6222 unitv_theta(3)=-dsin(x_sph(2))
6223 unitv_phi(1)=-dsin(x_sph(3))
6224 unitv_phi(2)=dcos(x_sph(3))
6229 subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
6233 integer,
intent(in) :: qunit,nXO1,nXO2,nWO
6234 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6235 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6236 double precision,
intent(inout) :: wO(nXO1,nXO2,nWO)
6237 character(20),
intent(in) :: datatype
6239 integer :: nPiece,nP1,nP2,nC1,nC2,nWC
6240 integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
6241 double precision :: uniform_tol
6242 double precision,
allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
6249 if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
6256 if (datatype==
'image_euv' .or. datatype==
'image_sxr')
then
6258 piece_nmax1=block_nx2
6259 piece_nmax2=block_nx3
6261 piece_nmax1=block_nx3
6262 piece_nmax2=block_nx1
6264 piece_nmax1=block_nx1
6265 piece_nmax2=block_nx2
6267 else if (datatype==
'spectrum_euv')
then
6270 piece_nmax2=block_nx1
6272 piece_nmax2=block_nx2
6274 piece_nmax2=block_nx3
6281 loopn1:
do j=piece_nmax1,1,-1
6282 if(mod(nxo1,j)==0)
then
6287 loopn2:
do j=piece_nmax2,1,-1
6288 if(mod(nxo2,j)==0)
then
6301 case(
'EIvtuCCmpi',
'ESvtuCCmpi',
'SIvtuCCmpi',
'WIvtuCCmpi')
6303 allocate(xc(npiece,nc1,nc2,2))
6304 allocate(dxc(npiece,nc1,nc2,2))
6305 allocate(wc(npiece,nc1,nc2,nwo))
6309 ix1=mod(ipc-1,np1)*nc1+ixc1
6310 ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
6311 xc(ipc,ixc1,ixc2,1)=xo1(ix1)
6312 xc(ipc,ixc1,ixc2,2)=xo2(ix2)
6313 dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
6314 dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
6316 wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
6323 deallocate(xc,dxc,wc)
6324 case(
'EIvtiCCmpi',
'ESvtiCCmpi',
'SIvtiCCmpi',
'WIvtiCCmpi')
6326 (maxval(abs(dxo1(:)-dxo1(1)))>uniform_tol*max(one,abs(dxo1(1))) .or. &
6327 maxval(abs(dxo2(:)-dxo2(1)))>uniform_tol*max(one,abs(dxo2(1)))))
then
6328 call mpistop(
"vti needs uniform dat-resolution image grids")
6330 call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
6333 call mpistop(
"Error in synthesize emission: Unknown convert_type")
6339 subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
6343 integer,
intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
6344 double precision,
intent(in) :: xO1(nXO1),xO2(nxO2)
6345 double precision,
intent(in) :: dxO1(nxO1),dxO2(nxO2)
6346 double precision,
intent(in) :: wO(nXO1,nXO2,nWO)
6348 double precision :: origin(1:3), spacing(1:3)
6349 integer :: wholeExtent(1:6)
6351 integer :: ixC1,ixC2
6355 character (70) :: subname,wname,vname,nameL,nameS
6356 character (len=std_len) :: filename
6357 logical :: sph_datres_no_doppler
6360 origin(1)=xo1(1)-0.5d0*dxo1(1)
6361 origin(2)=xo2(1)-0.5d0*dxo2(1)
6372 inquire(qunit,opened=fileopen)
6373 if(.not.fileopen)
then
6378 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vti"
6380 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vti"
6386 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6390 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6391 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="ImageData"'
6392 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6393 write(qunit,
'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')
' <ImageData Origin="',&
6394 origin,
'" WholeExtent="',wholeextent,
'" Spacing="',spacing,
'">'
6396 write(qunit,
'(a)')
'<FieldData>'
6397 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6398 'NumberOfTuples="1" format="ascii">'
6400 write(qunit,
'(a)')
'</DataArray>'
6401 write(qunit,
'(a)')
'</FieldData>'
6403 write(qunit,
'(a,6(i10),a)')
'<Piece Extent="',wholeextent,
'">'
6404 write(qunit,
'(a)')
'<CellData>'
6415 if (trim(
emission_model)==
'pseudo_current' .and. iw==1) vname=
'pseudo_current'
6416 if (trim(
emission_model)==
'radio_ff' .and. iw==1) vname=
'radio_brightness_temperature'
6418 if (iw==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6424 ((
dat_resolution .and. ((sph_datres_no_doppler .and. iw==2) .or. &
6425 ((.not. sph_datres_no_doppler) .and. iw==3))) .or. &
6439 vname=
'absorption_fraction'
6450 if (iw==1)
write(vname,
'(a)')
'B'
6451 if (iw==2)
write(vname,
'(a)')
'pB'
6459 write(qunit,
'(a,a,a)')&
6460 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6461 write(qunit,
'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=1,nxo1),ixc2=1,nxo2)
6462 write(qunit,
'(a)')
'</DataArray>'
6464 write(qunit,
'(a)')
'</CellData>'
6465 write(qunit,
'(a)')
'</Piece>'
6467 write(qunit,
'(a)')
'</ImageData>'
6468 write(qunit,
'(a)')
'</VTKFile>'
6478 integer,
intent(in) :: qunit
6479 integer,
intent(in) :: nPiece,nC1,nC2,nWC
6480 double precision,
intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
6481 double precision,
intent(in) :: wC(nPiece,nC1,nC2,nWC)
6482 character(20),
intent(in) :: datatype
6485 double precision :: xP(nPiece,nC1+1,nC2+1,2)
6488 character (70) :: subname,wname,vname,nameL,nameS
6489 character (len=std_len) :: filename
6490 integer :: ixC1,ixC2,ixP,ix1,ix2,j
6491 integer :: nc,np,icel,VTK_type
6492 logical :: sph_datres_no_doppler
6503 if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
6504 if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
6505 if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
6506 if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
6511 inquire(qunit,opened=fileopen)
6512 if(.not.fileopen)
then
6516 if (datatype==
'image_euv')
then
6517 write(filename,
'(a,i4.4,a)') trim(
filename_euv),filenr,
".vtu"
6518 else if (datatype==
'image_sxr')
then
6519 write(filename,
'(a,i4.4,a)') trim(
filename_sxr),filenr,
".vtu"
6520 else if (datatype==
'image_whitelight')
then
6522 else if (datatype==
'spectrum_euv')
then
6525 open(qunit,file=filename,status=
'unknown',form=
'formatted')
6528 write(qunit,
'(a)')
'<?xml version="1.0"?>'
6529 write(qunit,
'(a)',advance=
'no')
'<VTKFile type="UnstructuredGrid"'
6530 write(qunit,
'(a)')
' version="0.1" byte_order="LittleEndian">'
6531 write(qunit,
'(a)')
'<UnstructuredGrid>'
6532 write(qunit,
'(a)')
'<FieldData>'
6533 write(qunit,
'(2a)')
'<DataArray type="Float32" Name="TIME" ',&
6534 'NumberOfTuples="1" format="ascii">'
6536 write(qunit,
'(a)')
'</DataArray>'
6537 write(qunit,
'(a)')
'</FieldData>'
6539 write(qunit,
'(a,i7,a,i7,a)') &
6540 '<Piece NumberOfPoints="',np,
'" NumberOfCells="',nc,
'">'
6541 write(qunit,
'(a)')
'<CellData>'
6543 if (datatype==
'image_euv')
then
6552 if (trim(
emission_model)==
'pseudo_current') vname=
'pseudo_current'
6553 if (trim(
emission_model)==
'radio_ff') vname=
'radio_brightness_temperature'
6556 if (j==2 .and.
dat_resolution .and. (.not. sph_datres_no_doppler) .and. &
6562 ((
dat_resolution .and. ((sph_datres_no_doppler .and. j==2) .or. &
6563 ((.not. sph_datres_no_doppler) .and. j==3))) .or. &
6577 vname=
'absorption_fraction'
6579 else if (datatype==
'image_sxr')
then
6587 else if (datatype==
'image_whitelight')
then
6588 write(vname,
'(a)')
'whitelight'
6589 else if (datatype==
'spectrum_euv')
then
6596 write(qunit,
'(a,a,a)')&
6597 '<DataArray type="Float64" Name="',trim(vname),
'" format="ascii">'
6598 write(qunit,
'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
6599 write(qunit,
'(a)')
'</DataArray>'
6601 write(qunit,
'(a)')
'</CellData>'
6602 write(qunit,
'(a)')
'<Points>'
6603 write(qunit,
'(a)')
'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
6608 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6610 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6612 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6616 write(qunit,
'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
6618 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
6620 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6623 write(qunit,
'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
6627 write(qunit,
'(a)')
'</DataArray>'
6628 write(qunit,
'(a)')
'</Points>'
6630 write(qunit,
'(a)')
'<Cells>'
6631 write(qunit,
'(a)')
'<DataArray type="Int32" Name="connectivity" format="ascii">'
6634 write(qunit,
'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
6635 ix1-1+ix2*np1,ix1+ix2*np1
6638 write(qunit,
'(a)')
'</DataArray>'
6640 write(qunit,
'(a)')
'<DataArray type="Int32" Name="offsets" format="ascii">'
6642 write(qunit,
'(i7)') icel*(2**2)
6644 write(qunit,
'(a)')
'</DataArray>'
6646 write(qunit,
'(a)')
'<DataArray type="Int32" Name="types" format="ascii">'
6650 write(qunit,
'(i2)') vtk_type
6652 write(qunit,
'(a)')
'</DataArray>'
6653 write(qunit,
'(a)')
'</Cells>'
6654 write(qunit,
'(a)')
'</Piece>'
6656 write(qunit,
'(a)')
'</UnstructuredGrid>'
6657 write(qunit,
'(a)')
'</VTKFile>'
6663 double precision,
intent(in) :: vec1(1:3),vec2(1:3)
6664 double precision,
intent(out) :: res
6666 res=vec1(1)*vec2(1)+vec1(2)*vec2(2)+vec1(3)*vec2(3)
6671 double precision,
intent(in) :: vec_in1(1:3),vec_in2(1:3)
6672 double precision,
intent(out) :: vec_out(1:3)
6674 vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
6675 vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
6676 vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
6682 double precision :: LOS_psi
6683 double precision :: vec_car(1:3),vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6684 double precision :: vec_LOS_sph(1:3),vec_xI1_sph(1:3),vec_xI2_sph(1:3)
6688 vec_los(2)=dpi*los_theta/180.d0
6699 if (los_theta==zero)
then
6701 vec_temp1(2)=dpi/2.d0
6702 vec_temp1(3)=dpi*los_phi/180.d0
6716 los_psi=dpi*image_rotate/180.d0
6717 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6718 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6729 vec_los_sph(2:3)=vec_los_sph(2:3)*180.d0/dpi
6730 vec_xi1_sph(2:3)=vec_xi1_sph(2:3)*180.d0/dpi
6731 vec_xi2_sph(2:3)=vec_xi2_sph(2:3)*180.d0/dpi
6733 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' ray direction (spherical): [',vec_los_sph(1),vec_los_sph(2),vec_los_sph(3),
']'
6734 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' xI1 direction (spherical): [',vec_xi1_sph(1),vec_xi1_sph(2),vec_xi1_sph(3),
']'
6735 if (mype==0)
write(*,
'(a,f3.1,f6.1,f6.1,a)')
' xI2 direction (spherical): [',vec_xi2_sph(1),vec_xi2_sph(2),vec_xi2_sph(3),
']'
6741 double precision,
intent(in) :: vec_sph(1:3)
6742 double precision,
intent(inout) :: vec_car(1:3)
6744 vec_car(1)=vec_sph(1)*dsin(vec_sph(2))*dcos(vec_sph(3))
6745 vec_car(2)=vec_sph(1)*dsin(vec_sph(2))*dsin(vec_sph(3))
6746 vec_car(3)=vec_sph(1)*dcos(vec_sph(2))
6752 double precision,
intent(in) :: vec_car(1:3)
6753 double precision,
intent(inout) :: vec_sph(1:3)
6755 vec_sph(1)=dsqrt(vec_car(1)**2+vec_car(2)**2+vec_car(3)**2)
6756 vec_sph(2)=dacos(vec_car(3)/vec_sph(1))
6757 vec_sph(3)=atan2(vec_car(2),vec_car(1))
6763 double precision :: LOS_psi
6764 double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
6767 vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
6768 vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
6769 vec_los(3)=-cos(dpi*los_theta/180.d0)
6775 if (los_theta==zero)
then
6776 vec_xi1(1)=cos(dpi*los_phi/180.d0)
6777 vec_xi1(2)=sin(dpi*los_phi/180.d0)
6785 los_psi=dpi*image_rotate/180.d0
6786 vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
6787 vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
6801 double precision,
intent(in) :: x_3D_sph(1:3)
6802 double precision,
intent(inout) :: x_image(1:2)
6803 double precision :: res,res_origin
6804 double precision :: x_3D(1:3)
6815 double precision,
intent(in) :: x_3D(1:3)
6816 double precision,
intent(inout) :: x_image(1:2)
6817 double precision :: res,res_origin
6821 x_image(1)=res-res_origin
6824 x_image(2)=res-res_origin
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter const_rsun
double precision, parameter kb_cgs
Boltzmann constant in cgs.
double precision, parameter half
double precision, parameter one
double precision, parameter dpi
Pi.
double precision, parameter zero
some frequently used numbers
double precision, parameter smalldouble
double precision, parameter mp_cgs
Proton mass in cgs.
double precision, parameter const_c
universal constants as specified in cgs units
Equation of state for AMRVAC, handled through a single eos_container object.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cartesian
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
integer spectrum_wl
wave length for spectrum
integer ixghi
Upper index of grid block arrays.
logical activate_unit_arcsec
use arcsec as length unit of images/spectra
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
logical output_absorption_fraction
output absorption fraction for thick/thin EUV synthesis when available
double precision radio_beam_fwhm
Gaussian radio beam full width at half maximum in arcsec.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
logical instrument_postprocess
Post-process dat-resolution EUV images onto the instrument pixel grid.
character(len=std_len) filename_whitelight
Base file name for synthetic white light.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter rpxmin
double precision unit_length
Physical scaling factor for length.
double precision location_slit
location of the slit
double precision time_convert_factor
Conversion factor for time unit.
integer icomm
The MPI communicator.
character(len=std_len) whitelight_instrument
white light observation instrument
integer mype
The rank of the current MPI task.
double precision radio_frequency
Observing frequency for radio free-free synthesis in Hz.
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision, dimension(:), allocatable, parameter d
logical slab
Cartesian geometry or not.
double precision radio_beam_pixel_size
Output pixel size for radio beam post-processing in arcsec; <=0 uses FWHM/3.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
logical dat_resolution
resolution of the images
double precision r_occultor
the white light emission below it (unit=Rsun) is not visible
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
integer npe
The number of MPI tasks.
logical output_tau
output optical-depth map for synthetic emission when available
double precision, dimension(^nd) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
double precision unit_velocity
Physical scaling factor for velocity.
integer radsyn_segment_batch_factor
Maximum ray segments per pixel batch, as a factor of radsyn_pixel_batch; <=0 uses memory budget....
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision los_theta
direction of the line of sight (LOS)
character(len=std_len) radiation_transfer
Synthetic emission transfer mode: thin or thick.
double precision spectrum_window_max
integer wavelength
wavelength for output
integer, dimension(ndim) stretch_type
What kind of stretching is used per dimension.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
integer, parameter rpxmax
integer radsyn_pixel_batch
Number of image pixels processed in one ray-segment MPI batch.
logical radsyn_verbose
Print synthetic-emission ray-tracing profiling counters.
logical big_image
big image
double precision instrument_resolution_factor
times for enhancing spatial resolution for EUV image/spectra
double precision radsyn_segment_memory_mb
Approximate per-rank temporary memory budget, in MiB, for automatic ray-segment batch sizing.
double precision spectrum_window_min
spectral window
integer refine_max_level
Maximal number of AMR levels.
character(len=std_len) ray_method
Synthetic emission ray traversal method.
integer direction_slit
direction of the slit (for dat resolution only)
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
character(len=std_len) emission_model
Synthetic emission physical model selector.
integer, dimension(:,:), allocatable node
integer radsyn_segment_comm_factor
Maximum ray segments per segmented MPI all-to-all round, as a factor of radsyn_pixel_batch.
integer, parameter ixglo
Lower index of grid block arrays (always 1)
This module defines the procedures of a physics module. It contains function pointers for the various...
double precision, dimension(1:3) vec_los
subroutine get_goes_flux_grid(ixil, ixol, w, x, dv, xboxl, xbl, fl, eflux_grid)
subroutine integrate_spectra_cartesian(igrid, wl, dwlg, xs, dxsg, spectra, numwl, numxs, fl)
subroutine sph_cart_to_coord(pos, sph)
subroutine get_spectrum_datresol(qunit, datatype, fl)
double precision, dimension(1:101) f_304
subroutine get_sph_intersection_datresol_spacing(dxi)
double precision, dimension(1:101) f_193
double precision, dimension(1:60) f_264
subroutine normalize_euv_doppler(ni1, ni2, euv, dpl, unitv)
subroutine get_sph_intersection_image_bounds(ximin1, ximax1, ximin2, ximax2)
double precision, dimension(1:60) f_263
integer function sph_locate_index(value, faces, imin, imax)
subroutine get_euv_image(qunit, fl)
double precision, dimension(1:60) t_eis1
subroutine postprocess_radio_beam_image(nsrc1, nsrc2, xsrc1, xsrc2, dxsrc1, dxsrc2, bright, nout1, nout2, xout1, xout2, dxout1, dxout2, wout, numwout, tau, brightthin)
subroutine get_sxr(ixil, ixol, w, x, fl, flux, el, eu)
integer function radsyn_euv_num_outputs(has_doppler, has_thick)
subroutine collect_euv_cart_dda_segments(ixil, ixol, source, opacity, sourcev, pixel_id, ray_origin, xface1, xface2, xface3, t_enter, t_exit, segments, nseg, capacity)
subroutine get_unit_vector_spherical(x_sph, unitv_r, unitv_theta, unitv_phi)
subroutine collect_euv_sph_intersection_segments(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface, thetaface, phiface, rface2, theta_cos, phi_sin, phi_cos, segments, nseg, capacity)
double precision, dimension(1:101) f_131
recursive subroutine quicksort_segment_indices(segments, idx, ilo, ihi)
subroutine get_line_info(wl, ion, mass, logte, line_center, spatial_px, spectral_px, sigma_psf, width_slit)
double precision, dimension(1:60) f_255
subroutine sph_block_pixel_range(rface, thetaface, phiface, ixol, nxi1, nxi2, xi1, xi2, dxi, ixpmin1, ixpmax1, ixpmin2, ixpmax2, has_pixels)
subroutine cart_dda_advance_axis(ray_origin_axis, ray_dir_axis, faces, imin, imax, idx, step, tmax, done)
subroutine get_pseudo_current(igrid, ixil, ixol, w, source)
logical function segment_is_valid(segments, is, nvars)
double precision function transfer_attenuation(tau)
double precision, dimension(1:101) t_aia
double precision function exp_clamped(argument)
subroutine write_image_vtucc(qunit, xc, wc, dxc, npiece, nc1, nc2, nwc, datatype)
logical function sph_segment_visible(pos, ximg1, ximg2)
subroutine get_cor_image(x_3d, x_image)
subroutine get_thomson_parameters(rl, a, b, c, d)
subroutine integrate_euv_datresol(igrid, nxif1, nxif2, xif1, xif2, dxif1, dxif2, fl, euv, dpl)
subroutine cart_dda_block_pixel_range(box_min, box_max, nxif1, nxif2, xif1, xif2, ixpmin1, ixpmax1, ixpmin2, ixpmax2, has_pixels)
subroutine sph_add_t_fixed(tvals, nt, t)
double precision, dimension(1:60) t_eis2
subroutine sph_add_theta_intersections(ray_origin, ray_dir, thetaface, tvals, nt, capacity)
subroutine integrate_euv_cart_dda_thick_datresol(nxif1, nxif2, xif1, xif2, fl, euv, dpl, tau, euvthin)
double precision, dimension(1:101) f_171
double precision, dimension(1:101) f_94
subroutine integrate_spectra_datresol(igrid, wl, dwl, spectra, numwl, numxs, dir_loc, fl)
subroutine acc_euv_cart_dda(ixil, ixol, source, sourcev, ray_origin, xface1, xface2, xface3, t_enter, t_exit, euvp, dplp)
subroutine sph_add_phi_intersection(ray_origin, ray_dir, phiface, tvals, nt, capacity)
double precision, dimension(1:3) vec_xi1
subroutine append_cart_dda_segment(segments, nseg, capacity, pixel_id, tseg, jds, kds, jvds)
subroutine radsyn_get_segment_batch_limits(pixel_batch_target, segment_batch_target, segment_comm_target)
subroutine get_spectrum(qunit, datatype, fl)
subroutine ray_box_intersection_cart(ray_origin, ray_dir, box_min, box_max, hit, t_enter, t_exit)
subroutine cartesian_to_spherical(vec_car, vec_sph)
subroutine get_cor_image_spherical(x_3d_sph, x_image)
subroutine dot_product_loc(vec1, vec2, res)
subroutine integrate_emission_spherical(igrid, numxi1, numxi2, xi1, xi2, dxi, fl, datatype, em)
subroutine get_image(qunit, datatype, fl)
subroutine integrate_emission_cartesian(igrid, numxi1, numxi2, xi1, xi2, dxi, fl, datatype, em)
logical function radsyn_euv_has_doppler_output()
integer function sph_locate_index_desc(value, faces, imin, imax)
subroutine init_vectors_spherical()
subroutine insertion_sort_segment_indices(segments, idx, ilo, ihi)
subroutine write_image_vticc(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, nc1, nc2)
subroutine get_sxr_image(qunit, fl)
subroutine check_synthetic_emission_options(datatype)
subroutine fill_euv_absorption_fraction(ni1, ni2, euv, euvthin, smallflux, absorption, cap_to_one)
subroutine sph_sort_unique_t(tvals, nt)
subroutine init_vectors_cartesian()
double precision, dimension(1:60) f_192
subroutine sort_segment_indices_near_to_far(segments, idx, nidx)
integer function cart_dda_locate_index(pos, faces, imin, imax)
subroutine pack_euv_image_outputs(ni1, ni2, euv, wi, smallflux, has_doppler, has_thick, dpl, tau, euvthin, cap_absorption)
integer function segment_pixel_owner(pixel_id)
subroutine get_whitelight_thomson(rl, rin, ne, a, b, c, d, fluxtb, fluxpb)
subroutine get_image_datresol(qunit, datatype, fl)
subroutine collect_euv_sph_dda_interval(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface2, theta_cos, phiface, phi_sin, phi_cos, t_enter, t_exit, segments, nseg, capacity, ok)
subroutine get_goes_sxr_flux(xboxl, fl, eflux)
double precision function interpolate_response_value(temperature, t_table, f_table, n_table, log_temperature, log_response)
subroutine integrate_whitelight_spherical(igrid, numxi1, numxi2, numwi, xi1, xi2, dxi, fl, datatype, wlb)
double precision, dimension(1:3) vec_xi2
double precision, dimension(1:41) f_1354
subroutine cross_product_loc(vec_in1, vec_in2, vec_out)
character(len=std_len) ray_method_active
subroutine integrate_sxr_datresol(igrid, nxif1, nxif2, xif1, xif2, dxif1, dxif2, fl, sxr)
subroutine apply_temperature_response(ixil, ixol, te, flux, t_table, f_table, n_table, log_temperature, log_response)
double precision, dimension(1:101) f_335
subroutine sph_try_theta_exit_candidate(t, theta_face_cos, ray_origin, tnow, texit, epsray, tnext, found)
subroutine get_radio_ff_source_opacity(ixil, ixol, w, x, fl, source, kappa)
subroutine sph_next_cell_exit(ray_origin, rface2, theta_cos, phiface, phi_sin, phi_cos, ixol, ix1, ix2, ix3, tnow, texit, epsray, tnext, found)
subroutine integrate_euv_thick_datresol(nxif1, nxif2, fl, euv, dpl, tau, euvthin)
subroutine output_data(qunit, xo1, xo2, dxo1, dxo2, wo, nxo1, nxo2, nwo, datatype)
subroutine integrate_euv_sph_intersection_thick(numxi1, numxi2, xi1, xi2, dxi, fl, euv, tau, euvthin)
subroutine sph_try_phi_exit_candidate(t, phi_face_sin, phi_face_cos, ray_origin, tnow, texit, epsray, tnext, found)
double precision, dimension(1:41) t_iris
subroutine sph_locate_cell_fast(pos, rface2, theta_cos, phiface, ixol, ix1, ix2, ix3, inside)
subroutine get_euv_spectrum(qunit, fl)
double precision, dimension(1:101) f_211
subroutine build_sph_intersection_faces(ixil, ixol, x, dx, rface, thetaface, phiface)
subroutine collect_euv_sph_dda_segments(ixil, ixol, source, opacity, pixel_id, ray_origin, ximg1, ximg2, rface, thetaface, phiface, rface2, theta_cos, phi_sin, phi_cos, segments, nseg, capacity, fallback)
subroutine cart_dda_init_axis(ray_origin_axis, ray_dir_axis, faces, imin, imax, idx, step, tmax)
subroutine get_euv_hhe_opacity(wl, ixil, ixol, w, x, fl, kappa)
subroutine integrate_euv_sph_intersection_thin(numxi1, numxi2, xi1, xi2, dxi, fl, em)
subroutine integrate_euv_cart_dda_datresol(nxif1, nxif2, xif1, xif2, fl, euv, dpl)
subroutine get_whitelight_image(qunit, fl)
subroutine get_euv(wl, ixil, ixol, w, x, fl, flux)
subroutine sph_add_sphere_intersections(ray_origin, ray_dir, rface, tvals, nt, capacity)
subroutine build_cart_dda_faces(ixil, ixol, x, dx, xface1, xface2, xface3)
subroutine integrate_transfer_step_first_order(emissivity, opacity, path_length, intensity, tau)
double precision function pow10_clamped(exponent)
subroutine acc_euv_sph_intersection(ixil, ixol, source, ray_origin, ximg1, ximg2, rface, thetaface, phiface, euvp)
subroutine sph_locate_cell(pos, rface, thetaface, phiface, ixol, ix1, ix2, ix3, inside)
subroutine sph_try_exit_candidate(t, tnow, texit, epsray, tnext, found)
subroutine spherical_to_cartesian(vec_sph, vec_car)
subroutine sph_add_t(tvals, nt, capacity, t)
subroutine postprocess_euv_instrument_image(nsrc1, nsrc2, xsrc1, xsrc2, dxsrc1, dxsrc2, euv, dpl, nout1, nout2, xout1, xout2, dxout1, dxout2, wout, numwout, tau, euvthin)