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     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, i.e. scale DEPTH coincides with DM
C            > 0  -  polynomial interpolation of the (INTRPL-1)th order
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      PARAMETER (MINPUT=MLEVEL+4)
      DIMENSION ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),POPLTE(MLEVEL),
     *          TOTN(MDEPTH),PLTE(MLEVEL,MDEPTH)
      COMMON ESEMAT,BESE,POPLTE,POPUL0(MLEVEL,MDEPTH),X(MINPUT),
     *       TEMP0(MDEPTH),ELEC0(MDEPTH),DENS0(MDEPTH),PPL0(MDEPTH),
     *       PPL(MDEPTH),DEPTH(MDEPTH),DM0(MDEPTH),DP(MDEPTH)
      COMMON/NLTPOP/PNLT(MATOM,MION,MDEPTH)
      common/quasex/iexpl(mlevel),iltot(mlevel)
C
      NUMLT=3
      IF(INMOD.EQ.2) NUMLT=4
      READ(8,*) NDPTH,NUMPAR
      READ(8,*) (DEPTH(I),I=1,NDPTH)
      ND=NDPTH
      NUMP=ABS(NUMPAR)
      DO 30 ID=1,NDPTH
         READ(8,*) (X(I),I=1,NUMP)
         TEMP(ID)=X(1)
         ELEC(ID)=X(2)
         DENS(ID)=X(3)
         TOTN(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
         CALL WNSTOR(ID)
         CALL SABOLF(ID)
         IP=NUMLT
         IF(NUMPAR.LT.0) THEN
            IP=IP+1
            TOTN(ID)=X(IP)
         END IF
         IF(INMOD.EQ.2) IP=IP+1
c
c        first compute LTE level populations for all levels,
c        i.e. explicit, semi-explisit, and quasi-explicit
c
         NLEV0=NLEVEL
         TEMP(ID)=X(1)
         ELEC(ID)=X(2)
         DENS(ID)=X(3)
         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,ipri)
          else
            if(imode.gt.-2) then
               DO IAT=1,NATOM
                  ATTOT(IAT,ID)=DENS(ID)/WMM(ID)/YTOT(ID)*ABUND(IAT,ID)
               END DO
             else
               DO IAT=1,NATOM
                  ATTOT(IAT,ID)=DENS(ID)/WMM(1)/YTOT(1)*ABUND(IAT,1)
               END DO
            end if
         end if
         CALL WNSTOR(ID)
         CALL SABOLF(ID)
         CALL RATMAT(ID,ESEMAT,BESE)
         CALL LEVSOL(ESEMAT,BESE,POPLTE,NLEV0) 
         DO I=1,NLEV0
            POPUL(I,ID)=POPLTE(I)
            PLTE(I,ID)=POPLTE(I)
c           if(id.eq.1) write(6,651) i,ip,popul(i,id),plte(i,id)
         END DO
c
c        if the input file fort.8 contains also NLTE level populations
c        of b-factors, replace the LTE populations by those
c
         IF(NUMP.GT.IP) THEN
            NLEV0=NUMP-IP    
            DO I=1,NLEV0
               j=iltot(i)
               POPUL(J,ID)=X(IP+I)*RELAB(IATM(I),ID)
c              if(id.eq.1) write(6,651) i,j,x(ip+i),popul(i,id)
c 651          format('in',2i4,1p2e12.4)
            END DO
c           DO I=1,NLEV0
c              j=iltot(i)
c              if(popul(j,id).le.0.) then
c                 IE=IEL(I)
c                 N0I=NFIRST(IE)
c                 NKI=NNEXT(IE)
c                 POPUL(J,ID)=ELEC(ID)*POPUL(iltot(NKI),ID)*SBF(I)
c              end if
c           END DO
c
c           in the case the input "NLTE populations are in fact b-factors,
c           compute the real populations
c
            if(ibfac.eq.1) then
               do i=1,nlev0
                  j=iltot(i)
                  popul(j,id)=popul(j,id)*plte(j,id)
               end do
            end if
         END IF
   30 CONTINUE
C
      close(8)
c
      write(6,600)
  600 format(/' INPUT TLUSTY MODEL'/
     *        ' ------------------'/
     *         1H ,8X,'MASS',9X,'T',9X,'NE',9X,'DENS'//)
         nd=ndpth
         DO 40 ID=1,ND
            DM(ID)=DEPTH(ID)
            write(6,601) id,dm(id),temp(id),elec(id),dens(id),
     *       popul(1,id)
  601       format(i6,1pe10.3,0pf10.1,1p4e12.3)
   40    CONTINUE
C
      DO 100 ID=1,ND
         BCON=ELEC(ID)/TEMP(ID)/SQRT(TEMP(ID))*2.0706E-16
         DO 100 IONE=1,NION
            ION=IZ(IONE)
            IAT=NUMAT(IATM(NFIRST(IONE)))
            NKI=NNEXT(IONE)
            IF(ION.GT.0) PNLT(IAT,ION,ID)=POPUL(NKI,ID)/G(NKI)*BCON
  100 CONTINUE
c
c     check abundances
c
c     CALL CHCKAB
      RETURN
      END