MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_thermal_emission.t
Go to the documentation of this file.
1 ! module mod_thermal_emission -- synthesize emission flux of some
2 ! thermal lines
3 ! EUV lines database:
4 ! 'He_II_304' 'Fe_IX_171' 'Fe_XXIV_193' 'Fe_XIV_211' 'Fe_XVI_335'
5 ! 'Fe_XVIII_94' 'Fe_XXI_131'
6 ! subroutines:
7 ! get_EUV: get local EUV emission intensity (for 1d, 2d and 3d)
8 ! get_SXR: get local Soft X-ray emission intensity (for 1d, 2d and 3d)
9 
12  use mod_physics
13 
14  implicit none
15 
16  integer :: n_aia
17  double precision :: t_aia(1:101)
18  double precision :: f_94(1:101),f_131(1:101),f_171(1:101)
19  double precision :: f_193(1:101),f_211(1:101),f_304(1:101)
20  double precision :: f_335(1:101)
21  integer :: n_iris
22  double precision :: t_iris(1:41)
23  double precision :: f_1354(1:41)
24  integer :: n_eis
25  double precision :: t_eis1(1:60),t_eis2(1:60)
26  double precision :: f_263(1:60),f_264(1:60),f_192(1:60),f_255(1:60)
27 
28 
29  double precision :: vec_xi1(1:3),vec_xi2(1:3),vec_los(1:3)
30 
31  data n_aia / 101 /
32 
33  data t_aia / 4. , 4.05, 4.1, 4.15, 4.2, 4.25, 4.3, 4.35, &
34  4.4, 4.45, 4.5, 4.55, 4.6, 4.65, 4.7, 4.75, &
35  4.8, 4.85, 4.9, 4.95, 5. , 5.05, 5.1, 5.15, &
36  5.2, 5.25, 5.3, 5.35, 5.4, 5.45, 5.5, 5.55, &
37  5.6, 5.65, 5.7, 5.75, 5.8, 5.85, 5.9, 5.95, &
38  6. , 6.05, 6.1, 6.15, 6.2, 6.25, 6.3, 6.35, &
39  6.4, 6.45, 6.5, 6.55, 6.6, 6.65, 6.7, 6.75, &
40  6.8, 6.85, 6.9, 6.95, 7. , 7.05, 7.1, 7.15, &
41  7.2, 7.25, 7.3, 7.35, 7.4, 7.45, 7.5, 7.55, &
42  7.6, 7.65, 7.7, 7.75, 7.8, 7.85, 7.9, 7.95, &
43  8. , 8.05, 8.1, 8.15, 8.2, 8.25, 8.3, 8.35, &
44  8.4, 8.45, 8.5, 8.55, 8.6, 8.65, 8.7, 8.75, &
45  8.8, 8.85, 8.9, 8.95, 9. /
46 
47  data f_94 / 4.25022959d-37, 4.35880298d-36, 3.57054296d-35, 2.18175426d-34, &
48  8.97592571d-34, 2.68512961d-33, 7.49559346d-33, 2.11603751d-32, &
49  5.39752853d-32, 1.02935904d-31, 1.33822307d-31, 1.40884290d-31, &
50  1.54933156d-31, 2.07543102d-31, 3.42026227d-31, 6.31171444d-31, &
51  1.16559416d-30, 1.95360497d-30, 2.77818735d-30, 3.43552578d-30, &
52  4.04061803d-30, 4.75470982d-30, 5.65553769d-30, 6.70595782d-30, &
53  7.80680354d-30, 8.93247715d-30, 1.02618156d-29, 1.25979030d-29, &
54  1.88526483d-29, 3.62448572d-29, 7.50553279d-29, 1.42337571d-28, &
55  2.37912813d-28, 3.55232305d-28, 4.84985757d-28, 6.20662827d-28, &
56  7.66193687d-28, 9.30403645d-28, 1.10519802d-27, 1.25786927d-27, &
57  1.34362634d-27, 1.33185242d-27, 1.22302081d-27, 1.05677973d-27, &
58  9.23064720d-28, 8.78570994d-28, 8.02397416d-28, 5.87681142d-28, &
59  3.82272695d-28, 3.11492649d-28, 3.85736090d-28, 5.98893519d-28, &
60  9.57553548d-28, 1.46650267d-27, 2.10365847d-27, 2.79406671d-27, &
61  3.39420087d-27, 3.71077520d-27, 3.57296767d-27, 2.95114380d-27, &
62  2.02913103d-27, 1.13361825d-27, 5.13405629d-28, 2.01305089d-28, &
63  8.15781482d-29, 4.28366817d-29, 3.08701543d-29, 2.68693906d-29, &
64  2.51764203d-29, 2.41773103d-29, 2.33996083d-29, 2.26997246d-29, &
65  2.20316143d-29, 2.13810001d-29, 2.07424438d-29, 2.01149189d-29, &
66  1.94980213d-29, 1.88917920d-29, 1.82963583d-29, 1.77116920d-29, &
67  1.71374392d-29, 1.65740593d-29, 1.60214447d-29, 1.54803205d-29, &
68  1.49510777d-29, 1.44346818d-29, 1.39322305d-29, 1.34441897d-29, &
69  1.29713709d-29, 1.25132618d-29, 1.20686068d-29, 1.14226584d-29, &
70  1.09866413d-29, 1.05635524d-29, 1.01532444d-29, 9.75577134d-30, &
71  9.37102736d-30, 8.99873335d-30, 8.63860172d-30, 8.29051944d-30, &
72  7.95414793d-30 /
73 
74  data f_131 / 3.18403601d-37, 3.22254703d-36, 2.61657920d-35, &
75  1.59575286d-34, 6.65779556d-34, 2.07015132d-33, &
76  6.05768615d-33, 1.76074833d-32, 4.52633001d-32, &
77  8.57121883d-32, 1.09184271d-31, 1.10207963d-31, &
78  1.11371658d-31, 1.29105226d-31, 1.80385897d-31, &
79  3.27295431d-31, 8.92002136d-31, 3.15214579d-30, &
80  9.73440787d-30, 2.22709702d-29, 4.01788984d-29, &
81  6.27471832d-29, 8.91764995d-29, 1.18725647d-28, &
82  1.52888040d-28, 2.05082946d-28, 3.47651873d-28, &
83  8.80482184d-28, 2.66533063d-27, 7.05805149d-27, &
84  1.46072515d-26, 2.45282476d-26, 3.55303726d-26, &
85  4.59075911d-26, 5.36503515d-26, 5.68444094d-26, &
86  5.47222296d-26, 4.81119761d-26, 3.85959059d-26, &
87  2.80383406d-26, 1.83977650d-26, 1.11182849d-26, &
88  6.50748885d-27, 3.96843481d-27, 2.61876319d-27, &
89  1.85525324d-27, 1.39717024d-27, 1.11504283d-27, &
90  9.38169611d-28, 8.24801234d-28, 7.43331919d-28, &
91  6.74537063d-28, 6.14495760d-28, 5.70805277d-28, &
92  5.61219786d-28, 6.31981777d-28, 9.19747307d-28, &
93  1.76795732d-27, 3.77985446d-27, 7.43166191d-27, &
94  1.19785603d-26, 1.48234676d-26, 1.36673114d-26, &
95  9.61047146d-27, 5.61209353d-27, 3.04779780d-27, &
96  1.69378976d-27, 1.02113491d-27, 6.82223774d-28, &
97  5.02099099d-28, 3.99377760d-28, 3.36279037d-28, &
98  2.94767378d-28, 2.65740865d-28, 2.44396277d-28, &
99  2.28003967d-28, 2.14941419d-28, 2.04178995d-28, &
100  1.95031045d-28, 1.87011994d-28, 1.79777869d-28, &
101  1.73093957d-28, 1.66795789d-28, 1.60785455d-28, &
102  1.55002399d-28, 1.49418229d-28, 1.44022426d-28, &
103  1.38807103d-28, 1.33772767d-28, 1.28908404d-28, &
104  1.24196208d-28, 1.17437501d-28, 1.12854330d-28, &
105  1.08410498d-28, 1.04112003d-28, 9.99529904d-29, &
106  9.59358806d-29, 9.20512291d-29, 8.83009123d-29, &
107  8.46817043d-29, 8.11921928d-29 /
108 
109  data f_171 / 2.98015581d-42, 1.24696230d-40, 3.37614652d-39, 5.64103034d-38, &
110  5.20550266d-37, 2.77785939d-36, 1.16283616d-35, 6.50007689d-35, &
111  9.96177399d-34, 1.89586076d-32, 2.10982799d-31, 1.36946479d-30, &
112  6.27396553d-30, 2.29955134d-29, 7.13430211d-29, 1.91024282d-28, &
113  4.35358848d-28, 7.94807808d-28, 1.07431875d-27, 1.08399488d-27, &
114  9.16212938d-28, 7.34715770d-28, 6.59246382d-28, 9.13541375d-28, &
115  2.05939035d-27, 5.08206555d-27, 1.10148083d-26, 2.01884662d-26, &
116  3.13578384d-26, 4.14367719d-26, 5.36067711d-26, 8.74170213d-26, &
117  1.64161233d-25, 2.94587860d-25, 4.76298332d-25, 6.91765639d-25, &
118  9.08825111d-25, 1.08496183d-24, 1.17440114d-24, 1.13943939d-24, &
119  9.71696981d-25, 7.09593688d-25, 4.31376399d-25, 2.12708486d-25, &
120  8.47429567d-26, 3.17608104d-26, 1.95898842d-26, 1.98064242d-26, &
121  1.67706555d-26, 8.99126003d-27, 3.29773878d-27, 1.28896127d-27, &
122  8.51169698d-28, 7.53520167d-28, 6.18268143d-28, 4.30034650d-28, &
123  2.78152409d-28, 1.95437088d-28, 1.65896278d-28, 1.68740181d-28, &
124  1.76054383d-28, 1.63978419d-28, 1.32880591d-28, 1.00833205d-28, &
125  7.82252806d-29, 6.36181741d-29, 5.34633869d-29, 4.58013864d-29, &
126  3.97833422d-29, 3.49414760d-29, 3.09790940d-29, 2.76786227d-29, &
127  2.48806269d-29, 2.24823367d-29, 2.04016653d-29, 1.85977413d-29, &
128  1.70367499d-29, 1.56966125d-29, 1.45570643d-29, 1.35964565d-29, &
129  1.27879263d-29, 1.21016980d-29, 1.15132499d-29, 1.09959628d-29, &
130  1.05307482d-29, 1.01040261d-29, 9.70657096d-30, 9.33214234d-30, &
131  8.97689427d-30, 8.63761192d-30, 8.31149879d-30, 7.85162401d-30, &
132  7.53828281d-30, 7.23559452d-30, 6.94341530d-30, 6.66137038d-30, &
133  6.38929156d-30, 6.12669083d-30, 5.87346434d-30, 5.62943622d-30, &
134  5.39435202d-30 /
135 
136  data f_193 / 6.40066486d-32, 4.92737300d-31, 2.95342934d-30, 1.28061594d-29, &
137  3.47747667d-29, 5.88554792d-29, 7.72171179d-29, 9.75609282d-29, &
138  1.34318963d-28, 1.96252638d-28, 2.70163878d-28, 3.63192965d-28, &
139  5.28087341d-28, 8.37821446d-28, 1.39089159d-27, 2.31749718d-27, &
140  3.77510689d-27, 5.85198594d-27, 8.26021568d-27, 1.04870405d-26, &
141  1.25209374d-26, 1.47406787d-26, 1.77174067d-26, 2.24098537d-26, &
142  3.05926105d-26, 4.50018853d-26, 6.84720216d-26, 1.00595861d-25, &
143  1.30759222d-25, 1.36481773d-25, 1.15943558d-25, 1.01467304d-25, &
144  1.04092532d-25, 1.15071251d-25, 1.27416033d-25, 1.38463476d-25, &
145  1.47882726d-25, 1.57041238d-25, 1.69786224d-25, 1.94970397d-25, &
146  2.50332918d-25, 3.58321431d-25, 5.18061550d-25, 6.60405549d-25, &
147  6.64085365d-25, 4.83825816d-25, 2.40545020d-25, 8.59534098d-26, &
148  2.90920638d-26, 1.33204845d-26, 9.03933926d-27, 7.78910836d-27, &
149  7.29342321d-27, 7.40267022d-27, 8.05279981d-27, 8.13829291d-27, &
150  6.92634262d-27, 5.12521880d-27, 3.59527615d-27, 2.69617560d-27, &
151  2.84432713d-27, 5.06697306d-27, 1.01281903d-26, 1.63526978d-26, &
152  2.06759342d-26, 2.19482312d-26, 2.10050611d-26, 1.89837248d-26, &
153  1.66347131d-26, 1.43071097d-26, 1.21518419d-26, 1.02078343d-26, &
154  8.46936184d-27, 6.93015742d-27, 5.56973237d-27, 4.38951754d-27, &
155  3.38456457d-27, 2.55309556d-27, 1.88904224d-27, 1.38057546d-27, &
156  1.00718330d-27, 7.43581116d-28, 5.63562931d-28, 4.43359435d-28, &
157  3.63923535d-28, 3.11248143d-28, 2.75586846d-28, 2.50672237d-28, &
158  2.32419348d-28, 2.18325682d-28, 2.06834486d-28, 1.93497044d-28, &
159  1.84540751d-28, 1.76356504d-28, 1.68741425d-28, 1.61566157d-28, &
160  1.54754523d-28, 1.48249410d-28, 1.42020176d-28, 1.36045230d-28, &
161  1.30307965d-28 /
162 
163  data f_211 / 4.74439912d-42, 1.95251522d-40, 5.19700194d-39, 8.53120166d-38, &
164  7.72745727d-37, 4.04158559d-36, 1.64853511d-35, 8.56295439d-35, &
165  1.17529722d-33, 2.16867729d-32, 2.40472264d-31, 1.56418133d-30, &
166  7.20032889d-30, 2.65838271d-29, 8.33196904d-29, 2.26128236d-28, &
167  5.24295811d-28, 9.77791121d-28, 1.35913489d-27, 1.43957785d-27, &
168  1.37591544d-27, 1.49029886d-27, 2.06183401d-27, 3.31440622d-27, &
169  5.42497318d-27, 8.41100374d-27, 1.17941366d-26, 1.49269794d-26, &
170  1.71506074d-26, 1.71266353d-26, 1.51434781d-26, 1.36766622d-26, &
171  1.33483562d-26, 1.36834518d-26, 1.45829002d-26, 1.62575306d-26, &
172  1.88773347d-26, 2.22026986d-26, 2.54930499d-26, 2.80758138d-26, &
173  3.06176409d-26, 3.62799792d-26, 5.13226109d-26, 8.46260744d-26, &
174  1.38486586d-25, 1.86192535d-25, 1.78007934d-25, 1.16548409d-25, &
175  5.89293257d-26, 2.69952884d-26, 1.24891081d-26, 6.41273176d-27, &
176  4.08282914d-27, 3.26463328d-27, 2.76230280d-27, 2.08986882d-27, &
177  1.37658470d-27, 8.48489381d-28, 5.19304217d-28, 3.19312514d-28, &
178  2.02968197d-28, 1.50171666d-28, 1.39164218d-28, 1.42448821d-28, &
179  1.41714519d-28, 1.33341059d-28, 1.20759270d-28, 1.07259692d-28, &
180  9.44895400d-29, 8.29030041d-29, 7.25440631d-29, 6.33479483d-29, &
181  5.51563757d-29, 4.79002469d-29, 4.14990482d-29, 3.59384972d-29, &
182  3.12010860d-29, 2.72624742d-29, 2.40734791d-29, 2.15543565d-29, &
183  1.95921688d-29, 1.80682882d-29, 1.68695662d-29, 1.59020936d-29, &
184  1.50940886d-29, 1.43956179d-29, 1.37731622d-29, 1.32049043d-29, &
185  1.26771875d-29, 1.21803879d-29, 1.17074716d-29, 1.10507836d-29, &
186  1.06022834d-29, 1.01703080d-29, 9.75436986d-30, 9.35349257d-30, &
187  8.96744546d-30, 8.59527489d-30, 8.23678940d-30, 7.89144480d-30, &
188  7.55891138d-30 /
189 
190  data f_304 / 3.62695850d-32, 2.79969087d-31, 1.68340584d-30, 7.32681440d-30, &
191  1.99967770d-29, 3.41296785d-29, 4.55409104d-29, 5.94994635d-29, &
192  8.59864963d-29, 1.39787633d-28, 3.17701965d-28, 1.14474920d-27, &
193  4.44845958d-27, 1.54785841d-26, 4.70265345d-26, 1.24524365d-25, &
194  2.81535352d-25, 5.10093666d-25, 6.83545307d-25, 6.82110329d-25, &
195  5.66886188d-25, 4.36205513d-25, 3.29265688d-25, 2.49802368d-25, &
196  1.92527113d-25, 1.51058572d-25, 1.20596047d-25, 9.76884267d-26, &
197  7.89979266d-26, 6.18224289d-26, 4.67298332d-26, 3.57934505d-26, &
198  2.84535785d-26, 2.32853022d-26, 1.95228514d-26, 1.67880071d-26, &
199  1.47608785d-26, 1.32199691d-26, 1.20070960d-26, 1.09378177d-26, &
200  1.00031730d-26, 9.62434001d-27, 1.05063954d-26, 1.27267143d-26, &
201  1.45923057d-26, 1.36746707d-26, 1.03466970d-26, 6.97647829d-27, &
202  4.63141039d-27, 3.19031994d-27, 2.33373613d-27, 1.81589079d-27, &
203  1.48446917d-27, 1.26611478d-27, 1.12617468d-27, 1.03625148d-27, &
204  9.61400595d-28, 8.79016231d-28, 7.82612130d-28, 6.73762960d-28, &
205  5.59717956d-28, 4.53010243d-28, 3.65712196d-28, 3.00958686d-28, &
206  2.54011502d-28, 2.18102277d-28, 1.88736437d-28, 1.63817539d-28, &
207  1.42283147d-28, 1.23631916d-28, 1.07526003d-28, 9.36797928d-29, &
208  8.18565660d-29, 7.18152734d-29, 6.32523238d-29, 5.59513985d-29, &
209  4.96614048d-29, 4.42518826d-29, 3.95487628d-29, 3.54690294d-29, &
210  3.18953930d-29, 2.87720933d-29, 2.60186750d-29, 2.36011522d-29, &
211  2.14717806d-29, 1.95905217d-29, 1.79287981d-29, 1.64562262d-29, &
212  1.51489425d-29, 1.39876064d-29, 1.29496850d-29, 1.18665438d-29, &
213  1.10240474d-29, 1.02643099d-29, 9.57780996d-30, 8.95465151d-30, &
214  8.38950190d-30, 7.87283711d-30, 7.40136507d-30, 6.96804279d-30, &
215  6.56945323d-30 /
216 
217  data f_335 / 2.46882661d-32, 1.89476632d-31, 1.13216502d-30, 4.89532008d-30, &
218  1.32745970d-29, 2.25390335d-29, 3.00511672d-29, 3.96035934d-29, &
219  5.77977656d-29, 8.58600736d-29, 1.14083000d-28, 1.48644411d-28, &
220  2.15788823d-28, 3.51628877d-28, 6.12200698d-28, 1.08184987d-27, &
221  1.85590697d-27, 2.91679107d-27, 3.94405396d-27, 4.63610680d-27, &
222  5.13824456d-27, 5.66602209d-27, 6.30009232d-27, 7.03422868d-27, &
223  7.77973918d-27, 8.32371831d-27, 8.56724316d-27, 8.62601374d-27, &
224  8.13308844d-27, 6.53188216d-27, 4.55197029d-27, 3.57590087d-27, &
225  3.59571707d-27, 4.03502770d-27, 4.54366411d-27, 4.96914990d-27, &
226  5.24601170d-27, 5.39979250d-27, 5.43023669d-27, 5.26235042d-27, &
227  4.91585495d-27, 4.52628362d-27, 4.13385020d-27, 3.67538967d-27, &
228  3.39939742d-27, 3.81284533d-27, 5.02332701d-27, 6.19438602d-27, &
229  6.49613071d-27, 6.04010475d-27, 5.24664275d-27, 4.37225997d-27, &
230  3.52957182d-27, 2.76212276d-27, 2.08473158d-27, 1.50850518d-27, &
231  1.04602472d-27, 7.13091243d-28, 5.34289645d-28, 5.21079581d-28, &
232  6.22246365d-28, 6.99555864d-28, 6.29665489d-28, 4.45077026d-28, &
233  2.67046793d-28, 1.52774686d-28, 9.18061770d-29, 6.09116074d-29, &
234  4.48562572d-29, 3.59463696d-29, 3.05820218d-29, 2.70766652d-29, &
235  2.46144034d-29, 2.27758450d-29, 2.13331183d-29, 2.01537836d-29, &
236  1.91566180d-29, 1.82893912d-29, 1.75167748d-29, 1.68136168d-29, &
237  1.61615595d-29, 1.55481846d-29, 1.49643236d-29, 1.44046656d-29, &
238  1.38657085d-29, 1.33459068d-29, 1.28447380d-29, 1.23615682d-29, &
239  1.18963296d-29, 1.14478976d-29, 1.10146637d-29, 1.04039479d-29, &
240  9.98611410d-30, 9.58205147d-30, 9.19202009d-30, 8.81551313d-30, &
241  8.45252127d-30, 8.10224764d-30, 7.76469090d-30, 7.43954323d-30, &
242  7.12653873d-30 /
243 
244 
245  data n_iris / 41 /
246 
247  data t_iris / 4. , 4.1 , 4.2 , 4.3 , 4.40000001, &
248  4.50000001, 4.60000001, 4.70000001, 4.80000001, 4.90000001, &
249  5.00000001, 5.10000002, 5.20000002, 5.30000002, 5.40000002, &
250  5.50000002, 5.60000002, 5.70000003, 5.80000003, 5.90000003, &
251  6.00000003, 6.10000003, 6.20000003, 6.30000003, 6.40000004, &
252  6.50000004, 6.60000004, 6.70000004, 6.80000004, 6.90000004, &
253  7.00000004, 7.10000005, 7.20000005, 7.30000005, 7.40000005, &
254  7.50000005, 7.60000005, 7.70000006, 7.80000006, 7.90000006, &
255  8.00000006 /
256 
257  data f_1354 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
258  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
259  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
260  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
261  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
262  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.09503647d-39, &
263  5.47214550d-36, 2.42433983d-33, 2.75295034d-31, 1.21929718d-29, &
264  2.48392125d-28, 2.33268145d-27, 8.68623633d-27, 1.00166284d-26, &
265  3.63126633d-27, 7.45174807d-28, 1.38224064d-28, 2.69270994d-29, &
266  5.53314977d-30, 1.15313092d-30, 2.34195788d-31, 4.48242942d-32, &
267  7.94976380d-33 /
268 
269 
270  data n_eis / 60 /
271 
272  data t_eis1 / 1.99526231d+05, 2.23872114d+05, 2.51188643d+05, 2.81838293d+05, &
273  3.16227766d+05, 3.54813389d+05, 3.98107171d+05, 4.46683592d+05, &
274  5.01187234d+05, 5.62341325d+05, 6.30957344d+05, 7.07945784d+05, &
275  7.94328235d+05, 8.91250938d+05, 1.00000000d+06, 1.12201845d+06, &
276  1.25892541d+06, 1.41253754d+06, 1.58489319d+06, 1.77827941d+06, &
277  1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, &
278  3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
279  5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
280  7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
281  1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
282  1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
283  3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
284  5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
285  7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
286  1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08 /
287 
288  data t_eis2 / 1.99526231d+06, 2.23872114d+06, 2.51188643d+06, 2.81838293d+06, &
289  3.16227766d+06, 3.54813389d+06, 3.98107171d+06, 4.46683592d+06, &
290  5.01187234d+06, 5.62341325d+06, 6.30957344d+06, 7.07945784d+06, &
291  7.94328235d+06, 8.91250938d+06, 1.00000000d+07, 1.12201845d+07, &
292  1.25892541d+07, 1.41253754d+07, 1.58489319d+07, 1.77827941d+07, &
293  1.99526231d+07, 2.23872114d+07, 2.51188643d+07, 2.81838293d+07, &
294  3.16227766d+07, 3.54813389d+07, 3.98107171d+07, 4.46683592d+07, &
295  5.01187234d+07, 5.62341325d+07, 6.30957344d+07, 7.07945784d+07, &
296  7.94328235d+07, 8.91250938d+07, 1.00000000d+08, 1.12201845d+08, &
297  1.25892541d+08, 1.41253754d+08, 1.58489319d+08, 1.77827941d+08, &
298  1.99526231d+08, 2.23872114d+08, 2.51188643d+08, 2.81838293d+08, &
299  3.16227766d+08, 3.54813389d+08, 3.98107171d+08, 4.46683592d+08, &
300  5.01187234d+08, 5.62341325d+08, 6.30957344d+08, 7.07945784d+08, &
301  7.94328235d+08, 8.91250938d+08, 1.00000000d+09, 1.12201845d+09, &
302  1.25892541d+09, 1.41253754d+09, 1.58489319d+09, 1.77827941d+09 /
303 
304  data f_263 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
305  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, &
306  0.00000000d+00, 4.46454917d-45, 3.26774829d-42, 1.25292566d-39, &
307  2.66922338d-37, 3.28497742d-35, 2.38677554d-33, 1.03937729d-31, &
308  2.75075687d-30, 4.47961733d-29, 4.46729177d-28, 2.64862689d-27, &
309  8.90863800d-27, 1.72437548d-26, 2.22217752d-26, 2.27999477d-26, &
310  2.08264363d-26, 1.78226687d-26, 1.45821699d-26, 1.14675379d-26, &
311  8.63082492d-27, 6.15925429d-27, 4.11252514d-27, 2.51530564d-27, &
312  1.37090986d-27, 6.42443134d-28, 2.48392636d-28, 7.59187874d-29, &
313  1.77852938d-29, 3.23945221d-30, 4.90533903d-31, 6.75458158d-32, &
314  9.06878868d-33, 1.23927474d-33, 1.75769395d-34, 2.60710914d-35, &
315  4.04318030d-36, 6.53500581d-37, 1.09365022d-37, 1.88383322d-38, &
316  3.31425233d-39, 5.90964084d-40, 1.06147549d-40, 1.90706170d-41, &
317  3.41331584d-42, 6.07310718d-43, 1.07364738d-43, 1.89085498d-44, &
318  3.32598922d-45, 5.87125640d-46, 0.00000000d+00, 0.00000000d+00 /
319 
320  data f_264 / 0.00000000d+00, 2.81670057d-46, 1.28007268d-43, 2.54586603d-41, &
321  2.67887256d-39, 1.68413285d-37, 6.85702304d-36, 1.91797284d-34, &
322  3.84675839d-33, 5.69939170d-32, 6.36224608d-31, 5.39176489d-30, &
323  3.45478458d-29, 1.64848693d-28, 5.71476364d-28, 1.39909997d-27, &
324  2.37743056d-27, 2.86712530d-27, 2.65206348d-27, 2.07175767d-27, &
325  1.47866767d-27, 1.01087374d-27, 6.79605811d-28, 4.54746770d-28, &
326  3.04351751d-28, 2.03639149d-28, 1.35940991d-28, 9.01451939d-29, &
327  5.91289972d-29, 3.81821178d-29, 2.41434696d-29, 1.48871220d-29, &
328  8.93362094d-30, 5.21097445d-30, 2.95964719d-30, 1.64278748d-30, &
329  8.95571660d-31, 4.82096011d-31, 2.57390991d-31, 1.36821781d-31, &
330  7.27136350d-32, 3.87019426d-32, 2.06883430d-32, 1.11228884d-32, &
331  6.01883313d-33, 3.27790676d-33, 1.79805012d-33, 9.93085346d-34, &
332  5.52139556d-34, 3.08881387d-34, 1.73890315d-34, 9.84434964d-35, &
333  5.60603378d-35, 3.20626492d-35, 1.84111068d-35, 0.00000000d+00, &
334  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
335 
336  data f_192 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 4.35772105d-44, &
337  1.26162319d-41, 1.97471205d-39, 1.83409019d-37, 1.08206288d-35, &
338  4.27914363d-34, 1.17943846d-32, 2.32565755d-31, 3.33087991d-30, &
339  3.47013260d-29, 2.60375866d-28, 1.37737127d-27, 5.01053913d-27, &
340  1.23479810d-26, 2.11310542d-26, 2.71831513d-26, 2.89851163d-26, &
341  2.77312376d-26, 2.50025229d-26, 2.18323661d-26, 1.86980322d-26, &
342  1.58035034d-26, 1.31985651d-26, 1.08733133d-26, 8.81804906d-27, &
343  7.00417973d-27, 5.43356567d-27, 4.09857884d-27, 2.99651764d-27, &
344  2.11902962d-27, 1.45014127d-27, 9.62291023d-28, 6.21548647d-28, &
345  3.92807578d-28, 2.44230375d-28, 1.50167782d-28, 9.17611405d-29, &
346  5.58707641d-29, 3.40570915d-29, 2.08030862d-29, 1.27588676d-29, &
347  7.86535588d-30, 4.87646151d-30, 3.03888897d-30, 1.90578649d-30, &
348  1.20195947d-30, 7.61955060d-31, 4.85602199d-31, 3.11049969d-31, &
349  2.00087065d-31, 1.29223740d-31, 8.37422008d-32, 0.00000000d+00, &
350  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
351 
352  data f_255 / 0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 1.76014287d-44, &
353  5.07057938d-42, 7.90473970d-40, 7.31852999d-38, 4.30709255d-36, &
354  1.70009061d-34, 4.67925160d-33, 9.21703546d-32, 1.31918676d-30, &
355  1.37393161d-29, 1.03102379d-28, 5.45694018d-28, 1.98699648d-27, &
356  4.90346776d-27, 8.40524725d-27, 1.08321456d-26, 1.15714525d-26, &
357  1.10905152d-26, 1.00155023d-26, 8.75799161d-27, 7.50935839d-27, &
358  6.35253533d-27, 5.30919268d-27, 4.37669455d-27, 3.55185164d-27, &
359  2.82347055d-27, 2.19257595d-27, 1.65589541d-27, 1.21224987d-27, &
360  8.58395132d-28, 5.88163935d-28, 3.90721447d-28, 2.52593407d-28, &
361  1.59739995d-28, 9.93802874d-29, 6.11343388d-29, 3.73711135d-29, &
362  2.27618743d-29, 1.38793199d-29, 8.48060787d-30, 5.20305940d-30, &
363  3.20867365d-30, 1.99011277d-30, 1.24064551d-30, 7.78310544d-31, &
364  4.91013681d-31, 3.11338381d-31, 1.98451675d-31, 1.27135460d-31, &
365  8.17917486d-32, 5.28280497d-32, 3.42357159d-32, 0.00000000d+00, &
366  0.00000000d+00, 0.00000000d+00, 0.00000000d+00, 0.00000000d+00 /
367 
368  abstract interface
369  subroutine get_subr1(w,x,ixI^L,ixO^L,res)
371  integer, intent(in) :: ixI^L, ixO^L
372  double precision, intent(in) :: w(ixI^S,nw)
373  double precision, intent(in) :: x(ixI^S,1:ndim)
374  double precision, intent(out):: res(ixI^S)
375  end subroutine get_subr1
376 
377  end interface
378 
379  type te_fluid
380 
381  procedure(get_subr1), pointer, nopass :: get_rho => null()
382  procedure(get_subr1), pointer, nopass :: get_pthermal => null()
383  procedure(get_subr1), pointer, nopass :: get_var_rfactor => null()
384 
385  ! factor in eq of state p = Rfactor * rho * T
386  ! used for getting temperature
387  double precision :: rfactor = 1d0
388 
389  end type te_fluid
390 
391 
392  contains
393 
394  subroutine get_line_info(wl,ion,mass,logTe,line_center,spatial_px,spectral_px,sigma_PSF,width_slit)
395  ! get information of the spectral line
396  ! wl: wavelength
397  ! mass: ion mass, unit -- proton mass
398  ! logTe: peak temperature of emission line in logarithm
399  ! line_center: center wavelength of emission line, unit -- Angstrom (0.1 nm)
400  ! spatial_px: pixel size in space of instrument (for image), unit -- arcsec
401  ! spectral_px: pixel size in wagelength of instrument (for spectrum), unit -- Angstrom
402  ! sigma_PSF: width of point spread function core (for instrument), unit -- pixel
403  ! width_slit: width of slit for spectrograph, unit -- arcsec
405 
406  integer, intent(in) :: wl
407  integer, intent(out) :: mass
408  character(len=30), intent(out) :: ion
409  double precision, intent(out) :: logTe,line_center,spatial_px,spectral_px
410  double precision, intent(out) :: sigma_PSF,width_slit
411 
412  select case(wl)
413  case(304)
414  ion='He II'
415  mass=4
416  logte=4.7d0
417  line_center=303.8d0
418  spatial_px=0.6d0
419  spectral_px=0.02d0
420  sigma_psf=0.895d0
421  width_slit=0.6d0
422  case(171)
423  ion='Fe IX'
424  mass=56
425  logte=5.8d0
426  line_center=171.1d0
427  spatial_px=0.6d0
428  spectral_px=0.02d0
429  sigma_psf=1.019d0
430  width_slit=0.6d0
431  case(193)
432  ion='Fe XXIV'
433  mass=56
434  logte=7.3d0
435  line_center=193.5d0
436  spatial_px=0.6d0
437  spectral_px=0.02d0
438  sigma_psf=0.813d0
439  width_slit=0.6d0
440  case(211)
441  ion='Fe XIV'
442  mass=56
443  logte=6.3d0
444  line_center=211.3d0
445  spatial_px=0.6d0
446  spectral_px=0.02d0
447  sigma_psf=0.913d0
448  width_slit=0.6d0
449  case(335)
450  ion='Fe XVI'
451  mass=56
452  logte=6.4d0
453  line_center=335.4d0
454  spatial_px=0.6d0
455  spectral_px=0.02d0
456  sigma_psf=1.019d0
457  width_slit=0.6d0
458  case(94)
459  ion='Fe XVIII'
460  mass=56
461  logte=6.8d0
462  line_center=93.9d0
463  spatial_px=0.6d0
464  spectral_px=0.02d0
465  sigma_psf=1.025d0
466  width_slit=0.6d0
467  case(131)
468  ion='Fe XXI'
469  mass=56
470  logte=7.0d0
471  line_center=131.0d0
472  spatial_px=0.6d0
473  spectral_px=0.02d0
474  sigma_psf=0.984d0
475  width_slit=0.6d0
476  case(1354)
477  ion='Fe XXI'
478  mass=56
479  logte=7.0d0
480  line_center=1354.1d0
481  spatial_px=0.1663d0
482  spectral_px=12.98d-3
483  sigma_psf=1.d0
484  width_slit=0.33d0
485  case(263)
486  ion='Fe XVI'
487  mass=56
488  logte=6.4d0
489  line_center=262.976d0
490  spatial_px=1.d0
491  spectral_px=22.d-3
492  sigma_psf=1.d0
493  width_slit=2.d0
494  case(264)
495  ion='Fe XXIII'
496  mass=56
497  logte=7.1d0
498  line_center=263.765d0
499  spatial_px=1.d0
500  spectral_px=22.d-3
501  sigma_psf=1.d0
502  width_slit=2.d0
503  case(192)
504  ion='Fe XXIV'
505  mass=56
506  logte=7.2d0
507  line_center=192.028d0
508  spatial_px=1.d0
509  spectral_px=22.d-3
510  sigma_psf=1.d0
511  width_slit=2.d0
512  case(255)
513  ion='Fe XXIV'
514  mass=56
515  logte=7.2d0
516  line_center=255.113d0
517  spatial_px=1.d0
518  spectral_px=22.d-3
519  sigma_psf=1.d0
520  width_slit=2.d0
521  case default
522  call mpistop("No information about this line")
523  end select
524  end subroutine get_line_info
525 
526  subroutine get_euv(wl,ixI^L,ixO^L,w,x,fl,flux)
527  ! calculate the local emission intensity of given EUV line (optically thin)
528  ! wavelength is the wave length of the emission line
529  ! unit [DN cm^-1 s^-1 pixel^-1]
530  ! ingrate flux along line of sight: DN s^-1 pixel^-1
532 
533  integer, intent(in) :: wl
534  integer, intent(in) :: ixI^L, ixO^L
535  double precision, intent(in) :: x(ixI^S,1:ndim)
536  double precision, intent(in) :: w(ixI^S,1:nw)
537  type(te_fluid), intent(in) :: fl
538  double precision, intent(out) :: flux(ixI^S)
539 
540  integer :: n_table
541  double precision, allocatable :: t_table(:),f_table(:)
542  integer :: ix^D,iTt,i
543  double precision :: pth(ixI^S),Te(ixI^S),Ne(ixI^S)
544  double precision :: logT,logGT,GT
545 
546  ! selecting emission table
547  select case(wl)
548  case(94)
549  n_table=n_aia
550  allocate(t_table(1:n_table))
551  allocate(f_table(1:n_table))
552  t_table(1:n_table)=t_aia(1:n_aia)
553  f_table(1:n_table)=f_94(1:n_aia)
554  case(131)
555  n_table=n_aia
556  allocate(t_table(1:n_table))
557  allocate(f_table(1:n_table))
558  t_table(1:n_table)=t_aia(1:n_aia)
559  f_table(1:n_table)=f_131(1:n_aia)
560  case(171)
561  n_table=n_aia
562  allocate(t_table(1:n_table))
563  allocate(f_table(1:n_table))
564  t_table(1:n_table)=t_aia(1:n_aia)
565  f_table(1:n_table)=f_171(1:n_aia)
566  case(193)
567  n_table=n_aia
568  allocate(t_table(1:n_table))
569  allocate(f_table(1:n_table))
570  t_table(1:n_table)=t_aia(1:n_aia)
571  f_table(1:n_table)=f_193(1:n_aia)
572  case(211)
573  n_table=n_aia
574  allocate(t_table(1:n_table))
575  allocate(f_table(1:n_table))
576  t_table(1:n_table)=t_aia(1:n_aia)
577  f_table(1:n_table)=f_211(1:n_aia)
578  case(304)
579  n_table=n_aia
580  allocate(t_table(1:n_table))
581  allocate(f_table(1:n_table))
582  t_table(1:n_table)=t_aia(1:n_aia)
583  f_table(1:n_table)=f_304(1:n_aia)
584  case(335)
585  n_table=n_aia
586  allocate(t_table(1:n_table))
587  allocate(f_table(1:n_table))
588  t_table(1:n_table)=t_aia(1:n_aia)
589  f_table(1:n_table)=f_335(1:n_aia)
590  case(1354)
591  n_table=n_iris
592  allocate(t_table(1:n_table))
593  allocate(f_table(1:n_table))
594  t_table(1:n_table)=t_iris(1:n_iris)
595  f_table(1:n_table)=f_1354(1:n_iris)
596  case(263)
597  n_table=n_eis
598  allocate(t_table(1:n_table))
599  allocate(f_table(1:n_table))
600  t_table(1:n_table)=t_eis1(1:n_eis)
601  f_table(1:n_table)=f_263(1:n_eis)
602  case(264)
603  n_table=n_eis
604  allocate(t_table(1:n_table))
605  allocate(f_table(1:n_table))
606  t_table(1:n_table)=t_eis2(1:n_eis)
607  f_table(1:n_table)=f_264(1:n_eis)
608  case(192)
609  n_table=n_eis
610  allocate(t_table(1:n_table))
611  allocate(f_table(1:n_table))
612  t_table(1:n_table)=t_eis2(1:n_eis)
613  f_table(1:n_table)=f_192(1:n_eis)
614  case(255)
615  n_table=n_eis
616  allocate(t_table(1:n_table))
617  allocate(f_table(1:n_table))
618  t_table(1:n_table)=t_eis2(1:n_eis)
619  f_table(1:n_table)=f_255(1:n_eis)
620  case default
621  allocate(t_table(1))
622  allocate(f_table(1))
623  call mpistop("Unknown wavelength")
624  end select
625  call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
626  call fl%get_rho(w,x,ixi^l,ixo^l,ne)
627  if(associated(fl%get_var_Rfactor)) then
628  call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
629  te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
630  else
631  te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*fl%Rfactor)*unit_temperature
632  endif
633  if (si_unit) then
634  ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
635  flux(ixo^s)=ne(ixo^s)**2
636  else
637  ne(ixo^s)=ne(ixo^s)*unit_numberdensity
638  flux(ixo^s)=ne(ixo^s)**2
639  endif
640 
641  select case(wl)
642  case(94,131,171,193,211,304,335,1354)
643  ! temperature table in log
644  do i=1,n_table
645  if(f_table(i) .gt. 1.d-99) then
646  f_table(i)=log10(f_table(i))
647  else
648  f_table(i)=-99.d0
649  endif
650  enddo
651  loggt=zero
652  {do ix^db=ixomin^db,ixomax^db\}
653  logt=log10(te(ix^d))
654  if (logt>=t_table(1) .and. logt<=t_table(n_table)) then
655  do itt=1,n_table-1
656  if (logt>=t_table(itt) .and. logt<t_table(itt+1)) then
657  loggt=f_table(itt)*(logt-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
658  f_table(itt+1)*(logt-t_table(itt))/(t_table(itt+1)-t_table(itt))
659  endif
660  enddo
661  flux(ix^d)=flux(ix^d)*(10**(loggt))
662  if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
663  else
664  flux(ix^d)=zero
665  endif
666  {enddo\}
667  case default
668  ! temperature table linear
669  {do ix^db=ixomin^db,ixomax^db\}
670  if (te(ix^d)>=t_table(1) .and. te(ix^d)<=t_table(n_table)) then
671  do itt=1,n_table-1
672  if (te(ix^d)>=t_table(itt) .and. te(ix^d)<t_table(itt+1)) then
673  gt=f_table(itt)*(te(ix^d)-t_table(itt+1))/(t_table(itt)-t_table(itt+1))+&
674  f_table(itt+1)*(te(ix^d)-t_table(itt))/(t_table(itt+1)-t_table(itt))
675  endif
676  enddo
677  flux(ix^d)=flux(ix^d)*gt
678  if(flux(ix^d)<smalldouble) flux(ix^d)=0.d0
679  else
680  flux(ix^d)=zero
681  endif
682  {enddo\}
683  end select
684 
685  deallocate(t_table,f_table)
686  end subroutine get_euv
687 
688  subroutine get_sxr(ixI^L,ixO^L,w,x,fl,flux,El,Eu)
689  !synthesize thermal SXR from El keV to Eu keV released by cm^-3/m^-3
690  ! volume of plasma in 1 s
691  !flux (cgs): photons cm^-3 s^-1
692  !flux (SI): photons m^-3 s^-1
694 
695  integer, intent(in) :: ixI^L,ixO^L
696  integer, intent(in) :: El,Eu
697  double precision, intent(in) :: x(ixI^S,1:ndim)
698  double precision, intent(in) :: w(ixI^S,nw)
699  type(te_fluid), intent(in) :: fl
700  double precision, intent(out) :: flux(ixI^S)
701 
702  integer :: ix^D,ixO^D
703  integer :: iE,numE
704  double precision :: I0,kb,keV,dE,Ei
705  double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
706  double precision :: Ne(ixI^S),gff(ixI^S),fi(ixI^S)
707  double precision :: EM(ixI^S)
708 
709  i0=3.01d-15 ! I0*4*pi*AU**2, I0 from Pinto (2015)
710  kb=const_kb
711  kev=1.0d3*const_ev
712  de=0.1
713  nume=floor((eu-el)/de)
714  call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
715  call fl%get_rho(w,x,ixi^l,ixo^l,ne)
716 
717  if(associated(fl%get_var_Rfactor)) then
718  call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,te)
719  te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*te(ixo^s))*unit_temperature
720  else
721  te(ixo^s)=pth(ixo^s)/(ne(ixo^s)*fl%Rfactor)*unit_temperature
722  endif
723  if (si_unit) then
724  ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
725  em(ixo^s)=(ne(ixo^s))**2*1.d6 ! cm^-3 m^-3
726  else
727  ne(ixo^s)=ne(ixo^s)*unit_numberdensity
728  em(ixo^s)=(ne(ixo^s))**2
729  endif
730  kbt(ixo^s)=kb*te(ixo^s)/kev
731  flux(ixo^s)=0.0d0
732  do ie=0,nume-1
733  ei=de*ie+el*1.d0
734  gff(ixo^s)=1.d0
735  {do ix^db=ixomin^db,ixomax^db\}
736  if (kbt(ix^d)>0.01*ei) then
737  if(kbt(ix^d)<ei) gff(ix^d)=(kbt(ix^d)/ei)**0.4
738  fi(ix^d)=(em(ix^d)*gff(ix^d))*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
739  else
740  fi(ix^d)=zero
741  endif
742  {enddo\}
743  flux(ixo^s)=flux(ixo^s)+fi(ixo^s)*de
744  enddo
745  flux(ixo^s)=flux(ixo^s)*i0
746  end subroutine get_sxr
747 
748  subroutine get_goes_sxr_flux(xbox^L,fl,eflux)
749  !get GOES SXR 1-8A flux observing at 1AU from given box [w/m^2]
751 
752  double precision, intent(in) :: xbox^L
753  type(te_fluid), intent(in) :: fl
754  double precision, intent(out) :: eflux
755 
756  double precision :: dxb^D,xb^L
757  integer :: iigrid,igrid,j
758  integer :: ixO^L,ixI^L,ix^D
759  double precision :: eflux_grid,eflux_pe
760 
761  ^d&iximin^d=ixglo^d;
762  ^d&iximax^d=ixghi^d;
763  ^d&ixomin^d=ixmlo^d;
764  ^d&ixomax^d=ixmhi^d;
765  eflux_pe=zero
766  do iigrid=1,igridstail; igrid=igrids(iigrid);
767  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
768  ^d&xbmin^d=rnode(rpxmin^d_,igrid);
769  ^d&xbmax^d=rnode(rpxmax^d_,igrid);
770  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)
771  eflux_pe=eflux_pe+eflux_grid
772  enddo
773  call mpi_allreduce(eflux_pe,eflux,1,mpi_double_precision,mpi_sum,icomm,ierrmpi)
774  end subroutine get_goes_sxr_flux
775 
776  subroutine get_goes_flux_grid(ixI^L,ixO^L,w,x,dV,xbox^L,xb^L,fl,eflux_grid)
778 
779  integer, intent(in) :: ixI^L,ixO^L
780  double precision, intent(in) :: x(ixI^S,1:ndim),dV(ixI^S)
781  double precision, intent(in) :: w(ixI^S,nw)
782  double precision, intent(in) :: xbox^L,xb^L
783  type(te_fluid), intent(in) :: fl
784  double precision, intent(out) :: eflux_grid
785 
786  integer :: ix^D,ixO^D,ixb^L
787  integer :: iE,numE,j,inbox
788  double precision :: I0,kb,keV,dE,Ei,El,Eu,A_cgs
789  double precision :: pth(ixI^S),Te(ixI^S),kbT(ixI^S)
790  double precision :: Ne(ixI^S),EM(ixI^S)
791  double precision :: gff,fi,erg_SI
792 
793  ! check whether the grid is inside given box
794  inbox=0
795  {if (xbmin^d<xboxmax^d .and. xbmax^d>xboxmin^d) inbox=inbox+1\}
796 
797  if (inbox==ndim) then
798  ! indexes for cells inside given box
799  ^d&ixbmin^d=ixomin^d;
800  ^d&ixbmax^d=ixomax^d;
801  {if (xbmax^d>xboxmax^d) ixbmax^d=ixomax^d-ceiling((xbmax^d-xboxmax^d)/dxlevel(^d))\}
802  {if (xbmin^d<xboxmin^d) ixbmin^d=ceiling((xboxmin^d-xbmin^d)/dxlevel(^d))+ixomin^d\}
803 
804  i0=1.07d-38 ! photon flux index for observed at 1AU [photon cm^3 m^-2 s^-1 keV^-1]
805  kb=const_kb
806  kev=1.0d3*const_ev
807  erg_si=1.d-7
808  a_cgs=1.d-8 ! Angstrom
809  el=const_h*const_c/(8.d0*a_cgs)/kev ! 8 A
810  eu=const_h*const_c/(1.d0*a_cgs)/kev ! 1 A
811  de=0.1 ! keV
812  nume=floor((eu-el)/de)
813  call fl%get_pthermal(w,x,ixi^l,ixb^l,pth)
814  call fl%get_rho(w,x,ixi^l,ixb^l,ne)
815  if(associated(fl%get_var_Rfactor)) then
816  call fl%get_var_Rfactor(w,x,ixi^l,ixb^l,te)
817  te(ixb^s)=pth(ixb^s)/(ne(ixb^s)*te(ixb^s))*unit_temperature
818  else
819  te(ixb^s)=pth(ixb^s)/(ne(ixb^s)*fl%Rfactor)*unit_temperature
820  endif
821  if (si_unit) then
822  ne(ixo^s)=ne(ixo^s)*unit_numberdensity/1.d6 ! m^-3 -> cm-3
823  em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*(unit_length*1.d2)**3 ! cm^-3
824  else
825  ne(ixo^s)=ne(ixo^s)*unit_numberdensity
826  em(ixb^s)=(i0*(ne(ixb^s))**2)*dv(ixb^s)*unit_length**3
827  endif
828  kbt(ixb^s)=kb*te(ixb^s)/kev
829  eflux_grid=0.0d0
830 
831  do ie=0,nume-1
832  ei=de*ie+el
833  {do ix^db=ixbmin^db,ixbmax^db\}
834  if (kbt(ix^d)>1.d-2*ei) then
835  if(kbt(ix^d)<ei) then
836  gff=(kbt(ix^d)/ei)**0.4
837  else
838  gff=1.d0
839  endif
840  fi=(em(ix^d)*gff)*dexp(-ei/(kbt(ix^d)))/(ei*dsqrt(kbt(ix^d)))
841  eflux_grid=eflux_grid+fi*de*ei
842  endif
843  {enddo\}
844  enddo
845  eflux_grid=eflux_grid*kev*erg_si
846  endif
847 
848  end subroutine get_goes_flux_grid
849 
850  {^ifthreed
851  subroutine get_euv_spectrum(qunit,fl)
853 
854  integer, intent(in) :: qunit
855  type(te_fluid), intent(in) :: fl
856  character(20) :: datatype
857 
858  integer :: mass
859  character (30) :: ion
860  double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
861  double precision :: xslit,arcsec
862 
863  arcsec=7.25d7/unit_length
864 
865  if (mype==0) print *, '###################################################'
866  select case(spectrum_wl)
867  case (1354)
868  if (mype==0) print *, 'Systhesizing EUV spectrum (observed by IRIS).'
869  case (263,264,192,255)
870  if (mype==0) print *, 'Systhesizing EUV spectrum (observed by Hinode/EIS).'
871  case default
872  call mpistop('Wrong wavelength!')
873  end select
874 
875  call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
876  if (mype==0) write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
877  if (mype==0) print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
878  if (mype==0) write(*,'(a,f5.3,a,f5.1,a)') ' Pixel: ',wlrsl,' Angstrom x ',spacersl*725.0, ' km'
879 
881  call mpistop('Wrong spectrum window!')
882  endif
883 
884  datatype='spectrum_euv'
885 
886  if (resolution_spectrum=='data') then
887  if (mype==0) print *, 'Unit of wavelength: Angstrom (0.1 nm) '
888  if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
889  if (si_unit) then
890  if (mype==0) write(*,'(a,f8.1,a)') ' Location of slit: ',location_slit*unit_length/1.d6,' Mm'
891  else
892  if (mype==0) write(*,'(a,f8.1,a)') ' Location of slit: ',location_slit*unit_length/1.d8,' Mm'
893  endif
894  if (mype==0) write(*,'(a,f8.1,a)') ' Width of slit: ',wslit*725.0,' km'
895  call get_spectrum_data_resol(qunit,datatype,fl)
896  else if (resolution_spectrum=='instrument') then
897  if (mype==0) print *, 'Unit of wavelength: Angstrom (0.1 nm) '
898  if (mype==0) print *, 'Unit of length: arcsec (~725 km)'
899  if (mype==0) print *, 'Direction of the slit: parallel to xI2 vector'
900  if (mype==0) write(*,'(a,f8.1,a)') ' Location of slit: xI1 = ',location_slit,' arcsec'
901  if (mype==0) write(*,'(a,f8.1,a)') ' Width of slit: ',wslit,' arcsec'
902  call get_spectrum_inst_resol(qunit,datatype,fl)
903  else
904  call mpistop('Wrong resolution for resolution_spectrum!')
905  endif
906 
907  if (mype==0) print *, '###################################################'
908 
909  end subroutine get_euv_spectrum
910 
911  subroutine get_spectrum_inst_resol(qunit,datatype,fl)
912 
913  integer, intent(in) :: qunit
914  character(20), intent(in) :: datatype
915  type(te_fluid), intent(in) :: fl
916 
917  integer :: numWL,numXS,iwL,ixS,numWI,ix^D
918  double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
919  double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
920  double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
921  double precision :: vec_cor(1:3),xI_cor(1:2)
922  double precision :: res,r_loc,r_max
923 
924  integer :: mass
925  character (30) :: ion
926  double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
927  double precision :: unitv,arcsec,RHESSI_rsl,pixel
928  integer :: iigrid,igrid,i,j,numS
929  double precision :: xLmin,xLmax,xslit
930 
931  call init_vectors()
932 
933  ! calculate domain in space
934  do ix1=1,2
935  if (ix1==1) vec_cor(1)=xprobmin1
936  if (ix1==2) vec_cor(1)=xprobmax1
937  do ix2=1,2
938  if (ix2==1) vec_cor(2)=xprobmin2
939  if (ix2==2) vec_cor(2)=xprobmax2
940  do ix3=1,2
941  if (ix3==1) vec_cor(3)=xprobmin3
942  if (ix3==2) vec_cor(3)=xprobmax3
943  if (big_image) then
944  r_loc=(vec_cor(1)-x_origin(1))**2
945  r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
946  r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
947  r_loc=sqrt(r_loc)
948  if (ix1==1 .and. ix2==1 .and. ix3==1) then
949  r_max=r_loc
950  else
951  r_max=max(r_max,r_loc)
952  endif
953  else
954  call get_cor_image(vec_cor,xi_cor)
955  if (ix1==1 .and. ix2==1 .and. ix3==1) then
956  xsmin=xi_cor(2)
957  xsmax=xi_cor(2)
958  else
959  xsmin=min(xsmin,xi_cor(2))
960  xsmax=max(xsmax,xi_cor(2))
961  endif
962  endif
963  enddo
964  enddo
965  enddo
966  if (big_image) then
967  xsmin=-r_max
968  xsmax=r_max
969  endif
970  xscent=(xsmin+xsmax)/2.d0
971 
972  ! tables for storing spectra data
973  if (si_unit) then
974  arcsec=7.25d5/unit_length
975  else
976  arcsec=7.25d7/unit_length
977  endif
978  call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
979  dxsg=spacersl*arcsec
980  numxs=ceiling((xsmax-xscent)/dxsg)
981  xsmin=xscent-numxs*dxsg
982  xsmax=xscent+numxs*dxsg
983  numxs=numxs*2
984  dwlg=wlrsl
985  numwl=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwlg))
986  wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
987  wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
988  allocate(wl(numwl),dwl(numwl),xs(numxs),dxs(numxs))
989  numwi=1
990  allocate(wi(numwl,numxs,numwi),spectra(numwl,numxs),spectra_rc(numwl,numxs))
991  do iwl=1,numwl
992  wl(iwl)=wlmin+iwl*dwlg-half*dwlg
993  dwl=dwlg
994  enddo
995  do ixs=1,numxs
996  xs(ixs)=xsmin+dxsg*(ixs-half)
997  dxs(ixs)=dxsg
998  enddo
999 
1000  ! find slit and do integration
1001  spectra=zero
1002  do iigrid=1,igridstail; igrid=igrids(iigrid);
1003  do ix1=1,2
1004  if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
1005  if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
1006  do ix2=1,2
1007  if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
1008  if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
1009  do ix3=1,2
1010  if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
1011  if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
1012  call get_cor_image(vec_cor,xi_cor)
1013  if (ix1==1 .and. ix2==1 .and. ix3==1) then
1014  xlmin=xi_cor(1)
1015  xlmax=xi_cor(1)
1016  else
1017  xlmin=min(xlmin,xi_cor(1))
1018  xlmax=max(xlmax,xi_cor(1))
1019  endif
1020  enddo
1021  enddo
1022  enddo
1023 
1024  xslit=location_slit*arcsec
1025  if (xslit>=xlmin-wslit*arcsec .and. xslit<=xlmax+wslit*arcsec) then
1026  call integrate_spectra_inst_resol(igrid,wl,dwlg,xs,dxsg,spectra,numwl,numxs,fl)
1027  endif
1028  enddo
1029 
1030  nums=numwl*numxs
1031  call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1032  mpi_sum,icomm,ierrmpi)
1033  do iwl=1,numwl
1034  do ixs=1,numxs
1035  if (spectra_rc(iwl,ixs)>smalldouble) then
1036  wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1037  else
1038  wi(iwl,ixs,1)=zero
1039  endif
1040  enddo
1041  enddo
1042 
1043  xs=xs/arcsec
1044  dxs=dxs/arcsec
1045 
1046  call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1047 
1048  deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1049 
1050  end subroutine get_spectrum_inst_resol
1051 
1052  subroutine integrate_spectra_inst_resol(igrid,wL,dwLg,xS,dxSg,spectra,numWL,numXS,fl)
1053 
1054  integer, intent(in) :: igrid,numWL,numXS
1055  double precision, intent(in) :: wL(numWL),xS(numXS)
1056  double precision, intent(in) :: dwLg,dxSg
1057  double precision, intent(inout) :: spectra(numWL,numXS)
1058  type(te_fluid), intent(in) :: fl
1059 
1060  integer :: ixO^L,ixI^L,ix^D,ixOnew,j
1061  double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1062  double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
1063  double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),dxIloc(1:2)
1064  integer :: nSubC^D,iSubC^D,iwL,ixS,ixSmin,ixSmax,iwLmin,iwLmax,nwL
1065  double precision :: slit_width,dxSubC^D,xerf^L,fluxSubC
1066  double precision :: xSubC(1:3),xCent(1:2)
1067 
1068  integer :: mass
1069  double precision :: logTe,lineCent
1070  character (30) :: ion
1071  double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1072  double precision :: sigma_wl,sigma_xs,factor
1073 
1074  if (si_unit) then
1075  arcsec=7.25d5/unit_length
1076  else
1077  arcsec=7.25d7/unit_length
1078  endif
1079  xslit=location_slit*arcsec
1080 
1081  call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1082 
1083  ^d&ixomin^d=ixmlo^d\
1084  ^d&ixomax^d=ixmhi^d\
1085  ^d&iximin^d=ixglo^d\
1086  ^d&iximax^d=ixghi^d\
1087  allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1088  ! get local EUV flux and velocity
1089  call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1090  call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1091  call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1092  if(associated(fl%get_var_Rfactor)) then
1093  call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1094  te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1095  else
1096  te(ixo^s)=pth(ixo^s)/(fl%Rfactor*rho(ixo^s))
1097  endif
1098  {do ix^d=ixomin^d,ixomax^d\}
1099  do j=1,3
1100  vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1101  enddo
1102  call dot_product_loc(vloc,vec_los,res)
1103  v(ix^d)=res
1104  {enddo\}
1105 
1106  deallocate(rho)
1107 
1108  slit_width=wslit*arcsec
1109  sigma_wl=sigma_psf*dwlg
1110  sigma_xs=sigma_psf*dxsg
1111  {do ix^d=ixomin^d,ixomax^d\}
1112  xloc(1:3)=ps(igrid)%x(ix^d,1:3)
1113  dxloc(1:3)=ps(igrid)%dx(ix^d,1:3)
1114  call get_cor_image(xloc,xiloc)
1115  call dot_product_loc(dxloc,vec_xi1,res)
1116  dxiloc(1)=abs(res)
1117  if (xiloc(1)>=xslit-half*(slit_width+dxiloc(1)) .and. &
1118  xiloc(1)<=xslit+half*(slit_width+dxiloc(1))) then
1119  ^d&nsubc^d=1;
1120  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(slit_width/16.d0)));
1121  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxsg/4.d0)));
1122  ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1123  ! local line center and line width
1124  if (si_unit) then
1125  fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxsg/dxsg ! DN s^-1
1126  wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1127  else
1128  fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxsg/dxsg ! DN s^-1
1129  wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1130  endif
1131  wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1132  wlwd=wlwd*linecent/const_c
1133  ! dividing a cell to several parts to get more accurate integrating values
1134  {do isubc^d=1,nsubc^d\}
1135  ^d&xsubc(^d)=xloc(^d)-half*dxloc(^d)+(isubc^d-half)*dxsubc^d;
1136  call get_cor_image(xsubc,xcent)
1137  dst_slit=abs(xcent(1)-xslit) ! space distance to slit center
1138  if (dst_slit<=half*slit_width) then
1139  ixs=floor((xcent(2)-(xs(1)-half*dxsg))/dxsg)+1
1140  ixsmin=max(1,ixs-3)
1141  ixsmax=min(ixs+3,numxs)
1142  iwl=floor((wlc-(wl(1)-half*dwlg))/dwlg)+1
1143  nwl=3*ceiling(wlwd/dwlg+1)
1144  iwlmin=max(1,iwl-nwl)
1145  iwlmax=min(iwl+nwl,numwl)
1146  ! calculate the contribution to nearby pixels
1147  do iwl=iwlmin,iwlmax
1148  do ixs=ixsmin,ixsmax
1149  xerfmin1=(wl(iwl)-half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1150  xerfmax1=(wl(iwl)+half*dwlg-wlc)/sqrt(2.d0*(sigma_wl**2+wlwd**2))
1151  xerfmin2=(xs(ixs)-half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1152  xerfmax2=(xs(ixs)+half*dxsg-xcent(2))/(sqrt(2.d0)*sigma_xs)
1153  factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1154  spectra(iwl,ixs)=spectra(iwl,ixs)+fluxsubc*factor
1155  enddo
1156  enddo
1157  ! nearby pixels
1158  endif
1159  {enddo\}
1160  endif
1161  {enddo\}
1162 
1163  deallocate(flux,v,pth,te)
1164  end subroutine integrate_spectra_inst_resol
1165 
1166  subroutine get_spectrum_data_resol(qunit,datatype,fl)
1167 
1168  integer, intent(in) :: qunit
1169  character(20), intent(in) :: datatype
1170  type(te_fluid), intent(in) :: fl
1171 
1172  integer :: numWL,numXS,iwL,ixS,numWI,numS
1173  double precision :: dwLg,xSmin,xSmax,wLmin,wLmax
1174  double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
1175  double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
1176  integer :: strtype,nstrb,nbb,nuni,nstr,bnx
1177  double precision :: qs,dxfirst,dxmid,lenstr
1178 
1179  integer :: iigrid,igrid,j,dir_loc
1180  double precision :: xbmin(1:ndim),xbmax(1:ndim)
1181 
1182  dwlg=1.d-3
1183  numwl=4*int((spectrum_window_max-spectrum_window_min)/(4.d0*dwlg))
1184  wlmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwlg*numwl/2
1185  wlmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwlg*numwl/2
1186  allocate(wl(numwl),dwl(numwl))
1187  dwl(:)=dwlg
1188  do iwl=1,numwl
1189  wl(iwl)=wlmin+iwl*dwlg-half*dwlg
1190  enddo
1191 
1192  select case(direction_slit)
1193  case (1)
1194  numxs=domain_nx1*2**(refine_max_level-1)
1195  xsmin=xprobmin1
1196  xsmax=xprobmax1
1197  bnx=block_nx1
1198  nbb=domain_nx1
1199  strtype=stretch_type(1)
1201  qs=qstretch_baselevel(1)
1202  case (2)
1203  numxs=domain_nx2*2**(refine_max_level-1)
1204  xsmin=xprobmin2
1205  xsmax=xprobmax2
1206  bnx=block_nx2
1207  nbb=domain_nx2
1208  strtype=stretch_type(2)
1210  qs=qstretch_baselevel(2)
1211  case (3)
1212  numxs=domain_nx3*2**(refine_max_level-1)
1213  xsmin=xprobmin3
1214  xsmax=xprobmax3
1215  bnx=block_nx3
1216  nbb=domain_nx3
1217  strtype=stretch_type(3)
1219  qs=qstretch_baselevel(3)
1220  case default
1221  call mpistop('Wrong direction_slit')
1222  end select
1223 
1224  allocate(xs(numxs),dxs(numxs),spectra(numwl,numxs),spectra_rc(numwl,numxs))
1225  numwi=1
1226  allocate(wi(numwl,numxs,numwi))
1227 
1228  select case(strtype)
1229  case(0) ! uniform
1230  dxs(:)=(xsmax-xsmin)/numxs
1231  do ixs=1,numxs
1232  xs(ixs)=xsmin+dxs(ixs)*(ixs-half)
1233  enddo
1234  case(1) ! uni stretch
1235  qs=qs**(one/2**(refine_max_level-1))
1236  dxfirst=(xsmax-xsmin)*(one-qs)/(one-qs**numxs)
1237  dxs(1)=dxfirst
1238  do ixs=2,numxs
1239  dxs(ixs)=dxfirst*qs**(ixs-1)
1240  xs(ixs)=dxs(1)/(one-qs)*(one-qs**(ixs-1))+half*dxs(ixs)
1241  enddo
1242  case(2) ! symm stretch
1243  ! base level, nbb = nstr + nuni + nstr
1244  nstr=nstrb*bnx/2
1245  nuni=nbb-nstrb*bnx
1246  lenstr=(xsmax-xsmin)/(2.d0+nuni*(one-qs)/(one-qs**nstr))
1247  dxfirst=(xsmax-xsmin)/(dble(nuni)+2.d0/(one-qs)*(one-qs**nstr))
1248  dxmid=dxfirst
1249  ! refine_max level, numXI = nstr + nuni + nstr
1250  nstr=nstr*2**(refine_max_level-1)
1251  nuni=nuni*2**(refine_max_level-1)
1252  qs=qs**(one/2**(refine_max_level-1))
1253  dxfirst=lenstr*(one-qs)/(one-qs**nstr)
1254  dxmid=dxmid/2**(refine_max_level-1)
1255  ! uniform center
1256  if(nuni .gt. 0) then
1257  do ixs=nstr+1,nstr+nuni
1258  dxs(ixs)=dxmid
1259  xs(ixs)=lenstr+(dble(ixs)-0.5d0-nstr)*dxs(ixs)+xsmin
1260  enddo
1261  endif
1262  ! left half
1263  do ixs=nstr,1,-1
1264  dxs(ixs)=dxfirst*qs**(nstr-ixs)
1265  xs(ixs)=xsmin+lenstr-dxs(ixs)*half-dxfirst*(one-qs**(nstr-ixs))/(one-qs)
1266  enddo
1267  ! right half
1268  do ixs=nstr+nuni+1,numxs
1269  dxs(ixs)=dxfirst*qs**(ixs-nstr-nuni-1)
1270  xs(ixs)=xsmax-lenstr+dxs(ixs)*half+dxfirst*(one-qs**(ixs-nstr-nuni-1))/(one-qs)
1271  enddo
1272  case default
1273  call mpistop("unknown stretch type")
1274  end select
1275 
1276 
1277  if (los_phi==0 .and. los_theta==90 .and. direction_slit==2) then
1278  ! LOS->x slit->y
1279  dir_loc=3
1280  else if (los_phi==0 .and. los_theta==90 .and. direction_slit==3) then
1281  ! LOS->x slit->z
1282  dir_loc=2
1283  else if (los_phi==90 .and. los_theta==90 .and. direction_slit==1) then
1284  ! LOS->y slit->x
1285  dir_loc=3
1286  else if (los_phi==90 .and. los_theta==90 .and. direction_slit==3) then
1287  ! LOS->y slit->z
1288  dir_loc=1
1289  else if (los_theta==0 .and. direction_slit==1) then
1290  ! LOS->z slit->x
1291  dir_loc=2
1292  else if (los_theta==0 .and. direction_slit==2) then
1293  ! LOS->z slit->y
1294  dir_loc=1
1295  else
1296  call mpistop('Wrong combination of LOS and slit direction!')
1297  endif
1298 
1299  if (dir_loc==1) then
1300  if (location_slit>xprobmax1 .or. location_slit<xprobmin1) then
1301  call mpistop('Wrong value for location_slit!')
1302  endif
1303  else if (dir_loc==2) then
1304  if (location_slit>xprobmax2 .or. location_slit<xprobmin2) then
1305  call mpistop('Wrong value for location_slit!')
1306  endif
1307  else
1308  if (location_slit>xprobmax3 .or. location_slit<xprobmin3) then
1309  call mpistop('Wrong value for location_slit!')
1310  endif
1311  endif
1312 
1313  ! find slit and do integration
1314  spectra=zero
1315  do iigrid=1,igridstail; igrid=igrids(iigrid);
1316  ^d&xbmin(^d)=rnode(rpxmin^d_,igrid);
1317  ^d&xbmax(^d)=rnode(rpxmax^d_,igrid);
1318  if (location_slit>=xbmin(dir_loc) .and. location_slit<xbmax(dir_loc)) then
1319  call integrate_spectra_data_resol(igrid,wl,dwl,spectra,numwl,numxs,dir_loc,fl)
1320  endif
1321  enddo
1322 
1323 
1324  nums=numwl*numxs
1325  call mpi_allreduce(spectra,spectra_rc,nums,mpi_double_precision, &
1326  mpi_sum,icomm,ierrmpi)
1327  do iwl=1,numwl
1328  do ixs=1,numxs
1329  if (spectra_rc(iwl,ixs)>smalldouble) then
1330  wi(iwl,ixs,1)=spectra_rc(iwl,ixs)
1331  else
1332  wi(iwl,ixs,1)=zero
1333  endif
1334  enddo
1335  enddo
1336 
1337  call output_data(qunit,wl,xs,dwl,dxs,wi,numwl,numxs,numwi,datatype)
1338 
1339  deallocate(wl,xs,dwl,dxs,spectra,spectra_rc,wi)
1340 
1341  end subroutine get_spectrum_data_resol
1342 
1343  subroutine integrate_spectra_data_resol(igrid,wL,dwL,spectra,numWL,numXS,dir_loc,fl)
1344  use mod_constants
1345 
1346  integer, intent(in) :: igrid,numWL,numXS,dir_loc
1347  type(te_fluid), intent(in) :: fl
1348  double precision, intent(in) :: wL(numWL),dwL(numWL)
1349  double precision, intent(inout) :: spectra(numWL,numXS)
1350 
1351  integer :: direction_LOS
1352  integer :: ixO^L,ixI^L,ix^D,ixOnew
1353  double precision, allocatable :: flux(:^D&),v(:^D&),pth(:^D&),Te(:^D&),rho(:^D&)
1354  double precision :: wlc,wlwd
1355 
1356  integer :: mass
1357  double precision :: logTe,lineCent
1358  character (30) :: ion
1359  double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
1360 
1361  integer :: levelg,rft,ixSmin,ixSmax,iwL
1362  double precision :: flux_pix,dL
1363 
1364  call get_line_info(spectrum_wl,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1365 
1366  if (los_phi==0 .and. los_theta==90) then
1367  direction_los=1
1368  else if (los_phi==90 .and. los_theta==90) then
1369  direction_los=2
1370  else
1371  direction_los=3
1372  endif
1373 
1374  ^d&ixomin^d=ixmlo^d\
1375  ^d&ixomax^d=ixmhi^d\
1376  ^d&iximin^d=ixglo^d\
1377  ^d&iximax^d=ixghi^d\
1378  allocate(flux(ixi^s),v(ixi^s),pth(ixi^s),te(ixi^s),rho(ixi^s))
1379 
1380  ^d&ix^d=ixomin^d;
1381  if (dir_loc==1) then
1382  do ix1=ixomin1,ixomax1
1383  if (location_slit>=(ps(igrid)%x(ix^d,1)-half*ps(igrid)%dx(ix^d,1)) .and. &
1384  location_slit<(ps(igrid)%x(ix^d,1)+half*ps(igrid)%dx(ix^d,1))) then
1385  ixonew=ix1
1386  endif
1387  enddo
1388  ixomin1=ixonew
1389  ixomax1=ixonew
1390  else if (dir_loc==2) then
1391  do ix2=ixomin2,ixomax2
1392  if (location_slit>=(ps(igrid)%x(ix^d,2)-half*ps(igrid)%dx(ix^d,2)) .and. &
1393  location_slit<(ps(igrid)%x(ix^d,2)+half*ps(igrid)%dx(ix^d,2))) then
1394  ixonew=ix2
1395  endif
1396  enddo
1397  ixomin2=ixonew
1398  ixomax2=ixonew
1399  else
1400  do ix3=ixomin3,ixomax3
1401  if (location_slit>=(ps(igrid)%x(ix^d,3)-half*ps(igrid)%dx(ix^d,3)) .and. &
1402  location_slit<(ps(igrid)%x(ix^d,3)+half*ps(igrid)%dx(ix^d,3))) then
1403  ixonew=ix3
1404  endif
1405  enddo
1406  ixomin3=ixonew
1407  ixomax3=ixonew
1408  endif
1409 
1410  call get_euv(spectrum_wl,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1411  call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1412  v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
1413  call fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,pth)
1414  if(associated(fl%get_var_Rfactor)) then
1415  call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,te)
1416  te(ixo^s)=pth(ixo^s)/(te(ixo^s)*rho(ixo^s))
1417  else
1418  te(ixo^s)=pth(ixo^s)/(fl%Rfactor*rho(ixo^s))
1419  endif
1420 
1421  ! grid parameters
1422  levelg=ps(igrid)%level
1423  rft=2**(refine_max_level-levelg)
1424 
1425  {do ix^d=ixomin^d,ixomax^d\}
1426  if (si_unit) then
1427  wlc=linecent*(1.d0+v(ix^d)*unit_velocity*1.d2/const_c)
1428  else
1429  wlc=linecent*(1.d0+v(ix^d)*unit_velocity/const_c)
1430  endif
1431  wlwd=sqrt(kb_cgs*te(ix^d)*unit_temperature/(mass*mp_cgs))
1432  wlwd=wlwd*linecent/const_c
1433 
1434  select case(direction_slit)
1435  case(1)
1436  ixsmin=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1))*rft+1
1437  ixsmax=(block_nx1*(node(pig1_,igrid)-1)+(ix1-ixomin1+1))*rft
1438  case(2)
1439  ixsmin=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2))*rft+1
1440  ixsmax=(block_nx2*(node(pig2_,igrid)-1)+(ix2-ixomin2+1))*rft
1441  case(3)
1442  ixsmin=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3))*rft+1
1443  ixsmax=(block_nx3*(node(pig3_,igrid)-1)+(ix3-ixomin3+1))*rft
1444  end select
1445 
1446  select case(direction_los)
1447  case(1)
1448  dl=ps(igrid)%dx(ix^d,1)*unit_length
1449  case(2)
1450  dl=ps(igrid)%dx(ix^d,2)*unit_length
1451  case default
1452  dl=ps(igrid)%dx(ix^d,3)*unit_length
1453  end select
1454  if (si_unit) dl=dl*1.d2
1455 
1456  do iwl=1,numwl
1457  flux_pix=flux(ix^d)*wlrsl*dl*exp(-(wl(iwl)-wlc)**2/(2*wlwd**2))/(sqrt(2*dpi)*wlwd)
1458  flux_pix=flux_pix*wslit/spacersl
1459  spectra(iwl,ixsmin:ixsmax)=spectra(iwl,ixsmin:ixsmax)+flux_pix
1460  enddo
1461 
1462  {enddo\}
1463 
1464  deallocate(flux,v,pth,te,rho)
1465 
1466  end subroutine integrate_spectra_data_resol
1467 
1468  subroutine get_euv_image(qunit,fl)
1470 
1471  integer, intent(in) :: qunit
1472  type(te_fluid), intent(in) :: fl
1473  character(20) :: datatype
1474 
1475  integer :: mass
1476  character (30) :: ion
1477  double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1478 
1479  call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1480  if (mype==0) print *, '###################################################'
1481  if (mype==0) print *, 'Systhesizing EUV image'
1482  if (mype==0) write(*,'(a,f8.3,a)') ' Wavelength: ',linecent,' Angstrom'
1483  if (mype==0) print *, 'Unit of EUV flux: DN s^-1 pixel^-1'
1484  if (mype==0) write(*,'(a,f5.1,a,f5.1,a)') ' Pixel: ',spacersl*725.0,' km x ',spacersl*725.0, ' km'
1485  if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
1486  '] of the simulation box is located at [X=0,Y=0] of the image'
1487  if (mype==0) print *, 'Negative Doppler velocity indicates blueshift'
1488 
1489  datatype='image_euv'
1490 
1491  if (resolution_euv=='data') then
1492  if (si_unit) then
1493  if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d6,' Mm'
1494  else
1495  if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d8,' Mm'
1496  endif
1497  if (los_phi==0 .and. los_theta==90) then
1498  call get_image_data_resol(qunit,datatype,fl)
1499  else if (los_phi==90 .and. los_theta==90) then
1500  call get_image_data_resol(qunit,datatype,fl)
1501  else if (los_theta==0) then
1502  call get_image_data_resol(qunit,datatype,fl)
1503  else
1504  call mpistop('ERROR: Wrong LOS for synthesizing emission!')
1505  endif
1506  else if (resolution_euv=='instrument') then
1507  if (mype==0) print *, 'Unit of length: arcsec (~725 km)'
1508  call get_image_inst_resol(qunit,datatype,fl)
1509  else
1510  call mpistop('ERROR: Wrong resolution type')
1511  endif
1512 
1513  if (mype==0) print *, '###################################################'
1514 
1515  end subroutine get_euv_image
1516 
1517  subroutine get_sxr_image(qunit,fl)
1519 
1520  integer, intent(in) :: qunit
1521  type(te_fluid), intent(in) :: fl
1522  character(20) :: datatype
1523 
1524  if (mype==0) print *, '###################################################'
1525  if (mype==0) print *, 'Systhesizing SXR image (observed at 1 AU).'
1526  if (mype==0) write(*,'(a,i2,a,i2,a)') ' Passband: ',emin_sxr,' - ',emax_sxr,' keV'
1527  if (mype==0) print *, 'Unit of SXR flux: photons cm^-2 s^-1 pixel^-1'
1528  if (mype==0) write(*,'(a,f6.1,a,f6.1,a)') ' Pixel: ',2.3*725.0,' km x ',2.3*725.0, ' km'
1529  if (mype==0) write(*,'(a,f6.3,f8.3,f8.3,a)') ' Mapping: [',x_origin(1),x_origin(2),x_origin(3), &
1530  '] of the simulation box is located at [X=0,Y=0] of the image'
1531 
1532  datatype='image_sxr'
1533 
1534  if (resolution_sxr=='data') then
1535  if (si_unit) then
1536  if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d3,' km'
1537  else
1538  if (mype==0) write(*,'(a,f8.1,a)') ' Unit of length: ',unit_length/1.d5,' km'
1539  endif
1540  if (los_phi==0 .and. los_theta==90) then
1541  call get_image_data_resol(qunit,datatype,fl)
1542  else if (los_phi==90 .and. los_theta==90) then
1543  call get_image_data_resol(qunit,datatype,fl)
1544  else if (los_theta==0) then
1545  call get_image_data_resol(qunit,datatype,fl)
1546  else
1547  call mpistop('ERROR: Wrong LOS for synthesizing emission!')
1548  endif
1549  else if (resolution_sxr=='instrument') then
1550  if (mype==0) print *, 'Unit of length: arcsec (~725 km)'
1551  call get_image_inst_resol(qunit,datatype,fl)
1552  else
1553  call mpistop('ERROR: Wrong resolution type')
1554  endif
1555 
1556  if (mype==0) print *, '###################################################'
1557 
1558  end subroutine get_sxr_image
1559 
1560  subroutine get_image_inst_resol(qunit,datatype,fl)
1561  ! integrate emission flux along line of sight (LOS)
1562  ! in a 3D simulation box and get a 2D EUV image
1564  use mod_constants
1565 
1566  integer, intent(in) :: qunit
1567  type(te_fluid), intent(in) :: fl
1568  character(20), intent(in) :: datatype
1569 
1570  integer :: ix^D,numXI1,numXI2,numWI
1571  double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
1572  double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),wI(:,:,:)
1573  double precision, allocatable :: EUVs(:,:),EUV(:,:),Dpls(:,:),Dpl(:,:)
1574  double precision, allocatable :: SXRs(:,:),SXR(:,:)
1575  double precision :: vec_temp1(1:3),vec_temp2(1:3)
1576  double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
1577  double precision :: res,LOS_psi,r_max,r_loc
1578 
1579  integer :: mass
1580  character (30) :: ion
1581  double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
1582  double precision :: unitv,arcsec,RHESSI_rsl,pixel
1583  integer :: iigrid,igrid,i,j,numSI
1584 
1585  call init_vectors()
1586 
1587  ! calculate domain of the image
1588  do ix1=1,2
1589  if (ix1==1) vec_cor(1)=xprobmin1
1590  if (ix1==2) vec_cor(1)=xprobmax1
1591  do ix2=1,2
1592  if (ix2==1) vec_cor(2)=xprobmin2
1593  if (ix2==2) vec_cor(2)=xprobmax2
1594  do ix3=1,2
1595  if (ix3==1) vec_cor(3)=xprobmin3
1596  if (ix3==2) vec_cor(3)=xprobmax3
1597  if (big_image) then
1598  r_loc=(vec_cor(1)-x_origin(1))**2
1599  r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
1600  r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
1601  r_loc=sqrt(r_loc)
1602  if (ix1==1 .and. ix2==1 .and. ix3==1) then
1603  r_max=r_loc
1604  else
1605  r_max=max(r_max,r_loc)
1606  endif
1607  else
1608  call get_cor_image(vec_cor,xi_cor)
1609  if (ix1==1 .and. ix2==1 .and. ix3==1) then
1610  ximin1=xi_cor(1)
1611  ximax1=xi_cor(1)
1612  ximin2=xi_cor(2)
1613  ximax2=xi_cor(2)
1614  else
1615  ximin1=min(ximin1,xi_cor(1))
1616  ximax1=max(ximax1,xi_cor(1))
1617  ximin2=min(ximin2,xi_cor(2))
1618  ximax2=max(ximax2,xi_cor(2))
1619  endif
1620  endif
1621  enddo
1622  enddo
1623  enddo
1624  if (big_image) then
1625  ximin1=-r_max
1626  ximin2=-r_max
1627  ximax1=r_max
1628  ximax2=r_max
1629  endif
1630  xicent1=(ximin1+ximax1)/2.d0
1631  xicent2=(ximin2+ximax2)/2.d0
1632 
1633  ! tables for image
1634  if (si_unit) then
1635  arcsec=7.25d5/unit_length
1636  else
1637  arcsec=7.25d7/unit_length
1638  endif
1639  if (datatype=='image_euv') then
1640  call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1641  dxi=spacersl*arcsec ! intrument resolution of image
1642  else
1643  rhessi_rsl=2.3d0
1644  dxi=rhessi_rsl*arcsec
1645  endif
1646  numxi1=2*ceiling((ximax1-xicent1)/dxi/2.d0)
1647  ximin1=xicent1-numxi1*dxi
1648  ximax1=xicent1+numxi1*dxi
1649  numxi1=numxi1*2
1650  numxi2=2*ceiling((ximax2-xicent2)/dxi/2.d0)
1651  ximin2=xicent2-numxi2*dxi
1652  ximax2=xicent2+numxi2*dxi
1653  numxi2=numxi2*2
1654  allocate(xi1(numxi1),xi2(numxi2),dxi1(numxi1),dxi2(numxi2))
1655  do ix1=1,numxi1
1656  xi1(ix1)=ximin1+dxi*(ix1-half)
1657  dxi1(ix1)=dxi
1658  enddo
1659  do ix2=1,numxi2
1660  xi2(ix2)=ximin2+dxi*(ix2-half)
1661  dxi2(ix2)=dxi
1662  enddo
1663 
1664  ! calculate emission
1665  if (datatype=='image_euv') then
1666  if (si_unit) then
1667  unitv=unit_velocity/1.0e3 ! km/s
1668  else
1669  unitv=unit_velocity/1.0e5 ! km/s
1670  endif
1671  numwi=2
1672  allocate(wi(numxi1,numxi2,numwi))
1673  allocate(euv(numxi1,numxi2),euvs(numxi1,numxi2))
1674  allocate(dpl(numxi1,numxi2),dpls(numxi1,numxi2))
1675  wi=zero
1676  euv=zero
1677  euvs=zero
1678  dpl=zero
1679  dpls=zero
1680  do iigrid=1,igridstail; igrid=igrids(iigrid);
1681  call integrate_euv_inst_resol(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,euvs,dpls)
1682  enddo
1683  numsi=numxi1*numxi2
1684  call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
1685  mpi_sum,icomm,ierrmpi)
1686  call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
1687  mpi_sum,icomm,ierrmpi)
1688  do ix1=1,numxi1
1689  do ix2=1,numxi2
1690  if (euv(ix1,ix2)<smalldouble) euv(ix1,ix2)=zero
1691  if(euv(ix1,ix2)/=0) then
1692  dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
1693  else
1694  dpl(ix1,ix2)=0.d0
1695  endif
1696  if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
1697  wi(ix1,ix2,1)=euv(ix1,ix2)
1698  wi(ix1,ix2,2)=dpl(ix1,ix2)
1699  enddo
1700  enddo
1701 
1702  xi1=xi1/arcsec
1703  dxi1=dxi1/arcsec
1704  xi2=xi2/arcsec
1705  dxi2=dxi2/arcsec
1706  call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
1707  deallocate(wi,euv,euvs,dpl,dpls)
1708  endif
1709 
1710  if (datatype=='image_sxr') then
1711  numwi=1
1712  allocate(wi(numxi1,numxi2,numwi))
1713  allocate(sxr(numxi1,numxi2),sxrs(numxi1,numxi2))
1714  wi=zero
1715  sxr=zero
1716  sxrs=zero
1717  do iigrid=1,igridstail; igrid=igrids(iigrid);
1718  call integrate_sxr_inst_resol(igrid,numxi1,numxi2,xi1,xi2,dxi,fl,sxrs)
1719  enddo
1720  numsi=numxi1*numxi2
1721  call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
1722  mpi_sum,icomm,ierrmpi)
1723 
1724  do ix1=1,numxi1
1725  do ix2=1,numxi2
1726  if (sxr(ix1,ix2)<smalldouble) sxr(ix1,ix2)=zero
1727  enddo
1728  enddo
1729  wi(:,:,1)=sxr(:,:)
1730 
1731  xi1=xi1/arcsec
1732  dxi1=dxi1/arcsec
1733  xi2=xi2/arcsec
1734  dxi2=dxi2/arcsec
1735  call output_data(qunit,xi1,xi2,dxi1,dxi2,wi,numxi1,numxi2,numwi,datatype)
1736  deallocate(wi,sxr,sxrs)
1737  endif
1738 
1739  deallocate(xi1,xi2,dxi1,dxi2)
1740 
1741  end subroutine get_image_inst_resol
1742 
1743  subroutine integrate_sxr_inst_resol(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,SXR)
1744  integer, intent(in) :: igrid,numXI1,numXI2
1745  double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
1746  double precision, intent(in) :: dxI
1747  type(te_fluid), intent(in) :: fl
1748  double precision, intent(out) :: SXR(numXI1,numXI2)
1749 
1750  integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
1751  double precision :: xb^L,xd^D
1752  double precision, allocatable :: flux(:^D&)
1753  double precision :: vloc(1:3),res
1754  integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
1755  double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
1756  double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
1757 
1758  double precision :: sigma_PSF,RHESSI_rsl,sigma0,factor
1759  double precision :: arcsec,pixel,area_1AU
1760 
1761  ^d&ixomin^d=ixmlo^d\
1762  ^d&ixomax^d=ixmhi^d\
1763  ^d&iximin^d=ixglo^d\
1764  ^d&iximax^d=ixghi^d\
1765  ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
1766  ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
1767 
1768  allocate(flux(ixi^s))
1769  ! get local SXR flux
1770  call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
1771 
1772  ! integrate emission
1773  if (si_unit) then
1774  arcsec=7.25d5/unit_length
1775  else
1776  arcsec=7.25d7/unit_length
1777  endif
1778  rhessi_rsl=2.3d0
1779  sigma_psf=1.d0
1780  pixel=rhessi_rsl*arcsec
1781  sigma0=sigma_psf*pixel
1782  area_1au=2.81d27
1783  {do ix^d=ixomin^d,ixomax^d\}
1784  ^d&nsubc^d=1;
1785  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(dxi/4.d0)));
1786  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxi/4.d0)));
1787  ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1788  ! dividing a cell to several parts to get more accurate integrating values
1789  {do isubc^d=1,nsubc^d\}
1790  ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
1791  fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length**3/area_1au
1792  ! mapping the 3D coordinate to location at the image
1793  call get_cor_image(xsubc,xcent)
1794  ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
1795  ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
1796  ixpmin1=max(1,ixp1-3)
1797  ixpmax1=min(ixp1+3,numxi1)
1798  ixpmin2=max(1,ixp2-3)
1799  ixpmax2=min(ixp2+3,numxi2)
1800  do ixp1=ixpmin1,ixpmax1
1801  do ixp2=ixpmin2,ixpmax2
1802  xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1803  xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1804  xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1805  xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1806  factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1807  sxr(ixp1,ixp2)=sxr(ixp1,ixp2)+fluxsubc*factor
1808  enddo !ixP2
1809  enddo !ixP1
1810  {enddo\} !iSubC
1811  {enddo\} !ix
1812 
1813  deallocate(flux)
1814  end subroutine integrate_sxr_inst_resol
1815 
1816  subroutine integrate_euv_inst_resol(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,EUV,Dpl)
1817  integer, intent(in) :: igrid,numXI1,numXI2
1818  double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
1819  double precision, intent(in) :: dxI
1820  type(te_fluid), intent(in) :: fl
1821  double precision, intent(out) :: EUV(numXI1,numXI2),Dpl(numXI1,numXI2)
1822 
1823  integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
1824  double precision :: xb^L,xd^D
1825  double precision, allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
1826  double precision :: vloc(1:3),res
1827  integer :: ixP^L,ixP^D,nSubC^D,iSubC^D
1828  double precision :: xSubP1,xSubP2,dxSubP,xerf^L,fluxsubC
1829  double precision :: xSubC(1:3),dxSubC^D,xCent(1:2)
1830 
1831  integer :: mass
1832  double precision :: logTe
1833  character (30) :: ion
1834  double precision :: lineCent
1835  double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
1836  double precision :: unitv,arcsec,pixel
1837 
1838  ^d&ixomin^d=ixmlo^d\
1839  ^d&ixomax^d=ixmhi^d\
1840  ^d&iximin^d=ixglo^d\
1841  ^d&iximax^d=ixghi^d\
1842  ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
1843  ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
1844 
1845  allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
1846  ! get local EUV flux and velocity
1847  call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
1848  call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
1849  {do ix^d=ixomin^d,ixomax^d\}
1850  do j=1,3
1851  vloc(j)=ps(igrid)%w(ix^d,iw_mom(j))/rho(ix^d)
1852  enddo
1853  call dot_product_loc(vloc,vec_los,res)
1854  v(ix^d)=res
1855  {enddo\}
1856  deallocate(rho)
1857 
1858  ! integrate emission
1859  call get_line_info(wavelength,ion,mass,logte,linecent,spacersl,wlrsl,sigma_psf,wslit)
1860  if (si_unit) then
1861  arcsec=7.25d5/unit_length
1862  else
1863  arcsec=7.25d7/unit_length
1864  endif
1865  pixel=spacersl*arcsec
1866  sigma0=sigma_psf*pixel
1867  {do ix^d=ixomin^d,ixomax^d\}
1868  ^d&nsubc^d=1;
1869  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi1(^d))/(dxi/2.d0)));
1870  ^d&nsubc^d=max(nsubc^d,ceiling(ps(igrid)%dx(ix^dd,^d)*abs(vec_xi2(^d))/(dxi/2.d0)));
1871  ^d&dxsubc^d=ps(igrid)%dx(ix^dd,^d)/nsubc^d;
1872  if (si_unit) then
1873  fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length*1.d2/dxi/dxi ! DN s^-1
1874  else
1875  fluxsubc=flux(ix^d)*dxsubc1*dxsubc2*dxsubc3*unit_length/dxi/dxi ! DN s^-1
1876  endif
1877  ! dividing a cell to several parts to get more accurate integrating values
1878  {do isubc^d=1,nsubc^d\}
1879  ^d&xsubc(^d)=ps(igrid)%x(ix^dd,^d)-half*ps(igrid)%dx(ix^dd,^d)+(isubc^d-half)*dxsubc^d;
1880  ! mapping the 3D coordinate to location at the image
1881  call get_cor_image(xsubc,xcent)
1882  ! distribution at nearby pixels
1883  ixp1=floor((xcent(1)-(xi1(1)-half*dxi))/dxi)+1
1884  ixp2=floor((xcent(2)-(xi2(1)-half*dxi))/dxi)+1
1885  ixpmin1=max(1,ixp1-3)
1886  ixpmax1=min(ixp1+3,numxi1)
1887  ixpmin2=max(1,ixp2-3)
1888  ixpmax2=min(ixp2+3,numxi2)
1889  do ixp1=ixpmin1,ixpmax1
1890  do ixp2=ixpmin2,ixpmax2
1891  xerfmin1=((xi1(ixp1)-half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1892  xerfmax1=((xi1(ixp1)+half*dxi)-xcent(1))/(sqrt(2.d0)*sigma0)
1893  xerfmin2=((xi2(ixp2)-half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1894  xerfmax2=((xi2(ixp2)+half*dxi)-xcent(2))/(sqrt(2.d0)*sigma0)
1895  factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-erfc(xerfmax2))/4.d0
1896  euv(ixp1,ixp2)=euv(ixp1,ixp2)+fluxsubc*factor
1897  dpl(ixp1,ixp2)=dpl(ixp1,ixp2)+fluxsubc*factor*v(ix^d)
1898  enddo !ixP2
1899  enddo !ixP1
1900  {enddo\} !iSubC
1901  {enddo\} !ix
1902 
1903  deallocate(flux,v)
1904  end subroutine integrate_euv_inst_resol
1905 
1906  subroutine get_image_data_resol(qunit,datatype,fl)
1907  ! integrate emission flux along line of sight (LOS)
1908  ! in a 3D simulation box and get a 2D EUV image
1910  use mod_constants
1911 
1912  integer, intent(in) :: qunit
1913  character(20), intent(in) :: datatype
1914  type(te_fluid), intent(in) :: fl
1915 
1916  double precision :: dx^D
1917  integer :: numX^D,ix^D
1918  double precision, allocatable :: EUV(:,:),EUVs(:,:),Dpl(:,:),Dpls(:,:)
1919  double precision, allocatable :: SXR(:,:),SXRs(:,:),wI(:,:,:)
1920  double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:),dxIi
1921  integer :: numXI1,numXI2,numSI,numWI
1922  double precision :: xI^L
1923  integer :: iigrid,igrid,i,j
1924  double precision, allocatable :: xIF1(:),xIF2(:),dxIF1(:),dxIF2(:)
1925  integer :: nXIF1,nXIF2
1926  double precision :: xIF^L
1927 
1928  double precision :: unitv,arcsec,RHESSI_rsl
1929  integer :: strtype^D,nstrb^D,nbb^D,nuni^D,nstr^D,bnx^D
1930  double precision :: qs^D,dxfirst^D,dxmid^D,lenstr^D
1931 
1932  numx1=domain_nx1*2**(refine_max_level-1)
1933  numx2=domain_nx2*2**(refine_max_level-1)
1934  numx3=domain_nx3*2**(refine_max_level-1)
1935 
1936  ! parameters for creating table
1937  if (los_phi==0 .and. los_theta==90) then
1938  nxif1=domain_nx2*2**(refine_max_level-1)
1939  nxif2=domain_nx3*2**(refine_max_level-1)
1940  xifmin1=xprobmin2
1941  xifmax1=xprobmax2
1942  xifmin2=xprobmin3
1943  xifmax2=xprobmax3
1944  bnx1=block_nx2
1945  bnx2=block_nx3
1946  nbb1=domain_nx2
1947  nbb2=domain_nx3
1948  strtype1=stretch_type(2)
1949  strtype2=stretch_type(3)
1950  nstrb1=nstretchedblocks_baselevel(2)
1951  nstrb2=nstretchedblocks_baselevel(3)
1952  qs1=qstretch_baselevel(2)
1953  qs2=qstretch_baselevel(3)
1954  if (mype==0) write(*,'(a)') ' LOS vector: [-1.00 0.00 0.00]'
1955  if (mype==0) write(*,'(a)') ' xI1 vector: [ 0.00 1.00 0.00]'
1956  if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
1957  else if (los_phi==90 .and. los_theta==90) then
1958  nxif1=domain_nx3*2**(refine_max_level-1)
1959  nxif2=domain_nx1*2**(refine_max_level-1)
1960  xifmin1=xprobmin3
1961  xifmax1=xprobmax3
1962  xifmin2=xprobmin1
1963  xifmax2=xprobmax1
1964  bnx1=block_nx3
1965  bnx2=block_nx1
1966  nbb1=domain_nx3
1967  nbb2=domain_nx1
1968  strtype1=stretch_type(3)
1969  strtype2=stretch_type(1)
1970  nstrb1=nstretchedblocks_baselevel(3)
1971  nstrb2=nstretchedblocks_baselevel(1)
1972  qs1=qstretch_baselevel(3)
1973  qs2=qstretch_baselevel(1)
1974  if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 -1.00 0.00]'
1975  if (mype==0) write(*,'(a)') ' xI1 vector: [-1.00 0.00 0.00]'
1976  if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 0.00 1.00]'
1977  else
1978  nxif1=domain_nx1*2**(refine_max_level-1)
1979  nxif2=domain_nx2*2**(refine_max_level-1)
1980  xifmin1=xprobmin1
1981  xifmax1=xprobmax1
1982  xifmin2=xprobmin2
1983  xifmax2=xprobmax2
1984  bnx1=block_nx1
1985  bnx2=block_nx2
1986  nbb1=domain_nx1
1987  nbb2=domain_nx2
1988  strtype1=stretch_type(1)
1989  strtype2=stretch_type(2)
1990  nstrb1=nstretchedblocks_baselevel(1)
1991  nstrb2=nstretchedblocks_baselevel(2)
1992  qs1=qstretch_baselevel(1)
1993  qs2=qstretch_baselevel(2)
1994  if (mype==0) write(*,'(a)') ' LOS vector: [ 0.00 0.00 -1.00]'
1995  if (mype==0) write(*,'(a)') ' xI1 vector: [ 1.00 0.00 0.00]'
1996  if (mype==0) write(*,'(a)') ' xI2 vector: [ 0.00 1.00 0.00]'
1997  endif
1998  allocate(xif1(nxif1),xif2(nxif2),dxif1(nxif1),dxif2(nxif2))
1999 
2000  ! initialize image coordinate
2001  select case(strtype1)
2002  case(0) ! uniform
2003  dxif1(:)=(xifmax1-xifmin1)/nxif1
2004  do ix1=1,nxif1
2005  xif1(ix1)=xifmin1+dxif1(ix1)*(ix1-half)
2006  enddo
2007  case(1) ! uni stretch
2008  qs1=qs1**(one/2**(refine_max_level-1))
2009  dxfirst1=(xifmax1-xifmin1)*(one-qs1)/(one-qs1**nxif1)
2010  dxif1(1)=dxfirst1
2011  do ix1=2,nxif1
2012  dxif1(ix1)=dxfirst1*qs1**(ix1-1)
2013  xif1(ix1)=dxif1(1)/(one-qs1)*(one-qs1**(ix1-1))+half*dxif1(ix1)
2014  enddo
2015  case(2) ! symm stretch
2016  ! base level, nbb = nstr + nuni + nstr
2017  nstr1=nstrb1*bnx1/2
2018  nuni1=nbb1-nstrb1*bnx1
2019  lenstr1=(xifmax1-xifmin1)/(2.d0+nuni1*(one-qs1)/(one-qs1**nstr1))
2020  dxfirst1=(xifmax1-xifmin1)/(dble(nuni1)+2.d0/(one-qs1)*(one-qs1**nstr1))
2021  dxmid1=dxfirst1
2022  ! refine_max level, numXI = nstr + nuni + nstr
2023  nstr1=nstr1*2**(refine_max_level-1)
2024  nuni1=nuni1*2**(refine_max_level-1)
2025  qs1=qs1**(one/2**(refine_max_level-1))
2026  dxfirst1=lenstr1*(one-qs1)/(one-qs1**nstr1)
2027  dxmid1=dxmid1/2**(refine_max_level-1)
2028  ! uniform center
2029  if(nuni1 .gt. 0) then
2030  do ix1=nstr1+1,nstr1+nuni1
2031  dxif1(ix1)=dxmid1
2032  xif1(ix1)=lenstr1+(dble(ix1)-0.5d0-nstr1)*dxif1(ix1)+xifmin1
2033  enddo
2034  endif
2035  ! left half
2036  do ix1=nstr1,1,-1
2037  dxif1(ix1)=dxfirst1*qs1**(nstr1-ix1)
2038  xif1(ix1)=xifmin1+lenstr1-dxif1(ix1)*half-dxfirst1*(one-qs1**(nstr1-ix1))/(one-qs1)
2039  enddo
2040  ! right half
2041  do ix1=nstr1+nuni1+1,nxif1
2042  dxif1(ix1)=dxfirst1*qs1**(ix1-nstr1-nuni1-1)
2043  xif1(ix1)=xifmax1-lenstr1+dxif1(ix1)*half+dxfirst1*(one-qs1**(ix1-nstr1-nuni1-1))/(one-qs1)
2044  enddo
2045  case default
2046  call mpistop("unknown stretch type")
2047  end select
2048 
2049  select case(strtype2)
2050  case(0) ! uniform
2051  dxif2(:)=(xifmax2-xifmin2)/nxif2
2052  do ix2=1,nxif2
2053  xif2(ix2)=xifmin2+dxif2(ix2)*(ix2-half)
2054  enddo
2055  case(1) ! uni stretch
2056  qs2=qs2**(one/2**(refine_max_level-1))
2057  dxfirst2=(xifmax2-xifmin2)*(one-qs2)/(one-qs2**nxif2)
2058  dxif2(1)=dxfirst2
2059  do ix2=2,nxif1
2060  dxif2(ix2)=dxfirst2*qs2**(ix2-1)
2061  xif2(ix2)=dxif2(1)/(one-qs2)*(one-qs2**(ix2-1))+half*dxif2(ix2)
2062  enddo
2063  case(2) ! symm stretch
2064  ! base level, nbb = nstr + nuni + nstr
2065  nstr2=nstrb2*bnx2/2
2066  nuni2=nbb2-nstrb2*bnx2
2067  lenstr2=(xifmax2-xifmin2)/(2.d0+nuni2*(one-qs2)/(one-qs2**nstr2))
2068  dxfirst2=(xifmax2-xifmin2)/(dble(nuni2)+2.d0/(one-qs2)*(one-qs2**nstr2))
2069  dxmid2=dxfirst2
2070  ! refine_max level, numXI = nstr + nuni + nstr
2071  nstr2=nstr2*2**(refine_max_level-1)
2072  nuni2=nuni2*2**(refine_max_level-1)
2073  qs2=qs2**(one/2**(refine_max_level-1))
2074  dxfirst2=lenstr2*(one-qs2)/(one-qs2**nstr2)
2075  dxmid2=dxmid2/2**(refine_max_level-1)
2076  ! uniform center
2077  if(nuni2 .gt. 0) then
2078  do ix2=nstr2+1,nstr2+nuni2
2079  dxif2(ix2)=dxmid2
2080  xif2(ix2)=lenstr2+(dble(ix2)-0.5d0-nstr2)*dxif2(ix2)+xifmin2
2081  enddo
2082  endif
2083  ! left half
2084  do ix2=nstr2,1,-1
2085  dxif2(ix2)=dxfirst2*qs2**(nstr2-ix2)
2086  xif2(ix2)=xifmin2+lenstr2-dxif2(ix2)*half-dxfirst2*(one-qs2**(nstr2-ix2))/(one-qs2)
2087  enddo
2088  ! right half
2089  do ix2=nstr2+nuni2+1,nxif2
2090  dxif2(ix2)=dxfirst2*qs2**(ix2-nstr2-nuni2-1)
2091  xif2(ix2)=xifmax2-lenstr2+dxif2(ix2)*half+dxfirst2*(one-qs2**(ix2-nstr2-nuni2-1))/(one-qs2)
2092  enddo
2093  case default
2094  call mpistop("unknown stretch type")
2095  end select
2096 
2097  ! integrate EUV flux and get cell average flux for image
2098  if (datatype=='image_euv') then
2099  if (si_unit) then
2100  unitv=unit_velocity/1.0e3 ! km/s
2101  else
2102  unitv=unit_velocity/1.0e5 ! km/s
2103  endif
2104  numwi=2
2105  allocate(wi(nxif1,nxif2,numwi))
2106  allocate(euvs(nxif1,nxif2),euv(nxif1,nxif2))
2107  allocate(dpl(nxif1,nxif2),dpls(nxif1,nxif2))
2108  euvs=0.0d0
2109  euv=0.0d0
2110  dpl=0.d0
2111  dpls=0.d0
2112  do iigrid=1,igridstail; igrid=igrids(iigrid);
2113  call integrate_euv_data_resol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,euvs,dpls)
2114  enddo
2115  numsi=nxif1*nxif2
2116  call mpi_allreduce(euvs,euv,numsi,mpi_double_precision, &
2117  mpi_sum,icomm,ierrmpi)
2118  call mpi_allreduce(dpls,dpl,numsi,mpi_double_precision, &
2119  mpi_sum,icomm,ierrmpi)
2120  do ix1=1,nxif1
2121  do ix2=1,nxif2
2122  if (euv(ix1,ix2)<smalldouble) euv(ix1,ix2)=zero
2123  if(euv(ix1,ix2)/=0) then
2124  dpl(ix1,ix2)=(dpl(ix1,ix2)/euv(ix1,ix2))*unitv
2125  else
2126  dpl(ix1,ix2)=0.d0
2127  endif
2128  if (abs(dpl(ix1,ix2))<smalldouble) dpl(ix1,ix2)=zero
2129  enddo
2130  enddo
2131  wi(:,:,1)=euv(:,:)
2132  wi(:,:,2)=dpl(:,:)
2133 
2134  call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2135  deallocate(wi,euv,euvs,dpl,dpls)
2136  endif
2137 
2138  ! integrate EUV flux and get cell average flux for image
2139  if (datatype=='image_sxr') then
2140  if (si_unit) then
2141  arcsec=7.25d5
2142  else
2143  arcsec=7.25d7
2144  endif
2145  rhessi_rsl=2.3d0
2146  numwi=1
2147  allocate(wi(nxif1,nxif2,numwi))
2148  allocate(sxrs(nxif1,nxif2),sxr(nxif1,nxif2))
2149  sxrs=0.0d0
2150  sxr=0.0d0
2151  do iigrid=1,igridstail; igrid=igrids(iigrid);
2152  call integrate_sxr_data_resol(igrid,nxif1,nxif2,xif1,xif2,dxif1,dxif2,fl,sxrs)
2153  enddo
2154  numsi=nxif1*nxif2
2155  call mpi_allreduce(sxrs,sxr,numsi,mpi_double_precision, &
2156  mpi_sum,icomm,ierrmpi)
2157 
2158  sxr=sxr*(rhessi_rsl*arcsec)**2
2159  do ix1=1,nxif1
2160  do ix2=1,nxif2
2161  if (sxr(ix1,ix2)<smalldouble) sxr(ix1,ix2)=zero
2162  enddo
2163  enddo
2164  wi(:,:,1)=sxr(:,:)
2165 
2166  call output_data(qunit,xif1,xif2,dxif1,dxif2,wi,nxif1,nxif2,numwi,datatype)
2167  deallocate(wi,sxr,sxrs)
2168  endif
2169 
2170  deallocate(xif1,xif2,dxif1,dxif2)
2171 
2172  end subroutine get_image_data_resol
2173 
2174  subroutine integrate_euv_data_resol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,EUV,Dpl)
2176 
2177  integer, intent(in) :: igrid,nXIF1,nXIF2
2178  double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2179  double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2180  type(te_fluid), intent(in) :: fl
2181  double precision, intent(out) :: EUV(nXIF1,nXIF2),Dpl(nXIF1,nXIF2)
2182 
2183  integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2184  double precision :: xb^L,xd^D
2185  double precision, allocatable :: flux(:^D&),v(:^D&),rho(:^D&)
2186  double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2187  double precision, allocatable :: EUVg(:,:),Fvg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2188  integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2189  double precision :: EUVt,Fvt,xc^L,xg^L,r2
2190  integer :: ixP^L,ixP^D
2191  integer :: direction_LOS
2192 
2193  if (los_phi==0 .and. los_theta==90) then
2194  direction_los=1
2195  else if (los_phi==90 .and. los_theta==90) then
2196  direction_los=2
2197  else
2198  direction_los=3
2199  endif
2200 
2201  ^d&ixomin^d=ixmlo^d\
2202  ^d&ixomax^d=ixmhi^d\
2203  ^d&iximin^d=ixglo^d\
2204  ^d&iximax^d=ixghi^d\
2205  ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
2206  ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
2207 
2208  allocate(flux(ixi^s),v(ixi^s),rho(ixi^s))
2209  allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2210  dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2211  dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2212  dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2213  ! get local EUV flux and velocity
2214  call get_euv(wavelength,ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux)
2215  call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixi^l,ixo^l,rho)
2216  v(ixo^s)=-ps(igrid)%w(ixo^s,iw_mom(direction_los))/rho(ixo^s)
2217  deallocate(rho)
2218 
2219  ! grid parameters
2220  levelg=ps(igrid)%level
2221  rft=2**(refine_max_level-levelg)
2222 
2223  ! fine table for storing EUV flux of current grid
2224  select case(direction_los)
2225  case(1)
2226  nxg1=iximax2*rft
2227  nxg2=iximax3*rft
2228  case(2)
2229  nxg1=iximax3*rft
2230  nxg2=iximax1*rft
2231  case(3)
2232  nxg1=iximax1*rft
2233  nxg2=iximax2*rft
2234  end select
2235  allocate(euvg(nxg1,nxg2),fvg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2236  euvg=zero
2237  fvg=zero
2238  xg1=zero
2239  xg2=zero
2240 
2241  ! integrate for different direction
2242  select case(direction_los)
2243  case(1)
2244  do ix2=ixomin2,ixomax2
2245  ixgmin1=(ix2-1)*rft+1
2246  ixgmax1=ix2*rft
2247  do ix3=ixomin3,ixomax3
2248  ixgmin2=(ix3-1)*rft+1
2249  ixgmax2=ix3*rft
2250  euvt=0.d0
2251  fvt=0.d0
2252  do ix1=ixomin1,ixomax1
2253  euvt=euvt+flux(ix^d)*dxb1(ix^d)*unit_length
2254  fvt=fvt+flux(ix^d)*dxb1(ix^d)*unit_length*v(ix^d)
2255  enddo
2256  euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2257  fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2258  enddo
2259  enddo
2260  case(2)
2261  do ix3=ixomin3,ixomax3
2262  ixgmin1=(ix3-1)*rft+1
2263  ixgmax1=ix3*rft
2264  do ix1=ixomin1,ixomax1
2265  ixgmin2=(ix1-1)*rft+1
2266  ixgmax2=ix1*rft
2267  euvt=0.d0
2268  fvt=0.d0
2269  do ix2=ixomin2,ixomax2
2270  euvt=euvt+flux(ix^d)*dxb2(ix^d)*unit_length
2271  fvt=fvt+flux(ix^d)*dxb2(ix^d)*unit_length*v(ix^d)
2272  enddo
2273  euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2274  fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2275  enddo
2276  enddo
2277  case(3)
2278  do ix1=ixomin1,ixomax1
2279  ixgmin1=(ix1-1)*rft+1
2280  ixgmax1=ix1*rft
2281  do ix2=ixomin2,ixomax2
2282  ixgmin2=(ix2-1)*rft+1
2283  ixgmax2=ix2*rft
2284  euvt=0.d0
2285  fvt=0.d0
2286  do ix3=ixomin3,ixomax3
2287  euvt=euvt+flux(ix^d)*dxb3(ix^d)*unit_length
2288  fvt=fvt+flux(ix^d)*dxb3(ix^d)*unit_length*v(ix^d)
2289  enddo
2290  euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=euvt
2291  fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=fvt
2292  enddo
2293  enddo
2294  end select
2295  if (si_unit) then
2296  euvg=euvg*1.d2
2297  fvg=fvg*1.d2
2298  endif
2299 
2300  ! mapping grid data to global table
2301  ! index ranges in local table
2302  select case(direction_los)
2303  case(1)
2304  ixgmin1=(ixomin2-1)*rft+1
2305  ixgmax1=ixomax2*rft
2306  ixgmin2=(ixomin3-1)*rft+1
2307  ixgmax2=ixomax3*rft
2308  case(2)
2309  ixgmin1=(ixomin3-1)*rft+1
2310  ixgmax1=ixomax3*rft
2311  ixgmin2=(ixomin1-1)*rft+1
2312  ixgmax2=ixomax1*rft
2313  case(3)
2314  ixgmin1=(ixomin1-1)*rft+1
2315  ixgmax1=ixomax1*rft
2316  ixgmin2=(ixomin2-1)*rft+1
2317  ixgmax2=ixomax2*rft
2318  end select
2319  ! index ranges in global table & mapping
2320  select case(direction_los)
2321  case(1)
2322  ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
2323  ixpmax1=node(pig2_,igrid)*rft*block_nx2
2324  ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
2325  ixpmax2=node(pig3_,igrid)*rft*block_nx3
2326  case(2)
2327  ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
2328  ixpmax1=node(pig3_,igrid)*rft*block_nx3
2329  ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
2330  ixpmax2=node(pig1_,igrid)*rft*block_nx1
2331  case(3)
2332  ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
2333  ixpmax1=node(pig1_,igrid)*rft*block_nx1
2334  ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
2335  ixpmax2=node(pig2_,igrid)*rft*block_nx2
2336  end select
2337  xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2338  xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2339  dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2340  dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2341  euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=euv(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2342  euvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2343  dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=dpl(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2344  fvg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2345 
2346  deallocate(flux,v,dxb1,dxb2,dxb3,euvg,fvg,xg1,xg2,dxg1,dxg2)
2347 
2348  end subroutine integrate_euv_data_resol
2349 
2350  subroutine integrate_sxr_data_resol(igrid,nXIF1,nXIF2,xIF1,xIF2,dxIF1,dxIF2,fl,SXR)
2352 
2353  integer, intent(in) :: igrid,nXIF1,nXIF2
2354  double precision, intent(in) :: xIF1(nXIF1),xIF2(nXIF2)
2355  double precision, intent(in) :: dxIF1(nXIF1),dxIF2(nXIF2)
2356  type(te_fluid), intent(in) :: fl
2357  double precision, intent(out) :: SXR(nXIF1,nXIF2)
2358 
2359  integer :: ixO^L,ixO^D,ixI^L,ix^D,i,j
2360  double precision :: xb^L,xd^D
2361  double precision, allocatable :: flux(:^D&)
2362  double precision, allocatable :: dxb1(:^D&),dxb2(:^D&),dxb3(:^D&)
2363  double precision, allocatable :: SXRg(:,:),xg1(:),xg2(:),dxg1(:),dxg2(:)
2364  integer :: levelg,nXg1,nXg2,iXgmin1,iXgmax1,iXgmin2,iXgmax2,rft,iXg^D
2365  double precision :: SXRt,xc^L,xg^L,r2,area_1AU
2366  integer :: ixP^L,ixP^D
2367  integer :: direction_LOS
2368 
2369  if (los_phi==0 .and. los_theta==90) then
2370  direction_los=1
2371  else if (los_phi==90 .and. los_theta==90) then
2372  direction_los=2
2373  else
2374  direction_los=3
2375  endif
2376 
2377  ^d&ixomin^d=ixmlo^d\
2378  ^d&ixomax^d=ixmhi^d\
2379  ^d&iximin^d=ixglo^d\
2380  ^d&iximax^d=ixghi^d\
2381  ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
2382  ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
2383 
2384  allocate(flux(ixi^s))
2385  allocate(dxb1(ixi^s),dxb2(ixi^s),dxb3(ixi^s))
2386  dxb1(ixo^s)=ps(igrid)%dx(ixo^s,1)
2387  dxb2(ixo^s)=ps(igrid)%dx(ixo^s,2)
2388  dxb3(ixo^s)=ps(igrid)%dx(ixo^s,3)
2389  ! get local SXR flux
2390  call get_sxr(ixi^l,ixo^l,ps(igrid)%w,ps(igrid)%x,fl,flux,emin_sxr,emax_sxr)
2391 
2392  ! grid parameters
2393  levelg=ps(igrid)%level
2394  rft=2**(refine_max_level-levelg)
2395 
2396  ! fine table for storing EUV flux of current grid
2397  select case(direction_los)
2398  case(1)
2399  nxg1=iximax2*rft
2400  nxg2=iximax3*rft
2401  case(2)
2402  nxg1=iximax3*rft
2403  nxg2=iximax1*rft
2404  case(3)
2405  nxg1=iximax1*rft
2406  nxg2=iximax2*rft
2407  end select
2408  allocate(sxrg(nxg1,nxg2),xg1(nxg1),xg2(nxg2),dxg1(nxg1),dxg2(nxg2))
2409  sxrg=zero
2410  xg1=zero
2411  xg2=zero
2412 
2413  ! integrate for different direction
2414  select case(direction_los)
2415  case(1)
2416  do ix2=ixomin2,ixomax2
2417  ixgmin1=(ix2-1)*rft+1
2418  ixgmax1=ix2*rft
2419  do ix3=ixomin3,ixomax3
2420  ixgmin2=(ix3-1)*rft+1
2421  ixgmax2=ix3*rft
2422  sxrt=0.d0
2423  do ix1=ixomin1,ixomax1
2424  sxrt=sxrt+flux(ix^d)*dxb1(ix^d)*unit_length
2425  enddo
2426  sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2427  enddo
2428  enddo
2429  case(2)
2430  do ix3=ixomin3,ixomax3
2431  ixgmin1=(ix3-1)*rft+1
2432  ixgmax1=ix3*rft
2433  do ix1=ixomin1,ixomax1
2434  ixgmin2=(ix1-1)*rft+1
2435  ixgmax2=ix1*rft
2436  sxrt=0.d0
2437  do ix2=ixomin2,ixomax2
2438  sxrt=sxrt+flux(ix^d)*dxb2(ix^d)*unit_length
2439  enddo
2440  sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2441  enddo
2442  enddo
2443  case(3)
2444  do ix1=ixomin1,ixomax1
2445  ixgmin1=(ix1-1)*rft+1
2446  ixgmax1=ix1*rft
2447  do ix2=ixomin2,ixomax2
2448  ixgmin2=(ix2-1)*rft+1
2449  ixgmax2=ix2*rft
2450  sxrt=0.d0
2451  do ix3=ixomin3,ixomax3
2452  sxrt=sxrt+flux(ix^d)*dxb3(ix^d)*unit_length
2453  enddo
2454  sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)=sxrt
2455  enddo
2456  enddo
2457  end select
2458 
2459  area_1au=2.81d27
2460  sxrg=sxrg/area_1au
2461 
2462  ! mapping grid data to global table
2463  ! index ranges in local table
2464  select case(direction_los)
2465  case(1)
2466  ixgmin1=(ixomin2-1)*rft+1
2467  ixgmax1=ixomax2*rft
2468  ixgmin2=(ixomin3-1)*rft+1
2469  ixgmax2=ixomax3*rft
2470  case(2)
2471  ixgmin1=(ixomin3-1)*rft+1
2472  ixgmax1=ixomax3*rft
2473  ixgmin2=(ixomin1-1)*rft+1
2474  ixgmax2=ixomax1*rft
2475  case(3)
2476  ixgmin1=(ixomin1-1)*rft+1
2477  ixgmax1=ixomax1*rft
2478  ixgmin2=(ixomin2-1)*rft+1
2479  ixgmax2=ixomax2*rft
2480  end select
2481  ! index ranges in global table & mapping
2482  select case(direction_los)
2483  case(1)
2484  ixpmin1=(node(pig2_,igrid)-1)*rft*block_nx2+1
2485  ixpmax1=node(pig2_,igrid)*rft*block_nx2
2486  ixpmin2=(node(pig3_,igrid)-1)*rft*block_nx3+1
2487  ixpmax2=node(pig3_,igrid)*rft*block_nx3
2488  case(2)
2489  ixpmin1=(node(pig3_,igrid)-1)*rft*block_nx3+1
2490  ixpmax1=node(pig3_,igrid)*rft*block_nx3
2491  ixpmin2=(node(pig1_,igrid)-1)*rft*block_nx1+1
2492  ixpmax2=node(pig1_,igrid)*rft*block_nx1
2493  case(3)
2494  ixpmin1=(node(pig1_,igrid)-1)*rft*block_nx1+1
2495  ixpmax1=node(pig1_,igrid)*rft*block_nx1
2496  ixpmin2=(node(pig2_,igrid)-1)*rft*block_nx2+1
2497  ixpmax2=node(pig2_,igrid)*rft*block_nx2
2498  end select
2499  xg1(ixgmin1:ixgmax1)=xif1(ixpmin1:ixpmax1)
2500  xg2(ixgmin2:ixgmax2)=xif2(ixpmin2:ixpmax2)
2501  dxg1(ixgmin1:ixgmax1)=dxif1(ixpmin1:ixpmax1)
2502  dxg2(ixgmin2:ixgmax2)=dxif2(ixpmin2:ixpmax2)
2503  sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)=sxr(ixpmin1:ixpmax1,ixpmin2:ixpmax2)+&
2504  sxrg(ixgmin1:ixgmax1,ixgmin2:ixgmax2)
2505 
2506  deallocate(flux,dxb1,dxb2,dxb3,sxrg,xg1,xg2,dxg1,dxg2)
2507 
2508  end subroutine integrate_sxr_data_resol
2509 
2510  subroutine output_data(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,datatype)
2511  ! change the format of data and write data
2513 
2514  integer, intent(in) :: qunit,nXO1,nXO2,nWO
2515  double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
2516  double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
2517  double precision, intent(inout) :: wO(nXO1,nXO2,nWO)
2518  character(20), intent(in) :: datatype
2519 
2520  integer :: nPiece,nP1,nP2,nC1,nC2,nWC
2521  integer :: piece_nmax1,piece_nmax2,ix1,ix2,j,ipc,ixc1,ixc2
2522  double precision, allocatable :: xC(:,:,:,:),wC(:,:,:,:),dxC(:,:,:,:)
2523  character(len=std_len) :: resolution_type
2524 
2525  ! clean small values
2526  do ix1=1,nxo1
2527  do ix2=1,nxo2
2528  do j=1,nwo
2529  if (abs(wo(ix1,ix2,j))<smalldouble) wo(ix1,ix2,j)=zero
2530  enddo
2531  enddo
2532  enddo
2533 
2534  ! how many cells in each grid
2535  if (datatype=='image_euv' .and. resolution_euv=='data') then
2536  if (los_phi==0 .and. los_theta==90) then
2537  piece_nmax1=block_nx2
2538  piece_nmax2=block_nx3
2539  else if (los_phi==90 .and. los_theta==90) then
2540  piece_nmax1=block_nx3
2541  piece_nmax2=block_nx1
2542  else
2543  piece_nmax1=block_nx1
2544  piece_nmax2=block_nx2
2545  endif
2546  else if (datatype=='image_sxr' .and. resolution_sxr=='data') then
2547  if (los_phi==0 .and. los_theta==90) then
2548  piece_nmax1=block_nx2
2549  piece_nmax2=block_nx3
2550  else if (los_phi==90 .and. los_theta==90) then
2551  piece_nmax1=block_nx3
2552  piece_nmax2=block_nx1
2553  else
2554  piece_nmax1=block_nx1
2555  piece_nmax2=block_nx2
2556  endif
2557  else if (datatype=='spectrum_euv' .and. resolution_spectrum=='data') then
2558  piece_nmax1=16
2559  if (direction_slit==1) then
2560  piece_nmax2=block_nx1
2561  else if (direction_slit==2) then
2562  piece_nmax2=block_nx2
2563  else
2564  piece_nmax2=block_nx3
2565  endif
2566  else
2567  piece_nmax1=20
2568  piece_nmax2=20
2569  endif
2570 
2571  loopn1: do j=piece_nmax1,1,-1
2572  if(mod(nxo1,j)==0) then
2573  nc1=j
2574  exit loopn1
2575  endif
2576  enddo loopn1
2577  loopn2: do j=piece_nmax2,1,-1
2578  if(mod(nxo2,j)==0) then
2579  nc2=j
2580  exit loopn2
2581  endif
2582  enddo loopn2
2583  ! how many grids
2584  np1=nxo1/nc1
2585  np2=nxo2/nc2
2586  npiece=np1*np2
2587  nwc=nwo
2588 
2589  select case(convert_type)
2590  case('EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi')
2591  ! put data into grids
2592  allocate(xc(npiece,nc1,nc2,2))
2593  allocate(dxc(npiece,nc1,nc2,2))
2594  allocate(wc(npiece,nc1,nc2,nwo))
2595  do ipc=1,npiece
2596  do ixc1=1,nc1
2597  do ixc2=1,nc2
2598  ix1=mod(ipc-1,np1)*nc1+ixc1
2599  ix2=floor(1.0*(ipc-1)/np1)*nc2+ixc2
2600  xc(ipc,ixc1,ixc2,1)=xo1(ix1)
2601  xc(ipc,ixc1,ixc2,2)=xo2(ix2)
2602  dxc(ipc,ixc1,ixc2,1)=dxo1(ix1)
2603  dxc(ipc,ixc1,ixc2,2)=dxo2(ix2)
2604  do j=1,nwc
2605  wc(ipc,ixc1,ixc2,j)=wo(ix1,ix2,j)
2606  enddo
2607  enddo
2608  enddo
2609  enddo
2610  ! write data into vtu file
2611  call write_image_vtucc(qunit,xc,wc,dxc,npiece,nc1,nc2,nwc,datatype)
2612  deallocate(xc,dxc,wc)
2613  case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi')
2614  if (convert_type=='EIvtiCCmpi') resolution_type=resolution_euv
2615  if (convert_type=='ESvtiCCmpi') resolution_type=resolution_spectrum
2616  if (convert_type=='SIvtiCCmpi') resolution_type=resolution_sxr
2617  if (sum(stretch_type(:))>0 .and. resolution_type=='data') then
2618  call mpistop("Error in synthesize emission: vti is not supported for data resolution")
2619  else
2620  call write_image_vticc(qunit,xo1,xo2,dxo1,dxo2,wo,nxo1,nxo2,nwo,nc1,nc2)
2621  endif
2622  case default
2623  call mpistop("Error in synthesize emission: Unknown convert_type")
2624  end select
2625 
2626  end subroutine output_data
2627  }
2628 
2629  subroutine write_image_vticc(qunit,xO1,xO2,dxO1,dxO2,wO,nXO1,nXO2,nWO,nC1,nC2)
2630  ! write image data to vti
2632 
2633  integer, intent(in) :: qunit,nXO1,nXO2,nWO,nC1,nC2
2634  double precision, intent(in) :: xO1(nXO1),xO2(nxO2)
2635  double precision, intent(in) :: dxO1(nxO1),dxO2(nxO2)
2636  double precision, intent(in) :: wO(nXO1,nXO2,nWO)
2637 
2638  double precision :: origin(1:3), spacing(1:3)
2639  integer :: wholeExtent(1:6),extent(1:6)
2640  integer :: nP1,nP2,iP1,iP2,iw
2641  integer :: ixC1,ixC2,ixCmin1,ixCmax1,ixCmin2,ixCmax2
2642 
2643  integer :: filenr
2644  logical :: fileopen
2645  character (70) :: subname,wname,vname,nameL,nameS
2646  character (len=std_len) :: filename
2647  integer :: mass
2648  double precision :: logTe
2649  character (30) :: ion
2650  double precision :: line_center
2651  double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
2652 
2653 
2654  origin(1)=xo1(1)-0.5d0*dxo1(1)
2655  origin(2)=xo2(1)-0.5d0*dxo2(1)
2656  origin(3)=zero
2657  spacing(1)=dxo1(1)
2658  spacing(2)=dxo2(1)
2659  spacing(3)=zero
2660  wholeextent=zero
2661  wholeextent(2)=nxo1
2662  wholeextent(4)=nxo2
2663  np1=nxo1/nc1
2664  np2=nxo2/nc2
2665 
2666  ! get information of emission line
2667  if (convert_type=='EIvtiCCmpi') then
2668  call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
2669  else if (convert_type=='ESvtiCCmpi') then
2670  call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
2671  endif
2672 
2673  if (mype==0) then
2674  inquire(qunit,opened=fileopen)
2675  if(.not.fileopen)then
2676  ! generate filename
2677  filenr=snapshotini
2678  if (autoconvert) filenr=snapshotnext
2679  if (convert_type=='EIvtiCCmpi') then
2680  write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vti"
2681  else if (convert_type=='SIvtiCCmpi') then
2682  write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vti"
2683  else if (convert_type=='ESvtiCCmpi') then
2684  write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vti"
2685  endif
2686  open(qunit,file=filename,status='unknown',form='formatted')
2687  endif
2688 
2689  ! generate xml header
2690  write(qunit,'(a)')'<?xml version="1.0"?>'
2691  write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
2692  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2693  write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
2694  origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
2695  ! file info
2696  write(qunit,'(a)')'<FieldData>'
2697  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2698  'NumberOfTuples="1" format="ascii">'
2699  write(qunit,*) real(global_time*time_convert_factor)
2700  write(qunit,'(a)')'</DataArray>'
2701  if (convert_type=='EIvtiCCmpi' .or. convert_type=='ESvtiCCmpi') then
2702  write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
2703  'NumberOfTuples="1" format="ascii">'
2704  write(qunit,*) real(logte)
2705  write(qunit,'(a)')'</DataArray>'
2706  endif
2707  write(qunit,'(a)')'</FieldData>'
2708  ! pixel/cell data
2709  do ip1=1,np1
2710  do ip2=1,np2
2711  extent=zero
2712  extent(1)=(ip1-1)*nc1
2713  extent(2)=ip1*nc1
2714  extent(3)=(ip2-1)*nc2
2715  extent(4)=ip2*nc2
2716  ixcmin1=extent(1)+1
2717  ixcmax1=extent(2)
2718  ixcmin2=extent(3)+1
2719  ixcmax2=extent(4)
2720  write(qunit,'(a,6(i10),a)') &
2721  '<Piece Extent="',extent,'">'
2722  write(qunit,'(a)')'<CellData>'
2723  do iw=1,nwo
2724  ! variable name
2725  if (convert_type=='EIvtiCCmpi') then
2726  if (iw==1) then
2727  if (wavelength<100) then
2728  write(vname,'(a,i2)') "AIA",wavelength
2729  else if (wavelength<1000) then
2730  write(vname,'(a,i3)') "AIA",wavelength
2731  else
2732  write(vname,'(a,i4)') "IRIS",wavelength
2733  endif
2734  endif
2735  if (iw==2) vname='Doppler_velocity'
2736  else if (convert_type=='SIvtiCCmpi') then
2737  if (emin_sxr<10 .and. emax_sxr<10) then
2738  write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2739  else if (emin_sxr<10 .and. emax_sxr>=10) then
2740  write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2741  else
2742  write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2743  endif
2744  else if (convert_type=='ESvtiCCmpi') then
2745  if (spectrum_wl==1354) then
2746  write(vname,'(a,i4)') "SG",spectrum_wl
2747  else
2748  write(vname,'(a,i3)') "EIS",spectrum_wl
2749  endif
2750  endif
2751  write(qunit,'(a,a,a)')&
2752  '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
2753  write(qunit,'(200(1pe14.6))') ((wo(ixc1,ixc2,iw),ixc1=ixcmin1,ixcmax1),ixc2=ixcmin2,ixcmax2)
2754  write(qunit,'(a)')'</DataArray>'
2755  enddo
2756  write(qunit,'(a)')'</CellData>'
2757  write(qunit,'(a)')'</Piece>'
2758  enddo
2759  enddo
2760  ! end
2761  write(qunit,'(a)')'</ImageData>'
2762  write(qunit,'(a)')'</VTKFile>'
2763  close(qunit)
2764  endif
2765 
2766  end subroutine write_image_vticc
2767 
2768  subroutine write_image_vtucc(qunit,xC,wC,dxC,nPiece,nC1,nC2,nWC,datatype)
2769  ! write image data to vtu
2771 
2772  integer, intent(in) :: qunit
2773  integer, intent(in) :: nPiece,nC1,nC2,nWC
2774  double precision, intent(in) :: xC(nPiece,nC1,nC2,2),dxC(nPiece,nc1,nc2,2)
2775  double precision, intent(in) :: wC(nPiece,nC1,nC2,nWC)
2776  character(20), intent(in) :: datatype
2777 
2778  integer :: nP1,nP2
2779  double precision :: xP(nPiece,nC1+1,nC2+1,2)
2780  integer :: filenr
2781  logical :: fileopen
2782  character (70) :: subname,wname,vname,nameL,nameS
2783  character (len=std_len) :: filename
2784  integer :: ixC1,ixC2,ixP,ix1,ix2,j
2785  integer :: nc,np,icel,VTK_type
2786  integer :: mass
2787  double precision :: logTe
2788  character (30) :: ion
2789  double precision :: line_center
2790  double precision :: spatial_rsl,spectral_rsl,sigma_PSF,wslit
2791 
2792  np1=nc1+1
2793  np2=nc2+1
2794  np=np1*np2
2795  nc=nc1*nc2
2796  ! cell corner location
2797  do ixp=1,npiece
2798  do ix1=1,np1
2799  do ix2=1,np2
2800  if (ix1<np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1,1,1)-0.5d0*dxc(ixp,ix1,1,1)
2801  if (ix1==np1) xp(ixp,ix1,ix2,1)=xc(ixp,ix1-1,1,1)+0.5d0*dxc(ixp,ix1-1,1,1)
2802  if (ix2<np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2,2)-0.5d0*dxc(ixp,1,ix2,2)
2803  if (ix2==np2) xp(ixp,ix1,ix2,2)=xc(ixp,1,ix2-1,2)+0.5d0*dxc(ixp,1,ix2-1,2)
2804  enddo
2805  enddo
2806  enddo
2807  ! get information of emission line
2808  if (datatype=='image_euv') then
2809  call get_line_info(wavelength,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
2810  else if (datatype=='spectrum_euv') then
2811  call get_line_info(spectrum_wl,ion,mass,logte,line_center,spatial_rsl,spectral_rsl,sigma_psf,wslit)
2812  endif
2813 
2814  if (mype==0) then
2815  inquire(qunit,opened=fileopen)
2816  if(.not.fileopen)then
2817  ! generate filename
2818  filenr=snapshotini
2819  if (autoconvert) filenr=snapshotnext
2820  if (datatype=='image_euv') then
2821  write(filename,'(a,i4.4,a)') trim(filename_euv),filenr,".vtu"
2822  else if (datatype=='image_sxr') then
2823  write(filename,'(a,i4.4,a)') trim(filename_sxr),filenr,".vtu"
2824  else if (datatype=='spectrum_euv') then
2825  write(filename,'(a,i4.4,a)') trim(filename_spectrum),filenr,".vtu"
2826  endif
2827  open(qunit,file=filename,status='unknown',form='formatted')
2828  endif
2829  ! generate xml header
2830  write(qunit,'(a)')'<?xml version="1.0"?>'
2831  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2832  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2833  write(qunit,'(a)')'<UnstructuredGrid>'
2834  write(qunit,'(a)')'<FieldData>'
2835  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2836  'NumberOfTuples="1" format="ascii">'
2837  write(qunit,*) real(global_time*time_convert_factor)
2838  write(qunit,'(a)')'</DataArray>'
2839  if (datatype=='image_euv' .or. datatype=='spectrum_euv') then
2840  write(qunit,'(2a)')'<DataArray type="Float32" Name="logT" ',&
2841  'NumberOfTuples="1" format="ascii">'
2842  write(qunit,*) real(logte)
2843  write(qunit,'(a)')'</DataArray>'
2844  endif
2845  write(qunit,'(a)')'</FieldData>'
2846  do ixp=1,npiece
2847  write(qunit,'(a,i7,a,i7,a)') &
2848  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2849  write(qunit,'(a)')'<CellData>'
2850  do j=1,nwc
2851  if (datatype=='image_euv') then
2852  if (j==1) then
2853  if (wavelength<100) then
2854  write(vname,'(a,i2)') "AIA",wavelength
2855  else if (wavelength<1000) then
2856  write(vname,'(a,i3)') "AIA",wavelength
2857  else
2858  write(vname,'(a,i4)') "IRIS",wavelength
2859  endif
2860  endif
2861  if (j==2) vname='Doppler_velocity'
2862  else if (datatype=='image_sxr') then
2863  if (emin_sxr<10 .and. emax_sxr<10) then
2864  write(vname,'(a,i1,a,i1,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2865  else if (emin_sxr<10 .and. emax_sxr>=10) then
2866  write(vname,'(a,i1,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2867  else
2868  write(vname,'(a,i2,a,i2,a)') "SXR",emin_sxr,"-",emax_sxr,"keV"
2869  endif
2870  else if (datatype=='spectrum_euv') then
2871  if (spectrum_wl==1354) then
2872  write(vname,'(a,i4)') "SG",spectrum_wl
2873  else
2874  write(vname,'(a,i3)') "EIS",spectrum_wl
2875  endif
2876  endif
2877  write(qunit,'(a,a,a)')&
2878  '<DataArray type="Float64" Name="',trim(vname),'" format="ascii">'
2879  write(qunit,'(200(1pe14.6))') ((wc(ixp,ixc1,ixc2,j),ixc1=1,nc1),ixc2=1,nc2)
2880  write(qunit,'(a)')'</DataArray>'
2881  enddo
2882  write(qunit,'(a)')'</CellData>'
2883  write(qunit,'(a)')'<Points>'
2884  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2885  do ix2=1,np2
2886  do ix1=1,np1
2887  if (datatype=='image_euv' .and. resolution_euv=='data') then
2888  if (los_phi==0 .and. los_theta==90) then
2889  write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
2890  else if (los_phi==90 .and. los_theta==90) then
2891  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
2892  else
2893  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2894  endif
2895  else if (datatype=='image_sxr' .and. resolution_sxr=='data') then
2896  if (los_phi==0 .and. los_theta==90) then
2897  write(qunit,'(3(1pe14.6))') 0.d0,xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2)
2898  else if (los_phi==90 .and. los_theta==90) then
2899  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,2),0.d0,xp(ixp,ix1,ix2,1)
2900  else
2901  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2902  endif
2903  else
2904  write(qunit,'(3(1pe14.6))') xp(ixp,ix1,ix2,1),xp(ixp,ix1,ix2,2),0.d0
2905  endif
2906  enddo
2907  enddo
2908  write(qunit,'(a)')'</DataArray>'
2909  write(qunit,'(a)')'</Points>'
2910  ! connetivity part
2911  write(qunit,'(a)')'<Cells>'
2912  write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
2913  do ix2=1,nc2
2914  do ix1=1,nc1
2915  write(qunit,'(4(i7))') ix1-1+(ix2-1)*np1,ix1+(ix2-1)*np1,&
2916  ix1-1+ix2*np1,ix1+ix2*np1
2917  enddo
2918  enddo
2919  write(qunit,'(a)')'</DataArray>'
2920  ! offsets data array
2921  write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
2922  do icel=1,nc
2923  write(qunit,'(i7)') icel*(2**2)
2924  enddo
2925  write(qunit,'(a)')'</DataArray>'
2926  ! VTK cell type data array
2927  write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
2928  ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
2929  vtk_type=8
2930  do icel=1,nc
2931  write(qunit,'(i2)') vtk_type
2932  enddo
2933  write(qunit,'(a)')'</DataArray>'
2934  write(qunit,'(a)')'</Cells>'
2935  write(qunit,'(a)')'</Piece>'
2936  enddo
2937  write(qunit,'(a)')'</UnstructuredGrid>'
2938  write(qunit,'(a)')'</VTKFile>'
2939  close(qunit)
2940  endif
2941  end subroutine write_image_vtucc
2942 
2943  subroutine dot_product_loc(vec_in1,vec_in2,res)
2944  double precision, intent(in) :: vec_in1(1:3),vec_in2(1:3)
2945  double precision, intent(out) :: res
2946  integer :: j
2947 
2948  res=zero
2949  do j=1,3
2950  res=res+vec_in1(j)*vec_in2(j)
2951  enddo
2952 
2953  end subroutine dot_product_loc
2954 
2955  subroutine cross_product_loc(vec_in1,vec_in2,vec_out)
2956  double precision, intent(in) :: vec_in1(1:3),vec_in2(1:3)
2957  double precision, intent(out) :: vec_out(1:3)
2958 
2959  vec_out(1)=vec_in1(2)*vec_in2(3)-vec_in1(3)*vec_in2(2)
2960  vec_out(2)=vec_in1(3)*vec_in2(1)-vec_in1(1)*vec_in2(3)
2961  vec_out(3)=vec_in1(1)*vec_in2(2)-vec_in1(2)*vec_in2(1)
2962 
2963  end subroutine cross_product_loc
2964 
2965  subroutine init_vectors()
2966  integer :: j
2967  double precision :: LOS_psi
2968  double precision :: vec_z(1:3),vec_temp1(1:3),vec_temp2(1:3)
2969 
2970  ! vectors for image coordinate
2971  vec_los(1)=-cos(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
2972  vec_los(2)=-sin(dpi*los_phi/180.d0)*sin(dpi*los_theta/180.d0)
2973  vec_los(3)=-cos(dpi*los_theta/180.d0)
2974  do j=1,3
2975  if (abs(vec_los(j))<=smalldouble) vec_los(j)=zero
2976  enddo
2977  vec_z(:)=zero
2978  vec_z(3)=1.d0
2979  if (los_theta==zero) then
2980  vec_xi1=zero
2981  vec_xi2=zero
2982  vec_xi1(1)=1.d0
2983  vec_xi2(2)=1.d0
2984  else
2985  call cross_product_loc(vec_los,vec_z,vec_xi1)
2987  endif
2988  vec_temp1=vec_xi1/sqrt(vec_xi1(1)**2+vec_xi1(2)**2+vec_xi1(3)**2)
2989  vec_temp2=vec_xi2/sqrt(vec_xi2(1)**2+vec_xi2(2)**2+vec_xi2(3)**2)
2990  los_psi=dpi*image_rotate/180.d0
2991  vec_xi1=vec_temp1*cos(los_psi)-vec_temp2*sin(los_psi)
2992  vec_xi2=vec_temp2*cos(los_psi)+vec_temp1*sin(los_psi)
2993 
2994  do j=1,3
2995  if (abs(vec_xi1(j))<smalldouble) vec_xi1(j)=zero
2996  if (abs(vec_xi2(j))<smalldouble) vec_xi2(j)=zero
2997  enddo
2998 
2999  if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' LOS vector: [',vec_los(1),vec_los(2),vec_los(3),']'
3000  if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI1 vector: [',vec_xi1(1),vec_xi1(2),vec_xi1(3),']'
3001  if (mype==0) write(*,'(a,f5.2,f6.2,f6.2,a)') ' xI2 vector: [',vec_xi2(1),vec_xi2(2),vec_xi2(3),']'
3002 
3003  end subroutine init_vectors
3004 
3005  subroutine get_cor_image(x_3D,x_image)
3006  double precision :: x_3D(1:3),x_image(1:2)
3007  double precision :: res,res_origin
3008 
3009  call dot_product_loc(x_3d,vec_xi1,res)
3010  call dot_product_loc(x_origin,vec_xi1,res_origin)
3011  x_image(1)=res-res_origin
3012  call dot_product_loc(x_3d,vec_xi2,res)
3013  call dot_product_loc(x_origin,vec_xi2,res_origin)
3014  x_image(2)=res-res_origin
3015 
3016  end subroutine get_cor_image
3017 
3018 end module mod_thermal_emission
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
Module for physical and numeric constants.
Definition: mod_constants.t:2
double precision, parameter kb_cgs
Boltzmann constant in cgs.
Definition: mod_constants.t:34
double precision, parameter half
Definition: mod_constants.t:20
double precision, parameter one
Definition: mod_constants.t:18
double precision, parameter dpi
Pi.
Definition: mod_constants.t:25
double precision, parameter zero
some frequently used numbers
Definition: mod_constants.t:17
double precision, parameter smalldouble
Definition: mod_constants.t:8
double precision, parameter mp_cgs
Proton mass in cgs.
Definition: mod_constants.t:28
double precision, parameter const_c
Definition: mod_constants.t:48
This module contains definitions of global parameters and variables and some generic functions/subrou...
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.
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
double precision global_time
The global simulation time.
integer snapshotini
Resume from the snapshot with this index.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
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.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
double precision image_rotate
rotation of image
character(len=std_len) resolution_sxr
resolution of the SXR image
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
double precision unit_velocity
Physical scaling factor for velocity.
double precision, dimension(ndim) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
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)
double precision spectrum_window_max
character(len=std_len) resolution_euv
resolution of the EUV image
integer wavelength
wavelength for output
integer, dimension(ndim) stretch_type
What kind of stretching is used per dimension.
character(len=std_len) resolution_spectrum
resolution of the spectrum
double precision spectrum_window_min
spectral window
integer refine_max_level
Maximal number of AMR levels.
integer direction_slit
direction of the slit
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
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...
Definition: mod_physics.t:4
double precision, dimension(1:3) vec_los
subroutine integrate_euv_data_resol(igrid, nXIF1, nXIF2, xIF1, xIF2, dxIF1, dxIF2, fl, EUV, Dpl)
double precision, dimension(1:101) f_304
double precision, dimension(1:101) f_193
double precision, dimension(1:60) f_264
double precision, dimension(1:60) f_263
subroutine get_euv_image(qunit, fl)
subroutine get_line_info(wl, ion, mass, logTe, line_center, spatial_px, spectral_px, sigma_PSF, width_slit)
double precision, dimension(1:60) t_eis1
subroutine integrate_spectra_data_resol(igrid, wL, dwL, spectra, numWL, numXS, dir_loc, fl)
double precision, dimension(1:101) f_131
subroutine get_image_inst_resol(qunit, datatype, fl)
double precision, dimension(1:60) f_255
subroutine integrate_sxr_inst_resol(igrid, numXI1, numXI2, xI1, xI2, dxI, fl, SXR)
subroutine write_image_vtucc(qunit, xC, wC, dxC, nPiece, nC1, nC2, nWC, datatype)
subroutine get_euv(wl, ixIL, ixOL, w, x, fl, flux)
subroutine get_cor_image(x_3D, x_image)
subroutine write_image_vticc(qunit, xO1, xO2, dxO1, dxO2, wO, nXO1, nXO2, nWO, nC1, nC2)
double precision, dimension(1:101) t_aia
subroutine get_spectrum_inst_resol(qunit, datatype, fl)
subroutine integrate_euv_inst_resol(igrid, numXI1, numXI2, xI1, xI2, dxI, fl, EUV, Dpl)
double precision, dimension(1:60) t_eis2
subroutine get_spectrum_data_resol(qunit, datatype, fl)
double precision, dimension(1:101) f_171
subroutine output_data(qunit, xO1, xO2, dxO1, dxO2, wO, nXO1, nXO2, nWO, datatype)
double precision, dimension(1:101) f_94
double precision, dimension(1:3) vec_xi1
subroutine get_image_data_resol(qunit, datatype, fl)
subroutine dot_product_loc(vec_in1, vec_in2, res)
subroutine get_sxr_image(qunit, fl)
double precision, dimension(1:60) f_192
subroutine integrate_sxr_data_resol(igrid, nXIF1, nXIF2, xIF1, xIF2, dxIF1, dxIF2, fl, SXR)
double precision, dimension(1:3) vec_xi2
double precision, dimension(1:41) f_1354
subroutine cross_product_loc(vec_in1, vec_in2, vec_out)
double precision, dimension(1:101) f_335
subroutine get_sxr(ixIL, ixOL, w, x, fl, flux, El, Eu)
double precision, dimension(1:41) t_iris
subroutine get_euv_spectrum(qunit, fl)
double precision, dimension(1:101) f_211
subroutine integrate_spectra_inst_resol(igrid, wL, dwLg, xS, dxSg, spectra, numWL, numXS, fl)
subroutine get_goes_flux_grid(ixIL, ixOL, w, x, dV, xboxL, xbL, fl, eflux_grid)
subroutine get_goes_sxr_flux(xboxL, fl, eflux)