subroutine difeq(k,k1,k2,jsf,is1,isf,indexv, 1 ne,s,nsi,nsj,y,nyj,nyk) implicit double precision (a-h,o-z) c c Subroutine to supply coefficients for SOLVDE for the c differential equations of stellar structure equations c k is the point at which the difference equations c are to be evaluated c k1 and k2 label the first and last points c in the mesh. c is1 : starting row (equation) that needs to be filled c isf : final row (equation) that needs to be filled c jsf : column containing the r.h.s. of the difference c equations; c y : y(j,k)value of variable j at meshpoint k c s : s(i,j), where i labels the equation, and j is the c index of the variable. j=1,n refers to the values c at meshpoint k-1, j=n+1,2n refers to the values c at meshpoint k. c c Note: variables postscripted with a 0 are for zone k-1 c variables postscripted with a 1 are for zone k c c Arrays for difference equation solution: parameter(jmax=1999) dimension y(nyj,nyk) dimension indexv(nyj),s(nsi,nsj) c Common block CTRL contains most of the control parameters for passing common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block BETA contains offsets for difference equations common/BETA/br,bp,bt,bl c Common block TORHS contains things for use in computing DYDX c things exported are compostion variables cmp, sp and st are changes c in P and T from last model to this one common/TORHS/cmp(15),sp,st,ishell c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn c dimension yy(4),dydx0(4),dydx1(4),dcdt(15) c c h(j,1)=lnR c h(j,2)=lnP c h(j,3)=lnT c h(j,4)=Fl c h(j,5)=ln(rho) c h(j,6)=delad c h(j,7)=delrad c h(j,8)=del c h(j,9)=q =-dln(rho)/dln(T)|P=Xt/Xd c h(j,10)=Xt =dln(P)/dln(T)|rho c h(j,11)=Xd =dln(P)/dln(rho)|T c h(j,12)=kappa c h(j,13)=dln(kappa)/dln(T)|rho c h(j,14)=dln(kappa)/dln(rho)|T c h(j,15)=epsilon c h(j,16)=dln(epsilon)/dln(T)|rho c h(j,17)=dln(epsilon)/dln(rho)|T c h(j,18)=mu c h(j,19)=delta(lnP) from previous model to this model c h(j,20)=delta(lnT) " " " " " " c c CENTRAL BOUNDARY CONDITIONS: really trivial, as most of the work done c for central convergence is at zone 1 c if (k.eq.k1) then call DIFCBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) c print '(i5,1p4e13.5)',1,0.d0,0.d0,s(3,jsf),s(4,jsf) return endif c c SURFACE BOUNDARY CONDITIONS c if (k.gt.k2) then call DIFOBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) c print '(i5,1p4e13.5)',k,s(1,jsf),s(2,jsf),0.d0,0.d0 return endif c c EQUATIONS FOR THE INTERIOR POINTS c sr1=sr(k) sr0=sr(k-1) c fmr is Mr fmr1=fms-sr1 fmr0=fms-sr0 r1=dexp(y(1,k)) r0=dexp(y(1,k-1)) p1=dexp(y(2,k)) p0=dexp(y(2,k-1)) t1=dexp(y(3,k)) t0=dexp(y(3,k-1)) fl1=y(4,k) fl0=y(4,k-1) c careful here, if fl1 or fl0 is exactly zero then things blow up... if zero c then set them to something very small! if (fl1.eq.0.d0) fl1=1.d-99 if (fl0.eq.0.d0) fl0=1.d-99 cc Note that the H, COMP, SR arrays include the center cc as Zone 1, but the Y array for numerics doesn't. cc Hence the value of k is one greater when dealing cc with H, COMP, SR than with Y sp1=h(k,19) sp0=h(k-1,19) st1=h(k,20) st0=h(k-1,20) x1=comp(k,1)+comp(k,2) x0=comp(k-1,1)+comp(k-1,2) y1=comp(k,4)+comp(k,3) y0=comp(k-1,4)+comp(k-1,3) XC121=comp(k,8) XC120=comp(k-1,8) XO161=comp(k,12) XO160=comp(k-1,12) c c Find right hand sides of the differential equations c c load common block TORHS at zone k, then call RHS do 10 i=1,15 cmp(i)=comp(k,i) 10 continue YY(1)=y(1,k) YY(2)=y(2,k) YY(3)=y(3,k) YY(4)=y(4,k) sp=sp1 st=st1 c set shell number for RHS call so that you don't need to c call EOS within RHS ishell=k call RHS(sr1,yy,dydx1) c constant for figuring radiative gradient cgr=3.d0*cls/(16.d0*cp4*cg*csb*cms) c Now extract some useful things for zone K (subscript 1) q1=h(k,9) ga1=h(k,6) gr1=h(k,7) g1=h(k,8) d1=dexp(h(k,5)) xd1=h(k,11) xt1=h(k,10) fk1=h(k,12) dlklt1=h(k,13) dlkld1=h(k,14) eps1=h(k,15) dlelt1=h(k,16) dleld1=h(k,17) cp1=p1*q1/d1/t1/ga1 c gr1=cgr*fl1*fk1*p1/(fmr1*t1**4) c find the partials of kappa w.r.t. P at constant T, and T at constant P dlklp1=dlkld1/xd1 dlklt1=dlklt1-q1*dlkld1 c c load common block TORHS at zone k-1, then call RHS c do 11 i=1,15 cmp(i)=comp(k-1,i) 11 continue YY(1)=y(1,k-1) YY(2)=y(2,k-1) YY(3)=y(3,k-1) YY(4)=y(4,k-1) sp=sp0 st=st0 c set shell number for RHS call so that you don't need to c call EOS within RHS c ishell=k ishell=k-1 call RHS(sr0,yy,dydx0) c Now extract some useful things for zone K-1 (subscript 0) q0=h(k-1,9) ga0=h(k-1,6) gr0=h(k-1,7) g0=h(k-1,8) d0=dexp(h(k-1,5)) xd0=h(k-1,11) xt0=h(k-1,10) fk0=h(k-1,12) dlklt0=h(k-1,13) dlkld0=h(k-1,14) eps0=h(k-1,15) dlelt0=h(k-1,16) dleld0=h(k-1,17) cp0=p0*q0/d0/t0/ga0 c gr0=cgr*fl0*fk0*p0/(fmr0*t0**4) c find the partials of kappa w.r.t. P at constant T, and T at constant P dlklp0=dlkld0/xd0 dlklt0=dlklt0-q0*dlkld0 c c Pressure and temperature derivatives of Xt and Ga at zone k c also derivatives of gradient when convective... c p1del=dexp(dlog(P1)+0.001d0) c tcut=5.5d0 call EOS(p1del,t1,x1,y1,XC121,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalp1=(gad-ga1)/0.001d0/ga1 dlqlp1=(qd-q1)/0.001d0/q1 c if convective, compute d(ln del)/d(ln P) if (ga1.lt.gr1) then cpd=p1del*qd/dd/t1/gad fkd=dexp(dlog(fk1)+0.001d0*dlklp1) grd=cgr*fl1*fkd*p1del/(fmr1*t1**4) if (gad.lt.grd) then call TGRAD(fmr1,r1,t1,p1del,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglp1=(gd-g1)/0.001d0/g1 endif c t1del=dexp(dlog(T1)+0.001d0) call EOS(p1,t1del,x1,y1,XC121,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalt1=(gad-ga1)/0.001d0/ga1 dlqlt1=(qd-q1)/0.001d0/q1 c if convective, compute d(ln del)/d(ln T) if (ga1.lt.gr1) then cpd=p1*qd/dd/t1del/gad fkd=dexp(dlog(fk1)+0.001d0*dlklt1) grd=cgr*fl1*fkd*p1/(fmr1*t1del**4) if (gad.lt.grd) then call TGRAD(fmr1,r1,t1del,p1,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglt1=(gd-g1)/0.001d0/g1 c also fix up dlgfl1 (would be 0 if adiabatic, 1 if radiative) fac=(g1-ga1)/(gr1-ga1) dlgfl1=fac c write (31,101) k, dlglp1,1.d0+dlklp1,dlglt1,-4+dlklt1,dlgfl1 c write (31,102) k, dlgalp1, dlgalt1 c write (31,101) k, ga1,gr1,g1,del1 101 format(i4,1p2e11.3,3x,2e11.3,3x,2e11.3) 102 format(i4,11x,1pe11.3,14x,e11.3) else c not convective, use derivitaves of delrad for derivatives of del dlglp1=1.d0+dlklp1 dlglt1=-4.d0+dlklt1 dlgfl1=1.d0 endif c c and at zone k-1 p0del=dexp(dlog(P0)+0.001d0) call EOS(p0del,t0,x0,y0,XC120,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalp0=(gad-ga0)/0.001d0/ga0 dlqlp0=(qd-q0)/0.001d0/q0 c if convective, compute d(ln del)/d(ln P) at zone k-1 if (ga0.lt.gr0) then c use a fresh value for the gradient cpd=p0del*qd/dd/t0/gad fkd=dexp(dlog(fk0)+0.001d0*dlklp0) grd=cgr*fl0*fkd*p0del/(fmr0*t0**4) if (gad.lt.grd) then call TGRAD(fmr0,r0,t0,p0del,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglp0=(gd-g0)/0.001d0/g0 endif c t0del=dexp(dlog(T0)+0.001d0) call EOS(p0,t0del,x0,y0,XC120,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalt0=(gad-ga0)/0.001d0/ga0 dlqlt0=(qd-q0)/0.001d0/q0 c if convective, compute d(ln del)/d(ln T) at zone k-1 if (ga0.lt.gr0) then cpd=p0*qd/dd/t0del/gad fkd=dexp(dlog(fk0)+0.001d0*dlklt0) grd=cgr*fl0*fkd*p0/(fmr0*t0del**4) if (gad.lt.grd) then call TGRAD(fmr0,r0,t0del,p0,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglt0=(gd-g0)/0.001d0/g0 c also fix up dlgfl0 (would be 0 if adiabatic, 1 if radiative) fac=(g0-ga0)/(gr0-ga0) dlgfl0=fac c write (31,101) k-1, dlglp0,1.d0+dlklp0,dlglt0,-4+dlklt0,dlgfl0 c write (31,102) k-1, dlgalp0, dlgalt0 c write (31,101) k-1, ga0,gr0,g0,del0 c write (31,101) else c not convective, use derivitaves of delrad for derivatives of del dlglp0=1.d0+dlklp0 dlglt0=-4.d0+dlklt0 dlgfl0=1.d0 endif c ================================================================== c And now for the difference equations... c the step in surface mass fraction: c ds=sr(k+1)-sr(k) ds=sr(k)-sr(k-1) c c Continuity equation (equation for lnR): c drnds1=dydx1(1) drnds0=dydx0(1) s(1,4+indexv(1))=1.d0+3.d0*ds*br*drnds1 s(1,indexv(1))=-1.d0+3.d0*ds*(1.d0-br)*drnds0 s(1,4+indexv(2))=ds*br*drnds1/Xd1 s(1,indexv(2))=ds*(1.d0-br)*drnds0/Xd0 s(1,4+indexv(3))=-ds*br*drnds1*Q1 s(1,indexv(3))=-ds*(1.d0-br)*drnds0*Q0 s(1,4+indexv(4))=0.d0 s(1,indexv(4))=0.d0 s(1,jsf)=dlog(r1)-dlog(r0)-ds*(br*drnds1+(1.d0-br)*drnds0) c c Hydrostatic equilibrium (equation for lnP) c dpnds1=dydx1(2) dpnds0=dydx0(2) s(2,4+indexv(1))=4.d0*ds*bp*dpnds1 s(2,indexv(1))=4.d0*ds*(1.d0-bp)*dpnds0 s(2,4+indexv(2))=1.d0+ds*bp*dpnds1 s(2,indexv(2))=-1.d0+ds*(1.d0-bp)*dpnds0 s(2,4+indexv(3))=0.d0 s(2,indexv(3))=0.d0 s(2,4+indexv(4))=0.d0 s(2,indexv(4))=0.d0 s(2,jsf)=dlog(p1)-dlog(p0)-ds*(bp*dpnds1+(1.d0-bp)*dpnds0) c c Thermal equilibrium (equation for lnT) c c gradients of gradients computed earlier whether convective or radiative... c dTnds1=dydx1(3) dTnds0=dydx0(3) s(3,4+indexv(1))=4.d0*ds*bt*dTnds1 s(3,indexv(1))=4.d0*ds*(1.d0-bt)*dTnds0 s(3,4+indexv(2))=ds*bt*dTnds1*(1.d0-dlglp1) s(3,indexv(2))=ds*(1.d0-bt)*dTnds0*(1.d0-dlglp0) s(3,4+indexv(3))=1.d0-ds*bt*dTnds1*dlglt1 s(3,indexv(3))=-1.d0-ds*(1.d0-bt)*dTnds0*dlglt0 s(3,4+indexv(4))=-ds*bt*dTnds1/fl1*dlgfl1 s(3,indexv(4))=-ds*(1.d0-bt)*dTnds0/fl0*dlgfl0 s(3,jsf)=dlog(T1)-dlog(T0)-ds*(bt*dTnds1+(1.d0-bt)*dTnds0) c c Energy conservation c c first some preliminaries, a-la Prather's thesis, for zone K. cl=-cms/cls sbar1=P1*q1/d1/ga1*(st1-ga1*sp1) dFlds1=dydx1(4) c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsbTn1=sbar1*(q1+dlqlt1)+q1*p1/d1/ga1* c 1 (-dlgalt1*st1) c dsbPn1=sbar1*(1.d0+dlqlp1-1.d0/Xd1)- c 1 q1*p1/d1/ga1*(st1*dlgalp1) dsbTn1=sbar1*(q1+dlqlt1)+q1*p1/d1/ga1* 1 (1.d0-dlgalt1*st1) dsbPn1=sbar1*(1.d0+dlqlp1-1.d0/Xd1)- 1 q1*p1/d1/ga1*(ga1+st1*dlgalp1) c and some preliminaries, a-la Prather's thesis, for zone K-1. sbar0=P0*q0/d0/ga0*(st0-ga0*sp0) dFlds0=dydx0(4) c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsbTn0=sbar0*(q0+dlqlt0)+q0*p0/d0/ga0* c 1 (-dlgalt0*st0) c dsbPn0=sbar0*(1.d0+dlqlp0-1.d0/Xd0)- c 1 q0*p0/d0/ga0*(st0*dlgalp0) dsbTn0=sbar0*(q0+dlqlt0)+q0*p0/d0/ga0* 1 (1.d0-dlgalt0*st0) dsbPn0=sbar0*(1.d0+dlqlp0-1.d0/Xd0)- 1 q0*p0/d0/ga0*(ga0+st0*dlgalp0) c s(4,4+indexv(1))=0.d0 s(4,indexv(1))=0.d0 s(4,4+indexv(2))=-cl*ds*bl*(eps1*(dleld1/Xd1)-dsbPn1/dtime) s(4,indexv(2))=-cl*ds*(1.d0-bl)* 1 (eps0*(dleld0/Xd0)-dsbPn0/dtime) s(4,4+indexv(3))=-cl*ds*bl*(eps1*(dlelt1-dleld1*q1) 1 -dsbTn1/dtime) s(4,indexv(3))=-cl*ds*(1.d0-bl)*(eps0*(dlelt0-dleld0*q0) 1 -dsbTn0/dtime) s(4,4+indexv(4))=1.d0 s(4,indexv(4))=-1.d0 s(4,jsf)=fl1-fl0-ds*(bl*dFlds1+(1.d0-bl)*dFlds0) c c Done with equations for the interior! c c print 100,k,(s(i,jsf),i=1,4) 100 format(i5,1p4e13.5) return end c subroutine DIFCBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) implicit double precision (a-h,o-z) parameter(jmax=1999) dimension y(nyj,nyk),indexv(nyj),s(nsi,nsj) c Common block CTRL contains most of the control parameters for passing common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn dimension cmp(15),dcdt(15) c c Central Boundary conditions c c some preliminaries: note that CBC's applied at zone 2. Thermo c quantities derived from pcl,tcl c fac=1.d0/dlog(10.d0) fmr=fms-sr(1) r=dexp(y(1,1)) fl=y(4,1) p=dexp(y(2,1)) t=dexp(y(3,1)) pc=dexp(pcl/fac) tc=dexp(tcl/fac) c p=pc t=tc c sp=h(1,19) st=h(1,20) x=comp(1,1)+comp(1,2) YY=comp(1,3)+comp(1,4) XC12=comp(1,8) XO16=comp(1,12) do 10 i=1,15 cmp(i)=comp(1,i) 10 continue c tcut=5.5d0 c tcut=6.d0 call EOS(P,T,x,YY,XC12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu) call BNUKE(dlog10(d),dlog10(T),cmp,dcdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu,1,dlelt,dleld) c Derivatives of del-ad and q w.r.t. p and t Pdel=dexp(dlog(p)+0.001d0) call EOS(Pdel,T,x,YY,XC12,tcut,q1,cp1,ga1,d1,xd1,xt1,u1, 1 cv1,eta1,fmu1) dlgalp=(ga1-ga)/0.001d0/ga dlqlp=(q1-q)/0.001d0/q Tdel=dexp(dlog(t)+0.001d0) call EOS(P,Tdel,x,YY,XC12,tcut,q1,cp1,ga1,d1,xd1,xt1,u1, 1 cv1,eta1,fmu1) dlgalt=(ga1-ga)/0.001d0/ga dlqlt=(q1-q)/0.001d0/q c c Continuity boundary condition: Relates R with Rho(P,T) through Mr c ccont=3.d0/cp4*cms/crs**3 s(3,4+indexv(1))=1.d0 s(3,indexv(1))=0.d0 s(3,4+indexv(2))=1.d0/3.d0/xd s(3,indexv(2))=0.d0 s(3,4+indexv(3))=-1.d0/3.d0*q s(3,indexv(3))=0.d0 s(3,4+indexv(4))=0.d0 s(3,indexv(4))=0.d0 s(3,jsf)=dlog(r)-1.d0/3.d0*(dlog(ccont)+dlog(fmr)-dlog(d)) c c Luminosity boundary condition: c clum=cms/cls dtnds=0.d0 dpnds=0.d0 c sbar=P*q/d/ga*(st-ga*sp) sbar=0.d0 c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsblt=sbar*(q+dlqlt)+q*P/d/ga*(-dlgalt*st) c dsblp=sbar*(1.d0+dlqlp-1.d0/Xd)-q*P/d*(st/ga*dlgalp) c* dsblt=sbar*(q+dlqlt)+q*P/d/ga*(1.d0-dlgalt*st) c* dsblp=sbar*(1.d0+dlqlp-1.d0/Xd)-q*P/d*(1.d0+st/ga*dlgalp) dsblt=0.d0 dsblp=0.d0 s(4,4+indexv(1))=0.d0 s(4,indexv(1))=0.d0 s(4,4+indexv(2))=-clum*fmr*(eps*(dleld/xd)-dsblp/dtime) s(4,indexv(2))=0.d0 s(4,4+indexv(3))=-clum*fmr*(eps*(dlelt-dleld*q)-dsblt/dtime) s(4,indexv(3))=0.d0 s(4,4+indexv(4))=1.d0 s(4,indexv(4))=0.d0 s(4,jsf)=fl-clum*fmr*(eps-sbar/dtime) return end c subroutine DIFOBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) implicit double precision (a-h,o-z) parameter(jmax=1999) dimension indexv(nyj),s(nsi,nsj),cmp(15),y(nyj,nyk) c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR c Common block KIPP contains positions for the corners of the OBC triangle COMMON/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot c c Surface boundary conditions for relaxation scheme c c rn=y(1,nj-1) c pn=y(2,nj-1) c tn=y(3,nj-1) c fl=y(4,nj-1) rn=y(1,nj) pn=y(2,nj) tn=y(3,nj) fl=y(4,nj) fac=1.d0/dlog(10.d0) c call OBCS(dlog10(fl),teffl,sr(nj),rl,pl,tl, c 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) call OBCS(dlog10(fls),teffl,sr(nj),rl,pl,tl, 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) c Derivatves of Te,L w.r.t. P,T at base dlTelP=1.d0/(dlPlTe-dlPlL*(dlTlTe/dlTlL)) dlTelT=1.d0/(dlTlTe-dlTlL*(dlPlTe/dlPlL)) dlLlP=1.d0/(dlPlL-dlPlTe*(dlTlL/dlTlTe)) dlLlT=1.d0/(dlTlL-dlTlTe*(dlPlL/dlPlTe)) c Derivative of R w.r.t. P,T at base dlRlP=dlRlTe*dlTelP+dlRlL*dlLlP dlRlT=dlRlTe*dlTelT+dlRlL*dlLlT c c Radius boundary condition c s(1,4+indexv(1))=fac s(1,indexv(1))=0.d0 s(1,4+indexv(2))=-fac*dlRlP s(1,indexv(2))=0.d0 s(1,4+indexv(3))=-fac*dlRlT s(1,indexv(3))=0.d0 s(1,4+indexv(4))=0.d0 s(1,indexv(4))=0.d0 rbasel=brl(1)+(pn*fac-bpl(1))*dlRlP+(tn*fac-btl(1))*dlRlT s(1,jsf)=rn*fac-rbasel c c Luminosity boundary condition c s(2,4+indexv(1))=0.d0 s(2,indexv(1))=0.d0 s(2,4+indexv(2))=-fac*dlLlP s(2,indexv(2))=0.d0 s(2,4+indexv(3))=-fac*dlLlT s(2,indexv(3))=0.d0 s(2,4+indexv(4))=fac/fl s(2,indexv(4))=0.d0 fll=cl(1)+(pn*fac-bpl(1))*dlLlP+(tn*fac-btl(1))*dlLlT s(2,jsf)=dlog10(fl)-fll c return end