SUBROUTINE INTHYD(W0,X0,Z0,IWL,ILINE)
C
C     Interpolation in temperature and electron density from the
C     hydrogen odening tables to the actual valus of
C     temperature and electron density
C
      INCLUDE 'PARAMS.FOR'
      PARAMETER (TWO=2.)
      DIMENSION ZZ(3),XX(3),WX(3),WZ(3)
C
      NX=3
      NZ=3
      NT=NTH(ILINE)
      NE=NEH(ILINE)
      BETA=WL(IWL,ILINE)/FXK
      IF(ILEMKE.EQ.1) THEN
         BETA=WL(IWL,ILINE)/XK
         NX=2
         NZ=2
      END IF
C
C     for values lower than the lowest grid value of electron density
C     the profiles are determined by the approximate expression
C     (see STARKA); not by an extrapolation in the HYD tables which may
C     be very inaccurate
C
      IF(Z0.LT.XNE(1,ILINE)*0.99.OR.Z0.GT.XNE(NE,ILINE)*1.01) THEN
         CALL DIVSTR(A,DIV)
         W0=STARKA(BETA,A,DIV,TWO)*DBETA
         W0=LOG10(W0)
         GO TO 500
      END IF
C
C     Otherwise, one interpolates (or extrapolates for higher than the
C     highes grid value of electron density) in the HYD tables
C
      DO IZZ=1,NE-1
         IPZ=IZZ
         IF(Z0.LE.XNE(IZZ+1,ILINE)) GO TO 20
      END DO
   20 N0Z=IPZ-NZ/2+1
      IF(N0Z.LT.1) N0Z=1
      IF(N0Z.GT.NE-NZ+1) N0Z=NE-NZ+1
      N1Z=N0Z+NZ-1
C
      DO 300 IZZ=N0Z,N1Z
         I0Z=IZZ-N0Z+1
         ZZ(I0Z)=XNE(IZZ,ILINE)
C
C     Likewise, the approximate expression instead of extrapolation
C     is used for higher that the highest grid value of temperature,
C     if the Doppler width expressed in beta units (BETAD) is
C     sufficiently large (> 10)
C
         IF(X0.GT.1.01*XT(NT,ILINE).AND.BETAD.GT.10.) THEN
            CALL DIVSTR(A,DIV)
            W0=STARKA(BETA,A,DIV,TWO)*DBETA
            W0=LOG10(W0)
            GO TO 500
         END IF
C
C     Otherwise, normal inter- or extrapolation
C
C     Both interpolations (in T as well as in electron density) are
C     by default the quadratic interpolations in logarithms
C
         DO IX=1,NT-1
            IPX=IX
            IF(X0.LE.XT(IX+1,ILINE)) GO TO 40
         END DO
   40    N0X=IPX-NX/2+1
         IF(N0X.LT.1) N0X=1
         IF(N0X.GT.NT-NX+1) N0X=NT-NX+1
         N1X=N0X+NX-1
         DO IX=N0X,N1X
            I0=IX-N0X+1
            XX(I0)=XT(IX,ILINE)
            WX(I0)=PRF(IWL,IX,IZZ,ILINE)
         END DO
         IF(WX(1).LT.-99..OR.WX(2).LT.-99..OR.WX(3).LT.-99.) THEN
            CALL DIVSTR(A,DIV)
            W0=STARKA(BETA,A,DIV,TWO)*DBETA
            W0=LOG10(W0)
            GO TO 500
          ELSE
            WZ(I0Z)=YINT(XX,WX,X0)
         END IF
  300 CONTINUE
      W0=YINT(ZZ,WZ,Z0)
  500 CONTINUE
      RETURN
      END