SUBROUTINE INPMOD
C     =================
C
C     Read an initial model atmosphere from unit 8
C     File 8 contains:
C      1. NDPTH -  number of depth points in which the initial model is
C                  given (if not equal to ND, routine interpolates
C                  automatically to the set DM by linear interpolation
C                  in log(DM)
C         NUMPAR - number of input model parameters in each depth
C                  = 3 for LTE model - ie. N, T, N(electron);
C                  > 3 for NLTE model)
C      2. DEPTH(ID),ID=1,NDPTH - mass-depth points for the input model
C      3. for each depth:
C                 T   - temperature
C                 ANE - electron density
C                 RHO - mass density
C                 level populations - only for NLTE input model
C                       Number of input level populations need not be
C                       equal to NLEVEL; in that case the procedure
C                       CHANGE is called from START to calculate the
C                       remaining level populations
C
C     Note: The output file 7, which is created by this program
C           (procedure OUTPUT) has the same structure as file 8
C           and may thus be used as input to another run of the
C           program
C
C     INTRPL - switch indicating whether (and, if so, how) interpolate
C              the initial model if the depth scales for the input model
C              and the present depth scale are different
C            = 0  -  no interpolation, scale DEPTH replaces DM;
C                    i.e. DEPTH will be used as DM, regardless of which
C                    values of DM were read in subroutine START
C            > 0  -  polynomial interpolation of the (INTRPL-1)th order
C
C            < 0  -  reads different initial models (eg. -1 : Kurucz model)
C
C     If INTRPL > 0, there is an additional input from unit 8, namely
C     new depth scale DM, the one which will be used in the present run
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      PARAMETER (MINPUT=MLEVEL+5,
     *           MDEPTI=MDEPTH)
      common/eospar/anmol(600,mdepth),
     *              anato(100,mdepth),
     *              anion(100,mdepth)
      COMMON POPUL0(MLEVEL,MDEPTI),ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),
     *       TEMP0(MDEPTI),ELEC0(MDEPTI),DENS0(MDEPTI),PPL0(MDEPTI),
     *       DEPTH(MDEPTI),PPL(MDEPTH),POPLTE(MLEVEL),X(MINPUT),
     *       ZD0(MDEPTH)
      dimension a(mlevel,mlevel),b(mlevel),iifor0(mlevel)
C
      LCHC0=LCHC
      LCHC=.TRUE.
      LTE0=LTE
      LTE=.TRUE.
      NUMLT=3
      IF(IDISK.EQ.1) NUMLT=4
      IF(IFMOL.GT.0) NUMLT=NUMLT+1
      RRDIL=1.
      TEMPBD=0.
      IF(INTRPL.GE.0) THEN
         READ(8,*) NDPTH,NUMPAR
         ND=NDPTH
         IF(NDPTH.LE.0) CALL QUIT('NDPTH.LE.0 in Unit 8',ndpth,0)
         IF(NDPTH.GT.MDEPTI) 
     *      CALL QUIT('NDPTH.GT.MDEPTI in Unit 8',ndpth,mdepti)
         READ(8,*) (DM(ID),ID=1,ND)
         IDSTD=0
         NUMP=ABS(NUMPAR)
         DO ID=1,NDPTH
            READ(8,*) (X(I),I=1,NUMP)
            do i=1,nump
               if(x(i).lt.0.) x(i)=0.
            end do
            TEMP(ID)=X(1)
            ELEC(ID)=X(2)
            DENS(ID)=X(3)
            TOTN(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
            IP=3
            IF(NUMPAR.LT.0) THEN
               IP=IP+1
               TOTN(ID)=X(IP)
            END IF
            IF(IDISK.EQ.1) THEN
               IP=IP+1
               ZD(ID)=X(IP)
            END IF
            IF(TEMP(ID).LT.TEFF) IDSTD=ID
c           IF(NUMP.GT.IP.AND..NOT.LTE0.AND.ICHANG.NE.-2) THEN
            t=temp(id)
            if(ifmol.gt.0.and.t.lt.tmolim) then
               ipri=1
               aein=elec(id)
               an=totn(id)
               call moleq(id,t,an,aein,ane,enrg,entt,wm,ipri)
            end if
            IF(NUMP.GT.IP) THEN
               NLEV0=NUMP-IP
               DO I=1,NLEV0
                  POPUL(I,ID)=X(NUMLT+I)
               END DO
             ELSE
               NLEV0=NLEVEL
               CALL WNSTOR(ID)
               CALL SABOLF(ID)
               DO I=1,NLEV0
                  IIFOR0(I)=I
               END DO
               CALL RATMAT(ID,IIFOR0,-1,A,B)
               CALL LEVSOL(A,B,POPLTE,IIFOR0,NLEV0,1) 
               DO I=1,NLEV0
                  POPUL(I,ID)=POPLTE(I)
               END DO
               if(ifmol.le.0.or.t.ge.tmolim) then
                  if(n0hn.gt.0) then
                     anato(1,id)=popul(n0hn,id)
                   else
                     anato(1,id)=dens(id)/wmm(id)/ytot(id)
                  end if
                  if(iathe.gt.0) then
                     anato(2,id)=popul(n0a(iathe),id)
                   else
                     anato(2,id)=dens(id)/wmm(id)/ytot(id)*abndd(2,id)
                  end if
               end if
c

            END IF
         END DO
         ELSTD=ELEC(IDSTD)
         READ(8,*,END=10,ERR=10) INTRPL
       ELSE IF(INTRPL.GT.-10) THEN
         CALL KURUCZ(NDPTH)
         NUMPAR=3
         IF(ND.NE.NDPTH .AND. INTRPL.EQ.0)
     *      CALL QUIT('ND.NE.NDPTH in KURUCZ',nd,ndpth)
       ELSE
         CALL INCLDY(NDPTH)
         NUMPAR=3
         IF(ND.NE.NDPTH .AND. INTRPL.EQ.0)
     *      CALL QUIT('ND.NE.NDPTH in INCLDY',nd,ndpth)
      END IF
   10 LCHC=LCHC0
      LTE=LTE0
C
      IF(IDISK.EQ.1) THEN
         ZND=ZD(ND)
         IF(ZND.GT.0.) IFZ0=-1
      END IF
c
c !!!!! attention - temporary fix
c
      if(idisk.eq.1) zd(nd)=0.
C      
c
      if(ioptab.ge.0) return
      IF(IOPTAB.LT.0) RETURN
      do id=1,nd
         ptotal(id)=dm(id)*grav
         an=ptotal(id)/(bolk*temp(id))
         elec(id)=1.e-16*an
         wmm(id)=dens(id)/an
      end do
      RETURN
      END