c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/intsumpol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine intsum
c
c     last modified 03/07/95...dave flittner
c     purpose: Spherical correction to the integration of the out-going
c     beam for each view angle.  Lines changed denoted by !def
c
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn.
c
c     last mod. Mar. 14, 2000...def
c     purpose: use definition of itmax as # of multiple scattering events.
c     Also change sb and sbp to single elements.
c
      implicit none
      include 'parameter.inc'
      include 'consts.cmn'
      include 'switches.cmn'
      include 'contrl.cmn'
      include 'totals.cmn'
      include 'out.cmn'
      include 'in.cmn'
      include 'depolt.cmn'
      include 'buff1.cmn'
      include 'buff2.cmn'
      include 'buff3.cmn'
      include 'buff4.cmn'
      include 'crefdir.cmn'
c local
      real*8 tt(max_scan)
      real*8 ttl(max_scan),ttr(max_scan),temp(max_scan)
      real*8 xizero,xi1,xi2,tfunc,gistl,gistr
      integer*4 i,j,it,k
c
c
c
c
c**** compute total reflected intensity at top of atmosphere ristar
c**** loop over albedos and polar look anglesc
      do 100 i=1,nalb
c
c**** compute ig, ground reflected radiation
c
         tensig(i)=(fs+gg)*qr(i)
         do 200 j=1,imu
c
            do 124 it=1,itmax+2
               ttz(it,j)=(fs+ggz(it))*(tnstrz(it,j)+refdir(j))       !def
124         continue
            tt(j)=(fs+gg)*(tnstr(j)+refdir(j))                    !def
            ttl(j)=(fs+gg)*(tnstrl(j)+0.50d0*refdir(j))           !def     
            ttr(j)=(fs+gg)*(tnstrr(j)+0.50d0*refdir(j))           !def
            temp(j)=ttl(j)+ttr(j)    
            if(jprint(1).eq.1)write(33,123)fs,gg,refdir(j),tnstr(j),
     &                tt(j),temp(j)
 123  format('intsum...fs,gg,refdir,tnstr(j),tt,temp',6f8.4)
c**** compute sum of direct (refdir) and diffuse (tnstr) reflected inten
c
            ristar(j,i)=tensig(i)*(tnstr(j)+refdir(j))            !def
  200    continue
  100 continue
c
c**** sum ristar and totint to calculate the total radiation (scattered
c**** and reflected (izero+i1+i2+istar) at top of atmosphere.
c
      do 500 i=1,imu
         do 400 j=1,iazmth
            do 300 k=1,nalb
               total(i,j,k)=totint(i,j)+ristar(i,k)
  300       continue
  400    continue
  500 continue
c
c     compute degree of polarization
c
      do i=1,imu
         do j=1,iazmth
            do k=1,nalb
               gistl=qr(k)*ttl(i)
               gistr=qr(k)*ttr(i)
               eittl(i,j,k)=eitl(i,j)+gistl
               eittr(i,j,k)=eitr(i,j)+gistr
               eitotz(i,j,k)=eittl(i,j,k)+eittr(i,j,k)
               polz(i,j,k)=dsqrt((eittl(i,j,k)-eittr(i,j,k))**2
     1                     +eitu(i,j)**2)/eitotz(i,j,k)
               if (eittl(i,j,k) .gt. eittr(i,j,k))
     *                             polz(i,j,k) = -polz(i,j,k)
               if(jprint(1).eq.1)then
                 write(33,'(''intsum...eittl,eittr,eitu'',3f9.5)')
     *             eittl(i,j,k), eittr(i,j,k), eitu(i,j)
                 write(33,135)i,j,k,tnstrl(i),tnstrr(i),qr(k),ttl(i),
     1               ttr(i),eitotz(i,j,k),polz(i,j,k)
               endif
            enddo
        enddo
      enddo
  135 format('im,jaz,kalb,tnstrl,tnstrr,qr,ttl,ttr,eitotz,polz'/
     1        3i4,7f9.5)

      if(jprint(1).eq.1)then
c     print izero, i1,i2 etc as func. of sza and view angle
c
        write(33,132)
132   format('theta0','theta',t15,'i0',t25,'i1',t35,'i2',t45,'t',
     1        t55,'sb',t65,'rho',t73,'ipol'/)
        do i=1,imu
           tfunc=tt(i)/pi
           xizero=eizero(i,1)/pi
           xi1=eiaz1(i,1)/pi
           xi2=eiaz2(i,1)/pi
           write(33,133)thnot(imuz),thta(i),xizero,xi1,xi2,tfunc,
     1sb,rhon,ipol
133        format(2f5.0,1x,6f10.5,1x,i3)
        enddo
        write(33,134)
134   format(/ )
      endif
      return
      end