SUBROUTINE INITIA C ================= C C driver for input and initializations - "new" routine C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ITERAT.FOR' INCLUDE 'ODFPAR.FOR' INCLUDE 'ALIPAR.FOR' PARAMETER (T15=1.D-15) parameter (xcon=8.0935d-21) CHARACTER*20 FINSTD CHARACTER*40 FILEI CHARACTER*4 TYPIOI DIMENSION IGLE(18),IGMN(25),IGFE(26),IGNI(28) dimension fst(mfreq),extrd0(mfreq) COMMON/STRPAR/IMER,ITR,IC,IL,IP,NLASTE,NHOD,LASV COMMON/INUNIT/IUNIT common/freqcl/frmin,frmax,nfrecl DATA IGLE/2,1,2,1,6,9,4,9,6,1,2,1,6,9,4,9,6,1/ DATA IGMN/2,1,2,1,6,9,4,9,6,1,2,1,6,9,4,9,6,1, * 10,21,28,25,6,7,6/ DATA IGFE/2,1,2,1,6,9,4,9,6,1,2,1,6,9,4,9,6,1, * 10,21,28,25,6,25,30,25/ DATA IGNI/2,1,2,1,6,9,4,9,6,1,2,1,6,9,4,9,6,1, * 10,21,28,25,6,25,28,21,10,21/ C C ---------------------- C Basic input parameters C ---------------------- C CALL READBF(0) FRTABM=0. C C for a stellar atmosphere C IF(IDISK.EQ.0) THEN READ(IBUFF,*) TEFF,GRAV WRITE(6,601) TEFF,GRAV GRAV=EXP(2.3025851*GRAV) C C for an accretion disk C ELSE READ(IBUFF,*) XMSTAR,XMDOT,RSTAR,RELDST END IF c c the rest of input is the same for atmospheres and disks C READ(IBUFF,*) LTE,LTGREY READ(IBUFF,*) FINSTD CALL NSTPAR(FINSTD) C C ------------------------ C Specific input for disks C ------------------------ C IF(IDISK.EQ.1) CALL INPDIS C C ------------------------ c Initialize opacity table C ------------------------ c if(ioptab.ne.0) call tabini C C ---------------------------- C Frequency points and weights C ---------------------------- C C NFREAD - number of "continuum" frequency points which are C read. May be less than NFREQ - the actual number of C frequencies, because the frequency points in lines C can be calculated. C FREQ(IJ) - frequency (in sec**-1) of the IJ-th freq. point C W(IJ) - corresponding frequency quadrature weight C READ(IBUFF,*) NFREAD NJREAD=NFREAD if(njread.gt.mfreq) CALL QUIT('njread.gt.mfreq',njread,mfreq) C IF(NJREAD.LE.0) THEN frmin = FRCMIN frmax = FRCMAX NFREQ=-NJREAD if(ifrset.ge.0) then nfreq=ifrset frmin=frtab(numfreq) frmax=frtab(1) end if if(ifrset.lt.0) nfreqc=-ifrset nfreqc=nfreq if(nfreq.eq.1) then freq(1)=frmin w(1)=1. ijali(1)=1 else do ij=1,nfreq fr=log(frmin)+(log(frmax)-log(frmin))*(ij-1)/(nfreq-1) freq(nfreq-ij+1)=exp(fr) ijali(ij)=1 end do end if C C setting simple quadrature weights - trapezoidal integration C if(nfreq.gt.1) then w(1)=0.5*(freq(1)-freq(2)) w(nfreq)=0.5*(freq(nfreq-1)-freq(nfreq)) do ij=2,nfreq-1 w(ij)=0.5*(freq(ij-1)-freq(ij+1)) end do end if ELSE NFREQC=NJREAD END IF C C ---------------------------------------------------- C Initialize 1/i*i, 1/i*i*i; and turbulent velocities C ---------------------------------------------------- C DO I=1,NLMX X=I XI2(I)=UN/(X*X) XI3(I)=XI2(I)/X END DO C IF(ABS(VTB).LT.1.E3) VTB=VTB*1.E5 DO ID=1,ND IF(VTB.GT.0.) VTURB(ID)=VTB IF(IPTURB.EQ.0) VTURB(ID)=0. VTURBS(ID)=ABS(VTB) END DO C C ---------------------------------------------------- C Input parameters for explicit and non-explicit atoms C ---------------------------------------------------- C C Input parameters are read by procedure STATE C (see description there) C CALL STATE(0,ID,X1,X2) ID=1 WRITE(6,602) YTOT(ID),WMY(ID),WMM(ID) C C check consistency of the opacity table C c if(ioptab.gt.0.and.iopold.eq.0) call chctab c DO I=1,MLEVEL ILK(I)=0 DO J=1,MLEVEL ITRA(J,I)=0 END DO END DO c if(nfrecl.gt.nfreq) nfrecl=nfreq do ij=1,nfrecl ijali(ij)=0 end do C IF(NFFIX.EQ.2) THEN DO IJ=1,NFREQC IJALI(IJ)=1 END DO END IF C C -------------------------------------------------------------- C Input of parameters for explicit ions, levels, and transitions C -------------------------------------------------------------- C ILEV=0 IATLST=0 ION=0 IA=0 IUNIT=20 NATOM=0 NLEVEL=0 10 CONTINUE READ(IBUFF,*,END=20,ERR=20) IATII,IZII,NLEVSI,ILASTI,ILVLIN, * NONSTD,TYPIOI,FILEI IF(ILASTI.EQ.0) THEN ION=ION+1 if(ion.gt.mion) CALL QUIT('ion.gt.mion',ion,mion) IATI(ION)=IATII IZI(ION)=IZII ILTION(ION)=0 if(nlevsi.lt.0) then iltion(ion)=2 nlevsi=-nlevsi end if NLEVS(ION)=NLEVSI TYPION(ION)=TYPIOI FIDATA(ION)=FILEI NLLIM(ION)=ILVLIN INODF1(ION)=0 INODF2(ION)=0 IUPSUM(ION)=0 ICUP(ION)=16 FF(ION)=0. MODEFF=1 NFF=0 IF(IATI(ION).EQ.1.AND.IZI(ION).EQ.0) THEN IUPSUM(ION)=-100 MODEFF=2 END IF IF(IATI(ION).EQ.1.AND.IZI(ION).EQ.-1) THEN IUPSUM(ION)=0 ICUP(ION)=0 MODEFF=3 END IF IF(IATI(ION).EQ.2.AND.IZI(ION).EQ.1) THEN MODEFF=2 ICUP(ION)=32 END IF IF(NONSTD.GT.0) THEN READ(IBUFF,*) IUPSUM(ION),ICUP(ION),MODEFF,NFF ELSE IF(NONSTD.LT.0) THEN READ(IBUFF,*) INODF1(ION),INODF2(ION),FIODF1(ION), * FIODF2(ION),FIBFCS(ION) IKOBS(ION)=IABS(NONSTD) if(ispodf.ge.1) then IF(INODF1(ION).EQ.0) THEN IUNIT=IUNIT+1 INODF1(ION)=IUNIT END IF IF(INODF2(ION).EQ.0) THEN IUNIT=IUNIT+1 INODF2(ION)=IUNIT END IF end if IF(FIBFCS(ION).NE.' ') THEN IUNIT=IUNIT+1 INBFCS(ION)=IUNIT END IF IUPSUM(ION)=1 ICUP(ION)=0 END IF C IF(IATI(ION).EQ.IATLST) THEN NFIRST(ION)=ILEV ELSE NFIRST(ION)=ILEV+1 IATLST=IATI(ION) IA=IATEX(IATLST) N0A(IA)=NFIRST(ION) NATOM=MAX(NATOM,IA) END IF NLAST(ION)=NFIRST(ION)+NLEVS(ION)-1 NNEXT(ION)=NLAST(ION)+1 ILEV=NNEXT(ION) if(ilev.gt.mlevel) CALL QUIT('ilev.gt.mlevel',ilev,mlevel) IZ(ION)=IZI(ION)+1 CHARG2(ION)=IZ(ION)*IZ(ION) IF(IFMOFF.gt.0) THEN NFF=NQUANT(NLAST(ION))+1 IF(NFF.GT.0) FF(ION)=EH/H*IZ(ION)*IZ(ION)/NFF/NFF END IF N0I=NFIRST(ION) N1I=NLAST(ION) NKI=NNEXT(ION) ITRA(NKI,NKI)=MODEFF DO II=N0I,N1I IEL(II)=ION IATM(II)=IA END DO ILK(NKI)=ION IATM(NKI)=IA C IF(NUMAT(IA).EQ.1) THEN IATH=IA IF(IZ(ION).EQ.1) IELH=ION IF(IZ(ION).EQ.0) THEN IELHM=ION IOPHMI=0 END IF END IF IF(NUMAT(IA).EQ.2) THEN IATHE=IA IF(IZ(ION).EQ.1) IELHE1=ION IF(IZ(ION).EQ.2) IELHE2=ION END IF C IF(ION.EQ.1) WRITE(6,603) WRITE(6,604) ION,TYPION(ION),N0I,N1I,NKI,IZ(ION), * IUPSUM(ION),ICUP(ION),FF(ION) C ELSE IF(ILASTI.GT.0) THEN ENION(ILEV)=0. G(ILEV)=ILASTI NQUANT(ILEV)=1 TYPLEV(ILEV)=TYPIOI IMODL(ILEV)=0 IFWOP(ILEV)=0 IEL(ILEV)=ION NKA(IA)=NNEXT(ION) if(modref.ge.0) nref(ia)=nka(ia) IF(ILASTI.EQ.1.AND.IATII.GT.IZII) THEN IF(IATII.LT.25) THEN G(ILEV)=IGLE(IATII-IZII) ELSE IF(IATII.EQ.25) THEN G(ILEV)=IGMN(IATII-IZII) ELSE IF(IATII.EQ.26) THEN G(ILEV)=IGFE(IATII-IZII) ELSE IF(IATII.EQ.28) THEN G(ILEV)=IGNI(IATII-IZII) ENDIF ENDIF ELSE GO TO 20 END IF GO TO 10 20 CONTINUE NION=ION NLEVEL=NKI c if(ioptab.gt.0.and.iopold.eq.0) call chctab write(6,631) * iopadd,irsct,irsche,irsch2,iophmi,ioph2p,iopoh,iopch 631 format(/'iopadd,irsct,irsche,irsch2,iophmi,ioph2p,iopoh,iopch', * 8i6) C IMER=0 ITR=0 ITRX=0 IC=0 IL=0 IP=0 IF(NFREAD.LE.0) THEN NLASTE=NFREQC NFREQ=NFREQC ELSE NLASTE=0 NFREQ=0 END IF LBPFX=.TRUE. C lasv=.false. DO ION=1,NION CALL RDATA(ION) IF(IFMOFF.gt.0) THEN NFF=NQUANT(NLAST(ION))+1 IF(NFF.GT.0) FF(ION)=EH/H*IZ(ION)*IZ(ION)/NFF/NFF END IF END DO C NTRANS=ITR NFREQ=NLASTE NCON=IC if(ntrans.gt.mtrans) CALL QUIT('ntrans.gt.mtrans',ntrans,mtrans) C CALL LEVSET C c WRITE(6,605) DO I=1,NLEVEL IF(I.EQ.1) WRITE(6,605) WRITE(6,606) I,TYPLEV(I),TYPION(IEL(I)),ENION(I)/H,G(I), * NQUANT(I),IEL(I),ILK(I),IATM(I), * IMODL(I),ILTLEV(I),IIEXP(I),IIFOR(I) END DO C C ----------------------------------------------------------- C useful quantities for later use C ----------------------------------------------------------- C DO 125 IT=1,NTRANS FRQMX(IT)=0. LINEXP(IT)=(.NOT.LINE(IT).OR.INDEXP(IT).EQ.0) IF(LINEXP(IT)) GO TO 125 FFMX=0. DO IJ=IFR0(IT),IFR1(IT) FFMX=MAX(FFMX,FREQ(IJ)) END DO FRQMX(IT)=FFMX 125 CONTINUE C C ------------------------------------------------------ C setup continuum frequencies (if they are not read) C ------------------------------------------------------ C if(nfreq.gt.mfreq) CALL QUIT('nfreq.gt.mfreq',nfreq,mfreq) if(nfreqc.gt.mfreqc) CALL QUIT('nfreqc.gt.mfreqc',nfreqc,mfreqc) IF(NFREAD.GT.0.AND.ISPODF.EQ.0) THEN IF(IOPTAB.EQ.0) THEN CALL INIFRC(0) ELSE IF(IOPTAB.GT.0) THEN CALL INIFRT END IF END IF IF(ISPODF.GE.1) CALL INIFRS IF(FRLMAX.LE.0.) FRLMAX=FREQ(1) C C --------------------------------------------------------- C setup depth-independent line profiles (sampling mode) C --------------------------------------------------------- C IF(ISPODF.GE.1) THEN TSTD=0.75*TEFF DO 80 IT=1,NTRANS IF(LINEXP(IT)) GOTO 80 INDXPA=IABS(INDEXP(IT)) IF(INDXPA.GE.2 .AND. INDXPA.LE.4) GO TO 80 CALL DOPGAM(IT,1,TSTD,DOP,AGAM) CALL LINSPL(IT,DOP,AGAM) 80 CONTINUE END IF C call rdatax(0,ic,0) C C check the dimensions C if(nd.gt.mdepth) CALL QUIT('nd.gt.mdepth',nd,mdepth) if(nfreq.gt.mfreq) CALL QUIT('nfreq.gt.mfreq',nfreq,mfreq) if(nfreqc.gt.mfreqc) CALL QUIT('nfreqc.gt.mfreqc',nfreqc,mfreqc) if(nfreqe.gt.mfrex) CALL QUIT('nfreqe.gt.mfrex',nfreqe,mfrex) if(nlevel.gt.mlevel) CALL QUIT('nlevel.gt.mlevel',nlevel,mlevel) if(ntrans.gt.mtrans) CALL QUIT('ntrans.gt.mtrans',ntrans,mtrans) if(natom.gt.matom) CALL QUIT('natom.gt.matom',natom,matom) if(nion.gt.mion) CALL QUIT('nion.gt.mion',nion,mion) if(ibfint.eq.0.and.nfreq.gt.mfreqc) * CALL QUIT('nfreq.gt.mfreqc for ibfint.eq.0',nfreq,mfreqc) if(ntrans.gt.32767) * call quit(' Too many transitions to define ITRLIN as INTEGER*2', * ntrans,ntrans) nncdw=0 do ibft=1,ntranc itr=itrbf(ibft) icdw=mcdw(itr) if(icdw.ge.1) nncdw=nncdw+1 end do if(nncdw.gt.mmcdw) & call quit(' Too many pseudo-continua',nncdw,mmcdw) c C ----------------------------------------------------------- C read the input model C ----------------------------------------------------------- C IF(.NOT.LTGREY) THEN CALL INPMOD IF(ICHANG.NE.0) CALL CHANGE END IF C CALL LINSET(0,0,0,0,0.d0,0.d0,0.d0) C C ----------------------------------------------------------- C There is no additional input for ODF transitions (i.e. ABS(MODE)=3) C from standard input; all information is taken from corresponding C input files (subroutine ODFSET) C C Input parameters for ODF lines or for iron line sampling C ----------------------------------------------------------- C DOPSTD=SQRT(TWO*BOLK*TEFF/HMASS+VTB*VTB) IF(NHOD.GT.0) CALL ODFHYS(DOPSTD) IF(ISPODF.EQ.0) THEN CALL ODFSET ELSE IF(ISPODF.GE.1) THEN CALL IROSET END IF C C ----------------------------------------------------------- C select explicit continuum frequencies C ----------------------------------------------------------- C c IF(NFREAD.GT.0.AND.ISPODF.EQ.0.and.ioptab.eq.0) CALL INIFRC(1) IF(NFREAD.GT.0.AND.ISPODF.GE.0) CALL INIFRC(1) C CALL TRAINI C c C -------------------------------------------------------- c Interpolate the opacity table to current frequencies C -------------------------------------------------------- c if(ioptab.ne.0) then call tabint call rayini do ij=1,nfreqc ijx(ij)=1 end do end if C C ----------------------------------------------------------- C sorting frequencies & new weights C ----------------------------------------------------------- C DO IJ=NFREQC+1,NFREQ IJX(IJ)=1 END DO if(ioptab.ge.0) then CALL SRTFRQ else do ij=1,nfreq ijfr(ij)=ij jik(ij)=ij end do end if C C ----------------------------------------------------------- C other frequency-dependent quantities C ----------------------------------------------------------- C if(ifryb.gt.0.or.nfrecl.ge.nfreq) then do ij=1,nfreq ijali(ij)=0 ijfr(ij)=ij end do end if c CALL CORRWM C do ij=1,nfreq if(icompt.eq.0) then sigec(ij)=sige else c c first-order expression c if(knish.eq.0) then SIGEC(IJ)=SIGE*(un-two*freq(ij)*xcon) c C Use full Klein-Nishina cross section c (Rybicki & Lightman 1975): c else xf=xcon*freq(ij) if(xf.lt.1.d-1) then SIGEC(IJ)=SIGE*(1.-xf*(2.-xf*(26./5.-xf*(13.3 * -xf*(1144./35.-xf*(544./7.-xf*(3784./21. * -xf*(6148./15.-xf*(151552./165. * -xf*111872./55.))))))))) else if(xf.gt.1.d3) then SIGEC(IJ)=SIGE*3./8./xf*(log(2.*xf)+0.5) else SIGEC(IJ)=SIGE*0.75*((1.+xf)/xf**3*(2.*xf*(1.+xf)/ * (1.+2.*xf)-log(1.+2.*xf)) * +0.5*log(1.+2.*xf)/xf * -(1.+3.*xf)/(1.+2.*xf)**2) endif endif endif end do c CALL RTEANG C C ----------------------------------------------------------- C print out some important transition parameters C ----------------------------------------------------------- C if(iptran.gt.0) then WRITE(6,607) NTRANS c cone=3.28805e15*6.6256e-27/109678.758 cone=3.2880869e15*6.6256e-27/109678.758 DO IT=1,NTRANS IF(FR0(IT).GT.0.) THEN ALAM=2.997925D18/FR0(IT) ELSE ALAM=0. END IF if(indexp(it).ne.0) then ii=ilow(it) io=iel(ii) n1=nfirst(io)-1 WRITE(6,608) IT,ILOW(IT),IUP(IT), * typion(io),ii-n1,iup(it)-n1, * INDEXP(IT),ICOL(IT), * IFR0(IT),IFR1(IT),OSC0(IT),FR0(IT),ALAM itr=it if(ioptab.gt.0.and.line(itr)) then jj=iup(itr) iaa=numat(iatm(ii)) atn=float(iaa)+(float(iz(io))-un)*0.01 gf=log10(osc0(itr)*g(ii)) eel=-(enion(ii)-enion(n1+1))/cone xl=(g(ii)-1.)/2. eeu=-(enion(jj)-enion(n1+1))/cone xu=(g(jj)-1.)/2. wlam=alam if(wlam.gt.2000.) then ALM=1.E8/(alam*alam) xn1=64.328+29498.1/(146.-ALM)+255.4/(41.-ALM) wlam=alam/(xn1*1.d-6+un) end if write(66,666) alam,wlam,atn,gf,eel,xl,eeu,xu 666 format(2f11.3,f7.2,f7.3,f12.3,f6.1,f12.3,f6.1) end if end if END DO end if C C -------------------------------------------------------- C Evaluation of the photoionization cross-sections in both C explicit and fixed-option continuum transitions C C array CROSS(transition,frequency) C -------------------------------------------------------- C NFREQB=NFREQ IF(IBFINT.GT.0) NFREQB=NFREQC if(lasv) call sigave DO 220 ITR=1,NTRANS IF(LINE(ITR).OR.INDEXP(ITR).EQ.0) GO TO 220 IC=ITRCON(ITR) if(ibf(ic).gt.49.and.ibf(ic).lt.100) go to 220 ISIK=0 MODW=IABS(INDEXP(ITR)) IF(MODW.EQ.5 .OR. MODW.EQ.15) ISIK=1 DO IJ=1,NFREQB FR=FREQ(IJ) IF(ISPODF.GE.1) FR=FREQ(IFREQB(IJ)) CS=SIGK(FR,ITR,ISIK) BFCS(IC,IJ)=real(CS) IF(FR.LT.FR0PC(ITR)) BFCS(IC,IJ)=0. IF(IFWOP(ILOW(ITR)).LT.0) BFCS(IC,IJ)=1.E-30 END DO 220 CONTINUE C call rdatax(-1,ic,0) c call gomini C C ----------------------------------------- C Input parameters for additional opacities C ----------------------------------------- C IF(IOPADD.GT.0) THEN C C procedure OPADD0 now sets up array CROSS that C contain relevant cross-sections C DO IJ=1,NFREQB CALL OPADD0(IJ) END DO END IF C IF(IOPHL1.NE.0 .OR. IOPHL2.NE.0) then if(ielh.le.0) * CALL QUIT('ielh.le.0 for iophl1.gt.0',ielh,iophl1) CALL OPAHST end if if(iophe1.gt.0.and.ielhe1.le.0) * CALL QUIT('iophe1.gt.0.and.ielhe1.le.0',iophe1,ielhe1) if(iophe2.gt.0.and.ielhe2.le.0) * CALL QUIT('iophe2.gt.0.and.ielhe2.le.0',iophe2,ielhe2) if(iphe2c.gt.0.and.ielhe2.le.0) * CALL QUIT('iophe2c.gt.0.and.ielhe2.le.0',iophe2c,ielhe2) C C -------------------------------------------------------------------- C Parameters specifying the overall organization and structure of the C global matrices of complete linearization C -------------------------------------------------------------------- C if(hmix0.gt.0..and.iconv.eq.0) iconv=1 IF(HMIX0.GT.0.and.inpc.gt.0.and.inse.gt.0) THEN INDL = 3 INPC = 4 INSE = 5 INZD = 0 END IF c if(ifryb.gt.0) then inhe=0 inre=1 indl=0 inpc=0 inmp=0 inse=0 inzd=0 do ij=1,nfreq ijali(ij)=0 end do end if c nn=nfreqe if(inhe.gt.0) nn=nn+1 if(inre.gt.0) nn=nn+1 if(inpc.gt.0) nn=nn+1 if(inmp.gt.0) nn=nn+1 if(indl.gt.0) nn=nn+1 if(inzd.gt.0) nn=nn+1 nngg=nn if(inse.gt.0) nn=nn+NLVEXP LCHMAT=.FALSE. c C change LCHC and IRSPLT in case that INPC=0 C IF(INPC.EQ.0) THEN IRSPLT=1 LCHC=.TRUE. END IF C C -------------------------------------------------------------------- C Acceleration parameters C -------------------------------------------------------------------- C C if ITEK > 0, full iteration after Ng accelerations C < 0, only first ITEK iterations are full iterations; C KANT=0 - signifies a full complete linearization iteration; C KANT=1 - a Kantorovich iteration C DO IKT=1,NITER KANT(IKT)=0 END DO IF(ORELAX.EQ.0.) ORELAX=UN IACC0=IACC-3 LAC2=.FALSE. LRES2=.TRUE. IF(ITEK.GT.0) THEN DO IKT=ITEK+1,NITER KANT(IKT)=1 END DO DO IKT=IACC,18,IACD KANT(IKT)=0 END DO ELSE IF(ITEK.LT.0) THEN ITEK=ABS(ITEK) DO IKT=ITEK+1,NITER KANT(IKT)=1 END DO END IF C IF(KSNG.EQ.0) THEN DO I=1,NN LSNG(I)=.TRUE. END DO ELSE DO I=1,NN LSNG(I)=.FALSE. END DO END IF IF(KSNG.EQ.1) THEN DO I=1,NFREQE LSNG(I)=.TRUE. END DO IF(INSE.GT.0) THEN DO ION=1,NION N0=NFIRST(ION) N0G=NNGG+ABS(IIEXP(N0)) LSNG(N0G)=.TRUE. NKE=NNEXT(ION) NKG=NNGG+ABS(IIEXP(NKE)) LSNG(NKG)=.TRUE. END DO END IF END IF IF(INHE.GT.0) LSNG(NFREQE+INHE)=.TRUE. IF(INRE.GT.0) LSNG(NFREQE+INRE)=.TRUE. IF(INPC.GT.0) LSNG(NFREQE+INPC)=.TRUE. IF(INMP.GT.0) LSNG(NFREQE+INMP)=.TRUE. IF(INDL.GT.0) LSNG(NFREQE+INDL)=.TRUE. C C ------------------------------------ C Wind-blanketing parameters - albedos; or irradiation intensity C ------------------------------------ C IF(IWINBL.EQ.1) THEN READ(IBUFF,*) (ALBE(I),I=1,NFREQ) ELSE IF(IWINBL.EQ.2) THEN READ(IBUFF,*) (ALBE(I),I=1,NFREQ) ELSE IF(IWINBL.EQ.3) THEN READ(IBUFF,*) ALCON0,ALCHE2,ALLIN0 DO IJ=1,NFREQC ALBE(IJ)=ALCON0 IF(FREQ(IJ).GT.1.315D16) ALBE(IJ)=ALCHE2 END DO DO IJ=NFREQC+1,NFREQ ALBE(IJ)=ALLIN0 END DO READ(IBUFF,*) NUMMOD IF(NUMMOD.GT.0) THEN DO I=1,NUMMOD READ(IBUFF,*) IJ1,IJ2,ALB DO IJ=IJ1,IJ2 ALBE(IJ)=ALB END DO END DO END IF END IF C C ------------------------------------ C External irradiation intensity C ------------------------------------ C if(wdil.le.0..and.adist.gt.0.) * wdil=4.*(rstar/adist)**2 EXTOT=0. IF(TRAD.EQ.0) THEN DO IJ=1,NFREQ EXTRAD(IJ)=0. END DO ELSE IF(TRAD.GT.0) THEN DO IJ=1,NFREQ EXTRAD(IJ)=BNUE(IJ)/(EXP(HK*FREQ(IJ)/TRAD)-UN)*WDIL EXTOT=EXTOT+W(IJ)*EXTRAD(IJ) END DO ELSE open(48,file='stellarspectrum.dat',status='old') nfreq0=-ifix(real(trad)) DO IJ=NFREQ0,1,-1 read(48,*) fst(ij), EXTRD0(IJ) end do close(48) call interp(fst,extrd0,freq,extrad,nfreq0,nfreq,2,0,0) DO IJ=1,NFREQ EXTRAD(IJ)=EXTRAD(IJ)*WDIL EXTOT=EXTOT+W(IJ)*EXTRAD(IJ) END DO END IF DO IJ=1,NFREQ DO IMU=1,NMU EXTINT(IJ,IMU)=EXTRAD(IJ) END DO HEXTRD(IJ)=HALF*EXTRAD(IJ) END DO TSTAR=TRAD EXTOT=EXTOT/SIG4P/4. EXTOT0=WDIL*TRAD*TRAD*TRAD*TRAD write(6,496) extot0,extot 496 format(/' EXTERNAL IRRADIATION - EXTOT0, EXTOT: ',1p2e12.3/) C if(nd.gt.mdepth) CALL QUIT('nd.gt.mdepth',nd,mdepth) if(nfreq.gt.mfreq) CALL QUIT('nfreq.gt.mfreq',nfreq,mfreq) if(nfreqe.gt.mfrex) CALL QUIT('nfreqe.gt.mfrex',nfreqe,mfrex) if(nlevel.gt.mlevel) CALL QUIT('nlevel.gt.mlevel',nlevel,mlevel) if(ntrans.gt.mtrans) CALL QUIT('ntrans.gt.mtrans',ntrans,mtrans) if(natom.gt.matom) CALL QUIT('natom.gt.matom',natom,matom) if(nion.gt.mion) CALL QUIT('nion.gt.mion',nion,mion) if(nn.gt.mtot) CALL QUIT('nn.gt.mtot',nn,mtot) if(nlambd.gt.mlambd) CALL QUIT('nlambd.gt.mlambd',nlambd,mlambd) if(nhod.gt.mhod) CALL QUIT('nhod.gt.mhod',nhod,mhod) if(ntrans.gt.32767) * call quit(' Too many transitions to define ITRLIN as INTEGER*2', * ntrans,ntrans) C C ----------------------------------------------------------- C calculate the input LTE-grey model, if needed C ----------------------------------------------------------- C IF(LTGREY.AND.IDISK.EQ.0) CALL LTEGR IF(LTGREY.AND.IDISK.GT.0) CALL LTEGRD C CALL NSTOUT CALL DMDER C C 601 FORMAT(1H1,'*******************************'// * ' M O D E L A T M O S P H E R E'// * ' *******************************'// * ' TEFF =',F10.0/' LOG G =',F10.2/) c 602 FORMAT(1H0//' YTOT =',F11.5/ c * ' WMY =',1PD15.5/' WMM =',D15.5/) 602 FORMAT(' YTOT WMY WMM ',F11.5,1P2D15.5) 603 FORMAT(1H0//' EXPLICIT IONS INCLUDED'/ * ' ----------------------'// * ' NO. ION N0 N1 NK IZ IUPSUM ICUP FF'/) 604 FORMAT(1H ,I3,2X,A4,6I6,1PD15.3) 605 FORMAT(1H0//' EXPLICIT ENERGY LEVELS INCLUDED'/ * ' -------------------------------'// * ' NO. LEVEL ION ION.FREQ.(s^-1) G NQ', * ' IEL ILK IAT IMOD ILT IIE IIF'/) 606 FORMAT(1H ,I4,1X,A10,A4,1PD15.7,0PF10.2,8I5) 607 FORMAT(1H0//' TRANSITION PARAMETERS (TOTAL OF',I6,' )'/ * ' ---------------------'// * ' ITR ILOW IUP',14x,' INDEXP ICOL IFR0 IFR1 OSC', * ' FR0',8X,' LAMBDA'/) 608 FORMAT(1H ,I6,2I5,a6,2i4,2i5,2I7,1P2D12.3,0PF12.3) c 609 FORMAT(1H0//' FREQUENCY POINTS AND WEIGHTS - EXPLICIT'/ c * ' ---------------------------------------'// c * ' IJ ',7X,'FREQ',13X,'WEIGHT',11X,'PROF'/) c 610 FORMAT(1H ,I7,1P2D17.8,D15.5,I5,D17.8) C RETURN END