MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_solar_atmosphere.t
Go to the documentation of this file.
1!> User can use subroutine get_atm_para to generate 1D solar stmosphere.
2!> User should provide heights (h), number density at h=0,
3!> number of points (nh), and the gravity (grav) at each point.
4!> User can select temperature profile.
5!> This subroutine will return density and pressure at each point.
6module mod_solar_atmosphere
7 use mod_global_parameters
8 use mod_comm_lib, only: mpistop
9 implicit none
10
11 double precision, parameter :: solar_gravity=2.74d4 ! cm s^-2
12 double precision :: h_valc(1:52),t_valc(1:52)
13 double precision :: h_hong(1:300),t_hong(1:300)
14 double precision :: h_fontenla(1:56),t_fontenla(1:56)
15 double precision :: h_alc7(1:120),t_alc7(1:120)
16 integer :: n_valc,n_hong,n_fontenla,n_alc7
17
18 data n_alc7 / 120 /
19
20 data h_alc7 / 0.0, 10.0, 20.0, 35.0, 50.0, &
21 75.0, 100.0, 125.0, 150.0, 175.0, &
22 200.0, 250.0, 300.0, 350.0, 400.0, &
23 450.0, 490.0, 525.0, 560.0, 615.0, &
24 660.0, 700.0, 750.0, 800.0, 854.0, &
25 900.0, 946.0, 971.0, 1003.0, 1032.0, &
26 1065.0, 1101.0, 1143.0, 1214.0, 1299.0, &
27 1398.0, 1520.0, 1617.0, 1722.0, 1820.0, &
28 1894.0, 1946.0, 1989.0, 2024.0, 2055.0, &
29 2083.0, 2098.0, 2110.0, 2120.0, 2126.0, &
30 2130.0, 2132.0, 2134.0, 2136.0, 2138.0, &
31 2141.9, 2145.2, 2147.1, 2148.7, 2150.1, &
32 2151.2, 2152.0, 2152.9, 2153.5, 2154.1, &
33 2154.7, 2155.5, 2156.2, 2156.8, 2157.6, &
34 2158.2, 2158.8, 2159.4, 2160.0, 2160.6, &
35 2161.3, 2162.3, 2163.3, 2164.6, 2166.0, &
36 2167.7, 2170.0, 2172.3, 2175.3, 2178.7, &
37 2182.4, 2186.1, 2191.5, 2196.5, 2201.9, &
38 2207.7, 2213.7, 2220.0, 2227.7, 2235.3, &
39 2243.3, 2251.6, 2259.4, 2269.7, 2282.4, &
40 2294.4, 2308.6, 2325.8, 2346.1, 2378.2, &
41 2418.7, 2457.9, 2495.6, 2535.5, 2578.2, &
42 2628.3, 2688.1, 2761.5, 2844.4, 2970.8, &
43 3133.8, 3425.0, 3752.7, 4295.9, 4968.6/
44
45 data t_alc7 / 6583, 6397, 6231, 6006, 5826, &
46 5607, 5431, 5288, 5165, 5080, &
47 5010, 4907, 4805, 4700, 4590, &
48 4485, 4435, 4410, 4400, 4435, &
49 4510, 4640, 4840, 5090, 5430, &
50 5720, 5969, 6100, 6225, 6315, &
51 6400, 6474, 6531, 6576, 6598, &
52 6610, 6623, 6633, 6643, 6652, &
53 6660, 6667, 6674, 6680, 6686, &
54 6694, 6700, 6706, 6718, 6740, &
55 6768, 6800, 6870, 6992, 7248, &
56 7950, 9115, 10980, 13200, 15760, &
57 18140, 20510, 23100, 25120, 27130, &
58 29500, 32260, 34580, 36870, 39400, &
59 41450, 43400, 45140, 46800, 48490, &
60 50140, 52690, 55020, 57790, 60790, &
61 63950, 68000, 71810, 76330, 81220, &
62 86120, 90640, 96860, 102300, 107800, &
63 113200, 118700, 124000, 130200, 135800, &
64 141600, 147200, 152300, 158600, 166000, &
65 172600, 180100, 188600, 198100, 212000, &
66 227800, 241900, 254400, 266700, 279000, &
67 292500, 307300, 324000, 341300, 365000, &
68 392000, 432900, 471600, 524300, 577400/
69
70
71
72 data n_fontenla / 56 /
73
74 data h_fontenla / -100.0, -90.0, -80.0, -70.0, -60.0, &
75 -50.0, -40.0, -30.0, -20.0, -10.0, &
76 -0.7, 9.5, 20.0, 35.0, 50.0, &
77 75.0, 100.0, 125.0, 150.0, 175.0, &
78 200.0, 250.0, 300.0, 347.0, 398.0, &
79 447.5, 485.0, 519.0, 570.0, 625.0, &
80 696.0, 772.0, 812.0, 838.0, 874.5, &
81 899.5, 926.0, 947.0, 965.0, 1020.0, &
82 1060.0, 1120.0, 1180.0, 1245.0, 1319.0, &
83 1447.0, 1571.0, 1687.0, 1801.0, 1911.0, &
84 2013.0, 2142.0, 2263.0, 2357.0, 2425.0, &
85 2472.0 /
86
87 data t_fontenla / 9460, 9220, 8940, 8620, 8280, &
88 7940, 7600, 7270, 6960, 6690, &
89 6490, 6310, 6150, 5950, 5780, &
90 5550, 5380, 5245, 5130, 5035, &
91 4960, 4835, 4740, 4660, 4580, &
92 4510, 4450, 4390, 4300, 4200, &
93 4080, 3955, 3890, 3850, 3800, &
94 3810, 3950, 4150, 4400, 5350, &
95 5770, 6120, 6270, 6390, 6450, &
96 6490, 6530, 6570, 6620, 6670, &
97 6720, 6770, 6830, 6910, 7030, &
98 7210 /
99
100
101
102 data n_valc / 52 /
103
104 data h_valc / -75, -50, -25, 0, 50, &
105 100, 150, 250, 350, 450, &
106 515, 555, 605, 655, 705, &
107 755, 855, 905, 980, 1065, &
108 1180, 1280, 1380, 1515, 1605, &
109 1785, 1925, 1990, 2016, 2050, &
110 2070, 2080, 2090, 2104, 2107, &
111 2109, 2113, 2115, 2120, 2129, &
112 2160, 2200, 2230, 2255, 2263, &
113 2267, 2271, 2274, 2280, 2290, &
114 2298, 2543 /
115
116 data t_valc / 8320, 7610, 6910, 6420, 5840, &
117 5455, 5180, 4780, 4465, 4220, &
118 4170, 4230, 4420, 4730, 5030, &
119 5280, 5650, 5755, 5925, 6040, &
120 6150, 6220, 6280, 6370, 6440, &
121 6630, 6940, 7160, 7360, 7660, &
122 7940, 8180, 8440, 9500, 10700, &
123 12300, 18500, 21000, 22500, 23000, &
124 23500, 24000, 24200, 24500, 25500, &
125 28000, 32000, 37000, 50000, 89100, &
126 141000, 447000 /
127
128
129 data n_hong / 300 /
130
131 data h_hong / -3.7616d+07, -3.7384d+07, -3.7045d+07, -3.6561d+07, -3.5880d+07,&
132 -3.4950d+07, -3.3718d+07, -3.2153d+07, -3.0250d+07, -2.8050d+07,&
133 -2.5627d+07, -2.3088d+07, -2.0556d+07, -1.8151d+07, -1.5972d+07,&
134 -1.4004d+07, -1.2234d+07, -1.0625d+07, -9.1510d+06, -7.7773d+06,&
135 -6.4599d+06, -5.1584d+06, -3.8455d+06, -2.4946d+06, -1.0577d+06,&
136 4.7767d+05, 2.1041d+06, 3.8237d+06, 5.6198d+06, 7.4817d+06,&
137 9.3930d+06, 1.1339d+07, 1.3304d+07, 1.5278d+07, 1.7250d+07,&
138 1.9213d+07, 2.1163d+07, 2.3099d+07, 2.5019d+07, 2.6922d+07,&
139 2.8809d+07, 3.0678d+07, 3.2532d+07, 3.4371d+07, 3.6196d+07,&
140 3.8009d+07, 3.9811d+07, 4.1605d+07, 4.3396d+07, 4.5186d+07,&
141 4.6977d+07, 4.8773d+07, 5.0576d+07, 5.2390d+07, 5.4218d+07,&
142 5.6065d+07, 5.7934d+07, 5.9833d+07, 6.1766d+07, 6.3736d+07,&
143 6.5744d+07, 6.7791d+07, 6.9880d+07, 7.2011d+07, 7.4185d+07,&
144 7.6403d+07, 7.8660d+07, 8.0954d+07, 8.3284d+07, 8.5647d+07,&
145 8.8038d+07, 9.0456d+07, 9.2896d+07, 9.5357d+07, 9.7838d+07,&
146 1.0033d+08, 1.0285d+08, 1.0537d+08, 1.0791d+08, 1.1046d+08,&
147 1.1302d+08, 1.1559d+08, 1.1818d+08, 1.2077d+08, 1.2337d+08,&
148 1.2599d+08, 1.2862d+08, 1.3126d+08, 1.3392d+08, 1.3659d+08,&
149 1.3927d+08, 1.4198d+08, 1.4470d+08, 1.4744d+08, 1.5019d+08,&
150 1.5295d+08, 1.5571d+08, 1.5847d+08, 1.6121d+08, 1.6389d+08,&
151 1.6648d+08, 1.6895d+08, 1.7122d+08, 1.7325d+08, 1.7499d+08,&
152 1.7642d+08, 1.7754d+08, 1.7838d+08, 1.7900d+08, 1.7944d+08,&
153 1.7974d+08, 1.7996d+08, 1.8010d+08, 1.8020d+08, 1.8027d+08,&
154 1.8032d+08, 1.8035d+08, 1.8037d+08, 1.8039d+08, 1.8041d+08,&
155 1.8042d+08, 1.8044d+08, 1.8045d+08, 1.8046d+08, 1.8047d+08,&
156 1.8049d+08, 1.8051d+08, 1.8053d+08, 1.8055d+08, 1.8059d+08,&
157 1.8063d+08, 1.8068d+08, 1.8074d+08, 1.8082d+08, 1.8092d+08,&
158 1.8103d+08, 1.8118d+08, 1.8135d+08, 1.8157d+08, 1.8182d+08,&
159 1.8214d+08, 1.8252d+08, 1.8299d+08, 1.8357d+08, 1.8428d+08,&
160 1.8515d+08, 1.8621d+08, 1.8749d+08, 1.8904d+08, 1.9089d+08,&
161 1.9307d+08, 1.9561d+08, 1.9852d+08, 2.0180d+08, 2.0544d+08,&
162 2.0941d+08, 2.1368d+08, 2.1819d+08, 2.2291d+08, 2.2781d+08,&
163 2.3284d+08, 2.3798d+08, 2.4320d+08, 2.4849d+08, 2.5383d+08,&
164 2.5922d+08, 2.6463d+08, 2.7008d+08, 2.7554d+08, 2.8102d+08,&
165 2.8651d+08, 2.9202d+08, 2.9753d+08, 3.0305d+08, 3.0858d+08,&
166 3.1412d+08, 3.1966d+08, 3.2520d+08, 3.3075d+08, 3.3630d+08,&
167 3.4186d+08, 3.4741d+08, 3.5297d+08, 3.5854d+08, 3.6410d+08,&
168 3.6967d+08, 3.7523d+08, 3.8080d+08, 3.8637d+08, 3.9195d+08,&
169 3.9752d+08, 4.0309d+08, 4.0867d+08, 4.1424d+08, 4.1982d+08,&
170 4.2540d+08, 4.3098d+08, 4.3656d+08, 4.4214d+08, 4.4772d+08,&
171 4.5330d+08, 4.5888d+08, 4.6446d+08, 4.7004d+08, 4.7563d+08,&
172 4.8121d+08, 4.8679d+08, 4.9238d+08, 4.9796d+08, 5.0355d+08,&
173 5.0913d+08, 5.1472d+08, 5.2030d+08, 5.2589d+08, 5.3148d+08,&
174 5.3706d+08, 5.4265d+08, 5.4824d+08, 5.5382d+08, 5.5941d+08,&
175 5.6500d+08, 5.7059d+08, 5.7617d+08, 5.8176d+08, 5.8735d+08,&
176 5.9294d+08, 5.9853d+08, 6.0412d+08, 6.0971d+08, 6.1530d+08,&
177 6.2089d+08, 6.2648d+08, 6.3207d+08, 6.3766d+08, 6.4325d+08,&
178 6.4884d+08, 6.5443d+08, 6.6002d+08, 6.6561d+08, 6.7120d+08,&
179 6.7679d+08, 6.8238d+08, 6.8798d+08, 6.9357d+08, 6.9916d+08,&
180 7.0475d+08, 7.1034d+08, 7.1593d+08, 7.2152d+08, 7.2712d+08,&
181 7.3271d+08, 7.3830d+08, 7.4389d+08, 7.4948d+08, 7.5508d+08,&
182 7.6067d+08, 7.6626d+08, 7.7185d+08, 7.7744d+08, 7.8304d+08,&
183 7.8863d+08, 7.9422d+08, 7.9981d+08, 8.0540d+08, 8.1100d+08,&
184 8.1659d+08, 8.2218d+08, 8.2777d+08, 8.3337d+08, 8.3896d+08,&
185 8.4455d+08, 8.5014d+08, 8.5574d+08, 8.6133d+08, 8.6692d+08,&
186 8.7252d+08, 8.7811d+08, 8.8370d+08, 8.8929d+08, 8.9489d+08,&
187 9.0048d+08, 9.0607d+08, 9.1166d+08, 9.1726d+08, 9.2285d+08,&
188 9.2844d+08, 9.3404d+08, 9.3963d+08, 9.4522d+08, 9.5082d+08,&
189 9.5641d+08, 9.6200d+08, 9.6759d+08, 9.7319d+08, 9.7878d+08,&
190 9.8437d+08, 9.8997d+08, 9.9556d+08, 1.0012d+09, 1.0067d+09 /
191
192 data t_hong / 1.0000d+04, 1.0000d+04, 1.0001d+04, 1.0002d+04, 1.0003d+04,&
193 1.0005d+04, 1.0008d+04, 1.0012d+04, 1.0018d+04, 1.0026d+04,&
194 1.0036d+04, 1.0047d+04, 1.0058d+04, 1.0011d+04, 9.6700d+03,&
195 9.3468d+03, 8.9643d+03, 8.6116d+03, 8.2357d+03, 7.8466d+03,&
196 7.4734d+03, 7.1441d+03, 6.8397d+03, 6.5263d+03, 6.2978d+03,&
197 6.1165d+03, 5.9274d+03, 5.7787d+03, 5.6233d+03, 5.4920d+03,&
198 5.3684d+03, 5.2602d+03, 5.1621d+03, 5.0768d+03, 4.9988d+03,&
199 4.9202d+03, 4.8435d+03, 4.7743d+03, 4.7107d+03, 4.6492d+03,&
200 4.5903d+03, 4.5333d+03, 4.4783d+03, 4.4272d+03, 4.3832d+03,&
201 4.3391d+03, 4.2948d+03, 4.2540d+03, 4.2277d+03, 4.2067d+03,&
202 4.1940d+03, 4.1851d+03, 4.1956d+03, 4.2262d+03, 4.2723d+03,&
203 4.3381d+03, 4.4204d+03, 4.5278d+03, 4.6445d+03, 4.7643d+03,&
204 4.8856d+03, 5.0037d+03, 5.1176d+03, 5.2250d+03, 5.3211d+03,&
205 5.4063d+03, 5.4888d+03, 5.5778d+03, 5.6443d+03, 5.7069d+03,&
206 5.7603d+03, 5.8132d+03, 5.8700d+03, 5.9167d+03, 5.9592d+03,&
207 5.9938d+03, 6.0256d+03, 6.0561d+03, 6.0805d+03, 6.1051d+03,&
208 6.1298d+03, 6.1515d+03, 6.1720d+03, 6.1904d+03, 6.2096d+03,&
209 6.2281d+03, 6.2458d+03, 6.2654d+03, 6.2877d+03, 6.3134d+03,&
210 6.3429d+03, 6.3754d+03, 6.4108d+03, 6.4497d+03, 6.4949d+03,&
211 6.5456d+03, 6.6066d+03, 6.6837d+03, 6.7794d+03, 6.8983d+03,&
212 7.0383d+03, 7.1943d+03, 7.3537d+03, 7.4943d+03, 7.6133d+03,&
213 7.7131d+03, 7.8049d+03, 7.8909d+03, 7.9677d+03, 8.0333d+03,&
214 8.0851d+03, 8.1219d+03, 8.1534d+03, 8.2225d+03, 8.4087d+03,&
215 8.7744d+03, 9.3187d+03, 1.0097d+04, 1.1128d+04, 1.2398d+04,&
216 1.3941d+04, 1.5857d+04, 1.8286d+04, 2.1320d+04, 2.4676d+04,&
217 2.7878d+04, 3.0947d+04, 3.4019d+04, 3.7182d+04, 4.0498d+04,&
218 4.4021d+04, 4.7801d+04, 5.1890d+04, 5.6344d+04, 6.1230d+04,&
219 6.6619d+04, 7.2587d+04, 7.9209d+04, 8.6564d+04, 9.4709d+04,&
220 1.0364d+05, 1.1329d+05, 1.2363d+05, 1.3469d+05, 1.4653d+05,&
221 1.5920d+05, 1.7275d+05, 1.8719d+05, 2.0250d+05, 2.1859d+05,&
222 2.3535d+05, 2.5260d+05, 2.7013d+05, 2.8767d+05, 3.0498d+05,&
223 3.2184d+05, 3.3811d+05, 3.5372d+05, 3.6862d+05, 3.8283d+05,&
224 3.9636d+05, 4.0926d+05, 4.2157d+05, 4.3333d+05, 4.4459d+05,&
225 4.5538d+05, 4.6575d+05, 4.7573d+05, 4.8535d+05, 4.9464d+05,&
226 5.0362d+05, 5.1232d+05, 5.2076d+05, 5.2895d+05, 5.3691d+05,&
227 5.4467d+05, 5.5222d+05, 5.5958d+05, 5.6677d+05, 5.7379d+05,&
228 5.8066d+05, 5.8738d+05, 5.9395d+05, 6.0040d+05, 6.0672d+05,&
229 6.1291d+05, 6.1900d+05, 6.2497d+05, 6.3084d+05, 6.3661d+05,&
230 6.4228d+05, 6.4786d+05, 6.5335d+05, 6.5876d+05, 6.6408d+05,&
231 6.6933d+05, 6.7450d+05, 6.7960d+05, 6.8463d+05, 6.8959d+05,&
232 6.9449d+05, 6.9932d+05, 7.0409d+05, 7.0880d+05, 7.1346d+05,&
233 7.1806d+05, 7.2260d+05, 7.2709d+05, 7.3154d+05, 7.3593d+05,&
234 7.4027d+05, 7.4457d+05, 7.4882d+05, 7.5303d+05, 7.5719d+05,&
235 7.6131d+05, 7.6539d+05, 7.6943d+05, 7.7343d+05, 7.7739d+05,&
236 7.8131d+05, 7.8518d+05, 7.8898d+05, 7.9272d+05, 7.9641d+05,&
237 8.0004d+05, 8.0362d+05, 8.0715d+05, 8.1063d+05, 8.1406d+05,&
238 8.1745d+05, 8.2079d+05, 8.2409d+05, 8.2734d+05, 8.3056d+05,&
239 8.3373d+05, 8.3687d+05, 8.3997d+05, 8.4303d+05, 8.4606d+05,&
240 8.4906d+05, 8.5205d+05, 8.5502d+05, 8.5796d+05, 8.6090d+05,&
241 8.6381d+05, 8.6670d+05, 8.6958d+05, 8.7244d+05, 8.7528d+05,&
242 8.7810d+05, 8.8091d+05, 8.8371d+05, 8.8648d+05, 8.8924d+05,&
243 8.9199d+05, 8.9472d+05, 8.9743d+05, 9.0013d+05, 9.0281d+05,&
244 9.0548d+05, 9.0814d+05, 9.1078d+05, 9.1341d+05, 9.1602d+05,&
245 9.1862d+05, 9.2121d+05, 9.2378d+05, 9.2634d+05, 9.2889d+05,&
246 9.3142d+05, 9.3394d+05, 9.3645d+05, 9.3895d+05, 9.4143d+05,&
247 9.4391d+05, 9.4637d+05, 9.4882d+05, 9.5126d+05, 9.5368d+05,&
248 9.5610d+05, 9.5850d+05, 9.6089d+05, 9.6328d+05, 9.6565d+05,&
249 9.6801d+05, 9.7036d+05, 9.7270d+05, 9.7503d+05, 9.7735d+05,&
250 9.7965d+05, 9.8195d+05, 9.8424d+05, 9.8652d+05, 9.8879d+05,&
251 9.9105d+05, 9.9330d+05, 9.9554d+05, 9.9778d+05, 9.9778d+05 /
252
253
254contains
255
256 subroutine get_atm_para(h,rho,pth,grav,nh,Tcurve,hc,rhohc,Tem)
257 use mod_physics, only: phys_partial_ionization
258 use mod_ionization_degree
259 ! input:h,grav,nh,rho0,Tcurve; output:rho,pth (dimensionless units)
260 ! nh -- number of points
261 ! rho0 -- number density at h=0
262 ! Tcurve -- 'VAL-C' | 'Hong2017' | 'SPRM305' | 'AL-C7'
263
264 integer, intent(in) :: nh
265 double precision, intent(in) :: h(nh),grav(nh)
266 double precision, intent(out) :: rho(nh),pth(nh)
267 double precision, intent(out),optional :: Tem(nh)
268 double precision :: rhohc,hc
269 character(*) :: Tcurve
270
271 double precision :: h_cgs(nh),Te_cgs(nh),Te(nh)
272 double precision :: invT,dh,rhot,dht,ratio
273 integer :: j
274
275 h_cgs=h*unit_length
276
277 select case(tcurve)
278 case('VAL-C')
279 call get_te_valc(h_cgs,te_cgs,nh)
280 if (mype==0) print *, 'Temperature curve from Vernazza, Avrett & Loeser, 1981, ApJS, 45, 635'
281
282 case('Hong2017')
283 call get_te_hong(h_cgs,te_cgs,nh)
284 if (mype==0) print *, 'Temperature curve from Hong et al. 2017, ApJ, 845, 144'
285
286 case('Fontenla')
287 call get_te_sprm(h_cgs,te_cgs,nh)
288 if (mype==0) print *, 'Temperature curve from Fontenla et al. 2007, ApJ, 667, 1243'
289
290 case('AL-C7')
291 call get_te_alc7(h_cgs,te_cgs,nh)
292 if (mype==0) print *, 'Temperature curve from Avrett & Loeser 2008, ApJS, 175, 229'
293
294 case default
295 call mpistop("Unknown temperature curve")
296
297 end select
298
299 te=te_cgs/unit_temperature
300 if(present(tem)) tem=te
301
302 if(phys_partial_ionization) then
303 do j=1,nh
304 te(j)=te(j)*r_ideal_gas_law_partial_ionization(te(j))
305 end do
306 end if
307
308 ! density and pressure profiles
309 rho(1)=1.d5
310 pth(1)=rho(1)*te(1)
311
312 rhot=1.d0
313 invt=0.d0
314 do j=2,nh
315 dh=h(j)-h(j-1)
316 invt=invt+dh*(grav(j)/te(j)+grav(j-1)/te(j-1))*0.5d0
317 pth(j)=pth(1)*dexp(invt)
318 rho(j)=pth(j)/te(j)
319
320 if (h(j-1)<=hc .and. h(j)>hc) then
321 dht=hc-h(j-1)
322 rhot=rho(j-1)+dht*(rho(j)-rho(j-1))/(h(j)-h(j-1))
323 end if
324 end do
325
326 ratio=rhohc/rhot
327 rho=rho*ratio
328 pth=pth*ratio
329
330 end subroutine get_atm_para
331
332 subroutine get_te_alc7(h,Te,nh)
333 use mod_interpolation
334
335 integer :: nh
336 double precision :: h(nh),Te(nh)
337
338 integer :: ih,j,imin,imax,n_table
339 double precision :: h_table(n_alc7),T_table(n_alc7)
340 double precision :: unit_h,unit_T,htra,Ttra,Fc,kappa,dTdh
341
342 ! temperature profile
343 unit_h=1.d5 ! km -> cm
344 unit_t=1.0 ! K -> K
345 fc=1.72d5
346 kappa=8.d-7
347
348 n_table=n_alc7
349 h_table=h_alc7*unit_h
350 t_table=t_alc7*unit_t
351 htra=4968.6d0*unit_h
352 ttra=577400.d0
353 ! adiabatic temperature gradient below solar surface (Toriumi and Takasao 2017 ApJ 850 39)
354 dtdh=-0.4d0*mh_cgs*solar_gravity/kb_cgs
355
356 do ih=1,nh
357 if (h(ih)>h_table(n_table)) then
358 ! above transition region
359 te(ih)=(3.5d0*fc*(h(ih)-htra)/kappa+ttra**3.5d0)**(2.d0/7.d0)
360 endif
361
362 if (h(ih)<=h_table(1)) then
363 ! below photosphere
364 te(ih)=(h(ih)-h_table(1))*dtdh+t_table(1)
365 endif
366 enddo
367
368 ! inside the table
369 imin=nh
370 imax=nh-1
371 if (h(1)>=h_table(1) .and. h(1)<=h_table(n_table)) then
372 imin=1
373 else
374 do ih=2,nh
375 if (h(ih-1)<h_table(1) .and. h(ih)>=h_table(1) .and. h(ih)<=h_table(n_table)) imin=ih
376 enddo
377 endif
378 if (h(nh)>=h_table(1) .and. h(nh)<=h_table(n_table)) then
379 imax=nh
380 else
381 do ih=1,nh-1
382 if (h(ih)<=h_table(n_table) .and. h(ih+1)>h_table(n_table) .and. h(ih)>=h_table(1)) imax=ih
383 enddo
384 endif
385
386 if (imin<=imax) then
387 call interp_linear(h_table,t_table,n_table,h(imin:imax),te(imin:imax),imax-imin+1)
388 endif
389
390 end subroutine get_te_alc7
391
392 subroutine get_te_sprm(h,Te,nh)
393 use mod_interpolation
394
395 integer :: nh
396 double precision :: h(nh),Te(nh)
397
398 double precision :: h_table(n_fontenla),T_table(n_fontenla)
399 double precision :: unit_h,unit_T,dTdh
400 double precision :: h1,h2,h3
401 double precision :: Tpho,Ttop,htanh,wtra
402 double precision :: htra,Ttra,Fc,kappa
403 integer :: ih,j,imin,imax,n_table
404
405 unit_h=1.d5 ! km -> cm
406 unit_t=1.0 ! K -> K
407
408 n_table=n_fontenla
409 h_table=h_fontenla*unit_h
410 t_table=t_fontenla*unit_t
411
412 ! height for shift curve table/function
413 h1=h_table(1)
414 h2=h_table(n_table)
415 h3=3271.d0*unit_h
416
417 ! parameter for T curve in (h2,h3)
418 htanh=3464.d0*unit_h
419 tpho=6726.d0
420 ttop=1.5d6
421 wtra=246.d0*unit_h
422
423 ! parameter for T curve above h3
424 htra=3271.d0*unit_h
425 ttra=2.642d5
426 kappa=8.d-7
427 fc=4.791d5
428
429 ! parameter for T curve below h1
430 dtdh=(t_table(2)-t_table(1))/(h_table(2)-h_table(1))
431
432
433 do ih=1,nh
434 ! below photosphere
435 if (h(ih)<=h1) then
436 te(ih)=(h(ih)-h_table(1))*dtdh+t_table(1)
437 endif
438
439 ! high chromosphere and low transition region
440 if (h(ih)>h2 .and. h(ih)<h3) then
441 te(ih)=tpho+0.5d0*(ttop-tpho)*(tanh((h(ih)-htanh)/wtra)+1.d0)
442 endif
443
444 ! high transition region and corona
445 if (h(ih)>=h3) then
446 te(ih)=(3.5d0*fc*(h(ih)-htra)/kappa+ttra**3.5)**(2.d0/7.d0)
447 endif
448 enddo
449
450
451 ! inside the table
452 imin=nh
453 imax=nh-1
454 if (h(1)>=h_table(1) .and. h(1)<=h_table(n_table)) then
455 imin=1
456 else
457 do ih=2,nh
458 if (h(ih-1)<h_table(1) .and. h(ih)>=h_table(1) .and. h(ih)<=h_table(n_table)) imin=ih
459 enddo
460 endif
461 if (h(nh)>=h_table(1) .and. h(nh)<=h_table(n_table)) then
462 imax=nh
463 else
464 do ih=1,nh-1
465 if (h(ih)<=h_table(n_table) .and. h(ih+1)>h_table(n_table) .and. h(ih)>=h_table(1)) imax=ih
466 enddo
467 endif
468
469 if (imin<=imax) then
470 call interp_linear(h_table,t_table,n_fontenla,h(imin:imax),te(imin:imax),imax-imin+1)
471 endif
472
473 end subroutine
474
475 subroutine get_te_valc(h,Te,nh)
476 use mod_interpolation
477
478 integer :: nh
479 double precision :: h(nh),Te(nh)
480
481 double precision :: h_table(n_valc),T_table(n_valc)
482 double precision :: unit_h,unit_T,htra,Ttra,Fc,kappa,dTdh
483 integer :: ih,j,imin,imax,n_table
484
485 ! temperature profile
486 unit_h=1.d5 ! km -> cm
487 unit_t=1.0 ! K -> K
488 fc=5.38d5
489 kappa=8.d-7
490
491 n_table=n_valc
492 h_table=h_valc*unit_h
493 t_table=t_valc*unit_t
494 htra=h_table(n_table)
495 ttra=t_table(n_table)
496 dtdh=(t_table(2)-t_table(1))/(h_table(2)-h_table(1))
497
498 do ih=1,nh
499 if (h(ih)>h_table(n_table)) then
500 ! above transition region
501 te(ih)=(3.5d0*fc*(h(ih)-htra)/kappa+ttra**3.5d0)**(2.d0/7.d0)
502 endif
503
504 if (h(ih)<=h_table(1)) then
505 ! below photosphere
506 te(ih)=(h(ih)-h_table(1))*dtdh+t_table(1)
507 endif
508 enddo
509
510
511 ! inside the table
512 imin=nh
513 imax=nh-1
514 if (h(1)>=h_table(1) .and. h(1)<=h_table(n_table)) then
515 imin=1
516 else
517 do ih=2,nh
518 if (h(ih-1)<h_table(1) .and. h(ih)>=h_table(1) .and. h(ih)<=h_table(n_table)) imin=ih
519 enddo
520 endif
521 if (h(nh)>=h_table(1) .and. h(nh)<=h_table(n_table)) then
522 imax=nh
523 else
524 do ih=1,nh-1
525 if (h(ih)<=h_table(n_table) .and. h(ih+1)>h_table(n_table) .and. h(ih)>=h_table(1)) imax=ih
526 enddo
527 endif
528
529 if (imin<=imax) then
530 call interp_linear(h_table,t_table,n_valc,h(imin:imax),te(imin:imax),imax-imin+1)
531 endif
532
533 end subroutine get_te_valc
534
535 subroutine get_te_hong(h,Te,nh)
536 use mod_interpolation
537
538 integer :: nh
539 double precision :: h(nh),Te(nh)
540
541 double precision :: h_table(n_hong),T_table(n_hong)
542 double precision :: dTdh
543 double precision :: htra,Ttra,Fc,kappa
544 integer :: ih,j,imin,imax,n_table
545
546 n_table=n_hong
547 h_table=h_hong
548 t_table=t_hong
549
550 htra=h_table(150)
551 ttra=t_table(150)
552 fc=2.76d5 ! heat flux in high corona
553 kappa=8.d-7
554
555 do ih=1,nh
556 if (h(ih)>=h_table(n_table)) then
557 ! above max height of the table
558 te(ih)=(3.5d0*fc*(h(ih)-htra)/kappa+ttra**3.5d0)**(2.d0/7.d0)
559 endif
560
561 if (h(ih)<=h_table(1)) then
562 ! below min height of the table
563 dtdh=(t_table(2)-t_table(1))/(h_table(2)-h_table(1))
564 te(ih)=(h(ih)-h_table(1))*dtdh+t_table(1)
565 endif
566 enddo
567
568
569 ! inside the table
570 imin=nh
571 imax=nh-1
572 if (h(1)>=h_table(1) .and. h(1)<=h_table(n_table)) then
573 imin=1
574 else
575 do ih=2,nh
576 if (h(ih-1)<h_table(1) .and. h(ih)>=h_table(1) .and. h(ih)<=h_table(n_table)) imin=ih
577 enddo
578 endif
579 if (h(nh)>=h_table(1) .and. h(nh)<=h_table(n_table)) then
580 imax=nh
581 else
582 do ih=1,nh-1
583 if (h(ih)<=h_table(n_table) .and. h(ih+1)>h_table(n_table) .and. h(ih)>=h_table(1)) imax=ih
584 enddo
585 endif
586
587 if (imin<=imax) then
588 call interp_linear(h_table,t_table,n_table,h(imin:imax),te(imin:imax),imax-imin+1)
589 endif
590
591 end subroutine get_te_hong
592
593end module mod_solar_atmosphere
mod_comm_lib
Definition mod_comm_lib.t:1
mod_comm_lib::mpistop
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition mod_comm_lib.t:208
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::unit_length
double precision unit_length
Physical scaling factor for length.
Definition mod_global_parameters.t:72
mod_global_parameters::mype
integer mype
The rank of the current MPI task.
Definition mod_global_parameters.t:223
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_global_parameters::unit_temperature
double precision unit_temperature
Physical scaling factor for temperature.
Definition mod_global_parameters.t:84
mod_interpolation
Definition mod_interpolation.t:1
mod_interpolation::interp_linear
subroutine interp_linear(x_table, y_table, n_table, x_itp, y_itp, n_itp)
Definition mod_interpolation.t:9
mod_ionization_degree
module ionization degree - get ionization degree for given temperature
Definition mod_ionization_degree.t:2
mod_ionization_degree::r_ideal_gas_law_partial_ionization
double precision function r_ideal_gas_law_partial_ionization(te)
Definition mod_ionization_degree.t:177
mod_physics
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
mod_physics::phys_partial_ionization
logical phys_partial_ionization
Solve partially ionized one-fluid plasma.
Definition mod_physics.t:44
mod_solar_atmosphere
User can use subroutine get_atm_para to generate 1D solar stmosphere. User should provide heights (h)...
Definition mod_solar_atmosphere.t:6
mod_solar_atmosphere::t_alc7
double precision, dimension(1:120) t_alc7
Definition mod_solar_atmosphere.t:15
mod_solar_atmosphere::n_fontenla
integer n_fontenla
Definition mod_solar_atmosphere.t:16
mod_solar_atmosphere::get_te_sprm
subroutine get_te_sprm(h, te, nh)
Definition mod_solar_atmosphere.t:393
mod_solar_atmosphere::t_hong
double precision, dimension(1:300) t_hong
Definition mod_solar_atmosphere.t:13
mod_solar_atmosphere::t_fontenla
double precision, dimension(1:56) t_fontenla
Definition mod_solar_atmosphere.t:14
mod_solar_atmosphere::get_atm_para
subroutine get_atm_para(h, rho, pth, grav, nh, tcurve, hc, rhohc, tem)
Definition mod_solar_atmosphere.t:257
mod_solar_atmosphere::t_valc
double precision, dimension(1:52) t_valc
Definition mod_solar_atmosphere.t:12
mod_solar_atmosphere::n_hong
integer n_hong
Definition mod_solar_atmosphere.t:16
mod_solar_atmosphere::get_te_valc
subroutine get_te_valc(h, te, nh)
Definition mod_solar_atmosphere.t:476
mod_solar_atmosphere::h_valc
double precision, dimension(1:52) h_valc
Definition mod_solar_atmosphere.t:12
mod_solar_atmosphere::h_hong
double precision, dimension(1:300) h_hong
Definition mod_solar_atmosphere.t:13
mod_solar_atmosphere::h_fontenla
double precision, dimension(1:56) h_fontenla
Definition mod_solar_atmosphere.t:14
mod_solar_atmosphere::n_valc
integer n_valc
Definition mod_solar_atmosphere.t:16
mod_solar_atmosphere::h_alc7
double precision, dimension(1:120) h_alc7
Definition mod_solar_atmosphere.t:15
mod_solar_atmosphere::get_te_hong
subroutine get_te_hong(h, te, nh)
Definition mod_solar_atmosphere.t:536
mod_solar_atmosphere::n_alc7
integer n_alc7
Definition mod_solar_atmosphere.t:16
mod_solar_atmosphere::solar_gravity
double precision, parameter solar_gravity
Definition mod_solar_atmosphere.t:11
mod_solar_atmosphere::get_te_alc7
subroutine get_te_alc7(h, te, nh)
Definition mod_solar_atmosphere.t:333
Generated by doxygen 1.9.8