SUBROUTINE RADTEM
C     =================
C
C     determination of the radiation temperatures
C     after Schmutz (1991); inversion done by Newton-Raphson
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'WINCOM.FOR'
      common/velaux/velmax,iemoff,nltoff,itrad
      PARAMETER (CON=2.0706D-16, un=1.)
      parameter (nterad=3)
C
      DO ID=1,ND
         rx=RD(ND)/RD(ID)
c        WDIL(ID)=0.5*(1.-sqrt(1.-rx*rx))
         wdil(id)=un-sqrt(un-rx*rx)
      END DO
      DO ITRD=1,NTERAD
         if(itrad.eq.0) then
         do id=1,nd
            trad(itrd,id)=temp(id)
         end do
         else
         II=0
         JJ=0
         IF(ITRD.LE.NION) II=NFIRST(ITRD)
         IF(ITRD.LE.NION) JJ=NNEXT(ITRD)
         DO ID=1,ND
            TRAD(ITRD,ID)=TEMP(ID)
            IF(II.GT.0) THEN
c           IF(II.GT.100000) THEN
            AA=POPUL(JJ,ID)/POPUL(II,ID)*ELEC(ID)*CON
            AA=AA*G(II)/G(JJ)/WDIL(ID)/SQRT(TEMP(ID))
            TR=TEMP(ID)
            ITER=0
   10       ITER=ITER+1
            XX=ENION(II)/BOLK/TR
            DTR=(AA*EXP(XX)-TR)/(1.+XX)
            DTRR=DTR/TR
            TR=TR+DTR
            IF(ABS(DTRR).GT.1.E-3.AND.ITER.LT.100) GO TO 10
            TRAD(ITRD,ID)=TR
            END IF
         END DO
         end if
      END DO
      write(6,600)
      do id=1,nd
         write(6,601) id,temp(id),trad(1,id),trad(2,id),trad(3,id)
      end do
  600 format(/' radiation temperatures/')
  601 format(i5,4f10.1)
      RETURN
      END