151 lines
4.8 KiB
Fortran
151 lines
4.8 KiB
Fortran
SUBROUTINE READPH
|
|
C =================
|
|
C
|
|
C Auxiliary routine for LINSET - read table of detailed
|
|
C photoinization cross-section from unit IPHT1,
|
|
C and interpolate to the set of current wavelengths (WLAM)
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
INCLUDE 'SYNTHP.FOR'
|
|
INCLUDE 'LINDAT.FOR'
|
|
COMMON/PHOTCS/PHOT(MFRQ,MPHOT),WPHT0,WPHT1,APHT(MPHOT),
|
|
* EPHT(MPHOT),GPHT(MPHOT),JPHT(MPHOT),
|
|
* NPHT
|
|
DIMENSION PHT0(MPHOT),PHT1(MPHOT),IPHT(MPHOT),IEND(MPHOT),
|
|
* IFILE(MPHOT),NELEM(MPHOT),INDEX(MPHOT,MPHOT)
|
|
PARAMETER (IPHT0=57)
|
|
SAVE IPHT,IEND,NELEM,INDEX,NUMFIL
|
|
C
|
|
C initialization - read basic information about files where the
|
|
C cross-sections are stored,
|
|
C and basic parameters for starting levels
|
|
C
|
|
IF(IBLANK.LE.1) THEN
|
|
NPHT=0
|
|
IPHT1=0
|
|
NUMFIL=0
|
|
DO 10 IJ=1,MFRQ
|
|
DO 10 I=1,MPHOT
|
|
10 PHOT(IJ,I)=0.
|
|
READ(IPHT0,*,END=50,err=50) NPHT
|
|
IF(NPHT.LE.0) RETURN
|
|
npht1=npht
|
|
READ(IPHT0,*,END=50) (IPHT(I),I=1,NPHT)
|
|
READ(IPHT0,*,END=50) (APHT(I),I=1,NPHT)
|
|
READ(IPHT0,*,END=50) (EPHT(I),I=1,NPHT)
|
|
READ(IPHT0,*,END=50) (GPHT(I),I=1,NPHT)
|
|
READ(IPHT0,*,END=50) (JPHT(I),I=1,NPHT)
|
|
C
|
|
C determination of the number of files (NFILE) and the
|
|
C partitioning of the individual cross-section to the corresponding
|
|
C files
|
|
C
|
|
NUMFIL=1
|
|
IFILE(1)=1
|
|
NELEM(1)=1
|
|
INDEX(1,1)=1
|
|
IF(NPHT.GT.1) THEN
|
|
DO 30 I=2,NPHT
|
|
DO 20 J=1,I-1
|
|
IF(IPHT(I).EQ.IPHT(J)) THEN
|
|
IFILE(I)=IFILE(J)
|
|
NELEM(IFILE(I))=NELEM(IFILE(I))+1
|
|
INDEX(IFILE(I),NELEM(IFILE(I)))=I
|
|
GO TO 30
|
|
END IF
|
|
20 CONTINUE
|
|
NUMFIL=NUMFIL+1
|
|
IFILE(I)=NUMFIL
|
|
NELEM(NUMFIL)=1
|
|
INDEX(NUMFIL,1)=I
|
|
30 CONTINUE
|
|
END IF
|
|
DO 40 IFIL=1,NUMFIL
|
|
IEND(IFIL)=0
|
|
40 CONTINUE
|
|
END IF
|
|
50 IF(NUMFIL.LE.0) RETURN
|
|
c
|
|
C loop over individual files containing the photoionization data
|
|
C
|
|
DO 300 IFIL=1,NUMFIL
|
|
IF(IEND(IFIL).EQ.1) GO TO 200
|
|
IF(IEND(IFIL).EQ.2) GO TO 300
|
|
NPHT1=NELEM(IFIL)
|
|
IPHT1=IPHT(INDEX(IFIL,1))
|
|
IF(IBLANK.LE.1) THEN
|
|
110 READ(IPHT1,*,END=200) WPHT1,(PHT1(I),I=1,NPHT1)
|
|
IF(WPHT1.LT.WLAM(1)) GO TO 110
|
|
BACKSPACE(IPHT1)
|
|
BACKSPACE(IPHT1)
|
|
READ(IPHT1,*,END=200) WPHT0,(PHT0(I),I=1,NPHT1)
|
|
ELSE
|
|
BACKSPACE(IPHT1)
|
|
BACKSPACE(IPHT1)
|
|
READ(IPHT1,*,END=200) WPHT0,(PHT0(I),I=1,NPHT1)
|
|
READ(IPHT1,*,END=200) WPHT1,(PHT1(I),I=1,NPHT1)
|
|
END IF
|
|
DW=WPHT1-WPHT0
|
|
A1=(WPHT1-WLAM(3))/DW
|
|
A2=(WLAM(3)-WPHT0)/DW
|
|
DO 130 I=1,NPHT1
|
|
INDX=INDEX(IFIL,I)
|
|
PHOT(1,INDX)=0.
|
|
PHOT(2,INDX)=0.
|
|
PHOT(3,INDX)=(A1*PHT0(I)+A2*PHT1(I))*1.E-18
|
|
DO 130 IJ=4,MFRQ
|
|
PHOT(IJ,INDX)=0.
|
|
130 CONTINUE
|
|
DO 190 IJ=4,MFRQ
|
|
IF(WLAM(IJ).LE.WPHT1) THEN
|
|
A1=(WPHT1-WLAM(IJ))/DW
|
|
A2=(WLAM(IJ)-WPHT0)/DW
|
|
DO 140 I=1,NPHT1
|
|
INDX=INDEX(IFIL,I)
|
|
PHOT(IJ,INDX)=(A1*PHT0(I)+A2*PHT1(I))*1.E-18
|
|
140 CONTINUE
|
|
ELSE
|
|
WPHT0=WPHT1
|
|
DO 150 I=1,NPHT1
|
|
150 PHT0(I)=PHT1(I)
|
|
IFSML=0
|
|
160 READ(IPHT1,*,END=180) WPHT1,(PHT1(I),I=1,NPHT1)
|
|
IF(WPHT1.LT.WLAM(IJ)) THEN
|
|
IFSML=1
|
|
GO TO 160
|
|
END IF
|
|
IF(IFSML.EQ.1) THEN
|
|
BACKSPACE(IPHT1)
|
|
BACKSPACE(IPHT1)
|
|
READ(IPHT1,*,END=180) WPHT0,(PHT0(I),I=1,NPHT1)
|
|
READ(IPHT1,*,END=180) WPHT1,(PHT1(I),I=1,NPHT1)
|
|
END IF
|
|
DW=WPHT1-WPHT0
|
|
A1=(WPHT1-WLAM(IJ))/DW
|
|
A2=(WLAM(IJ)-WPHT0)/DW
|
|
DO 170 I=1,NPHT1
|
|
INDX=INDEX(IFIL,I)
|
|
PHOT(IJ,INDX)=(A1*PHT0(I)+A2*PHT1(I))*1.E-18
|
|
170 CONTINUE
|
|
END IF
|
|
GO TO 190
|
|
180 IEND(IFIL)=1
|
|
DO 185 I=1,NPHT1
|
|
INDX=INDEX(IFIL,I)
|
|
PHOT(IJ,INDX)=0.
|
|
185 CONTINUE
|
|
190 CONTINUE
|
|
PHOT(1,INDX)=PHOT(3,INDX)
|
|
PHOT(2,INDX)=PHOT(MFRQ,INDX)
|
|
GO TO 300
|
|
200 IEND(IFIL)=2
|
|
DO 210 IJ=1,MFREQ
|
|
DO 210 I=1,NELEM(IFIL)
|
|
INDX=INDEX(IFIL,I)
|
|
PHOT(IJ,INDX)=0.
|
|
210 CONTINUE
|
|
300 CONTINUE
|
|
RETURN
|
|
END
|