SUBROUTINE RYBCHN(CHANGT)
C     =========================
C
C     handling relative changes in the Rybicki formalism
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
      common/rybpgs/CS(MDEPTH),PRAD2D(MDEPTH),F1HE
      DIMENSION CHANGT(MDEPTH),TMPOLD(MDEPTH)
C
      CHMX=0.
      CHN=0.
      JJR=0
      DPLP=DPSILT-UN
      DPLM=UN/DPSILT-UN
      NRE=NFREQE+1
      IF(ITER.EQ.1) WRITE(9,800)
      DO ID=ND,1,-1
         CHT=CHANGT(ID)/TEMP(ID)
         CHAN=CHT
         IF(CHAN.LE.DPLM) CHAN=DPLM
         IF(CHAN.GT.DPLP) CHAN=DPLP
         TMPOLD(ID)=TEMP(ID)
         TEMP(ID)=TEMP(ID)*(CHAN+UN)
         PSI0(NRE)=TEMP(ID)
         PSY0(NRE,ID)=PSI0(NRE)
         WRITE(9,801) ITER,ID,CHT,chn,chN,CHN,chT,jjr,jjr
         IF(ABS(CHT).GT.CHMX) CHMX=ABS(CHT)
      END DO
c
      if(ioptab.le.-2) go to 10
      if(nretc.lt.0) then
         do id=-nretc,1,-1
            temp(id)=temp(id+1)
            PSY0(NRE,ID)=PSY0(NRE,ID+1)
         end do
      end if
c   
      if(idisk.eq.0) then
         if(ifprad.gt.0) then
            DO ID=1,ND
               T=TEMP(ID)
               dtod=temp(id)/tmpold(id)-un
               if(dtod.gt.0.2) dtod=0.2
               if(dtod.lt.-0.2) dtod=-0.2
               gfac=un+4.*dtod
               if(idisk.eq.0) then
                  if(id.eq.1) then
                     pgs(id)=dm(id)*(grav-grd(id)*gfac)
                   else
                     pgs(id)=pgs(id-1)+grav*(dm(id)-dm(id-1))-
     *                       pck*grd(id)*gfac
                  end if
               end if
c              if(iprybc.gt.0) 
c              write(6,603) iter,id,grd(id),gfac,pgs(id),dtod,temp(id)
c 603          format(2i5,1p6e13.5)
c
               AN=PGS(ID)/BOLK/T
               CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
               RHO=WMM(ID)*(AN-ANE)
               DENS(ID)=RHO
               ELEC(ID)=ANE
            END DO
          else
            DO ID=1,ND
               T=TEMP(ID)
               PGS(ID)=DM(ID)*GRAV
               AN=PGS(ID)/BOLK/T
               CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
               RHO=WMM(ID)*(AN-ANE)
               DENS(ID)=RHO
               ELEC(ID)=ANE
            END DO
         end if
       else
          pgpre=pgs(1)
          id=1
         do id=2,nd
           pgpre=pgs(id)
            dtod=temp(id)/tmpold(id)-un
            if(dtod.gt.0.2) dtod=0.2
            if(dtod.lt.-0.2) dtod=-0.2
            gfac=un+4.*dtod
            grv=(dm(id)-dm(id-1))*qgrav*(zd(id)+zd(id-1))*half
            pgs(id)=pgs(id-1)-pck*grd(id)*gfac+grv
            pgs0(id)=pgs(id)
         end do
         pgs0(1)=pgs(1)
c
         itpg=0
    5    continue
         z1=zd(1)
         itpg=itpg+1
         call pgset(1)
c
         DO ID=1,ND
            T=TEMP(ID)
            AN=ANTP(ID)
            CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
            RHO=WMM(ID)*(AN-ANE)
            DENS(ID)=RHO
            ELEC(ID)=ANE
            PGS(ID)=BOLK*T*AN
            pgs0(id)=pgs(id)
         END DO
c
c        recomputing the z-distance
c
         do id=nd-1,1,-1
            ddp=(dm(id+1)-dm(id))*half
            zd(id)=zd(id+1)+ddp/dens(id+1)+ddp/dens(id)
         end do
c
         if(abs((zd(1)-z1)/z1).lt.1.e-3.or.itpg.gt.5) go to 8
c
         do id=1,nd
            cs(id)=pgs(id)/dens(id)/temp(id)
         end do
            
         hr1=grd(1)/qgrav
         hg1=sqrt(two*cs(1)*temp(1)/qgrav)
         x=(zd(1)-hr1)/hg1
         IF(X.LT.3.) THEN
            IF(X.LT.0.) X=0.
            F1HE=8.86226925D-1*EXP(X*X)*ERFCX(X)
          ELSE
            F1HE=HALF*(UN-HALF/X/X)/X
         END IF
         if(abs((zd(1)-z1)/z1).lt.1.e-4.or.itpg.gt.5) go to 8
         go to 5
    8    continue
      end if
   10 continue
C
C     STOP if changes become too large
C
      IF(ITER.NE.1 .AND. ABS(CHMX).GT.1.D16) THEN
         WRITE(6,610) ITER,CHMX
         WRITE(10,610) ITER,CHMX
         STOP
      END IF
C
C     Finally, set up quantity LFIN that indicates whether or not
C     this iteration of complete linearization is the last one
C
      LFIN=ABS(CHMX).LE.CHMAX.OR.ITER.GE.NITER
c
  610 FORMAT(' **** STOP in RYBSOL after ITER',I4,/,
     * ' Max change:',1PE12.2)
  800    FORMAT(' RELATIVE CHANGES OF VECTOR PSI'/
     *   '  ITER  ID    TEMP        NE       POP     RAD      MAXIMUM',
     *   '   ilev   ifr',/)
  801    FORMAT(2I5,1P5D10.2,2I5)
      RETURN
      END