SUBROUTINE RADTOT
C     =================
C
C     Evaluation of integrated radiative intensities and moments
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ODFPAR.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ITERAT.FOR'
      COMMON/OPTDPT/DT(MDEPTH)
      COMMON/SURFEX/EXTJ(MFREQ),EXTH(MFREQ)
      COMMON/TOTJHK/TOTJ(MDEPTH),TOTH(MDEPTH),TOTK(MDEPTH),
     *              RDOPAC(MDEPTH),FLOPAC(MDEPTH)
      DIMENSION SUMPL(MDEPTH)
C
C     zero the quantities
C
      DO ID=1,ND
         ABROSD(ID)=0.
         SUMDPL(ID)=0.
         ABPLAD(ID)=0.
         SUMPL(ID)=0.
         TOTJ(ID)=0.
         TOTH(ID)=0.
         TOTK(ID)=0.
         RDOPAC(ID)=0.
         FLOPAC(ID)=0.
         if(id.lt.nd) THEN
            DELDM(ID)=HALF*(DM(ID+1)-DM(ID))
            deldmz(id)=deldm(id)
            if(izscal.eq.1) deldmz(id)=half*(zd(id)-zd(id+1)) 
         end if
      END DO
      DEDM1=DM(1)/DENS(1)
C
C     loop over frequencies
C
      CALL TDPINI
      CALL OPAINI(1)
      DO IJ=1,NFREQ
         FR=FREQ(IJ)
         CALL OPACF1(IJ)
         CALL RTEFR1(IJ)
         WW=W(IJ)
         DO ID=1,ND
            PLAN=XKFB(ID)/XKF1(ID)*WW
            DPLAN=PLAN/XKF1(ID)*FREQ(IJ)*HKT21(ID)
            if(ioptab.ge.0) then
            ABROSD(ID)=ABROSD(ID)+DPLAN/ABSO1(ID)
            ABPLAD(ID)=ABPLAD(ID)+PLAN*(ABSO1(ID)-SCAT1(ID))
            RDOPAC(ID)=RDOPAC(ID)+WW*RAD1(ID)*(ABSO1(ID)-SCAT1(ID))
            else
            ar=(ABSO1(ID)-SCAT1(ID))*dens(id)
            ABROSD(ID)=ABROSD(ID)+DPLAN/(ABSO1(ID)*dens(id))
            ABPLAD(ID)=ABPLAD(ID)+PLAN*AR
            RDOPAC(ID)=RDOPAC(ID)+WW*RAD1(ID)*AR
            end if
            SUMDPL(ID)=SUMDPL(ID)+DPLAN
            SUMPL(ID)=SUMPL(ID)+PLAN
            TOTJ(ID)=TOTJ(ID)+WW*RAD1(ID)
            TOTK(ID)=TOTK(ID)+WW*RAD1(ID)*FAK1(ID)
            IF(ID.LT.ND) THEN
               FLUX1=RAD1(ID+1)*FAK1(ID+1)-RAD1(ID)*FAK1(ID)
               TOTH(ID+1)=TOTH(ID+1)+WW*FLUX1/DT(ID)
            END IF
         END DO
         WF=WW*(FH(IJ)*RAD1(1)-HEXTRD(IJ))
         TOTH(1)=TOTH(1)+WF
         if(ioptab.ge.0) then
            FLOPAC(1)=FLOPAC(1)+WF*ABSO1(1)/DENS(1)
          else
            FLOPAC(1)=FLOPAC(1)+WF*ABSO1(1)
         end if
      END DO  
c
c     Rosseland and Planck mean opacities
C
      DO ID=1,ND
         ABROSD(ID)=SUMDPL(ID)/ABROSD(ID)
         ABPLAD(ID)=ABPLAD(ID)/SUMPL(ID)
      END DO
C
c     Rosseland optical depth scale; flux mean
c
      ID=1
      TAUROS(ID)=HALF*DEDM1*ABROSD(ID)
      DO ID=2,ND         
         DTAUR=DELDM(ID-1)*(ABROSD(ID)*DENS1(ID)+
     *                      ABROSD(ID-1)*DENS1(ID-1))
         TAUROS(ID)=TAUROS(ID-1)+DTAUR
         FLOPAC(ID)=(TOTK(ID)-TOTK(ID-1))/(DM(ID)-DM(ID-1))
      END DO
c
c     final Rosseland and Planck mean opacities
C
      DO ID=1,ND
         ABROSD(ID)=ABROSD(ID)/DENS(ID)
         ABPLAD(ID)=ABPLAD(ID)/DENS(ID)
      END DO
C
      RETURN
      END