SUBROUTINE HYDLIN(ID,I0,I1,ABSOH,EMISH) C ======================================= C C opacity and emissivity of hydrogen lines C INCLUDE 'PARAMS.FOR' INCLUDE 'MODELP.FOR' INCLUDE 'SYNTHP.FOR' PARAMETER (FRH1=3.28805E15,FRH2=FRH1/4.,UN=1.,SIXTH=1./6.) PARAMETER (CPP=4.1412E-16,CPJ=157803.) PARAMETER (C00=1.25E-9,CDOP=1.284523E12,CID=0.02654,TWO=2.) PARAMETER (CPJ4=CPJ/4.,AL10=2.3025851,CINV=UN/2.997925E18) PARAMETER (CID1=0.01497) common/quasun/nunalp,nunbet,nungam,nunbal common/hhebrd/sthe,nunhhe common/gompar/hglim,ihgom DIMENSION PJ(40),PRF0(54),OSCH(4,22), * ABSO(MFREQ),EMIS(MFREQ),ABSOH(MFREQ),EMISH(MFREQ) dimension wlir(15),irlow(15),irupp(15) DATA FRH /3.289017E15/ data wlir/ * 123680., 75005., 59066., 51273.,190570.,113060., * 87577., 75061.,277960.,162050.,123840.,105010., * 223340.,168760.,141790./ data irlow/4*6, 4*7, 4*8, 3*9/ data irupp/7,8,9,10,8,9,10,11,9,10,11,12,11,12,13/ data nlinir/15/ c DATA INIT /0/ C DO IJ=I0,I1 ABSOH(IJ)=0. EMISH(IJ)=0. END DO c if(iath.le.0.or.rrr(1,1,1).eq.0.) return izz=1 C IF(INIT.EQ.0) THEN DO I=1,4 DO J=I+1,22 CALL STARK0(I,J,IZZ,XK,WL0,FIJ,FIJ0) WLINE(I,J)=WL0 OSCH(I,J)=FIJ+FIJ0 END DO END DO INIT=1 END IF DO IJ=I0,I1 ABSO(IJ)=0. EMIS(IJ)=0. END DO c if(ilowh.le.0) return c T=TEMP(ID) T1=UN/T SQT=SQRT(T) ANE=ELEC(ID) ANES=EXP(SIXTH*LOG(ANE)) TL=LOG10(T) ANEL=LOG10(ANE) C C populations of the first 40 levels of hydrogen C ANP=POPUL(NKH,ID) PP=CPP*ANE*ANP*T1/SQT NLH=N1H-N0HN+1 c if(ifwop(n1h).lt.0) nlh=nlh-1 nlh=nlh-1 DO IL=1,50 X=IL*IL IF(IL.LE.NLH) PJ(IL)=POPUL(N0HN+IL-1,ID) IF(IL.GT.NLH) PJ(IL)=PP*EXP(CPJ/X*T1)*X*wnhint(il,id) END DO p2=pp*exp(cpj4*t1)*4.*wnhint(2,id) c C Frequency- and line-independent parameters for evaluating the C asymptotic Stark profile C F00=C00*ANES*ANES*ANES*ANES DOP0=1.E8*SQRT(1.65E8*T+VTURB(ID)) C C ------------------------------------------------------------------- C overall loop over spectral series (only in the infrared region) C ------------------------------------------------------------------- C ISERL=ILOWH ISERU=ILOWH c if(wlam(i0).gt.14000.) iseru=4 if(wlam(i0).gt.22700.) iseru=5 if(wlam(i0).gt.32800.) iseru=6 if(wlam(i0).gt.44660.) iseru=7 if(wlam(i0).gt.60000.) iserl=4 c if(iserl.eq.3.and.iseru.eq.3.and.nunbal.gt.0) iserl=2 DO IJ=I0,I1 ABSO(IJ)=0. EMIS(IJ)=0. END DO C c ======================== c loop over spectral series c ======================== c DO I=ISERL,ISERU c c skip the following calculations if one uses the Gomez tables c if(ihgom.gt.0.and.elec(id).gt.hglim) then if(i.ge.1.and.i.le.ihgom) then call ghydop(id,i0,i1,pj,absoh,emish) go to 200 end if end if c II=I*I XII=UN/II POPI=PJ(I) IF(I.EQ.1) FRH=3.28805E15 C C determination of which hydrogen lines contribute in a current C frequency region C M1=M10 IF(I.LT.ILOWH) M1=ILOWH-1 M2=M1+1 M1=M1-1 M2=M20+3 IF(M1.LT.I+1) M1=I+1 if(grav.gt.3.) then m2=m2+5 m1=m1-3 if(m1.gt.i+6) m1=m1-3 end if c new! if(i.ge.3) then m1=i+1 m2=i+40 end if if(i.ge.4) m2=i+20 if(i.ge.6) m2=i+10 C C loop over lines which contribute at given wavelength region C m1=min(m1,40) m2=min(m2,40) m1=max(m1,i+1) m2=max(m2,i+2) DO J=M1,M2 ILINE=0 JJ=J*J XJJ=UN/JJ ABTRA=PJ(I)*WNHINT(J,ID) EMTRA=PJ(J)*WNHINT(I,ID)*II*XJJ*EXP(CPJ*(XII-XJJ)*T1) if(i.le.2.and.j.le.i+2) then abtra=pj(i) emtra=pj(j)*wnhint(i,id)/wnhint(j,id)* * ii*xjj*exp(cpj*(xii-xjj)*t1) end if IF(I.LE.4.AND.J.LE.22) ILINE=ILIN0(I,J) c c quasi-molecular opacity for Lyman-alpha and beta satellites c lquasi=i.eq.1.and.j.eq.2.and.nunalp.gt.0 lquasi=lquasi.or.i.eq.1.and.j.eq.3.and.nunbet.gt.0 lquasi=lquasi.or.i.eq.1.and.j.eq.4.and.nungam.gt.0 lquasi=lquasi.or.i.eq.2.and.j.eq.3.and.nunbal.gt.0 lalhhe=i.eq.1.and.j.eq.2.and.nunhhe.gt.0 if(lquasi) then DO IJ=I0,I1 call allard(wlam(ij),popi,anp,sg,i,j) ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA END DO end if ahe=0. if(iathe.gt.0) ahe=popul(n0a(iathe),id) if(lalhhe.and.ahe.gt.0.) then rel=1./6.2831855 do ij=i0,i1 call lyahhe(wlam(ij),ahe,sg0) sg=sg0*rel abso(ij)=abso(ij)+sg*abtra emis(ij)=emis(ij)+sg*emtra end do end if c c lines with special Stark broadening tables c IF(ILINE.GT.0) THEN FID=CID*OSCH(I,J) c c switch to either original Lemke/Tremblay of Xenomorph c if(ilxen(i,j).eq.0.or.anel.lt.xnemin) then c c original Lemke/Tremblay c NWL=NWLHYD(ILINE) DO IWL=1,NWL PRF0(IWL)=PRFHYD(ILINE,ID,IWL) END DO DO IJ=I0,I1 AL=ABS(WLAM(IJ)-WLINE(I,J)) IF(AL.LT.1.E-4) AL=1.E-4 IF(ILEMKE.EQ.1) AL=AL/F00 AL=LOG10(AL) DO 30 IWL=1,NWL-1 IW0=IWL IF(AL.LE.WLHYD(ILINE,IWL+1)) GO TO 40 30 CONTINUE 40 IW1=IW0+1 PRFF=(PRF0(IW0)*(WLHYD(ILINE,IW1)-AL)+PRF0(IW1)* * (AL-WLHYD(ILINE,IW0)))/ * (WLHYD(ILINE,IW1)-WLHYD(ILINE,IW0)) SG=EXP(PRFF*AL10)*FID sg0=EXP(PRFF*AL10) IF(ILEMKE.EQ.1) SG=SG*WLINE(I,J)**2*CINV/F00 ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA END DO c c XENOMORPH data for selected lines c else ixn=ilxen(i,j) nwl=nwlxen(ixn) fr0l=2.997925e18/wline(i,j) do ij=i0,i1 al=(freq(ij)-fr0l)/f00 if(abs(al).lt.1.e-4) al=1.e-4 all=log10(abs(al)) do 51 iwl=1,nwl-1 iw0=iwl if(all.le.alxen(ixn,iwl+1)) go to 52 51 continue 52 iw1=iw0+1 if(al.gt.0.) then prff=(prfb(ixn,id,iw0)*(alxen(ixn,iw1)-all)+ * prfb(ixn,id,iw1)*(all-alxen(ixn,iw0)))/ * (alxen(ixn,iw1)-alxen(ixn,iw0)) else prff=(prfr(ixn,id,iw0)*(alxen(ixn,iw1)-all)+ * prfr(ixn,id,iw1)*(all-alxen(ixn,iw0)))/ * (alxen(ixn,iw1)-alxen(ixn,iw0)) end if sg=exp(prff*al10)*fid/f00 ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA end do END IF c c lines without special Stark broadening tables c ELSE CALL STARK0(I,J,izz,XKIJ,WL0,FIJ,FIJ0) if((wl0.le.wlam(i1).and.1.25*wl0.gt.wlam(i0)). or. * (wl0.ge.wlam(i0).and.0.75*wl0.lt.wlam(i1))) then FXK=F00*XKIJ FXK1=UN/FXK DOP=DOP0/WL0 DBETA=WL0*WL0*CINV*FXK1 BETAD=DOP*DBETA FID=CID*FIJ*DBETA c FID0=CID1*FIJ0/DOP CALL DIVSTR(AD,DIV) fac=two if(lquasi) fac=un DO IJ=I0,I1 fr=freq(ij) BETA=ABS(WLAM(IJ)-WL0)*FXK1 IF(I.LT.5) THEN SG=STARKA(BETA,AD,DIV,fac)*FID if(iophli.eq.2.and.i.eq.1.and.j.eq.2) * sg=sg*feautr(fr,id) ELSE SG=STARKIR(II,JJ,T,ANE,BETA)*FID END IF ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA END DO END IF END IF END DO END DO C C far infrared hydrogen lines C if(wlam(i1).gt.70000.) then DO I=8,13 II=I*I XII=UN/II DO J=I+1,I+4 JJ=J*J XJJ=UN/JJ CALL STARK0(I,J,izz,XKIJ,WL0,FIJ,FIJ0) if((wl0.le.wlam(i1).and.1.5*wl0.gt.wlam(i0)). or. * (wl0.ge.wlam(i0).and.0.5*wl0.lt.wlam(i1))) then FXK=F00*XKIJ FXK1=UN/FXK DOP=DOP0/WL0 DBETA=WL0*WL0*CINV*FXK1 BETAD=DOP*DBETA FID=CID*FIJ*DBETA CALL DIVSTR(AD,DIV) fac=two DO IJ=I0,I1 fr=freq(ij) BETA=ABS(WLAM(IJ)-WL0)*FXK1 SG=STARKIR(II,JJ,T,ANE,BETA)*FID ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA END DO END IF END DO END DO END IF 200 continue c if(wlam(i1).gt.5.e5) then do ij=i0,i1 fr=freq(ij) do ilir=1,nlinir if(wlam(ij).gt.wlir(ilir)*0.95.and. * wlam(ij).lt.wlir(ilir)*1.05) then j=irupp(ilir) JJ=J*J i=irlow(ilir) II=I*I XII=UN/II XJJ=UN/JJ ABTRA=PJ(I)*WNHINT(J,ID) EMTRA=PJ(J)*WNHINT(I,ID)*II*XJJ*EXP(CPJ*(XII-XJJ)*T1) CALL STARK0(I,J,izz,XKIJ,WL0,FIJ,FIJ0) FXK=F00*XKIJ FXK1=UN/FXK DOP=DOP0/WL0 DBETA=WL0*WL0*CINV*FXK1 BETAD=DOP*DBETA FID=CID*FIJ*DBETA CALL DIVSTR(AD,DIV) fac=two BETA=ABS(WLAM(IJ)-WL0)*FXK1 SG=STARKA(BETA,AD,DIV,fac)*FID ABSO(IJ)=ABSO(IJ)+SG*ABTRA EMIS(IJ)=EMIS(IJ)+SG*EMTRA end if end do end do end if C C ---------------------------- C total opacity and emissivity C ---------------------------- C DO IJ=I0,I1 F=FREQ(IJ) F15=F*1.E-15 XKF=EXP(-4.79928e-11*F*T1) XKFB=XKF*1.4743E-2*F15*F15*F15 ABSOH(IJ)=ABSO(IJ)-XKF*EMIS(IJ) EMISH(IJ)=XKFB*EMIS(IJ) END DO RETURN END