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