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
|
subroutine blends
c******************************************************************************
c This program does abundance derivations from blended spectral features;
c only one element will have its abundance derived per run.
c******************************************************************************
implicit real*8 (a-h,o-z)
include 'Atmos.com'
include 'Linex.com'
include 'Factor.com'
include 'Mol.com'
include 'Pstuff.com'
include 'Dummy.com'
c*****examine the parameter file
call params
c*****open the files for standard output and summary abundances
nf1out = 20
lscreen = 4
array = 'STANDARD OUTPUT'
nchars = 15
call infile ('output ',nf1out,'formatted ',0,nchars,
. f1out,lscreen)
nf2out = 21
lscreen = lscreen + 2
array = 'SUMMARY ABUNDANCE OUTPUT'
nchars = 24
call infile ('output ',nf2out,'formatted ',0,nchars,
. f2out,lscreen)
c*****open and read the model atmosphere
nfmodel = 30
lscreen = lscreen + 2
array = 'THE MODEL ATMOSPHERE'
nchars = 20
call infile ('input ',nfmodel,'formatted ',0,nchars,
. fmodel,lscreen)
call inmodel
c*****initialize some variables
isynth = 1
isorun = 1
iatom = nint(cogatom)
pec(iatom) = 1
numpecatom = 1
pecabund(iatom,1) = 0.
c*****open and read the line list file; get ready for the line calculations
nflines = 31
lscreen = lscreen + 2
array = 'THE LINE LIST'
nchars = 13
call infile ('input ',nflines,'formatted ',0,nchars,
. flines,lscreen)
100 call inlines (1)
call eqlib
call nearly (1)
c*****start the large loop that will go through each blended feature
ewsynthopt = -1
mode = 4
30 do lll=1,1000
if (lim2line .eq. nlines) exit
c*****define the set of lines responsible for a blended feature
call linlimit
lim1 = lim1line
lim2 = lim2line
write (99,1007) iatom, wave1(lim1), wave1(lim2)
c*****make sure that the element whose abundance is to be fit has
c a representative line of the blend
ifind = 0
do j=lim1,lim2
if (atom1(j) .lt. 100.) then
if (iatom .eq. int(atom1(j))) then
ifind = 1
exit
endif
else
call sunder (atom1(j),ia,ib)
if (iatom.eq.ia .or. iatom.eq.ib) then
ifind = 1
exit
endif
endif
enddo
if (ifind .eq. 0) then
do j=lim1,lim2
abundout(j) = 999.99
enddo
write (nf1out,1002)
lim1line = lim2line + 1
if (lim1line .le. nlines+nstrong) cycle
endif
c*****do the syntheses, forcing each abundance to predict the
c feature equivalent width; a limit of 30 iterations is imposed
c arbitrarily
ncurve = 1
start = wave1(lim1) - delwave
sstop = wave1(lim2) + delwave
delta = wave1(lim2) - wave1(lim1) + delwave
oldstart = start
oldstop = sstop
oldstep = step
olddelta = delta
rwlgobs = dlog10(width(lim1)/wave1(lim1))
gf1(ncurve) = 1.
do k=1,30
call synspec
call total
error = (w(ncurve)-width(lim1))/width(lim1)
ratio = width(lim1)/w(ncurve)
ncurve = ncurve + 1
c*****here we go for another iteration
if (dabs(error) .ge. 0.0075) then
rwlcomp = dlog10(w(ncurve-1)/wave1(lim1))
if (rwlcomp.lt.-5.2 .and. rwlgobs.lt.-5.2) then
ratio = ratio
elseif (rwlcomp.ge.-5.2 .and. rwlgobs.ge.-5.2) then
ratio = ratio**2.0
else
ratio = ratio**1.5
endif
gf1(ncurve) = gf1(ncurve-1)*ratio
do j=lim1,lim2
if (atom1(j) .gt. 100.) then
call sunder (atom1(j),ia,ib)
if (ia.eq.iatom .or. ib.eq.iatom) then
do i=1,ntau
kapnu0(j,i) = kapnu0(j,i)*ratio
enddo
endif
elseif (int(atom1(j)) .eq. iatom) then
do i=1,ntau
kapnu0(j,i) = kapnu0(j,i)*ratio
enddo
endif
enddo
if (k .eq. 20) then
write (*,1008)
stop
endif
else
exit
endif
enddo
c*****here we do the final calculation when the predicted and observed are close
gf1(ncurve) = gf1(ncurve-1)*ratio
do j=lim1,lim2
if (atom1(j) .gt. 100.) then
call sunder (atom1(j),ia,ib)
if (ia.eq.iatom .or. ib.eq.iatom) then
do i=1,ntau
kapnu0(j,i) = kapnu0(j,i)*ratio
enddo
endif
elseif (int(atom1(j)) .eq. iatom) then
do i=1,ntau
kapnu0(j,i) = kapnu0(j,i)*ratio
enddo
endif
enddo
call synspec
call total
widout(lim1) = w(ncurve)
diff = dlog10(gf1(ncurve))
abundout(lim1) = dlog10(xabund(iatom)) + 12.0 + diff
if (ncurve .ne. 1) then
write (nf1out,1001) ncurve
endif
c*****here is where some auxiliary things like mean depth are computed
if (lim1.eq.lim2 .and. linprintopt.ge.3) then
wave = wave1(lim1)
call taukap
call cdcalc(2)
first = 0.4343*cd(1)
d(1) = rinteg(xref,cd,dummy1,ntau,first)
do i=1,ntau
dummy1(i) = xref(i)*cd(i)
enddo
first = 0.
cdmean = rinteg(xref,dummy1,dummy2,ntau,first)/
. rinteg(xref,cd,dummy2,ntau,first)
do i=1,ntau
if (cdmean .lt. cd(i)) exit
enddo
write (nf1out,1005) lim1, cdmean, i, xref(i)
do i=1,ntau
if (taunu(i)+taulam(i) .ge. 1.) exit
enddo
write (nf1out,1006) lim1, i, dlog10(tauref(i)),
. dlog10(taulam(i)), dlog10(taunu(i))
endif
c*****assign the abundance to the strongest line of the blend that
c contains cogatom; put 999.99's as the abundances of all but
c this line; go back for another blended feature
if (lim2 .gt. lim1) then
abunblend = abundout(lim1)
widblend = widout(lim1)
strongest = 0.
linstrongest = 0
do j=lim1,lim2
abundout(j) = 999.99
enddo
do j=lim1,lim2
if (atom1(j) .lt. 100.) then
if (dint(atom1(j)) .eq. cogatom) then
if (kapnu0(j,jtau5) .gt. strongest) then
strongest = kapnu0(j,jtau5)
linstrongest = j
endif
endif
else
call sunder (atom1(j),ia,ib)
if (dble(ia).eq.cogatom .or. dble(ib).eq.cogatom) then
if (kapnu0(j,jtau5) .gt. strongest) then
strongest = kapnu0(j,jtau5)
linstrongest = j
endif
endif
endif
enddo
abundout(linstrongest) = abunblend
widout(linstrongest) = widblend
endif
lim1line = lim2line + 1
if (lim1line .le. nlines+nstrong) cycle
enddo
c*****do abundance statistics; print out a summary of the abundances
lim1obs = 1
lim2obs = nlines+nstrong
call stats
rewind nf2out
call lineinfo (3)
c*****here a plot may be made on the terminal (and paper) if there
c are enough lines; then the user will be prompted on some
c options concerning what is seen on the plot
if (plotopt .ne. 0) then
call pltabun
if (choice.eq.'v' .or. choice.eq.'m') go to 100
endif
c*****now the option will be to redo the molecular equilibrium and redo
c the last species, at the user's option.
if (neq .ne. 0) then
do n=1,neq
if (iatom .eq. iorder(n)) then
write (array,1004)
ikount = min0(nlines+11,maxline)
nchars = 70
35 call getasci (nchars,ikount)
choice = chinfo(1:1)
if (choice.eq.'y' .or. nchars.le.0) then
xabund(iatom) = 10.**(average-12.)
call eqlib
call nearly (1)
lim1line = 0
lim2line = 0
rewind nf2out
go to 30
elseif (choice .eq. 'n') then
call finish (0)
else
go to 35
endif
endif
enddo
endif
c*****exit the program
call finish (0)
c*****format statements
1001 format (' This fit required ',i2,' iterations'/)
1002 format ('WARNING: NO ELEMENT LINE TO VARY IN BLEND')
1004 format ('SPECIES USES MOLECULAR EQUILIBRIUM!',
. ' REDO WITH NEW ABUNDANCE ([y]/n)? ')
1005 format (/'LINE ', i5, ':',
. ' weighted mean line contribution function C_d =',
. f6.2/ ' which occurs near level ', i3,
. ' with log tauref = ', f6.2)
1006 format (/'LINE ', i5, ':',
. ' tau(total) is greater than 1 at level',i3/
. ' logs of tauref, taulam, taunu =', 3f6.2)
1007 format (/'VARYING THIS ELEMENT: ', i8/
. 'USING THE LINE GROUP IN THE RANGE: ', 2f10.3)
1008 format (' MAX OF 30 ITERATIONS REACHED; I QUIT!')
end
|
