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
      implicit real*8(a-h,o-z),integer*4(i-n)
c
      include "contrl.inc"
      include "out.inc"
      include "totals.inc"
      include "emm.inc"
      include "atmos.inc"
      include "in.inc"
      include "depolt.inc"
      include "buff1.inc"
      include "buff2.inc"
      include "buff3.inc"
      include "buff4.inc"
      include "crefdir.inc"
c
      real*8 tt(15)
      real*8 totn(15,8,11)
      real*8 temp(15)
c
c**** compute total reflected intensity at top of atmosphere ristar
c**** loop over albedos and polar look angles
c
      itmaxp=itmax+1
      maxpp=itmaxp+1
      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,maxpp
      ttz(it,j)=(fs+ggz(it))*(tnstrz(it,j)+refdir(j))       !def
124   continue
      tt(j)=(fs+gg)*(tnstr(j)+refdir(j))                    !def
      ttlp(j)=(fs+gg)*(tnstrl(j)+0.50d0*refdir(j))           !def     
      ttrp(j)=(fs+gg)*(tnstrr(j)+0.50d0*refdir(j))           !def
      temp(j)=ttlp(j)+ttrp(j)    
c      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
c
      total(i,j,k)=totint(i,j)+ristar(i,k)
      totn(i,j,k)=-100.0*dlog10(total(i,j,k))
c     if(i.eq.1 .and. j.eq.1 .and. k.eq.1)then
c       write(33,666)i,j,k,thnot(imuz),thta(i),azmth(j),alb(k),
c    1               total(i,j,k),totn(i,j,k)
666     format('i,j,k,total,totn',3i3,3f7.1,f7.3,1p2e11.3)
c     endif
  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)*ttlp(i)
               gistr=qr(k)*ttrp(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)
c
c        write(33,'(''intsum...eittl,eittr,eitu'',3f9.5)')
c     *             eittl(i,j,k), eittr(i,j,k), eitu(i,j)
c        write(33,135)i,j,k,tnstrl(i),tnstrr(i),qr(k),ttlp(i),
c     1               ttrp(i),eitotz(i,j,k),polz(i,j,k)
            enddo
        enddo
      enddo
  135 format('im,jaz,kalb,tnstrl,tnstrr,qr,ttlp,ttrp,eitotz,polz'/
     1        3i4,7f9.5)

               
      pi=dacos(-1.0d0)
c     print izero, i1,i2 etc as func. of sza and view angle
c
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)/pi
         xi1=eiaz1(i,1)/pi
         xi2=eiaz2(i,1)/pi
c         write(33,133)thnot(imuz),thta(i),xizero,xi1,xi2,tfunc,sb(4),
c     1                rhon,ipol
133      format(2f5.0,1x,6f10.5,1x,i3)
      enddo
c      write(33,134)
134   format(/ )
c
      return
      end