diff options
Diffstat (limited to 'dfsynthe/dfsynthe.for')
| -rw-r--r-- | dfsynthe/dfsynthe.for | 3733 |
1 files changed, 3733 insertions, 0 deletions
diff --git a/dfsynthe/dfsynthe.for b/dfsynthe/dfsynthe.for new file mode 100644 index 0000000..dbaba37 --- /dev/null +++ b/dfsynthe/dfsynthe.for @@ -0,0 +1,3733 @@ + PROGRAM DFSYNTHE +c revised 28jul00 +C THIS PROGRAM IS REAL*4 EXCEPT FOR WAVELENGTHS AND ENERGY LEVELS +C +C TAPE5 INPUT reads IFTEFF +C TAPE6 OUTPUT +C TAPE10 input from XNFDF +C TAPE11 input line data from LOWLINESDF +C TAPE21 input line data from HIGHLINESDF +C TAPE31 input line data from DIATOMICSDF +C TAPE41 input line data from TIOLINESDF +C TAPE43 input line data from H2OLINESDF +C TAPE51 input line data from NLTELINESA for XLINOP +C TAPE61 input line data from NLTELINESB simple lines +C TAPE12 line data selected for each temperature from 11,21,31,41 +C TAPE14 temporary file of the 0 km/s opacity spectrum +C TAPE15 output file of distribution functions +C TAPE19 line data for each temperature selected from 51 +C TAPE22 input from XNFDF + PARAMETER (kw=25) + PARAMETER (LENBUFF=3510000,MAXPROF=10000) + IMPLICIT REAL*4 (A-H,O-Z) + COMMON /BUFFER/BUFFER(LENBUFF),PROFILE(MAXPROF),B10000(10000) + COMMON /CONTIN/CONTINUUM(LENBUFF) + 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(10001),E1TAB(2000) + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + COMMON /NBIGLIT/NBEG(1540),NEND(1540),NSTEPS(13,1540),NN(1540) + 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 /VTPROF/VTPROF(120,20),VTCENTER(20),NVT(20) + COMMON /XNFDOP/XNFPEL(594),DOPPLE(594),XNFDOP(594) + DIMENSION ABLOG(3,377),IFTP(kw),ABMIN(kw),MLINEJ(kw) + DIMENSION TITLE(74),IDMOL(60),MOMASS(60) + REAL*8 FREQEDGE(377),WLEDGE(377),CMEDGE(377),CONFRQ(1131) + REAL*8 DELEDGE(377),HALFEDGE(377) + REAL*8 AEDGE(377),BEDGE(377),CEDGE(377) + real*8 ablog8(3,377) + REAL*4 KAPPA0 + REAL*4 FREQ1000(4000) + REAL*8 TITLE,TEFF,GLOG,IDMOL,MOMASS + CHARACTER*9 NAME15(95) + COMMON /TAPE1112/ITAPE12(6,10000) + REAL*4 TAPE12(6,10000) + EQUIVALENCE (TAPE12(1,1),ITAPE12(1,1)) + REAL*8 RESOLU,RATIO,RATIOLG,SIGMA2,WLBEG,WLEND,WBEGIN,WLBEG0 + REAL*8 EMERGE + REAL*8 WL,E,EP,WLVAC,CENTER,CONCEN,WAVE + REAL*8 LABEL,LABELP,OTHER1,OTHER2 + INTEGER IFVT(20),IFTEFF(57) +C 0 KM/S + 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20 KM/S +C DATA IFVT/1,1,0,1,0,0,0,1,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/ + DATA IFVT/1,1,0,1,0,0,0,1,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/ + DATA IFTEFF/57*1/ +C DATA IFTEFF/57*0/ + DATA NAME15/ + 1'DFT01.DAT','DFT02.DAT','DFT03.DAT','DFT04.DAT','DFT05.DAT', + 2'DFT06.DAT','DFT07.DAT','DFT08.DAT','DFT09.DAT','DFT10.DAT', + 3'DFT11.DAT','DFT12.DAT','DFT13.DAT','DFT14.DAT','DFT15.DAT', + 4'DFT16.DAT','DFT17.DAT','DFT18.DAT','DFT19.DAT','DFT20.DAT', + 5'DFT21.DAT','DFT22.DAT','DFT23.DAT','DFT24.DAT','DFT25.DAT', + 6'DFT26.DAT','DFT27.DAT','DFT28.DAT','DFT29.DAT','DFT30.DAT', + 7'DFT31.DAT','DFT32.DAT','DFT33.DAT','DFT34.DAT','DFT35.DAT', + 8'DFT36.DAT','DFT37.DAT','DFT38.DAT','DFT39.DAT','DFT40.DAT', + 9'DFT41.DAT','DFT42.DAT','DFT43.DAT','DFT44.DAT','DFT45.DAT', + A'DFT46.DAT','DFT47.DAT','DFT48.DAT','DFT49.DAT','DFT50.DAT', + 1'DFT51.DAT','DFT52.DAT','DFT53.DAT','DFT54.DAT','DFT55.DAT', + 2'DFT56.DAT','DFT57.DAT','DFT58.DAT','DFT59.DAT','DFT60.DAT', + 3'DFT61.DAT','DFT62.DAT','DFT63.DAT','DFT64.DAT','DFT65.DAT', + 4'DFT66.DAT','DFT67.DAT','DFT68.DAT','DFT69.DAT','DFT70.DAT', + 5'DFT71.DAT','DFT72.DAT','DFT73.DAT','DFT74.DAT','DFT75.DAT', + 6'DFT76.DAT','DFT77.DAT','DFT78.DAT','DFT79.DAT','DFT80.DAT', + 7'DFT81.DAT','DFT82.DAT','DFT83.DAT','DFT84.DAT','DFT85.DAT', + 8'DFT86.DAT','DFT87.DAT','DFT88.DAT','DFT89.DAT','DFT90.DAT', + 9'DFT91.DAT','DFT92.DAT','DFT93.DAT','DFT94.DAT','DFT95.DAT'/ +C + READ(5,5)IFTEFF + 5 FORMAT(57I1) + WRITE(6,6)IFTEFF + 6 FORMAT(1X,57I1) +C CALL BEGTIME + RESOLU=500000.D0 + RATIO=1.D0+1.D0/RESOLU + RATIOLG=LOG(RATIO) + WLBEG=8.97666 + WLEND=10000. + WLBEG0=WLBEG/RATIO + IXWLBEG=DLOG(WLBEG)/RATIOLG + WBEGIN=DEXP(IXWLBEG*RATIOLG) + IF(WBEGIN.LT.WLBEG)THEN + IXWLBEG=IXWLBEG+1 + WBEGIN=DEXP(IXWLBEG*RATIOLG) + ENDIF + DO 17 I=1,3510 + 17 FREQ1000(I)=2.99792458E17/WLBEG/RATIO**(I*1000-999) + CUTOFF=.001 +C + OPEN(UNIT=10,STATUS='OLD',FORM='UNFORMATTED',READONLY,SHARED) + OPEN(UNIT=22,STATUS='OLD',FORM='UNFORMATTED',READONLY,SHARED) +C + DO 3456 I=1,10001 + 3456 EXTAB(I)=EXP(-(I-1)*.01) + 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) + CALL DFINTERVALS + CALL VTTAB +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,(FREQEDGE(IEDGE),WLEDGE(IEDGE),CMEDGE(IEDGE), + 1IEDGE=1,NEDGE),IDMOL,MOMASS + READ(10)NCON,(CONFRQ(NU),NU=1,NCON) + REWIND 10 + 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 CALL ENDTIME + OPEN(UNIT=15,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) + OPEN(UNIT=16,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) + OPEN(UNIT=17,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) + OPEN(UNIT=18,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) + OPEN(UNIT=20,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) +C + ITEMP=0 + DO 950 ITEFF=1,57 + ITEMP=ITEMP+1 +C OPEN(UNIT=15,TYPE='NEW',FORM='UNFORMATTED',NAME=NAME15(ITEFF), +C 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8) + READ(10) + READ(10) + READ(10) + IF(IFTEFF(ITEFF).EQ.0)THEN + READ(10)T(1) + TLOG10=LOG10(T(1)) + WRITE(6,202)ITEFF,TLOG10 + 202 FORMAT(//I3,' LOGT',F5.2) + DO 201 J=1,NRHOX + 201 READ(10) + READ(22) + READ(22) + READ(22) + READ(22) + GO TO 950 + ENDIF + READ(10)T,TKEV,TK,HKT,TLOG,HCKT,P,XNE,XNATOM,RHO,RHOX,VTURB, + 1XNFH,XNFHE,XNFH2 +C CALL BEGTIME + CALL SELECTLINES(N12,N19) +C CALL ENDTIME +C + IREC=0 + DO 500 J=1,NRHOX +C CALL BEGTIME + TLOG10=LOG10(T(J)) + PLOG10=LOG10(P(J)) + IVT=0 + WRITE(6,209)TLOG10,PLOG10,IVT + 209 FORMAT(//' LOGT',F5.2,' LOGP',F5.2,' VT',I3) + REWIND 12 +C INITIALIZE BUFFER + LENGTH=3507859 + DO 210 NBUFF=1,LENBUFF + 210 BUFFER(NBUFF)=0. + READ(10)XNFPEL,DOPPLE,ABLOG + DO 2002 IEDGE=1,NEDGE-1 + ablog8(1,iedge)=ablog(1,iedge) + ablog8(2,iedge)=ablog(2,iedge) + ablog8(3,iedge)=ablog(3,iedge) + IF(DELEDGE(IEDGE).EQ.0.)GO TO 2002 + AEDGE(IEDGE)=(ABLOG8(1,IEDGE)-2.*ABLOG8(2,IEDGE)+ABLOG8(3,IEDGE))/ + 1 DELEDGE(IEDGE) + BEDGE(IEDGE)=(-(HALFEDGE(IEDGE)+WLEDGE(IEDGE+1))*ABLOG8(1,IEDGE)+ + 1 (WLEDGE(IEDGE)+WLEDGE(IEDGE+1))*2.*ABLOG8(2,IEDGE)- + 2(WLEDGE(IEDGE)+HALFEDGE(IEDGE))*ABLOG8(3,IEDGE))/DELEDGE(IEDGE) + CEDGE(IEDGE)=(HALFEDGE(IEDGE)*WLEDGE(IEDGE+1)*ABLOG8(1,IEDGE)- + 1 WLEDGE(IEDGE)*WLEDGE(IEDGE+1)*2.*ABLOG8(2,IEDGE)+ + 2 WLEDGE(IEDGE)*HALFEDGE(IEDGE)*ABLOG8(3,IEDGE))/DELEDGE(IEDGE) + 2002 CONTINUE + NU=0 + 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)=CUTOFF*10.**((AEDGE(IEDGE)*WAVE+BEDGE(IEDGE))* + 1 WAVE+CEDGE(IEDGE)) +C + XNFPH(J,1)=XNFPEL(1) + XNFPH(J,2)=XNFPEL(2) + XNFPHE(J,1)=XNFPEL(7) + XNFPHE(J,2)=XNFPEL(8) + XNFPHE(J,3)=XNFPEL(9) + DO 205 NELION=1,594 +C PROBLEMS WITH OVERFLOW + XNFPEL(NELION)=XNFPEL(NELION)/RHO(J) +C QXNFDOP(NELION)=QXNFPEL(NELION)/QRHO(J)/QDOPPLE(NELION) + 205 XNFDOP(NELION)=XNFPEL(NELION)/DOPPLE(NELION) + TXNXN(J)=(XNFH(J)+.42*XNFHE(J,1)+.85*XNFH2(J))*(T(J)/10000.)**.3 +C +C ADD LINES TO BUFFER + MLINES=0 + IF(N19.GT.0)CALL XLINOP(J,N19) + IF(N12.EQ.0)GO TO 400 + L=10000 +C DO 350 ILINE=1,N12,10000 +C READ(12)NBUFF,CONGF,NELION,GAMRF,GAMSF,GAMWF + 301 L=L+1 + IF(L.EQ.10001)THEN + READ(12)TAPE12 + L=1 + ENDIF + NELION=ITAPE12(3,L) + IF(NELION.EQ.0)GO TO 400 + CONTMIN=CONTINUUM(ITAPE12(1,L)) + KAPPA0=TAPE12(2,L)*XNFDOP(NELION) +C KAPPA0=CONGF*XNFDOP(NELION) +C IF(KAPPA0.LT.CONTINUUM(NBUFF))GO TO 301 + IF(KAPPA0.LT.CONTMIN)GO TO 301 + MLINES=MLINES+1 +C ADAMP=(GAMRF+GAMSF*XNE(J)+GAMWF*TXNXN(J))/DOPPLE(NELION) + ADAMP=(TAPE12(4,L)+TAPE12(5,L)*XNE(J)+TAPE12(6,L)*TXNXN(J))/ + 1 DOPPLE(NELION) +C ADAMP=MIN(ADAMP,.2) + NBUFF=ITAPE12(1,L) +C +C BUFFER(NBUFF)=BUFFER(NBUFF)+KAPPA0*VOIGT(0.,ADAMP) +C DVOIGT=1./DOPPLE(NELION)/RESOLU +CC PROFILE INSIDE 10 DOPPLER WIDTHS +C N10DOP=10.*(DOPPLE(NELION)*RESOLU) +C DO 321 NSTEP=1,N10DOP +C PROFILE(NSTEP)=KAPPA0*VOIGT(FLOAT(NSTEP)*DVOIGT,ADAMP) +C IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 +C 321 CONTINUE +C +C COMPUTE VOIGT FUNCTION INLINE + IF(ADAMP.LE..2)THEN + BUFFER(NBUFF)=BUFFER(NBUFF)+KAPPA0*(1.-1.128*ADAMP) + TABSTEP=VSTEPS/(DOPPLE(NELION)*RESOLU) + TABI=1.5 +C PROFILE INSIDE 10 DOPPLER WIDTHS + N10DOP=10.*DOPPLE(NELION)*RESOLU + DO 1321 NSTEP=1,N10DOP + TABI=TABI+TABSTEP + PROFILE(NSTEP)=KAPPA0*(H0TAB(IFIX(TABI))+ADAMP*H1TAB(IFIX(TABI))) + IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 + 1321 CONTINUE +C FAR WINGS + X=PROFILE(N10DOP)*FLOAT(N10DOP)**2 + MAXSTEP=SQRT(X/CONTMIN)+1. + MAXSTEP=MIN(MAXSTEP,MAXPROF) +C CONGF HAS BEEN DIVIDED BY FREQ +C IF CONGF IS SMALL, THE WINGS ARE TRUNCATED AT 100 DOPPLER WIDTHS +C IF((CONGF*FREQ1000(NBUFF/1000+1).LT.0.0001) + IF(TAPE12(2,L)*FREQ1000(NBUFF/1000+1).LT.0.0001)MAXSTEP= + 1MIN(MAXSTEP,N10DOP*10) +C IF THE LINE IS MOLECULAR, THE WINGS ARE TRUNCATED AT 100 DOPPLER WIDTHS + IF(MOD(NELION,6).EQ.0)MAXSTEP=MIN(MAXSTEP,N10DOP*10) + N1=N10DOP+1 + DO 322 NSTEP=N1,MAXSTEP + 322 PROFILE(NSTEP)=X/FLOAT(NSTEP)**2 + NSTEP=MAXSTEP + GO TO 323 +C + ELSE IF(ADAMP.GT.100.)THEN + PROFILECENTER=KAPPA0/ADAMP*.79788/1.4142 + BUFFER(NBUFF)=BUFFER(NBUFF)+PROFILECENTER + AA=ADAMP*ADAMP + DVOIGT=1./DOPPLE(NELION)/RESOLU +C PROFILE INSIDE 10 DOPPLER WIDTHS +C N10DOP=10.*DOPPLE(NELION)*RESOLU + DO 1322 NSTEP=1,MAXPROF + PROFILE(NSTEP)=PROFILECENTER*AA/(AA+(NSTEP*DVOIGT)**2) + IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 + 1322 CONTINUE + NSTEP=MAXPROF + GO TO 323 +C + ELSE IF(ADAMP.GT.1.4)THEN + AA=ADAMP*ADAMP + PROFILEC=KAPPA0*ADAMP*.79788 + BUFFER(NBUFF)=BUFFER(NBUFF)+KAPPA0/ADAMP*.79788/1.4142* + 1 (1.+(.75/AA-.5)/AA) + DVOIGT=1./DOPPLE(NELION)/RESOLU + DO 1323 NSTEP=1,MAXPROF + VV=(NSTEP*DVOIGT)**2 + U=(AA+VV)*1.4142 + AAU=AA/U + VVU=VV/U + PROFILE(NSTEP)=PROFILEC/U* + 1 ((((AAU-10.*VVU)*AAU+5.*VVU*VVU)/U+VVU-AAU/3.)*3./U+1.) + IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 + 1323 CONTINUE + NSTEP=MAXPROF + GO TO 323 +C + ELSE +C ADAMP.GT.0.2.AND.ADAMP.LT.1.4 + IA=ADAMP*200.+1.5 + PROFILECENTER=KAPPA0*ATAB(IA,1) + BUFFER(NBUFF)=BUFFER(NBUFF)+PROFILECENTER + DVOIGT=1./DOPPLE(NELION)/RESOLU +C PROFILE INSIDE 5 DOPPLER WIDTHS + N5DOP=5.*DOPPLE(NELION)*RESOLU + DO 1324 NSTEP=1,N5DOP + IV=NSTEP*DVOIGT*200.+1.5 + PROFILE(NSTEP)=KAPPA0*ATAB(IA,IV) + IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 + 1324 CONTINUE + PROFILEC=KAPPA0*ADAMP*.79788 + 2324 DO 1325 NSTEP=N5DOP+1,MAXPROF + VV=(NSTEP*DVOIGT)**2 + U=(AA+VV)*1.4142 + AAU=AA/U + VVU=VV/U + PROFILE(NSTEP)=PROFILEC/U* + 1 ((((AAU-10.*VVU)*AAU+5.*VVU*VVU)/U+VVU-AAU/3.)*3/U+1.) + IF(PROFILE(NSTEP).LT.CONTMIN)GO TO 323 + 1325 CONTINUE + NSTEP=MAXPROF + GO TO 323 + ENDIF +C + 323 IF(NBUFF+NSTEP.LT.1.OR.NBUFF-NSTEP.GT.LENGTH)GO TO 301 + 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 301 +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) + GO TO 301 + 350 CONTINUE +C + 400 CONTINUE +C + DO 401 NBUFF=1,LENGTH + 401 BUFFER(NBUFF)=MAX(BUFFER(NBUFF),CONTINUUM(NBUFF)) + REWIND 14 + DO 402 K=1,LENBUFF,10000 + 402 WRITE(14)(BUFFER(I),I=K,K+9999) + IVT=0 + WRITE(6,404)TLOG10,PLOG10,IVT,MLINES + 404 FORMAT(' LOGT',F5.2,' LOGP',F5.2,' VT',I3, + 1' SPECTRUM COMPLETED WITH',I8,' LINES') +C CALL ENDTIME +C CALL BEGTIME +C CALL ROSSCALC(T(J),P(J),0.) +C CALL ENDTIME +C CALL BEGTIME + CALL DFCALC(ITEFF,J,IVT,T(J),P(J),0.) +C CALL ENDTIME +C + DO 430 IVT=1,20 + IF(IFVT(IVT).EQ.0)GO TO 430 +C CALL BEGTIME + REWIND 14 + DO 411 I=1,LENBUFF + 411 BUFFER(I)=0. + NBUFF=-10000 + LASTREC=LENGTH-(LENBUFF-10000) + DO 420 K=1,LENBUFF,10000 + READ(14)B10000 + NBUFF=NBUFF+10000 + DO 415 I=1,10000 + 415 BUFFER(NBUFF+I)=BUFFER(NBUFF+I)+B10000(I)*VTCENTER(IVT) + IF(K.EQ.1)THEN + DO 416 I=1,10000 + DO 416 IPROF=1,NVT(IVT) + 416 BUFFER(NBUFF+I+IPROF)=BUFFER(NBUFF+I+IPROF)+B10000(I)* + 1 VTPROF(IPROF,IVT) + DO 417 I=1,10000 + DO 417 IPROF=1,NVT(IVT) + NBACK=NBUFF+I-IPROF + 417 IF(NBACK.GT.0)BUFFER(NBACK)=BUFFER(NBACK)+B10000(I)* + 1 VTPROF(IPROF,IVT) + GO TO 420 + ENDIF + LASTREC=10000 + IF(K.GT.LENBUFF-10000)LASTREC=LENGTH-(LENBUFF-10000) + DO 418 I=1,LASTREC + DO 418 IPROF=1,NVT(IVT) + BUFFER(NBUFF+I+IPROF)=BUFFER(NBUFF+I+IPROF)+B10000(I)* + 1 VTPROF(IPROF,IVT) + 418 BUFFER(NBUFF+I-IPROF)=BUFFER(NBUFF+I-IPROF)+B10000(I)* + 1 VTPROF(IPROF,IVT) + 420 CONTINUE + VT=IVT + WRITE(6,405)TLOG10,PLOG10,IVT + 405 FORMAT(' LOGT',F5.2,' LOGP',F5.2,' VT',I3, + 1' SPECTRUM BROADENED') +C CALL ENDTIME +C CALL BEGTIME +C CALL ROSSCALC(T(J),P(J),VT) +C CALL ENDTIME +C CALL BEGTIME + CALL DFCALC(ITEFF,J,IVT,T(J),P(J),VT) +C WRITE(6,*)' DISTRIBUTION FUNCTIONS COMPLETED' +C CALL ENDTIME + 430 CONTINUE + 500 CONTINUE + 950 CONTINUE + CLOSE(UNIT=15) + CLOSE(UNIT=12,DISP='DELETE') + CLOSE(UNIT=14,DISP='DELETE') + CLOSE(UNIT=19,DISP='DELETE') + CALL EXIT + END + SUBROUTINE XLINOP(J,N19) + PARAMETER (kw=25) + PARAMETER (LENBUFF=3510000,MAXPROF=10000) + COMMON /BUFFER/BUFFER(LENBUFF),PROFILE(MAXPROF) + COMMON /CONTIN/CONTINUUM(LENBUFF) + 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(10001),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) + REAL*4 NSTARK,NDOPP,NMERGE,INGLIS + EQUIVALENCE (GAMMAS,ASHORE),(GAMMAW,BSHORE) + EQUIVALENCE (GF,G,CGF),(TYPE,NLAST),(GAMMAR,XSECT,GAUNT) + INTEGER TYPE + REAL*4 KAPPA,KAPMIN,KAPPA0,KAPCEN + REAL*4 KAPPA0RED,KAPPARED,KAPPA0BLUE,KAPPABLUE + REAL*4 OLDKAPPA,NEWKAPPA,DKAPPA +c REAL*8 WAVE,WCON,EMERGE,WMERGE,WSHIFT,CONTX(25,17),WTAIL + REAL*8 WAVE,WCON,EMERGE,WMERGE,WSHIFT,CONTX(26,17),WTAIL + REAL*8 RESOLU,RATIO,RATIOLG,SIGMA2,WLBEG,WLEND + REAL*8 WL,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.764,27419.659,12186.462,6854.871,4387.113, 1.00 + 2 3046.604,2238.320,1713.711,1354.044,1096.776, + 3 906.426,761.650,648.980,559.579,487.456/ + DATA CONTX/ + 1 109678.764,27419.659,12186.462,6854.871,4387.113, 1.00 + 2 3046.604,2238.320,1713.711,1354.044,1096.776,16*0., + 3 198310.760,38454.691,32033.214,29223.753,27175.760, 2.00 + 4 15073.868,0.,0.,0.,0.,16*0., + 5 438908.850,109726.529,48766.491,27430.925,17555.715, 2.01 + 6 12191.437,0.,0.,0.,0.,16*0., + 7 90883.840,90867.420,90840.420,90820.420,90804.000, 6.00 + 8 90777.000,80691.180,80627.760,69235.820,69172.400,16*0., + 9 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,16*0., 6.01 + A 61671.020,39820.615,39800.556,39759.842,22*0., 12.00 + 1 26*0., 12.01 + 2 48278.370,48166.309,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.000,65957.885,65811.843,65747.550,65670.435, 14.00 + 5 65524.393,59736.150,59448.700,50640.630,50553.180,16*0., + 6 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 +C + IF(ITEMP.EQ.ITEMP1)GO TO 95 + EHYD(1)=0. + EHYD(2)=82259.105 + EHYD(3)=97492.302 + EHYD(4)=102823.893 + EHYD(5)=105291.651 + EHYD(6)=106632.160 + EHYD(7)=107440.444 + EHYD(8)=107965.051 + 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)) +C WRITE(6,90) +C 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.312/NMERGE**2 + EMERGEH(K)=109677.576/NMERGE**2 +C 91 WRITE(6,92)K,NMERGE,EMERGE(K),EMERGEH(K) +C 92 FORMAT(I3,4F10.3) + 91 CONTINUE + ITEMP1=ITEMP +C + 95 REWIND 19 + DO 900 ILINE=1,N19 + READ(19)WL,GF,NBLO,NBUP,NELION,TYPE,NCON,NELIONX, + 1GAMMAR,GAMMAS,GAMMAW,NBUFF + 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) + 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 + KAPPA=XSECTG*XNFPEL(NELION) + 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 + 200 KAPPA0=CGF*XNFDOP(NELION) + IF(KAPPA0.LT.CONTINUUM(NBUFF))GO TO 900 + IFSTRONG=0 +C SINCE CONT=CONT/1000, CONT*500 IS 1/2 THE ACTUAL CONTINUUM. LINES ARE +C STRONG IF THE LINE CENTER IS 1/2 THE CONTINUUM. + IF(KAPPA0.GT.500.*CONTINUUM(NBUFF))IFSTRONG=1 + MLINES=MLINES+1 + WCON=0. + WTAIL=0. + IF(NCON.GT.0)THEN + WCON=1.E7/(CONTX(NCON,NELIONX)-EMERGE(J)) + WTAIL=1.E7/(CONTX(NCON,NELIONX)-EMERGE(J)-500.) + ENDIF + ADAMP=(GAMMAR+GAMMAS*XNE(J)+GAMMAW*TXNXN(J))/DOPPLE(NELION) + KAPCEN=KAPPA0*VOIGT(0.,ADAMP) + DOPWL=DOPPLE(NELION)*WL + IF(WL.GT.WLEND)GO TO 213 +C RED WING + MINRED=MAX(1,NBUFF) + MAXRED=MIN(NBUFF+1000,LENGTH) + IF(IFSTRONG.EQ.1)MAXRED=LENGTH + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 211 IBUFF=MINRED,MAXRED + 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) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 212 + 211 WAVE=WAVE*RATIO + 212 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WLBEG)GO TO 900 +C BLUE WING + 213 IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + IF(IFSTRONG.EQ.0)MAXBLUE=MIN(MAXBLUE,1000) + 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 +C KAPPA=KAPPA0*VOIGT(ABS(WAVE-WL)/DOPWL,ADAMP) + KAPPA=KAPPA0*VOIGT(VVOIGT,ADAMP) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))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 + 600 KAPPA0=CGF*XNFDOP(1) + IF(KAPPA0.LT.CONTINUUM(NBUFF))GO TO 900 + DOPPH(J)=DOPPLE(1) +C DEUTERIUM + IF(TYPE.EQ.-2)DOPPH(J)=DOPPH(J)/1.4142 + 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.) + WCON=MIN(WSHIFT+WSHIFT,WCON) + IF(WTAIL.LT.0.)WTAIL=WCON+WCON + WTAIL=MIN(WCON+WCON,WTAIL) + IF(WL.GT.WLEND)GO TO 613 +CC RED WING +C IF(WLBEG.GT.ALPHAHYD(NBLO))GO TO 613 +C REDCUT=1.D7/(109678.764D0-109677.576D0/(NBUP-0.8D0)**2-EHYD(NBLO)) +C WLMINUS1=1.D7/(EHYD(NBUP-1)-EHYD(NBLO)) +C WLMINUS2=1.D7/(EHYD(NBUP-2)-EHYD(NBLO)) +C KAPPA0RED=KAPPA0*HFNM(NBLO,NBUP-2)/HFNM(NBLO,NBUP)/ +C 1 (EHYD(NBUP-2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) +C MINRED=MAX0(1,NBUFF) +C WAVE=WBEGIN*RATIO**(MINRED-1) +C DO 611 IBUFF=MINRED,LENGTH +C IF(WAVE.LT.WCON)GO TO 611 +C IF(WAVE.GT.WLMINUS1)GO TO 612 +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) +C IF(WAVE.GT.REDCUT)THEN +C KAPPARED=KAPPA0RED*HPROF4(NBLO,NBUP-2,J,wl,WAVE-WLMINUS2,DOPPH) +C IF(WAVE.LT.WTAIL)KAPPARED=KAPPARED*(WAVE-WCON)/(WTAIL-WCON) +C IF(KAPPARED.GE.KAPPA)GO TO 612 +C ENDIF +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 612 +C 611 WAVE=WAVE*RATIO +C 612 IF(MINRED.EQ.1)GO TO 900 +C IF(WL.LT.WLBEG)GO TO 900 +CC BLUE WING +C 613 BLUECUT=1.D7/(109678.764D0-109677.576D0/(NBUP+0.8D0)**2- +C 1 EHYD(NBLO)) +C WLPLUS1=1.D7/(EHYD(NBUP+1)-EHYD(NBLO)) +C WLPLUS2=1.D7/(EHYD(NBUP+2)-EHYD(NBLO)) +C KAPPA0BLUE=KAPPA0*HFNM(NBLO,NBUP+2)/HFNM(NBLO,NBUP)/ +C 1 (EHYD(NBUP+2)-EHYD(NBLO))*(EHYD(NBUP)-EHYD(NBLO)) +C IBUFF=MIN0(LENGTH+1,NBUFF) +C MAXBLUE=IBUFF-1 +C WAVE=WBEGIN*RATIO**(IBUFF-1) +C DO 614 I=1,MAXBLUE +C IBUFF=IBUFF-1 +C WAVE=WAVE/RATIO +C IF(WAVE.LT.WCON)GO TO 900 +C IF(WAVE.LT.WLPLUS1)GO TO 900 +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) +C IF(WAVE.LT.BLUECUT)THEN +C KAPPABLUE=KAPPA0BLUE*HPROF4(NBLO,NBUP+2,J,wl,WAVE-WLPLUS2,DOPPH) +C IF(WAVE.LT.WTAIL)KAPPABLUE=KAPPABLUE*(WAVE-WCON)/(WTAIL-WCON) +C IF(KAPPABLUE.GT.KAPPA)GO TO 900 +C ENDIF +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 900 +C 614 CONTINUE +C GO TO 900 +CC ALPHA PRESUMED ISOLATED LINES AND BETA BLUE WINGS +C 620 IF(WL.GT.WLEND)GO TO 623 +CC RED WING +C MINRED=MAX0(1,NBUFF) +C WAVE=WBEGIN*RATIO**(MINRED-1) +C DO 621 IBUFF=MINRED,LENGTH +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 622 +C 621 WAVE=WAVE*RATIO +C 622 IF(MINRED.EQ.1)GO TO 900 +C IF(WL.LT.WLBEG)GO TO 900 +CC BLUE WING +C 623 IBUFF=MIN0(LENGTH+1,NBUFF) +C MAXBLUE=IBUFF-1 +C WAVE=WBEGIN*RATIO**(IBUFF-1) +C DO 624 I=1,MAXBLUE +C IBUFF=IBUFF-1 +C WAVE=WAVE/RATIO +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 900 +C 624 CONTINUE +C GO TO 900 +CC BETA LINES RED WING +C 630 IF(WL.GT.WLEND)GO TO 623 +CC RED WING +C MINRED=MAX0(1,NBUFF) +C WAVE=WBEGIN*RATIO**(MINRED-1) +C DO 631 IBUFF=MINRED,LENGTH +C IF(WAVE.GT.ALPHAHYD(NBLO))GO TO 633 +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 622 +C 631 WAVE=WAVE*RATIO +C GO TO 622 +C 633 IBUFFALPHA=IBUFF +C IF(NBUFF.LT.1)GO TO 900 +C IBUFF=NBUFF +C MAXBLUE=IBUFF-1 +C WAVE=WBEGIN*RATIO**(IBUFF-1) +C DO 634 I=1,MAXBLUE +C IBUFF=IBUFF-1 +C WAVE=WAVE/RATIO +C KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) +C BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA +C IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 900 +C 634 CONTINUE +C GO TO 900 +C HYDROGEN LINE +C 600 KAPPA0=CGF*XNFPEL(1)/DOPPLE(1) +C +C COMPUTES THE FIRST 20 POINTS, THEN EVERY 20TH POINT AND LOG INTERPOLATES +C THE 19 MISSING POINTS +C +C RED WING + IF(WLBEG.GT.ALPHAHYD(NBLO))GO TO 613 + REDCUT=1.D7/(109678.764D0-109677.576D0/(NBUP-0.8D0)**2-EHYD(NBLO)) + WLMINUS1=1.D7/(EHYD(NBUP-1)-EHYD(NBLO)) + WLMINUS2=1.D7/(EHYD(NBUP-2)-EHYD(NBLO)) + 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) + IB=0 + DO 611 IBUFF=MINRED,LENGTH + IF(WAVE.LT.WCON)GO TO 611 + IF(WAVE.GT.WLMINUS1)GO TO 612 + IB=IB+1 + IF(IB.GT.21)THEN + IF(MOD(IB,20).NE.1)GO TO 611 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + IF(WAVE.GT.REDCUT)THEN + KAPPARED=KAPPA0RED*HPROF4M(NBLO,NBUP-2,J,wl,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 + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(IB.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1611 K=1,19 + 1611 BUFFER(IBUFF-K)=BUFFER(IBUFF-K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 612 + 611 WAVE=WAVE*RATIO + 612 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WLBEG)GO TO 900 +C BLUE WING + 613 BLUECUT=1.D7/(109678.764D0-109677.576D0/(NBUP+0.8D0)**2- + 1 EHYD(NBLO)) + WLPLUS1=1.D7/(EHYD(NBUP+1)-EHYD(NBLO)) + WLPLUS2=1.D7/(EHYD(NBUP+2)-EHYD(NBLO)) + 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 + IF(I.GT.21)THEN + IF(MOD(I,20).NE.1)GO TO 614 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + IF(WAVE.LT.WTAIL)KAPPA=KAPPA*(WAVE-WCON)/(WTAIL-WCON) + IF(WAVE.LT.BLUECUT)THEN + KAPPABLUE=KAPPA0BLUE*HPROF4P(NBLO,NBUP+2,J,wl,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 + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(I.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1614 K=1,19 + 1614 BUFFER(IBUFF+K)=BUFFER(IBUFF+K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))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) + IB=0 + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 621 IBUFF=MINRED,LENGTH + IB=IB+1 + IF(IB.GT.21)THEN + IF(MOD(IB,20).NE.1)GO TO 621 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(IB.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1621 K=1,19 + 1621 BUFFER(IBUFF-K)=BUFFER(IBUFF-K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 622 + 621 WAVE=WAVE*RATIO + 622 IF(MINRED.EQ.1)GO TO 900 + IF(WL.LT.WLBEG)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 + IF(I.GT.21)THEN + IF(MOD(I,20).NE.1)GO TO 624 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(I.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1624 K=1,19 + 1624 BUFFER(IBUFF+K)=BUFFER(IBUFF+K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))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) + IB=0 + WAVE=WBEGIN*RATIO**(MINRED-1) + DO 631 IBUFF=MINRED,LENGTH + IF(WAVE.GT.ALPHAHYD(NBLO))GO TO 633 + IB=IB+1 + IF(IB.GT.21)THEN + IF(MOD(IB,20).NE.1)GO TO 631 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(IB.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1631 K=1,19 + 1631 BUFFER(IBUFF-K)=BUFFER(IBUFF-K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))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 + IF(I.GT.21)THEN + IF(MOD(I,20).NE.1)GO TO 634 + ENDIF + KAPPA=KAPPA0*HPROF4(NBLO,NBUP,J,wl,WAVE-WL,DOPPH) + BUFFER(IBUFF)=BUFFER(IBUFF)+KAPPA + NEWKAPPA=LOG(MAX(KAPPA,1.E-20)) + IF(I.GT.21)THEN + DKAPPA=(NEWKAPPA-OLDKAPPA)/20. + DO 1634 K=1,19 + 1634 BUFFER(IBUFF+K)=BUFFER(IBUFF+K)+EXP(NEWKAPPA-K*DKAPPA) + ENDIF + OLDKAPPA=NEWKAPPA + IF(KAPPA.LT.CONTINUUM(IBUFF))GO TO 900 + 634 CONTINUE + GO TO 900 +C +C +C AUTOIONIZING LINE + 700 KAPPA0=BSHORE*G*XNFPEL(NELION) + IF(KAPPA0.LT.CONTINUUM(NBUFF))GO TO 900 + MLINES=MLINES+1 + FRELIN=2.99792458E17/WL + IF(WL.GT.WLEND)GO TO 713 +C RED WING + MINRED=MAX0(1,NBUFF) + FREQ=2.99792458E17/(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))GO TO 712 + 711 FREQ=FREQ/RATIO + 712 IF(NBUFF.EQ.1)GO TO 900 +C BLUE WING + 713 IBUFF=MIN0(LENGTH+1,NBUFF) + MAXBLUE=IBUFF-1 + FREQ=2.99792458E17/(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))GO TO 900 + 714 CONTINUE + GO TO 900 +C + 900 CONTINUE + RETURN + END + FUNCTION HPROF4(N,M,J,wl0,DELW,DOPPH) +C FASTER AND STRIPPED FROM SYNTHE VERSION +C FROM DEANE PETERSON + PARAMETER (kw=25) + REAL*8 DELW,EMERGE,wl0 + 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) +C change from SYNTHE + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + 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(10001),E1TAB(2000) +C change from SYNTHE +c COMMON /EXTAB/EXTAB(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 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/ +C not used by dfsynthe +c FASTEX(X)=EXTAB(IFIX(X)+1)* +c 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 + 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.99792458E18/FREQNM**2/XKNM + WAVENM=2.99792458E18/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 +CC RESONT=HFNM(1,M)/XM/(1.-1./XM2) +CC IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2) + RESONT=HFNM(1,M)/(1.-1./XM2) + IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/(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 +CC RESONT=RESONT*3.92E-24/GNM +C 0.5773 is (ga/ge)**1/2 for H(n=1 (S))-H(n=2,3,(P)) +C =(2/6)**0.5 + RESONT=RESONT*3.579E-24*0.5773/GNM + VDW=4.45E-26/GNM*(XM2*(7.*XM2+5.))**.4 +CCC GUESS THAT H2 IS TWICE AS STRONG AS HE 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 +C 30 WL=WAVENM+DELW*10. + 30 WL=WL0*10.+DELW*10. + FREQ=2.99792458E18/WL + freq0=2.99792458E17/wl0 +c DEL=ABS(FREQ-FREQNM) + DEL=ABS(FREQ-FREQ0) +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 +c DOP=FREQNM*HWDOP + DOP=FREQ0*HWDOP + IF(IFCORE.EQ.1)GO TO (32,40,50),NWID +C DO DOPPLER +C PUT FINE STRUCTURE IN DOPPLER CORE + 32 DO 33 I=1,IFINS +c D=ABS(FREQ-FREQNM-FINEST(I))/DOP + D=ABS(FREQ-FREQ0-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 +C IF(D.LE.7.)HPROF4=HPROF4+FASTEX(D*D)*FINSWT(I) +C change in DFSYNTHE for speed + IF(D.LE.7.)HPROF4=HPROF4+H0TAB(INT(200.*D+1.5))*FINSWT(I) + 33 CONTINUE + IF(IFCORE.EQ.1)RETURN +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 +C if (lambda<1277A) + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)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.99792458E10)GO TO 43 +C TABULATED AT 200 CM-1 SPACING + SPACING=200.*2.99792458E10 + FREQ22000=(82259.105-22000.)*2.99792458E10 + IF(FREQ.LT.FREQ22000)THEN + CUTOFF=(CUTOFFH2(2)-CUTOFFH2(1))/SPACING*(FREQ-FREQ22000)+ + 1CUTOFFH2(1) + ELSE + ICUT=(FREQ-FREQ22000)/SPACING + ICUT=MIN(ICUT,89) + 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.99792458E10 + 43 HPROFRES=CUTOFF*1.77245*DOP +C + 44 HPROFRAD=0. +C RAYLEIGHSCATTERING EXCEPT NEAR DOPPLER CORE +C CROSSOVER FROM ABSORPTION TO RAYLEIGH SCATTERING IN HRAYOP +C LAMBDA>/=1217A +C +C IF(FREQ.LE.2.463E15)HPROFRAD=0. +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..755*3.28805E15) +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..77*3.28805E15) + IF(FREQ.GT.2.4190611E15.AND.FREQ.LT..77*3.28805E15) + 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 WAAL CUTOFF FOR HE AND FOR H2 +C GUESS BOTH 60000 CM-1=1.8E15 HZ + IF(FREQ.LT.1.8E15)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.288051E15.AND. + 1 FREQ.LT..890*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman gamma + IF(M.EQ.4.AND.FREQ.GT..936*3.288051E15.AND. + 1 FREQ.LT..938*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman delta + IF(M.EQ.5.AND.FREQ.GT..959*3.288051E15.AND. + 1 FREQ.LT..961*3.288051E15)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 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))* +C GAM=G1*(.5*FASTEX(AMIN1(80.,Y1))+VCSE1F(Y1)-.5*VCSE1F(Y2))* +C 1(1.-GCON1(J)/(1.+(90.*Y1)**3)-GCON2(J)/(1.+2000.*Y1)) +C change for DFSYNTHE for speed + GAM=G1*(.5*EXTAB(NINT(100*AMIN1(80.,Y1)))+ + 1VCSE1F(Y1)-.5*VCSE1F(Y2))* + 2(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.105-20000.)*2.99792458E10)GO TO 53 + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)GO TO 52 +C TABULATED AT 100 CM-1 SPACING + FREQ15000=(82259.105-15000.)*2.99792458E10 + SPACING=100.*2.99792458E10 + IF(FREQ.LT.FREQ15000)THEN + CUTOFF=(CUTOFFH2PLUS(2)-CUTOFFH2PLUS(1))/SPACING* + 1(FREQ-FREQ15000)+CUTOFFH2PLUS(1) + ELSE + ICUT=(FREQ-FREQ15000)/SPACING + icut=min(icut,109) + 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.99792458E10*XNFPH(J,2) + HPROF4=HPROF4+CUTOFF*1.77245*DOP + GO TO 53 + 52 BETA4000=4000.*2.99792458E10/FO(J)*DBETA + PRQSP4000=SOFBET(BETA4000,PP(J),N,M)*.5/FO(J)*DBETA + CUTOFF4000=(10.**(-11.07-14.))/2.99792458E10*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 HPROF4m(N,M,J,wl0,DELW,DOPPH) +C FASTER AND STRIPPED FROM SYNTHE VERSION +C FROM DEANE PETERSON + PARAMETER (kw=25) + REAL*8 DELW,EMERGE,wl0 + DIMENSION 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 /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + 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(10001),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 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/ +c FASTEX(X)=EXTAB(IFIX(X)+1)* +c 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.99792458E18/FREQNM**2/XKNM + WAVENM=2.99792458E18/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 +CC RESONT=HFNM(1,M)/XM/(1.-1./XM2) +CC IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2) + RESONT=HFNM(1,M)/(1.-1./XM2) + IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/(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 +CC RESONT=RESONT*3.92E-24/GNM +C 0.5773 is (ga/ge)**1/2 for H(n=1 (S))-H(n=2,3,(P)) +C =(2/6)**0.5 + RESONT=RESONT*3.579E-24*0.5773/GNM + VDW=4.45E-26/GNM*(XM2*(7.*XM2+5.))**.4 +C GUESS THAT H2 IS TWICE AS STRONG AS HE 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 +C 30 WL=WAVENM+DELW*10. + 30 WL=WL0*10.+DELW*10. + FREQ=2.99792458E18/WL + freq0=2.99792458E17/wl0 +c DEL=ABS(FREQ-FREQNM) + DEL=ABS(FREQ-FREQ0) +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. + HPROF4m=0. + IFCORE=0 + IF(ABS(DEL).LE.HFWID)IFCORE=1 +c DOP=FREQNM*HWDOP + DOP=FREQ0*HWDOP + IF(IFCORE.EQ.1)GO TO (32,40,50),NWID +C DO DOPPLER +C PUT FINE STRUCTURE IN DOPPLER CORE + 32 DO 33 I=1,IFINS +c D=ABS(FREQ-FREQNM-FINEST(I))/DOP + D=ABS(FREQ-FREQ0-FINEST(I))/DOP +C IF(D.LE.7.)HPROF4m=HPROF4m+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 +C IF(D.LE.7.)HPROF4m=HPROF4m+FASTEX(D*D)*FINSWT(I) +C change in DFSYNTHE for speed + IF(D.LE.7.)HPROF4m=HPROF4m+H0TAB(INT(200.*D+1.5))*FINSWT(I) + 33 CONTINUE + IF(IFCORE.EQ.1)RETURN +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 +C if (lambda<1277A) + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)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.99792458E10)GO TO 43 +C TABULATED AT 200 CM-1 SPACING + SPACING=200.*2.99792458E10 + FREQ22000=(82259.105-22000.)*2.99792458E10 + IF(FREQ.LT.FREQ22000)THEN + CUTOFF=(CUTOFFH2(2)-CUTOFFH2(1))/SPACING*(FREQ-FREQ22000)+ + 1CUTOFFH2(1) + ELSE + ICUT=(FREQ-FREQ22000)/SPACING + ICUT=MIN(ICUT,89) + 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.99792458E10 + 43 HPROFRES=CUTOFF*1.77245*DOP +C + 44 HPROFRAD=0. +C RAYLEIGHSCATTERING EXCEPT NEAR DOPPLER CORE +C CROSSOVER FROM ABSORPTION TO RAYLEIGH SCATTERING IN HRAYOP +C LAMBDA>/=1217A +C +C IF(FREQ.LE.2.463E15)HPROFRAD=0. +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..755*3.28805E15) +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..77*3.28805E15) + IF(FREQ.GT.2.4190611E15.AND.FREQ.LT..77*3.28805E15) + 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 WAAL CUTOFF FOR HE AND FOR H2 +C GUESS BOTH 60000 CM-1=1.8E15 HZ + IF(FREQ.LT.1.8E15)HPROFVDW=0. + HPROFLOR=HPROFRES+HPROFRAD+HPROFVDW + HPROF4m=HPROF4m+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.288051E15.AND. + 1 FREQ.LT..890*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman gamma + IF(M.EQ.4.AND.FREQ.GT..936*3.288051E15.AND. + 1 FREQ.LT..938*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman delta + IF(M.EQ.5.AND.FREQ.GT..959*3.288051E15.AND. + 1 FREQ.LT..961*3.288051E15)TOP=HHW-FREQNM*HWRAD + 49 HPROFLOR=TOP/3.14159/(DEL**2+HHW**2)*1.77245*DOP + HPROF4m=HPROF4m+HPROFLOR + IF(IFCORE.EQ.1)RETURN +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))* +C GAM=G1*(.5*FASTEX(AMIN1(80.,Y1))+VCSE1F(Y1)-.5*VCSE1F(Y2))* +C 1(1.-GCON1(J)/(1.+(90.*Y1)**3)-GCON2(J)/(1.+2000.*Y1)) +C change for DFSYNTHE for speed + GAM=G1*(.5*EXTAB(NINT(100.*AMIN1(80.,Y1)))+ + 1VCSE1F(Y1)-.5*VCSE1F(Y2))* + 2(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.105-20000.)*2.99792458E10)GO TO 53 + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)GO TO 52 +C TABULATED AT 100 CM-1 SPACING + FREQ15000=(82259.105-15000.)*2.99792458E10 + SPACING=100.*2.99792458E10 + IF(FREQ.LT.FREQ15000)THEN + CUTOFF=(CUTOFFH2PLUS(2)-CUTOFFH2PLUS(1))/SPACING* + 1(FREQ-FREQ15000)+CUTOFFH2PLUS(1) + ELSE + ICUT=(FREQ-FREQ15000)/SPACING + icut=min(icut,109) + 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.99792458E10*XNFPH(J,2) + HPROF4m=HPROF4m+CUTOFF*1.77245*DOP + GO TO 53 + 52 BETA4000=4000.*2.99792458E10/FO(J)*DBETA + PRQSP4000=SOFBET(BETA4000,PP(J),N,M)*.5/FO(J)*DBETA + CUTOFF4000=(10.**(-11.07-14.))/2.99792458E10*XNFPH(J,2) + HPROF4m=HPROF4m+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 + HPROF4m=HPROF4m+(PRQS*(1.+FNS)+F)/FO(J)*DBETA*1.77245*DOP + RETURN + END + FUNCTION HPROF4p(N,M,J,wl0,DELW,DOPPH) +C FASTER AND STRIPPED FROM SYNTHE VERSION +C FROM DEANE PETERSON + PARAMETER (kw=25) + REAL*8 DELW,EMERGE,wl0 + DIMENSION 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 /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + 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(10001),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 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/ +c FASTEX(X)=EXTAB(IFIX(X)+1)* +c 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.99792458E18/FREQNM**2/XKNM + WAVENM=2.99792458E18/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 +CC RESONT=HFNM(1,M)/XM/(1.-1./XM2) +CC IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/XN/(1.-1./XN2) + RESONT=HFNM(1,M)/(1.-1./XM2) + IF(N.NE.1)RESONT=RESONT+HFNM(1,N)/(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 +CC RESONT=RESONT*3.92E-24/GNM +C 0.5773 is (ga/ge)**1/2 for H(n=1 (S))-H(n=2,3,(P)) +C =(2/6)**0.5 + RESONT=RESONT*3.579E-24*0.5773/GNM + VDW=4.45E-26/GNM*(XM2*(7.*XM2+5.))**.4 +C GUESS THAT H2 IS TWICE AS STRONG AS HE 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 +C 30 WL=WAVENM+DELW*10. + 30 WL=WL0*10.+DELW*10. + FREQ=2.99792458E18/WL + freq0=2.99792458E17/wl0 +c DEL=ABS(FREQ-FREQNM) + DEL=ABS(FREQ-FREQ0) +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. + HPROF4p=0. + IFCORE=0 + IF(ABS(DEL).LE.HFWID)IFCORE=1 +c DOP=FREQNM*HWDOP + DOP=FREQ0*HWDOP + IF(IFCORE.EQ.1)GO TO (32,40,50),NWID +C DO DOPPLER +C PUT FINE STRUCTURE IN DOPPLER CORE + 32 DO 33 I=1,IFINS +c D=ABS(FREQ-FREQNM-FINEST(I))/DOP + D=ABS(FREQ-FREQ0-FINEST(I))/DOP +C IF(D.LE.7.)HPROF4p=HPROF4p+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 +C IF(D.LE.7.)HPROF4p=HPROF4p+FASTEX(D*D)*FINSWT(I) +C change in DFSYNTHE for speed + IF(D.LE.7.)HPROF4p=HPROF4p+H0TAB(INT(200.*D+1.5))*FINSWT(I) + 33 CONTINUE + IF(IFCORE.EQ.1)RETURN +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 +C if (lambda<1277A) + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)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.99792458E10)GO TO 43 +C TABULATED AT 200 CM-1 SPACING + SPACING=200.*2.99792458E10 + FREQ22000=(82259.105-22000.)*2.99792458E10 + IF(FREQ.LT.FREQ22000)THEN + CUTOFF=(CUTOFFH2(2)-CUTOFFH2(1))/SPACING*(FREQ-FREQ22000)+ + 1CUTOFFH2(1) + ELSE + ICUT=(FREQ-FREQ22000)/SPACING + ICUT=MIN(ICUT,89) + 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.99792458E10 + 43 HPROFRES=CUTOFF*1.77245*DOP +C + 44 HPROFRAD=0. +C RAYLEIGHSCATTERING EXCEPT NEAR DOPPLER CORE +C CROSSOVER FROM ABSORPTION TO RAYLEIGH SCATTERING IN HRAYOP +C LAMBDA>/=1217A +C +C IF(FREQ.LE.2.463E15)HPROFRAD=0. +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..755*3.28805E15) +c IF(FREQ.GT..74*3.288051E15.AND.FREQ.LT..77*3.28805E15) + IF(FREQ.GT.2.4190611E15.AND.FREQ.LT..77*3.28805E15) + 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 WAAL CUTOFF FOR HE AND FOR H2 +C GUESS BOTH 60000 CM-1=1.8E15 HZ + IF(FREQ.LT.1.8E15)HPROFVDW=0. + HPROFLOR=HPROFRES+HPROFRAD+HPROFVDW + HPROF4p=HPROF4p+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.288051E15.AND. + 1 FREQ.LT..890*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman gamma + IF(M.EQ.4.AND.FREQ.GT..936*3.288051E15.AND. + 1 FREQ.LT..938*3.288051E15)TOP=HHW-FREQNM*HWRAD +C Lyman delta + IF(M.EQ.5.AND.FREQ.GT..959*3.288051E15.AND. + 1 FREQ.LT..961*3.288051E15)TOP=HHW-FREQNM*HWRAD + 49 HPROFLOR=TOP/3.14159/(DEL**2+HHW**2)*1.77245*DOP + HPROF4p=HPROF4p+HPROFLOR + IF(IFCORE.EQ.1)RETURN +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))* +C GAM=G1*(.5*FASTEX(AMIN1(80.,Y1))+VCSE1F(Y1)-.5*VCSE1F(Y2))* +C 1(1.-GCON1(J)/(1.+(90.*Y1)**3)-GCON2(J)/(1.+2000.*Y1)) +C change for DFSYNTHE for speed + GAM=G1*(.5*EXTAB(NINT(100.*AMIN1(80.,Y1)))+ + 1VCSE1F(Y1)-.5*VCSE1F(Y2))* + 2(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.105-20000.)*2.99792458E10)GO TO 53 + IF(FREQ.GT.(82259.105-4000.)*2.99792458E10)GO TO 52 +C TABULATED AT 100 CM-1 SPACING + FREQ15000=(82259.105-15000.)*2.99792458E10 + SPACING=100.*2.99792458E10 + IF(FREQ.LT.FREQ15000)THEN + CUTOFF=(CUTOFFH2PLUS(2)-CUTOFFH2PLUS(1))/SPACING* + 1(FREQ-FREQ15000)+CUTOFFH2PLUS(1) + ELSE + ICUT=(FREQ-FREQ15000)/SPACING + icut=min(icut,109) + 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.99792458E10*XNFPH(J,2) + HPROF4p=HPROF4p+CUTOFF*1.77245*DOP + GO TO 53 + 52 BETA4000=4000.*2.99792458E10/FO(J)*DBETA + PRQSP4000=SOFBET(BETA4000,PP(J),N,M)*.5/FO(J)*DBETA + CUTOFF4000=(10.**(-11.07-14.))/2.99792458E10*XNFPH(J,2) + HPROF4p=HPROF4p+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 + HPROF4p=HPROF4p+(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(10001),E1TAB(2000) + 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)-.57721566+X*(.99999193+X*(-.24991055+X*(.05519968+ + 1X*(-.00976004+X*.00107857)))) + RETURN + 20 IF(X.GT.30.)RETURN + VCSE1F=(X*(X+2.334733)+.25062)/(X*(X+3.330657)+1.681534)/X*EXP(-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),PROBETA(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/ + DATA MSAVE,NSAVE,PSAVE,INDEXSAVE/0,0,0.,0./ + IF(B.LT.500.)GO TO 1 + SOFBET=(1.5/SQRT(B)+27./B**2)/B**2 + RETURN + 1 IF(M.EQ.MSAVE.AND.N.EQ.NSAVE.AND.P.EQ.PSAVE)GO TO 5 + MSAVE=M + NSAVE=N + PSAVE=P + INDEX=7 + IF(N.LE.3.AND.M-N.LE.2)INDEX=N+M-2 + IF(INDEX.EQ.INDEXSAVE)GO TO 5 + IM=MIN0(INT(5.*P)+1,4) + IP=IM+1 + WTPP=5.*(P-PP(IM)) + WTPM=1.-WTPP + CC=C(IP,INDEX)*WTPP+C(IM,INDEX)*WTPM + DD=D(IP,INDEX)*WTPP+D(IM,INDEX)*WTPM + DO 2 J=1,15 + 2 PROBETA(J)=PROPBM(IP,J,INDEX)*WTPP+PROPBM(IM,J,INDEX)*WTPM + INDEXSAVE=INDEX + 5 B2=B*B + SB=SQRT(B) + IF(B.LT.25.12)GO TO 7 + CORR=1.+DD/(CC+B*SB) + SOFBET=(1.5/SB+27./B2)/B2*CORR + RETURN + 7 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 + CORR=1.+PROBETA(JP)*WTBP+PROBETA(JM)*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 + END + FUNCTION SOFBETM(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),PROBETA(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/ + DATA MSAVE,NSAVE,PSAVE,INDEXSAVE/0,0,0.,0./ + IF(B.LT.500.)GO TO 1 + SOFBETM=(1.5/SQRT(B)+27./B**2)/B**2 + RETURN + 1 IF(M.EQ.MSAVE.AND.N.EQ.NSAVE.AND.P.EQ.PSAVE)GO TO 5 + MSAVE=M + NSAVE=N + PSAVE=P + INDEX=7 + IF(N.LE.3.AND.M-N.LE.2)INDEX=N+M-2 + IF(INDEX.EQ.INDEXSAVE)GO TO 5 + IM=MIN0(INT(5.*P)+1,4) + IP=IM+1 + WTPP=5.*(P-PP(IM)) + WTPM=1.-WTPP + CC=C(IP,INDEX)*WTPP+C(IM,INDEX)*WTPM + DD=D(IP,INDEX)*WTPP+D(IM,INDEX)*WTPM + DO 2 J=1,15 + 2 PROBETA(J)=PROPBM(IP,J,INDEX)*WTPP+PROPBM(IM,J,INDEX)*WTPM + INDEXSAVE=INDEX + 5 B2=B*B + SB=SQRT(B) + IF(B.LT.25.12)GO TO 7 + CORR=1.+DD/(CC+B*SB) + SOFBETM=(1.5/SB+27./B2)/B2*CORR + RETURN + 7 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 + CORR=1.+PROBETA(JP)*WTBP+PROBETA(JM)*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 + SOFBETM=(PR1*WT+PR2*(1.-WT))*CORR + RETURN + END + FUNCTION SOFBETP(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),PROBETA(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/ + DATA MSAVE,NSAVE,PSAVE,INDEXSAVE/0,0,0.,0./ + IF(B.LT.500.)GO TO 1 + SOFBETP=(1.5/SQRT(B)+27./B**2)/B**2 + RETURN + 1 IF(M.EQ.MSAVE.AND.N.EQ.NSAVE.AND.P.EQ.PSAVE)GO TO 5 + MSAVE=M + NSAVE=N + PSAVE=P + INDEX=7 + IF(N.LE.3.AND.M-N.LE.2)INDEX=N+M-2 + IF(INDEX.EQ.INDEXSAVE)GO TO 5 + IM=MIN0(INT(5.*P)+1,4) + IP=IM+1 + WTPP=5.*(P-PP(IM)) + WTPM=1.-WTPP + CC=C(IP,INDEX)*WTPP+C(IM,INDEX)*WTPM + DD=D(IP,INDEX)*WTPP+D(IM,INDEX)*WTPM + DO 2 J=1,15 + 2 PROBETA(J)=PROPBM(IP,J,INDEX)*WTPP+PROPBM(IM,J,INDEX)*WTPM + INDEXSAVE=INDEX + 5 B2=B*B + SB=SQRT(B) + IF(B.LT.25.12)GO TO 7 + CORR=1.+DD/(CC+B*SB) + SOFBETP=(1.5/SB+27./B2)/B2*CORR + RETURN + 7 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 + CORR=1.+PROBETA(JP)*WTBP+PROBETA(JM)*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 + SOFBETP=(PR1*WT+PR2*(1.-WT))*CORR + RETURN + END + SUBROUTINE HE2LIN + RETURN + END + SUBROUTINE TABVOIGT(VSTEPS,N) + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + 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) + DO 5 IA=41,281 + A=(IA-1)/200. + AA=A*A + DO 4 IV=1,1001 + V=(IV-1)/200. + VV=V*V + IF(A+V.GT.3.2)GO TO 3 + 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 + ATAB(IA,IV)=((((HH4*A+HH3)*A+HH2)*A+HH1)*A+H0TAB(IV))* + 1 (((-.122727278*A+.532770573)*A-.96284325)*A+.979895032) + GO TO 4 + 3 U=(AA+VV)*1.4142 + AAU=AA/U + VVU=VV/U + ATAB(IA,IV)=A*.79788/U* + 1((((AAU-10.*VVU)*AAU+5.*VVU*VVU)/U+VVU-AAU/3.)*3./U+1.) + 4 CONTINUE + 5 CONTINUE + RETURN + END + FUNCTION VOIGT(V,A) +C FAST VOIGT + COMMON /H1TAB/H0TAB(2001),H1TAB(2001),H2TAB(2001),ATAB(281,1001) + IF(A.LT..2)GO TO 10 + IF(A.GT.1.4)GO TO 2 + IF(A+V.GT.3.2)GO TO 2 + IV=V*200.+1.5 + IA=A*200.+1.5 + VOIGT=ATAB(IA,IV) +C VV=V*V +C HH1=H1TAB(IV)+H0TAB(IV)*1.12838 +C HH2=H2TAB(IV)+HH1*1.12838-H0TAB(IV) +C HH3=(1.-H2TAB(IV))*.37613-HH1*.66667*VV+HH2*1.12838 +C HH4=(3.*HH3-HH1)*.37613+H0TAB(IV)*.66667*VV*VV +C VOIGT=((((HH4*A+HH3)*A+HH2)*A+HH1)*A+H0TAB(IV))* +C 1 (((-.122727278*A+.532770573)*A-.96284325)*A+.979895032) + 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 + VOIGT=VOIGT* + 1((((AAU-10.*VVU)*AAU+5.*VVU*VVU)/U+VVU-AAU/3.)*3./U+1.) + RETURN + 10 IF(V.GT.10.)GO TO 12 + IV=V*200.+1.5 + 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.43668255,4.42054938,3.16274620,7.68641124, + 1 5.65655216/ + DATA C0,C1,C2,D1,D2/.0012102205,.98147989,.75339742,1.6198645, + 1 .29135151/ + DATA E0,E1,F1/-.9969698,-.4257849,2.318261/ + 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(10001),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 +C SUBROUTINE ABORT +C PARAMETER (SS$_ABORT='002C'X) +C CALL LIB$STOP(%VAL(SS$_ABORT)) +C END +C SUBROUTINE BEGTIME +C WRITES OUT ELAPSED SYSTEM ACCOUNTING MESSAGE +C PARAMETER (SS$_ABORT='002C'X) +C IF(.NOT.LIB$INIT_TIMER())CALL LIB$STOP(%VAL(SS$_ABORT)) +C RETURN +C ENTRY ENDTIME +C IF(.NOT.LIB$SHOW_TIMER())CALL LIB$STOP(%VAL(SS$_ABORT)) +C RETURN +C END + SUBROUTINE SELECTLINES(N12,N19) + PARAMETER (LENBUFF=3510000,MAXPROF=10000) + IMPLICIT REAL*4 (A-H,O-Z) + COMMON /CONTIN/CONTINUUM(LENBUFF) + REAL*8 RESOLU,RATIO,RATIOLG,WLBEG,WLEND,RATIOLOG + REAL*8 WL,E,EP,WLVAC,START,STOP + REAL*8 TEFF,GLOG + INTEGER*2 IELION,IELO,IGFLOG,IGR,IGS,IGW +C COMMON /IIIIIII/IWL,IELION,IELO,IGFLOG,IGR,IGS,IGW +C INTEGER*4 IIIIIII(4) +C EQUIVALENCE (IIIIIII(1),IWL) + COMMON /TAPE1112/ITAPE12(6,10000) + REAL*4 TAPE12(6,10000) + INTEGER*4 ITAPE11(4,2000) + INTEGER*2 I2TAPE11(8,2000) + EQUIVALENCE (TAPE12(1,1),ITAPE12(1,1)) + EQUIVALENCE (ITAPE11(1,1),I2TAPE11(1,1),CONTINUUM(1)) + INTEGER*4 GNNNGGG(11) + EQUIVALENCE (GNNNGGG(1),GF) + REAL*8 FREQEDGE(377),WLEDGE(377),CMEDGE(377),FREQSET(1131) + REAL*8 WAVETAB(377) + INTEGER*4 IWLTAB(377) + REAL*4 TABLOG(32768),TABBOLT(32768) + REAL*4 XNFDOPMAX(377,594) + EQUIVALENCE (XNFDOPMAX(1,1),CONTINUUM(8001)) +C + DO 1 I=1,32768 + 1 TABLOG(I)=10.**((I-16384)*.001) + RATIO=1.D0+1.D0/500000.D0 + RATIOLG=LOG(RATIO) + RATIOLG4=RATIOLG + WLBEG=8.97666 + WLEND=10000. + WLBEG0=WLBEG/RATIO + WLBEG04=WLBEG0 + IXWLBEG=DLOG(WLBEG)/RATIOLG + IF(DEXP(IXWLBEG*RATIOLG).LT.WLBEG)IXWLBEG=IXWLBEG+1 +C + READ(22)NT,TEFF,GLOG + T=TEFF + HKT=6.6256E-27/1.38054E-16/T + HCKT=HKT*2.99792458E10 + HKT100=HKT*100. + HCKT100=HCKT*100. + READ(22)NEDGE,(FREQEDGE(I),WLEDGE(I),CMEDGE(I),I=1,NEDGE) + READ(22)NUMNU,(FREQSET(NU),NU=1,NUMNU) + READ(22)XNFDOPMAX +C PRINT 3333,(XNFDOPMAX(200,NELION),NELION=1,594) +C PRINT 3333,(XNFDOPMAX(250,NELION),NELION=1,594) +C PRINT 3333,(XNFDOPMAX(300,NELION),NELION=1,594) + 3333 FORMAT(1P6E12.3) + DO 2 I=1,NEDGE + WAVETAB(I)=ABS(WLEDGE(I)) + IXWL=DLOG(WAVETAB(I))/RATIOLG + IXWL=IXWL+1 +C THIS IWL IS ACTUALLY NBUFF + 2 IWLTAB(I)=IXWL-IXWLBEG+1 + DO 3 I=1,32768 + 3 TABBOLT(I)=EXP(-(I-1)*HCKT) +C + N12=0 + N22=0 + N32=0 + N42=0 + N44=0 + N52=0 + N62=0 + NU=1 + NBUFF=0 +C +C LOWLINES + OPEN(UNIT=11,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=669) + REWIND 12 + L=0 + L11=2000 + DO 580 LINE=1,31108567 + L11=L11+1 + IF(L11.EQ.2001)THEN + READ(11)ITAPE11 + L11=1 + ENDIF +C READ(11)IIIIIII +C READ(11)IWL,IELION,IELO,IGFLOG,IGR,IGS,IGW + IWL=ITAPE11(1,L11) + 570 IF(IWL.LT.IWLTAB(NU))GO TO 571 + NU=NU+1 + WRITE(6,*) NU,IWLTAB(NU) + GO TO 570 +C LINECEN=.026538/1.77245*GF*EXP(-ELO*HCKT)*STIM*XNFP/RHO/DOPPLE/FREQ +C 571 NELION=IELION + 571 NELION=I2TAPE11(3,L11) + IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 580 +C CONGF=.026538/1.77245*TABLOG(IGFLOG) + CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 580 +C ELO=TABLOG(IELO) + ELO=TABLOG(I2TAPE11(4,L11)) + CONGF=CONGF*EXP(-ELO*HCKT) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 580 + IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 + NBUFF=IWL + L=L+1 + ITAPE12(1,L)=NBUFF + FREQ4=2.99792458E17/WLVAC4 +C CONGF=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + ITAPE12(3,L)=NELION + FRQ4PI=FREQ4*12.5664 +C GAMRF=TABLOG(IGR)/FRQ4PI +C GAMSF=TABLOG(IGS)/FRQ4PI +C GAMWF=TABLOG(IGW)/FRQ4PI + TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI + TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI + TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI +C WRITE(12)NBUFF,CONGF,NELION,GAMRF,GAMSF,GAMWF + N12=N12+1 + IF(L.EQ.10000)THEN + WRITE(12)TAPE12 + L=0 + ENDIF + 580 CONTINUE + WRITE(6,8881)N12 + 8881 FORMAT(I10,' LINES FROM LOWLINES') + CLOSE(UNIT=11) +C +C HILINES + 669 IF(T.LT.20000.)GO TO 769 + OPEN(UNIT=21,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=769) + NU=1 + NBUFF=0 + L11=2000 + DO 680 LINE=1,10057574 + L11=L11+1 + IF(L11.EQ.2001)THEN + READ(21)ITAPE11 + L11=1 + ENDIF + IWL=ITAPE11(1,L11) + if(iwl.eq.0)go to 680 + 670 IF(IWL.LT.IWLTAB(NU))GO TO 671 + NU=NU+1 + GO TO 670 + 671 NELION=I2TAPE11(3,L11) + IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 680 + CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 680 + ELO=TABLOG(I2TAPE11(4,L11)) + CONGF=CONGF*EXP(-ELO*HCKT) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 680 + IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 + NBUFF=IWL + L=L+1 + ITAPE12(1,L)=NBUFF + FREQ4=2.99792458E17/WLVAC4 + TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + ITAPE12(3,L)=NELION + FRQ4PI=FREQ4*12.5664 + TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI + TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI + TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI + N22=N22+1 + IF(L.EQ.10000)THEN + WRITE(12)TAPE12 + L=0 + ENDIF + 680 CONTINUE + WRITE(6,8882)N22 + 8882 FORMAT(I10,' LINES FROM HILINES') + CLOSE(UNIT=21) +C +C DIATOMICS + 769 IF(T.GE.10000.)GO TO 869 + OPEN(UNIT=31,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=869) + NU=1 + NBUFF=0 + L11=2000 + DO 780 LINE=1,7771848 + L11=L11+1 + IF(L11.EQ.2001)THEN + READ(31)ITAPE11 + L11=1 + ENDIF + IWL=ITAPE11(1,L11) + 770 IF(IWL.LT.IWLTAB(NU))GO TO 771 + NU=NU+1 + GO TO 770 + 771 NELION=I2TAPE11(3,L11) + IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 780 + CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 780 + ELO=TABLOG(I2TAPE11(4,L11)) + CONGF=CONGF*EXP(-ELO*HCKT) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 780 + IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 + NBUFF=IWL + L=L+1 + ITAPE12(1,L)=NBUFF + FREQ4=2.99792458E17/WLVAC4 + TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + ITAPE12(3,L)=NELION + FRQ4PI=FREQ4*12.5664 + TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI + TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI + TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI + N32=N32+1 + IF(L.EQ.10000)THEN + WRITE(12)TAPE12 + L=0 + ENDIF + 780 CONTINUE + WRITE(6,8883)N32 + 8883 FORMAT(I10,' LINES FROM DIATOMICS') + CLOSE(UNIT=31) +C +C TIO + 869 IF(T.GE.5000.)GO TO 969 + OPEN(UNIT=41,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=969) + NU=1 + NBUFF=0 + L11=2000 + DO 880 LINE=1,36979284 + L11=L11+1 + IF(L11.EQ.2001)THEN + READ(41)ITAPE11 + L11=1 + ENDIF + IWL=ITAPE11(1,L11) + 870 IF(IWL.LT.IWLTAB(NU))GO TO 871 + NU=NU+1 + GO TO 870 + 871 NELION=I2TAPE11(3,L11) + IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 880 + CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 880 + ELO=TABLOG(I2TAPE11(4,L11)) + CONGF=CONGF*EXP(-ELO*HCKT) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 880 + IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 + NBUFF=IWL + L=L+1 + ITAPE12(1,L)=NBUFF + FREQ4=2.99792458E17/WLVAC4 + TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + ITAPE12(3,L)=NELION + FRQ4PI=FREQ4*12.5664 + TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI + TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI + TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI + N42=N42+1 + IF(L.EQ.10000)THEN + WRITE(12)TAPE12 + L=0 + ENDIF + 880 CONTINUE + WRITE(6,8884)N42 + 8884 FORMAT(I10,' LINES FROM TIOLINES') + CLOSE(UNIT=41) +C +C +C H2O + 969 IF(T.GE.5000.)GO TO 979 + OPEN(UNIT=43,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=979) + NU=1 + NBUFF=0 + L11=2000 +C DO 1880 LINE=1,48999850 +C LAST LINE IS OUTSIDE RANGE + DO 1880 LINE=1,48999843 + L11=L11+1 + IF(L11.EQ.2001)THEN + READ(43)ITAPE11 + L11=1 + ENDIF + IWL=ITAPE11(1,L11) + 1870 IF(IWL.LT.IWLTAB(NU))GO TO 1871 + NU=NU+1 + GO TO 1870 + 1871 NELION=I2TAPE11(3,L11) + IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 1880 + CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 1880 +C ELO=TABLOG(I2TAPE11(4,L11)) +C CONGF=CONGF*EXP(-ELO*HCKT) + CONGF=CONGF*TABBOLT(I2TAPE11(4,L11)+1) + IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 1880 + IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 + NBUFF=IWL + L=L+1 + ITAPE12(1,L)=NBUFF + FREQ4=2.99792458E17/WLVAC4 + TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) + ITAPE12(3,L)=NELION + FRQ4PI=FREQ4*12.5664 + TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI + TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI + TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI + N44=N44+1 + IF(L.EQ.10000)THEN + WRITE(12)TAPE12 + L=0 + ENDIF + 1880 CONTINUE + WRITE(6,1884)N44 + 1884 FORMAT(I10,' LINES FROM H2OLINES') + CLOSE(UNIT=43) +C +C NLTELINESA + 979 OPEN(UNIT=51,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, + 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=16,ERR=369) + OPEN(UNIT=19,TYPE='NEW',FORM='UNFORMATTED', + 1RECORDTYPE='FIXED',BLOCKSIZE=28000,RECL=14) + NU=1 + DO 980 LINE=1,999999 +C READ(51)WL,ELO,GF,NBLO,NBUP,NELION,TYPE,NCON,NELIONX, +C 1GAMMAR,GAMMAS,GAMMAW,NBUFF + READ(51,END=981)WL,ELO,GNNNGGG + NELION=GNNNGGG(4) + NBUFF=GNNNGGG(11) + 970 IF(NBUFF.LT.IWLTAB(NU))GO TO 971 + NU=NU+1 + GO TO 970 + 971 IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 980 + BOLTSTIM=EXP(-ELO*HCKT)*(1.-EXP(-1.E7/WL*HCKT)) + IF(BOLTSTIM*XNFDOPMAX(NU,NELION).LT.1.)GO TO 980 + GF=GF*BOLTSTIM +C WRITE(19)WL,GF,NBLO,NBUP,NELION,TYPE,NCON,NELIONX, +C 1GAMMAR,GAMMAS,GAMMAW,NBUFF + WRITE(19)WL,GNNNGGG + N52=N52+1 + 980 CONTINUE + 981 WRITE(6,8885)N52 + 8885 FORMAT(I10,' LINES FROM NLTELINESA') + N19=N52 + CLOSE(UNIT=51) + + 369 CONTINUE +CC +CC NLTELINESB +C 369 OPEN(UNIT=61,TYPE='OLD',FORM='UNFORMATTED',SHARED,READONLY, +C 1RECORDTYPE='FIXED',BLOCKSIZE=32000,RECL=8000,ERR=882) +C L=0 +C L11=2000 +C DO 380 LINE=1,29226 +C L11=L11+1 +C IF(L11.EQ.2001)THEN +C READ(61)ITAPE11 +C L11=1 +C ENDIF +C IWL=ITAPE11(1,L11) +C 370 IF(IWL.LT.IWLTAB(NU))GO TO 371 +C NU=NU+1 +C GO TO 370 +C 371 NELION=I2TAPE11(3,L11) +C IF(XNFDOPMAX(NU,NELION).LT.1.)GO TO 380 +C CONGF=.026538/1.77245*TABLOG(I2TAPE11(5,L11)) +C IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 380 +C ELO=TABLOG(I2TAPE11(4,L11)) +C CONGF=CONGF*EXP(-ELO*HCKT) +C IF(CONGF*XNFDOPMAX(NU,NELION).LT.1.)GO TO 380 +C IF(IWL.NE.NBUFF)WLVAC4=EXP(IWL*RATIOLG4)*WLBEG04 +C NBUFF=IWL +C L=L+1 +C ITAPE12(1,L)=NBUFF +C FREQ4=2.99792458E17/WLVAC4 +C TAPE12(2,L)=CONGF/FREQ4*(1.-EXP(-FREQ4*HKT)) +C ITAPE12(3,L)=NELION +C FRQ4PI=FREQ4*12.5664 +C TAPE12(4,L)=TABLOG(I2TAPE11(6,L11))/FRQ4PI +C TAPE12(5,L)=TABLOG(I2TAPE11(7,L11))/FRQ4PI +C TAPE12(6,L)=TABLOG(I2TAPE11(8,L11))/FRQ4PI +C N62=N62+1 +C IF(L.EQ.10000)THEN +C WRITE(12)TAPE12 +C L=0 +C ENDIF +C 380 CONTINUE +C 381 WRITE(6,8886)N62 +C 8886 FORMAT(I10,' LINES FROM NLTELINESB') +C CLOSE(UNIT=61) + + 882 N12=N12+N22+N32+N42+N44+N62 + NLINES=N12+N19 + WRITE(6,8888)NLINES + 8888 FORMAT(I10,' LINES TOTAL') + L=L+1 + ITAPE12(3,L)=0 + WRITE(12)TAPE12 + RETURN + END + SUBROUTINE DFCALC(ITEFF,J,IVT,T,P,VT) + PARAMETER (LENBUFF=3510000,MAXPROF=10000) + IMPLICIT REAL*4 (A-H,O-Z) + COMMON /NBIGLIT/NBEG(1540),NEND(1540),NSTEPS(13,1540),NN(1540) + COMMON /BUFFER/BUFFER(LENBUFF),PROFILE(MAXPROF),B10000(10000) + INTEGER*4 IBUFFER(LENBUFF) + EQUIVALENCE (BUFFER(1),IBUFFER(1)) + COMMON /XNFDOPMAX/XNFDOPMAX(377,594) + INTEGER*4 IBUFFSORT(223938) + EQUIVALENCE (XNFDOPMAX(1,1),IBUFFSORT(1)) + DIMENSION STEPS(12,1540),STEPSBIG(12,328),STEPSLIT(12,1212) + EQUIVALENCE (STEPS(1,1),STEPSBIG(1,1)) + EQUIVALENCE (STEPS(1,329),STEPSLIT(1,1)) + INTEGER*4 IDFOUT(7) + INTEGER*2 ISTEPS(12) + EQUIVALENCE (IDFOUT(2),ISTEPS(1)) + TLOG=ALOG10(T) + PLOG=ALOG10(P) + LENGTH=3507859 + DO 1 NBUFF=1,LENGTH + 1 IBUFFER(NBUFF)=NINT(LOG10(BUFFER(NBUFF))*1000.) + DO 5 IDF=1,1540 + NBEGIN=NBEG(IDF) + NSTOP=NEND(IDF) + NTOT=NSTOP-NBEGIN+1 +C PRINT 88,IDF,NBEGIN,NSTOP,NTOT + 88 FORMAT(10I10) + DO 6 I=NBEGIN,NSTOP + 6 IBUFFSORT(I-NBEGIN+1)=IBUFFER(I) + CALL BISORT(IBUFFSORT,NTOT,IBUFFSORT(100000)) + DO 3 ISTEP=1,12 + ISUM=0 + DO 2 I=NSTEPS(ISTEP,IDF),NSTEPS(ISTEP+1,IDF)-1 + 2 ISUM=ISUM+IBUFFSORT(I) + 3 ISTEPS(ISTEP)=(ISUM+(NSTEPS(ISTEP+1,IDF)-NSTEPS(ISTEP,IDF))/2)/ + 1 (NSTEPS(ISTEP+1,IDF)-NSTEPS(ISTEP,IDF)) + IDFOUT(1)=IVT*100000000+IDF*10000+ITEFF*100+J +C WRITE(15)IDFOUT + IF(IVT.EQ.0)WRITE(15)IDFOUT + IF(IVT.EQ.1)WRITE(16)IDFOUT + IF(IVT.EQ.2)WRITE(17)IDFOUT + IF(IVT.EQ.4)WRITE(18)IDFOUT + IF(IVT.EQ.8)WRITE(20)IDFOUT +C WRITE(8,77)IDFOUT(1),isteps +C PRINT 77,IDFOUT(1),ISTEPS + 77 FORMAT(13I10) + 5 CONTINUE + RETURN + END + SUBROUTINE BISORT(K,N,K1) + DIMENSION K(1),K1(1) + DO 111 J=1,N + 111 K(J)=K(J)+30000 + MASK=1 + MAX=K(1) + DO 1 J=2,N + 1 MAX=MAX0(MAX,K(J)) + DO 10 I=1,59 + IF(I.EQ.1)GO TO 6 + MAX=MAX/2 + MASK=MASK*2 + 6 IF(MAX.EQ.0)THEN + DO 112 J=1,N + 112 K(J)=K(J)-30000 + RETURN + ENDIF + J1=0 + J0=0 + DO 4 J=1,N + IF(K(J).AND.MASK)3,3,2 + 2 J1=J1+1 + K1(J1)=K(J) + GO TO 4 + 3 J0=J0+1 + K(J0)=K(J) + 4 CONTINUE + IF(J1.EQ.0)GO TO 10 + DO 5 J=1,J1 + JJ0=J+J0 + 5 K(JJ0)=K1(J) + 10 CONTINUE + END + SUBROUTINE DFINTERVALS + IMPLICIT REAL*8 (A-H,O-Z) + COMMON /NBIGLIT/NBIGBEG(328),NLITBEG(1212),NBIGEND(328), + 1 NLITEND(1212),NSTEPS(13,1540),NN(1540) + COMMON /BIGLIT/WAVEBIG(329),WAVELIT(1213) + REAL*8 BIGA(95),BIGB(125),BIGC(109) + REAL*8 LITA(95),LITB(190),LITC(190),LITD(190),LITE(190) + REAL*8 LITF(190),LITG(168) + EQUIVALENCE (WAVEBIG(1),BIGA(1)),(WAVEBIG(96),BIGB(1)) + EQUIVALENCE (WAVEBIG(221),BIGC(1)) + EQUIVALENCE (WAVELIT(1),LITA(1)),(WAVELIT(96),LITB(1)) + EQUIVALENCE (WAVELIT(286),LITC(1)),(WAVELIT(476),LITD(1)) + EQUIVALENCE (WAVELIT(666),LITE(1)),(WAVELIT(856),LITF(1)) + EQUIVALENCE (WAVELIT(1046),LITG(1)) + DATA BIGA/ + 1 8.97666, 9.2, 9.5, 9.71730, 9.81590, + 2 10.10901, 10.3, 10.46451, 10.65, 10.88545, + 3 11.2, 11.6, 11.95562, 12.3, 12.66565, + 4 12.74676, 12.93504, 13.16033, 13.31723, 13.46700, + 5 13.86253, 14.10462, 14.54782, 14.9, 15.3, + 6 15.74874, 16.00346, 16.4, 16.8, 17.25088, + 7 17.43258, 17.93660, 18.10114, 18.24438, 18.97796, + 8 19.22422, 19.55098, 20.11987, 20.23592, 20.8, + 9 21.3, 21.94783, 22.01983, 22.58070, 22.78377, + A 23.2, 23.65221, 24.3, 25.00744, 25.46614, + 1 25.89239, 26.14316, 26.6, 27.1, 27.6, + 2 28.09263, 28.65871, 29.3, 29.95849, 30.26260, + 3 31.3, 32.3, 33.3, 34.3, 35.27931, + 4 36.14100, 37., 38.01600, 38.96642, 40., + 5 41., 41.87992, 42.5, 43.5, 44.73260, + 6 45.5, 46.5, 47.5, 48.50178, 49.5, + 7 50.42590, 50.87630, 52., 54., 56., + 8 57.40614, 59.5, 61.5, 63.5, 65.53786, + 9 67., 69., 71., 72.24011, 74./ + DATA BIGB/ + 1 76., 78., 80., 82., 84., + 2 86., 88., 90., 91.17535, 94., + 3 98., 102., 106., 110.03056, 113., + 4 116., 120., 123.92928, 128., 132., + 5 136., 140., 144.43391, 148., 151.43485, + 6 157., 162.15072, 167.40282, 171., 175., + 7 180., 184., 188., 193., 197.46990, + 8 202., 207.13210, 210., 215., 220., + 9 225., 230., 235., 240., 245., + A 251.12621, 255., 260., 265., 270., + 1 275., 280., 285., 290., 300., + 2 310., 320., 330., 340., 350., + 3 360., 364.70183, 370., 380., 390., + 4 400., 410., 420., 430., 450., 460., 470., 480., 490., 500., + 5 510., 520., 530., 540., 550., 560., 570., 580., 590., 600., + 6 610., 620., 630., 640., 650., 660., 670., 680., 690., 700., + 7 710., 720., 730., 740., 750., 760., 770., 780., 790., 800., + 8 810.,820.58271,830.,840., 850., 860., 870., 880., 890., 900., + 9 910., 920., 930., 940., 950., 960., 970., 980., 990.,1000./ + DATA BIGC/ + 1 1025.,1050.,1075.,1100.,1125.,1150.,1175.,1200.,1225.,1250., + 2 1275.,1300.,1325.,1350.,1375.,1400.,1425.,1458.81670,1475.,1500., + 3 1525.,1550.,1575.,1600.,1640.,1680.,1720.,1760.,1800.,1840., + 4 1880.,1920.,1960.,2000.,2050.,2100.,2150.,2200.,2250.,2279.40330, + 5 2300.,2350.,2400.,2450.,2500.,2550.,2600.,2650.,2700.,2750., + 6 2800.,2850.,2900.,2950.,3000.,3050.,3100.,3150.,3200.,3282.34320, + 7 3400.,3500.,3600.,3700.,3800.,3900.,4000.,4100.,4200.,4300., + 8 4400.,4500.,4600.,4700.,4800.,4900.,5000.,5100.,5200.,5300., + 9 5400.,5500.,5600.,5700.,5800.,5900.,6000.,6100.,6200.,6300., + A 6400.,6600.,6800.,7000.,7200.,7400.,7600.,7800.,8000.,8200., + 1 8400.,8600.,8800.,9000.,9200.,9400.,9600.,9800.,10000./ + + DATA LITA/ + 1 8.97666, 9.2, 9.5, 9.71730, 9.81590, + 2 10.10901, 10.3, 10.46451, 10.65000, 10.88545, + 3 11.2, 11.6, 11.95562, 12.3, 12.66565, + 4 12.74676, 12.93504, 13.16033, 13.31723, 13.46700, + 5 13.86253, 14.10462, 14.54782, 14.9, 15.3, + 6 15.74874, 16.00346, 16.4, 16.8, 17.25088, + 7 17.43258, 17.93660, 18.10114, 18.24438, 18.97796, + 8 19.22422, 19.55098, 20.11987, 20.23592, 20.8, + 9 21.3, 21.94783, 22.01983, 22.58070, 22.78377, + A 23.2, 23.65221, 24.3, 25.00744, 25.46614, + 1 25.89239, 26.14316, 26.6, 27.1, 27.6, + 2 28.09263, 28.65871, 29.3, 29.95849, 30.26260, + 3 31.3, 32.3, 33.3, 34.3, 35.27931, + 4 36.14100, 37., 38.01600, 38.96642, 40., + 5 41., 41.87992, 42.5, 43.5, 44.73260, + 6 45.5, 46.5, 47.5, 48.50178, 49.5, + 7 50.42590, 50.87630, 51.5, 52.5, 53.5, + 8 54.5, 55.5, 56.5, 57.40614, 58.5, + 9 59.5, 60.5, 61.5, 62.5, 63.5/ + DATA LITB/ + 1 64.5, 65.53786, 66., 67., 68., 69., 70., 71., 72.24011, 73., + 2 74., 75., 76., 77., 78., 79., 80., 81., 82., 83., + 3 84., 85., 86., 87., 88., 89., 90., 91.17535, 92., 93., + 4 94., 95., 96., 97., 98., 99., 100., 101., 102., 103., + 5 104., 105., 106., 107., 108., 109., 110.03056, 111., 112., 113., + 6 114., 115., 116., 117., 118., 119., 120., 121., 122., 123., + 7 123.92928, 125., 126., 127., 128., 129., 130., 131., 132., 133., + 8 134., 135., 136., 137., 138., 139., 140., 141., 142., 143., + 9 144.,144.43391,145.,146.,147.,148.,149.,150.,151.,151.43485, + A 152., 153., 154., 155., 156., 157., 158., 159., 160., 161., + 1 162.15072,163.,164.,165.,166.,167.,167.40282,168.,169.,170., + 2 171., 172., 173., 174., 175., 176., 177., 178., 179., 180., + 3 181., 182., 183., 184., 185., 186., 187., 188., 189., 190., + 4 191.,192.,193.,194.,195.,196.,197.46990,197.81149,199.,200., + 5 201.,202.,203.,204.,205.,206.,207.13210,207.61400,208.,209., + 6 210., 211., 212., 213., 214., 215., 216., 217., 218., 219., + 7 220., 221., 222., 223., 224., 225., 226., 227., 228., 229., + 8 230., 231., 232., 233., 234., 235., 236., 237., 238., 239., + 9 240., 241., 242., 243., 244., 245., 246., 247., 248., 249./ + DATA LITC/ + 1 250.,251.12621,251.51005,252.,253.,254.,255.,256.,257.,258., + 2 259., 260., 261., 262., 263., 264., 265., 266., 267., 268., + 3 269., 270., 271., 272., 273., 274., 275., 276., 277., 278., + 4 279., 280., 281., 282., 283., 284., 285., 286., 287., 288., + 5 289., 290., 292., 294., 296., 298., 300., 302., 304., 306., + 6 308., 310., 312., 314., 316., 318., 320., 322., 324., 326., + 7 328., 330., 332., 334., 336., 338., 340., 342., 344., 346., + 8 348., 350., 352., 354., 356., 358., 360., 362., 364., 364.70183, + 9 366., 368., 370., 372., 374., 376., 378., 380., 382., 384., + A 386., 388., 390., 392., 394., 396., 398., 400., 402., 404., + 1 406., 408., 410., 412., 414., 416., 418., 420., 422., 424., + 2 426., 428., 430., 432., 434., 436., 438., 440., 442., 444., + 3 446., 448., 450., 452., 454., 456., 458., 460., 462., 464., + 4 466., 468., 470., 472., 474., 476., 478., 480., 482., 484., + 5 486., 488., 490., 492., 494., 496., 498., 500., 502., 504., + 6 506., 508., 510., 512., 514., 516., 518., 520., 522., 524., + 7 526., 528., 530., 532., 534., 536., 538., 540., 542., 544., + 8 546., 548., 550., 552., 554., 556., 558., 560., 562., 564., + 9 566., 568., 570., 572., 574., 576., 578., 580., 582., 584./ + DATA LITD/ + 1 586., 588., 590., 592., 594., 596., 598., 600., 602., 604., + 2 606., 608., 610., 612., 614., 616., 618., 620., 622., 624., + 3 626., 628., 630., 632., 634., 636., 638., 640., 642., 644., + 4 646., 648., 650., 652., 654., 656., 658., 660., 662., 664., + 5 666., 668., 670., 672., 674., 676., 678., 680., 682., 684., + 6 686., 688., 690., 692., 694., 696., 698., 700., 702., 704., + 7 706., 708., 710., 712., 714., 716., 718., 720., 722., 724., + 8 726., 728., 730., 732., 734., 736., 738., 740., 742., 744., + 9 746., 748., 750., 752., 754., 756., 758., 760., 762., 764., + A 766., 768., 770., 772., 774., 776., 778., 780., 782., 784., + 1 786., 788., 790., 792., 794., 796., 798., 800., 802., 804., + 2 806., 808., 810., 812., 814., 816., 818., 820.58271, 822., 824., + 3 826., 828., 830., 832., 834., 836., 838., 840., 842., 844., + 4 846., 848., 850., 852., 854., 856., 858., 860., 862., 864., + 5 866., 868., 870., 872., 874., 876., 878., 880., 882., 884., + 6 886., 888., 890., 892., 894., 896., 898., 900., 902., 904., + 7 906., 908., 910., 912., 914., 916., 918., 920., 922., 924., + 8 926., 928., 930., 932., 934., 936., 938., 940., 942., 944., + 9 946., 948., 950., 952., 954., 956., 958., 960., 962., 964./ + DATA LITE/ + 1 966., 968., 970., 972., 974., 976., 978., 980., 982., 984., + 2 986., 988., 990., 992., 994., 996., 998.,1000.,1005.,1010., + 3 1015.,1020.,1025.,1030.,1035.,1040.,1045.,1050.,1055.,1060., + 4 1065.,1070.,1075.,1080.,1085.,1090.,1095.,1100.,1105.,1110., + 5 1115.,1120.,1125.,1130.,1135.,1140.,1145.,1150.,1155.,1160., + 6 1165.,1170.,1175.,1180.,1185.,1190.,1195.,1200.,1205.,1210., + 7 1215.,1220.,1225.,1230.,1235.,1240.,1245.,1250.,1255.,1260., + 8 1265.,1270.,1275.,1280.,1285.,1290.,1295.,1300.,1305.,1310., + 9 1315.,1320.,1325.,1330.,1335.,1340.,1345.,1350.,1355.,1360., + A 1365.,1370.,1375.,1380.,1385.,1390.,1395.,1400.,1405.,1410., + 1 1415.,1420.,1425.,1430.,1435.,1440.,1445.,1450.,1455.,1458.81670, + 2 1465.,1470.,1475.,1480.,1485.,1490.,1495.,1500.,1505.,1510., + 3 1515.,1520.,1525.,1530.,1535.,1540.,1545.,1550.,1555.,1560., + 4 1565.,1570.,1575.,1580.,1585.,1590.,1595.,1600.,1610.,1620., + 5 1630.,1640.,1650.,1660.,1670.,1680.,1690.,1700.,1710.,1720., + 6 1730.,1740.,1750.,1760.,1770.,1780.,1790.,1800.,1810.,1820., + 7 1830.,1840.,1850.,1860.,1870.,1880.,1890.,1900.,1910.,1920., + 8 1930.,1940.,1950.,1960.,1970.,1980.,1990.,2000.,2010.,2020., + 9 2030.,2040.,2050.,2060.,2070.,2080.,2090.,2100.,2110.,2120./ + DATA LITF/ + 1 2130.,2140.,2150.,2160.,2170.,2180.,2190.,2200.,2210.,2220., + 2 2230.,2240.,2250.,2260.,2270.,2279.40330,2290.,2300.,2310.,2320., + 3 2330.,2340.,2350.,2360.,2370.,2380.,2390.,2400.,2410.,2420., + 4 2430.,2440.,2450.,2460.,2470.,2480.,2490.,2500.,2510.,2520., + 5 2530.,2540.,2550.,2560.,2570.,2580.,2590.,2600.,2610.,2620., + 6 2630.,2640.,2650.,2660.,2670.,2680.,2690.,2700.,2710.,2720., + 7 2730.,2740.,2750.,2760.,2770.,2780.,2790.,2800.,2810.,2820., + 8 2830.,2840.,2850.,2860.,2870.,2880.,2890.,2900.,2910.,2920., + 9 2930.,2940.,2950.,2960.,2970.,2980.,2990.,3000.,3010.,3020., + A 3030.,3040.,3050.,3060.,3070.,3080.,3090.,3100.,3110.,3120., + 1 3130.,3140.,3150.,3160.,3170.,3180.,3190.,3200.,3220.,3240., + 2 3260.,3282.34820,3300.,3320.,3340.,3360.,3380.,3400.,3420.,3440., + 3 3460.,3480.,3500.,3520.,3540.,3560.,3580.,3600.,3620.,3640., + 4 3660.,3680.,3700.,3720.,3740.,3760.,3780.,3800.,3820.,3840., + 5 3860.,3880.,3900.,3920.,3940.,3960.,3980.,4000.,4020.,4040., + 6 4060.,4080.,4100.,4120.,4140.,4160.,4180.,4200.,4220.,4240., + 7 4260.,4280.,4300.,4320.,4340.,4360.,4380.,4400.,4420.,4440., + 8 4460.,4480.,4500.,4520.,4540.,4560.,4580.,4600.,4620.,4640., + 9 4660.,4680.,4700.,4720.,4740.,4760.,4780.,4800.,4820.,4840./ + DATA LITG/ + 1 4860.,4880.,4900.,4920.,4940.,4960.,4980.,5000.,5020.,5040., + 2 5060.,5080.,5100.,5120.,5140.,5160.,5180.,5200.,5220.,5240., + 3 5260.,5280.,5300.,5320.,5340.,5360.,5380.,5400.,5420.,5440., + 4 5460.,5480.,5500.,5520.,5540.,5560.,5580.,5600.,5620.,5640., + 5 5660.,5680.,5700.,5720.,5740.,5760.,5780.,5800.,5820.,5840., + 6 5860.,5880.,5900.,5920.,5940.,5960.,5980.,6000.,6020.,6040., + 7 6060.,6080.,6100.,6120.,6140.,6160.,6180.,6200.,6220.,6240., + 8 6260.,6280.,6300.,6320.,6340.,6360.,6380.,6400.,6440.,6480., + 9 6520.,6560.,6600.,6640.,6680.,6720.,6760.,6800.,6840.,6880., + A 6920.,6960.,7000.,7040.,7080.,7120.,7160.,7200.,7240.,7280., + 1 7320.,7360.,7400.,7440.,7480.,7520.,7560.,7600.,7640.,7680., + 2 7720.,7760.,7800.,7840.,7880.,7920.,7960.,8000.,8040.,8080., + 3 8120.,8160.,8200.,8240.,8280.,8320.,8360.,8400.,8440.,8480., + 4 8520.,8560.,8600.,8640.,8680.,8720.,8760.,8800.,8840.,8880., + 5 8920.,8960.,9000.,9040.,9080.,9120.,9160.,9200.,9240.,9280., + 6 9320.,9360.,9400.,9440.,9480.,9520.,9560.,9600.,9640.,9680., + 7 9720.,9760.,9800.,9840.,9880.,9920.,9960.,10000./ + DO 3 I=1,328 + 3 NBIGEND(I)=DLOG(WAVEBIG(I+1))/DLOG(1.D0+1.D0/500000.D0) + NZERO=DLOG(WAVEBIG(1))/DLOG(1.D0+1.D0/500000.D0) + NZERO=NZERO+1 + WBEGIN=(1.D0+1.D0/500000.D0)**NZERO + DO 4 I=1,328 + 4 NBIGEND(I)=NBIGEND(I)-NZERO+1 + DO 5 I=2,328 + 5 NBIGBEG(I)=NBIGEND(I-1)+1 + WRITE(6,*)WBEGIN + NBIGBEG(1)=1 + WEND=(1.D0+1.D0/500000.D0)**(NBIGEND(328)+NZERO-1) + WRITE(6,*)WEND + DO 7 I=1,328 + NN(I)=NBIGEND(I)-NBIGBEG(I) + 7 WRITE(6,6)I,NBIGBEG(I),NBIGEND(I),NN(I) + 6 FORMAT(4I10) + DO 13 I=1,1212 + 13 NLITEND(I)=DLOG(WAVELIT(I+1))/DLOG(1.D0+1.D0/500000.D0) + NZERO=DLOG(WAVELIT(1))/DLOG(1.D0+1.D0/500000.D0) + NZERO=NZERO+1 + WBEGIN=(1.D0+1.D0/500000.D0)**NZERO + DO 14 I=1,1212 + 14 NLITEND(I)=NLITEND(I)-NZERO+1 + DO 15 I=2,1212 + 15 NLITBEG(I)=NLITEND(I-1)+1 + WRITE(6,*)WBEGIN + NLITBEG(1)=1 + WEND=(1.D0+1.D0/500000.D0)**(NLITEND(1212)+NZERO-1) + WRITE(6,*)WEND + DO 17 I=1,1212 + NN(I+328)=NLITEND(I)-NLITBEG(I) + 17 WRITE(6,6)I,NLITBEG(I),NLITEND(I),NN(I) + DO 18 I=1,1540 +C 1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/20.1/30,1/60 +C 1/10,2/10,3/10,4/10,5/10,6/10.7/10.8/10,9/10,19/20,59/60,60/60 +C + NSTEPS(1,I)=1. + NSTEPS(2,I)=(NN(I)*6)/60+1 + NSTEPS(3,I)=(NN(I)*12)/60+1 + NSTEPS(4,I)=(NN(I)*18)/60+1 + NSTEPS(5,I)=(NN(I)*24)/60+1 + NSTEPS(6,I)=(NN(I)*30)/60+1 + NSTEPS(7,I)=(NN(I)*36)/60+1 + NSTEPS(8,I)=(NN(I)*42)/60+1 + NSTEPS(9,I)=(NN(I)*48)/60+1 + NSTEPS(10,I)=(NN(I)*54)/60+1 + NSTEPS(11,I)=(NN(I)*57)/60+1 + NSTEPS(12,I)=(NN(I)*59)/60+1 + NSTEPS(13,I)=NN(I)+1 + 18 CONTINUE + RETURN + END + SUBROUTINE VTTAB + COMMON /VTPROF/VTPROF(120,20),VTCENTER(20),NVT(20) + DO 9 IV=1,20 + V=IV + SUM=1. + DO 3 I=1,120 + VTPROF(I,IV)=EXP(-(I/V*2.99792458E5/500000.)**2) + SUM=SUM+VTPROF(I,IV)*2. + IF(VTPROF(I,IV).LT.1.E-5)GO TO 4 + 3 CONTINUE + 4 NVT(IV)=I + DO 5 I=1,NVT(IV) + 5 VTPROF(I,IV)=VTPROF(I,IV)/SUM + VTCENTER(IV)=1./SUM + WRITE(6,7)IV,VTCENTER(IV),(VTPROF(I,IV),I=1,NVT(IV)) + 7 FORMAT(I3,' KM/S',10F10.7/(8X,10F10.7)) + 9 CONTINUE + RETURN + END + SUBROUTINE ROSSCALC(T,P,VT) + PARAMETER (LENBUFF=3510000,MAXPROF=10000) + COMMON /BUFFER/BUFFER(LENBUFF),PROFILE(MAXPROF),B10000(10000) + COMMON /CONTIN/CONTINUUM(LENBUFF) + COMMON /NLINES/WLBEG,WLEND,RESOLU,RATIO,RATIOLG,WBEGIN, + 1 LENGTH,MLINES,IXWLBEG + ABROSS=0. +C HAVE TO ADD IN UV AND IR TAILS + LENGTH=3507859 + HKT=6.6256E-27/1.38054E-16/T + FREQ0=2.99792458E17/WBEGIN*RATIO + DO 5 I=1,LENGTH + FREQ=FREQ0*RATIO**I + EHVKT=EXP(-FREQ*HKT) + STIM=1.-EHVKT + FREQ15=FREQ*1.E-15 + BNU=1.47439E-2*FREQ15**3*EHVKT/STIM + DBDT=BNU*FREQ*HKT/T/STIM +C 5 ABROSS=ABROSS+DBDT/(BUFFER(I)+CONTINUUM(I)*1000.)*FREQ/500000. + 5 ABROSS=ABROSS+DBDT/(BUFFER(I)+CONTINUUM(I)*1000.)*FREQ + ABROSS=(4.*5.6697E-5/3.14159)*T**3/ABROSS*500000. + RETURN + END |