subroutine instruct c-----VERSION of November 20, 1995----------------------------------------- c----------------------------------------------------------------------- c c This subroutine contains instructions for using the subroutine c OPAC( z, xh, xxc, xxo, t6, r ). c The purpose of the subroutine OPAC is to perform 4 or 5 variable c interpolation on log10(kappa). The opacity tables to be interpolated c are known to have somewhat random numerical errors of a few percent. c Consequently, adjusting the data prior to performing the interpolation c is justified at this level. The code first reads the original(unsmoothed) c tabular data, this data is then passed through a smoothing filter; using c a set of routines developed by Mike Seaton (see M.J. Seaton,MNRAS 265, c L25,1993). It is the adjusted data that is actually used in carrying out the c interpolations in OPAC. Steps are taken, as described below to insure c that the interpolated data is also smooth. The initial adjustment step c helps improve the smoothness of the OPAC output,particularly at the c smallest values of R. The medium to large R output is only slightly c effected by this step. It takes only a few seconds to carryout the initial c data smoothing step. This step can be skipped by setting the parameter c ismdata=1 in subroutine readco. c The interpolation variables are : c xh The hydrogen mass fraction, X c xxc The enhanced carbon mass fraction, delta Xc. c The total C mass fraction, Xc, is the sum of the initial c amount included in the metal mass fraction, Z, and c delta Xc c xxo The enhanced oxygen mass fraction, delta Xo c t6 The temperature in millions of degrees Kelvin, T6 c r =rho(gm/cc)/T6**3, R c c Additional input to OPAC is: c c z The metallicity, Z c c An interpolation between overlapping quadratics is used to obtain c smoothed results. A 4x4 grid in logT6 and logR is used to interpolate c in four different 3x3 sub-grids. Linear interpolation between quadratic c fits in these different sub-grids gives smoothed results in both log T6 c and Log R directions. This procedure produces results that are similar c to bicubic spline interpolation, but require storage of only local c information. c Each of the individual tables in a file Gx**x**z covers 70 temperatures c in the range logT=3.75[T6=0.0056341325]( referred to as temperature 1) to c logT=8.7[T6=501.187] and 19 values of log R in the range -8 (referred to as 1) c to +1; at half-integer steps. (NOTE: earlier tables were explicitly in c terms of T6. For convenience the present tables tabulate log Kappa vs logT. The c interpolation however still uses T6 for the temperature variable) c For specialized problems, if storage space is a problem, a reduced set of c data can be input . This requires a recompilation with altered parameter c values. In order to limit the range of R, set the parameter nrb= index of c lowest value of log R to use(count from log R=-8). Then set the parameter c nre to the index of the largest value of log R to use. (NOTE: nre-nrb must c be at least 5). To ignore the first ntb-1 values of T6 (starting from c temperature 1) set the parameter ntb to desired value. No other parameters c should be modified. c A five variable interpolation is done when X is greater than 0. In this c case it is assumed that variable values of X are also needed. We have provided c sets of tables for X=0, 0.03, 0.1, 0.35 and 0.7, for each of the metallicities c 0.0,0.001, 0.004,0.01,0.02,0.03,05 and .1. The five sets of tables associated with c a particular value of Z (eg. Gx0z01, Gx03x01,Gx10z01, Gx35z01, c Gx70z01) should be placed in files named codataa, codatab, codatac, codatad, c and codatae, respectively. To create storage for this data the routines must c be recompiled with the parameter mx=5. Again if storage is a problem the T6 c and log R ranges can be restricted. This version of the interpolation routine c does not interpolate in Z. Values of Z not in the table can be obtained by c interpolating the existing tables to produce similar tables for the Z of interest. c Interpolation in xh,xxo,xxc,t6, and r can be obtained as just described c A 4 variable interpolation in xxc, xxo, t6 and r is performed in the special c case X=0. The set of 60 data tables in (xxc,xxo) for a given Z that have c been provided for X=0 should be placed in a file called 'codataa'. This file c will be read from unit 2 in the subroutine readco. In this special case the c set of routines provided should be compiled with the parameter mx=1 (there are c 4 occurances). (NOTE: The version of the code described above, intended c for X> 0, also handles this case, but takes more storage space since mx=5). c If you want to work with a single value of X which is not zero, i.e., c X=0.03, 0.1, 0.35, or 0.70, then compile the code with mx=1 but change the statement c c data (xa(i), i=1,5 )/0.0, 0.03, 0.1, 0.35, 0.7/ c c If for example you want to use only the table for X=.70, then c c data (xa(i), i=1,5 )/0.70, 0.03, 0.1, 0.35, 0.0/ c c You also need to place the tables for X=0.7 into a file named c codataa. c ***CAUTION*** c As a result of the mixing procedure used to calculate the data a few c X=0.0, low T-small R, table values fell outside the range of T and R accessible c from the X=0.35 data directly calculated for this purpose. These T-R c locations are filled in with 9.99 (or for diagnostic purposes in some cases c larger values. At the same locations the derivatives are set to 99.9. When c T-R falls in this region a message is issued by the interpolation code to c inform the user of this situation. Presumable very few users will have c applications that take them into this region. c c Your routine that performs the call to OPAC should include the c statement: c c common/e/ opact,dopact,dopacr,dopactd c c These variables have the following meanings: c c OPACT Is the Log of the Rosseland mean opacity: Log(kappa) c DOPACT Is Dlog(kappa)/Dlog(T6) ! at constant R c DOPACTD Is Dlog(kappa)/Dlog(T6) ! at constant density c DOPACR Is Dlog(kappa)/Dlog(R), ! at constant T6 c dum=0.0 return end c c*********************************************************************** subroutine opac (z,xh,xxc,xxo,t6,r) c..... The purpose of this subroutine is to interpolate log kappa c and obtain smooth derivatives. c in C/O abundance and in T6, R,i.e. (xc,xo,T6,R) c xxc=Xc=carbon mass fraction c xxo=Xo=oxygen mass abundance c t6=T6=temperature in millions of degrees kelvin c r=R=density(g/cm**3)/T6**3 c..... to use opac insert common/e/ in the calling routine. c This common contains interpolated values for kappa and its c first derivities. c save integer w parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) common/aaa/ oxf(mx,mc),cxf(mx,mc),xcdf(mx,mc),xodf(mx,mc), . opl(mx,nt,nr),itime(mx),cxdf(mx,mc),oxdf(mx,mc) common/aa/ q(4), h(4), xcd(mc),xod(mc), xc(mc),xo(mo) . ,xcs(mc),xos(mo), cxd(mc),oxd(mo),cx(mc),ox(mo),zzz,xxh,xx(mx) . ,nc,no common/a/ co(mx,mc,mo,nt,nr), diag(mx,mc,nt,nr), index(101), . t6list(nt),alr(nr),n(mx,mc),alt(nt),diago(mx,mo,nt,nr),opk(nt,nr) . ,dfs(nt),dfsr(nr),a(3,mx),b(3),m,mf,xa(8),alrf(nrm),cof(ntm,nrm) .,t6listf(ntm),opk2(nt,nr),dfsx(mx) common/b/ itab(mx,ntabs),nta(nrm),x(mx,ntabs),y(mx,ntabs), . zz(mx,ntabs),xca(mx,ntabs),xoa(mx,ntabs) common/d/dkap common/bb/l1,l2,l3,l4,k1,k2,k3,k4,ip,iq,xodp,xcdp,xxco,cxx,oxx c..... OPACT- opacity obtained from a quadraric interpolation at c fixed log T6 at three values of log R; followed by quadratic c interpolation along log T6. Results smoothed by mixing c overlapping quadratics. c..... DOPACT- is Dlog(k)/Dlog(T6) smoothed by mixing quadratics. c..... DOPACR- is Dlog(k)/Dlog(R) smoothed by mixing quadratics. common/e/ opact,dopact,dopacr,dopactd common/recoin/ itimeco,mxzero c c..... INDEX refers to the index,i, of xc(i) or xo(i); xc(i) and xo(i) c are the abundance grid points. iop=1 ! provides smoothed interpolations; iop=0 gives no smoothing if(nr .lt. 6) go to 65 if((xh .gt. 1.e-6) .and. (mx .lt.4)) go to 64 c c..... set-up C/O axis points xxco=xxc+xxo if(z+xh+xxco-1.e-6 .gt. 1 ) go to 61 zzz=z+0.001 xxh=xh xxci=xxc xxoi=xxo xxi=xh t6i=t6 ri=r c c..... convert xxc and xxo to logarithmic shifted by Z cxx=log10(zzz+xxc) oxx=log10(zzz+xxo) xxx=log10(0.005+xh) slt=log10(t6) slr=log10(r) c..... set X indices ilo=2 ihi=mx 8 if(ihi-ilo .gt. 1) then imd=(ihi+ilo)/2 if(xh .le. xa(imd)+1.e-7) then ihi=imd else ilo=imd endif go to 8 endif i=ihi mf=i-2 mg=i-1 mh=i mi=i+1 mf2=mi istep1=1 if (mx .gt. 1) then istep1=mx-1 if((xh .le. xa(2)+1.e-7) .or. (xh .ge. xa(istep1)-1.e-7)) mf2=mh endif c if ((mx .eq. 1) .or. (xh .lt. 1.e-6)) then mf=1 mg=1 mh=1 mi=2 mf2=1 endif if (itime(1) .ne. 12345678) then alr(1)=-8.+(nrb-1)*0.5 do i=2,nr alr(i)=alr(i-1)+0.5 enddo alt(1)=-2.25+(ntb-1)*0.05 do i=ntb+1,46 alt(i)=alt(i-1)+.05 enddo ntd=47 if ((ntb +1) .gt. 47) ntd=ntb+1 do i=ntd,68 alt(i)=alt(i-1)+.1 enddo do i=68,70 alt(i)=alt(i-1)+.2 enddo do i=1,nt t6list(i)=10.**alt(i) enddo endif ilo=2 ihi=nr 12 if(ihi-ilo .gt. 1) then imd=(ihi+ilo)/2 if(slr .le. alr(imd)+1.e-7) then ihi=imd else ilo=imd endif go to 12 endif i=ihi l1=i-2 l2=i-1 l3=i l4=l3+1 c ilo=2 ihi=nt 11 if(ihi-ilo .gt. 1) then imd=(ihi+ilo)/2 if(t6 .le. t6list(imd)+1.e-7) then ihi=imd else ilo=imd endif go to 11 endif i=ihi k1=i-2 k2=i-1 k3=i k4=k3+1 k3s=k3+ntb-1 l3s=l3+nrb-1 c-----set-up indices to skip when low T&R data missing for X=0. kmin=0 k1in=k1 iadvance=0 mfin=mf if ((mfin .eq. 1) .and. (co(1,1,1,k1,l1) .gt. 9.)) then! data missing do i=1,6 if (co(1,1,1,i,l1) .gt. 9.) then if (xh .lt. .1) then kmin=i+1 else if (iadvance .eq. 0) then ! sfift X index to avoid X=0. iadvance=iadvance+1 mf=mf+1 mg=mg+1 mh=mh+1 mi=mi+1 mf2=mf2+1 endif endif endif enddo if ((iadvance .eq. 0) .and. (k1 .le. kmin) .and. x (slt .le. alt(kmin))) then k1=kmin if ((co(1,1,1,kmin,l1+1) .lt. 9.) .and. x ((slr+.01) .gt. alr(l1+1))) then l1=l1+1 kmin=0 k1=k1in do i=1,6 if (co(1,1,1,i,l1) .gt. 9.) kmin=i+1 enddo if ((kmin .ne. 0) .and. (k1in .lt. kmin)) k1=kmin endif endif if ((slt+.001) .lt. alt(k1)) then write (*,'("OPAL data not available for X=", f7.5," logT6=", f7.3, x " logR=",f7.3)') xh,slt,slr opact=30. dopact=99. dopacr=99. dopactd=99. return endif l2=l1+1 l3=l2+1 l4=l3+1 l3s=l3+nrb-1 k2=k1+1 k3=k2+1 k4=k3+1 k3s=k3+ntb-1 endif c-----end of check for missing data do 123 m=mf,mf2 if(mx .ge. 4) then c..... C and O fractions determined by the ray through the origin that c also passes through the point (Xc,Xo). Specific interpolation c values determined by tabulated X values;i.e. xa(m). Inter- c polation along the ray gives log (kappa(Xc,Xo)). (Advantage c of method: keeps indices within table boundaries) c Subtract Z to prevent out-of-range C+O values for small X if(1.-xh-z.gt.1.e-6)then cmod=(1.-xa(m)-z)/(1.-xh-z) else cmod=0. endif xxc=cmod*xxci xxo=cmod*xxoi cxx=log10(zzz+xxc) oxx=log10(zzz+xxo) endif c ythism=z+xa(m)+xxc+xxo do i=1,mc xhe=1.-xa(m)-z nc=i no=i xc(i)=xcs(i) xo(i)=xos(i) c if(xcs(i) .ge. xhe-1.e-6) then if(xcs(i) .gt. xhe) then xc(i)=xhe xo(i)=xhe go to 3 endif enddo 3 continue c if(itime(m) .ne. 12345678) then itime(m)=12345678 mxzero=0 do i=1,mx xx(i)=log10(0.005+xa(i)) if(xa(i) .eq. 0.0) mxzero=i enddo c..... this is the first time throught this m. Calculate the decadic c log of the perimeter points shifted by Z+0.001(to avoid divergence c at origin); m refers to xa(m); the hydrogen table value. c c note that the nc-th elements are sometimes needed! do i=1,nc oxf(m,i)=log10(zzz+xo(i)) cxf(m,i)=log10(zzz+xc(i)) xcdf(m,i)=-xo(i)+xo(no) xodf(m,i)=-xc(i)+xc(nc) cxdf(m,i)=log10(zzz+xcdf(m,i)) oxdf(m,i)=log10(zzz+xodf(m,i)) enddo c c note that the nc-th elements are sometimes needed! do i=1,nc ox(i)=oxf(m,i) cx(i)=cxf(m,i) xcd(i)=xcdf(m,i) xod(i)=xodf(m,i) cxd(i)=cxdf(m,i) oxd(i)=oxdf(m,i) enddo c c.....read the data files call readco endif c do i=1,nc ox(i)=oxf(m,i) cx(i)=cxf(m,i) xcd(i)=xcdf(m,i) xod(i)=xodf(m,i) cxd(i)=cxdf(m,i) oxd(i)=oxdf(m,i) enddo c c..... Determine log R and log T6 grid points to use in the c interpolation. if((slt .lt. alt(1)).or.(slt .gt. alt(nt))) go to 62 if((slr .lt. alr (1)).or.(slr .gt. alr(nr))) go to 62 c if (m .eq. mf) then ! calculate table indices c if((mf2 .ne. mxzero) .and. (k3s .gt. ntm)) go to 62 do i=14,18 if((l3s .gt. i) .and. (k3s .gt. nta(i+1))) go to 62 enddo ip=3 iq=3 ntlimit=nta(l3s) if((k3s .eq. ntlimit) .or. (iop .eq. 0)) then ip=2 iq=2 endif if(t6 .le. t6list(2)+1.e-7 .or. iop .eq. 0) ip=2 if((l3 .eq. nr) .or. (iop .eq. 0)) then ! right edge of full table iq=2 ip=2 endif if(slr .le. alr(2)+1.e-7 .or. iop .eq. 0) iq=2 endif c xodp=max(-xxc+xc(nc),0.) xcdp=max(-xxo+xo(no),0.) is=0 c call cointerp(xxc,xxo) 123 continue if((zz(mg,1) .ne. zz(mf,1)) .or. (zz(mh,1) .ne. zz(mf,1))) then write(*,'("Z does not match Z in codata* files you are" x ," using")') stop endif if(z .ne. zz(mf,1)) go to 66 xxc=xxci ! restores input value; necessary if stop replaced c with return xxo=xxoi ! restores input value is=0 iw=1 do 45 ir=l1,l1+iq do 46 it=k1,k1+ip if((mx .eq. 1) .or. (mf2 .eq. 1)) then opk(it,ir)=opl(mf,it,ir) go to 46 endif opk(it,ir)=quad(is,iw,xxx,opl(mf,it,ir),opl(mg,it,ir) x ,opl(mh,it,ir),xx(mf),xx(mg),xx(mh)) is=1 46 continue 45 continue if (mi .eq. mf2) then ! interpolate between quadratics is=0 iw=1 dixr=(xx(mh)-xxx)*dfsx(mh) do 47 ir=l1,l1+iq do it=k1,k1+ip opk2(it,ir)=quad(is,iw,xxx,opl(mg,it,ir),opl(mh,it,ir) x ,opl(mi,it,ir),xx(mg),xx(mh),xx(mi)) opk(it,ir)=opk(it,ir)*dixr+opk2(it,ir)*(1.-dixr) is=1 enddo 47 continue c interpolate X between two overlapping quadratics endif is=0 c c..... completed H,C,O interpolation. Now interpolate T6 and log R on a c 4x4 grid. (log(T6(i)),i=i1,i1+3),log(R(j)),j=j1,j1+3)).Procedure c mixes overlapping quadratics to obtain smoothed derivatives. c call t6rinterp(slr,slt) return c 61 write(*,'(" Mass fractions exceed unity")') xxc=xxci ! restores input value; required if stop replaced c with a return xxo=xxoi ! restores imput value stop 62 write(*,'(" T6/LogR outside of table range")') xxc=xxci ! restores input value; required if stop replaced c with a return xxo=xxoi ! restores input value stop 64 write(*,'(" X not equal to zero: To run this case it .is necessary"/ "to recompile with parameter (mx=5)")') stop 65 write(*,'("Too few R values; NRE+1-NRB < 6")') stop 66 write(*,'(" Z does not match Z in codata* files you are", . " using")') stop end c************************************************************************ subroutine cointerp(xxc,xxo) c The purpose of this subroutine is to interpolate in C and O abund- c ances. save integer w parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) common/aaa/ oxf(mx,mc),cxf(mx,mc),xcdf(mx,mc),xodf(mx,mc), . opl(mx,nt,nr),itime(mx),cxdf(mx,mc),oxdf(mx,mc) common/aa/ q(4), h(4), xcd(mc),xod(mc), xc(mc),xo(mo) . ,xcs(mc),xos(mo), cxd(mc),oxd(mo),cx(mc),ox(mo),zzz,xxh,xx(mx) . ,nc,no common/a/ co(mx,mc,mo,nt,nr), diag(mx,mc,nt,nr), index(101), . t6list(nt),alr(nr),n(mx,mc),alt(nt),diago(mx,mo,nt,nr),opk(nt,nr) . ,dfs(nt),dfsr(nr),a(3,mx),b(3),m,mf,xa(8),alrf(nrm),cof(ntm,nrm) .,t6listf(ntm),opk2(nt,nr),dfsx(mx) common/bb/l1,l2,l3,l4,k1,k2,k3,k4,ip,iq,xodp,xcdp,xxco,cxx,oxx is=0 if(xxco .lt. 1.e-6) then do ir=l1,l1+iq do it=k1,k1+ip opl(m,it,ir)=co(m,1,1,it,ir) enddo enddo is=1 go to 123 endif c include boundaries that could later cause division by 0! if(xxc .gt. xcd(3)-1.e-6) then c__________ oxdp=log10(zzz+xodp) c handle possibility that xodp=0 fac=max(min((oxx-ox(1))/max(oxdp-ox(1),1.e-6),1.),0.) do 40 ir=l1,l1+iq do 41 it=k1,k1+ip c c interpolation in region c1 c c include boundaries that could later cause division by 0! if(xxc .gt. xcd(2)-1.e-6) then iw=1 a(1,m)=quad(is,iw,cxx,co(m,nc-2,1,it,ir),co(m,nc-1,1,it,ir), . diag(m,1,it,ir),cx(nc-2),cx(nc-1),cx(nc)) iw=iw+1 a(2,m)=quad(is,iw,cxx,diag(m,1,it,ir),diag(m,2,it,ir), . diag(m,3,it,ir),cxd(1),cxd(2),cxd(3)) do w=1,2 b(w)=a(w,m) enddo c handle possibility that xodp=0 opl(m,it,ir)=b(1)+(b(2)-b(1))*fac is=1 go to 41 endif c interpolation in region c2 c iw=1 a(1,m)=quad(is,iw,cxx,co(m,nc-2,1,it,ir),co(m,nc-1,1,it,ir), . diag(m,1,it,ir),cx(nc-2),cx(nc-1),cx(nc)) iw=iw+1 a(2,m)=quad(is,iw,cxx,co(m,n(m,2)-2,2,it,ir),co(m,n(m,2)-1,2,it, . ir),diag(m,2,it,ir),cx(n(m,2)-2),cx(n(m,2)-1),cxd(2)) iw=iw+1 a(3,m)=quad(is,iw,cxx,diag(m,1,it,ir),diag(m,2,it,ir) .,diag(m,3,it,ir),cxd(1),cxd(2),cxd(3)) do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,oxx,b(1),b(2),b(3),ox(1),ox(2),oxdp) is=1 41 continue 40 continue if(is .eq. 1) go to 123 c__________ endif c c interpolation in region c3 to c6 is=0 c if(nc .ge. 5) then c__________ do 44 i=4,nc-1 c do not go beyond middle (where c3-c6 overlaps o3-o6), and if((xxc .gt. xcd(i)-1.e-6) .and. (xxo .gt. xo(i-1)-1.e-6) .and. $ (xcd(i-1) .gt. xc(i-1))) then do 42 ir=l1,l1+iq do 43 it=k1,k1+ip oxdp=log10(zzz+xodp) iw=1 m1=i-1 m2=i-2 a(1,m)=quad(is,iw,cxx,co(m,n(m,m2)-2,m2,it,ir),co(m,n(m,m2)-1,m2 x ,it,ir),diag(m,m2,it,ir),cx(n(m,m2)-2),cx(n(m,m2)-1),cxd(m2)) iw=iw+1 a(2,m)=quad(is,iw,cxx,co(m,n(m,m1)-2,m1,it,ir),co(m,n(m,m1)-1,m1 x ,it,ir),diag(m,m1,it,ir),cx(n(m,m1)-2),cx(n(m,m1)-1),cxd(m1)) iw=iw+1 a(3,m)=quad(is,iw,cxx,diag(m,m2,it,ir),diag(m,m1,it,ir), x diag(m,i,it,ir),cxd(m2),cxd(m1),cxd(i)) do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,oxx,b(1),b(2),b(3),ox(i-2),ox(i-1),oxdp) is=1 43 continue 42 continue if (is .eq. 1) go to 123 endif 44 continue c__________ endif c if (is .eq. 1) go to 123 c c include boundaries that could later cause division by 0! if(xxo .gt. xod(3)-1.e-6) then c__________ cxdp=log10(zzz+xcdp) c handle possibility that xcdp=0 fac=max(min((cxx-cx(1))/max(cxdp-cx(1),1.e-6),1.),0.) do 140 ir=l1,l1+iq do 141 it=k1,k1+ip c c interpolation in region o1 c c include boundaries that could later cause division by 0! if(xxo .gt. xod(2)-1.e-6) then iw=1 a(1,m)=quad(is,iw,oxx,co(m,1,no-2,it,ir),co(m,1,no-1,it,ir), . diago(m,no-1,it,ir),ox(no-2),ox(no-1),ox(no)) iw=iw+1 a(2,m)=quad(is,iw,oxx,diago(m,no-1,it,ir),diago(m,no-2,it,ir), .diago(m,no-3,it,ir),oxd(1),oxd(2),oxd(3)) do w=1,2 b(w)=a(w,m) enddo c handle possibility that xcdp=0 opl(m,it,ir)=b(1)+(b(2)-b(1))*fac is=1 go to 141 endif c interpolation in region o2 c iw=1 a(1,m)=quad(is,iw,oxx,co(m,1,no-2,it,ir),co(m,1,no-1,it,ir), . diago(m,no-1,it,ir),ox(no-2),ox(no-1),ox(no)) iw=iw+1 a(2,m)=quad(is,iw,oxx,co(m,2,n(m,2)-2,it,ir),co(m,2,n(m,2)-1,it, . ir),diago(m,no-2,it,ir),ox(n(m,2)-2),ox(n(m,2)-1),oxd(2)) iw=iw+1 a(3,m)=quad(is,iw,oxx,diago(m,no-1,it,ir),diago(m,no-2,it,ir), .diago(m,nc-3,it,ir),oxd(1),oxd(2),oxd(3)) do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,cxx,b(1),b(2),b(3),cx(1),cx(2),cxdp) is=1 141 continue 140 continue if(is .eq. 1) go to 123 c__________ endif c c interpolation in region o3 to o6 is=0 if(no .ge. 5) then c__________ do 144 i=4,no-1 c do not go beyond middle (where o3-o6 overlaps c3-c6), and if((xxo .gt. xod(i)-1.e-6) .and. (xxc .gt. xc(i-1)-1.e-6) .and. $ (xod(i-1) .gt. xo(i-1)-1.e-6)) then do 142 ir=l1,l1+iq do 143 it=k1,k1+ip cxdp=log10(zzz+xcdp) iw=1 m2=i-2 m1=i-1 a(1,m)=quad(is,iw,oxx,co(m,m2,n(m,m2)-2,it,ir),co(m,m2,n(m,m2)- .1,it,ir),diago(m,no-m2,it,ir),ox(n(m,m2)-2),ox(n(m,m2)-1),oxd(m2)) iw=iw+1 a(2,m)=quad(is,iw,oxx,co(m,m1,n(m,m1)-2,it,ir),co(m,m1,n(m,m1)-1, . it,ir),diago(m,no-m1,it,ir),ox(n(m,m1)-2),ox(n(m,m1)-1),oxd(m1)) iw=iw+1 a(3,m)=quad(is,iw,oxx,diago(m,no-m2,it,ir),diago(m,no-m1,it,ir), .diago(m,no-i,it,ir),oxd(m2),oxd(m1),oxd(i)) do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,cxx,b(1),b(2),b(3),cx(m2),cx(m1),cxdp) is=1 143 continue 142 continue if (is .eq. 1) go to 123 endif 144 continue c__________ endif c if (is .eq. 1) go to 123 c c.....find index of C grid. 52 ie=100*xxc+1 iei=index(ie)+1 c must also allow index = nc, to avoid extrapolation if (iei .gt. nc) iei=nc c if(iei .gt. 3) then i1=iei-2 i2=iei-1 i3=iei else i1=1 i2=2 i3=3 endif c c.....find index of O grid ie=100*xxo+1 iej=index(ie)+1 c must also allow index = no, to avoid extrapolation if (iej .gt. no) iej=no c if(iej .gt. 3) then j1=iej-2 j2=iej-1 j3=iej else j1=1 j2=2 j3=3 endif c c lower-O part of grid: interpolate C before O if(j3.lt.no .and. i3.le.n(m,j3) .and. $ (xxc.lt.xcd(j3)+1.e-6 .or. xxc.ge.xxo))then do 20 ir=l1,l1+iq do 21 it=k1,k1+ip iw=0 do jx=j1,j1+2 iw=iw+1 c if i3=n(m,jx), then must replace cx(i3) with cxd(jx) a(iw,m)=quad(is,iw,cxx,co(m,i1,jx,it,ir),co(m,i2,jx,it,ir), x co(m,i3,jx,it,ir),cx(i1),cx(i2),min(cx(i3),cxd(jx))) enddo do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,oxx,b(1),b(2),b(3),ox(j1),ox(j2),ox(j3)) is=1 21 continue 20 continue c else: high-O part of grid: must interpolate O before C else do ir=l1,l1+iq do it=k1,k1+ip iw=0 do ix=i1,i1+2 iw=iw+1 if(j3.lt.n(m,ix))then a(iw,m)=quad(is,iw,oxx,co(m,ix,j1,it,ir),co(m,ix,j2,it,ir), $ co(m,ix,j3,it,ir),ox(j1),ox(j2),ox(j3)) else a(iw,m)=quad(is,iw,oxx,co(m,ix,j1,it,ir),co(m,ix,j2,it,ir), $ diago(m,no-ix,it,ir),ox(j1),ox(j2),oxd(ix)) endif enddo do w=1,3 b(w)=a(w,m) enddo iw=iw+1 opl(m,it,ir)=quad(is,iw,cxx,b(1),b(2),b(3),cx(i1),cx(i2),cx(i3)) is=1 enddo enddo endif 123 continue return end c*********************************************************************** subroutine t6rinterp(slr,slt) c The purpose of this subroutine is to interpolate in logT6 and logR save parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) common/aaa/ oxf(mx,mc),cxf(mx,mc),xcdf(mx,mc),xodf(mx,mc), . opl(mx,nt,nr),itime(mx),cxdf(mx,mc),oxdf(mx,mc) common/aa/ q(4), h(4), xcd(mc),xod(mc), xc(mc),xo(mo) . ,xcs(mc),xos(mo), cxd(mc),oxd(mo),cx(mc),ox(mo),zzz,xxh,xx(mx) . ,nc,no common/a/ co(mx,mc,mo,nt,nr), diag(mx,mc,nt,nr), index(101), . t6list(nt),alr(nr),n(mx,mc),alt(nt),diago(mx,mo,nt,nr),opk(nt,nr) . ,dfs(nt),dfsr(nr),a(3,mx),b(3),m,mf,xa(8),alrf(nrm),cof(ntm,nrm) .,t6listf(ntm),opk2(nt,nr),dfsx(mx) common/d/dkap common/bb/l1,l2,l3,l4,k1,k2,k3,k4,ip,iq,xodp,xcdp,xxco,cxx,oxx common/e/ opact,dopact,dopacr,dopactd c is=0 iu=0 do kx=k1,k1+ip iw=1 iu=iu+1 h(iu)=quad(is,iw,slr,opk(kx,l1),opk(kx,l2),opk(kx,l3), x alr(l1),alr(l2),alr(l3)) if(iq. eq. 3) then iw=2 q(iu)=quad(is,iw,slr,opk(kx,l2),opk(kx,l3),opk(kx,l4), x alr(l2),alr(l3),alr(l4)) endif is=1 enddo c is=0 iw=1 c..... k and Dlog(k)/dlog(T6) in lower-right 3x3(i=i1,i1+2 j=j1,j1+2) opact=quad(is,iw,slt,h(1),h(2),h(3),alt(k1),alt(k2),alt(k3)) dopact=dkap dkap1=dkap if (iq. eq. 3) then c..... k and Dlog(k)/Dlog(T6) upper-right 3x3(i=i1+1,i1+3 j=j1,j1+2) opactq=quad(is,iw,slt,q(1),q(2),q(3),alt(k1),alt(k2),alt(k3)) dkapq1=dkap endif if(ip .eq. 3) then c..... k and Dlog(k)/Dlog(T6) in lower-left 3x3. opact2=quad(is,iw,slt,h(2),h(3),h(4),alt(k2),alt(k3),alt(k4)) dkap2=dkap c..... k and Dlog(k)/Dlog(T6) smoothed in left 3x4 dix=(alt(k3)-slt)*dfs(k3) dopact=dkap1*dix+dkap2*(1.-dix) opact=opact*dix+opact2*(1.-dix) if(iq .eq. 3) then c..... k and Dlog(k)/Dlog(T6) in upper-right 3x3. opactq2=quad(is,iw,slt,q(2),q(3),q(4),alt(k2),alt(k3),alt(k4)) dkapq2=dkap dopactq=dkapq1*dix+dkapq2*(1.-dix) opactq=opactq*dix+opactq2*(1.-dix) endif endif c iu=0 do lx=l1,l1+iq iw=1 iu=iu+1 h(iu)=quad(is,iw,slt,opk(k1,lx),opk(k2,lx),opk(k3,lx), x alt(k1),alt(k2),alt(k3)) if(ip .eq. 3) then iw=2 q(iu)=quad(is,iw,slt,opk(k2,lx),opk(k3,lx),opk(k4,lx), x alt(k2),alt(k3),alt(k4)) endif is=1 enddo c is=0 iw=1 c..... k and Dlog(k)/Dlog(R) in lower-left 3x3 opacr=quad(is,iw,slr,h(1),h(2),h(3),alr(l1),alr(l2),alr(l3)) dopacr=dkap if(ip .eq. 3) then opacrq=quad(is,iw,slr,q(1),q(2),q(3),alr(l1),alr(l2),alr(l3)) c..... k and Dlog(k)/Dlog(R) in upper-left 3x3. dopacrq=dkap endif if(iq .eq. 3) then c..... k and Dlog(k)/Dlog(R) in lower-right 3x3. opact2=quad(is,iw,slr,h(2),h(3),h(4),alr(l2),alr(l3),alr(l4)) dix2=(alr(l3)-slr)*dfsr(l3) dopacr=dopacr*dix2+dkap*(1.-dix2) c..... k and Dlog(k)/Dlog(T6) smoothed in both log(T6) and log(R) dopact=dopact*dix2+dopactq*(1.-dix2) opact=opact*dix2+opactq*(1.-dix2) endif if(ip .eq. 3) then if(iq .eq. 3) then c..... k and Dlog(k)/Dlog(R) in upper-right 3x3. opacrq=quad(is,iw,slr,q(2),q(3),q(4),alr(l2),alr(l3),alr(l4)) c..... Dlog(k)/Dlog(R) smoothed in both log(T6) and Log(R). dopacrq=dopacrq*dix2+dkap*(1.-dix2) endif dopacr=dopacr*dix+dopacrq*(1.-dix) endif dopactd=dopact-3.*dopacr if (opact .gt. 1.e+15) then write(*,'("Interpolation indices out of range", x ";please report conditions.")') stop endif if (opact .gt. 9) then opact=30. dopact=99. dopactr=99. dopactd=99. endif return end c************************************************************************ subroutine readco c..... The purpose of this subroutine is to read the data tables save parameter (ismdata=0) ! modified parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) common/aa/ q(4), h(4), xcd(mc),xod(mc), xc(mc),xo(mo) . ,xcs(mc),xos(mo), cxd(mc),oxd(mo),cx(mc),ox(mo),zzz,xxh,xx(mx) . ,nc,no common/a/ co(mx,mc,mo,nt,nr), diag(mx,mc,nt,nr), index(101), . t6list(nt),alr(nr),n(mx,mc),alt(nt),diago(mx,mo,nt,nr),opk(nt,nr) . ,dfs(nt),dfsr(nr),a(3,mx),b(3),m,mf,xa(8),alrf(nrm),cof(ntm,nrm) .,t6listf(ntm),opk2(nt,nr),dfsx(mx) common/b/ itab(mx,ntabs),nta(nrm),x(mx,ntabs),y(mx,ntabs), . zz(mx,ntabs),xca(mx,ntabs),xoa(mx,ntabs) common/alink/ NTEMP,NSM,nrlow,nrhigh,RLE,t6arr(100),coff(100,nr) COMMON/CST/NRL,RLS,nset,tmax ! modified common/e/ opact,dopact,dopacr,dopacrd character*1 dumarra(250) common/recoin/ itimeco,mxzero c if (itimeco .ne. 12345678) then do i=1,mx do j=1,mc do k=1,mo do l=1,nt do mq=1,nr co(i,j,k,l,mq)=1.e+35 enddo enddo enddo enddo enddo do i=1,mx do j=1,mc do l=1,nt do mq=1,nr diag(i,j,l,mq)=1.e+35 diago(i,j,l,mq)=1.e+35 enddo enddo enddo enddo itimeco=12345678 endif do 20 j=1,nc-1 do 21 i=1,nc if(xcd(j) .ge. xc(i)) then n(m,j)=i+1 if(xcd(j) .lt. xc(i)+1.e-6) n(m,j)=i endif 21 continue 20 continue n(m,nc)=0 c close (2) c..... read X=0.0 tables if(m .eq. 1) open(2, FILE='codataa') c..... read X=0.03 tables if(m .eq. 2) open(2, FILE='codatab') c..... read X=0.10 tables if(m .eq. 3) open(2, FILE='codatac') c..... read X=0.35 tables if(m .eq. 4) open(2, FILE='codatad') c.....read X=0.70 tables if(m .eq. 5) open(2, FILE='codatae') c c read header read(2,'(a)') (dumarra(i),i=1,240) c int=0 do 1 j=1,no-1 do 2 i=1,n(m,j) int=int+1 c read(2,'(f10.5)') dum read (2,'(7x,i3,26x,f6.4,3x,f6.4,3x,f6.4,5x,f6.4,5x,f6.4)') x itab(m,int),x(m,int),y(m,int),zz(m,int),xca(m,int),xoa(m,int) xca(m,int)=min(xca(m,int),1.-x(m,int)-zz(m,int)-xoa(m,int)) read(2,'(f10.5)') dum,dum,dum read(2,'(4x,f6.1,18f7.1)') (alrf(kk),kk=1,nrm) read(2,'(f10.5)') dum do k=1,ntm read(2,'(f4.2,19f7.3)') altin,(cof(k,l), l=1,nrm) do ll=1,nrm ! modified coff(k,ll)=cof(k,ll) enddo enddo if (isett6 .ne. 1234567) then do k=1,ntm t6arr(k)=t6list(k) enddo endif isett6=1234567 if (ismdata .eq. 0) then if ((nrm .ne. nr) .or. (ntm .ne. nt)) then write (*,'("Not set up to smooth data with reduced ", x "T-Rho grid")') stop endif tmax=10. ! modified nset=67 ! 65 in earlier version RLS=-8. nsm=1 RLE=1.0 nrlow=1 nrhigh=2*(RLE-RLS)+1 call opaltab !modified endif ll=1 do 110 kk=nrb,nre alr(ll)=alrf(kk) do k=1,nt if (ismdata .eq. 0) then if ((m .eq. 1) .and. (k .le. 9)) then co(m,i,j,k,ll)=cof(k+ntb-1,kk) else co(m,i,j,k,ll)=coff(k+ntb-1,kk) endif else co(m,i,j,k,ll)=coff(k+ntb-1,kk) endif enddo 110 ll=ll+1 2 continue 1 continue if(x(m,1) .ne. xa(m)) then write(*,'(" X in the codata? file does not match xa(m)")') stop endif c do i=1, nc-1 do k=1,nt do l=1,nr diag(m,i,k,l)=co(m,n(m,i),i,k,l) enddo enddo enddo c do 6 j=1,no-1 int=int+1 read(2,'(f10.5)') dum read (2,'(7x,i3,26x,f6.4,3x,f6.4,3x,f6.4,5x,f6.4,5x,f6.4)') x itab(m,int),x(m,int),y(m,int),zz(m,int),xca(m,int),xoa(m,int) read(2,'(f10.5)') dum,dum,dum read(2,'(4x,f6.1,18f7.1)') (alrf(kk),kk=1,nrm) read(2,'(f10.5)') dum do k=1,ntm read (2,'(f4.2,19f7.3)') dum, (cof(k,l), l=1,nrm) c set up to smooth final "diago" opacity tables do l=1,nrm coff(k,l)=cof(k,l) enddo enddo c smooth final "diago" opacity tables too! if (ismdata .eq. 0) then tmax=10. ! modified nset=67 !65 in earlier version RLS=-8. nsm=1 RLE=1.0 nrlow=1 nrhigh=2*(RLE-RLS)+1 call opaltab !modified do k=3,NTEMP-2 do ll=nrlow,nrhigh c Following skip required because, due to missing data, c the X=0 low T data cannot be smoothed if ((m .eq. 1) .and. (k .le. 9)) then cof(k,ll)=cof(k,ll) else cof(k,ll)=coff(k,ll) endif enddo enddo ll=1 do kk=nrb,nre do k=1,nt diago(m,j,k,ll)=cof(k+ntb-1,kk) enddo ll=ll+1 enddo endif 6 continue do i=2,nt dfs(i)=1./(alt(i)-alt(i-1)) enddo do i=2,nr dfsr(i)=1./(alr(i)-alr(i-1)) enddo istep=-1 if (mx .gt. 1 ) then istep=1 do i=2,mx,istep dfsx(i)=1./(xx(i)-xx(i-1)) enddo endif return end c c************************************************************************ function quad(ic,i,x,y1,y2,y3,x1,x2,x3) c..... this function performs a quadratic interpolation. save common/d/dkap common/coquad/ xx(3),yy(3),xx12(30),xx13(30),xx23(30),xx1sq(30) . ,xx1pxx2(30) xx(1)=x1 xx(2)=x2 xx(3)=x3 yy(1)=y1 yy(2)=y2 yy(3)=y3 if(ic .eq. 0) then xx12(i)=1./(xx(1)-xx(2)) xx13(i)=1./(xx(1)-xx(3)) xx23(i)=1./(xx(2)-xx(3)) xx1sq(i)=xx(1)*xx(1) xx1pxx2(i)=xx(1)+xx(2) endif c3=(yy(1)-yy(2))*xx12(i) c3=c3-(yy(2)-yy(3))*xx23(i) c3=c3*xx13(i) c2=(yy(1)-yy(2))*xx12(i)-(xx1pxx2(i))*c3 c1=yy(1)-xx(1)*c2-xx1sq(i)*c3 dkap=c2+(x+x)*c3 quad=c1+x*(c2+x*c3) return end c c************************************************************************ block data parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) common/aa/ q(4), h(4), xcd(mc),xod(mc), xc(mc),xo(mo) . ,xcs(mc),xos(mo), cxd(mc),oxd(mo),cx(mc),ox(mo),zzz,xxh,xx(mx) . ,nc,no common/a/ co(mx,mc,mo,nt,nr), diag(mx,mc,nt,nr), index(101), . t6list(nt),alr(nr),n(mx,mc),alt(nt),diago(mx,mo,nt,nr),opk(nt,nr) . ,dfs(nt),dfsr(nr),a(3,mx),b(3),m,mf,xa(8),alrf(nrm),cof(ntm,nrm) .,t6listf(ntm),opk2(nt,nr),dfsx(mx) common/b/ itab(mx,ntabs),nta(nrm),x(mx,ntabs),y(mx,ntabs), . zz(mx,ntabs),xca(mx,ntabs),xoa(mx,ntabs) common/aaa/ oxf(mx,mc),cxf(mx,mc),xcdf(mx,mc),xodf(mx,mc), . opl(mx,nt,nr),itime(mx),cxdf(mx,mc),oxdf(mx,mc) common/recoin/ itimeco,mxzero data itime/mx*0/,itimeco/0/ data ( index(i),i=1,101)/1,2,2,3,3,3,3,3,3,3,4,4,4,4,4,4, . 4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6, . 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, . 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7/ data (xcs(i),i= 1,mc)/ 0.0,0.01,0.03,0.1,0.2,0.4,0.6,1.0/ data (xos(i),i= 1,mc)/0.0,0.01,0.03,0.1,0.2,0.4,0.6,1.0/ data (xa(i),i=1,5)/0.0,0.03,0.1,0.35,0.7/ data (nta(i),i=1,nrm)/14*70,69,64,60,58,57/ end C********************************************************************* SUBROUTINE SPLINE(X,Y,N,Y2) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(N),U(NMAX) C C FIRST DERIVATIVES AT END POINTS USING CUBIC FIT YP1=((Y(3)-Y(1))*(X(2)-X(1))**2 + -(Y(2)-Y(1))*(X(3)-X(1))**2)/ + ((X(3)-X(1))*(X(2)-X(1))*(X(2)-X(3))) YPN=((Y(N-2)-Y(N))*(X(N-1)-X(N))**2 + -(Y(N-1)-Y(N))*(X(N-2)-X(N))**2)/ + ((X(N-2)-X(N))*(X(N-1)-X(N))*(X(N-1)-X(N-2))) C Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END C********************************************************************* SUBROUTINE SPLINT(XA,YA,N,Y2A,X,Y,YP) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) IF (H.EQ.0.) PAUSE 'Bad XA input.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. YP=0.05* ( (-YA(KLO)+YA(KHI))/H + + ( -(3*A**2-1)*Y2A(KLO) + +(3*B**2-1)*Y2A(KHI) )*H/6. ) RETURN END C********************************************************************* SUBROUTINE FITY C C THIS ROUTINE MAKES SPLINE FITS FOR F AND FX, AND OBTAINS C FY AND FXY C C COMMON/CST/NRL,RLS,nset,tmax ! modified C PARAMETER(IPR=20) DIMENSION A(IPR),B(IPR),AD(IPR),BD(IPR) COMMON/CF/F(85,IPR),FX(85,IPR),FY(85,IPR),FXY(85,IPR) C DO 30 I=1,nset ! modified DO 10 J=1,NRL A(J)=F(I,J) B(J)=FX(I,J) 10 CONTINUE C CALL GETD(A,NRL,AD,AP1,APN) CALL GETD(B,NRL,BD,BP1,BPN) C FY(I,1)=AP1 FY(I,NRL)=APN FXY(I,1)=BP1 FXY(I,NRL)=BPN DO 20 J=2,NRL-1 FY(I,J)= -A(J)+A(J+1)-2.*AD(J)-AD(J+1) FXY(I,J)=-B(J)+B(J+1)-2.*BD(J)-BD(J+1) 20 CONTINUE 30 CONTINUE C RETURN END C********************************************************************* SUBROUTINE FITX C C THIS ROUTINE IS USED ONLY AFTER SMOOTHING. C ITS FUNCTION IS TO RECOMPUTE FX USING SMOOTHED F. C C PARAMETER(IPR=20) DIMENSION A(85),D(85) C COMMON/CST/NRL,RLS,nset,tmax ! modified COMMON/CF/F(85,IPR),FX(85,IPR),FY(85,IPR),FXY(85,IPR) C DO 30 J=1,NRL DO 10 I=1,nset ! modified A(I)=F(I,J) 10 CONTINUE CALL GETD(A,nset,D,AP1,APN) ! modified FX(1,J)=AP1 FX(nset,J)=APN ! modified DO 20 I=2,nset-1 ! modified FX(I,J)=-A(I)+A(I+1)-2.*D(I)-D(I+1) 20 CONTINUE 30 CONTINUE C RETURN END C C************************************************************************ SUBROUTINE GETD(F,N,D,FP1,FPN) C C SIMPLIFIED CODE FOR SPLINE COEFFICIENTS, FOR CASE OF INTERVALS C OF UNITY. C C DIMENSION F(N),D(N),T(85) C FP1=(-11.*F(1)+18.*F(2)-9.*F(3)+2.*F(4))/6. FPN=(11.*F(N)-18.*F(N-1)+9.*F(N-2)-2.*F(N-3))/6. C D(1)=-.5 T(1)=.5*(-F(1)+F(2)-FP1) C DO 10 J=2,N-1 D(J)=-1./(4.+D(J-1)) T(J)=-D(J)*(F(J-1)-2.*F(J)+F(J+1)-T(J-1)) 10 CONTINUE C D(N)=(FPN+F(N-1)-F(N)-T(N-1))/(2.+D(N-1)) C DO 20 J=N-1,1,-1 D(J)=D(J)*D(J+1)+T(J) 20 CONTINUE C RETURN END C C********************************************************************* SUBROUTINE INTERP(FLT,FLRHO,G,DGDT,DGDRHO,IERR) C C GIVEN F,FX,FY AND FXY ON THE GRID POINTS, THIS ROUTINE C DOES BI-CUBIC INTERPOLATIONS USING METHODS DESCRIBED IN C Numerical Recipes, PP. 118 TO 120 C C PARAMETER(IPR=20) COMMON/CF/ F(85,IPR),FX(85,IPR),FY(85,IPR),FXY(85,IPR) DIMENSION B(16) LOGICAL IERR C COMMON/CST/NRL,RLS,nset,tmax ! modified C C EXTREME LIMITS ALLOWED ARE:- C (3.800-0.0125) TO (8.000+0.0125) FOR LOG10(T) C (RLS-0.125) TO (RLE+0.1254) FOR LOG10(R) C (ALLOWING FOR SMALL EXTRAPOLATIONS BEYOND TABULAR VALUES) C C FUNCTION DEFINITIONS FOR CUBIC EXPANSION C FF(S,T)= B( 1)+T*(B( 2)+T*(B( 3)+T*B( 4))) + +S*( B( 5)+T*(B( 6)+T*(B( 7)+T*B( 8))) + +S*( B( 9)+T*(B(10)+T*(B(11)+T*B(12))) + +S*( B(13)+T*(B(14)+T*(B(15)+T*B(16))) ))) C FFX(S,T)= B( 5)+T*(B( 6)+T*(B( 7)+T*B( 8))) + +S*( 2*(B( 9)+T*(B(10)+T*(B(11)+T*B(12)))) + +S*( 3*(B(13)+T*(B(14)+T*(B(15)+T*B(16)))) )) C FFY(S,T)= B( 2)+S*(B( 6)+S*(B(10)+S*B(14))) + +T*( 2*(B( 3)+S*(B( 7)+S*(B(11)+S*B(15)))) + +T*( 3*(B( 4)+S*(B( 8)+S*(B(12)+S*B(16)))) )) C FFXY(S,T)= B( 6)+T*(2*B( 7)+3*T*B( 8)) + +S*( 2*B(10)+T*(4*B(11)+6*T*B(12)) + +S*( 3*B(14)+T*(6*B(15)+9*T*B(16)) )) C C IERR=.FALSE. C X=20.*(FLT-3.800)+1 FLR=FLRHO+18.-3.*FLT Y=2*( FLR - RLS )+1 C IF(X.LT.2.)THEN IF(X.LT.0.75)THEN IERR=.TRUE. ELSE I=1 ENDIF ELSEIF(X.GT.84)THEN IF(X.GT.85.25)THEN IERR=.TRUE. ELSE I=84 ENDIF ELSE I=X ENDIF U=X-I C IF(Y.LT.2.)THEN IF(Y.LT.0.75)THEN IERR=.TRUE. ELSE J=1 ENDIF ELSEIF(Y.GT.NRL-1)THEN IF(Y.GT.NRL+.25)THEN IERR=.TRUE. ELSE J=NRL-1 ENDIF ELSE J=Y ENDIF V=Y-J C IF(IERR)THEN G=9.999 DGDT=9.999 DGDRHO=9.999 RETURN ENDIF C C C GIVEN FUNCTIONS AND DERIVATIVES AT GRID POINTS, COMPUTE COEFFICIENTS. B(1)=F(I,J) B(2)=FY(I,J) B(3)=3*(-F(I,J)+F(I,J+1))-2*FY(I,J)-FY(I,J+1) B(4)=2*(F(I,J)-F(I,J+1))+FY(I,J)+FY(I,J+1) C B(5)=FX(I,J) B(6)=FXY(I,J) B(7)=3*(-FX(I,J)+FX(I,J+1))-2*FXY(I,J)-FXY(I,J+1) B(8)=2*(FX(I,J)-FX(I,J+1))+FXY(I,J)+FXY(I,J+1) C B(9)=3*(-F(I,J)+F(I+1,J))-2*FX(I,J)-FX(I+1,J) B(10)=3*(-FY(I,J)+FY(I+1,J))-2*FXY(I,J)-FXY(I+1,J) B(11)=9*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1)) + +6*(FX(I,J)-FX(I,J+1)+FY(I,J)-FY(I+1,J)) + +4*FXY(I,J) + +3*(FX(I+1,J)-FX(I+1,J+1)-FY(I+1,J+1)+FY(I,J+1)) + +2*(FXY(I,J+1)+FXY(I+1,J)) + +FXY(I+1,J+1) B(12)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1)) + +4*(-FX(I,J)+FX(I,J+1)) + +3*(-FY(I,J)+FY(I+1,J)+FY(I+1,J+1)-FY(I,J+1)) + +2*(-FX(I+1,J)+FX(I+1,J+1)-FXY(I,J)-FXY(I,J+1)) + -FXY(I+1,J)-FXY(I+1,J+1) C B(13)=2*(F(I,J)-F(I+1,J))+FX(I,J)+FX(I+1,J) B(14)=2*(FY(I,J)-FY(I+1,J))+FXY(I,J)+FXY(I+1,J) B(15)=6*(-F(I,J)+F(I+1,J)-F(I+1,J+1)+F(I,J+1)) + +4*(-FY(I,J)+FY(I+1,J)) + +3*(-FX(I,J)-FX(I+1,J)+FX(I+1,J+1)+FX(I,J+1)) + +2*(FY(I+1,J+1)-FY(I,J+1)-FXY(I,J)-FXY(I+1,J)) + -FXY(I+1,J+1)-FXY(I,J+1) B(16)=4*(F(I,J)-F(I+1,J)+F(I+1,J+1)-F(I,J+1)) + +2*(FX(I,J)+FX(I+1,J)-FX(I+1,J+1)-FX(I,J+1) + +FY(I,J)-FY(I+1,J)-FY(I+1,J+1)+FY(I,J+1)) + +FXY(I,J)+FXY(I+1,J)+FXY(I+1,J+1)+FXY(I,J+1) C C GET G=LOG10(ROSS), DGDT=d LOG10(ROSS)/d LOG10(T), C DGDRHO=d LOG10(ROSS)/d LOG10(RHO) G=FF(U,V) DGDT=20.*FFX(U,V)-6.*FFY(U,V) DGDRHO=2.*FFY(U,V) C C RETURN END C********************************************************************* SUBROUTINE SMOOTH C C THIS SUBROUTINE USES A 2-DIMENSIONAL GENERALISATION OF THE SMOOTHING C TECHNIQUES DESCRIBED ON PP. 644 TO 649 OF Numerical Recipes. C C CONSIDER THE 25 POINTS DEFINED BY C I+n, n=-2,-1,0,1,2 AND J+m, m=-2,-1,0,1,2. C THE FUNCTION TO BE SMOOTHED IS FITTED TO A BI-CUBIC, INVOLVING C 16 COEFFICIENTS, USING TECHNIQUES OF LEAST-SQUARES. THE SMOOTHED C FUNCTION (TEMPORARILY STORED IN FXY) IS GIVEN BY THE FITTED VALUE C AT THE POINT I AND J. C C THE FITTING IS SHIFTED FOR POINTS CLOSE TO BOUNDARIES. C C PARAMETER(IPR=20) C COMMON/CST/NRL,RLS,nset,tmax ! modified COMMON/CF/F(85,IPR),FX(85,IPR),FY(85,IPR),FXY(85,IPR) C DIMENSION GAM(6) DATA GAM/+0.0073469388,-0.0293877551,-0.0416326531, + +0.1175510204,+0.1665306122,+0.2359183673/ DIMENSION BET(11) DATA BET/ + -0.0048979592,-0.0661224490,-0.0293877551,+0.0195918367, + 0.2644897959,+0.1175510204,-0.0783673469,+0.0277551020, + 0.3746938776,+0.1665306122,-0.1110204082/ DIMENSION ALP(11) DATA ALP/ + -0.0844897959,-0.0048979592,+0.0073469388,+0.0012244898, + 0.3379591837,+0.0195918367,-0.0293877551,+0.4787755102, + 0.0277551020,-0.0416326531,-0.0069387755/ C C DO 20 I=3,nset-2 C J=1 FXY(I,J)= + ALP(1)*( F(I-2,J )+F(I+2,J ) ) + +ALP(2)*( F(I-2,J+1)+F(I+2,J+1)+F(I-2,J+3)+F(I+2,J+3) + +F(I-1,J+4)+F(I+1,J+4) ) + +ALP(3)*( F(I-2,J+2)+F(I+2,J+2) ) + +ALP(4)*( F(I-2,J+4)+F(I+2,J+4) ) + +ALP(5)*( F(I-1,J )+F(I+1,J ) ) + +ALP(6)*( F(I-1,J+1)+F(I+1,J+1)+F(I-1,J+3)+F(I+1,J+3) ) + +ALP(7)*( F(I-1,J+2)+F(I+1,J+2) ) + +ALP(8)* F(I ,J ) + +ALP(9)*( F(I ,J+1)+F(I ,J+3) ) + +ALP(10)* F(I ,J+2) +ALP(11)*F(I ,J+4) C J=2 FXY(I,J)= + BET(1)*( F(I-2,J-1)+F(I+2,J-1)+F(I-2,J+3)+F(I+2,J+3) ) + +BET(2)*( F(I-2,J )+F(I+2,J ) ) + +BET(3)*( F(I-2,J+1)+F(I+2,J+1) ) + +BET(4)*( F(I-2,J+2)+F(I+2,J+2)+F(I-1,J-1)+F(I+1,J-1) + +F(I-1,J+3)+F(I+1,J+3) ) + +BET(5)*( F(I-1,J )+F(I+1,J ) ) + +BET(6)*( F(I-1,J+1)+F(I+1,J+1) ) + +BET(7)*( F(I-1,J+2)+F(I+1,J+2) ) + +BET(8)*( F(I ,J-1)+F(I ,J+3) ) + +BET(9)*F(I ,J ) +BET(10)*F(I ,J+1) +BET(11)*F(I ,J+2) C DO 10 J=3,NRL-2 FXY(I,J)= + GAM(1)*( F(I-2,J-2)+F(I-2,J+2)+F(I+2,J-2)+F(I+2,J+2) ) + +GAM(2)*( F(I-2,J+1)+F(I-2,J-1)+F(I-1,J-2)+F(I-1,J+2) + +F(I+1,J-2)+F(I+1,J+2)+F(I+2,J-1)+F(I+2,J+1) ) + +GAM(3)*( F(I-2,J )+F(I+2,J )+F(I ,J-2)+F(I ,J+2) ) + +GAM(4)*( F(I-1,J-1)+F(I-1,J+1)+F(I+1,J-1)+F(I+1,J+1) ) + +GAM(5)*( F(I-1,J )+F(I ,J-1)+F(I ,J+1)+F(I+1,J ) ) + +GAM(6)* F(I ,J ) 10 CONTINUE C J=NRL-1 FXY(I,J)= + BET(1)*( F(I-2,J+1)+F(I+2,J+1)+F(I-2,J-3)+F(I+2,J-3) ) + +BET(2)*( F(I-2,J )+F(I+2,J ) ) + +BET(3)*( F(I-2,J-1)+F(I+2,J-1) ) + +BET(4)*( F(I-2,J-2)+F(I+2,J-2)+F(I-1,J+1)+F(I+1,J+1) + +F(I-1,J-3)+F(I+1,J-3) ) + +BET(5)*( F(I-1,J )+F(I+1,J ) ) + +BET(6)*( F(I-1,J-1)+F(I+1,J-1) ) + +BET(7)*( F(I-1,J-2)+F(I+1,J-2) ) + +BET(8)*( F(I ,J+1)+F(I ,J-3) ) + +BET(9)*F(I ,J ) +BET(10)*F(I ,J-1) +BET(11)*F(I ,J-2) C J=NRL FXY(I,J)= + ALP(1)*( F(I-2,J )+F(I+2,J ) ) + +ALP(2)*( F(I-2,J-1)+F(I+2,J-1)+F(I-2,J-3)+F(I+2,J-3) + +F(I-1,J-4)+F(I+1,J-4) ) + +ALP(3)*( F(I-2,J-2)+F(I+2,J-2) ) + +ALP(4)*( F(I-2,J-4)+F(I+2,J-4) ) + +ALP(5)*( F(I-1,J )+F(I+1,J ) ) + +ALP(6)*( F(I-1,J-1)+F(I+1,J-1)+F(I-1,J-3)+F(I+1,J-3) ) + +ALP(7)*( F(I-1,J-2)+F(I+1,J-2) ) + +ALP(8)* F(I ,J ) + +ALP(9)*( F(I ,J-1)+F(I ,J-3) ) + +ALP(10)* F(I ,J-2) +ALP(11)*F(I ,J-4) C 20 CONTINUE C DO 40 I=3,nset-2 ! modified DO 30 J=1,NRL F(I,J)=FXY(I,J) 30 CONTINUE 40 CONTINUE C RETURN END C C********************************************************************* subroutine opaltab C C CODE FOR FITTING AND SMOOTHING OPAL DATA. ADAPTED FROM A CODE C WRITTEN BY MIKE SEATON(obtained june 1993) C C OPAL DATA. C ASSUMES FIRST T6=0.006, LAST T6=10.OR 0.04). Depending on position C in the table. C USES RECTANGULAR ARRAY FOR VARIABLES T6 AND LOG10(R) C C (1) NSM=NUMBER OF PASSES THROUGH SMOOTHING FILTER. C USE OF NSM=1 OR 2 IS RECOMMENDED. C NO SMOOTHING WITH NSM=0 C (2) RANGE FOR LOG10(R), C RLS=FIRST VALUE, RLE=LAST VALE C (RLS MUST BE FIRST VALUYE IN TABLE) C C SUBROUTINE INTERP C AFTER PROCESSING, DATA ARE IN A FORM FOR USE OF C SUBROUTINE INTERP C WHICH GIVES LOG(ROSS) AND TWO FIRST DERIVATIVES FOR ANY C VALUES OF LOG(T) AND LOG(RHO). SEE BELOW FOR FURTHER C EXPLANATION. C C OUTPUT FOR THE CASE OF NSM.GT.0. C INTERP IS USED TO OBTAIN SMOOTHED DATA INTERPOLATED C BACK TO THE ORIGINAL OPAL MESH. TWO FILES ARE WRITTEN. C C C THE SUBROUTINES SPLINE AND SPLINT ARE ADAPTED FROM THOSE GIVE BY C W.H. Press, S.A. Teulolsky, W.T. Vettering and B.P. Flannery, C "Numerical Recipes in FORTRAN", 2nd edn., 1992, C.U.P. C OTHER REFERENCES ARE MADE TO METHODS DESCRIBED IN THAT BOOK. C PARAMETER(IP=100,IPR=20) DIMENSION U(IP),ROSSL(IP,IPR),V(IP),V2(IP) parameter (mx=5,mc=8,mo=8,nrm=19,nrb=1,nre=19,nr=nre+1-nrb . ,ntabs=60,ntm=70,ntb=1,nt=ntm+1-ntb) COMMON/CF/F(85,IPR),FX(85,IPR),FY(85,IPR),FXY(85,IPR) CHARACTER*1 HEAD(100) COMMON/CST/NRL,RLS,nset,tmax ! modified common/alink/ N,NSM,nrlow,nrhigh,RLE,t6arr(100),coff(100,nr) LOGICAL IERR C NRL=2*(RLE-RLS)+1 C C STORE LOG10(T) IN U AND LOG10(ROSS) IN ROSSL C CHECK FIRST VALUE OF T6 T6=t6arr(1) do j=1,NRL ROSSL(1,j)=coff(1,j) enddo if( abs(T6 -.0056341325) .lt. 1.e-8) then U(1)=6.+LOG10(T6) ENDIF C SET ROSSL UP TO T6=t6arr(nset) I=1 5 I=I+1 T6=t6arr(I) do j=1,NRL ROSSL(I,j)=coff(I,j) enddo U(I)=6+LOG10(T6) IF(T6.LT.tmax)GOTO 5 N=I IF(N.GT.IP)THEN PRINT*,' REQUIRE PARAMETER IP OF AT LEAST ',N STOP ENDIF C C C DEFINE VARIABLES C X=20.0*(LOG10(T)-3.80)+1 C Y=2.0*(LOG10(R)-RLS)+1 C USE INDICES I=1 TO nset AND J=1 TO NRL C X AND Y ARE SUCH THAT, ON MESH-POINT (I,J), X=I AND Y=J C OBTAIN:- C F(I,J)=LOG10(ROSS) C FX(I,J)=dF/dX C FY(I,J)=dF/dY C FXY(I,J)=ddF/dXdY C C C FIRST GET F AND FX, INTERPOLATING FROM OPAL T6 TO C INTERVAL OF 0.05 IN LOG10(T). DO 40 J=1,NRL C FOR EACH LOG10(R), STORE LOG10(ROSS) IN V(I) DO 20 I=1,N V(I)=ROSSL(I,J) 20 CONTINUE C C GET FIRST DERIVATIVES AT END POINTS C C GET SECOND DERIVATIVES FOR SPLINE FIT CALL SPLINE(U,V,N,V2) C C INTERPOLATE TO LOG10(T)=FLT, FLT=3.8(0.05)8.0 DO 30 I=1,nset ! modified FLT=3.75+0.05*I CALL SPLINT(U,V,N,V2,FLT,F(I,J),FX(I,J)) 30 CONTINUE C 40 CONTINUE C C C OPTION FOR SMOOTHING IF(NSM.GT.0)THEN DO 35 NS=1,NSM CALL SMOOTH 35 CONTINUE CALL FITX ENDIF C C C GET FY AND FXY CALL FITY C C THE ARRAYS F, FX, FY AND FXY ARE NOW STORED C C CAN NOW DO INTERPOLATIONS USING C CALL INTERP(FLT,FLRHO,G,DGDT,DGDRHO,IERR) C INPUT IS FLT=LOG10(T), FLRHO=LOG10(RHO) C OUTPUT IS G=LOG10(ROSS) C DGDT=dG/d(LOG10(T)) C DGDRHO=dG/d(LOG10(RHO)) C IERR=.TRUE. IF INPUT FLT, FLRHO ARE OUT-OF-RANGE, C ELSE IERR=.FALSE. C C INTERPOLATE BACK TO OPAL POINTS IF(NSM.GT.0)THEN do l=1,NRL coff(1,l)=ROSSL(1,l) enddo DO 70 K=2,N FLT=U(K) DO 50 L=nrlow,nrhigh FLR=RLS+.5*(L-1) FLRHO=FLR-18.+3.*FLT CALL INTERP(FLT,FLRHO,G,DGDT,DGDRHO,IERR) IF(IERR)THEN ENDIF V(L)=G 50 CONTINUE T6=t6arr(K) do l=nrlow,nrhigh coff(K,l)=V(l) enddo 70 CONTINUE ENDIF return C 1000 FORMAT(' SMOOTHED OPAL DATA') 1100 FORMAT(' OPAL DATA, (SMOOTHED-ORIGINAL)') 2000 FORMAT(100A1) 2222 FORMAT(F8.3,20F7.3) 6000 FORMAT(/' FIRST T6=',1P,E10.3,', SHOULD BE 0.006') 6003 FORMAT(/' !!! OUT-OF-RANGE !!!'/' FLT=',1P,E10.3,', FLRHO=',E10.3, + ', FLR=',E10.3) C END