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