SUBROUTINE HYDINI
C
C     Initializes necessary arrays for evaluating hydrogen line profiles
C     from the Lemke, Tremblay-Bergeron, or Schoening-Butler tables
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
c     DIMENSION WLINE(4,22)
      DIMENSION IILW(100),IIUP(100)
      CHARACTER*1 CHAR
      DATA INIT /0/
C
      IF(INIT.EQ.0) THEN
         DO I=1,4
            DO J=I+1,22
               CALL STARK0(I,J,IZZ,XK,WL0,FIJ,FIJ0)
               WLINE(I,J)=WL0
c              OSCH(I,J)=FIJ+FIJ0
            END DO
         END DO
         INIT=1
      END IF
      DO I=1,4
         DO J=1,22
            ILIN0(I,J)=0
         END DO
      END DO
C
C --------------------------------------------
C     Schoening-Butler tables - for IHYDPR < 0
C --------------------------------------------
C
      IF(IHYDPR.LT.0) THEN
      IHYDPR=67
      ILEMKE=0
      NLINE=12
c
      OPEN(UNIT=IHYDPR,FILE='./data/hydprf.dat',STATUS='OLD')
      write(6,*) ' reading Schoening-Butler tables'
C
      DO I=1,12
         READ(IHYDPR,500)
      END DO
      DO 100 ILINE=1,NLINE
C
C     read the tables, which have to be stored in file
C     unit IHYDPR (which is the input parameter in the progarm)
C
         READ(IHYDPR,501) I,J
         IF(ILINE.EQ.12) J=10
         WL0=WLINE(I,J)
         ILIN0(I,J)=ILINE
         READ(IHYDPR,*) CHAR,NWL,(WL(I,ILINE),I=1,NWL)
         READ(IHYDPR,*) CHAR,NT,(XT(I,ILINE),I=1,NT)
         READ(IHYDPR,*) CHAR,NE,(XNE(I,ILINE),I=1,NE)
         READ(IHYDPR,500)
         NWLH(ILINE)=NWL
         NWLHYD(ILINE)=NWL
         NTH(ILINE)=NT
         NEH(ILINE)=NE
C
         DO I=1,NWL
            IF(WL(I,ILINE).LT.1.E-4) WL(I,ILINE)=1.E-4
            WLHYD(ILINE,I)=LOG10(WL(I,ILINE))
         END DO
C
         DO IE=1,NE
            DO IT=1,NT
               READ(IHYDPR,500)
               READ(IHYDPR,*) (PRF(IWL,IT,IE,ILINE),IWL=1,NWL)
            END DO
         END DO
C
C        coefficient for the asymptotic profile is determined from
C        the input data
C
         XCLOG=PRF(NWL,1,1,ILINE)+2.5*LOG10(WL(NWL,ILINE))+31.5304-
     *         XNE(1,ILINE)-2.*LOG10(WL0)
         XKLOG=0.6666667*(XCLOG-0.176)
         XK=EXP(XKLOG*2.3025851)
C
         DO ID=1,ND
C
C           temperature is modified in order to account for the
C           effect of turbulent velocity on the Doppler width
C
            T=TEMP(ID)+6.06E-9*VTURB(ID)
            ANE=ELEC(ID)
            TL=LOG10(T)
            ANEL=LOG10(ANE)
            F00=1.25E-9*ANE**0.666666667
            FXK=F00*XK
            DOP=1.E8/WL0*SQRT(1.65E8*T)
            DBETA=WL0*WL0/2.997925E18/FXK
            BETAD=DBETA*DOP
C
C       interpolation to the actual values of temperature and electron
C       density. The result is stored at array PRFHYD, having indices
C       ILINE (line number: 1 for L-alpha,..., 4 for H-delta, etc.);
C                           5 for H-alpha,..., 8 for H-delta, etc.)
C       ID - depth index
C       IWL - wavelength index 
C
            DO IWL=1,NWL
               CALL INTHYD(PROF,TL,ANEL,IWL,ILINE)
               PRFHYD(ILINE,ID,IWL)=PROF
            END DO
         END DO
  100 CONTINUE
      CLOSE(IHYDPR)
C
  500 FORMAT(1X)
  501 FORMAT(12X,I1,9X,I1)
C
      IHYDPR=-IHYDPR
      RETURN
      END IF
C
C ---------------------------------
C     read Lemke or Tremblay tables
C ---------------------------------
C
      if(ihydpr.lt.20) ihydpr=ihydpr+20
      if(ihydpr.eq.21) then
         open(unit=ihydpr,file='./data/lemke.dat',status='old')
         write(6,641) ihydpr
       else if(ihydpr.eq.22) then
         open(unit=ihydpr,file='./data/tremblay.dat',status='old')
         write(6,642) ihydpr
      end if
  641 format(' -----------'/
     *       ' reading Lemke tables; ihydpr =',i3,/
     *       ' -----------')
  642 format(' -----------'/
     *       ' reading Tremblay tables; ihydpr =',i3,/
     *       ' -----------')
C
      ILEMKE=1
      READ(IHYDPR,*) NTAB
      write(6,611) ntab 
  611 format(' ntab',i4)
      DO ITAB=1,NTAB
         ILINEB=ILINE
         READ(IHYDPR,*) NLLY
         DO ILI=1,NLLY
            ILINE=ILINE+1
            READ(IHYDPR,*) I,J,ALMIN,ANEMIN,TMIN,DLA,DLE,DLT,
     *                     NWL,NE,NT
            WL0=WLINE(I,J)
            ILIN0(I,J)=ILINE
            NWLH(ILINE)=NWL
            NWLHYD(ILINE)=NWL
            NTH(ILINE)=NT
            NEH(ILINE)=NE
            iilw(iline)=i
            iiup(iline)=j
            DO IWL=1,NWL
               WL(IWL,ILINE)=ALMIN+(IWL-1)*DLA
               WLHYD(ILINE,IWL)=WL(IWL,ILINE)
               WL(IWL,ILINE)=EXP(2.3025851*WL(IWL,ILINE))
            END DO
            DO INE=1,NE
               XNE(INE,ILINE)=ANEMIN+(INE-1)*DLE
            END DO
            DO IT=1,NT
               XT(IT,ILINE)=TMIN+(IT-1)*DLT
            END DO 
         END DO
c
         DO ILI=1,NLLY         
            ILNE=ILINEB+ILI
            NWL=NWLH(ILNE)
            READ(IHYDPR,500)
            DO INE=1,NEH(ILNE)
               DO IT=1,NTH(ILNE)
                  READ(IHYDPR,*) QLT,(PRF(IWL,IT,INE,ILNE),IWL=1,NWL)
               END DO
            END DO
C
            i=iilw(ilne)
            j=iiup(ilne)
            DO ID=1,ND
               CALL HYDTAB(I,J,ID)
            END DO
         END DO
      END DO
      NLIHYD=ILNE
      CLOSE(IHYDPR)
C
      RETURN
      END