Nbody6
 All Files Functions Variables
zcnsts.f
Go to the documentation of this file.
1 ***
2  SUBROUTINE zcnsts(z,zpars)
3 *
4  implicit none
5  integer kw
6 *
7  include 'zdata.h'
8  real*8 z,zpars(20)
9  real*8 tm,tn,tscls(20),lums(10),gb(10)
10  real*8 lzs,dlzs,lz,lzd,dum1,m1,m2,rr,rb,mhefl,lhefl,thefl,lx
11  real*8 tbgbf,thef,lbagbf,lheif,lhef,lzahbf
12  real*8 rgbf,ragbf,rminf,mcgbf
13  external tbgbf,thef,lbagbf,lheif,lhef,lzahbf
14  external rgbf,ragbf,rminf,mcgbf
15 *
16  real*8 msp(200),gbp(200),c(5)
17  common /mscff/ msp
18  common /gbcff/ gbp
19  data c /3.040581d-01, 8.049509d-02, 8.967485d-02,
20  & 8.780198d-02, 2.219170d-02/
21 *
22 * ------------------------------------------------------------
23 *
24 * zpars: 1; M below which hook doesn't appear on MS, Mhook.
25 * 2; M above which He ignition occurs non-degenerately, Mhef.
26 * 3; M above which He ignition occurs on the HG, Mfgb.
27 * 4; M below which C/O ignition doesn't occur, Mup.
28 * 5; M above which C ignites in the centre, Mec.
29 * 6; value of log D for M<= zpars(3)
30 * 7; value of x for Rgb propto M^(-x)
31 * 8; value of x for tMS = MAX(tHOOK,x*tBGB)
32 * 9; constant for McHeIf when computing Mc,BGB, mchefl.
33 * 10; constant for McHeIf when computing Mc,HeI, mchefl.
34 * 11; hydrogen abundance.
35 * 12; helium abundance.
36 * 13; constant x in rmin = rgb*x**y used by LM CHeB.
37 * 14; z**0.4 to be used for WD L formula.
38 *
39 * ------------------------------------------------------------
40 *
41  lzs = log10(z/0.02d0)
42  dlzs = 1.d0/(z*log(10.d0))
43  lz = log10(z)
44  lzd = lzs + 1.d0
45 *
46  zpars(1) = 1.0185d0 + lzs*(0.16015d0 + lzs*0.0892d0)
47  zpars(2) = 1.995d0 + lzs*(0.25d0 + lzs*0.087d0)
48  zpars(3) = 16.5d0*z**0.06d0/(1.d0 + (1.0d-04/z)**1.27d0)
49  zpars(4) = max(6.11044d0 + 1.02167d0*lzs, 5.d0)
50  zpars(5) = zpars(4) + 1.8d0
51  zpars(6) = 5.37d0 + lzs*0.135d0
52  zpars(7) = c(1) + lzs*(c(2) + lzs*(c(3) + lzs*(c(4) + lzs*c(5))))
53  zpars(8) = max(0.95d0,max(0.95d0-(10.d0/3.d0)*(z-0.01d0),
54  & min(0.99d0,0.98d0-(100.d0/7.d0)*(z-0.001d0))))
55 ***
56 * Lzams
57  msp(1) = xz(1)+lzs*(xz(2)+lzs*(xz(3)+lzs*(xz(4)+lzs*xz(5))))
58  msp(2) = xz(6)+lzs*(xz(7)+lzs*(xz(8)+lzs*(xz(9)+lzs*xz(10))))
59  msp(3) = xz(11)+lzs*(xz(12)+lzs*(xz(13)+lzs*(xz(14)+lzs*xz(15))))
60  msp(4) = xz(16)+lzs*(xz(17)+lzs*(xz(18)+lzs*(xz(19)+lzs*xz(20))))
61  msp(5) = xz(21)+lzs*(xz(22)+lzs*(xz(23)+lzs*(xz(24)+lzs*xz(25))))
62  msp(6) = xz(26)+lzs*(xz(27)+lzs*(xz(28)+lzs*(xz(29)+lzs*xz(30))))
63  msp(7) = xz(31)+lzs*(xz(32)+lzs*(xz(33)+lzs*(xz(34)+lzs*xz(35))))
64 * Rzams
65  msp(8) = xz(36)+lzs*(xz(37)+lzs*(xz(38)+lzs*(xz(39)+lzs*xz(40))))
66  msp(9) = xz(41)+lzs*(xz(42)+lzs*(xz(43)+lzs*(xz(44)+lzs*xz(45))))
67  msp(10) = xz(46)+lzs*(xz(47)+lzs*(xz(48)+lzs*(xz(49)+lzs*xz(50))))
68  msp(11) = xz(51)+lzs*(xz(52)+lzs*(xz(53)+lzs*(xz(54)+lzs*xz(55))))
69  msp(12) = xz(56)+lzs*(xz(57)+lzs*(xz(58)+lzs*(xz(59)+lzs*xz(60))))
70  msp(13) = xz(61)
71  msp(14) = xz(62)+lzs*(xz(63)+lzs*(xz(64)+lzs*(xz(65)+lzs*xz(66))))
72  msp(15) = xz(67)+lzs*(xz(68)+lzs*(xz(69)+lzs*(xz(70)+lzs*xz(71))))
73  msp(16) = xz(72)+lzs*(xz(73)+lzs*(xz(74)+lzs*(xz(75)+lzs*xz(76))))
74 * Tbgb
75  msp(17) = xt(1)+lzs*(xt(2)+lzs*(xt(3)+lzs*xt(4)))
76  msp(18) = xt(5)+lzs*(xt(6)+lzs*(xt(7)+lzs*xt(8)))
77  msp(19) = xt(9)+lzs*(xt(10)+lzs*(xt(11)+lzs*xt(12)))
78  msp(20) = xt(13)+lzs*(xt(14)+lzs*(xt(15)+lzs*xt(16)))
79  msp(21) = xt(17)
80 * dTbgb/dz
81  msp(117) = dlzs*(xt(2)+lzs*(2.d0*xt(3)+3.d0*lzs*xt(4)))
82  msp(118) = dlzs*(xt(6)+lzs*(2.d0*xt(7)+3.d0*lzs*xt(8)))
83  msp(119) = dlzs*(xt(10)+lzs*(2.d0*xt(11)+3.d0*lzs*xt(12)))
84  msp(120) = dlzs*(xt(14)+lzs*(2.d0*xt(15)+3.d0*lzs*xt(16)))
85 * Thook
86  msp(22) = xt(18)+lzs*(xt(19)+lzs*(xt(20)+lzs*xt(21)))
87  msp(23) = xt(22)
88  msp(24) = xt(23)+lzs*(xt(24)+lzs*(xt(25)+lzs*xt(26)))
89  msp(25) = xt(27)+lzs*(xt(28)+lzs*(xt(29)+lzs*xt(30)))
90  msp(26) = xt(31)
91 * Ltms
92  msp(27) = xl(1)+lzs*(xl(2)+lzs*(xl(3)+lzs*(xl(4)+lzs*xl(5))))
93  msp(28) = xl(6)+lzs*(xl(7)+lzs*(xl(8)+lzs*(xl(9)+lzs*xl(10))))
94  msp(29) = xl(11)+lzs*(xl(12)+lzs*(xl(13)+lzs*xl(14)))
95  msp(30) = xl(15)+lzs*(xl(16)+lzs*(xl(17)+lzs*(xl(18)+lzs*xl(19))))
96  msp(27) = msp(27)*msp(30)
97  msp(28) = msp(28)*msp(30)
98  msp(31) = xl(20)+lzs*(xl(21)+lzs*(xl(22)+lzs*xl(23)))
99  msp(32) = xl(24)+lzs*(xl(25)+lzs*(xl(26)+lzs*xl(27)))
100 * Lalpha
101  m2 = 2.d0
102  msp(33) = xl(28)+lzs*(xl(29)+lzs*(xl(30)+lzs*xl(31)))
103  msp(34) = xl(32)+lzs*(xl(33)+lzs*(xl(34)+lzs*xl(35)))
104  msp(35) = xl(36)+lzs*(xl(37)+lzs*(xl(38)+lzs*xl(39)))
105  msp(36) = xl(40)+lzs*(xl(41)+lzs*(xl(42)+lzs*xl(43)))
106  msp(37) = max(0.9d0,1.1064d0+lzs*(0.415d0+0.18d0*lzs))
107  msp(38) = max(1.d0,1.19d0+lzs*(0.377d0+0.176d0*lzs))
108  if(z.gt.0.01d0)then
109  msp(37) = min(msp(37),1.d0)
110  msp(38) = min(msp(38),1.1d0)
111  endif
112  msp(39) = max(0.145d0,0.0977d0-lzs*(0.231d0+0.0753d0*lzs))
113  msp(40) = min(0.24d0+lzs*(0.18d0+0.595d0*lzs),0.306d0+0.053d0*lzs)
114  msp(41) = min(0.33d0+lzs*(0.132d0+0.218d0*lzs),
115  & 0.3625d0+0.062d0*lzs)
116  msp(42) = (msp(33)+msp(34)*m2**msp(36))/
117  & (m2**0.4d0+msp(35)*m2**1.9d0)
118 * Lbeta
119  msp(43) = xl(44)+lzs*(xl(45)+lzs*(xl(46)+lzs*(xl(47)+lzs*xl(48))))
120  msp(44) = xl(49)+lzs*(xl(50)+lzs*(xl(51)+lzs*(xl(52)+lzs*xl(53))))
121  msp(45) = xl(54)+lzs*(xl(55)+lzs*xl(56))
122  msp(46) = min(1.4d0,1.5135d0+0.3769d0*lzs)
123  msp(46) = max(0.6355d0-0.4192d0*lzs,max(1.25d0,msp(46)))
124 * Lhook
125  msp(47) = xl(57)+lzs*(xl(58)+lzs*(xl(59)+lzs*xl(60)))
126  msp(48) = xl(61)+lzs*(xl(62)+lzs*(xl(63)+lzs*xl(64)))
127  msp(49) = xl(65)+lzs*(xl(66)+lzs*(xl(67)+lzs*xl(68)))
128  msp(50) = xl(69)+lzs*(xl(70)+lzs*(xl(71)+lzs*xl(72)))
129  msp(51) = min(1.4d0,1.5135d0+0.3769d0*lzs)
130  msp(51) = max(0.6355d0-0.4192d0*lzs,max(1.25d0,msp(51)))
131 * Rtms
132  msp(52) = xr(1)+lzs*(xr(2)+lzs*(xr(3)+lzs*(xr(4)+lzs*xr(5))))
133  msp(53) = xr(6)+lzs*(xr(7)+lzs*(xr(8)+lzs*(xr(9)+lzs*xr(10))))
134  msp(54) = xr(11)+lzs*(xr(12)+lzs*(xr(13)+lzs*(xr(14)+lzs*xr(15))))
135  msp(55) = xr(16)+lzs*(xr(17)+lzs*(xr(18)+lzs*xr(19)))
136  msp(56) = xr(20)+lzs*(xr(21)+lzs*(xr(22)+lzs*xr(23)))
137  msp(52) = msp(52)*msp(54)
138  msp(53) = msp(53)*msp(54)
139  msp(57) = xr(24)
140  msp(58) = xr(25)+lzs*(xr(26)+lzs*(xr(27)+lzs*xr(28)))
141  msp(59) = xr(29)+lzs*(xr(30)+lzs*(xr(31)+lzs*xr(32)))
142  msp(60) = xr(33)+lzs*(xr(34)+lzs*(xr(35)+lzs*xr(36)))
143  msp(61) = xr(37)+lzs*(xr(38)+lzs*(xr(39)+lzs*xr(40)))
144 *
145  msp(62) = max(0.097d0-0.1072d0*(lz+3.d0),max(0.097d0,min(0.1461d0,
146  & 0.1461d0+0.1237d0*(lz+2.d0))))
147  msp(62) = 10.d0**msp(62)
148  m2 = msp(62) + 0.1d0
149  msp(63) = (msp(52)+msp(53)*msp(62)**msp(55))/
150  & (msp(54)+msp(62)**msp(56))
151  msp(64) = (msp(57)*m2**3+msp(58)*m2**msp(61)+
152  & msp(59)*m2**(msp(61)+1.5d0))/(msp(60)+m2**5)
153 * Ralpha
154  msp(65) = xr(41)+lzs*(xr(42)+lzs*(xr(43)+lzs*xr(44)))
155  msp(66) = xr(45)+lzs*(xr(46)+lzs*(xr(47)+lzs*xr(48)))
156  msp(67) = xr(49)+lzs*(xr(50)+lzs*(xr(51)+lzs*xr(52)))
157  msp(68) = xr(53)+lzs*(xr(54)+lzs*(xr(55)+lzs*xr(56)))
158  msp(69) = xr(57)+lzs*(xr(58)+lzs*(xr(59)+lzs*(xr(60)+lzs*xr(61))))
159  msp(70) = max(0.9d0,min(1.d0,1.116d0+0.166d0*lzs))
160  msp(71) = max(1.477d0+0.296d0*lzs,min(1.6d0,-0.308d0-1.046d0*lzs))
161  msp(71) = max(0.8d0,min(0.8d0-2.d0*lzs,msp(71)))
162  msp(72) = xr(62)+lzs*(xr(63)+lzs*xr(64))
163  msp(73) = max(0.065d0,0.0843d0-lzs*(0.0475d0+0.0352d0*lzs))
164  msp(74) = 0.0736d0+lzs*(0.0749d0+0.04426d0*lzs)
165  if(z.lt.0.004d0) msp(74) = min(0.055d0,msp(74))
166  msp(75) = max(0.091d0,min(0.121d0,0.136d0+0.0352d0*lzs))
167  msp(76) = (msp(65)*msp(71)**msp(67))/(msp(66) + msp(71)**msp(68))
168  if(msp(70).gt.msp(71))then
169  msp(70) = msp(71)
170  msp(75) = msp(76)
171  endif
172 * Rbeta
173  msp(77) = xr(65)+lzs*(xr(66)+lzs*(xr(67)+lzs*xr(68)))
174  msp(78) = xr(69)+lzs*(xr(70)+lzs*(xr(71)+lzs*xr(72)))
175  msp(79) = xr(73)+lzs*(xr(74)+lzs*(xr(75)+lzs*xr(76)))
176  msp(80) = xr(77)+lzs*(xr(78)+lzs*(xr(79)+lzs*xr(80)))
177  msp(81) = xr(81)+lzs*(xr(82)+lzs*lzs*xr(83))
178  if(z.gt.0.01d0) msp(81) = max(msp(81),0.95d0)
179  msp(82) = max(1.4d0,min(1.6d0,1.6d0+lzs*(0.764d0+0.3322d0*lzs)))
180 * Rgamma
181  msp(83) = max(xr(84)+lzs*(xr(85)+lzs*(xr(86)+lzs*xr(87))),
182  & xr(96)+lzs*(xr(97)+lzs*xr(98)))
183  msp(84) = min(0.d0,xr(88)+lzs*(xr(89)+lzs*(xr(90)+lzs*xr(91))))
184  msp(84) = max(msp(84),xr(99)+lzs*(xr(100)+lzs*xr(101)))
185  msp(85) = xr(92)+lzs*(xr(93)+lzs*(xr(94)+lzs*xr(95)))
186  msp(85) = max(0.d0,min(msp(85),7.454d0+9.046d0*lzs))
187  msp(86) = min(xr(102)+lzs*xr(103),max(2.d0,-13.3d0-18.6d0*lzs))
188  msp(87) = min(1.5d0,max(0.4d0,2.493d0+1.1475d0*lzs))
189  msp(88) = max(1.d0,min(1.27d0,0.8109d0-0.6282d0*lzs))
190  msp(88) = max(msp(88),0.6355d0-0.4192d0*lzs)
191  msp(89) = max(5.855420d-02,-0.2711d0-lzs*(0.5756d0+0.0838d0*lzs))
192 * Rhook
193  msp(90) = xr(104)+lzs*(xr(105)+lzs*(xr(106)+lzs*xr(107)))
194  msp(91) = xr(108)+lzs*(xr(109)+lzs*(xr(110)+lzs*xr(111)))
195  msp(92) = xr(112)+lzs*(xr(113)+lzs*(xr(114)+lzs*xr(115)))
196  msp(93) = xr(116)+lzs*(xr(117)+lzs*(xr(118)+lzs*xr(119)))
197  msp(94) = min(1.25d0,
198  & max(1.1d0,1.9848d0+lzs*(1.1386d0+0.3564d0*lzs)))
199  msp(95) = 0.063d0 + lzs*(0.0481d0 + 0.00984d0*lzs)
200  msp(96) = min(1.3d0,max(0.45d0,1.2d0+2.45d0*lzs))
201 * Lneta
202  if(z.gt.0.0009d0)then
203  msp(97) = 10.d0
204  else
205  msp(97) = 20.d0
206  endif
207 * Lbgb
208  gbp(1) = xg(1)+lzs*(xg(2)+lzs*(xg(3)+lzs*xg(4)))
209  gbp(2) = xg(5)+lzs*(xg(6)+lzs*(xg(7)+lzs*xg(8)))
210  gbp(3) = xg(9)+lzs*(xg(10)+lzs*(xg(11)+lzs*xg(12)))
211  gbp(4) = xg(13)+lzs*(xg(14)+lzs*(xg(15)+lzs*xg(16)))
212  gbp(5) = xg(17)+lzs*(xg(18)+lzs*xg(19))
213  gbp(6) = xg(20)+lzs*(xg(21)+lzs*xg(22))
214  gbp(3) = gbp(3)**gbp(6)
215  gbp(7) = xg(23)
216  gbp(8) = xg(24)
217 * Lbagb
218 * set gbp(16) = 1.d0 until it is reset later with an initial
219 * call to Lbagbf using mass = zpars(2) and mhefl = 0.0
220  gbp(9) = xg(25) + lzs*(xg(26) + lzs*xg(27))
221  gbp(10) = xg(28) + lzs*(xg(29) + lzs*xg(30))
222  gbp(11) = 15.d0
223  gbp(12) = xg(31)+lzs*(xg(32)+lzs*(xg(33)+lzs*xg(34)))
224  gbp(13) = xg(35)+lzs*(xg(36)+lzs*(xg(37)+lzs*xg(38)))
225  gbp(14) = xg(39)+lzs*(xg(40)+lzs*(xg(41)+lzs*xg(42)))
226  gbp(15) = xg(43)+lzs*xg(44)
227  gbp(12) = gbp(12)**gbp(15)
228  gbp(14) = gbp(14)**gbp(15)
229  gbp(16) = 1.d0
230 * Rgb
231  gbp(17) = -4.6739d0-0.9394d0*lz
232  gbp(17) = 10.d0**gbp(17)
233  gbp(17) = max(gbp(17),-0.04167d0+55.67d0*z)
234  gbp(17) = min(gbp(17),0.4771d0-9329.21d0*z**2.94d0)
235  gbp(18) = min(0.54d0,0.397d0+lzs*(0.28826d0+0.5293d0*lzs))
236  gbp(19) = max(-0.1451d0,-2.2794d0-lz*(1.5175d0+0.254d0*lz))
237  gbp(19) = 10.d0**gbp(19)
238  if(z.gt.0.004d0)then
239  gbp(19) = max(gbp(19),0.7307d0+14265.1d0*z**3.395d0)
240  endif
241  gbp(20) = xg(45)+lzs*(xg(46)+lzs*(xg(47)+lzs*(xg(48)+
242  & lzs*(xg(49)+lzs*xg(50)))))
243  gbp(21) = xg(51)+lzs*(xg(52)+lzs*(xg(53)+lzs*(xg(54)+lzs*xg(55))))
244  gbp(22) = xg(56)+lzs*(xg(57)+lzs*(xg(58)+lzs*(xg(59)+
245  & lzs*(xg(60)+lzs*xg(61)))))
246  gbp(23) = xg(62)+lzs*(xg(63)+lzs*(xg(64)+lzs*(xg(65)+lzs*xg(66))))
247 * Ragb
248  gbp(24) = min(0.99164d0-743.123d0*z**2.83d0,
249  & 1.0422d0+lzs*(0.13156d0+0.045d0*lzs))
250  gbp(25) = xg(67)+lzs*(xg(68)+lzs*(xg(69)+lzs*(xg(70)+
251  & lzs*(xg(71)+lzs*xg(72)))))
252  gbp(26) = xg(73)+lzs*(xg(74)+lzs*(xg(75)+lzs*(xg(76)+lzs*xg(77))))
253  gbp(27) = xg(78)+lzs*(xg(79)+lzs*(xg(80)+lzs*(xg(81)+
254  & lzs*(xg(82)+lzs*xg(83)))))
255  gbp(28) = xg(84)+lzs*(xg(85)+lzs*(xg(86)+lzs*(xg(87)+lzs*xg(88))))
256  gbp(29) = xg(89)+lzs*(xg(90)+lzs*(xg(91)+lzs*(xg(92)+
257  & lzs*(xg(93)+lzs*xg(94)))))
258  gbp(30) = xg(95)+lzs*(xg(96)+lzs*(xg(97)+lzs*(xg(98)+
259  & lzs*(xg(99)+lzs*xg(100)))))
260  m1 = zpars(2) - 0.2d0
261  gbp(31) = gbp(29) + gbp(30)*m1
262  gbp(32) = min(gbp(25)/zpars(2)**gbp(26),gbp(27)/zpars(2)**gbp(28))
263 * Mchei
264  gbp(33) = xg(101)**4
265  gbp(34) = xg(102)*4.d0
266 * Mcagb
267  gbp(35) = xg(103)+lzs*(xg(104)+lzs*(xg(105)+lzs*xg(106)))
268  gbp(36) = xg(107)+lzs*(xg(108)+lzs*(xg(109)+lzs*xg(110)))
269  gbp(37) = xg(111)+lzs*xg(112)
270  gbp(35) = gbp(35)**4
271  gbp(36) = gbp(36)*4.d0
272  gbp(37) = gbp(37)**4
273 * Lhei
274 * set gbp(41) = -1.d0 until it is reset later with an initial
275 * call to Lheif using mass = zpars(2) and mhefl = 0.0
276  gbp(38) = xh(1)+lzs*xh(2)
277  gbp(39) = xh(3)+lzs*xh(4)
278  gbp(40) = xh(5)
279  gbp(41) = -1.d0
280  gbp(42) = xh(6)+lzs*(xh(7)+lzs*xh(8))
281  gbp(43) = xh(9)+lzs*(xh(10)+lzs*xh(11))
282  gbp(44) = xh(12)+lzs*(xh(13)+lzs*xh(14))
283  gbp(42) = gbp(42)**2
284  gbp(44) = gbp(44)**2
285 * Lhe
286  gbp(45) = xh(15)+lzs*(xh(16)+lzs*xh(17))
287  if(lzs.gt.-1.d0)then
288  gbp(46) = 1.d0 - xh(19)*(lzs+1.d0)**xh(18)
289  else
290  gbp(46) = 1.d0
291  endif
292  gbp(47) = xh(20)+lzs*(xh(21)+lzs*xh(22))
293  gbp(48) = xh(23)+lzs*(xh(24)+lzs*xh(25))
294  gbp(45) = gbp(45)**gbp(48)
295  gbp(47) = gbp(47)**gbp(48)
296  gbp(46) = gbp(46)/zpars(3)**0.1d0+(gbp(46)*gbp(47)-gbp(45))/
297  & zpars(3)**(gbp(48)+0.1d0)
298 * Rmin
299  gbp(49) = xh(26)+lzs*(xh(27)+lzs*(xh(28)+lzs*xh(29)))
300  gbp(50) = xh(30)+lzs*(xh(31)+lzs*(xh(32)+lzs*xh(33)))
301  gbp(51) = xh(34)+lzs*(xh(35)+lzs*(xh(36)+lzs*xh(37)))
302  gbp(52) = 5.d0+xh(38)*z**xh(39)
303  gbp(53) = xh(40)+lzs*(xh(41)+lzs*(xh(42)+lzs*xh(43)))
304  gbp(49) = gbp(49)**gbp(53)
305  gbp(51) = gbp(51)**(2.d0*gbp(53))
306 * The
307 * set gbp(57) = -1.d0 until it is reset later with an initial
308 * call to Thef using mass = zpars(2), mc = 0.0 and mhefl = 0.0
309  gbp(54) = xh(44)+lzs*(xh(45)+lzs*(xh(46)+lzs*xh(47)))
310  gbp(55) = xh(48)+lzs*(xh(49)+lzs*xh(50))
311  gbp(55) = max(gbp(55),1.d0)
312  gbp(56) = xh(51)
313  gbp(57) = -1.d0
314  gbp(58) = xh(52)+lzs*(xh(53)+lzs*(xh(54)+lzs*xh(55)))
315  gbp(59) = xh(56)+lzs*(xh(57)+lzs*(xh(58)+lzs*xh(59)))
316  gbp(60) = xh(60)+lzs*(xh(61)+lzs*(xh(62)+lzs*xh(63)))
317  gbp(61) = xh(64)+lzs*xh(65)
318  gbp(58) = gbp(58)**gbp(61)
319  gbp(60) = gbp(60)**5
320 * Tbl
321  dum1 = zpars(2)/zpars(3)
322  gbp(62) = xh(66)+lzs*xh(67)
323  gbp(62) = -gbp(62)*log10(dum1)
324  gbp(63) = xh(68)
325  if(lzd.gt.0.d0) then
326  gbp(64) = 1.d0-lzd*(xh(69)+lzd*(xh(70)+lzd*xh(71)))
327  else
328  gbp(64) = 1.d0
329  end if
330  gbp(65) = 1.d0-gbp(64)*dum1**gbp(63)
331  gbp(66) = 1.d0 - lzd*(xh(77) + lzd*(xh(78) + lzd*xh(79)))
332  gbp(67) = xh(72) + lzs*(xh(73) + lzs*(xh(74) + lzs*xh(75)))
333  gbp(68) = xh(76)
334 * Lzahb
335  gbp(69) = xh(80) + lzs*(xh(81) + lzs*xh(82))
336  gbp(70) = xh(83) + lzs*(xh(84) + lzs*xh(85))
337  gbp(71) = 15.d0
338  gbp(72) = xh(86)
339  gbp(73) = xh(87)
340 * Rzahb
341  gbp(75) = xh(88) + lzs*(xh(89) + lzs*(xh(90) + lzs*xh(91)))
342  gbp(76) = xh(92) + lzs*(xh(93) + lzs*(xh(94) + lzs*xh(95)))
343  gbp(77) = xh(96) + lzs*(xh(97) + lzs*(xh(98) + lzs*xh(99)))
344 ***
345 * finish Lbagb
346  mhefl = 0.d0
347  lx = lbagbf(zpars(2),mhefl)
348  gbp(16) = lx
349 * finish LHeI
350  dum1 = 0.d0
351  lhefl = lheif(zpars(2),mhefl)
352  gbp(41) = (gbp(38)*zpars(2)**gbp(39)-lhefl)/
353  & (exp(zpars(2)*gbp(40))*lhefl)
354 * finish THe
355  thefl = thef(zpars(2),dum1,mhefl)*tbgbf(zpars(2))
356  gbp(57) = (thefl-gbp(54))/(gbp(54)*exp(gbp(56)*zpars(2)))
357 * finish Tblf
358  rb = ragbf(zpars(3),lheif(zpars(3),zpars(2)),mhefl)
359  rr = 1.d0 - rminf(zpars(3))/rb
360  rr = max(rr,1.0d-12)
361  gbp(66) = gbp(66)/(zpars(3)**gbp(67)*rr**gbp(68))
362 * finish Lzahb
363  gbp(74) = lhefl*lhef(zpars(2))
364 ***
365  kw = 0
366  tm = 0.d0
367  tn = 0.d0
368  CALL star(kw,zpars(2),zpars(2),tm,tn,tscls,lums,gb,zpars)
369  zpars(9) = mcgbf(lums(3),gb,lums(6))
370  zpars(10) = mcgbf(lums(4),gb,lums(6))
371 * set the hydrogen and helium abundances
372  zpars(11) = 0.76d0 - 3.d0*z
373  zpars(12) = 0.24d0 + 2.d0*z
374 * set constant for low-mass CHeB stars
375  zpars(13) = rminf(zpars(2))/
376  & rgbf(zpars(2),lzahbf(zpars(2),zpars(9),zpars(2)))
377 *
378  zpars(14) = z**0.4d0
379 *
380  return
381  end
382 ***