SUBROUTINE RTECMU
C     =================
C
C     Solution of the radiative transfer equation with Compton scattering
C     for one frequency at a time (assuming the radiation intensity in other 
C     frequencies given), for a number of specific intensities (Gaussian)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
      COMMON/OPTDPT/DT(MDEPTH)
      COMMON/AUXRTE/
     *          COMA(MDEPTH),COMB(MDEPTH),COMC(MDEPTH),VL(MDEPTH),
     *          COME(MDEPTH),U(MDEPTH),V(MDEPTH),BS(MDEPTH),
     *          AL(MDEPTH),BE(MDEPTH),GA(MDEPTH)
      DIMENSION RI(MDEPTH),DTAU(MDEPTH),ST0(MDEPTH),ali(mdepth)
c     PARAMETER (MW=10, nw=10, zero=0.)
      PARAMETER (MW=3, nw=3, zero=0.)
      DIMENSION rmu(MW),b(MW),rintmu(mw),rintpo(mw),rdwn(2*mw),
     *          rmmu(2*MW),wmmu(2*mw)
      DIMENSION RJTOT(MDEPTH),RJNUT(MDEPTH),RDJ1(MDEPTH),
     *          ABSCAD(MDEPTH),ABRAD(MDEPTH),PLTOT(MDEPTH),
     *          retot(mdepth),re1(mdepth),re2(mdepth),
     *          recm(mdepth),recm0(mdepth),scom(mdepth)
c
c     set-up Gaussian angle points
c
      call gauleg(zero,un,rmu,b,nw,mw)
      do i=1,nw
         rmmu(i)=-rmu(nw-i+1)
         rmmu(i+nw)=rmu(i)
         wmmu(i)=b(nw-i+1)
         wmmu(i+nw)=b(i)
      end do
C
      SUMW=0.
      DO ID=1,ND
         RJTOT(ID)=0.
         RJNUT(ID)=0.
         ABSCAD(ID)=0.
         ABPLAD(ID)=0.
         ABRAD(ID)=0.
         pltot(id)=0.
         retot(id)=0.
         re1(id)=0.
         re2(id)=0.
         recm(id)=0.
         recm0(id)=0.
      END DO
C
c     --------------------- loop over frequencies
c
      do ij=1,nfreq
      IJI=NFREQ-KIJ(IJ)+1
      FR=FREQ(IJ)
      xcomp=xcon*fr
      CALL OPACF1(IJ)
      if(icompt.gt.0) then
         CALL RTECF0(IJ)
       else
         do id=1,nd
            if(id.lt.nd) dt(id)=deldmz(id)*(absot(id+1)+absot(id))
            comb(id)=elec(id)*sige/abso1(id)
            vl(id)=emis1(id)/abso1(id)
            coma(id)=0.
            comc(id)=0.
           bs(id)=0.
         end do
      end if
c
      SUMW=SUMW+W(IJ)
      do id=1,nd
         x0=elec(id)*sige/abso1(id)
         vl(id)=emis1(id)/abso1(id)
         st0(id)=vl(id)+(comb(id)+bs(id))*rad(iji,id)
         RDJ1(ID)=0.
         ABSCAD(ID)=ABSCAD(ID)+SCAT1(ID)*W(IJ)
         scom(id)=(comb(id)-x0*(un-two*xcomp)+bs(id))*rad(iji,id)
      end do
      rdh1=0.
      rdhd=0.
c
      if(iji.gt.1) then
         do id=1,nd
            st0(id)=st0(id)+coma(id)*rad(iji-1,id)
            scom(id)=scom(id)+coma(id)*rad(iji-1,id)
         end do
      end if
      if(iji.lt.nfreq) then
         do id=1,nd
            st0(id)=st0(id)+comc(id)*rad(iji+1,id)
            scom(id)=scom(id)+comc(id)*rad(iji+1,id)
         end do
      end if
c
      if(ifz0.lt.0) then
         fr15=fr*1.d-15
         bnu=bn*fr15*fr15*fr15
         pland=0.
         x=hk*fr/temp(nd)
         ex=exp(-x)
         pland=bnu*ex/(un-ex)
         dplan=0.
         x=hk*fr/temp(nd-1)
         ex=exp(-x)
         dplan=bnu*ex/(un-ex)
         dplan=(pland-dplan)/dt(nd-1)
      end if
C
c     --------------------- loop over angles
c
      do i=1,2*nw
         do id=1,nd-1
            dtau(id)=dt(id)/abs(rmmu(i))
         end do
c
c        boundary conditions
c
         rup=extint(ij,i)
C
C        diffusion approximation for semi-infinite atmospheres
C
         if(ifz0.lt.0) rdown=pland+rmmu(i)*dplan
c
c        the case of finite slab - irradiation of the back side
c
         if(rmmu(i).gt.0.) rdown=rdwn(nw-i+1)
c
c        solution of the transfer equation
c
         call rtesol(dtau,st0,rup,rdown,rmmu(i),ri,ali)
c
         if(rmmu(i).gt.0.) then
            if(ri(1).lt.1.e-35) ri(1)=1.e-35
            rintmu(i-nw)=ri(1)
            rintpo(i-nw)=0.
            rdh1=rdh1+rmmu(i)*wmmu(i)*ri(1)*half
         end if
         rdwn(i)=ri(nd)
         rdhd=rdhd+abs(rmmu(i)*wmmu(i))*ri(nd)*half
         DO ID=1,ND
            RDJ1(ID)=RDJ1(ID)+WMMU(I)*RI(ID)*HALF
         END DO
      end do
C
c     --------------------- end of loop over angles
c
      BBN=1.4743E-2*(FR*1.E-15)**3
      DO ID=1,ND
         pla=0.
         x=hk*fr/temp(nd)
         ex=exp(-x)
         pla=bbn*ex/(un-ex)*w(ij)
         RJTOT(ID)=RJTOT(ID)+RDJ1(ID)*W(IJ)
         RJNUT(ID)=RJNUT(ID)+RDJ1(ID)*FREQ(IJ)*W(IJ)
         ABRAD(ID)=ABRAD(ID)+RDJ1(ID)*W(IJ)*(ABSO1(ID)-SCAT1(ID))
         ABPLAD(ID)=ABPLAD(ID)+PLA*(ABSO1(ID)-SCAT1(ID))
         PLTOT(ID)=PLTOT(ID)+PLA
         retot(id)=retot(id)+abso1(id)*(st0(id)-rdj1(id))*w(ij)
         re1(id)=re1(id)+(abso1(id)-scat1(id))*rdj1(id)*w(ij)
         re2(id)=re2(id)+emis1(id)*w(ij)
         recm(id)=recm(id)+(st0(id)-vl(id)-
     *            scat1(id)/abso1(id)*rdj1(id))*w(ij)
         recm0(id)=recm0(id)+scom(id)*w(ij)
      END DO
c
c     wll=2.997925e18/fr
c     WRITE(14,641) wll,rdh1,(RINTMU(I),RINTPO(I),I=1,NW)
      END DO  
C
c     --------------------- end of loop over frequencies
c
c 641 FORMAT(1H ,f15.3,1pe15.5/(1P5E15.5))
c
      tautot=dm(nd)*elec(nd)*sige/dens(nd)
      DO ID=1,ND
         ABSCAD(ID)=elec(ID)*sige/DENS(ID)
         ABRAD(ID)=ABRAD(ID)/DENS(ID)/RJTOT(ID)
         ABPLAD(ID)=ABPLAD(ID)/DENS(ID)/PLTOT(ID)
         XNU=RJNUT(ID)/RJTOT(ID)
         re1(id)=re1(id)/dens(id)
         re2(id)=re2(id)/dens(id)
         retot(id)=retot(id)/dens(id)
         taurr=dm(id)*abscad(id)
         xl=abplad(id)*(temp(id)/teff)**4
         xr1=0.75*(1./sqrt(3.)+taurr*(1.-0.5*taurr/tautot))
         xr3a=4.*temp(id)*ycon
         xr3b=xnu*xcon
         xr3=xr3a-xr3b
         xr4=abscad(id)*xr3
         xx1=xr1*(abrad(id)-xr4)
         xx2=0.25/dm(nd)
         xr=xx1+xx2
         xtj=sig4p*4.*teff**4*xr1
         XH1=ABPLAD(ID)*PLTOT(ID)
         XH2=ABRAD(ID)*RJTOT(ID)
         XH12=XH1-XH2
         XH3=XR4*RJTOT(ID)
         XH123=XH12+XH3
         XHR=SIG4P*TEFF**4/DM(ND)
      END DO
      RETURN
      END