SUBROUTINE RESOLV
C     =================
C
C     Control procedure for the formal solution, i.e. all calculations
C     between two consecutive iterations of complete linearization
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      common/icnrsp/iconrs
      DIMENSION PGR(MLVEXP)
C
C     Initialization - procedure INILAM:
C
      ilam=0
      CALL INILAM
      if(ioptab.lt.0.or.ioptab.gt.0) call rayset
      call prd(0)
      NLAMBD=NITLAM(ITER)
      IF(NLAMBD.LE.0) GO TO 80
      IF(LFIN.AND.NITER.GT.0) NLAMBD=1
      LAC2P=.FALSE.
      IACC0P=IACPP-3
C
C     solution of the transfer equation with Compton scattering
C
      if(icompt.ne.0.and.iter.eq.1) then
         CALL OPAINI(1)
         ilam=0
         do ij=1,nfreq
            call opacf1(ij)
            call rtefr1(ij)
         end do
         CALL RTECOM
      end if
C
      IF(ITER.LE.1.and.ioptab.eq.0) CALL LINSEL
C
C     Set of NLAMBD procedures, called overall "lambda" iterations
C     (not to be confused with ordinary lambda iterations)
C     Each "lambda" iteration contains:
C
      DO ILAM=1,NLAMBD
         CALL OPAINI(1)
         if(icompt.ne.0.and.ilam.gt.1) CALL RTECOM
C
C        Radiative rates in all transitions, in all depths
C
         IF(IFPREC.EQ.0) THEN
            CALL RATES1(0)
          ELSE
            CALL RATSP1
         END IF
c
         call prd(0)
C
C ****** evaluation of the new populations,
C        using all the previously calculated radiative rates
C
         DO ID=1,ND
c           CALL STEQEQ(ID,POP,0)
            CALL STEQEQ(ID,POP,1)
            CALL NEWPOP(ID,POP)
            IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
         END DO
C
         if(iprind.eq.2) call output
C
C ****** acceleration of convergence (if required)
C
         IF(IACPP.GT.0) THEN
            CALL ACCELP
         END IF
C
         call lucy
      END DO
C
   80 IF(ITER.EQ.1.OR.LFIN) CALL ROSSTD(0)
C
      CALL OUTPUT
C
C     in case of convection is considered: call pzeval to evaluate pressures
C     and the logarithmic gradient of the total pressure
C
      if(iter.le.nitzer) call pzert
C
      if((iheso6.ne.0.or.HMIX0.GT.0.).and.init.eq.1) CALL PZEVAL
      call radpre
      CALL TIMING(1,ITER)
      ipng=1
      if(iacd.gt.0) ipng=mod((iter-iacc),iacd)
      if(ipng.eq.0 .and. iter.ge.iacc .and. lres2) goto 90
c
c     call prnt
      if(hmix0.eq.0.) then
         write(6,611) iter-1
  611    format(/'**  CONVECTIVE FLUX: RESOLV; GLOBAL ITERATION =',I3/)
         call conout(1,ipconf)
       else if(hmix0.gt.0.) then
         if(iconre.gt.0.and.iter.le.iconre.and.iter.ge.iconrs) 
     *   call conref
         IF(ipconf.gt.0.or.(ipconf.eq.0.and.lfin)) then
            WRITE(6,611) ITER-1
            CALL CONOUT(1,1)
         END IF
      end if
C
C     evaluate necessary ALI aprameters
C
c     call prnt
      CALL OPAINI(0)
c     call prnt
      if(icompt.ne.0.and.ilam.gt.1) CALL RTECOM
      IF(KANT(ITER).EQ.1.OR.LFIN) THEN
         CALL ALISK2
       ELSE
         IF(IRDER.EQ.0) THEN
            CALL ALIST1
          ELSE
            CALL ALIST2
         END IF
      END IF
c
C     if IFPOPR=2 - evaluate new populations
C
      if(ifpopr.eq.2) then
         DO ID=1,ND
            CALL STEQEQ(ID,POP,1)
            IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
         END DO
      END IF
   90 CONTINUE
      id=1
      do ij=1,nfreqe
         absoe1(ij)=absoex(ij,id)
      end do
C
      if(ihecor.ge.-2.and.izscal.eq.0) then
      if(inzd.gt.0.or.(idisk.eq.1.and.ifryb.gt.0)) then
         if(iheso6.eq.0) then
            CALL PZEVLD
          else
            CALL HESOL6
         end if
      end if
      end if
      if(izscal.eq.1) call dmeval
c
      if(ifryb.gt.0) call rybheq
C
C     Output of condensed model atmosphere to file 7
C     This file can serve as input of initial model atmosphere for
C     another run of the program.
C
      CALL OUTPUT
c     call prnt
C
C     Output of computed model atmosphere - standard output file
C
      IF(LFIN) THEN
         IF(.NOT.LTE) CALL PRINC
         CALL OUTPRI
         IF(ICOOLP.NE.0.OR.IPOPAC.NE.0) CALL COOLRT
         CALL RECHCK
         IF(ICHCKP.NE.0) CALL CHCKSE
         IF(INTENS.GT.0) CALL RTEINT
         if(icompt.gt.0) then
            call rtecmu
            call opaini(0)
            do ij=1,nfreq
               call opacf1(ij)
               call taufr1(ij)
            end do
         end if
         RETURN
      END IF
C
C     The final part - to store previously calculated
C     mean intensities  - RADEX(IJ,ID),
C     and other model parameters (total number density, temperature,
C     electron density, and populations) for
C     further use in procedure SOLVE (actual complete linearization)
C
      IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            DO ID=1,ND
               PSY0(IJ,ID)=RADEX(IJ,ID)
            END DO
         END DO
      END IF
      DO ID=1,ND
         DO II=1,NLVEXP
            PSY0(NFREQE+INSE+II-1,ID)=0.
            PGR(II)=0.
         END DO
         DO I=1,NLEVEL
            II=IABS(IIEXP(I))
            IF(II.GT.0) PGR(II)=PGR(II)+POPUL(I,ID)
         END DO
         DO III=1,NLVEXZ
            PSY0(NFREQE+INSE+III-1,ID)=PGR(INDLGZ(III))
         END DO
         TOTN(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
         IF(INRE.NE.0) PSY0(NFREQE+INRE,ID)=TEMP(ID)
         IF(INPC.NE.0) PSY0(NFREQE+INPC,ID)=ELEC(ID)
         IF(INHE.NE.0) PSY0(NFREQE+INHE,ID)=TOTN(ID)
         IF(INZD.NE.0) PSY0(NFREQE+INZD,ID)=ZD(ID)
         IF(INMP.NE.0) PSY0(NFREQE+INMP,ID)=DENS(ID)/WMM(ID)
         IF(INDL.NE.0) PSY0(NFREQE+INDL,ID)=DELTA(ID)
      END DO
C
      IF(INIT.EQ.1.AND.NATOM.GT.0) THEN 
         WRITE(6,600)
         DO ID=1,ND
            WRITE(6,601) ID,(NREFS(I,ID),I=1,NATOM)
         END DO
  600    FORMAT(//' REFERENCE LEVEL INDICES AS FUNCTIONS OF DEPTH'/,
     *          'ITER =',i4)
  601    FORMAT(' ID=',I3,2X,15I4)
      END IF
C           
      RETURN
      END