SUBROUTINE TABINI C ================= C C Initialization and reading of the opacity table for thermal processe C and Rayleigh scattering c raytab: scattering opacities in cm^2/gm at 5.0872638d14 Hz (sodium D) c (NOTE: Quantities in rayleigh.tab are in log_e) C c tempvec: array of temperatures c rhomat: array of densities (gm/cm^3) c nu: array of frequencies c table: absorptive opacities in cm^2/gm c (NOTE: Quantities in absorption.tab are in log_e) C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ATOMIC.FOR' parameter (mtabto=100,mtabro=100) common/intcff/yint(mfreq),jint(mfreq) common/abntab/abunt(matom),abuno(matom),tmolit, * iophmt,ioph2t,iophet,iopcht,iopoht, * ioh2mt,ih2h2t,ih2het,ioh2ht,iohhet, * ifmolt common/eletab/elecgr(mtabt,mtabr) CHARACTER(len=100) :: DUM c real*4 absopa0 dimension frlt(mfrtab) dimension absopa0(mtabto,mtabro),indt(mtabto),indr(mtabro) dimension tempve0(100),rhove0(100,100), * elecg0(100,100) dimension rhov(mtabr),nden(mtabt) dimension typa(matom) character*(80) optable character*4 typa c istept=1 istepr=1 istepf=1 c read(15,*,err=10,end=10) optable,ibinop go to 20 10 optable='./data/absopac.dat' 20 if(ibinop.eq.0) then open(53,file=optable,status='old') else open(53,file=optable,form='unformatted',status='old') end if c read(15,*,err=30,end=30) istept,istepr,istepf 30 continue c write(6,601) optable,ibinop 601 format(/' ' * /' OPACITY TABLE: READ FROM THE FILE ',a70/' IBINOP=',i2/ * ' --------------'/) if(istept.gt.1.or.istepr.gt.1.or.istepf.gt.1) then write(6,*) ' BUT ONLY SELECTED DATA USED' write(6,*) end if c c reading the opacity table in the old mode c if(ioptab.lt.0.and.iopold.gt.0) then if(ibinop.eq.0) then read(53,*) numfre0,numtem0,numrh0 read(53,*) read(53,*) (tempve0(i),i=1,numtem0) read(53,*) read(53,*) (rhove0(1,j),j=1,numrh0) else read(53) numfre0,numtem0,numrh0 read(53) (tempve0(i),i=1,numtem0) read(53) (rhove0(1,j),j=1,numrh0) end if c c reading the opacity table in the new mode c else io2hmt=0 ih2h2t=0 ih2het=0 ioh2ht=0 iohhet=0 if(ibinop.eq.0) then read(53,*) read(53,*) do iat=1,matom read(53,'(a100)') dum read(dum,*,iostat=kst) typa(iat),abunt(iat),abuno(iat) if(kst.ne.0) go to 40 end do read(53,*) 40 continue read(53,*) read(53,*) ifmolt,tmolit read(53,*) read(53,*) read(53,'(a100)') dum read(dum,*,iostat=kst) iophmt,ioph2t,iophet,iopcht,iopoht, * ioh2mt,ih2h2t,ih2het,ioh2ht,iohhet if(kst.ne.0) read(dum,*) iophmt,ioph2t,iophet,iopcht,iopoht read(53,*) read(53,*) read(53,*) numfre0,numtem0,numrh0 if(numrh0.gt.0) then read(53,*) read(53,*) (tempve0(i),i=1,numtem0) read(53,*) read(53,*) (rhov(j),j=1,numrh0) nden(1)=numrh0 do j=1,numrh0 do i=1,numtem0 rhove0(i,j)=rhov(j) end do end do do i=1,numtem0 nden(i)=nden(1) end do read(53,*) read(53,*) ((elecg0(i,j),j=1,numrh0),i=1,numtem0) else read(53,*) (nden(i),i=1,numtem0) read(53,*) read(53,*) (tempve0(i),i=1,numtem0) read(53,*) do i=1,numtem0 read(53,*) (rhove0(i,j),j=1,nden(i)) end do read(53,*) do i=1,numtem0 read(53,*) (elecg0(i,j),j=1,nden(i)) end do end if else do iat=1,92 read(53) typa(iat),abunt(iat),abuno(iat) end do read(53) ifmolt,tmolit read(53) iophmt,ioph2t,iophet,iopcht,iopoht, * ioh2mt,ih2h2t,ih2het,ioh2ht,iohhet read(53) numfre0,numtem0,numrh0 if(numrh0.gt.0) then read(53) (tempve0(i),i=1,numtem0) read(53) (rhov(j),j=1,numrh0) read(53) ((elecg0(i,j),j=1,numrh0),i=1,numtem0) nden(1)=numrh0 do j=1,numrh0 do i=1,numtem0 rhove0(i,j)=rhov(j) end do end do do i=1,numtem0 nden(i)=nden(1) end do else read(53) (nden(i),i=1,numtem0) read(53) (tempve0(i),i=1,numtem0) do i=1,numtem0 read(53) (rhove0(i,j),j=1,nden(i)) end do do i=1,numtem0 read(53) (elecg0(i,j),j=1,nden(i)) end do end if end if c end if c C select only a part of tabular data (if required) C j=0 nrmax=0. do it=1,numtem0,istept j=j+1 tempvec(j)=tempve0(it) indt(j)=it k=0. numr=nden(it) nrmax=max(nrmax,numr) do ir=1,numr,istepr k=k+1 rhomat(j,k)=rhove0(it,ir) indr(k)=ir elecgr(j,k)=elecg0(it,ir) end do numrh(j)=k end do numtemp=j if(numrh0.gt.0) numrho=k if(numtemp.gt.mtabt) then write(*,*) 'number of temperatures in opac.table too large' write(*,*) 'numtemp,mtabt = ',numtemp,mtabt write(*,*) 'recompile with MTABT.ge.NUMTEMP in BASICS.FOR' stop end if c if(nrmax.gt.mtabr) then write(*,*) 'number of densities in opac,table too large' write(*,*) 'numrho,mtabr = ',numtemp,mtabr write(*,*) 'recompile with MTABR.ge.NUMRHO in BASICS.FOR' stop end if c ij=0 do k=1,numfre0 if(mod(k,istepf).eq.0) ij=ij+1 end do numfreq=ij if(numfreq.gt.mfreqc) then write(*,*) 'number of wavelengths in opac.table too large' write(*,*) 'numfreq,mfreqc = ',numfreq,mfreqc write(*,*) 'recompile with MFREQC.ge.NUMFREQ in BASICS.FOR' stop end if c write(6,602) numfre0,numtem0,numrh0 602 format(' original number of frequencies, temperatures', * ' densities: ',3i7/) C if(istept.gt.1.or.istepr.gt.1.or.istepf.gt.1) * write(6,603) numfreq,numtemp,numrho 603 format(' modified number of frequencies, temperatures', * ' densities: ',3i7/) C write(*,*) 'temperatures:' write(6,604) (exp(tempvec(i)),i=1,numtemp) write(*,*) 'densities:' do i=1,numtemp numr=numrh(i) write(6,605) exp(tempvec(i)),(exp(rhomat(i,j)),j=1,numr) end do 604 format(10f8.1) 605 format('for T=',f7.1,(1p10e10.2)) c c check the consistency (or a lack thereof) the parameters of c the opacity table c c if(ioptab.gt.0.or.iopold.eq.0) call chctab c RTAB1 = rhomat(1,1) RTAB2 = rhomat(1,numrho) TTAB1 = tempvec(1) TTAB2 = tempvec(numtemp) c if(ibinop.eq.0) then ij=0 nden0=numrh0 do k=1,numfre0 read(53,*) read(53,*) read(53,*) frta if(numrh0.gt.0) then do j = 1, numrh0 read(53,*) (absopa0(i,j),i=1,numtem0) end do else do i=1,numtem0 nden0=nden(i) read(53,*) (absopa0(i,j),j=1,nden0) end do end if if(mod(k,istepf).eq.0.and.frta.lt.frtlim) then ij=ij+1 frtab(ij)=frta frlt(ij)=log10(frta) do i=1,numtemp numr=numrh(i) do j=1,numr absopac(i,j,ij)=absopa0(indt(i),indr(j)) end do end do end if end do numfreq=ij close(53) else ij=0 do k=1,numfre0 read(53) frta if(numrh0.gt.0) then do j = 1, numrh0 read(53) (absopa0(i,j),i=1,numtem0) end do else do i=1,numtem0 nden0=nden(i) read(53,err=11) (absopa0(i,j),j=1,nden0) 11 continue end do end if if(mod(k,istepf).eq.0.and.frta.lt.frtlim) then ij=ij+1 frtab(ij)=frta frlt(ij)=log10(frta) do i=1,numtemp numr=numrh(i) do j=1,numr absopac(i,j,ij)=absopa0(indt(i),indr(j)) end do end do end if end do numfreq=ij c c write(*,*) 'final NUMFREQ, FREQ(1)',numfreq,freq(1) c DO K=1,NUMFREQ IF(FRTAB(K).GT.2.997925E13) THEN K0=K ELSE do i=1,numtemp numr=numrh(i) do j=1,numr absopac(i,j,k)=absopac(i,j,k0) end do end do END IF END DO close(53) end if c frtabm=max(frtab(1),frtab(numfreq)) RETURN END