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