subroutine Iplus_calc(ntau,tau,source,Iplus,optthin)
*
* Calculates Iplus (outwards directed intensities) along
* a ray at all depths in a model atmosphere, by integration of 
* the source function. Input: tau is the optical depth along
* the ray (i.e. tau/mu), source is the source function ,
* both provided at ntau depth points.  BPz 19/07-2001
*
* THIS VERSION ASSUMES SPHERICAL SYMMETRY
*
      implicit none
      logical optthin
      integer i,ntau
      real Iplus(*),source(*),tau(*),ex(1000)

      if (ntau.gt.1000) stop 'Iplus_calc.f : increase ex() dimension!'
      do i=2,ntau
        ex(i)=exp(tau(i-1)-tau(i))
      enddo
      if (optthin) then
* optically thin: compute incoming radiation at taumax, 
* by integrating intensity inwards from the surface
* ASSUMES SPHERICAL SYMMETRY

        Iplus(1)=source(1)*(1.-exp(-tau(1)))
        do i=2,ntau
          Iplus(i)=Iplus(i-1)*ex(i)+
     &             (source(i-1)+source(i))*0.5*
     &              (1.-ex(i))
        enddo
      else
* optically thick: diffusion approximation
        Iplus(ntau)=source(ntau)+
     &     (source(ntau)-source(ntau-1))/(tau(ntau)-tau(ntau-1))
      endif
      do i=ntau-1,1,-1
        Iplus(i)=Iplus(i+1)*ex(i+1)+
     &           (source(i+1)+source(i))*0.5*
     &                  (1.-ex(i+1))
      enddo

      return
      end