-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathInmodel.f
executable file
·462 lines (412 loc) · 16.7 KB
/
Inmodel.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
subroutine inmodel
c******************************************************************************
c This subroutine reads in the model
c******************************************************************************
implicit real*8 (a-h,o-z)
include 'Atmos.com'
include 'Linex.com'
include 'Mol.com'
include 'Quants.com'
include 'Factor.com'
include 'Dummy.com'
include 'Pstuff.com'
real*8 element(95), logepsilon(95)
real*8 kaprefmass(100)
real*8 bmol(110)
character list*80, list2*70
integer append
c*****Read in the key word to define the model type
modelnum = modelnum + 1
rewind nfmodel
read (nfmodel,2001) modtype
write (nf1out,1010) modtype
if (modtype .eq. 'begn ' .or. modtype .eq. 'BEGN ')
. write (nf1out,1011)
c*****Read a comment line (usually describing the model)
read (nfmodel,2002) moditle
c*****Read the number of depth points
read (nfmodel,2002) list
list2 = list(11:)
read (list2,*) ntau
if (ntau .gt. 100) then
write (array,1012)
call prinfo (10)
stop
endif
c*****EITHER: Read in a model from the output of the experimental new
c MARCS code. This modtype is called "NEWMARCS". On each line
c the numbers are:
c tau(5000), t, pe, pgas, rho, model microtrubulent velocity,
c and mean opacity (cm^2/gm) at the reference wavelength (5000A).
if (modtype .eq. 'NEWMARCS ') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) tauref(i),t(i),ne(i),pgas(i),rho(i),
. vturb(1),kaprefmass(i)
enddo
c*****OR: Read in a model from the output of the on-line new
c MARCS code. This modtype is called "WEBMARCS". On each line
c the numbers are:
c layer number (not needed), log{tau(Rosseland)} (not needed),
c log{tau(5000)}, depth, t, pe, pgas, prad (not read in) and
c pturb (not read in)
elseif (modtype .eq. 'WEBMARCS') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) k, dummy1(k), tauref(i), dummy2(k), t(i),
. ne(i), pgas(i)
enddo
c*****OR: Read in a model from an alternative form of on-line new
c MARCS code. This modtype is called "WEB2MARC". On each line
c the numbers are:
c atmospheric layer number (not needed), log{tau(5000)}, t,
c log(Pe), log(Pgas), rhox
elseif (modtype .eq. 'WEB2MARC') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) k,tauref(i),t(i),ne(i),pgas(i),rhox(i)
enddo
c OR: Read in a model from the output of the ATLAS code. This
c modtype is called "KURUCZ". On each line the numbers are:
c rhox, t, pgas, ne, and Rosseland mean opacity (cm^2/gm), and
c two numbers not used by MOOG.
elseif (modtype .eq. 'KURUCZ ') then
do i=1,ntau
read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i),kaprefmass(i)
enddo
c OR: Read in a model from the output of the NEXTGEN code. This
c modtype is called "NEXTGEN". These models have tau scaled at a
c specific wavelength that is read in before the model. MOOG will
c need to generate the opacities internally.On each line the numbers
c are: tau, t, pgas, pe, density, mean molecular weight, two numbers
c not used by MOOG, and Rosseland mean opacity (cm^2/gm).
elseif (modtype .eq. 'NEXTGEN ') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i), rho(i),
. molweight(i), x2, x3, kaprefmass(i)
enddo
c OR: Read in a model from the output of the MARCS code. This modtype
c type is called "BEGN". On each line the numbers are:
c tauross, t, log(pg), log(pe), mol weight, and kappaross.
elseif (modtype .eq. 'BEGN ') then
do i=1,ntau
read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i),
. molweight(i), kaprefmass(i)
enddo
c OR: Read in a model generated from ATLAS, but without accompanying
c opacities. MOOG will need to generate the opacities internally,
c using a reference wavelength that it reads in before the model.
elseif (modtype .eq. 'KURTYPE') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) rhox(i),t(i),pgas(i),ne(i)
enddo
c OR: Read in a model generated from ATLAS, with output generated
c in Padova. The columns are in somewhat different order than normal
elseif (modtype .eq. 'KUR-PADOVA') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) tauref(i),t(i),kaprefmass(i),
. ne(i),pgas(i),rho(i)
enddo
c OR: Read in a generic model that has a tau scale at a specific
c wavelength that is read in before the model.
c MOOG will need to generate the opacities internally.
elseif (modtype .eq. 'GENERIC ') then
read (nfmodel,*) wavref
do i=1,ntau
read (nfmodel,*) tauref(i),t(i),pgas(i),ne(i)
enddo
c OR: quit in utter confusion if those model types are not specified
else
write (*,1001)
stop
endif
c*****Compute other convenient forms of the temperatures
do i=1,ntau
theta(i) = 5040./t(i)
tkev(i) = 8.6171d-5*t(i)
tlog(i) = dlog(t(i))
enddo
c*****Convert from logarithmic Pgas scales, if needed
if (pgas(ntau)/pgas(1) .lt. 10.) then
do i=1,ntau
pgas(i) = 10.0**pgas(i)
enddo
endif
c*****Convert from logarithmic Ne scales, if needed
if(ne(ntau)/ne(1) .lt. 20.) then
do i=1,ntau
ne(i) = 10.0**ne(i)
enddo
endif
c*****Convert from Pe to Ne, if needed
if(ne(ntau) .lt. 1.0e7) then
do i=1,ntau
ne(i) = ne(i)/1.38054d-16/t(i)
enddo
endif
c*****compute the atomic partition functions
do j=1,95
elem(j) = dble(j)
call partfn (elem(j),j)
enddo
c*****Read the microturbulence (either a single value to apply to
c all layers, or a value for each of the ntau layers).
c Conversion to cm/sec from km/sec is done if needed
read (nfmodel,2003) (vturb(i),i=1,6)
if (vturb(2) .ne. 0.) then
read (nfmodel,2003) (vturb(i),i=7,ntau)
else
do i=2,ntau
vturb(i) = vturb(1)
enddo
endif
if (vturb(1) .lt. 100.) then
write (moditle(55:62),1008) vturb(1)
do i=1,ntau
vturb(i) = 1.0e5*vturb(i)
enddo
else
write (moditle(55:62),1008) vturb(1)/1.0e5
endif
c*****Read in the abundance data, storing the original abundances in xabu
c*****The abundances not read in explicity are taken from the default
c*****solar ones contained in array xsolar.
read (nfmodel,2002) list
list2 = list(11:)
read (list2,*) natoms,abscale
write (moditle(63:73),1009) abscale
if(natoms .ne. 0)
. read (nfmodel,*) (element(i),logepsilon(i),i=1,natoms)
xhyd = 10.0**xsolar(1)
xabund(1) = 1.0
xabund(2) = 10.0**xsolar(2)/xhyd
do i=3,95
xabund(i) = 10.0**(xsolar(i)+abscale)/xhyd
xabu(i) = xabund(i)
enddo
if (natoms .ne. 0) then
do i=1,natoms
xabund(idint(element(i))) = 10.0**logepsilon(i)/xhyd
xabu(idint(element(i))) = 10.0**logepsilon(i)/xhyd
enddo
endif
c*****Compute the mean molecular weight, ignoring molecule formation
c in this approximation (maybe make more general some day?)
wtnum = 0.
wtden = 0.
do i=1,95
wtnum = wtnum + xabund(i)*xam(i)
wtden = wtden + xabund(i)
enddo
wtmol = wtnum/(xam(1)*wtden)
nomolweight = 0
if (modtype .eq. 'BEGN ' .or. modtype .eq. 'NEXTGEN ') then
nomolweight = 1
endif
if (nomolweight .ne. 1) then
do i=1,ntau
molweight(i) = wtmol
enddo
endif
c*****Compute the density
if (modtype .ne. 'NEXTGEN ') then
do i=1,ntau
rho(i) = pgas(i)*molweight(i)*1.6606d-24/(1.38054d-16*t(i))
enddo
endif
c*****Calculate the fictitious number density of hydrogen
c Note: ph = (-b1 + dsqrt(b1*b1 - 4.0*a1*c1))/(2.0*a1)
iatom = 1
call partfn (dble(iatom),iatom)
do i=1,ntau
th = 5040.0/t(i)
ah2 = 10.0**(-(12.7422+(-5.1137+(0.1145-0.0091*th)*th)*th))
a1 = (1.0+2.0*xabund(2))*ah2
b1 = 1.0 + xabund(2)
c1 = -pgas(i)
ph = (-b1/2.0/a1)+dsqrt((b1**2/(4.0*a1*a1))-(c1/a1))
nhtot(i) = (ph+2.0*ph*ph*ah2)/(1.38054d-16*t(i))
enddo
c*****Molecular equilibrium called here.
c First, a default list of ions and molecules is considered. Then the
c user's list is read in. A check is performed to see if any of these
c species need to be added. If so, then they are appended to the
c default list. The molecular equilibrium routine is then called.
c Certain species are important for continuous opacities and damping
c calculations - these are read from the equilibrium solution and
c saved.
c*****Set up the default molecule list
if (molset .eq. 0) then
nmol = 21
else
nmol = 57
endif
if (molset .eq. 0) then
do i=1,110
amol(i) = smallmollist(i)
enddo
elseif (molset .eq. 1) then
do i=1,110
amol(i) = largemollist(i)
enddo
else
array = 'molset = 0 or 1 only; I quit!'
call putasci (29,6)
stop
endif
c*****Read in the names of additional molecules to be used in
c molecular equilibrium if needed.
read (nfmodel,2002,end=101) list
list2 = list(11:)
read (list2,*) moremol
if (moremol .ne. 0) then
read (nfmodel,*) (bmol(i),i=1,moremol)
append = 1
do k=1,moremol
do l=1,nmol
if (nint(bmol(k)) .eq. nint(amol(l)))
. append = 0
enddo
if (append .eq. 1) then
nmol = nmol + 1
amol(nmol) = bmol(k)
endif
append = 1
enddo
endif
c*****do the general molecular equilibrium
101 call eqlib
c In the number density array "numdens", the elements denoted by
c the first subscripts are named in at the ends of the assignment
c lines; at present these are the only ones needed for continuous
c opacities
c
do i=1,ntau
numdens(1,1,i) = xamol(1,i) H I
numdens(1,2,i) = xmol(1,i) H II
numdens(2,1,i) = xamol(2,i) He I
numdens(2,2,i) = xmol(2,i) Hi II
numdens(3,1,i) = xamol(3,i) C I
numdens(3,2,i) = xmol(3,i) C II
numdens(4,1,i) = xamol(6,i) Mg I
numdens(4,2,i) = xmol(6,i) Mg II
numdens(5,1,i) = xamol(7,i) Al I
numdens(5,2,i) = xmol(7,i) Al II
numdens(6,1,i) = xamol(8,i) Si I
numdens(6,2,i) = xmol(8,i) Si II
numdens(7,1,i) = xamol(16,i) Fe I
numdens(7,2,i) = xmol(16,i) Fe II
numdens(8,1,i) = xmol(17,i) H_2
enddo
c*****SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units
if (modtype .eq. 'NEWMARCS ') then
do i=1,ntau
kapref(i) = kaprefmass(i)*rho(i)
enddo
c SPECIAL NEEDS: for KURUCZ models, to create the optical depth array,
c and to convert kaprefs to our units
elseif (modtype .eq. 'KURUCZ ') then
first = rhox(1)*kaprefmass(1)
tottau = rinteg(rhox,kaprefmass,tauref,ntau,first)
tauref(1) = first
do i=2,ntau
tauref(i) = tauref(i-1) + tauref(i)
enddo
do i=1,ntau
kapref(i) = kaprefmass(i)*rho(i)
enddo
c SPECIAL NEEDS: for NEXTGEN models, to convert kaprefs to our units
elseif (modtype .eq. 'NEXTGEN ') then
do i=1,ntau
kapref(i) = kaprefmass(i)*rho(i)
enddo
c SPECIAL NEEDS: for BEGN models, to convert kaprefs to our units
elseif (modtype .eq. 'BEGN ') then
do i=1,ntau
kapref(i) = kaprefmass(i)*rho(i)
enddo
c SPECIAL NEEDS: for KURTYPE models, to create internal kaprefs,
c and to compute taurefs from the kaprefs converted to mass units
elseif (modtype .eq. 'KURTYPE ') then
call opacit (1,wavref)
do i=1,ntau
kaprefmass(i) = kapref(i)/rho(i)
enddo
first = rhox(1)*kaprefmass(1)
tottau = rinteg(rhox,kaprefmass,tauref,ntau,first)
tauref(1) = first
do i=2,ntau
tauref(i) = tauref(i-1) + tauref(i)
enddo
c SPECIAL NEEDS: for NEWMARCS models, to convert kaprefs to our units
elseif (modtype .eq. 'KUR-PADOVA') then
do i=1,ntau
kapref(i) = kaprefmass(i)*rho(i)
enddo
c SPECIAL NEEDS: for generic models, to create internal kaprefs,
elseif (modtype .eq. 'GENERIC ' .or.
. modtype .eq. 'WEBMARCS ' .or.
. modtype .eq. 'WEB2MARC ') then
call opacit (1,wavref)
endif
c*****Convert from logarithmic optical depth scales, or vice versa.
c xref will contain the log of the tauref
if(tauref(1) .lt. 0.) then
do i=1,ntau
xref(i) = tauref(i)
tauref(i) = 10.0**xref(i)
enddo
else
do i=1,ntau
xref(i) = dlog10(tauref(i))
enddo
endif
c*****Write information to output files
if (modprintopt .lt. 1) return
write (nf1out,1002) moditle
do i=1,ntau
dummy1(i) = dlog10(pgas(i))
dummy2(i) = dlog10(ne(i)*1.38054d-16*t(i))
enddo
write (nf1out,1003) wavref,(i,xref(i),tauref(i),t(i),dummy1(i),
. pgas(i),dummy2(i),ne(i),vturb(i),i=1,ntau)
write (nf1out,1004)
do i=1,95
dummy1(i) = dlog10(xabund(i)) + 12.0
enddo
write (nf1out,1005) (names(i),i,dummy1(i),i=1,95)
write (nf1out,1006) modprintopt, molopt, linprintopt, fluxintopt
write (nf1out,1007) (kapref(i),i=1,ntau)
return
c*****format statements
2001 format (a10)
2002 format (a80)
2003 format (6d13.0)
1001 format('permitted model types are:'/'KURUCZ, BEGN, ',
. 'KURTYPE, KUR-PADOVA, NEWMARCS, WEBMARCS, NEXTGEN, ',
. 'WEB2MARC, or GENERIC'/ 'MOOG quits!')
1002 format (/'MODEL ATMOSPHERE HEADER:'/a80/)
1003 format ('INPUT ATMOSPHERE QUANTITIES', 10x,
. '(reference wavelength =',f10.2,')'/3x, 'i', 2x, 'xref',
. 3x, 'tauref', 7x, 'T', 6x, 'logPg', 4x, 'Pgas',
. 6x, 'logPe', 5x, 'Ne', 9x, 'Vturb'/
. (i4, 0pf6.2, 1pd11.4, 0pf9.1, f8.3, 1pd11.4, 0pf8.3,
. 1pd11.4, d11.2))
1004 format (/'INPUT ABUNDANCES: (log10 number densities, log H=12)'/
. ' Default solar abundances: Asplund et al. 2009')
1005 format (5(3x,a2, '(',i2,')=', f5.2))
1006 format (/'OPTIONS: atmosphere = ', i1, 5x, 'molecules = ', i1/
. ' lines = ', i1, 5x, 'flux/int = ', i1)
1007 format (/'KAPREF ARRAY:'/(6(1pd12.4)))
1008 format ('vt=', f5.2)
1009 format (' M/H=', f5.2)
1010 format (13('-'),'MOOG OUTPUT FILE', 10('-'),
. '(MOOG version from 01/28/09)', 13('-')//
. 'THE MODEL TYPE: ', a10)
1011 format (' The Rosseland opacities and optical depths have ',
. 'been read in')
1012 format ('HOUSTON, WE HAVE MORE THAN 100 DEPTH POINTS! I QUIT!')
end