subroutine OBCS(fll,Teffl,sstop,rl,pl,tl, 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) implicit double precision (a-h,o-z) c c Checks to see if the input values of Fll, Te are within the current c OBC interpolation triangle. If not, it grinds away to find a triangle c that does surround the point. It then computes the conditions at the c base of the envelopes of the triangle, and computes the numerical c logarithmic derivatives of R, P, and T w.r.t. log Te and log L. c Finally, computes the value of logR, logP, and logT for the specified c log(L) and log(Te). c c The triangle has three corners at (cl(i),ct(i)) in the H-R diagram. c At the base of each of the three envelopes, the quanities logP, logT, c and log R are stored in bpl(i), btl(i), and brl(i) respectively c c Common block CTRL contains most of the control parameters for ISUEVO c its used to get dll and dlt common/CTRL/dll,dlt,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot common/KFIRST/newrun c if first step of run, compute three new corner positions if (newrun.ne.2) then cl(1)=fll-0.2d0/2.2415927d0 ct(1)=teffl-0.2d0/2.2415927d0 cl(2)=cl(1) ct(2)=ct(1)+dlt cl(3)=cl(1)+dll ct(3)=ct(1) newrun=2 krot=1 endif c icorn1=0 icorn2=0 icorn3=0 c call SETOBC(ICORN1,ICORN2,ICORN3,teffl,fll) c if (icorn1.eq.1) then call ENVELOP(cl(1),ct(1),sstop,rbotl,pbotl,tbotl) write (4,102) 1,rbotl,pbotl,tbotl,ct(1),cl(1) 102 format(' Env #',i2,'; Rbot= ',f7.4,' Pbot= ',f7.4, 1 ' Tbot=',f6.4,' Te=',f6.4,' L= ', f7.4) brl(1)=rbotl bpl(1)=pbotl btl(1)=tbotl endif if (icorn2.eq.1) then call ENVELOP(cl(2),ct(2),sstop,rbotl,pbotl,tbotl) write (4,102) 2,rbotl,pbotl,tbotl,ct(2),cl(2) brl(2)=rbotl bpl(2)=pbotl btl(2)=tbotl endif if (icorn3.eq.1) then call ENVELOP(cl(3),ct(3),sstop,rbotl,pbotl,tbotl) write (4,102) 3,rbotl,pbotl,tbotl,ct(3),cl(3) brl(3)=rbotl bpl(3)=pbotl btl(3)=tbotl endif c Care taken to ensure that the correct differences are taken w.r.t. c triangle parity and orientation. c logarithmic Derivatives with respect to log(L) at constant temperature if (dabs(cl(3)-cl(1)).gt.1d-8) then dlRlL=(brl(3)-brl(1))/(cl(3)-cl(1)) dlPlL=(bpl(3)-bpl(1))/(cl(3)-cl(1)) dlTlL=(btl(3)-btl(1))/(cl(3)-cl(1)) else dlRlL=(brl(2)-brl(1))/(cl(2)-cl(1)) dlPlL=(bpl(2)-bpl(1))/(cl(2)-cl(1)) dlTlL=(btl(2)-btl(1))/(cl(2)-cl(1)) endif c logarithmic Derivatives with respect to log(Te) at consant L if (dabs(ct(3)-ct(1)).gt.1d-8) then dlRlTe=(brl(3)-brl(1))/(ct(3)-ct(1)) dlPlTe=(bpl(3)-bpl(1))/(ct(3)-ct(1)) dlTlTe=(btl(3)-btl(1))/(ct(3)-ct(1)) else dlRlTe=(brl(2)-brl(1))/(ct(2)-ct(1)) dlPlTe=(bpl(2)-bpl(1))/(ct(2)-ct(1)) dlTlTe=(btl(2)-btl(1))/(ct(2)-ct(1)) endif c Interpolation calculation of quantities at base at the point dll=fll-cl(1) dlte=teffl-ct(1) rl=brl(1)+dll*dlRlL+dlte*dlRlTe pl=bpl(1)+dll*dlPlL+dlte*dlPlTe tl=btl(1)+dll*dlTlL+dlte*dlTlTe if (icorn1.eq.1.or.icorn2.eq.1.or.icorn3.eq.1) then write (4,104) dlRlL,dlPlL,dlTlL,dlRlTe,dlPlTe,dlTlTe 104 format (" Derivs: L:(",2(f8.4,","),f8.4,")", 1 " Te: (",2(f8.4,","),f8.4,")") endif c return end c subroutine ENVELOP(fllin,tefflin,sstop,rbotl,pbotl,tbotl) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Envelope integration down to specified mass sstop. Integrates an c envelope at log(L/Lo)=fllin, log(Te)=tefflin c 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 common/TORHS/cmp(15),sp,st,ishell dimension y1(4) COMMON/PATH/KMAX,KOUNT,DXSAV,XP(2000),YP(50,2000) external RHS,rkqs kmax=2000 dxsav=0.d0 c fl=10.d0**fllin teff=10.d0**tefflin x=comp(nj,1)+comp(nj,2) y=comp(nj,4)+comp(nj,3) c12=comp(nj,8) o16=comp(nj,12) c c First the atmospheric integration to obtain the pressure and mass c at the photosphere. call ATMINT(fl,teff,alfa,x,y,z,c12,o16,fms,fmenv,Ro,peff) reff=dsqrt(cls*fl/cp4/csb/teff**4)/crs y1(1)=dlog(reff) y1(2)=dlog(peff) y1(3)=dlog(Teff) y1(4)=fl do 10 i=1,15 cmp(i)=comp(nj,i) 10 continue c c Now integrate downwards to interior mass point with Sr=Sstop c using the RHS subroutine for the derivatives of the physical c quantities. c itmp=ishell ishell=0 sfirst=fmenv/10.d0 eps=3.d-6 call ODEINT(y1,4,fmenv,sstop,eps,sfirst,0.d0,nok,nbad, 1 RHS,rkqs) ishell=itmp rbotl=y1(1)/dlog(10.d0) pbotl=y1(2)/dlog(10.d0) tbotl=y1(3)/dlog(10.d0) c return end c subroutine SETOBC(Icorn1,Icorn2,Icorn3,teffl,fll) implicit double precision (a-h,o-z) c c Determine the new points in the OBC grid within the H-R diagram, and c set the Icorn flag to signal the integrator. Uses values of g1,g2,g3 c computed by IZITOUT to flip triangle once to try and enclose the point c in the triangle. Does NOT reset the Icorns so that multiple calls can c result in multiple resets of corner positions without forgetting what was c done before. c common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot c c Calculate new postions in L,Te space for corners of triangle c 100 continue g1=krot*((cl(2)-cl(3))*(Teffl-ct(2)) 1 +(ct(3)-ct(2))*(fll-cl(2))) g2=krot*((cl(3)-cl(1))*(Teffl-ct(3)) 2 +(ct(1)-ct(3))*(fll-cl(3))) g3=krot*((cl(1)-cl(2))*(Teffl-ct(1)) 3 +(ct(2)-ct(1))*(fll-cl(1))) c G1 is less than zero: Compute new position for point 1, change c sense of rotation, and set flag for new calculation of point 1 if (g1.lt.0.d0) then cl(1)=cl(3)+cl(2)-cl(1) ct(1)=ct(3)+ct(2)-ct(1) icorn1=1 krot=-krot go to 100 endif c G2 is less than zero: Compute new position for point 2, change c sense of rotation, and set flag for new calculation of point 2 if (g2.lt.0.d0) then cl(2)=2.d0*cl(1)-cl(2) ct(2)=2.d0*ct(1)-ct(2) icorn2=1 krot=-krot go to 100 endif c G3 is less than zero: Compute new position for point 3, change c sense of rotation, and set flag for new calculation of point 3 if (g3.lt.0.d0) then cl(3)=2.d0*cl(1)-cl(3) ct(3)=2.d0*ct(1)-ct(3) icorn3=1 krot=-krot go to 100 endif c return end c subroutine PTLTE(pl,tl,teobcl,flobcl) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Subroutine to compute the value of log(Te) for a model c from the values of ln(P),ln(T) at the last interior zone c 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 KIPP contains info about the OBC envelope triangle common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot 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 C Compute derivatives of P,T w.r.t. L,Te. Make sure that the correct c triangle corners are used in derivatives! c logarithmic Derivatives with respect to log(L) at constant temperature c if (dabs(cl(3)-cl(1)).gt.1d-8) then dlRlL=(brl(3)-brl(1))/(cl(3)-cl(1)) dlPlL=(bpl(3)-bpl(1))/(cl(3)-cl(1)) dlTlL=(btl(3)-btl(1))/(cl(3)-cl(1)) else dlRlL=(brl(2)-brl(1))/(cl(2)-cl(1)) dlPlL=(bpl(2)-bpl(1))/(cl(2)-cl(1)) dlTlL=(btl(2)-btl(1))/(cl(2)-cl(1)) endif c logarithmic Derivatives with respect to log(Te) at consant L if (dabs(ct(3)-ct(1)).gt.1d-8) then dlRlTe=(brl(3)-brl(1))/(ct(3)-ct(1)) dlPlTe=(bpl(3)-bpl(1))/(ct(3)-ct(1)) dlTlTe=(btl(3)-btl(1))/(ct(3)-ct(1)) else dlRlTe=(brl(2)-brl(1))/(ct(2)-ct(1)) dlPlTe=(bpl(2)-bpl(1))/(ct(2)-ct(1)) dlTlTe=(btl(2)-btl(1))/(ct(2)-ct(1)) endif c Derivatve of Te w.r.t. P,T at base dlTelP=1.d0/(dlPlTe-dlPlL*(dlTlTe/dlTlL)) dlTelT=1.d0/(dlTlTe-dlTlL*(dlPlTe/dlPlL)) c Derivatve of L w.r.t. P,T at base dlLlP=1.d0/(dlPlL-dlPlTe*(dlTlL/dlTlTe)) dlLlT=1.d0/(dlTlL-dlTlTe*(dlPlL/dlPlTe)) c Effective temperature teobcl=ct(1)+(tl-btl(1))*dlTelT+(pl-bpl(1))*dlTelP c Luminosity flobcl=cl(1)+(tl-btl(1))*dlLlT+(pl-bpl(1))*dlLlP c return end