-
Notifications
You must be signed in to change notification settings - Fork 0
/
Lineabund_SCAT.f
executable file
·103 lines (83 loc) · 3.15 KB
/
Lineabund_SCAT.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
subroutine lineabund_scat (abundin)
c******************************************************************************
c This routine iteratively determines the abundance for a single line;
c for ease and speed of computation the actual changes are done to
c gf-values, and then at the end the adjustment is made to input
c abundance by the amount that the gf had to be changed
c******************************************************************************
implicit real*8 (a-h,o-z)
include 'Atmos.com'
include 'Linex.com'
include 'Mol.com'
include 'Pstuff.com'
lim2 = lim1
ncurve = 1
c*****find the c-o-g gf that matches the observed RW
gf1(ncurve) = gf(lim1)
rwlgobs = dlog10(width(lim1)/wave1(lim1))
gfobs = gftab(ntabtot)
do i=2,ntabtot
if (rwtab(i) .gt. rwlgobs) then
gfobs = gftab(i-1) + (gftab(i)-gftab(i-1))*
. (rwlgobs-rwtab(i-1))/(rwtab(i)-rwtab(i-1))
go to 15
endif
enddo
c*****using the current gf, compute a line and find the c-o-g gf that matches
c*****the computed RW
15 call oneline_scat (1)
rwlgcal = dlog10(w(ncurve)/wave1(lim1))
gfcal = gftab(ntabtot)
do i=2,ntabtot
if (rwtab(i) .gt. rwlgcal) then
gfcal = gftab(i-1) + (gftab(i)-gftab(i-1))*
. (rwlgcal-rwtab(i-1))/(rwtab(i)-rwtab(i-1))
go to 20
endif
enddo
c*****are the observed and computed RWs too different? if so, iterate
c*****by adjusting the assumed gf, and recompute the line.
c*****for strong lines, the iterations are slowed down by using the
c*****square root of the proposed gf shift.
20 error = (w(ncurve)-width(lim1))/width(lim1)
ratio = 10.**(gfobs-gfcal)
ncurve = ncurve + 1
if (dabs(error) .ge. 0.0015 .and. ncurve .lt. 20) then
rwlcomp = dlog10(w(ncurve-1)/wave1(lim1))
if (rwlcomp .gt. -4.7) then
gf1(ncurve) = gf1(ncurve-1)*dsqrt(ratio)
do i=1,ntau
kapnu0(lim1,i) = kapnu0(lim1,i)*dsqrt(ratio)
enddo
else
gf1(ncurve) = gf1(ncurve-1)*ratio
do i=1,ntau
kapnu0(lim1,i) = kapnu0(lim1,i)*ratio
enddo
endif
go to 15
endif
c*****if the observed and computed RWs are close, do a final gf adjustment
c*****and finish with one more line recomputation
if (ncurve .eq. 30) then
write (nf1out,1001)
write (nf2out,1001)
endif
gf1(ncurve) = gf1(ncurve-1)*ratio
do i=1,ntau
kapnu0(lim1,i) = kapnu0(lim1,i)*ratio
enddo
call oneline_scat (2)
widout(lim1) = w(ncurve)
wid1comp(lim1) = w(1)
diff = dlog10(gf1(ncurve)/gf(lim1))
abundout(lim1) = abundin + diff
if (ncurve .ne. 1) then
write (nf1out,1002) ncurve
endif
return
c*****format statements
1001 format ('OH NO! ANOTHER FAILED ITERATION!')
1002 format (' This fit required ',i2,' iterations including ',
. 'a final small adjustment'/)
end