SUBROUTINE KURUCZ(NDPTH)
C     ========================
C
C     Read an initial model atmosphere from unit 8
C      in Kurucz ATLAS' format
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      PARAMETER (MINPUT=7)
      CHARACTER KUR*15
      DIMENSION KKFIX0(MLEVEL)
      COMMON POPUL0(MLEVEL,MDEPTH),ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),
     *       TEMP0(MDEPTH),ELEC0(MDEPTH),DENS0(MDEPTH),PPL0(MDEPTH),
     *       DEPTH(MDEPTH),PPL(MDEPTH),POPLTE(MLEVEL),X(MINPUT)
      common/temlim/tfloor
      dimension a(mlevel,mlevel),b(mlevel),iifor0(mlevel)
C
      do iat=1,natom
         kkfix0(iat)=iifix(iat)
         iifix(iat)=0
      end do
      nlev0=nlevel
      LCHC0=LCHC
      LCHC=.TRUE.
      LTE0=LTE
      LTE=.TRUE.
C
      if(ifixde.gt.0) then
      READ(8,551) TEF,GRAV
      WRITE(6,600) TEF,GRAV
      READ(8,552) ND
      ND=ND-1
c     write(6,553) nd
  551 FORMAT(4X,F8.0,9X,F8.5)
  552 FORMAT(/////////////////////10X,I3/)
c 553 format(' nd',i4)
C
  600 FORMAT(1H1,' INPUT KURUCZ MODEL FOR TEFF=',F7.0,'   LOG G =',
     *       F7.2//1H ,7X,'MASS',9X,'T',9X,'P',9X,'DENS',9X,'ELEC'//)
      DO ID=1,ND
            READ(8,*) DM(ID),TEMP(ID),P,ane0,a1,a2,a3,vel,rho
            if(temp(id).lt.tfloor) temp(id)=tfloor
            CALL RHONEN(ID,TEMP(ID),RHO,AN,ANE)
c           ELEC(ID)=ANE0
            ELEC(ID)=ANE
            DENS(ID)=RHO
            TOTN(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
            an=rho/wmm(id)+ane0
            rho0=WMM(ID)*(AN-ELEC(ID))
            WRITE(6,651) ID,DM(ID),TEMP(ID),P,DENS(ID),ELEC(ID),rho0,an
         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
      end do
c     write(6,553) nd
      go to 100
      end if
  651 FORMAT(1H ,I5,1PE10.3,0PF10.1,1P5D12.3)
c
      READ(8,801) KUR,GRAVK
      READ(KUR,802) TEFFK
  801 FORMAT(A15,6X,F8.5)
  802 FORMAT(4X,F8.0)
C
      IF(KUR(1:4).NE.'TEFF')
     *  CALL QUIT(' Unit 8 is NOT a Kurucz model as expected',0,0)
      IF(ABS(TEFFK-TEFF).GT.50.) then
        ieff=int(teff)
        ieffk=int(teffk)
c        CALL QUIT(' Teff not corresponding to Kurucz model',ieff,ieffk)
      END IF
      IF(ABS(GRAVK-LOG10(GRAV)).GT.0.02) then
        irav=int(log10(grav)+0.001)
        iravk=int(gravk)
c        CALL QUIT(' Gravity not corresponding to Kurucz model',
c     *             irav,iravk)
      END IF
C
      DO WHILE(KUR(1:9).NE.'READ DECK')
         READ(8,'(A15)') KUR
      END DO
      READ(KUR,803) NDPTH
  803 FORMAT(10X,I3)
      NDPTH=NDPTH-1
      NUMPAR=3
      NLEV0=NLEVEL
      READ(8,*) TTT
      IF(NDPTH.gt.mdepth.and.intrpl.eq.0) 
     *   CALL QUIT('ndpth.gt.mdepth in KURUCZ',ndpth,mdepth)
      DO ID=1,NDPTH
         READ(8,*) (X(I),I=1,MINPUT)
         DEPTH(ID)=X(1)
         TEMP0(ID)=X(2)
         ELEC0(ID)=X(4)
         AN=X(3)/BOLK/TEMP0(ID)
         DENS0(ID)=WMM(ID)*(AN-ELEC0(ID))
         TEMP(ID)=TEMP0(ID)
         ELEC(ID)=ELEC0(ID)
         DENS(ID)=DENS0(ID)
         ANMA(ID)=DENS(ID)/WMM(ID)
         ANTO(ID)=ANMA(ID)+ELEC(ID)
         TOTN(ID)=ANTO(ID)
         T=TEMP(ID)
         IF(IFMOL.GT.0.AND.T.LT.TMOLIM) THEN
            AN=TOTN(ID)
            AEIN=ELEC(ID)
            CALL MOLEQ(ID,T,AN,AEIN,ANE,ENR,ENT,WM,1)
         END IF
         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
            POPUL0(I,ID)=POPLTE(I)
         END DO
      END DO
C
      INTRPL=0
      DO WHILE(KUR(1:5).NE.'BEGIN')
         READ(8,'(A15)',END=100,ERR=100) KUR
      END DO
      READ(8,*,END=100,ERR=100) INTRPL
C
  100 do iat=1,natom
         iifix(iat)=kkfix0(iat)
      end do
      LCHC=LCHC0
      LTE=LTE0
C
      RETURN
      END