diff options
Diffstat (limited to 'synthe/synthe.for')
| -rw-r--r-- | synthe/synthe.for | 2978 |
1 files changed, 2978 insertions, 0 deletions
diff --git a/synthe/synthe.for b/synthe/synthe.for new file mode 100644 index 0000000..73a15bc --- /dev/null +++ b/synthe/synthe.for @@ -0,0 +1,2978 @@ + PROGRAM SYNTHE +c revised 11nov14 modified to agree with Fiorella Castelli's Linux version +c including He I lines and reverting to her HPROF +c revised 4nov14 constants given D exponents +c revised 5mar02 Jo Bruls found dimensions of CONTX too small, now 26 +c revised 20jul01 John Lester found dimensions of CONTX too small, now 25 +c revised 14nov99 +C THIS PROGRAM IS REAL*4 EXCEPT FOR WAVELENGTHS AND ENERGY LEVELS +C REAL*8 INPUT ON TAPE10 IS CONVERTED TO REAL*4 +C +C TAPE5 INPUT +C TAPE6 OUTPUT +C TAPE7 temporary file for line identifications +C TAPE8 file that passes the number of lines to SPECTRV +C TAPE9 output opacity vectors and line data +C TAPE10 input from XNFPELSYN +C TAPE12 input line data needed for calculation +C TAPE13 input line data, all data for each line +C TAPE14 temporary file of opacity spectra +C TAPE15 temporary file of opacity vectors for each line +C TAPE19 input line data from RNLTE +C TAPE93 parameters for this run from SYNBEG + PARAMETER (kw=99) + PARAMETER (LENREC=8000,MAXLEN=2000001,MAXPROF=10000, + 1 MAXBUFF=MAXLEN+MAXPROF,MAXLIN=MAXBUFF+MAXPROF*2) +C LENREC transposition is done in blocks of LENRECxkw +C MAXLEN number of points in the spectrum +C MAXPROF number of points in either wing of a line with a Voigt profile + COMMON /BUFFER/BUFFER(MAXBUFF),PROFILE(MAXPROF) + DIMENSION LINE(MAXLIN) + DIMENSION TRANSP(kw,LENREC),RECORD(LENREC),INDEXR(5000) + EQUIVALENCE (BUFFER(1),LINE(1),TRANSP(1)) +C DIMENSION TRANSP(55000,kw) + COMMON /CONTIN/CONTINUUM(MAXBUFF) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) + COMMON /EXTAB/EXTAB(1001),EXTABF(1001),E1TAB(2000) + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001) + COMMON /LINDAT/WL,E,EP,LABEL(2),LABELP(2),OTHER1(2),OTHER2(2), + 1 WLVAC,CENTER,CONCEN, NELION,GAMMAR,GAMMAS,GAMMAW,REF, + 2 NBLO,NBUP,ISO1,X1,ISO2,X2,GFLOG,XJ,XJP,CODE,ELO,GF,GS,GR,GW, + 3 DWL,DGFLOG,DGAMMAR,DGAMMAS,DGAMMAW,EXTRA1,EXTRA2,EXTRA3 + COMMON /NLINES/WLBEG,WLEND,RESOLU,RATIO,RATIOLG,WBEGIN, + 1 LENGTH,MLINES,IXWLBEG + COMMON /RHOX/RHOX(kw),NRHOX + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND, + 1 IFTURB + COMMON /TXNXN/EMERGE(kw),TXNXN(kw),BSTIM(kw),XNFH2(kw) + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + common /elem/abund(99),atmass(99),elem(99) + DIMENSION ABLOG(3,377),IFTP(kw),ABMIN(kw),MLINEJ(kw) + DIMENSION ASYNTH(kw),ALINEC(kw),ASYNCONT(kw),ALINECONT(kw) + DIMENSION DECKJ(7,kw) + DIMENSION VELSHIFT(kw),HFIELD(kw) + DIMENSION TITLE(74),IDMOL(60),MOMASS(60) + REAL*8 FRQEDG(377),WLEDGE(377),CMEDGE(377),CONFRQ(1131) + REAL*8 DELEDGE(377),HALFEDGE(377) + REAL KAPPA0,KAPPA,KAPCEN,KAPMIN + REAL*8 QT(kw),QTKEV(kw),QTK(kw),QHKT(kw),QTLOG(kw),QHCKT(kw) + REAL*8 QP(kw),QXNE(kw),QXNATOM(kw),QRHO(kw),QRHOX(kw),QVTURB(kw) + REAL*8 QXNFH(kw),QXNFHE(kw,2),QXNFH2(kw),QDOPPLE(594),QXNFPEL(594) + REAL*8 QABLOG(1131),QXNFDOP(594),QCONGF +C REAL*8 ASYNTH,ALINEC,TITLE,TEFF,GLOG,IDMOL,MOMASS + REAL*8 TITLE,TEFF,GLOG,IDMOL,MOMASS + dimension hfactor(kw),hefactor(kw),h2factor(kw) + REAL*8 CONWAV,CONCM,WAVE + REAL*8 LINDAT8(14) + REAL*4 LINDAT4(28) + EQUIVALENCE (LINDAT8(1),WL),(LINDAT4(1),NELION) + REAL*8 RESOLU,RATIO,RATIOLG,SIGMA2,WLBEG,WLEND,WBEGIN + REAL*8 EMERGE + REAL*8 WL,E,EP,WLVAC,CENTER,CONCEN + REAL*8 LABEL,LABELP,OTHER1,OTHER2 + real*8 wavel +C + FASTEX(X)=EXTAB(IFIX(X)+1)* + 1EXTABF(IFIX((X-FLOAT(IFIX(X)))*1000.+1.5)) +C + READ(93)NLINES,LENGTH,IFVAC,IFNLTE,N19,TURBV,DECKJ,IFPRED, + 1WLBEG,WLEND,RESOLU,RATIO,RATIOLG,CUTOFF,LINOUT + IXWLBEG=DLOG(WLBEG)/RATIOLG + WBEGIN=DEXP(IXWLBEG*RATIOLG) + IF(WBEGIN.LT.WLBEG)THEN + IXWLBEG=IXWLBEG+1 + WBEGIN=DEXP(IXWLBEG*RATIOLG) + ENDIF + CLOSE(UNIT=93,DISP='DELETE') + OPEN(UNIT=10,STATUS='OLD',FORM='UNFORMATTED',READONLY,SHARED) + OPEN(UNIT=12,STATUS='OLD',FORM='UNFORMATTED',ACCESS='APPEND') + OPEN(UNIT=13,STATUS='NEW',FORM='UNFORMATTED') + OPEN(UNIT=14,STATUS='OLD',FORM='UNFORMATTED',ACCESS='APPEND') + OPEN(UNIT=19,STATUS='OLD',FORM='UNFORMATTED',ACCESS='APPEND') + OPEN(UNIT=20,STATUS='OLD',FORM='UNFORMATTED',ACCESS='APPEND') + IF(LINOUT.LT.0.)GO TO 3442 + IF(N19.GT.0)THEN + REWIND 20 + DO 3440 I=1,N19 + READ(20)LINDAT8,LINDAT4 + 3440 WRITE(13)LINDAT8,LINDAT4 + ENDIF + IF(NLINES.GT.0)THEN + REWIND 14 + DO 3441 I=1,NLINES + READ(14)LINDAT8,LINDAT4 + 3441 WRITE(13)LINDAT8,LINDAT4 + ENDIF + 3442 CLOSE(UNIT=20,DISP='DELETE') + CLOSE(UNIT=14,DISP='DELETE') + OPEN(UNIT=14,TYPE='NEW',FORM='UNFORMATTED',RECORDTYPE='FIXED', + 1ACCESS='DIRECT',BLOCKSIZE=LENREC*4,RECORDSIZE=LENREC) +C + DO 3456 I=1,1001 + EXTAB(I)=EXP(-FLOAT(I-1)) + 3456 EXTABF(I)=EXP(-FLOAT(I-1)*.001) + DO 3457 I=1,2000 + 3457 E1TAB(I)=EXPI(1,FLOAT(I)*.01) +C PRETABULATE VOIGT FUNCTION +C 200 STEPS PER DOPPLER WIDTH + VSTEPS=200. + CALL TABVOIGT(VSTEPS,2001) +C + READ(10)NRHOX,TEFF,GLOG,TITLE + WRITE(6,2000)TEFF,GLOG,TITLE + 2000 FORMAT(6H TEFF=,F10.1,3X,5HGRAV=,F6.3,3X,74A1) + READ(10)NEDGE,(FRQEDG(IEDGE),WLEDGE(IEDGE),CMEDGE(IEDGE), + 1IEDGE=1,NEDGE),IDMOL,MOMASS + READ(10)NCON,(CONFRQ(NU),NU=1,NCON) + WLEDGE(1)=ABS(WLEDGE(1)) + DO 2001 IEDGE=2,NEDGE + WLEDGE(IEDGE)=ABS(WLEDGE(IEDGE)) + HALFEDGE(IEDGE-1)=(WLEDGE(IEDGE-1)+WLEDGE(IEDGE))*.5 + 2001 DELEDGE(IEDGE-1)=(WLEDGE(IEDGE)-WLEDGE(IEDGE-1))**2*.5 +C + ITEMP=1 +C READ(10)T,TKEV,TK,HKT,TLOG,HCKT,P,XNE,XNATOM,RHO,RHOX,VTURB, +C 1XNFH,XNFHE,XNFH2 +C ON VAX VARIABLES ARE READ REAL*8 AND CONVERTED TO REAL*4 + READ(10)QT,QTKEV,QTK,QHKT,QTLOG,QHCKT,QP,QXNE,QXNATOM,QRHO, + 1QRHOX,QVTURB,QXNFH,QXNFHE,QXNFH2 + DO 2010 J=1,NRHOX + T(J)=QT(J) + TKEV(J)=QTKEV(J) + TK(J)=QTK(J) + HKT(J)=QHKT(J) + TLOG(J)=QTLOG(J) + HCKT(J)=QHCKT(J) + P(J)=QP(J) + XNE(J)=QXNE(J) + XNATOM(J)=QXNATOM(J) + RHO(J)=QRHO(J) + RHOX(J)=QRHOX(J) + VTURB(J)=QVTURB(J) + XNFH(J)=QXNFH(J) + XNFHE(J,1)=QXNFHE(J,1) + XNFHE(J,2)=QXNFHE(J,2) + 2010 XNFH2(J)=QXNFH2(J) +C + DO 2011 J=1,kw + VELSHIFT(J)=DECKJ(1,J) + HFIELD(J)=DECKJ(2,J) + 2011 CONTINUE +C + WRITE(9)WLBEG,RESOLU,WLEND,LENGTH,NRHOX,LINOUT,TURBV,IFVAC + WRITE(9)NEDGE,(FRQEDG(IEDGE),WLEDGE(IEDGE),CMEDGE(IEDGE), + 1IEDGE=1,NEDGE),IDMOL,MOMASS +C + ILINES=0 + N12=NLINES + NLINES=NLINES+N19 + IREC=0 + DO 500 J=1,NRHOX + REWIND 12 +C INITIALIZE BUFFER + DO 210 NBUFF=1,LENGTH + 210 BUFFER(NBUFF)=0. + READ(10)QABLOG + READ(10) + READ(10) + NU=0 + DO 2002 IEDGE=1,NEDGE-1 + NU=NU+1 + ABLOG(1,IEDGE)=QABLOG(NU) + NU=NU+1 + ABLOG(2,IEDGE)=QABLOG(NU) + NU=NU+1 + 2002 ABLOG(3,IEDGE)=QABLOG(NU) + IEDGE=1 + DO 2005 NBUFF=1,LENGTH + WAVE=WBEGIN*RATIO**(NBUFF-1) + 2003 IF(WAVE.LT.WLEDGE(IEDGE+1))GO TO 2005 + IEDGE=IEDGE+1 + GO TO 2003 + 2005 CONTINUUM(NBUFF)= + 1((WAVE-HALFEDGE(IEDGE))*(WAVE-WLEDGE(IEDGE+1))*ABLOG(1,IEDGE)+ + 2 (WLEDGE(IEDGE)-WAVE)*(WAVE-WLEDGE(IEDGE+1))*2.*ABLOG(2,IEDGE)+ + 3 (WAVE-WLEDGE(IEDGE))*(WAVE-HALFEDGE(IEDGE))*ABLOG(3,IEDGE))/ + 4 DELEDGE(IEDGE) + DO 2006 NBUFF=1,LENGTH + 2006 CONTINUUM(NBUFF)=10.**CONTINUUM(NBUFF) +C + READ(10)QXNFPEL,QDOPPLE + WRITE(28,2828)J,QXNFPEL + 2828 FORMAT(I5/(1P6E21.14)) +C IF(IFTP(J).EQ.0)GO TO 400 + XNFPH(J,1)=QXNFPEL(1) + XNFPH(J,2)=QXNFPEL(2) + XNFPHE(J,1)=QXNFPEL(7) + XNFPHE(J,2)=QXNFPEL(8) + XNFPHE(J,3)=QXNFPEL(9) + DO 203 NELION=1,594 + XNFPEL(NELION)=0. +C PATCH FOR NLTE HELIUM WHERE CAN GET OVERFLOW ON VAX +C NLTE HELIUM IS NOT COMPUTED BY SYNTHE IN ANY CASE + IF(QXNFPEL(NELION).LT.1.D25)XNFPEL(NELION)=QXNFPEL(NELION) + 203 CONTINUE + DO 205 NELION=1,594 + QDOPPLE(NELION)=SQRT(QDOPPLE(NELION)**2+(TURBV/299792.458D0)**2) + DOPPLE(NELION)=QDOPPLE(NELION) +C PROBLEMS WITH OVERFLOW + XNFPEL(NELION)=XNFPEL(NELION)/RHO(J) + QXNFDOP(NELION)=QXNFPEL(NELION)/QRHO(J)/QDOPPLE(NELION) + 205 XNFDOP(NELION)=QXNFDOP(NELION) + TXNXN(J)=(XNFH(J)+.42D0*XNFHE(J,1)+.85D0*XNFH2(J))* + 1(T(J)/10000.D0)**.3 +C +C DOPPLER SHIFT IN POINT NUMBERS + NVSHIFT=RESOLU*VELSHIFT(J)/299792.458D0+.5D0 + WRITE(6,215)J,VELSHIFT(J),NVSHIFT + 215 FORMAT(I5,' VELOCITY SHIFT',F9.3,I7) +C ADD LINES TO BUFFER + MLINES=0 + IF(N19.GT.0)CALL XLINOP(J,N19,CUTOFF,VELSHIFT(J),IFVAC,LINOUT) + IF(N12.EQ.0)GO TO 400 + N191=N19+1 + alpha=0. + DO 350 ILINE=N191,NLINES + READ(12)NBUFF,CONGF,NELION,ELO,GAMRF,GAMSF,GAMWF +c +c include Barklem, Anstee, and O'Mara van der Waals +c READ(12)NBUFF,CONGF,NELION,ELO,GAMRF,GAMSF,GAMWF,alpha +c + QCONGF=CONGF + KAPPA0=CONGF*QXNFDOP(NELION) +C KAPPA0=CONGF*XNFDOP(NELION) +C PROBLEMS WITH OVERFLOW ON VAX +C KAPPA0=CONGF/DOPPLE(NELION)*XNFPEL(NELION) + KAPMIN=CONTINUUM(MIN(MAX(NBUFF,1),LENGTH))*CUTOFF + IF(KAPPA0.LT.KAPMIN)GO TO 350 + KAPPA0=KAPPA0*FASTEX(ELO*HCKT(J)) +C KAPPA0=KAPPA0*EXP(-ELO*HCKT(J)) + IF(KAPPA0.LT.KAPMIN)GO TO 350 +C +c Castelli for Barklem, Anstee, and O'mara van der Waals broadening + if(alpha.ne.0.)then + nelem=int(nelion/6)+1 + v2=(1.-alpha)/2 + wavel=wbegin*ratio**(nbuff-1) + type*,wavel,nelem,alpha,v2 + hfactor(j)=(t(j)/10000.)**v2 + hefactor(j)=0.628*(2.0991D-4*T(j)*(1/4+1.008/atmass(nelem)))**v2 + h2factor(j)=1.08*(2.0991D-4*T(j)*(1/2+1.008/atmass(nelem)))**v2 + txnxn(j)=xnfh(j)*hfactor(j)+xnfhe(J,1)*hefactor(j)+xnfh2(j)*h2factor(j) + endif +c +C VOIGT APPROXIMATION ACCURATE ONLY TO ADAMP**2 + ADAMP=(GAMRF+GAMSF*XNE(J)+GAMWF*TXNXN(J))/DOPPLE(NELION) + NBUFF=NBUFF+NVSHIFT + IF(NBUFF.LT.1.OR.NBUFF.GT.LENGTH)GO TO 320 + MLINES=MLINES+1 + IF(ADAMP.LT..2)THEN + KAPCEN=KAPPA0*(1.-1.128*ADAMP) + ELSE + KAPCEN=KAPPA0*VOIGT(0.,ADAMP) + ENDIF + IF(LINOUT.GE.0)WRITE(15)ILINE,KAPCEN + BUFFER(NBUFF)=BUFFER(NBUFF)+KAPCEN +C PROFILE INSIDE 10 DOPPLER WIDTHS + 320 N10DOP=10.*(DOPPLE(NELION)*RESOLU) + IF(ADAMP.LT..2)THEN + TABSTEP=VSTEPS/(DOPPLE(NELION)*RESOLU) + TABI=1.5 + DO 321 NSTEP=1,N10DOP + TABI=TABI+TABSTEP + PROFILE(NSTEP)=KAPPA0*(H0TAB(IFIX(TABI))+ADAMP*H1TAB(IFIX(TABI))) + IF(PROFILE(NSTEP).LT.KAPMIN)GO TO 323 + 321 CONTINUE + ELSE + DVOIGT=1./DOPPLE(NELION)/RESOLU + DO 1321 NSTEP=1,N10DOP + PROFILE(NSTEP)=KAPPA0*VOIGT(FLOAT(NSTEP)*DVOIGT,ADAMP) + IF(PROFILE(NSTEP).LT.KAPMIN)GO TO 323 + 1321 CONTINUE + ENDIF +C FAR WINGS + X=PROFILE(N10DOP)*FLOAT(N10DOP)**2 + MAXSTEP=SQRT(X/KAPMIN)+1. + MAXSTEP=MIN(MAXSTEP,MAXPROF) + N1=N10DOP+1 + DO 322 NSTEP=N1,MAXSTEP + 322 PROFILE(NSTEP)=X/FLOAT(NSTEP)**2 + NSTEP=MAXSTEP + 323 IF(NBUFF+NSTEP.LT.1.OR.NBUFF-NSTEP.GT.LENGTH)GO TO 350 + IF(NBUFF.GE.LENGTH)GO TO 325 +C RED WING + MAXRED=MIN0(LENGTH-NBUFF,NSTEP) + MINRED=MAX0(1,1-NBUFF) + DO 324 ISTEP=MINRED,MAXRED + 324 BUFFER(NBUFF+ISTEP)=BUFFER(NBUFF+ISTEP)+PROFILE(ISTEP) + IF(NBUFF.LE.1)GO TO 350 +C BLUE WING + 325 MAXBLUE=MIN0(NBUFF-1,NSTEP) + MINBLUE=MAX0(1,NBUFF-LENGTH) + DO 326 ISTEP=MINBLUE,MAXBLUE + 326 BUFFER(NBUFF-ISTEP)=BUFFER(NBUFF-ISTEP)+PROFILE(ISTEP) + 350 CONTINUE +C + 400 CONTINUE +C MEMORY TRANSPOSITION +C DO 26 I=1,LENGTH +C 26 TRANSP(I,J)=BUFFER(I) +C DIRECT IO TRANSPOSITION + NUMREC=(LENGTH+LENREC-1)/LENREC + DO 26 NBEG=1,LENGTH,LENREC + IREC=IREC+1 + NEND=NBEG+LENREC-1 + WRITE(14,REC=IREC)(BUFFER(I),I=NBEG,NEND) + 26 CONTINUE + WRITE(6,499)J,MLINES + 499 FORMAT(30X,2I10,' LINES USED') + MLINEJ(J)=MLINES + ILINES=ILINES+MLINES + 500 CONTINUE + WRITE(6,106)ILINES + 106 FORMAT(I10) +C +C DIRECT IO TRANSPOSITION + N9=0 + DO 99 N=1,NUMREC + WRITE(6,106)N + DO 93 J=1,NRHOX + IREC=J*NUMREC-NUMREC+N + READ(14,REC=IREC)(TRANSP(J,I),I=1,LENREC) + 93 CONTINUE + NOUT=LENREC + LASTREC=LENGTH-NUMREC*LENREC+LENREC + IF(N.EQ.NUMREC)NOUT=LASTREC + DO 95 I=1,NOUT + N9=N9+1 + WAVE=WBEGIN*RATIO**(N9-1) + FREQ=2.99792458D17/WAVE + DO 94 J=1,NRHOX + 94 ASYNTH(J)=TRANSP(J,I)*(1.-EXP(-FREQ*HKT(J))) + WRITE(29,2929)WAVE,(ASYNTH(J),J=1,NRHOX) + 2929 FORMAT(F15.8/(1P8E15.7)) + 95 WRITE(9)ASYNTH + 99 CONTINUE +C N9=LENGTH IS A CHECK THAT THE TRANSPOSITION WORKED PROPERLY + WRITE(6,96)N9 + 96 FORMAT(I10,27H OPACITY VECTORS ON TAPE 9 ) +C +C SAVE ALL LINES USED + N9=0 + IF(NLINES.EQ.0)GO TO 810 + IF(LINOUT.LT.0)GO TO 810 + DO 800 I=1,NLINES + 800 LINE(I)=0 + IF(ILINES.EQ.0)GO TO 810 + REWIND 15 + I=0 + 803 I=I+1 + READ(15)ILINE + LINE(ILINE)=1 + IF(I.LT.ILINES)GO TO 803 + REWIND 13 + REWIND 7 + DO 809 I=1,NLINES + READ(13)LINDAT8,LINDAT4 + IF(LINE(I).EQ.0)GO TO 809 + LINE(I)=0 + WRITE(7)LINDAT8,LINDAT4 + N9=N9+1 + LINE(I)=N9 + 809 CONTINUE + 810 WRITE(6,106)N9 + WRITE(8)N9 + WRITE(9)N9 + CLOSE(UNIT=12,DISP='DELETE') + CLOSE(UNIT=13,DISP='DELETE') +C +C SET UP LINE CENTER OPACITY FOR EACH LINE + IF(ILINES.EQ.0)CALL EXIT + IF(LINOUT.LT.0)CALL EXIT + REWIND 15 + IREC=0 + DO 830 J=1,NRHOX + MAXLINE=MLINEJ(J) +C WRITE(6,817)J,MAXLINE + 817 FORMAT(3I10) + K=0 + JOUT=0 + IF(MAXLINE.EQ.0)GO TO 820 + I9LAST=0 + DO 815 L=1,MAXLINE + READ(15)ILINE,KAPCEN + I9=LINE(ILINE) + IF(I9.EQ.0)GO TO 815 + NSKIP=I9-I9LAST-1 + IF(NSKIP.EQ.0)GO TO 811 + DO 812 ISKIP=1,NSKIP + K=K+1 + RECORD(K)=0. + JOUT=JOUT+1 + IF(K.LT.LENREC)GO TO 812 + IREC=IREC+1 + WRITE(14,REC=IREC)RECORD + K=0 + 812 CONTINUE + 811 I9LAST=I9 + K=K+1 + JOUT=JOUT+1 + RECORD(K)=KAPCEN + IF(K.LT.LENREC)GO TO 815 + IREC=IREC+1 + WRITE(14,REC=IREC)RECORD + K=0 + 815 CONTINUE + NSKIP=N9-I9LAST + IF(NSKIP.EQ.0)GO TO 825 + GO TO 821 + 820 NSKIP=N9 + 821 DO 822 ISKIP=1,NSKIP + K=K+1 + JOUT=JOUT+1 + RECORD(K)=0. + IF(K.LT.LENREC)GO TO 822 + IREC=IREC+1 + WRITE(14,REC=IREC)RECORD + K=0 + 822 CONTINUE + 825 IF(K.EQ.0)GO TO 830 + IREC=IREC+1 + WRITE(14,REC=IREC)RECORD + K=0 + 830 WRITE(6,817)J,MAXLINE,JOUT +C +C TRANSPOSE + IF(N9.LE.MAXLIN)GO TO 808 + WRITE(6,877) + 877 FORMAT(28H TOO MANY LINES TO TRANSPOSE) + CALL EXIT + 808 CONTINUE + REWIND 7 + NUMREC=(N9+LENREC-1)/LENREC + NLAST=N9-NUMREC*LENREC+LENREC + NCEN=0 + DO 899 N=1,NUMREC + WRITE(6,106)N + DO 893 J=1,NRHOX + IREC=J*NUMREC-NUMREC+N + READ(14,REC=IREC)(TRANSP(J,I),I=1,LENREC) + 893 CONTINUE + NOUT=LENREC + IF(N.EQ.NUMREC)NOUT=NLAST + DO 895 I=1,NOUT + NCEN=NCEN+1 + READ(7)LINDAT8,LINDAT4 + FREQ=2.99792458D17/WLVAC + DO 897 J=1,NRHOX + 897 ALINEC(J)=TRANSP(J,I)*(1.-EXP(-FREQ*HKT(J))) + 895 WRITE(9)LINDAT8,LINDAT4,ALINEC + 899 CONTINUE + CLOSE(UNIT=7,DISP='DELETE') + CLOSE(UNIT=15,DISP='DELETE') + CLOSE(UNIT=14,DISPOSE='DELETE') + WRITE(6,96)NCEN + CALL EXIT + END + SUBROUTINE XLINOP(J,N19,CUTOFF,VELSHIFT,IFVAC,LINOUT) + PARAMETER (kw=99) + PARAMETER (LENREC=8000,MAXLEN=2000001,MAXPROF=10000, + 1 MAXBUFF=MAXLEN+MAXPROF,MAXLIN=MAXBUFF+MAXPROF*2) + COMMON /BUFFER/BUFFER(MAXBUFF),PROFILE(MAXPROF) + COMMON /CONTIN/CONTINUUM(MAXBUFF) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) + COMMON /EXTAB/EXTAB(1001),EXTABF(1001),E1TAB(2000) + COMMON /NLINES/WLBEG,WLEND,RESOLU,RATIO,RATIOLG,WBEGIN, + 1 LENGTH,MLINES,IXWLBEG + COMMON /RHOX/RHOX(kw),NRHOX + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND, + 1 IFTURB + COMMON /TXNXN/EMERGE(kw),TXNXN(kw),BSTIM(kw),XNFH2(kw) + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + common /elem/abund(99),atmass(99),elem(99) + dimension hfactor(kw),hefactor(kw),h2factor(kw) + REAL NSTARK,NDOPP,NMERGE,INGLIS + EQUIVALENCE (GAMMAS,ASHORE),(GAMMAW,BSHORE) + EQUIVALENCE (GF,G,CGF),(TYPE,NLAST),(GAMMAR,XSECT,GAUNT) + INTEGER TYPE + REAL KAPPA,KAPMIN,KAPPA0,KAPCEN + REAL KAPPA0RED,KAPPARED,KAPPA0BLUE,KAPPABLUE + REAL*8 WAVE,WCON,EMERGE,WMERGE,WSHIFT,CONTX(26,17),WTAIL + REAL*8 RESOLU,RATIO,RATIOLG,SIGMA2,WLBEG,WLEND +c error 7jun2012 velshift is real*4 +C REAL*8 WL,DOPRATIO,VELSHIFT,WBEGIN,EMERGEH(kw) + REAL*8 WL,DOPRATIO,WBEGIN,EMERGEH(kw) + REAL*8 BLUECUT,WLPLUS1,WLPLUS2,REDCUT,WLMINUS1,WLMINUS2,VACAIR + REAL*8 EHYD(100),CONTH(15) ,ALPHAHYD(99) + DIMENSION DOPPH(kw) + DATA CONTH/ + 1 109678.764D0,27419.659D0,12186.462D0,6854.871,4387.113D0, 1.00 + 2 3046.604D0,2238.320D0,1713.711D0,1354.044D0,1096.776D0, + 3 906.426D0,761.650D0,648.980D0,559.579D0,487.456D0/ + DATA CONTX/ + 1 109678.764D0,27419.659D0,12186.462D0,6854.871D0,4387.113D0, 1.00 + 2 3046.604D0,2238.320D0,1713.711D0,1354.044D0,1096.776D0,16*0., + 3 198310.760D0,38454.691D0,32033.214D0,29223.753D0,27175.760D0, 2.00 + 4 15073.868D0,0.,0.,0.,0.,16*0., + 5 438908.850D0,109726.529D0,48766.491D0,27430.925D0,17555.715D0, 2.01 + 6 12191.437D0,0.,0.,0.,0.,16*0., + 7 90883.840D0,90867.420D0,90840.420D0,90820.420D0,90804.000D0, 6.00 + 8 90777.000D0,80691.180D0,80627.760D0,69235.820D0,69172.400D0, + 9 16*0., 26*0., 6.01 + A 61671.020D0,39820.615D0,39800.556D0,39759.842D0,22*0., 12.00 + 1 26*0., 12.01 + 2 48278.370D0,48166.309D0,0.,0.,0.,0.,0.,0.,0.,0.,16*0., 13.00 + 3 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,16*0., 13.01 + 4 66035.000D0,65957.885D0,65811.843D0,65747.550D0,65670.435D0, 14.00 + 5 65524.393D0,59736.150D0,59448.700D0,50640.630D0,50553.180D0, + 6 16*0., 26*0., 14.01 + 7 26*0.,26*0., 20.00,01 + 8 26*0.,26*0., 8.,11. + 9 26*0.,26*0./ 5.,19. + DATA ITEMP1/0/ +C + FASTEX(X)=EXTAB(IFIX(X)+1)* + 1EXTABF(IFIX((X-FLOAT(IFIX(X)))*1000.+1.5)) +C + IF(ITEMP.EQ.ITEMP1)GO TO 95 + EHYD(1)=0.D0 + EHYD(2)=82259.105D0 + EHYD(3)=97492.302D0 + EHYD(4)=102823.893D0 + EHYD(5)=105291.651D0 + EHYD(6)=106632.160D0 + EHYD(7)=107440.444D0 + EHYD(8)=107965.051D0 + DO 1 N=9,100 + 1 EHYD(N)=109678.764D0-109677.576D0/N**2 + DO 2 N=1,99 + 2 ALPHAHYD(N)=1.D7/(EHYD(N+1)-EHYD(N)) + WRITE(6,90) + 90 FORMAT(' NMERGE EMERGE EMERGEH') + DO 91 K=1,NRHOX +C FOR NEUTRALS FOR IONS NSTARK=NSTARK*Z**.25 NDOPP=NDOPP*Z**(2./3.) +C EMERGE=EMERGE*Z**2 +C INGLIS=1194./XNE(K)**.125 +C EMPIRICAL + INGLIS=1600./XNE(K)**(2./15.) +C NMERGE=INGLIS+.5 + NMERGE=INGLIS-1.5 + EMERGE(K)=109737.312D0/NMERGE**2 + EMERGEH(K)=109677.576D0/NMERGE**2 + 91 WRITE(6,92)K,NMERGE,EMERGE(K),EMERGEH(K) + 92 FORMAT(I3,4F10.3) + ITEMP1=ITEMP +C + 95 BOLT=1. + BOLTH=1. + OLDELO=1.E30*FLOAT(J*ITEMP) + OLDELOH=1.E30*FLOAT(J*ITEMP) + DOPRATIO=1.D0+VELSHIFT/299792.458D0 + alpha=0 + REWIND 19 + DO 900 ILINE=1,N19 + READ(19)WL,ELO,GF,NBLO,NBUP,NELION,TYPE,NCON,NELIONX, + 1GAMMAR,GAMMAS,GAMMAW,NBUFF +c +c include Barklem, Anstee, and O'Mara van der Waals +c READ(19)WL,ELO,GF,NBLO,NBUP,NELION,TYPE,NCON,NELIONX, +c 1GAMMAR,GAMMAS,GAMMAW,alpha,NBUFF +c + WL=WL*DOPRATIO + 97 IF(TYPE.EQ.2)GO TO 500 + IF(TYPE.EQ.0)GO TO 200 + IF(TYPE.EQ.-1)GO TO 600 + IF(TYPE.EQ.-2)GO TO 600 + IF(TYPE.EQ.1)GO TO 700 + IF(TYPE.EQ.3)GO TO 300 +C HE LINES DO NOT WORK YET + IF(TYPE.LT.-2)GO TO 200 +C +C MERGED CONTINUUM +C EDGE WAVELENGTHS ARE IN VACUUM + WSHIFT=1.D7/(1.D7/WL-109737.312D0/NLAST**2) + WMERGE=1.D7/(1.D7/WL-EMERGE(J)) + IF(NELION.EQ.1)THEN + WSHIFT=1.D7/(1.D7/WL-109677.576D0/NLAST**2) + WMERGE=1.D7/(1.D7/WL-EMERGEH(J)) + ENDIF + IF(WMERGE.LT.0.)WMERGE=WSHIFT+WSHIFT + WMERGE=MAX(WMERGE,WSHIFT) + WMERGE=MIN(WSHIFT+WSHIFT,WMERGE) + WTAIL=1.D7/(1.D7/WMERGE-500.) + IF(WTAIL.LT.0.)WTAIL=WMERGE+WMERGE + WTAIL=MIN(WMERGE+WMERGE,WTAIL) + IF(IFVAC.EQ.0)THEN + WMERGE=VACAIR(WMERGE)*DOPRATIO + WTAIL=VACAIR(WTAIL)*DOPRATIO + ENDIF + IXWL= LOG(WL)/RATIOLG + EDGEBLUE=EXP(IXWL*RATIOLG) + IF(EDGEBLUE.GT.WL)IXWL=IXWL-1 + NBUFF1=IXWL+1-IXWLBEG+1 + IXWL= LOG(WMERGE)/RATIOLG+.5 + NBUFF2=IXWL-IXWLBEG+1 + IXWL= LOG(WTAIL)/RATIOLG+.5 + NBUFF3=IXWL-IXWLBEG+1 + IF(NBUFF1.GT.LENGTH)GO TO 900 + IF(NBUFF3.LT.1)GO TO 900 + DNBUFF=NBUFF3-NBUFF2 + NBUFF1=MAX0(NBUFF1,1) + XSECTG=GF +C KAPPA=XSECTG*XNFPEL(NELION)*EXP(-ELO*HCKT(J)) + IF(NELION.GT.1)KAPPA=XSECTG*XNFPEL(NELION)*FASTEX(ELO*HCKT(J)) + IF(NELION.EQ.1)KAPPA=XSECTG*XNFPEL(NELION)* + 1FASTEX(ELO*HCKT(J)) + TAIL=1. + DO 190 IBUFF=NBUFF1,MIN0(NBUFF3,LENGTH) + IF(IBUFF.GT.NBUFF2)TAIL=(NBUFF3-IBUFF)/DNBUFF + 190 BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA*TAIL + GO TO 900 +C +C NORMAL LINE +C 200 KAPPA0=CGF*XNFPEL(NELION)/DOPPLE(NELION) + 200 KAPPA0=CGF*XNFDOP(NELION) + KAPMIN=CONTINUUM(MIN(MAX(NBUFF,1),LENGTH))*CUTOFF + IF(KAPPA0.LT.KAPMIN)GO TO 900 + IF(ELO.EQ.OLDELO)GO TO 210 +C BOLT=EXP(-ELO*HCKT(J)) + BOLT=FASTEX(ELO*HCKT(J)) + OLDELO=ELO + 210 KAPPA0=KAPPA0*BOLT + IF(KAPPA0.LT.KAPMIN)GO TO 900 + MLINES=MLINES+1 + WCON=0. + WTAIL=0. + IF(NCON.GT.0)THEN + WCON=1.E7/(CONTX(NCON,NELIONX)-EMERGE(J))*DOPRATIO + WTAIL=1.E7/(CONTX(NCON,NELIONX)-EMERGE(J)-500.)*DOPRATIO + ENDIF +c +c Castelli for Barklem, Anstee, and O'Mara van der Waals broadenint + if(alpha.ne.0.)then + nelem=int(nelion/6)+1 + v2=(1.-alpha)/2 + hfactor(j)=(t(j)/10000.)**v2 + hefactor(j)=0.628*(2.0991D-4*T(j)*(1/4+1.008/atmass(nelem)))**v2 + h2factor(j)=1.08*(2.0991D-4*T(j)*(1/2+1.008/atmass(nelem)))**v2 + txnxn(j)=xnfh(j)*hfactor(j)+xnfhe(J,1)*hefactor(j)+xnfh2(j)* + 1h2factor(j) + endif +c + ADAMP=(GAMMAR+GAMMAS*XNE(J)+GAMMAW*TXNXN(J))/DOPPLE(NELION) + KAPCEN=KAPPA0*VOIGT(0.,ADAMP) + IF(LINOUT.GE.0)WRITE(15)ILINE,KAPCEN + DOPWL=DOPPLE(NELION)*WL + IF(WL.GT.WLEND)GO TO 213 +C RED WING + MINRED=MAX0(1,NBUFF) + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 211 IBUFF=MINRED,LENGTH + IF(WAVE.LT.WCON)GO TO 211 + VVOIGT=ABS(WAVE-WL)/DOPWL + KAPPA=KAPPA0*VOIGT(VVOIGT,ADAMP) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) +C bug from john lester 23Jun03 +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 212 + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 212 + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + 211 WAVE=WAVE*RATIO + 212 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WBEGIN)GO TO 900 +C BLUE WING + 213 IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + WAVE=WBEGIN*RATIO**(IBUFF-1) + DO 214 I=1,MAXBLUE + IBUFF=IBUFF-1 + WAVE=WAVE/RATIO + IF(WAVE.LT.WCON)GO TO 214 + VVOIGT=ABS(WAVE-WL)/DOPWL + KAPPA=KAPPA0*VOIGT(VVOIGT,ADAMP) +C KAPPA=KAPPA0*VOIGT(ABS(WAVE-WL)/DOPWL,ADAMP) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 214 CONTINUE + GO TO 900 +C +C PRD LINE + 300 GO TO 200 +C +C CORONAL LINE + 500 GO TO 900 +C +C HYDROGEN LINE +C 600 KAPPA0=CGF*XNFPEL(1)/DOPPLE(1) + 600 KAPPA0=CGF*XNFDOP(1) + KAPMIN=CONTINUUM(MIN(MAX(NBUFF,1),LENGTH))*CUTOFF + IF(NBUP.EQ.2)KAPMIN=CONTINUUM(LENGTH)*CUTOFF + IF(KAPPA0.LT.KAPMIN)GO TO 900 +C IF(ELO.EQ.OLDELOH)GO TO 610 +C BOLTH=EXP(-ELO*HCKT(J)) + BOLTH=FASTEX(ELO*HCKT(J)) + OLDELOH=ELO + DOPPH(J)=DOPPLE(1) +C DEUTERIUM + IF(TYPE.EQ.-2)DOPPH(J)=DOPPH(J)/1.4142D0 + 610 KAPPA0=KAPPA0*BOLTH + IF(KAPPA0.LT.KAPMIN)GO TO 900 + IF(LINOUT.GE.0)WRITE(15)ILINE,KAPPA0 + MLINES=MLINES+1 + IF(NCON.EQ.0)GO TO 620 + IF(NBUP.EQ.NBLO+1)GO TO 620 + IF(NBUP.EQ.NBLO+2)GO TO 630 + WSHIFT=1.D7/(CONTH(NCON)-109677.576D0/81.D0**2) + WMERGE=1.D7/(CONTH(NCON)-EMERGEH(J)) + IF(WMERGE.LT.0.)WMERGE=WSHIFT+WSHIFT + WCON=MAX(WSHIFT,WMERGE) + WTAIL=1.D7/(1.D7/WCON-500.D0) + WCON=MIN(WSHIFT+WSHIFT,WCON) + IF(WTAIL.LT.0.)WTAIL=WCON+WCON + WTAIL=MIN(WCON+WCON,WTAIL) + IF(IFVAC.EQ.0)THEN + WCON=VACAIR(WCON) + WTAIL=VACAIR(WTAIL) + ENDIF + WCON=WCON*DOPRATIO + IF(IFVAC.EQ.0)WL=VACAIR(1.D7/(EHYD(NBUP)-EHYD(NBLO)))*DOPRATIO + IF(WL.GT.WLEND)GO TO 613 +C RED WING + IF(WBEGIN.GT.ALPHAHYD(NBLO))GO TO 613 + REDCUT=1.D7/(109678.764D0-109677.576D0/(NBUP-0.8D0)**2-EHYD(NBLO)) + IF(IFVAC.EQ.0)REDCUT=VACAIR(REDCUT) + REDCUT=REDCUT*DOPRATIO + WLMINUS1=1.D7/(EHYD(NBUP-1)-EHYD(NBLO)) + IF(IFVAC.EQ.0)WLMINUS1=VACAIR(WLMINUS1) + WLMINUS1=WLMINUS1*DOPRATIO + WLMINUS2=1.D7/(EHYD(NBUP-2)-EHYD(NBLO)) + IF(IFVAC.EQ.0)WLMINUS2=VACAIR(WLMINUS2) + WLMINUS2=WLMINUS2*DOPRATIO + KAPPA0RED=KAPPA0*HFNM(NBLO,NBUP-2)/HFNM(NBLO,NBUP)/ + 1 (EHYD(NBUP-2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) + MINRED=MAX0(1,NBUFF) + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 611 IBUFF=MINRED,LENGTH + IF(WAVE.LT.WCON)GO TO 611 + IF(WAVE.GT.WLMINUS1)GO TO 612 + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + IF(WAVE.GT.REDCUT)THEN + KAPPARED=KAPPA0RED*HPROF4(NBLO,NBUP-2,J,WAVE-WLMINUS2,DOPPH) + IF(WAVE.LT.WTAIL)KAPPARED=KAPPARED*(WAVE-WCON)/(WTAIL-WCON) + IF(KAPPARED.GE.KAPPA)GO TO 612 + ENDIF + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 612 + 611 WAVE=WAVE*RATIO + 612 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WBEGIN)GO TO 900 +C BLUE WING + 613 BLUECUT=1.D7/(109678.764D0-109677.576D0/(NBUP+0.8D0)**2- + 1 EHYD(NBLO)) + IF(IFVAC.EQ.0)BLUECUT=VACAIR(BLUECUT) + BLUECUT=BLUECUT*DOPRATIO + WLPLUS1=1.D7/(EHYD(NBUP+1)-EHYD(NBLO)) + IF(IFVAC.EQ.0)WLPLUS1=VACAIR(WLPLUS1) + WLPLUS1=WLPLUS1*DOPRATIO + WLPLUS2=1.D7/(EHYD(NBUP+2)-EHYD(NBLO)) + IF(IFVAC.EQ.0)WLPLUS2=VACAIR(WLPLUS2) + WLPLUS2=WLPLUS2*DOPRATIO + KAPPA0BLUE=KAPPA0*HFNM(NBLO,NBUP+2)/HFNM(NBLO,NBUP)/ + 1 (EHYD(NBUP+2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) + IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + WAVE=WBEGIN*RATIO**(IBUFF-1) + DO 614 I=1,MAXBLUE + IBUFF=IBUFF-1 + WAVE=WAVE/RATIO + IF(WAVE.LT.WCON)GO TO 900 + IF(WAVE.LT.WLPLUS1)GO TO 900 + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + IF(WAVE.LT.BLUECUT)THEN + KAPPABLUE=KAPPA0BLUE*HPROF4(NBLO,NBUP+2,J,WAVE-WLPLUS2,DOPPH) + IF(WAVE.LT.WTAIL)KAPPABLUE=KAPPABLUE*(WAVE-WCON)/(WTAIL-WCON) + IF(KAPPABLUE.GT.KAPPA)GO TO 900 + ENDIF + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 614 CONTINUE + GO TO 900 +C ALPHA PRESUMED ISOLATED LINES AND BETA BLUE WINGS + 620 IF(WL.GT.WLEND)GO TO 623 +C RED WING + MINRED=MAX0(1,NBUFF) + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 621 IBUFF=MINRED,LENGTH + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 622 + 621 WAVE=WAVE*RATIO + 622 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WBEGIN)GO TO 900 +C ALPHA OR BETA BLUE WING + 623 IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + WAVE=WBEGIN*RATIO**(IBUFF-1) + DO 624 I=1,MAXBLUE + IBUFF=IBUFF-1 + WAVE=WAVE/RATIO + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 624 CONTINUE + GO TO 900 +C BETA LINES RED WING + 630 IF(WL.GT.WLEND)GO TO 623 +C RED WING + MINRED=MAX0(1,NBUFF) + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 631 IBUFF=MINRED,LENGTH + IF(WAVE.GT.ALPHAHYD(NBLO))GO TO 633 + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 622 + 631 WAVE=WAVE*RATIO + GO TO 622 + 633 IBUFFALPHA=IBUFF + IF(NBUFF.LT.1)GO TO 900 + IBUFF=NBUFF + MAXBLUE=IBUFF-1 + WAVE=WBEGIN*RATIO**(IBUFF-1) + DO 634 I=1,MAXBLUE + IBUFF=IBUFF-1 + WAVE=WAVE/RATIO + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 634 CONTINUE + GO TO 900 +C +C AUTOIONIZING LINE + 700 KAPPA0=BSHORE*G*XNFPEL(NELION) + KAPMIN=CONTINUUM(MIN(MAX(NBUFF,1),LENGTH))*CUTOFF + IF(KAPPA0.LT.KAPMIN)GO TO 900 + KAPPA0=KAPPA0*FASTEX(ELO*HCKT(J)) +C KAPPA0=KAPPA0*EXP(-ELO*HCKT(J)) + IF(KAPPA0.LT.KAPMIN)GO TO 900 + IF(LINOUT.GE.0)WRITE(15)ILINE,KAPPA0 + MLINES=MLINES+1 + FRELIN=2.99792458E17/WL + IF(WL.GT.WLEND)GO TO 713 +C RED WING + MINRED=MAX0(1,NBUFF) + FREQ=2.99792458D17/(WBEGIN*RATIO**(MINRED-1)) + DO 711 IBUFF=MINRED,LENGTH + EPSIL=2.*(FREQ-FRELIN)/GAMMAR + KAPPA=KAPPA0*(ASHORE*EPSIL+BSHORE)/(EPSIL**2+1.)/BSHORE + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 712 + 711 FREQ=FREQ/RATIO + 712 IF(NBUFF.EQ.1)GO TO 900 + IF(WL.LT.WBEGIN)GO TO 900 +C BLUE WING + 713 IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + FREQ=2.99792458D17/(WBEGIN*RATIO**(IBUFF-1)) + DO 714 I=1,MAXBLUE + IBUFF=IBUFF-1 + FREQ=FREQ*RATIO + EPSIL=2.*(FREQ-FRELIN)/GAMMAR + KAPPA=KAPPA0*(ASHORE*EPSIL+BSHORE)/(EPSIL**2+1.)/BSHORE + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 714 CONTINUE + GO TO 900 +c +c He I +c +C LINE=1 4471 +C LINE=2 4026 +C LINE=3 4387 +C LINE=4 4921 +C LINE>4 OTHERS + 800 KAPPA0=CGF*XNFPEL(NELION)/DOPPLE(NELION) +C 3HE + IF(TYPE.EQ.-4)KAPPA0=KAPPA0/1.155 + KAPMIN=CONTINUUM(MIN(MAX(NBUFF,1),LENGTH))*CUTOFF + IF(KAPPA0.LT.KAPMIN)GO TO 900 + BOLT=FASTEX(ELO*HCKT(J)) + KAPPA0=KAPPA0*BOLT + IF(KAPPA0.LT.KAPMIN)GO TO 900 +c DOPWL=DOPPLE(7)*WL + DOPWL=DOPPLE(7)*WL4 + IF(TYPE.EQ.-4)DOPWL=DOPWL*1.155 +c KAPCEN=KAPPA0*HE1PROF(J,WL,WL,DOPWL,GAMMAR,GAMMAS) + KAPCEN=KAPPA0*HE1PROF(J,WL4,WL4,DOPWL,GAMMAR,GAMMAS) + WRITE(15)ILINE,KAPCEN + MLINES=MLINES+1 + IF(WL.GT.WLEND)GO TO 813 +C RED WING + REDCUT=WLEND + MINRED=MAX0(1,NBUFF) + WAVE=WBEGIN*RATIO**(MINRED-1) + wave4=wave + DO 811 IBUFF=MINRED,LENGTH + IF(WAVE.GT.REDCUT)GO TO 812 + KAPPA=KAPPA0*HE1PROF(J,WAVE4,WL4,DOPWL,GAMMAR,GAMMAS) + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 812 + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + WAVE=WAVE*RATIO + 811 wave4=wave + 812 IF(NBUFF.EQ.1)GO TO 900 + IF(WL.LT.WLBEG)GO TO 900 +C BLUE WING + 813 BLUECUT=WLBEG + IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + WAVE=WBEGIN*RATIO**(IBUFF-1) + DO 814 I=1,MAXBLUE + IBUFF=IBUFF-1 + WAVE=WAVE/RATIO + wave4=wave + IF(WAVE.LT.BLUECUT)GO TO 900 + KAPPA=KAPPA0*HE1PROF(J,WAVE4,WL4,DOPWL,GAMMAR,GAMMAS) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF)*CUTOFF)GO TO 900 + 814 CONTINUE + GO TO 900 + 850 CONTINUE + 900 CONTINUE + RETURN + END + FUNCTION HPROF4(N,M,J,DELW,DOPPH) +C FUNCTION HPROFL(N,M,J,DELW) +C VERSION FINE STRUCTURE LIKE GENERAL BUT APPROXIMATELY INCLUDES FINE +C STRUCTURE IN THE DOPPLER CORES. EXACT PATTERN +C IS USED FOR ALPHA LINES, M INFINITE PATTERN +C IS USED FOR ALL OTHER LINES. +C FROM DEANE PETERSON +C REQUIRES VCSE1F AND SOFBET + PARAMETER (kw=99) + REAL*8 DELW,EMERGE + REAL*4 DOPPH(kw) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) + COMMON /RHOX/RHOX(kw),NRHOX + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /TXNXN/EMERGE(kw),TXNXN(kw),BSTIM(kw),XNFH2(kw) + COMMON /EXTAB/EXTAB(1001),EXTABF(1001),E1TAB(2000) + DIMENSION PP(kw),FO(kw),GCON1(kw),GCON2(kw),Y1B(kw),Y1S(kw), + 1C1D(kw),C2D(kw),Y1WTM(2,2),XKNMTB(4,3) + DIMENSION T3NHE(kw),T3NH2(kw) + DIMENSION STCOMP(5,4),STALPH(34),ISTAL(4),LNGHAL(4),STWTAL(34), + 1STCPWT(5,4),LNCOMP(4),FINEST(14),FINSWT(14) + REAL*4 LORWING,ASUM(100),ASUMLYMAN(100) + REAL*4 CUTOFFH2PLUS(111),CUTOFFH2(91) +C + DATA XKNMTB/.0001716,.009019,.1001,.5820,.0005235,.01772,.171,.866 + 1,.0008912,.02507,.223,1.02/ + DATA Y1WTM/1.E18,1.E17,1.E16,1.E14/ + DATA ITEMP1/0/,N1/0/,M1/0/,RYDH/3.2880515E15/ +C FINE STRUCTURE COMPONENTS FOR ALPHA LINES IN FREQ*10**-7 + DATA STALPH/-730.,370.,188.,515.,327.,619.,-772.,-473.,-369.,120., + 1256.,162.,285.,-161.,-38.3,6.82,-174.,-147.,-101.,-77.5,55.,126., + 275.,139.,-60.,3.7,27.,-69.,-42.,-18.,-5.5,-9.1,-33.,-24./ +C ALPHA COMPONENT WEIGHTS + DATA STWTAL/1.,2.,1.,2.,1.,2.,1.,2.,3.,1.,2.,1.,2.,1.,4.,6.,1.,2., + 13.,4.,1.,2.,1.,2.,1.,4.,6.,1.,7.,6.,4.,4.,4.,5./ + DATA ISTAL/1,3,10,21/, LNGHAL/2,7,11,14/ +C FINE STRUCTURE FOR M.EQ.INFINITY IN FREQ*10**-7 + DATA STCOMP/0.,0.,0.,0.,0.,468.,576.,-522.,0.,0.,260.,290.,-33., + 1-140.,0.0,140.,150.,18.,-27.,-51./ +C WEIGHTS + DATA STCPWT/1.,0.,0.,0.,0.,1.,1.,2.,0.,0.,1.,1.,4.,3.,0.,1.,1., + 14.,6.,4./ + DATA LNCOMP/1,3,4,5/ +c revised by Fiorella Castelli +c Allard et al. (1998) A&A 372,260. Castelli & Kurucz (2001)a&A 372, 360. +C LYMAN ALPHA QUASI H2+ CUTOFF +C DELTA WAVENO = -15000+100*(N-1) N=1,111 UP TO -4000 + DATA CUTOFFH2PLUS/ + 1-15.14,-15.06,-14.97,-14.88,-14.80,-14.71,-14.62,-14.53, + 2-14.44,-14.36,-14.27,-14.18,-14.09,-14.01,-13.92,-13.83, + 3-13.74,-13.65,-13.57,-13.48,-13.39,-13.30,-13.21,-13.13, + 4-13.04,-12.95,-12.86,-12.77,-12.69,-12.60,-12.51,-12.40, + 5-12.29,-12.15,-12.02,-11.90,-11.76,-11.63,-11.53,-11.41, + 6-11.30,-11.22,-11.15,-11.09,-11.07,-11.06,-11.07,-11.09, + 7-11.12,-11.16,-11.19,-11.21,-11.24,-11.27,-11.30,-11.33, + 8-11.36,-11.39,-11.42,-11.45,-11.48,-11.48,-11.48,-11.48, + 9-11.48,-11.48,-11.48,-11.48,-11.48,-11.48,-11.48,-11.48, + A-11.48,-11.48,-11.48,-11.48,-11.41,-11.40,-11.39,-11.38, + 1-11.37,-11.36,-11.35,-11.34,-11.33,-11.32,-11.30,-11.29, + 2-11.28,-11.27,-11.27,-11.27,-11.26,-11.25,-11.24,-11.23, + 3-11.22,-11.21,-11.20,-11.19,-11.18,-11.17,-11.15,-11.14, + 4-11.13,-11.12,-11.11,-11.10,-11.09,-11.08,-11.07/ +C DATA CUTOFFH2PLUS/ +C 1-17.15,-17.00,-16.86,-16.71,-16.56,-16.41,-16.27,-16.12, +C 2-15.97,-15.82,-15.68,-15.53,-15.38,-15.23,-15.09,-14.94, +C 3-14.79,-14.64,-14.50,-14.35,-14.20,-14.05,-13.91,-13.76, +C 4-13.61,-13.46,-13.32,-13.17,-13.02,-12.87,-12.73,-12.58, +C 5-12.43,-12.28,-12.14,-11.99,-11.85,-11.72,-11.60,-11.49, +C 6-11.40,-11.36,-11.335,-11.33,-11.33,-11.335,-11.35,-11.37, +C 7-11.40,-11.42,-11.45,-11.48,-11.52,-11.55,-11.59,-11.61, +C 8-11.61,-11.61,-11.60,-11.59,-11.58,-11.57,-11.56,-11.55, +C 9-11.54,-11.53,-11.52,-11.51,-11.49,-11.48,-11.47,-11.46, +C A-11.45,-11.44,-11.43,-11.42,-11.41,-11.40,-11.39,-11.38, +C 1-11.37,-11.36,-11.35,-11.34,-11.33,-11.32,-11.30,-11.29, +C 2-11.28,-11.27,-11.26,-11.25,-11.24,-11.23,-11.22,-11.21, +C 3-11.20,-11.19,-11.18,-11.17,-11.16,-11.15,-11.13,-11.12, +C 4-11.11,-11.10,-11.09,-11.08,-11.07,-11.06,-11.05/ +C LYMAN ALPHA QUASI H2 CUTOFF +C DELTA WAVENO = -22000+200*(N-1) N=1,91 -4000 + DATA CUTOFFH2/ + 1-13.64,-13.52,-13.39,-13.27,-13.14,-13.01,-12.87,-12.74, + 2-12.63,-12.56,-12.51,-12.48,-12.47,-12.49,-12.52,-12.55, + 3-12.57,-12.61,-12.65,-12.69,-12.72,-12.76,-12.79,-12.82, + 4-12.84,-12.85,-12.87,-12.90,-12.93,-12.94,-12.93,-12.95, + 5-12.95,-12.96,-12.97,-12.96,-12.96,-12.95,-12.95,-12.96, + 6-12.98,-12.99,-12.95,-12.96,-13.00,-13.00,-12.98,-12.97, + 7-13.00,-13.00,-13.00,-13.00,-13.00,-13.00,-13.00,-13.00, + 8-13.00,-13.00,-13.00,-13.00,-13.00,-13.00,-13.00,-13.00, + 9-13.00,-13.00,-13.00,-13.00,-13.00,-13.00,-12.89,-12.88, + A-12.87,-12.86,-12.85,-12.84,-12.83,-12.81,-12.80,-12.79, + 1-12.78,-12.76,-12.74,-12.72,-12.70,-12.68,-12.65,-12.62, + 2-12.59,-12.56,-12.53/ +cc DATA CUTOFFH2/ +cc 1-13.43,-13.32,-13.21,-13.10,-12.98,-12.86,-12.79,-12.72, +cc 2-12.65,-12.58,-12.51,-12.47,-12.45,-12.45,-12.48,-12.51, +cc 3-12.53,-12.56,-12.59,-12.62,-12.65,-12.69,-12.73,-12.77, +cc 4-12.81,-12.85,-12.87,-12.89,-12.90,-12.90,-12.90,-12.90, +cc 5-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90, +cc 6-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90, +cc 7-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90, +cc 8-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.90, +cc 9-12.90,-12.90,-12.90,-12.90,-12.90,-12.90,-12.89,-12.88, +cc A-12.87,-12.86,-12.85,-12.84,-12.83,-12.81,-12.80,-12.79, +cc 1-12.78,-12.76,-12.74,-12.72,-12.70,-12.68,-12.65,-12.62, +cc 2-12.59,-12.56,-12.53/ +cc +C DATA CUTOFFH2/ +C 1-13.87,-13.73,-13.59,-13.45,-13.31,-13.17,-13.05,-12.93, +C 2-12.83,-12.74,-12.67,-12.63,-12.60,-12.60,-12.63,-12.65, +C 3-12.67,-12.69,-12.71,-12.73,-12.75,-12.77,-12.79,-12.81, +C 4-12.83,-12.85,-12.86,-12.87,-12.87,-12.86,-12.85,-12.84, +C 5-12.83,-12.82,-12.80,-12.79,-12.78,-12.77,-12.76,-12.75, +C 6-12.74,-12.72,-12.71,-12.70,-12.69,-12.68,-12.67,-12.66, +C 7-12.64,-12.63,-12.62,-12.61,-12.60,-12.59,-12.58,-12.57, +C 8-12.55,-12.54,-12.53,-12.52,-12.51,-12.50,-12.49,-12.47, +C 9-12.46,-12.45,-12.44,-12.43,-12.42,-12.41,-12.39,-12.38, +C A-12.37,-12.36,-12.35,-12.34,-12.33,-12.31,-12.30,-12.29, +C 1-12.28,-12.26,-12.24,-12.22,-12.20,-12.18,-12.15,-12.12, +C 2-12.09,-12.06,-12.03/ + DATA ASUMLYMAN/ + 1 0.000E+00, 6.265E+08, 1.897E+08, 8.126E+07, 4.203E+07, 2.450E+07, + 2 1.236E+07, 8.249E+06, 5.782E+06, 4.208E+06, 3.158E+06, 2.430E+06, + 3 1.910E+06, 1.567E+06, 1.274E+06, 1.050E+06, 8.752E+05, 7.373E+05, + 4 6.269E+05, 5.375E+05, 4.643E+05, 4.038E+05, 3.534E+05, 3.111E+05, + 5 2.752E+05, 2.447E+05, 2.185E+05, 1.959E+05, 1.763E+05, 1.593E+05, + 6 1.443E+05, 1.312E+05, 1.197E+05, 1.094E+05, 1.003E+05, 9.216E+04, + 7 8.489E+04, 7.836E+04, 7.249E+04, 6.719E+04, 6.239E+04, 5.804E+04, + 8 5.408E+04, 5.048E+04, 4.719E+04, 4.418E+04, 4.142E+04, 3.888E+04, + 9 3.655E+04, 3.440E+04, 3.242E+04, 3.058E+04, 2.888E+04, 2.731E+04, + A 2.585E+04, 2.449E+04, 2.322E+04, 2.204E+04, 2.094E+04, 1.991E+04, + 1 1.894E+04, 1.804E+04, 1.720E+04, 1.640E+04, 1.566E+04, 1.496E+04, + 2 1.430E+04, 1.368E+04, 1.309E+04, 1.254E+04, 1.201E+04, 1.152E+04, + 3 1.105E+04, 1.061E+04, 1.019E+04, 9.796E+03, 9.419E+03, 9.061E+03, + 4 8.721E+03, 8.398E+03, 8.091E+03, 7.799E+03, 7.520E+03, 7.255E+03, + 5 7.002E+03, 6.760E+03, 6.530E+03, 6.310E+03, 6.100E+03, 5.898E+03, + 6 5.706E+03, 5.522E+03, 5.346E+03, 5.177E+03, 5.015E+03, 4.860E+03, + 7 4.711E+03, 4.569E+03, 4.432E+03, 4.300E+03/ + DATA ASUM/ + 1 0.000E+00, 4.696E+08, 9.980E+07, 3.017E+07, 1.155E+07, 5.189E+06, + 2 2.616E+06, 1.437E+06, 8.444E+05, 5.234E+05, 3.389E+05, 2.275E+05, + 3 1.575E+05, 1.120E+05, 8.142E+04, 6.040E+04, 4.560E+04, 3.496E+04, + 4 2.719E+04, 2.141E+04, 1.711E+04, 1.377E+04, 1.119E+04, 9.166E+03, + 5 7.572E+03, 6.341E+03, 5.338E+03, 4.523E+03, 3.854E+03, 3.302E+03, + 6 2.844E+03, 2.460E+03, 2.138E+03, 1.866E+03, 1.635E+03, 1.438E+03, + 7 1.269E+03, 1.124E+03, 9.983E+02, 8.894E+02, 7.947E+02, 7.120E+02, + 8 6.396E+02, 5.759E+02, 5.198E+02, 4.703E+02, 4.263E+02, 3.873E+02, + 9 3.526E+02, 3.215E+02, 2.938E+02, 2.689E+02, 2.465E+02, 2.264E+02, + A 2.082E+02, 1.918E+02, 1.769E+02, 1.634E+02, 1.512E+02, 1.400E+02, + 1 1.298E+02, 1.206E+02, 1.121E+02, 1.043E+02, 9.720E+01, 9.066E+01, + 2 8.465E+01, 7.912E+01, 7.403E+01, 6.933E+01, 6.498E+01, 6.097E+01, + 3 5.725E+01, 5.381E+01, 5.061E+01, 4.765E+01, 4.489E+01, 4.232E+01, + 4 3.994E+01, 3.771E+01, 3.563E+01, 3.369E+01, 3.188E+01, 3.019E+01, + 5 2.860E+01, 2.712E+01, 2.572E+01, 2.442E+01, 2.319E+01, 2.204E+01, + 6 2.096E+01, 1.994E+01, 1.898E+01, 1.808E+01, 1.722E+01, 1.642E+01, + 7 1.566E+01, 1.495E+01, 1.427E+01, 1.363E+01/ + FASTEX(X)=EXTAB(IFIX(X)+1)* + 1EXTABF(IFIX((X-FLOAT(IFIX(X)))*1000.+1.5)) + IF(ITEMP.EQ.ITEMP1)GO TO 20 +C SET UP DEPTH VECTORS + ITEMP1=ITEMP + DO 10 K=1,NRHOX + XNE16=XNE(K)**.1666667 + PP(K)=XNE16*.08989/SQRT(T(K)) + FO(K)=XNE16**4*1.25E-9 + Y1B(K)=2./(1.+.012/T(K)*SQRT(XNE(K)/T(K))) + T4=T(K)/10000. + T43=T4**.3 + Y1S(K)=T43/XNE16 +C T3NHE(K)=T43*XNFPHE(K,1) + T3NHE(K)=T43*XNFHE(K,1) + T3NH2(K)=T43*XNFH2(K) + C1D(K)=FO(K)*78940./T(K) + C2D(K)=FO(K)**2/5.96E-23/XNE(K) + GCON1(K)=.2+.09*SQRT(T4)/(1.+XNE(K)/1.E13) + GCON2(K)=.2/(1.+XNE(K)/1.E15) + 10 CONTINUE +C SET UP FOR THIS LINE + 20 IF(N.EQ.N1.AND.M.EQ.M1)GO TO 30 + N1=N + M1=M + MMN=M-N + XN=N + XN2=XN*XN + XM=M + XM2=XM*XM + XMN2=XM2*XN2 + XM2MN2=XM2-XN2 + GNM=XM2MN2/XMN2 + IF(MMN.LE.3.AND.N.LE.4)XKNM=XKNMTB(N,MMN) + IF(MMN.GT.3.OR.N.GT.4)XKNM=5.5E-5/GNM*XMN2/(1.+.13/FLOAT(MMN)) + Y1NUM=320. + IF(M.EQ.2)Y1NUM=550. + IF(M.EQ.3)Y1NUM=380. + Y1WHT=1.E13 + IF(MMN.LE.3)Y1WHT=1.E14 + IF(MMN.LE.2.AND.N.LE.2)Y1WHT=Y1WTM(N,MMN) + FREQNM=RYDH*GNM + DBETA=2.99792458D18/FREQNM**2/XKNM + WAVENM=2.99792458D18/FREQNM + C1CON=XKNM/WAVENM*GNM*XM2MN2 + C2CON=(XKNM/WAVENM)**2 +C RADAMP=1.389E9/XM**4.53/(1.+5./XM2/XM) +C IF(N.NE.1)RADAMP=RADAMP+1.389E9/XN**4.53/(1.+5./XN2/XN) + RADAMP=ASUM(N)+ASUM(M) + IF(N.EQ.1)RADAMP=ASUMLYMAN(M) + RADAMP=RADAMP/12.5664 + RADAMP=RADAMP/FREQNM + RESONT=HFNM(1,M)/XM/(1.-1./XM2) + IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2) +C FUDGE TO BASCHEK*2 +C RESONT=HFNM(1,M)/XM/(1.-1./XM2)*XM/3.*.791*2. +C IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2)*XN/3.*.791*2. +C 2 IS FOR CONVERTING XNFPH TO XNFH +C RESONT=RESONT*5.593E-24/GNM*2. +C RESONT=RESONT*5.593E-24/GNM +C error in constant corrected 26nov95 +C KURUCZ +CC RESONT=RESONT*3.92E-24/GNM +C Castelli after Ali and Griem (1963) + RESONT=RESONT*3.579E-24/GNM + VDW=4.45E-26/GNM*(XM2*(7.*XM2+5.))**.4 +C GUESS THAT H2 IS TWICE AS STRONG AS HE AS IN TXNXN + HWVDW=VDW*T3NHE(J)+2.*VDW*T3NH2(J) + HWRAD=RADAMP + STARK=1.6678E-18*FREQNM*XKNM +C FINE STRUCTURE COMPONENTS +C +C IF(N.GT.4)THEN + IF(N.GT.4.OR.M.GT.10)THEN + IFINS=1 + FINEST(1)=0. + FINSWT(1)=1. + GO TO 30 + ENDIF +C + IF(MMN.EQ.1)GO TO 22 +C USE M.EQ.INF STRUCTURE + IFINS=LNCOMP(N) + DO 21 I=1,IFINS + FINEST(I)=STCOMP(I,N)*1.E7 + FINSWT(I)=STCPWT(I,N)/XN2 + 21 CONTINUE + GO TO 30 +C FOR ALPHA LINES + 22 IFINS=LNGHAL(N) + IPOS=ISTAL(N) + DO 23 I=1,IFINS + K=IPOS-1+I + FINEST(I)=STALPH(K)*1.E7 + FINSWT(I)=STWTAL(K)/XN2/3. + 23 CONTINUE +C NOW DO THIS DEPTH +C 30 DEL=-10.*DELW/WAVENM*FREQNM error +C FREQ=FREQNM+DEL +C WAVENM IN A DELW IN NM + 30 WL=WAVENM+DELW*10. + FREQ=2.99792458D18/WL + DEL=ABS(FREQ-FREQNM) +C WL IN NM + WL=WL/10. +C THESE HALF-WIDTHS ARE REALLY DNU/NU + HWSTK=STARK*FO(J) + HWVDW=VDW*T3NHE(J)+2.*VDW*T3NH2(J) + HWRAD=RADAMP +C XNFPH(J,1)*2 IS THE NUMBER IN THE GROUND STATE + HWRES=RESONT*XNFPH(J,1)*2. + HWLOR=HWRES+HWVDW+HWRAD + HWDOP=DOPPH(J) +C SPECIFY LARGEST HALF WIDTH IN CASE OF CORE CALC +C NWID=1, DOPPLER =2, LORENTZ =3, STARK + NWID=1 + IF(HWDOP.GE.HWSTK.AND.HWDOP.GE.HWLOR)GO TO 31 + NWID=2 + IF(HWLOR.GE.HWSTK)GO TO 31 + NWID=3 + 31 HFWID=FREQNM*AMAX1(HWDOP,HWLOR,HWSTK) +C SETS FLAG IF IN A LINE CORE +C HPROFL=0. + HPROF4=0. + IFCORE=0 + IF(ABS(DEL).LE.HFWID)IFCORE=1 + DOP=FREQNM*HWDOP + IF(IFCORE.EQ.1)GO TO (32,40,50),NWID +C +C DO DOPPLER +C PUT FINE STRUCTURE IN DOPPLER CORE + 32 DO 33 I=1,IFINS + D=ABS(FREQ-FREQNM-FINEST(I))/DOP +C IF(D.LE.7.)HPROF4=HPROF4+EXP(-D*D)/1.77245/DOP*FINSWT(I) +C IF(D.LE.7.)HPROFL=HPROFL+EXP(-D*D)/1.77245/DOP*FINSWT(I) +C SAME NORMALIZATION AS VOIGT FUNCTION + IF(D.LE.7.)HPROF4=HPROF4+FASTEX(D*D)*FINSWT(I) + 33 CONTINUE + IF(IFCORE.EQ.1)RETURN +C +C DO LORENTZ + 40 IF(N.NE.1)GO TO 48 + IF(M.NE.2)GO TO 48 +C Lyman alpha +C near center +c modify old resonance broadening to match at 4000 cm-1 + HWRES=HWRES*4. + HWLOR=HWRES+HWVDW+HWRAD + HHW=FREQNM*HWLOR + IF(FREQ.GT.(82259.10D0-4000.D0)*2.99792458D10)THEN + HPROFRES=HWRES*FREQNM/3.14159/(DEL**2+HHW**2)*1.77245*DOP + GO TO 44 + ENDIF +C only far red wing +C Data from N.F. Allard, March 1997 +C Insert Lyman alpha cutoff a la N.F. Allard and D. Koester, A&A, 258, +C 464-468. 1992. + CUTOFF=0. + IF(FREQ.LT.50000.*2.99792458D10)GO TO 43 +C TABULATED AT 200 CM-1 SPACING + SPACING=200.*2.99792458D10 + FREQ22000=(82259.10D0-22000.D0)*2.99792458D10 + IF(FREQ.LT.FREQ22000)THEN + CUTOFF=(CUTOFFH2(2)-CUTOFFH2(1))/SPACING*(FREQ-FREQ22000)+ + 1CUTOFFH2(1) + ELSE + ICUT=(FREQ-FREQ22000)/SPACING + CUTFREQ=ICUT*SPACING+FREQ22000 + CUTOFF=(CUTOFFH2(ICUT+2)-CUTOFFH2(ICUT+1))/SPACING* + 1(FREQ-CUTFREQ)+CUTOFFH2(ICUT+1) + ENDIF + CUTOFF=(10.**(CUTOFF-14.))*XNFPH(J,1)*2./2.99792458D10 + 43 HPROFRES=CUTOFF*1.77245*DOP +C + 44 HPROFRAD=0. +C RAYLEIGH SCATTERING EXCEPT NEAR DOPPLER CORE +C CROSSOVER FROM ABSORPTION TO RAYLEIGH SCATTTERING IN HRAYOP +C IF(FREQ.GT.2.463E15) +C IF(FREQ.GT..74*3.288051D15.AND.FREQ.LT..78*3.288051D15) +C USE FREQUENCY IN CONTINUA.DAT 2.419061115 AS CUTOFF INSTEAD + IF(FREQ.GT.2.4190611E15.AND.FREQ.LT..77*3.288051D15) + 1 HPROFRAD=HWRAD*FREQNM/3.14159/(DEL**2+HHW**2)*1.77245*DOP +C CORRECTION TO LORENTZ PROFILE ALLER P.164 NOT USED +C HPROFRAD=HPROFRAD*4.*FREQ**2/(FREQ**2+FREQNM**2) +C + HPROFVDW=HWVDW*FREQNM/3.14159/(DEL**2+HHW**2)*1.77245*DOP +C VAN DER WAALS CUTOFF FOR HE AND FOR H2 +C GUESS BOTH 60000 CM-1 = 1.8D15 HZ + IF(FREQ.LT.1.8D15)HPROFVDW=0. + HPROFLOR=HPROFRES+HPROFRAD+HPROFVDW + HPROF4=HPROF4+HPROFLOR + IF(IFCORE.EQ.1)RETURN + GO TO 50 +C +C not Lyman alpha + 48 HHW=FREQNM*HWLOR + TOP=HHW + IF(N.NE.1.OR.M.GT.5)GO TO 49 +C Lyman beta + IF(M.EQ.3.AND.FREQ.GT..885*3.288051D15.AND. + 1 FREQ.LT..890*3.288051D15)TOP=HHW-FREQNM*HWRAD +C Lyman gamma + IF(M.EQ.4.AND.FREQ.GT..936*3.288051D15.AND. + 1 FREQ.LT..938*3.288051D15)TOP=HHW-FREQNM*HWRAD +C Lyman delta + IF(M.EQ.5.AND.FREQ.GT..959*3.288051D15.AND. + 1 FREQ.LT..961*3.288051D15)TOP=HHW-FREQNM*HWRAD + 49 HPROFLOR=TOP/3.14159/(DEL**2+HHW**2)*1.77245*DOP + HPROF4=HPROF4+HPROFLOR + IF(IFCORE.EQ.1)RETURN +C +C DO STARK + 50 WTY1=1./(1.+XNE(J)/Y1WHT) + Y1SCAL=Y1NUM*Y1S(J)*WTY1+Y1B(J)*(1.-WTY1) + C1=C1D(J)*C1CON*Y1SCAL + C2=C2D(J)*C2CON + G1=6.77*SQRT(C1) + GNOT=G1*AMAX1(0.,.2114+ALOG(SQRT(C2)/C1))*(1.-GCON1(J)-GCON2(J)) + BETA=ABS(DEL)/FO(J)*DBETA + Y1=C1*BETA + Y2=C2*BETA**2 + GAM=GNOT +C IF(Y2.LE..001)GO TO 51 + IF(Y2.LE.1.E-4.AND.Y1.LE.1.E-5)GO TO 51 +C GAM=G1*(.5*FASTEX(AMIN1(80.,Y1))+FASTE1(Y1)-.5*FASTE1(Y2))* +C GAM=G1*(.5*EXP(-AMIN1(80.,Y1))+VCSE1F(Y1)-.5*VCSE1F(Y2))* + GAM=G1*(.5*FASTEX(AMIN1(80.,Y1))+VCSE1F(Y1)-.5*VCSE1F(Y2))* + 1(1.-GCON1(J)/(1.+(90.*Y1)**3)-GCON2(J)/(1.+2000.*Y1)) + IF(GAM.LE.1.E-20)GAM=0. + 51 PRQS=SOFBET(BETA,PP(J),N,M) + IF(M.GT.2)GO TO 53 +C ASSUME QUASISTATIC PROFILE IS HALF PROTONS, HALF ELECTRONS + PRQS=PRQS*.5 + CUTOFF=0. +C LYMAN ALPHA QUASI H2+ CUTOFF +C Data from N.F. Allard, March 1997 + IF(FREQ.LT.(82259.10D0-20000.)*2.99792458D10)GO TO 53 + IF(FREQ.GT.(82259.10D0-4000.)*2.99792458D10)GO TO 52 +C TABULATED AT 100 CM-1 SPACING + FREQ15000=(82259.10D0-15000.)*2.99792458D10 + SPACING=100.*2.99792458D10 + IF(FREQ.LT.FREQ15000)THEN + CUTOFF=(CUTOFFH2PLUS(2)-CUTOFFH2PLUS(1))/SPACING* + 1(FREQ-FREQ15000)+CUTOFFH2PLUS(1) + ELSE + ICUT=(FREQ-FREQ15000)/SPACING + CUTFREQ=ICUT*SPACING+FREQ15000 + CUTOFF=(CUTOFFH2PLUS(ICUT+2)-CUTOFFH2PLUS(ICUT+1))/ + 1SPACING*(FREQ-CUTFREQ)+CUTOFFH2PLUS(ICUT+1) + ENDIF +C XNFPH(J,2)=XNFH(J,2) + CUTOFF=(10.**(CUTOFF-14.))/2.99792458D10*XNFPH(J,2) + HPROF4=HPROF4+CUTOFF*1.77245*DOP + GO TO 53 + 52 BETA4000=4000.*2.99792458D10/FO(J)*DBETA + PRQSP4000=SOFBET(BETA4000,PP(J),N,M)*.5/FO(J)*DBETA + CUTOFF4000=(10.**(-11.07-14.))/2.99792458D10*XNFPH(J,2) + HPROF4=HPROF4+CUTOFF4000/PRQSP4000*PRQS/FO(J)*DBETA*1.77245*DOP + 53 F=0. + IF(GAM.GT.0.)F=GAM/3.14159/(GAM**2+BETA**2) + P1=(.9*Y1)**2 + FNS=(P1+.03*SQRT(Y1))/(P1+1.) +C SAME NORMALIZATION AS VOIGT FUNCTION + HPROF4=HPROF4+(PRQS*(1.+FNS)+F)/FO(J)*DBETA*1.77245*DOP + RETURN + END + FUNCTION HFNM(N,M) +C CALCULATES HYDROGEN OSCILLATOR STRENGTHS + DATA NSTR/0/,MSTR/0/ + HFNM=0. + IF(M.LE.N)RETURN + IF(N.EQ.NSTR)GO TO 10 + XN=N + GINF=.2027/XN**.71 + GCA=.124/XN + FKN=XN*1.9603 + WTC=.45-2.4/XN**3*(XN-1.) + NSTR=N + GO TO 15 + 10 IF(M.EQ.MSTR)GO TO 20 + 15 XM=M + XMN=M-N + FK=FKN*(XM/(XMN*(XM+XN)))**3 + XMN12=XMN**1.2 + WT=(XMN12-1.)/(XMN12+WTC) + FNM=FK*(1.-WT*GINF-(.222+GCA/XM)*(1.-WT)) + MSTR=M + 20 HFNM=FNM + RETURN + END + FUNCTION VCSE1F(X) +C ROUGH, BUT ARRANGED TO BE FAST. X.GE.0 + COMMON /EXTAB/EXTAB(1001),EXTABF(1001),E1TAB(2000) + FASTEX(X)=EXTAB(IFIX(X)+1)* + 1EXTABF(IFIX((X-FLOAT(IFIX(X)))*1000.+1.5)) + VCSE1F=0. + IF(X.LE.0.)RETURN + IF(X.GT..01)GO TO 10 + VCSE1F=-ALOG(X)-.577215+X + RETURN + 10 IF(X.GT.1.)GO TO 20 + VCSE1F=-ALOG(X)-.57721566D0+X*(.99999193D0+X*(-.24991055D0+X* + 1 (.05519968D0+X*(-.00976004D0+X*.00107857D0)))) + RETURN + 20 IF(X.GT.30.)RETURN +C VCSE1F=(X*(X+2.334733)+.25062)/(X*(X+3.330657)+1.681534)/X*EXP(-X) + VCSE1F=(X*(X+2.334733D0)+.25062D0)/(X*(X+3.330657D0)+1.681534D0) + 1 /X*FASTEX(X) + RETURN + END + FUNCTION SOFBET(B,P,N,M) +C GENERATES S(BETA,P) FOR HYDROGEN LINES. THE ALPHA AND BETA LINES +C OF THE FIRST THREE SERIES ARE EXPLICITLY INCLUDED AND THE H18 +C PROFILE IS USED FOR THE REST. +C +C THESE PROFILES HAVE BEEN RENORMALIZED TO FULL OSCILLATOR STRENGTH +C +C STORAGE FOR CORRECTIONS (P,BETA,IND),(P,IND),(P,IND) + DIMENSION PROPBM(5,15,7),C(5,7),D(5,7) + DIMENSION PP(5),BETA(15) + DIMENSION PROB1(75),PROB2(75),PROB3(75),PROB4(75),PROB5(75) + DIMENSION PROB6(75),PROB7(75) + DIMENSION C1(5),C2(5),C3(5),C4(5),C5(5),C6(5),C7(5) + DIMENSION D1(5),D2(5),D3(5),D4(5),D5(5),D6(5),D7(5) + EQUIVALENCE (PROPBM(1),PROB1(1)),(PROPBM(76),PROB2(1)) + EQUIVALENCE (PROPBM(151),PROB3(1)),(PROPBM(226),PROB4(1)) + EQUIVALENCE (PROPBM(301),PROB5(1)),(PROPBM(376),PROB6(1)) + EQUIVALENCE (PROPBM(451),PROB7(1)) + EQUIVALENCE (C(1),C1(1)),(C(6),C2(1)),(C(11),C3(1)),(C(16),C4(1)) + EQUIVALENCE (C(21),C5(1)),(C(26),C6(1)),(C(31),C7(1)) + EQUIVALENCE (D(1),D1(1)),(D(6),D2(1)),(D(11),D3(1)),(D(16),D4(1)) + EQUIVALENCE (D(21),D5(1)),(D(26),D6(1)),(D(31),D7(1)) +C LYMAN ALPHA + DATA PROB1/ + 1-.980,-.967,-.948,-.918,-.873,-.968,-.949,-.921,-.879,-.821, + 2-.950,-.922,-.883,-.830,-.764,-.922,-.881,-.830,-.770,-.706, + 3-.877,-.823,-.763,-.706,-.660,-.806,-.741,-.682,-.640,-.625, + 4-.691,-.628,-.588,-.577,-.599,-.511,-.482,-.484,-.514,-.568, + 5-.265,-.318,-.382,-.455,-.531,-.013,-.167,-.292,-.394,-.478, + 6 .166,-.056,-.216,-.332,-.415, .251, .035,-.122,-.237,-.320, + 7 .221, .059,-.068,-.168,-.247, .160, .055,-.037,-.118,-.189, + 8 .110, .043,-.022,-.085,-.147/ + DATA C1 /-18.396, 84.674,-96.273, 3.927, 55.191/ + DATA D1 / 11.801, 9.079, -.651,-11.071,-26.545/ +C LYMAN BETA + DATA PROB2/ + 1-.242, .060, .379, .671, .894, .022, .314, .569, .746, .818, + 2 .273, .473, .605, .651, .607, .432, .484, .489, .442, .343, + 3 .434, .366, .294, .204, .091, .304, .184, .079,-.025,-.135, + 4 .167, .035,-.082,-.189,-.290, .085,-.061,-.183,-.287,-.374, + 5 .032,-.127,-.249,-.344,-.418,-.024,-.167,-.275,-.357,-.420, + 6-.061,-.170,-.257,-.327,-.384,-.047,-.124,-.192,-.252,-.306, + 7-.043,-.092,-.142,-.190,-.238,-.038,-.070,-.107,-.146,-.187, + 8-.030,-.049,-.075,-.106,-.140/ + DATA C2 / 95.740, 18.489, 14.902, 24.466, 42.456/ + DATA D2 / -6.665, -7.136,-10.605,-15.882,-23.632/ +C BALMER ALPHA + DATA PROB3/ + 1-.484,-.336,-.206,-.111,-.058,-.364,-.264,-.192,-.154,-.144, + 2-.299,-.268,-.250,-.244,-.246,-.319,-.333,-.337,-.336,-.337, + 3-.397,-.414,-.415,-.413,-.420,-.456,-.455,-.451,-.456,-.478, + 4-.446,-.441,-.446,-.469,-.512,-.358,-.381,-.415,-.463,-.522, + 5-.214,-.288,-.360,-.432,-.503,-.063,-.196,-.304,-.394,-.468, + 6 .063,-.108,-.237,-.334,-.409, .151,-.019,-.148,-.245,-.319, + 7 .149, .016,-.091,-.177,-.246, .115, .023,-.056,-.126,-.189, + 8 .078, .021,-.036,-.091,-.145/ + DATA C3 /-25.088,145.882,-50.165, 7.902, 51.003/ + DATA D3 / 7.872, 5.592, -2.716,-12.180,-25.661/ +C BALMER BETA + DATA PROB4/ + 1-.082, .163, .417, .649, .829, .096, .316, .515, .660, .729, + 2 .242, .393, .505, .556, .534, .320, .373, .394, .369, .290, + 3 .308, .274, .226, .152, .048, .232, .141, .052,-.046,-.154, + 4 .148, .020,-.094,-.200,-.299, .083,-.070,-.195,-.299,-.385, + 5 .031,-.130,-.253,-.348,-.422,-.023,-.167,-.276,-.359,-.423, + 6-.053,-.165,-.254,-.326,-.384,-.038,-.119,-.190,-.251,-.306, + 7-.034,-.088,-.140,-.190,-.239,-.032,-.066,-.103,-.144,-.186, + 8-.027,-.048,-.075,-.106,-.142/ + DATA C4 / 93.783, 10.066, 9.224, 20.685, 40.136/ + DATA D4 / -5.918, -6.501,-10.130,-15.588,-23.570/ +C PASCHEN ALPHA + DATA PROB5/ + 1-.819,-.759,-.689,-.612,-.529,-.770,-.707,-.638,-.567,-.498, + 2-.721,-.659,-.595,-.537,-.488,-.671,-.617,-.566,-.524,-.497, + 3-.622,-.582,-.547,-.523,-.516,-.570,-.545,-.526,-.521,-.537, + 4-.503,-.495,-.496,-.514,-.551,-.397,-.418,-.448,-.492,-.547, + 5-.246,-.315,-.384,-.453,-.522,-.080,-.210,-.316,-.406,-.481, + 6 .068,-.107,-.239,-.340,-.418, .177,-.006,-.143,-.246,-.324, + 7 .184, .035,-.082,-.174,-.249, .146, .042,-.046,-.123,-.190, + 8 .103, .036,-.027,-.088,-.146/ + DATA C5 /-19.819, 94.981,-79.606, 3.159, 52.106/ + DATA D5 / 10.938, 8.028, -1.267,-11.375,-26.047/ +C PASCHEN BETA + DATA PROB6/ + 1-.073, .169, .415, .636, .809, .102, .311, .499, .639, .710, + 2 .232, .372, .479, .531, .514, .294, .349, .374, .354, .279, + 3 .278, .253, .212, .142, .040, .215, .130, .044,-.051,-.158, + 4 .141, .015,-.097,-.202,-.300, .080,-.072,-.196,-.299,-.385, + 5 .029,-.130,-.252,-.347,-.421,-.022,-.166,-.275,-.359,-.423, + 6-.050,-.164,-.253,-.325,-.384,-.035,-.118,-.189,-.252,-.306, + 7-.032,-.087,-.139,-.190,-.240,-.029,-.064,-.102,-.143,-.185, + 8-.025,-.046,-.074,-.106,-.142/ + DATA C6 /111.107, 11.910, 9.857, 21.371, 41.006/ + DATA D6 / -5.899, -6.381,-10.044,-15.574,-23.644/ +C BALMER 18 + DATA PROB7/ + 1 .005, .128, .260, .389, .504, .004, .109, .220, .318, .389, + 2-.007, .079, .162, .222, .244,-.018, .041, .089, .106, .080, + 3-.026,-.003, .003,-.023,-.086,-.025,-.048,-.087,-.148,-.234, + 4-.008,-.085,-.165,-.251,-.343, .018,-.111,-.223,-.321,-.407, + 5 .032,-.130,-.255,-.354,-.431, .014,-.148,-.269,-.359,-.427, + 6-.005,-.140,-.243,-.323,-.386, .005,-.095,-.178,-.248,-.307, + 7-.002,-.068,-.129,-.187,-.241,-.007,-.049,-.094,-.139,-.186, + 8-.010,-.036,-.067,-.103,-.143/ + DATA C7 /511.318, 1.532, 4.044, 19.266, 41.812/ + DATA D7 / -6.070, -4.528, -8.759,-14.984,-23.956/ + DATA PP/0.,.2,.4,.6,.8/ + DATA BETA/1.,1.259,1.585,1.995,2.512,3.162,3.981,5.012,6.310,7.943 + 1,10.,12.59,15.85,19.95,25.12/ + CORR=1. + B2=B*B + SB=SQRT(B) + IF(B.GT.500.)GO TO 40 + INDX=7 + MMN=M-N + IF(N.LE.3.AND.MMN.LE.2)INDX=2*(N-1)+MMN +C DETERMINE RELEVANT DEBYE RANGE + IM=MIN0(INT(5.*P)+1,4) + IP=IM+1 + WTPP=5.*(P-PP(IM)) + WTPM=1.-WTPP + IF(B.GT.25.12)GO TO 30 + DO 10 J=2,15 + IF(B.LE.BETA(J))GO TO 20 + 10 CONTINUE + 20 JM=J-1 + JP=J + WTBP=(B-BETA(JM))/(BETA(JP)-BETA(JM)) + WTBM=1.-WTBP + CBP=PROPBM(IP,JP,INDX)*WTPP+PROPBM(IM,JP,INDX)*WTPM + CBM=PROPBM(IP,JM,INDX)*WTPP+PROPBM(IM,JM,INDX)*WTPM + CORR=1.+CBP*WTBP+CBM*WTBM +C GET INNER APPROXIMATE PROFILE + PR1=0. + PR2=0. + WT=AMAX1(AMIN1(.5*(10.-B),1.),0.) + IF(B.LE.10.)PR1=8./(83.+(2.+.95*B2)*B) + IF(B.GE.8.)PR2=(1.5/SB+27./B2)/B2 + SOFBET=(PR1*WT+PR2*(1.-WT))*CORR + RETURN +C ASYMPTOTIC PARTS + 30 CC=C(IP,INDX)*WTPP+C(IM,INDX)*WTPM + DD=D(IP,INDX)*WTPP+D(IM,INDX)*WTPM + CORR=1.+DD/(CC+B*SB) + 40 SOFBET=(1.5/SB+27./B2)/B2*CORR + RETURN + END + SUBROUTINE TABVOIGT(VSTEPS,N) + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001) + DIMENSION TABVI(81),TABH1(81) + DATA TABVI/0.,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,1.1,1.2,1.3,1.4,1.5, + 11.6,1.7,1.8,1.9,2.,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.,3.1,3.2, + 2 3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0,4.2,4.4,4.6,4.8,5.0,5.2,5.4,5.6, + 3 5.8,6.0,6.2,6.4,6.6,6.8,7.0,7.2,7.4,7.6,7.8,8.0,8.2,8.4,8.6,8.8, + 4 9.0,9.2,9.4,9.6,9.8,10.0,10.2,10.4,10.6,10.8,11.0,11.2,11.4,11.6, + 5 11.8,12.0/ + DATA TABH1/-1.12838,-1.10596,-1.04048,-.93703,-.80346,-.64945, + 1-.48552,-.32192,-.16772,-.03012,.08594,.17789,.24537,.28981, + 2.31394,.32130,.31573,.30094,.28027,.25648,.231726,.207528,.184882, + 3.164341,.146128,.130236,.116515,.104739,.094653,.086005,.078565, + 4 .072129,.066526,.061615,.057281,.053430,.049988,.046894,.044098, + 5 .041561,.039250,.035195,.031762,.028824,.026288,.024081,.022146, + 6 .020441,.018929,.017582,.016375,.015291,.014312,.013426,.012620, + 7 .0118860,.0112145,.0105990,.0100332,.0095119,.0090306,.0085852, + 8 .0081722,.0077885,.0074314,.0070985,.0067875,.0064967,.0062243, + 9 .0059688,.0057287,.0055030,.0052903,.0050898,.0049006,.0047217, + T .0045526,.0043924,.0042405,.0040964,.0039595/ +C PRETABULATE VOIGT FUNCTION +C 100 STEPS PER DOPPLER WIDTH GIVES 2 PER CENT ACCURACY + DO 1 I=1,N + 1 H0TAB(I)=FLOAT(I-1)/VSTEPS + CALL MAP1(TABVI,TABH1,81,H0TAB,H1TAB,N) + DO 2 I=1,N + VV=(FLOAT(I-1)/VSTEPS)**2 + H0TAB(I)=EXP(-VV) + 2 H2TAB(I)=H0TAB(I)-(VV+VV)*H0TAB(I) + RETURN + END + FUNCTION VOIGT(V,A) +C FAST VOIGT + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001) + IV=V*200.+1.5 + IF(A.LT..2)GO TO 10 + IF(A.GT.1.4)GO TO 2 + IF(A+V.GT.3.2)GO TO 2 + VV=V*V + HH1=H1TAB(IV)+H0TAB(IV)*1.12838 + HH2=H2TAB(IV)+HH1*1.12838-H0TAB(IV) + HH3=(1.-H2TAB(IV))*.37613-HH1*.66667*VV+HH2*1.12838 + HH4=(3.*HH3-HH1)*.37613+H0TAB(IV)*.66667*VV*VV + VOIGT=((((HH4*A+HH3)*A+HH2)*A+HH1)*A+H0TAB(IV))* + 1 (((-.122727278D0*A+.532770573D0)*A-.96284325D0)*A+ + 2 .979895032D0) + RETURN + 2 AA=A*A + VV=V*V + U=(AA+VV)*1.4142 + VOIGT=A*.79788/U + IF(A.GT.100.)RETURN + AAU=AA/U + VVU=VV/U + UU=U*U + VOIGT= + A((((AAU-10.*VVU)*AAU*3.+15.*VVU*VVU)+3.*VV-AA)/UU+1.)*VOIGT +C A ((((AA-10.*VV)*AA*3.+15.*VV*VV)/UU+3.*VV-AA)/UU+1.)*A*.79788/U + RETURN + 10 IF(V.GT.10.)GO TO 12 + 11 VOIGT=(H2TAB(IV)*A+H1TAB(IV))*A+H0TAB(IV) + RETURN + 12 VOIGT=.5642*A/V**2 + RETURN + END + FUNCTION EXPI(N,X) +C LOW PRECISION VERSION 1.E-5 +C EXPONENTIAL INTEGRAL FOR POSITIVE ARGUMENTS AFTER CODY AND +C THACHER, MATH. OF COMP.,22,641(1968) + DATA X1/-1.E20/ + DATA A0,A1,A2,B0,B1/-4.43668255D0,4.42054938D0,3.16274620D0, + 1 7.68641124D0,5.65655216D0/ + DATA C0,C1,C2,D1,D2/.0012102205D0,.98147989D0,.75339742D0, + 1 1.6198645D0,.29135151D0/ + DATA E0,E1,F1/-.9969698D0,-.4257849D0,2.318261D0/ + IF(X.EQ.X1)GO TO 40 + EX=EXP(-X) + X1=X + IF(X.GT.4.)GO TO 10 + IF(X.GT.1.)GO TO 20 + IF(X.GT.0.)GO TO 30 + EX1=0. + GO TO 40 + 10 EX1=(EX+EX*(E0+E1/X)/(X+F1))/X + GO TO 40 + 20 EX1=EX*(C2+(C1+C0*X)*X)/(D2+(D1+X)*X) + GO TO 40 + 30 EX1=(A0+(A1+A2*X)*X)/(B0+(B1+X)*X)-ALOG(X) + 40 EXPI=EX1 + IF(N.EQ.1)RETURN + N1=N-1 + DO 41 I=1,N1 + 41 EXPI=(EX-X*EXPI)/FLOAT(I) + RETURN + END + FUNCTION FASTE1(X) + COMMON /EXTAB/EXTAB(1001),EXTABF(1001),E1TAB(2000) +C DO 3457 I=1,2000 +C3457 E1TAB(I)=EXPI(1,FLOAT(I)*.01) + FASTE1=0. + IF(X.GT.20)RETURN + IF(X.LT..5)GO TO 1 + FASTE1=E1TAB(IFIX(X*100.+.5)) + RETURN + 1 IF(X.LE.0.)RETURN + FASTE1=(1.-.22464*X)*X-ALOG(X)-.57721 + RETURN + END + SUBROUTINE MAP1(XOLD,FOLD,NOLD,XNEW,FNEW,NNEW) + DIMENSION XOLD(1),FOLD(1),XNEW(1),FNEW(1) + L=2 + LL=0 + DO 50 K=1,NNEW + 10 IF(XNEW(K).LT.XOLD(L))GO TO 20 + L=L+1 + IF(L.GT.NOLD)GO TO 30 + GO TO 10 + 20 IF(L.EQ.LL)GO TO 50 + IF(L.EQ.2)GO TO 30 + L1=L-1 + IF(L.GT.LL+1.OR.L.EQ.3)GO TO 21 + CBAC=CFOR + BBAC=BFOR + ABAC=AFOR + IF(L.EQ.NOLD)GO TO 22 + GO TO 25 + 21 L2=L-2 + D=(FOLD(L1)-FOLD(L2))/(XOLD(L1)-XOLD(L2)) + CBAC=FOLD(L)/((XOLD(L)-XOLD(L1))*(XOLD(L)-XOLD(L2)))+ + 1(FOLD(L2)/(XOLD(L)-XOLD(L2))-FOLD(L1)/(XOLD(L)-XOLD(L1)))/ + 2(XOLD(L1)-XOLD(L2)) + BBAC=D-(XOLD(L1)+XOLD(L2))*CBAC + ABAC=FOLD(L2)-XOLD(L2)*D+XOLD(L1)*XOLD(L2)*CBAC + IF(L.LT.NOLD)GO TO 25 + 22 C=CBAC + B=BBAC + A=ABAC + LL=L + GO TO 50 + 25 D=(FOLD(L)-FOLD(L1))/(XOLD(L)-XOLD(L1)) + CFOR=FOLD(L+1)/((XOLD(L+1)-XOLD(L))*(XOLD(L+1)-XOLD(L1)))+ + 1(FOLD(L1)/(XOLD(L+1)-XOLD(L1))-FOLD(L)/(XOLD(L+1)-XOLD(L)))/ + 2(XOLD(L)-XOLD(L1)) + BFOR=D-(XOLD(L)+XOLD(L1))*CFOR + AFOR=FOLD(L1)-XOLD(L1)*D+XOLD(L)*XOLD(L1)*CFOR + WT=0. + IF(ABS(CFOR).NE.0.)WT=ABS(CFOR)/(ABS(CFOR)+ABS(CBAC)) + A=AFOR+WT*(ABAC-AFOR) + B=BFOR+WT*(BBAC-BFOR) + C=CFOR+WT*(CBAC-CFOR) + LL=L + GO TO 50 + 30 IF(L.EQ.LL)GO TO 50 + L=AMIN0(NOLD,L) + C=0. + B=(FOLD(L)-FOLD(L-1))/(XOLD(L)-XOLD(L-1)) + A=FOLD(L)-XOLD(L)*B + LL=L + 50 FNEW(K)=A+(B+C*XNEW(K))*XNEW(K) +C MAP1=LL-1 + RETURN + END + FUNCTION AIRVAC(W) + IMPLICIT REAL*8 (A-H,O-Z) +C W IS AIR WAVELENGTH IN NM +C WAVEN IS AIR WAVENUMBER WHICH IS USUALLY GOOD ENOUGH +C MUST ITERATE FOR EXACT SOLUTION + WAVEN=1.D7/W + WNEW=W*(1.0000834213D0+ + 1 2406030.D0/(1.30D10-WAVEN**2)+15997.D0/(3.89D9-WAVEN**2)) + WAVEN=1.D7/WNEW + WNEW=W*(1.0000834213D0+ + 1 2406030.D0/(1.30D10-WAVEN**2)+15997.D0/(3.89D9-WAVEN**2)) + WAVEN=1.D7/WNEW + AIRVAC=W*(1.0000834213D0+ + 1 2406030.D0/(1.30D10-WAVEN**2)+15997.D0/(3.89D9-WAVEN**2)) + RETURN + END + FUNCTION VACAIR(W) + IMPLICIT REAL*8 (A-H,O-Z) +C W IS VACUUM WAVELENGTH IN NM + WAVEN=1.D7/W + VACAIR=W/(1.0000834213D0+ + 1 2406030.D0/(1.30D10-WAVEN**2)+15997.D0/(3.89D9-WAVEN**2)) + RETURN + END + FUNCTION HE1PROF(J,WAVE,WL,DOPWL,GAMMAR,GAMMAS) +c REAL*8 WAVE,WL,WLSAVE + LINE=WL+1. + IF(LINE.EQ.448)THEN + HE1PROF=HE4471(J,WAVE,WL,DOPWL) + RETURN + ENDIF + IF(LINE.EQ.403)THEN + LINE=WL+.4 + IF(LINE.EQ.402)then + line=wl+1. + go to 8 + endif + HE1PROF=HE4026(J,WAVE,WL,DOPWL) + RETURN + ENDIF + IF(LINE.EQ.439)THEN + HE1PROF=HE4387(J,WAVE,WL,DOPWL) + RETURN + ENDIF + IF(LINE.EQ.493)THEN + HE1PROF=HE4921(J,WAVE,WL,DOPWL) + RETURN + ENDIF + 8 continue +C HE1=381.96, 386.7, 392.6, 400.9, 402.3, 414.37, 416.8 + if(line.eq.382.OR.LINE.EQ.387.or.line.eq.393. + 1 or.line.eq.401.or.LINE.eq.403.or.LINE.EQ.415. + 2 or. LINE.EQ.417)then + he1prof=dimitri(j,wave,wl,DOPWL,GAMMAR,GAMMAS) + RETURN + endif + 1 HE1PROF=GRIEM(J,WAVE,WL,DOPWL,GAMMAR,GAMMAS) + RETURN + END + FUNCTION HE4471(J,WAVE,WL,DOPWL) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,wave +C STANDARD ISOLATED LINE BROADENING PARAMETERS +c From Barnard,Cooper and Smith J.Q.S.R.T. 14,1025,1974 +c ws=we(A),ds=de/we, alfs=alpha; data for Ne=10**13 cm(-3) + DATA WS/0.001460,0.001269,0.001079,0.000898/ + DATA DS/0.036,-0.005,-0.026,-0.034/ + DATA ALFS/0.107,0.119,0.134,0.154/ + DATA DEN/1.0E13/ +c data forb1/?,7.66E-3,9.04E-3,1.01E-2/ ???? + DATA FORB1/0.,0.,0.,0./ + DATA FORB2/0.,0.,0.,0./ +C NANOMETERS + DATA DLP/0.021/ + DATA DLF1/-0.150/ + DATA DLF2/0./ + DATA TS/ 5.0E3,1.0E4,2.0E4,4.0E4/ +c +c 447.1469 -2.198 2.0 169086.864 1.0 191444.603 2.002p 3P 4d 3D +c 447.1469 7 9.262-3.690 NBS 415 5 +c 447.1473 -1.028 2.0 169086.864 2.0 191444.585 2.002p 3P 4d 3D +c 447.1473 7 -3.690 NBS 415 4 +c 447.1473 -0.278 2.0 169086.864 3.0 191444.583 2.002p 3P 4d 3D +c 447.1473 7 -3.690 NBS 415 4 +c 447.1485 -1.028 1.0 169086.940 1.0 191444.603 2.002p 3P 4d 3D +c 447.1485 7 -3.690 NBS 415 4 +c 447.1488 -0.548 1.0 169086.940 2.0 191444.585 2.002p 3P 4d 3D +c 447.1488 7 -3.690 NBS 415 4 +c 447.1682 -0.898 0.0 169087.928 1.0 191444.603 2.002p 3P 4d 3D +c 447.1682 7 -3.690 NBS 415 4 +c +c 447.1498 0.052 4.0 169087.008 7.0 191444.588 2.002p 3P 4d 3D +c 447.1498 7 7.502 -3.69 -7.22NBS 415 4 +c + HE4471=0. + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVE-WL + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,4.0E4) + IFORB=0. + DO 10 IT=1,4 + IF(FORB1(IT).LE.0.)GO TO 10 + IFORB=1 + 10 CONTINUE + IF(E.LE.1.0E13)GO TO 499 + CALL READBCS(1,J,TEMP,XNFHP,XNFHEP,E,DL,PHIHE) + HE4471=1.772453*PHIHE*DOPWL*10. + RETURN + 499 CONTINUE + DO 2 I=2,4 + IT=I + IF(TS(I).GT.TEMP)GO TO 3 + 2 CONTINUE + 3 X=(TEMP-TS(IT-1))/(TS(IT)-TS(IT-1)) + XX=E/DEN +C DAMPING WIDTH + W=XX*(X*WS(IT)+(1.0-X)*WS(IT-1)) +C RATIO OF SHIFT TO WIDTH + D=X*DS(IT)+(1.0-X)*DS(IT-1) +C ION BROADENING PARAMETER + ALF=XX**0.25*(X*ALFS(IT)+(1.0-X)*ALFS(IT-1)) +C FORBIDDEN COMPONENTS INTENSITIES + F1=XX*(X*FORB1(IT)+(1.0-X)*FORB1(IT-1)) + F2=XX*(X*FORB2(IT)+(1.0-X)*FORB2(IT-1)) +C RECIPROCAL PERTURBER VELOCITY +C KM! + XX=XNFHP/E + VM1=8.78E0*(XX+2.0*(1.0-XX))/SQRT(TEMP) +C MEAN INTERPARTICLE DISTANCE + RHOM=1.0/(4.19*E)**(1.0/3.0) +C ION VELOCITY PARAMETER + SIGMA=1.885E14*W*RHOM*VM1/(WL*10.)**2 +C 1.885E14=2*PAI*C*1.E8 1.E8 TRASFORMA RHOM DA CM ad A + X=ALF**(8.0/9.0)/SIGMA**(1.0/3.0) +C TOTAL WIDTHIS ANGSTROMS, then in nm + WTOT=W*(1.0+1.36*X) + wtot=wtot*0.1 +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D*(1.0+2.36*X/ABS(D)) + dtot=dtot*0.1 +C VOIGT PARAMETERS + A=WTOT/DOPWL +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/(DOPPLE(7)*WL4) + con=1. + WWD=WAVE-WL-DTOT +C ALLOW FOR FINE STRUCTURE SPLITTING + HE4471=VOIGT(ABS(WWD-.0184)/DOPWL,A)/9.+ + 1 VOIGT(ABS(WWD+.0013)/DOPWL,A)/12.+ + 2 VOIGT(ABS(WWD+.0010)/DOPWL,A)/4.+ + 3 VOIGT(ABS(WWD+.0029)/DOPWL,A)/180.+ + 4 VOIGT(ABS(WWD+.0025)/DOPWL,A)*11./20. + IF(F1.LE.0) GO TO 4 +C ADD IN P TO F FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF1)/DOPWL + HE4471=HE4471+F1*VOIGT(V,A) + IF(F2.LE.0) GO TO 4 +C ADD P TO G FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF2)/DOPWL + HE4471=HE4471+F2*VOIGT(V,A) + 4 HE4471=HE4471*CON +C TYPE*,HE4471 + RETURN + END + FUNCTION HE4026(J,WAVE,WL,DOPWL) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,WAVE +C STANDARD ISOLATED LINE BROADENING PARAMETERS +C FROM GRIEM 1974-book-Spectral line broadening by plasmas +c ws=we(A), ds=de/we, alfs=alpha, data for Ne=10**16 cm(-3) + DATA WS/4.04,3.49,2.96,2.47/ + DATA DS/0.1339,0.0960,0.0780,0.0709/ + DATA ALFS/0.969,1.083,1.225,1.403/ + DATA DEN/1.0E16/ +C +c data forb1/?,7.66E-3,9.04E-3,1.01E-2/ ???? + DATA FORB1/0.,0.,0.,0./ + DATA FORB2/0.,0.,0.,0./ +C NANOMETERS + DATA DLP/0.017/ + DATA DLF1/-0.080/ + DATA DLF2/0./ + DATA TS/5000.,10000.,20000.,40000./ +c +c 402.6185 -2.620 2.0 169086.864 1.0 193917.253 2.002p 3P 5d 3D +c 402.6185 7 7.216-3.160-7.070NBS 423 4 +c 402.6187 -1.450 2.0 169086.864 2.0 193917.243 2.002p 3P 5d 3D +c 402.6187 7 -3.160 NBS 423 4 +c 402.6187 -0.700 2.0 169086.864 3.0 193917.243 2.002p 3P 5d 3D +c 402.6187 7 -3.160 NBS 423 4 +c 402.6198 -1.450 1.0 169086.940 1.0 193917.253 2.002p 3P 5d 3D +c 402.6198 7 -3.160 NBS 423 4 +c 402.6199 -0.970 1.0 169086.940 2.0 193917.243 2.002p 3P 5d 3D +c 402.6199 7 -3.160 NBS 423 4 +c 402.6358 -1.320 0.0 169087.928 1.0 193917.253 2.002p 3P 5d 3D +c 402.6358 7 -3.160 NBS 423 4 +c +c 402.6210 -0.370 4.0 169087.008 7.0 193917.245 2.002p 3P 5d 3D +c 402.6210 7 7.216 -3.16 -7.07NBS 415 4 +c + HE4026=0. + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVE-WL + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,4.0E4) + IFORB=0. + DO 10 IT=1,4 + IF(FORB1(IT).LE.0.)GO TO 10 + IFORB=1 + 10 CONTINUE + IF(E.LE.1.E14)GO TO 499 + CALL READBCS(2,J,TEMP,XNFHP,XNFHEP,E,DL,PHIHE) + HE4026=1.772453*PHIHE*DOPWL*10. + RETURN + 499 CONTINUE + DO 2 I=2,4 + IT=I + IF(TS(I).GT.TEMP)GO TO 3 + 2 CONTINUE + 3 X=(TEMP-TS(IT-1))/(TS(IT)-TS(IT-1)) + XX=E/DEN +C DAMPING WIDTH + W=XX*(X*WS(IT)+(1.0-X)*WS(IT-1)) +C RATIO OF SHIFT TO WIDTH + D=X*DS(IT)+(1.0-X)*DS(IT-1) +C ION BROADENING PARAMETER + ALF=XX**0.25*(X*ALFS(IT)+(1.0-X)*ALFS(IT-1)) +C FORBIDDEN COMPONENTS INTENSITIES + F1=XX*(X*FORB1(IT)+(1.0-X)*FORB1(IT-1)) + F2=XX*(X*FORB2(IT)+(1.0-X)*FORB2(IT-1)) +C RECIPROCAL PERTURBER VELOCITY +C KM! + XX=XNFHP/E + VM1=8.78E0*(XX+2.0*(1.0-XX))/SQRT(TEMP) +C MEAN INTERPARTICLE DISTANCE + RHOM=1.0/(4.19*E)**(1.0/3.0) +C ION VELOCITY PARAMETER + SIGMA=1.885E14*W*RHOM*VM1/(WL*10.)**2 +C 1.885E14=2*PAI*C*1.E8 1.E8 TRASFORMA RHOM DA CM ad A + X=ALF**(8.0/9.0)/SIGMA**(1.0/3.0) +C TOTAL WIDTHIS ANGSTROMS, then in nm + WTOT=W*(1.0+1.36*X) + wtot=wtot*0.1 +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D*(1.0+2.36*X/ABS(D)) + dtot=dtot*0.1 +C VOIGT PARAMETERS + A=WTOT/DOPWL +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/(DOPPLE(7)*WL4) + con=1. + WWD=WAVE-WL-DTOT +C ALLOW FOR FINE STRUCTURE SPLITTING + HE4026=VOIGT(ABS(WWD-.0148)/DOPWL,A)/9.+ + 1 VOIGT(ABS(WWD+.0012)/DOPWL,A)/12.+ + 2 VOIGT(ABS(WWD+.0011)/DOPWL,A)/4.+ + 3 VOIGT(ABS(WWD+.0025)/DOPWL,A)/180.+ + 4 VOIGT(ABS(WWD+.0023)/DOPWL,A)*11./20. + IF(F1.LE.0) GO TO 4 +C ADD IN P TO F FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF1)/DOPWL + HE4026=HE4026+F1*VOIGT(V,A) + IF(F2.LE.0) GO TO 4 +C ADD P TO G FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF2)/DOPWL + HE4026=HE4026+F2*VOIGT(V,A) + 4 HE4026=HE4026*CON +C TYPE*,HE4026 + RETURN + END + FUNCTION HE4387(J,WAVE,WL,DOPWL) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,WAVE +C STANDARD ISOLATED LINE BROADENING PARAMETERS +C FROM GRIEM + DATA WS/6.13,5.15,4.24,3.45/ + DATA DS/0.411,0.363,0.325,0.293/ + DATA ALFS/1.159,1.321,1.527,1.783/ + DATA DEN/1.0E16/ +C +c data forb1/?,7.66E-3,9.04E-3,1.01E-2/ ???? + DATA FORB1/0.,0.,0.,0./ + DATA FORB2/0.,0.,0.,0./ +C NANOMETERS + DATA DLP/0.017/ + DATA DLF1/-0.080/ + DATA DLF2/0./ + DATA TS/5000.,10000.,20000.,40000./ +c 438.7929 -0.883 1.0 171135.000 2.0 193918.391 2.002p 1P 5d 1D +c 438.7929 7 9.258-3.080 NBS 524 5 + HE4387=0. + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVE-WL + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,4.0E4) + IFORB=0. + DO 10 IT=1,4 + IF(FORB1(IT).LE.0.)GO TO 10 + IFORB=1 + 10 CONTINUE + IF(E.LE.1.E14)GO TO 499 + CALL READBCS(3,J,TEMP,XNFHP,XNFHEP,E,DL,PHIHE) + HE4387=1.772453*PHIHE*DOPWL*10. + RETURN + 499 CONTINUE + DO 2 I=2,4 + IT=I + IF(TS(I).GT.TEMP)GO TO 3 + 2 CONTINUE + 3 X=(TEMP-TS(IT-1))/(TS(IT)-TS(IT-1)) + XX=E/DEN +C DAMPING WIDTH + W=XX*(X*WS(IT)+(1.0-X)*WS(IT-1)) +C RATIO OF SHIFT TO WIDTH + D=X*DS(IT)+(1.0-X)*DS(IT-1) +C ION BROADENING PARAMETER + ALF=XX**0.25*(X*ALFS(IT)+(1.0-X)*ALFS(IT-1)) +C FORBIDDEN COMPONENTS INTENSITIES + F1=XX*(X*FORB1(IT)+(1.0-X)*FORB1(IT-1)) + F2=XX*(X*FORB2(IT)+(1.0-X)*FORB2(IT-1)) +C RECIPROCAL PERTURBER VELOCITY +C KM! + XX=XNFHP/E + VM1=8.78E0*(XX+2.0*(1.0-XX))/SQRT(TEMP) +C MEAN INTERPARTICLE DISTANCE + RHOM=1.0/(4.19*E)**(1.0/3.0) +C ION VELOCITY PARAMETER + SIGMA=1.885E14*W*RHOM*VM1/(WL*10.)**2 +C 1.885E14=2*PAI*C*1.E8 1.E8 TRASFORMA RHOM DA CM ad A + X=ALF**(8.0/9.0)/SIGMA**(1.0/3.0) +C TOTAL WIDTHIS ANGSTROMS, then in nm + WTOT=W*(1.0+1.36*X) + wtot=wtot*0.1 +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D*(1.0+2.36*X/ABS(D)) + dtot=dtot*0.1 +C VOIGT PARAMETERS + A=WTOT/DOPWL +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/(DOPWL) + con=1. + HE4387=VOIGT(ABS(WAVE-WL-DTOT)/DOPWL,A) + RETURN + END + FUNCTION HE4921(J,WAVE,WL,DOPWL) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,wave +C STANDARD ISOLATED LINE BROADENING PARAMETERS +c From Barnard,Cooper and Smith J.Q.S.R.T. 15, 429, 1975 +c ws=we(A),ds=de/we, alfs=alpha; data for Ne=10**13 cm(-3) + DATA WS/0.002312,0.001963,0.001624,0.001315/ + DATA DS/0.3932,0.3394,0.2950,0.2593/ + DATA ALFS/0.1207,0.1365,0.1564,0.1844/ + DATA DEN/1.0E13/ +c data forb1/?,7.66E-3,9.04E-3,1.01E-2/ ???? + DATA FORB1/0.,0.,0.,0./ + DATA FORB2/0.,0.,0.,0./ +C NANOMETERS + DATA DLP/0.021/ + DATA DLF1/-0.150/ + DATA DLF2/0./ + DATA TS/ 5.0E3,1.0E4,2.0E4,4.0E4/ +c +c 492.1931 -0.435 1.0 171135.000 2.0 191446.559 2.002p 1P 4d 1D +c 492.1931 7 9.262-3.600 NBS 516 5 +c + HE4921=0. + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVE-WL + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,4.0E4) + IFORB=0. + DO 10 IT=1,4 + IF(FORB1(IT).LE.0.)GO TO 10 + IFORB=1 + 10 CONTINUE + IF(E.LE.1.0E13)GO TO 499 + CALL READBCS(4,J,TEMP,XNFHP,XNFHEP,E,DL,PHIHE) + HE4921=1.772453*PHIHE*DOPWL*10. + RETURN + 499 CONTINUE + DO 2 I=2,4 + IT=I + IF(TS(I).GT.TEMP)GO TO 3 + 2 CONTINUE + 3 X=(TEMP-TS(IT-1))/(TS(IT)-TS(IT-1)) + XX=E/DEN +C DAMPING WIDTH + W=XX*(X*WS(IT)+(1.0-X)*WS(IT-1)) +C RATIO OF SHIFT TO WIDTH + D=X*DS(IT)+(1.0-X)*DS(IT-1) +C ION BROADENING PARAMETER + ALF=XX**0.25*(X*ALFS(IT)+(1.0-X)*ALFS(IT-1)) +C FORBIDDEN COMPONENTS INTENSITIES + F1=XX*(X*FORB1(IT)+(1.0-X)*FORB1(IT-1)) + F2=XX*(X*FORB2(IT)+(1.0-X)*FORB2(IT-1)) +C RECIPROCAL PERTURBER VELOCITY +C KM! + XX=XNFHP/E + VM1=8.78E0*(XX+2.0*(1.0-XX))/SQRT(TEMP) +C MEAN INTERPARTICLE DISTANCE + RHOM=1.0/(4.19*E)**(1.0/3.0) +C ION VELOCITY PARAMETER + SIGMA=1.885E14*W*RHOM*VM1/(WL*10.)**2 +C 1.885E14=2*PAI*C*1.E8 1.E8 TRASFORMA RHOM DA CM ad A + X=ALF**(8.0/9.0)/SIGMA**(1.0/3.0) +C TOTAL WIDTHIS ANGSTROMS, then in nm + WTOT=W*(1.0+1.36*X) + wtot=wtot*0.1 +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D*(1.0+2.36*X/ABS(D)) + dtot=dtot*0.1 +C VOIGT PARAMETERS + A=WTOT/DOPWL +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/(DOPPLE(7)*WL4) + con=1. + WWD=WAVE-WL-DTOT +C ALLOW FOR FINE STRUCTURE SPLITTING + HE4921=VOIGT(ABS(WWD)/DOPWL,A) + IF(F1.LE.0) GO TO 4 +C ADD IN P TO F FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF1)/DOPWL + HE4921=HE4921+F1*VOIGT(V,A) + IF(F2.LE.0) GO TO 4 +C ADD P TO G FORBIDDEN COMPONENT + V=ABS((WAVE-WL)-DLF2)/DOPWL + HE4921=HE4921+F2*VOIGT(V,A) + 4 HE4921=HE4921*CON +C TYPE*,HE4921 + RETURN + END + FUNCTION GRIEM(J,WAVESYN,WL,DOPWL,GAMMAR,GAMMAS) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,WAVESYN + DIMENSION WIDTH(4,999),SHIFT(4,999),ALPHA(4,999),BETA(4,999) + DIMENSION CODE(999),XNELOG(999),TTAB(4,999),WAVE(999) + CHARACTER*10 TRANS(999) + CHARACTER*60 GRIEM0200(210) + DATA (GRIEM0200(I),I= 1, 5)/ + 1' 2.00 52.2213 1s1S-4p1P 16.0 ', + 2' 5000. .0179 -.0112 .275 .000016 ', + 3' 10000. .0168 -.00872 .290 .000032 ', + 4' 20000. .0152 -.00647 .311 .000064 ', + 5' 40000. .0135 -.00460 .341 .00013 '/ + DATA (GRIEM0200(I),I= 6, 10)/ + 1' 2.00 53.7030 1s1S-3p1P 16.0 ', + 2' 5000. .00432 -.00277 .153 .000054 ', + 3' 10000. .00409 -.00216 .160 .00011 ', + 4' 20000. .00378 -.00159 .169 .00022 ', + 5' 40000. .00341 -.00111 .183 .00043 '/ + DATA (GRIEM0200(I),I= 11, 15)/ + 1' 2.00 58.4334 1s1S-2p1P 16.0 ', + 2' 5000. .000121 -.0000299 .012 .0024 ', + 3' 10000. .000158 -.00000567 .010 .0049 ', + 4' 20000. .000199 .0000220 .008 .0098 ', + 5' 40000. .000234 .0000460 .007 .020 '/ + DATA (GRIEM0200(I),I= 16, 20)/ + 1' 2.00 282.9076 2s3S-6p3P 16.0 ', + 2' 5000. 1.79 1.07 .285 .0014 ', + 3' 10000. 1.87 .829 .276 .0027 ', + 4' 20000. 1.84 .619 .279 .0054 ', + 5' 40000. 1.72 .455 .294 .011 '/ + DATA (GRIEM0200(I),I= 21, 25)/ + 1' 2.00 294.5104 2s3S-5p3P 16.0 ', + 2' 5000. .808 .522 .204 .0030 ', + 3' 10000. .857 .412 .195 .0059 ', + 4' 20000. .857 .311 .195 .012 ', + 5' 40000. .811 .231 .203 .024 '/ + DATA (GRIEM0200(I),I= 26, 30)/ + 1' 2.00 318.7746 2s3S-4p3P 16.0 ', + 2' 5000. .313 .219 .134 .0087 ', + 3' 10000. .338 .176 .127 .017 ', + 4' 20000. .344 .134 .125 .035 ', + 5' 40000. .332 .101 .128 .069 '/ + DATA (GRIEM0200(I),I= 31, 35)/ + 1' 2.00 388.8649 2s3S-3p3P 16.0 ', + 2' 5000. .102 .0744 .075 .050 ', + 3' 10000. .112 .0603 .070 .099 ', + 4' 20000. .117 .0464 .067 .20 ', + 5' 40000. .117 .0348 .067 .40 '/ + DATA (GRIEM0200(I),I= 36, 40)/ + 1' 2.00 396.4729 2s1S-4p1P 16.0 ', + 2' 5000. 1.03 -.697 .275 .0071 ', + 3' 10000. .996 -.504 .290 .014 ', + 4' 20000. .877 -.374 .311 .028 ', + 5' 40000. .776 -.266 .341 .056 '/ + DATA (GRIEM0200(I),I= 41, 45)/ + 1' 2.00 402.6187 2p3P-5d3D 16.0 ', + 2' 5000. 4.04 .541 .969 .0016 ', + 3' 10000. 3.49 .335 1.083 .0033 ', + 4' 20000. 2.96 .231 1.225 .0066 ', + 5' 40000. 2.47 .175 1.403 .013 '/ + DATA (GRIEM0200(I),I= 46, 50)/ + 1' 2.00 412.0811 2p3P-5s3S 16.0 ', + 2' 5000. .785 .897 .171 .013 ', + 3' 10000. .897 .894 .155 .026 ', + 4' 20000. .984 .808 .144 .052 ', + 5' 40000. 1.01 .670 .141 .10 '/ + DATA (GRIEM0200(I),I= 51, 55)/ + 1' 2.00 438.7929 2p1P-5d1D 16.0 ', + 2' 5000. 6.13 2.52 1.159 .0017 ', + 3' 10000. 5.15 1.87 1.321 .0033 ', + 4' 20000. 4.24 1.38 1.527 .0067 ', + 5' 40000. 3.45 1.01 1.783 .013 '/ + DATA (GRIEM0200(I),I= 56, 60)/ + 1' 2.00 443.7551 2p1P-5s1S 16.0 ', + 2' 5000. 1.41 1.51 .199 .012 ', + 3' 10000. 1.57 1.43 .184 .024 ', + 4' 20000. 1.65 1.24 .177 .047 ', + 5' 40000. 1.62 .996 .179 .094 '/ + DATA (GRIEM0200(I),I= 61, 65)/ + 1' 2.00 447.1488 2p3P-4d3D 16.0 ', + 2' 5000. 1.44 .136 .589 .0058 ', + 3' 10000. 1.26 .0804 .650 .012 ', + 4' 20000. 1.09 .0614 .726 .023 ', + 5' 40000. .927 .0549 .819 .047 '/ + DATA (GRIEM0200(I),I= 66, 70)/ + 1' 2.00 471.3139 2p3P-4s3S 16.0 ', + 2' 5000. .342 .402 .115 .044 ', + 3' 10000. .393 .416 .103 .088 ', + 4' 20000. .437 .390 .095 .18 ', + 5' 40000. .459 .335 .092 .35 '/ + DATA (GRIEM0200(I),I= 71, 75)/ + 1' 2.00 492.1931 2p1P-4d1D 16.0 ', + 2' 5000. 2.30 1.02 .683 .0061 ', + 3' 10000. 1.96 .773 .773 .012 ', + 4' 20000. 1.63 .584 .885 .024 ', + 5' 40000. 1.35 .440 1.023 .049 '/ + DATA (GRIEM0200(I),I= 76, 80)/ + 1' 2.00 501.5678 2s1S-3p1P 16.0 ', + 2' 5000. .378 -.250 .154 .044 ', + 3' 10000. .359 -.200 .160 .088 ', + 4' 20000. .334 -.152 .169 .18 ', + 5' 40000. .306 -.111 .180 .35 '/ + DATA (GRIEM0200(I),I= 81, 85)/ + 1' 2.00 504.7738 2p1P-4s1S 16.0 ', + 2' 5000. .625 .699 .134 .038 ', + 3' 10000. .705 .685 .123 .077 ', + 4' 20000. .756 .611 .117 .15 ', + 5' 40000. .760 .504 .116 .31 '/ + DATA (GRIEM0200(I),I= 86, 90)/ + 1' 2.00 587.5615 2p3P-3d3D 16.0 ', + 2' 5000. .159 -.0881 .064 .23 ', + 3' 10000. .170 -.0553 .061 .46 ', + 4' 20000. .176 -.0256 .059 .92 ', + 5' 40000. .177 -.00504 .059 1.8 '/ + DATA (GRIEM0200(I),I= 91, 95)/ + 1' 2.00 667.8154 2p1P-3d1D 16.0 ', + 2' 5000. .423 .275 .146 .14 ', + 3' 10000. .386 .233 .157 .27 ', + 4' 20000. .349 .196 .169 .54 ', + 5' 40000. .318 .161 .181 1.1 '/ + DATA (GRIEM0200(I),I= 96,100)/ + 1' 2.00 706.5176 2p3P-3s3S 16.0 ', + 2' 5000. .180 .215 .067 .44 ', + 3' 10000. .207 .231 .060 .87 ', + 4' 20000. .235 .227 .055 1.7 ', + 5' 40000. .254 .203 .052 3.5 '/ + DATA (GRIEM0200(I),I=101,105)/ + 1' 2.00 728.1349 2p1P-3s1S 16.0 ', + 2' 5000. .320 .374 .081 .33 ', + 3' 10000. .365 .382 .073 .65 ', + 4' 20000. .403 .355 .068 1.3 ', + 5' 40000. .419 .303 .066 2.6 '/ + DATA (GRIEM0200(I),I=106,110)/ + 1' 2.00 836.1694 3s3S-6p3P 16.0 ', + 2' 5000. 15.6 9.33 .285 .035 ', + 3' 10000. 16.3 7.23 .276 .070 ', + 4' 20000. 16.1 5.38 .279 .14 ', + 5' 40000. 15.0 3.93 .293 .28 '/ + DATA (GRIEM0200(I),I=111,115)/ + 1' 2.00 946.3596 3s3S-5p3P 16.0 ', + 2' 5000. 8.34 5.30 .203 .099 ', + 3' 10000. 8.86 4.13 .194 .20 ', + 4' 20000. 8.88 3.06 .194 .40 ', + 5' 40000. 8.46 2.21 .201 .79 '/ + DATA (GRIEM0200(I),I=116,120)/ + 1' 2.00 960.3418 3s1S-6p1P 16.0 ', + 2' 5000. 40.7 -24.0 .614 .023 ', + 3' 10000. 37.3 -18.5 .656 .045 ', + 4' 20000. 33.2 -13.7 .716 .091 ', + 5' 40000. 28.8 -9.94 .796 .18 '/ + DATA (GRIEM0200(I),I=121,125)/ + 1' 2.00 1066.7641 3p3P-6s3S 16.0 ', + 2' 5000. 12.9 14.2 .238 .12 ', + 3' 10000. 14.7 13.9 .217 .23 ', + 4' 20000. 16.0 12.4 .203 .46 ', + 5' 40000. 16.5 10.3 .198 .92 '/ + DATA (GRIEM0200(I),I=126,130)/ + 1' 2.00 1083.0336 2s3S-2p3P 16.0 ', + 2' 5000. .0493 -.0518 .028 8.3 ', + 3' 10000. .0601 -.0557 .024 17. ', + 4' 20000. .0755 -.0540 .020 33. ', + 5' 40000. .0931 -.0465 .017 66. '/ + DATA (GRIEM0200(I),I=131,135)/ + 1' 2.00 1099.6561 3d3D-6p3P 16.0 ', + 2' 5000. 27.0 16.3 .286 .08 ', + 3' 10000. 28.3 12.8 .276 .16 ', + 4' 20000. 27.9 9.59 .279 .32 ', + 5' 40000. 26.1 7.10 .294 .64 '/ + DATA (GRIEM0200(I),I=136,140)/ + 1' 2.00 1101.3070 3s1S-5p1P 16.0 ', + 2' 5000. 23.1 -14.3 .429 .066 ', + 3' 10000. 21.4 -11.2 .454 .13 ', + 4' 20000. 19.3 -8.49 .491 .26 ', + 5' 40000. 16.9 -6.24 .541 .53 '/ + DATA (GRIEM0200(I),I=141,145)/ + 1' 2.00 1104.5003 3p1P-6d1D 16.0 ', + 2' 5000. 94.5 37.7 1.742 .013 ', + 3' 10000. 79.0 27.8 1.992 .026 ', + 4' 20000. 64.8 20.4 2.310 .052 ', + 5' 40000. 52.5 14.9 2.706 .10 '/ + DATA (GRIEM0200(I),I=146,150)/ + 1' 2.00 1196.9060 3p3P-5d3D 16.0 ', + 2' 5000. 35.8 4.42 .967 .043 ', + 3' 10000. 31.0 2.56 1.077 .086 ', + 4' 20000. 26.4 1.65 1.215 .17 ', + 5' 40000. 22.1 1.20 1.387 .34 '/ + DATA (GRIEM0200(I),I=151,155)/ + 1' 2.00 1252.7537 3s3S-4p3P 16.0 ', + 2' 5000. 4.83 3.09 .130 .54 ', + 3' 10000. 5.26 2.34 .122 1.1 ', + 4' 20000. 5.45 1.65 .119 2.1 ', + 5' 40000. 5.37 1.12 .120 4.3 '/ + DATA (GRIEM0200(I),I=156,160)/ + 1' 2.00 1284.5935 3p3P-5s3S 16.0 ', + 2' 5000. 7.61 8.33 .165 .40 ', + 3' 10000. 8.85 8.26 .147 .81 ', + 4' 20000. 9.85 7.44 .136 1.6 ', + 5' 40000. 10.3 6.17 .132 3.2 '/ + DATA (GRIEM0200(I),I=161,165)/ + 1' 2.00 1296.8439 3p1P-5d1D 16.0 ', + 2' 5000. 54.3 23.1 1.149 .043 ', + 3' 10000. 45.9 17.3 1.304 .086 ', + 4' 20000. 38.1 12.8 1.501 .17 ', + 5' 40000. 31.1 9.45 1.745 .34 '/ + DATA (GRIEM0200(I),I=166,170)/ + 1' 2.00 1298.4882 3d3D-5p3P 16.0 ', + 2' 5000. 15.8 10.4 .205 .25 ', + 3' 10000. 16.8 8.24 .196 .51 ', + 4' 20000. 16.8 6.27 .195 1.0 ', + 5' 40000. 16.0 4.65 .203 2.0 '/ + DATA (GRIEM0200(I),I=171,175)/ + 1' 2.00 1508.3656 3s1S-4p1P 16.0 ', + 2' 5000. 15.1 -10.2 .277 .38 ', + 3' 10000. 14.3 -8.34 .289 .77 ', + 4' 20000. 13.2 -6.55 .307 1.5 ', + 5' 40000. 11.9 -4.95 .331 3.1 '/ + DATA (GRIEM0200(I),I=176,180)/ + 1' 2.00 1700.2364 3p3P-4d3D 16.0 ', + 2' 5000. 21.2 .993 .577 .32 ', + 3' 10000. 18.9 .216 .629 .64 ', + 4' 20000. 16.7 .0624 .692 1.3 ', + 5' 40000. 14.5 .129 .769 2.6 '/ + DATA (GRIEM0200(I),I=181,185)/ + 1' 2.00 1868.5313 3d3D-4f3F 16.0 ', + 2' 5000. 16.2 -2.77 .650 .50 ', + 3' 10000. 13.7 -1.22 .734 .99 ', + 4' 20000. 11.8 -.300 .823 2.0 ', + 5' 40000. 10.2 .131 .914 4.0 '/ + DATA (GRIEM0200(I),I=186,190)/ + 1' 2.00 1869.7233 3d1D-4f1F 16.0 ', + 2' 5000. 18.6 -5.39 .756 .42 ', + 3' 10000. 15.6 -3.35 .862 .84 ', + 4' 20000. 13.2 -1.95 .978 1.7 ', + 5' 40000. 11.3 -1.10 1.101 3.4 '/ + DATA (GRIEM0200(I),I=191,195)/ + 1' 2.00 1908.9369 3p1P-4d1D 16.0 ', + 2' 5000. 37.3 18.1 .657 .35 ', + 3' 10000. 32.3 14.0 .732 .71 ', + 4' 20000. 27.4 10.7 .828 1.4 ', + 5' 40000. 23.0 8.03 .945 2.8 '/ + DATA (GRIEM0200(I),I=196,200)/ + 1' 2.00 1954.3172 3d3D-4p3P 16.0 ', + 2' 5000. 12.2 8.91 .137 1.9 ', + 3' 10000. 13.2 7.25 .128 3.9 ', + 4' 20000. 13.6 5.56 .126 7.7 ', + 5' 40000. 13.3 4.11 .128 15. '/ + DATA (GRIEM0200(I),I=201,205)/ + 1' 2.00 2058.1299 2s1S-2p1P 16.0 ', + 2' 5000. .364 -.412 .038 33. ', + 3' 10000. .433 -.430 .033 65. ', + 4' 20000. .514 -.400 .029 130. ', + 5' 40000. .590 -.332 .026 260. '/ + DATA (GRIEM0200(I),I=206,210)/ + 1' 2.00 2112.0002 3p3P-4s3S 16.0 ', + 2' 5000. 7.17 6.48 .089 4.6 ', + 3' 10000. 8.76 6.78 .077 9.2 ', + 4' 20000. 10.1 6.47 .069 18. ', + 5' 40000. 10.9 5.62 .065 37. '/ + DATA IREAD/0/ + IF(IREAD.EQ.1)GO TO 10 +C DO 5 ILINE=1,999 +C READ(21,3,END=7)CODE(iline),WAVE(iline),TRANS(iline), +C 1 XNELOG(iline), (TTAB(I,iline),WIDTH(I,iline),SHIFT(I,iline), +C 2 ALPHA(I,iline), BETA(I,iline),I=1,4) +C PRINT 3,CODE(iline),WAVE(iline),TRANS(iline), +C 1 XNELOG(iline), (TTAB(I,iline),WIDTH(I,iline),SHIFT(I,iline), +C 2 ALPHA(I,iline), BETA(I,iline),I=1,4) +C 3 FORMAT(F5.2,F10.4,A10,F5.1,F10.0,4F10.6/30X,F10.0,4F10.6/ +C 2 30X,F10.0,4F10.6/30X,F10.0,4F10.6) +C +C 2.00 52.2213 1s1S-4p1P 16.0 5000. .0179 -.0112 .275 .000016 +C 10000. .0168 -.00872 .290 .000032 +C 20000. .0152 -.00647 .311 .000064 +C 40000. .0135 -.00460 .341 .00013 +C +C 5 CONTINUE +C 7 NLINE=ILINE-1 + NLINE=42 + DO 5 ILINE=1,NLINE + READ(GRIEM0200(ILINE*5-4),3)CODE(iline),WAVE(iline),TRANS(iline), + 1XNELOG(iline) + 3 FORMAT(F5.2,F10.4,A10,F5.1,F10.0,4F10.6) + DO 5 I=1,4 + READ(GRIEM0200(ILINE*5-4+I),4)TTAB(I,iline),WIDTH(I,iline), + 1 SHIFT(I,iline),ALPHA(I,iline), BETA(I,iline) + 5 SHIFT(I,ILINE)=SHIFT(I,ILINE)/WIDTH(I,ILINE) + 4 FORMAT(F10.0,4F10.6) + IREAD=1 + 10 CONTINUE + DO 11 ILINE=1,NLINE + IF(WL.GT.WAVE(ILINE)-.1.AND.WL.LT.WAVE(ILINE)+.1)GO TO 12 + 11 CONTINUE + V=ABS(WAVESYN-WL)/DOPWL +C WANT DOPPLE(7) BUT DOES NOT INCLUDE ISOTOPIC MASS DOPWL/WL DOES + A=(GAMMAR+GAMMAS*XNE(J))/(DOPWL/WL) + GRIEM=VOIGT(V,A) + RETURN + 12 IL=ILINE + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVESYN-WL + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,8.0E4) +c DO 2 I=2,5 + DO 2 I=2,4 + IT=I + IF(Ttab(It,il).GT.TEMP)GO TO 33 + 2 CONTINUE + 33 X=(TEMP-Ttab(IT-1,il))/(Ttab(IT,il)-Ttab(IT-1,il)) + XX=E/(10.**xnelog(il)) +C DAMPING WIDTH + W=XX*(X*Width(IT,il)+(1.0-X)*Width(IT-1,il)) +C RATIO OF SHIFT TO WIDTH + D=X*shift(IT,il)+(1.0-X)*shift(IT-1,il) +C ION BROADENING PARAMETER + ALF=XX**0.25*(X*ALpha(IT,il)+(1.0-X)*ALpha(IT-1,il)) +C RECIPROCAL PERTURBER VELOCITY +C KM! + XX=XNFHP/E + VM1=8.78E0*(XX+2.0*(1.0-XX))/SQRT(TEMP) +C MEAN INTERPARTICLE DISTANCE + RHOM=1.0/(4.19*E)**(1.0/3.0) +C ION VELOCITY PARAMETER + SIGMA=1.885E14*W*RHOM*VM1/(WL*10.)**2 +C 1.885E14=2*PAI*C*1.E8 1.E8 TRASFORMA RHOM DA CM ad A + X=ALF**(8.0/9.0)/SIGMA**(1.0/3.0) +C TOTAL WIDTHIS ANGSTROMS, then in nm + WTOT=W*(1.0+1.36*X) + wtot=wtot*0.1 +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D*(1.0+2.36*X/ABS(D)) + dtot=dtot*0.1 +C VOIGT PARAMETERS +C WANT DOPPLE(7) FOR GAMMAR BUT DOES NOT INCLUDE ISOTOPIC MASS DOPWL/WL DOES + A=WTOT/DOPWL+GAMMAR/(DOPWL/WL) +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/DOPWL + con=1. + WWD=WAVESYN-WL-DTOT + GRIEM=VOIGT(ABS(WWD)/DOPWL,A) + RETURN + END + FUNCTION DIMITRI(J,WAVEsyn,WL,DOPWL,GAMMAR,GAMMAS) + PARAMETER (kw=99) + DIMENSION WS(4),DS(4),ALFS(4),FORB1(4),FORB2(4),TS(4) + COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw),PTOTAL(kw) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + COMMON /BHE/BHE1(kw,29),AHE1(kw),SHE1(kw),BHE2(kw,6),AHE2(kw), + 1 SHE2(kw),AHEMIN(kw),SIGHE(kw),XNFPHE(kw,3),XNFHE(kw,2) + COMMON /BHYD/BHYD(kw,8),AHYD(kw),SHYD(kw),AH2P(kw),BMIN(kw), + 1 AHMIN(kw),SHMIN(kw),SIGH(kw),SIGH2(kw),AHLINE(kw), + 2 SHLINE(kw),XNFPH(kw,2),XNFH(kw) +c REAL*8 WL,wavesyn + DIMENSION WIDTH(4,999),SHIFT(4,999),WIDTHP(4,999),SHIFTP(4,999), + 1 WIDTHHE(4,999),SHIFTHE(4,999) + dimension CODE(999),XNELOG(999),TTAB(4,999),wave(999),gam(999) + DIMENSION XNETAB(7) + CHARACTER*10 TRANS(999) + CHARACTER*61 DIMITRI0200(35) + DATA (DIMITRI0200(I),I= 1, 5)/ + 1' 2.00 381.9624 2p3P-6d3D 13.0 ', + 2' 5000. 1.47E-02-2.04E-04 9.31E-03 7.87E-03 0.00E+00 0.00E+00', + 3' 10000. 1.29E-02-4.41E-05 1.06E-02 9.31E-03 8.90E-03 7.72E-03', + 4' 20000. 1.10E-02-1.37E-04 1.72E-02 1.13E-02 1.02E-02 8.99E-03', + 5' 40000. 9.33E-03-2.31E-05 1.08E-02 1.32E-02 1.13E-02 1.06E-02'/ +c DATA (DIMITRI0200(I),I= 6, 10)/ +c 1' 2.00 386.7494 2p3P-6s3S 13.0 ', +c 2' 5000. 2.74E-03 2.09E-03 6.15E-04 5.72E-04 5.29E-04 4.90E-04', +c 3' 10000. 2.91E-03 2.02E-03 6.91E-04 6.48E-04 5.93E-04 5.55E-04', +c 4' 20000. 2.85E-03 1.58E-03 7.75E-04 7.29E-04 6.66E-04 6.25E-04', +c 5' 40000. 3.10E-03 1.10E-03 8.71E-04 8.20E-04 7.48E-04 7.04E-04'/ +c DATA (DIMITRI0200(I),I= 6, 10)/ +c 1' 2.00 386.7494 2p3P-6s3S 13.0 ', +c 2' 5000. 2.74E-03 1.79E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', +c 3' 10000. 2.91E-03 1.86E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', +c 4' 20000. 2.85E-03 1.51E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', +c 5' 40000. 3.10E-03 1.10E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00'/ + DATA (DIMITRI0200(I),I= 6, 10)/ + 1' 2.00 386.7494 2p3P-6s3S 16.0 ', + 2' 5000. 2.74E+00 1.79E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 3' 10000. 2.91E+00 1.86E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 4' 20000. 2.85E+00 1.51E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 5' 40000. 3.10E+00 1.10E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00'/ + DATA (DIMITRI0200(I),I= 11, 15)/ + 1' 2.00 392.6544 2p1P-6d1D 13.0 ', + 2' 5000 7.07E-02 7.13E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 3' 10000. 6.07E-02 4.27E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 4' 20000. 5.09E-02 1.96E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 5' 40000. 4.19E-02 5.60E-04 0.00E+00 0.00E+00 0.00E+00 0.00E+00'/ + DATA (DIMITRI0200(I),I= 16, 20)/ + 1' 2.00 400.9256 2p1P-7d1D 13.0 ', + 2' 5000 3.96E-02 6.23E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 3' 10000. 3.42E-02 2.78E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 4' 20000. 2.88E-02 1.85E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 5' 40000. 2.38E-02 1.07E-03 3.05E-02 3.67E-02 0.00E+00 0.00E+00'/ + DATA (DIMITRI0200(I),I= 21, 25)/ + 1' 2.00 402.40 2p1P-7s1S 13.0 ', + 2' 5000 9.24E-03 6.25E-03 2.09E-03 1.91E-03 1.79E-03 1.63E-03', + 3' 10000. 8.80E-03 5.30E-03 2.34E-03 2.17E-03 2.01E-03 1.86E-03', + 4' 20000. 9.18E-03 3.77E-03 2.63E-03 2.46E-03 2.42E-03 2.72E-03', + 5' 40000. 9.19E-03 2.32E-03 2.96E-03 2.77E-03 2.54E-03 2.38E-03'/ + DATA (DIMITRI0200(I),I= 26, 30)/ + 1' 2.00 414.3759 2p1P-6d1D 13.0 ', + 2' 5000. 2.21E-02 3.63E-03 0.00E+00 0.00E+00 0.00E+00 0.00E+00', + 3' 10000. 1.90E-02 2.37E-03 1.85E-02 1.67E-02 0.00E+00 0.00E+00', + 4' 20000. 1.59E-02 1.12E-03 1.87E-02 1.96E-02 0.00E+00 0.00E+00', + 5' 40000. 1.31E-02 6.75E-04 1.65E-02 2.16E-02 1.89E-02 1.90E-02'/ + DATA (DIMITRI0200(I),I= 31, 35)/ + 1' 2.00 416.8972 2p1P-6s1S 13.0 ', + 2' 5000. 4.83E-03 3.41E-03 1.08E-03 1.00E-03 9.29E-04 8.57E-04', + 3' 10000. 4.72E-03 3.06E-03 1.21E-03 1.13E-03 1.04E-03 9.72E-04', + 4' 20000. 4.84E-03 2.27E-03 1.36E-03 1.28E-03 1.17E-03 1.10E-03', + 5' 40000. 4.97E-03 1.55E-03 1.53E-03 1.44E-03 1.32E-03 1.24E-03'/ + DATA IREAD/0/ + wl4=sngl(wl) + wavesyn4=sngl(wavesyn) + IF(IREAD.EQ.1)GO TO 10 +c shift=d +cc DO 5 ILINE=1,69 +cc read(20,134)n,wave(iline),gam(iline) +134 format(i5,f10.1,f10.3) +cc k=1 +cc read(20,133)xnetab(k) +c type*,k,xnetab(k) +133 FORMAT(10X,E7.1) +C READ(20,3)TRANS(i,iline),(TTAB(I,iline),WIDTH(I,iline), +C 1 SHIFT(I,iline),WIDTHP(I,iline),SHIFTP(I,iline), +C 2 WIDTHHE(I,iline), SHIFTHE(I,iline),i=1,6) +C 3 FORMAT(5x,A10,F8.0,6(1X,E10.2)/15x,f8.0,6(1x,e10.2)/15x, +C 1 f8.0,6(1x,e10.2)/15x,f8.0,6(1x,e10.2)/15x,f8.0, +C 2 6(1x,e10.2)/15x,f8.0,6(1x,e10.2)) +C xnelog(iline)=alog10(xnetab(k)) + NLINE=7 + DO 5 ILINE=1,NLINE + READ(DIMITRI0200(ILINE*5-4),3)CODE(ILINE),WAVE(ILINE), + 1 TRANS(ILINE), XNELOG(ILINE) + 3 FORMAT(F5.2,F10.4,A10,F5.1,31x) + DO 5 I=1,4 + READ(DIMITRI0200(ILINE*5-4+I),4)TTAB(I,ILINE),WIDTH(I,ILINE), + 1 SHIFT(I,ILINE),WIDTHP(I,ILINE),SHIFTP(I,ILINE), + 2 WIDTHHE(I,ILINE),SHIFTHE(I,ILINE) + 4 FORMAT(F7.0,6E9.2) + SHIFT(I,ILINE)=SHIFT(I,ILINE)/WIDTH(I,ILINE) + IF(WIDTHP(I,ILINE).NE.0.) + 1 SHIFTP(I,ILINE)=SHIFTP(I,ILINE)/WIDTHP(I,ILINE) + IF(WIDTHHE(I,ILINE).NE.0.) + 1 SHIFTHE(I,ILINE)=SHIFTHE(I,ILINE)/WIDTHHE(I,ILINE) + 5 continue + IREAD=1 + 10 continue + DO 11 ILINE=1,NLINE + IF(WL.GT.WAVE(ILINE)-.1.AND.WL.LT.WAVE(ILINE)+.1)GO TO 12 + 11 CONTINUE + V=ABS(WAVESYN-WL)/DOPWL +C WANT DOPPLE(7) BUT DOES NOT INCLUDE ISOTOPIC MASS DOPWL/WL DOES + A=(GAMMAR+GAMMAS*XNE(J))/(DOPWL/WL) + DIMITRI=VOIGT(V,A) + RETURN + 12 IL=ILINE + E=XNE(J) + TEMP=T(J) +C NUMERO PROTONI=N(H+)/U(H+) U(H+)=1 + XNFHP=XNFPH(J,2) + XNFHEP=XNFHE(J,2) + DL=WAVEsyn4-WL4 + TEMP=AMAX1(TEMP,5.0E3) + TEMP=AMIN1(TEMP,8.0E4) + DO 2 I=2,4 + IT=I + IF(Ttab(It,il).GT.TEMP)GO TO 33 + 2 CONTINUE + 33 X=(TEMP-Ttab(IT-1,il))/(Ttab(IT,il)-Ttab(IT-1,il)) + 333 XX=E/(10.**xnelog(il)) + XXH=XNFHP/(10.**XNELOG(IL)) + XXHE=XNFHEP/(10.**XNELOG(IL)) +C DAMPING WIDTHs + W=XX*(X*Width(IT,il)+(1.0-X)*Width(IT-1,il)) + Wp=XXH*(X*Widthp(IT,il)+(1.0-X)*Widthp(IT-1,il)) + Whe=XXHE*(X*Widthhe(IT,il)+(1.0-X)*Widthhe(IT-1,il)) +C RATIO SHIFTS TO WIDTH + D=shift(IT,il)+(1.0-X)*shift(IT-1,il) + Dp=X*shiftp(IT,il)+(1.0-X)*shiftp(IT-1,il) + Dhe=X*shifthe(IT,il)+(1.0-X)*shifthe(IT-1,il) +C TOTAL WIDTH IN ANGSTROMS, then in nm + WTOT=W+Wp+Whe + wtot=wtot*0.1 +c half-half-width + wtot=wtot/2. +C TOTAL SHIFT IN ANGSTROM, then in nm + DTOT=W*D+Wp*dp+Whe*dhe + dtot=dtot*0.1 +C VOIGT PARAMETERS +C WANT DOPPLE(7) BUT DOES NOT INCLUDE ISOTOPIC MASS DOPWL/WL DOES + A=WTOT/DOPWL+GAMMAR/(DOPWL/WL) +C NORMALIZATION CONSTANT +c CON=0.564189583547756E0/(1.0+F1+F2)/(DOPPLE(7)*WL4) + con=1. + WWD=WAVEsyn-WL-DTOT + dimitri=VOIGT(ABS(WWD)/DOPWL,A) + RETURN + END + SUBROUTINE READBCS(LINE,J,TEMP,XNFHP,XNFHEP,XNE,DLNM,PHIHE) +C LINE=1 4471 +C LINE=2 4026 +C LINE=3 4387 +C LINE=4 4921 + PARAMETER (kw=99) + COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),HCKT(kw),ITEMP + CHARACTER*8 TITLE1,TITLE2 + DIMENSION DLAM(204,4),PHI(8),PHIHP(4,7,142),PHIHEP(4,7,142) + DIMENSION PHILAM(204),PHI4026(4,8,196),PHI4387(4,8,204) + DIMENSION PHI4921(4,7,142) + DIMENSION NDLAM(4),NXNE(4),XNE1(4) + DATA NDLAM/142,196,204,142/ + DATA NXNE/7,8,8,7/ + DATA XNE1/13.,14.,14.,13./ + DATA ITEMP1/0/ + IF(DLAM(1,1).EQ.-150.)GO TO 10 + OPEN(UNIT=18,FORM='FORMATTED',STATUS='OLD',READONLY,SHARED) +C 4471 + READ(18,1) TITLE1 + READ(18,1) TITLE2 + 1 FORMAT(A80) + DO 22 IL=1,142 + DO 20 NE=1,7 + READ(18,34)FNE,DWL,(PHI(I),I=1,8) + 34 FORMAT(1X,F5.1,F8.2,8F7.3) +c TYPE*,FNE,DWL,(PHI(I),I=1,8) + DO 21 IT=1,4 + PHIHP(IT,NE,IL)=PHI(IT) + 21 PHIHEP(IT,NE,IL)=PHI(IT+4) + 20 CONTINUE + 22 DLAM(IL,1)=DWL-150. +C 4026 + READ(18,1)TITLE1 + READ(18,1)TITLE2 + DO 32 IL=1,196 + DO 30 NE=1,8 + READ(18,34)FNE,DWL,(PHI4026(IT,NE,IL),IT=1,4) +C TYPE*,FNE,DWL,(PHI4026(IT,NE,IL),IT=1,4) + 30 CONTINUE + 32 DLAM(IL,2)=DWL-150. +C 4387 + READ(18,1)TITLE1 + READ(18,1)TITLE2 + DO 35 IL=1,204 + DO 33 NE=1,8 + READ(18,34)FNE,DWL,(PHI4387(IT,NE,IL),IT=1,4) + 33 CONTINUE + 35 DLAM(IL,3)=DWL-150. +C 4921 + READ(18,1) TITLE1 + READ(18,1) TITLE2 + DO 45 IL=1,142 + DO 43 NE=1,7 + READ(18,34)FNE,DWL,(PHI(I),I=1,8) +c TYPE*,FNE,DWL,(PHI(I),I=1,8) + DO 46 IT=1,4 + PHIHP(IT,NE,IL)=PHI(IT) + 46 PHIHEP(IT,NE,IL)=PHI(IT+4) + 43 CONTINUE + 45 DLAM(IL,4)=DWL-150. + CLOSE(UNIT=18) + JSAVE=0 +C + 10 IF(J*LINE.EQ.JSAVE)GO TO 550 +C TEMPERATURE AND IONS DENSITY INTERPOLATION + AT=LOG10(TEMP) + BT=(AT-3.698970)/.3010300+1. + IT=BT+0.00001 + IT=MAX(MIN(IT,3),1) + WT=BT-IT + AP=LOG10(XNE) + AP=MAX(XNE1(LINE),AP) + BP=(AP-XNE1(LINE))/0.5+1. + IP=BP+0.00001 + IP=MAX(MIN(IP,NXNE(LINE)-1),1) + WP=BP-IP + C1W1W=(1.-WP)*(1.-WT) + C1WW=(1.-WP)*WT + CW1W=WP*(1.-WT) + CWW=WP*WT + GO TO (410,420,430,440),LINE + 410 XXH=XNFHP/XNE + XXHE=XNFHEP/XNE + DO 411 I=1,NDLAM(LINE) + 411 PHILAM(I)=XXH*10.**(C1W1W*PHIHP(IT ,IP ,I)+ + 1 C1WW*PHIHP(IT+1,IP ,I)+ + 2 CW1W*PHIHP(IT ,IP+1,I)+ + 3 CWW*PHIHP(IT+1,IP+1,I))+ + 4 XXHE*10.**(C1W1W*PHIHEP(IT ,IP ,I)+ + 5 C1WW*PHIHEP(IT+1,IP ,I)+ + 6 CW1W*PHIHEP(IT ,IP+1,I)+ + 7 CWW*PHIHEP(IT+1,IP+1,I)) + GO TO 502 + 420 DO 421 I=1,NDLAM(LINE) + 421 PHILAM(I)=10.**(C1W1W*PHI4026(IT ,IP ,I)+ + 1 C1WW*PHI4026(IT+1,IP ,I)+ + 2 CW1W*PHI4026(IT ,IP+1,I)+ + 3 CWW*PHI4026(IT+1,IP+1,I)) + GO TO 502 + 430 DO 431 I=1,NDLAM(LINE) + 431 PHILAM(I)=10.**(C1W1W*PHI4387(IT ,IP ,I)+ + 1 C1WW*PHI4387(IT+1,IP ,I)+ + 2 CW1W*PHI4387(IT ,IP+1,I)+ + 3 CWW*PHI4387(IT+1,IP+1,I)) + GO TO 502 + 440 XXH=XNFHP/XNE + XXHE=XNFHEP/XNE + DO 441 I=1,NDLAM(LINE) + 441 PHILAM(I)=XXH*10.**(C1W1W*PHIHP(IT ,IP ,I)+ + 1 C1WW*PHIHP(IT+1,IP ,I)+ + 2 CW1W*PHIHP(IT ,IP+1,I)+ + 3 CWW*PHIHP(IT+1,IP+1,I))+ + 4 XXHE*10.**(C1W1W*PHIHEP(IT ,IP ,I)+ + 5 C1WW*PHIHEP(IT+1,IP ,I)+ + 6 CW1W*PHIHEP(IT ,IP+1,I)+ + 7 CWW*PHIHEP(IT+1,IP+1,I)) + 502 CALL INTEG(DLAM(1,LINE),PHILAM,PHINORM,NDLAM(LINE),0.) + DO 503 I=1,NDLAM(LINE) + 503 PHILAM(I)=LOG10(PHILAM(I)/PHINORM) + JSAVE=J*LINE +C +C NOW DLAM INTERPOLATION + 550 DL=DLNM*10. + DO 600 I=2,NDLAM(LINE) + IF(DL.GT.DLAM(I,LINE).AND.I.LT.NDLAM(LINE))GO TO 600 + A=(DLAM(I,LINE)-DL)/(DLAM(I,LINE)-DLAM(I-1,LINE)) + B=(DL-DLAM(I-1,LINE))/(DLAM(I,LINE)-DLAM(I-1,LINE)) + PHIHE=10.**(A*PHILAM(I-1)+B*PHILAM(I)) + RETURN + 600 CONTINUE + RETURN + END + SUBROUTINE INTEG(X,F,FINT,N,START) +C SUBROUTINE INTEG(X,F,FINT,N) + DIMENSION X(1),F(1) + DIMENSION A(1000),B(1000),C(1000) + CALL PARCOE(F,X,A,B,C,N) + FINT=START +C FINT(1)=(A(1)+(B(1)/2.+C(1)/3.*X(1))*X(1))*X(1) +C FINT(2)=(A(1)+(B(1)/2.+C(1)/3.*X(2))*X(2))*X(2) +C IF(N.EQ.2)RETURN + N1=N-1 +C DO 10 I=2,N1 + DO 10 I=1,N1 + 10 FINT=FINT+(A(I)+B(I)/2.*(X(I+1)+X(I))+ + 1C(I)/3.*((X(I+1)+X(I))*X(I+1)+X(I)*X(I)))*(X(I+1)-X(I)) + RETURN + END + SUBROUTINE PARCOE(F,X,A,B,C,N) + DIMENSION F(1),X(1),A(1),B(1),C(1) + C(1)=0. + B(1)=(F(2)-F(1))/(X(2)-X(1)) + A(1)=F(1)-X(1)*B(1) + N1=N-1 + C(N)=0. + B(N)=(F(N)-F(N1))/(X(N)-X(N1)) + A(N)=F(N)-X(N)*B(N) + IF(N.EQ.2)RETURN + DO 1 J=2,N1 + J1=J-1 + D=(F(J)-F(J1))/(X(J)-X(J1)) + C(J)=F(J+1)/((X(J+1)-X(J))*(X(J+1)-X(J1)))-F(J)/((X(J)-X(J1))* + 1(X(J+1)-X(J)))+F(J1)/((X(J)-X(J1))*(X(J+1)-X(J1))) + B(J)=D-(X(J)+X(J1))*C(J) + 1 A(J)=F(J1)-X(J1)*D+X(J)*X(J1)*C(J) + C(2)=0. + B(2)=(F(3)-F(2))/(X(3)-X(2)) + A(2)=F(2)-X(2)*B(2) + C(3)=0. + B(3)=(F(4)-F(3))/(X(4)-X(3)) + A(3)=F(3)-X(3)*B(3) + DO 2 J=2,N1 + IF(C(J).EQ.0.)GO TO 2 + J1=J+1 + WT=ABS(C(J1))/(ABS(C(J1))+ABS(C(J))) + A(J)=A(J1)+WT*(A(J)-A(J1)) + B(J)=B(J1)+WT*(B(J)-B(J1)) + C(J)=C(J1)+WT*(C(J)-C(J1)) + 2 CONTINUE + A(N1)=A(N) + B(N1)=B(N) + C(N1)=C(N) + RETURN + END + SUBROUTINE BLOCKE +c IMPLICIT REAL*8 (A-H,O-Z) + PARAMETER (kw=99) + COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99) +C Grevesse,N. and Anders, E. 1988. presented at the workshop +C on the "Solar Interior and Atmosphere", Tucson, Nov 15-18. +C Anders, E. and Grevesse, N. 1989 Geochimica et Cosmochimica Acta, +C vol. 53, pp. 197-214. +C H has been defined to be -.04 instead of 12 +C 1H 2HE + DATA ABUND/ 0.911,0.089, +C 3LI 4BE 5B 6C 7N 8O 9F 10NE + 1-10.88,-10.89, -9.44, -3.48, -3.99, -3.11, -7.48, -3.95, +C 11NA 12MG 13AL 14SI 15P 16S 17CL 18AR + 2 -5.71, -4.46, -5.57, -4.49, -6.59, -4.83, -6.54, -5.48, +C 19K 20CA 21SC 22TI 23V 24CR 25MN 26FE + 3 -6.92, -5.68, -8.94, -7.05, -8.04, -6.37, -6.65, -4.37, +C 27CO 28NI 29CU 30ZN 31GA 33GE 33AS 34SE + 4 -7.12, -5.79, -7.83, -7.44, -9.16, -8.63, -9.67, -8.69, +C 35BR 36KR 37RB 38SR 39Y 40ZR 41NB 42MO + 5 -9.41, -8.81, -9.44, -9.14, -9.80, -9.44,-10.62,-10.12, +C 43TC 44RU 45RH 46PD 47AG 48CD 49IN 50SN + 6-20.00,-10.20,-10.92,-10.35,-11.10,-10.18,-10.38,-10.04, +C 51SB 52TE 53I 54XE 55CS 56BA 57LA 58CE + 7-11.04, -9.80,-10.53, -9.81,-10.92, -9.91,-10.82,-10.49, +C 59PR 60ND 61PM 62SM 63EU 64GD 65TB 66DY + 8-11.33,-10.54,-20.00,-11.04,-11.53,-10.92,-12.14,-10.94, +C 67HO 68ER 69TM 70YB 71LU 72HF 73TA 74W + 9-11.78,-11.11,-12.04,-10.96,-11.28,-11.16,-11.91,-10.93, +C 75RE 76OS 77IR 78PT 79AU 80HG 81TL 82PB + T-11.77,-10.59,-10.69,-10.24,-11.03,-10.95,-11.14,-10.19, +C 83BI 84PO 85AT 86RN 87FR 88RA 89AC 90TH + 1-11.33,-20.00,-20.00,-20.00,-20.00,-20.00,-20.00,-11.92, +C 91PA 92U 93NP 94PU 95AM 96CM 97BK 98CF 99ES + 2-20.00,-12.51,-20.00,-20.00,-20.00,-20.00,-20.00,-20.00,-20.00/ + DATA ATMASS/ 1.008,4.003, + 1 6.939,9.013,10.81,12.01,14.01,16.00,19.00,20.18,22.99,24.31, + 2 26.98,28.09,30.98,32.07,35.45,39.95,39.10,40.08,44.96,47.90, + 3 50.94,52.00,54.94,55.85,58.94,58.71,63.55,65.37,69.72,72.60, + 4 74.92,78.96,79.91,83.80,85.48,87.63,88.91,91.22,92.91,95.95, + 5 99.00,101.1,102.9,106.4,107.9,112.4,114.8,118.7,121.8,127.6, + 6 126.9,131.3,132.9,137.4,138.9,140.1,140.9,144.3,147.0,150.4, + 7 152.0,157.3,158.9,162.5,164.9,167.3,168.9,173.0,175.0,178.5, + 8 181.0,183.9,186.3,190.2,192.2,195.1,197.0,200.6,204.4,207.2, + 9 209.0,210.0,211.0,222.0,223.0,226.1,227.1,232.0,231.0,238.0, + T 237.0,244.0,243.0,247.0,247.0,251.0,254.0/ + DATA ELEM/ 2HH , 2HHE, + 1 2HLI, 2HBE, 2HB , 2HC , 2HN , 2HO , 2HF , 2HNE, 2HNA, 2HMG, + 2 2HAL, 2HSI, 2HP , 2HS , 2HCL, 2HAR, 2HK , 2HCA, 2HSC, 2HTI, + 3 2HV , 2HCR, 2HMN, 2HFE, 2HCO, 2HNI, 2HCU, 2HZN, 2HGA, 2HGE, + 4 2HAS, 2HSE, 2HBR, 2HKR, 2HRB, 2HSR, 2HY , 2HZR, 2HNB, 2HMO, + 5 2HTC, 2HRU, 2HRH, 2HPD, 2HAG, 2HCD, 2HIN, 2HSN, 2HSB, 2HTE, + 6 2HI , 2HXE, 2HCS, 2HBA, 2HLA, 2HCE, 2HPR, 2HND, 2HPM, 2HSM, + 7 2HEU, 2HGD, 2HTB, 2HDY, 2HHO, 2HER, 2HTM, 2HYB, 2HLU, 2HHF, + 8 2HTA, 2HW , 2HRE, 2HOS, 2HIR, 2HPT, 2HAU, 2HHG, 2HTL, 2HPB, + 9 2HBI, 2HPO, 2HAT, 2HRN, 2HFR, 2HRA, 2HAC, 2HTH, 2HPA, 2HU , + T 2HNP, 2HPU, 2HAM, 2HCM, 2HBK, 2HCF, 2HES/ + RETURN + END |