subroutine ingrid(mode,inext,igrd) C ================================== C c setting state parameters for the opacity grid calculations c c input: c temp1 - lowest value of T c temp2 - largest value of T c ntemp - number of temperature values c dens1 - lowest value of the density parameter c dens2 - largest value of the density parameter c ndens - number of the density parameter values c c isdens = 0 - density parameter is electron density c > 0 - density parameter is mass density c < 0 - density parameter is gas pressure c c INCLUDE 'PARAMS.FOR' INCLUDE 'MODELP.FOR' INCLUDE 'LINDAT.FOR' parameter (un=1.,ten15=1.e-15,c18=2.997925e18) real*4 absgrd(mttab,mrtab,mfgrid),dtim common/alsave/ALAM0s,ALASTs,CUTOF0s,CUTOFSs,RELOPs,SPACEs common/gridp0/tempg(mttab),densg(mttab,mrtab),elecgr(mttab,mrtab), * densg0(mttab),temp1,ntemp,ndens,nden(mttab) common/gridf0/wlgrid(mfgrid),nfgrid common/fintab/absgrd common/prfrgr/ipfreq,indext,indexn common/igrddd/igrdd,irelin common/initab/absop(msftab),wltab(msftab), * nfrtab(mttab,mrtab),inttab common/elecm0/elecm(mdepth) common/timeta/dtim common/relabu/relabn(matom),popul0(mlevel,1) dimension abgrd(mfgrid),xli(3) character*(80) tabname common/tabout/tabname,ibingr,idens dimension templ(mttab) c c -------------- c initialization c -------------- c igrdd=igrd if(mode.eq.0) then c read(2,*) ntemp,temp1,temp2 read(2,*) idens if(idens.lt.10) then read(2,*) ndens,dens1,dens2 else if(idens.lt.20) then read(2,*) ndens,densl1,densl2,densu1,densu2 else do it=1,ntemp read(2,*) ndens,densl,densu densg(it,1)=densl densg(it,ndens)=densu nden(it)=ndens end do end if if(idens.lt.20) then do it=1,ntemp nden(it)=ndens end do end if if(ifeos.le.0) then read(2,*) nfgrid,inttab,wlam1,wlam2 read(2,*) tabname,ibingr end if c irsct=0 irsche=0 irsch2=0 c wl1=log(wlam1) wl2=log(wlam2) dwl=(wl2-wl1)/(nfgrid-1) do i=1,nfgrid wlgrid(i)=exp(wl1+(i-1)*dwl) end do c if(temp1.gt.0.) then at1=log(temp1) at2=log(temp2) dt=0. if(ntemp.gt.1) dt=(at2-at1)/(ntemp-1) do i=1,ntemp templ(i)=at1+(i-1)*dt tempg(i)=exp(templ(i)) end do if(idens.lt.10) then at1=log(dens1) at2=log(dens2) dr=0. ndens=nden(1) if(ndens.gt.0) dr=(at2-at1)/(ndens-1) do i=1,ntemp do j=1,ndens densg(i,j)=exp(at1+(j-1)*dr) end do end do else if(idens.lt.20) then rhol1=log(densl1) rhol2=log(densl2) rhou1=log(densu1) rhou2=log(densu2) do i=1,ntemp ndens=nden(i) dens1=rhol1+(rhou1-rhol1)/(at2-at1)*(templ(i)-at1) dens2=rhol2+(rhou2-rhol2)/(at2-at1)*(templ(i)-at1) dr=0. if(ndens.gt.1) dr=(dens2-dens1)/(ndens-1) do j=1,ndens densg(i,j)=exp(dens1+(j-1)*dr) end do end do else do i=1,ntemp ndens=nden(i) at1=log(densg(i,1)) at2=log(densg(i,ndens)) dr=0. if(ndens.gt.0) dr=(at2-at1)/(ndens-1) do j=2,ndens-1 densg(i,j)=exp(at1+(j-1)*dr) end do end do end if c write(6,621) ntemp,nden(1) do i=1,ntemp ndens=nden(i) write(6,622) tempg(i),(log10(densg(i,j)),j=1,ndens) end do 621 format(/' COMPUTING AN OPACITY TABLE WITH GRID PARAMETERS:'/ * ' ===== ntemp, ndens ',2i4) 622 format(f10.1,20f8.2) else call inpmod ntemp=nd ndens=1 do it=1,ntemp tempg(it)=temp(it) densg0(it)=dens(it) densg(it,1)=dens(it) elecm(it)=elec(it) end do if(ifeos.le.0) then write(6,621) ntemp,ndens do i=1,ntemp write(6,622) tempg(i),densg0(i) end do end if ndens=1 idens=2 end if c nd=1 idstd=1 inext=1 frmx=0. frmn=1.e20 idens0=mod(idens,10) c indext=1 indexn=1 ipfreq=0 irelin=1 temp(1)=tempg(indext) c write(6,646) indext,temp(1), * indexn,densg(indext,indexn) 646 format(/' ************************************', * /' GRID POINT OF THE OPACITY TABLE WITH:'/ * ' INDEX TEMP, T ',i4,f10.1/ * ' INDEX DENS, DENS',I4,1PE10.1, * /' ************************************'/) c if(temp1.le.0.) elec(1)=elecm(indext) call densit(densg(indext,indexn),idens0) if(ntemp.eq.1.and.ndens.eq.1) inext=0 elecgr(indext,indexn)=elec(1) call abnchn(0) return c c --------------------------------------------- c after computing the table for one T-rho pair: c --------------------------------------------- c else if(mode.eq.1) then if(ifeos.le.0) then c call timing(1,igrd+1) c do i=1,3 xli(i)=0. end do do i=1,nmlist xli(i)=float(nlinmt(i))*1.e-3 end do c if(imode.ge.-5) then if(indext.eq.1.and.indexn.eq.1) * write(29,625) write(29,626) indext,indexn,temp(1),dens(1),elec(1), * float(nlin0)*1.e-3, * (xli(i),i=1,3),dtim 625 format(' it ir t rho elec',6x, * ' atomic molec1 molec2 molec3 time'/) 626 format(2i4,f9.2,1p2e10.2,2x,0pf8.1,2x,3f8.1,2x,f8.2) else alam0=alam0s if(alam0s.eq.0.) alam0=5.e7/temp(1)/10. if(alam0s.lt.0.) alam0=-5.e7/temp(1)/alam0s alast=alasts if(alasts.eq.0.) alast=5.e7/temp(1)*20. if(alasts.lt.0.) alast=-5.e7/temp(1)*alasts if(alast.gt.1.e5) alast=1.e5 write(29,629) temp(1),elec(1),dens(1), * alam0,alast end if 629 format(1p3e11.3,0pf9.3,0pf12.3) c c ------------------------------------------------ c interpolate and store previously computed table c ------------------------------------------------ c nfr=ipfreq nfrtab(indext,indexn)=ipfreq write(*,*) 'indext,indexn,nfreq',indext,indexn,ipfreq write(*,*) 'nfr,nfgrid',nfr,nfgrid c if(inttab.eq.1) then c call interp(wltab,absop,wlgrid,abgrd,nfr,nfgrid,2,0,0) call intrp(wltab,absop,wlgrid,abgrd,nfr,nfgrid) else ij=0 ijgrd=0 30 continue ijgrd=ijgrd+1 wlgr=0.5*(wlgrid(ijgrd)+wlgrid(ijgrd+1)) isum=0 sum=0. 40 continue ij=ij+1 if(ij.gt.nfr) go to 50 wlt=wltab(ij) abl=absop(ij) if(wlt.le.wlgr) then sum=sum+exp(abl) isum=isum+1 go to 40 end if if(isum.gt.0) then abgrd(ijgrd)=log(sum/float(isum)) else abg=abl+(absop(ij+1)-abl)/(wltab(ij+1)-wlt)*(wlgr-wlt) abgrd(ijgrd)=abg c write(*,*) 'grd',ij,absop(ij+1),abl,wltab(ij+1), c * wlt,wlgr,abg,abgrd(ijgrd),ijgrd end if if(ijgrd.lt.nfgrid) then ij=ij-1 go to 30 else if(ijgrd.eq.nfgrid) then wlgr=wlgrid(nfgrid) sum=0. isum=0 if(ij.lt.nfr) ij=ij-1 go to 40 end if end if 50 continue c do ij=1,nfgrid absgrd(indext,indexn,ij)=real(abgrd(ij)) end do absgrd(indext,indexn,nfgrid)=absgrd(indext,indexn,nfgrid-1) end if c c ------------------------------ c prepare values for a new table c ------------------------------ c ipfreq=0 ndens=nden(indext) if(indexn.lt.ndens) then indexn=indexn+1 rho=densg(indext,indexn) write(6,646) indext,tempg(indext), * indexn,densg(indext,indexn) call densit(rho,idens0) inext=1 else indexn=1 irelin=1 if(indext.lt.ntemp) then indext=indext+1 temp(1)=tempg(indext) if(temp1.le.0.) then densg(indext,indexn)=densg0(indext) elec(1)=elecm(indext) end if rho=densg(indext,indexn) write(6,646) indext,tempg(indext), * indexn,densg(indext,indexn) call densit(rho,idens0) inext=1 else inext=0 end if end if if(inext.eq.1) then rewind(19) if(inlist.lt.0) rewind(19) end if c elecgr(indext,indexn)=elec(1) c call abnchn(0) id=1 do i=1,4 do j=i+1,22 call hydtab(i,j,id) end do end do end if c return end