subroutine evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
c***********************************************************************
c
c     subroutine evalit
c
c     purpose
c        evalit is used to perform the integration of the reduced
c        source functions over the total optical thickness to each
c        level in the model atmosphere.
c        evalit is also used to transfer solar angle dependent data
c        to the archive and tables tapes.
c
c    method
c       the iterated reduced source functions as calculated by itrate
c       are integrated over the total cumulative optical thickness.the
c       resultant integrated quantities (ztop,ztp1,ztp2,gg) are extrapol
c       extrapolated using a geometric progression. finally certain
c       matrix multiplications are carried out to transform the
c       integrated quantities into intensity terms.
c
c       calling sequence
c        call evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
c
c        variable     type    i/o    description
c        --------     ----    ---    -----------
c
c        zma1(4,202)  r*8      i     first element of z matrix
c                                    for each source function
c                                    level for final three iters.
c                                    fourth row of this matrix is
c                                    reserved for extrapolated values.
c        zma2,3,4     r*8     i/o    second,third,fourth elements
c                                    of z matrix.
c        ej1(4,202)     r*8      o   source functions for cos(phi)
c                                    dependent terms for last three
c                                    iterations of itrate.
c        ej2(4,202)   r*8      i/o   source functions for cos(2*phi)
c                                    dependent terms.
c        jmuz         i*4       i    current index counter for zenith
c                                    angle.
c        ncp1         i*4       i    number of levels in model atmos.
c     common /out/ contains intensity terms computed by
c     evalrf,evalit and intsum routines.
c
c     common /out/
c
c       variables    type    i/o  description
c       ---------    ----    ---  -----------
c
c       sb(4)        r*8      i   fractional atmos. back scatter
c                                 factor.sb(4) is extrapolated value.
c       tnstr(15)     r*8     i   extrapolated values of istar/ig for
c                                 each polar look angle.
c       qr(10)       r*8      i   combined reflectivity -sb term for
c                                 all albedo values.
c       eizero(15)   r*8      o   azimuth independent intensity at each
c                                 polar look angle.
c       eiaz1(15,8)  r*8      o   cos(phi) azimuth dependent term for
c                                 each polar look angle and scan plane
c                                 angle.
c       eiaz2(15,8)  r*8      o   cos(2*phi) dependent term for each
c                                 polar look angle and scan plane angle.
c       gg           r*8      o   integrated downward diffuse flux at
c                                 ground for current zenith angle.
c       fs           r*8      i   direct solar flux at ground for curren
c                                 current solar zenith angle.
c       totint(15,8) r*8      o   total intensity at top of atmosphere
c                                 for each polar and azimuth angle.
c
c       analysis and programming   k. f. klenk, p. m. smith sasc aug. 77
c
c       modifications (date name purpose)
c
c       last modified 11/18/92 by zia ahmad
c       modified for single iteration
c
c       last modified 9/22/93 ... zia ahmad
c       purpose: to add the effect of molecular depolarization
c
c       last modified 10/12/94... zia ahmad
c       purpose: modified to compute and print polarization of the 
c       diffused radiation.  
c
c       last modified 03/08/95... dave flittner
c       purpose: store Z1func=I1/(-3/8*muo*sqrt(1-muo^2)*sqrt(1-mu^2)) and
c       Z2func=I2/(3/32*(1-muo^2)*(1-mu^2)/mu) so can be used as input into
c       interpolation tables.
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.
c***********************************************************************
c
      implicit integer*4 (i-n), real*8(a-h,o-z)
      real*8 ztop(4,4),a(4),b(4),c(4),d(4),ztp1(4),ztp2(4),top(4),
     &       dum(8),gr(4),gl(4),
     &       zma1(4,202),zma2(4,202),zma3(4,202),zma4(4,202),
     &       ej1(4,202),ej2(4,202),hold(4),xhold(4),eiold(15)
      real*8 eiaz1l(15,8),eiaz1r(15,8),eiaz1u(15,8),
     1       eiaz2l(15,8),eiaz2r(15,8),eiaz2u(15,8),pol(15,8),
     2       eizrol(15),eizror(15)
c
      include "tabbuf.inc"
      include "hedout.inc"
      include "out.inc"
      include "emm.inc"
      include "eks.inc"
      include "consts.inc"
      include "contrl.inc"
      include "prints.inc"
      include "in.inc"
      include "buff3.inc"
      include "czfunc.inc"
c     new statement 
      include "depolt.inc"
c     end of statement
c
      dimension raray(906)
      equivalence (rarray(1),raray(1))
c
      data c316/0.1875d0/
c     extrapolate source functions layer by layer
c
      do 300 i=1,ncp1
         if (itmax.eq.1) then
            zma1(4,i) = zma1(1,i)
            zma2(4,i) = zma2(1,i)
            zma3(4,i) = zma3(1,i)
            zma4(4,i) = zma4(1,i)
            ej1(4,i)  = ej1(1,i)
            ej2(4,i)  = ej2(1,i)
         else
            call geopro(zma1(1,i))
            call geopro(zma2(1,i))
            call geopro(zma3(1,i))
            call geopro(zma4(1,i))
            call geopro(ej1(1,i))
            call geopro(ej2(1,i))
         endif
  300 continue
c
c*****put adjoint matrix for current solar zenith angle into b
c
      b(1)=admatx(1,jmuz)
      b(2)=admatx(2,jmuz)
      b(3)=admatx(3,jmuz)
      b(4)=0.0d0
c
      if(jprint(1).eq.1) write(33,1000)
      do 105 im=1,imu
      it=4
      gr(it)=0.0d0
      gl(it)=0.0d0
      ztp1(it)=0.0d0
      ztp2(it)=0.0d0
      do 107 jd=1,4
  107 ztop(it,jd)=0.0d0
c
c*****do integration over optical thickness
c
      do 110 i=1,ncp1
      tmu=extmu(i,im)
      ztop(it,1)=ztop(it,1)+zma1(it,i)*tmu
      ztop(it,2)=ztop(it,2)+zma2(it,i)*tmu
      ztop(it,3)=ztop(it,3)+zma3(it,i)*tmu
      ztop(it,4)=ztop(it,4)+zma4(it,i)*tmu
      ztp1(it)=ztp1(it)+ej1(it,i)*tmu
      ztp2(it)=ztp2(it)+ej2(it,i)*tmu
      gr(it)=gr(it)+ek4(i)*zma1(it,i)+ek5(i)*zma3(it,i)
      gl(it)=gl(it)+ek4(i)*zma2(it,i)+ek5(i)*zma4(it,i)
110   continue
      if(ipol.eq.0)then                                !def
        Z1func(im)=ztp1(4)                             !def
        Z2func(im)=ztp2(4)                             !def
      else                                             !def
        Z1func(im)=q2*ztp1(4)                          !def
        Z2func(im)=q2*ztp2(4)                          !def
      endif                                            !def
  100 continue
c
c     write(33,132)ipol, ztp1(it),ztp2(it)
 132  format('debug evalit..ipol,ztp1,ztp2',i2,1x,1p2e12.4)
c
c*****calculate the azimuth -independent intensity eizero at current
c*****value of the polar angle. see eqs(3.7) and (4.16)
c
c*****multiply matrix top(1-4) times the adjoint matrix admatx(=b). cal
c*****call the product c
c
c
c*****multiply c by ematx
c
      a(1)=ematx(1,im)
      a(2)=ematx(2,im)
      a(3)=ematx(3,im)
      a(4)=0.0d0
c
c
c
      do 120 i=1,4
c
c**** store extrapolated z-matrix in hold(4)
c
      hold(i)=ztop(it,i)
  120 continue
      call matmul(hold,b,c)
      call matmul(a,c,d)
      xhold(4)=d(1)+d(2)+d(3)+d(4)
  130 continue
c     new statements added here
      if(ipol.eq.0) then
            eizero(im)=0.09375d0*xhold(4)/emu(im)
            eizrol(im)=0.09375d0*(d(1)+d(2))/emu(im)
            eizror(im)=0.09375d0*(d(3)+d(4))/emu(im)
      else
            eizero(im)=0.125d0*(q1/sdp)*xhold(4)/emu(im)
            eizrol(im)=0.125d0*(q1/sdp)*(d(1)+d(2))/emu(im)
            eizror(im)=0.125d0*(q1/sdp)*(d(3)+d(4))/emu(im)
      endif
c
c     if(im.eq.1)then
        if (jprint(1) .eq. 1) write(33,514)im,emu(im),
     *                           d(1),d(2),d(3),d(4)
514     format('im,emu,d1,d2,d3,d4',i3,1x,5f8.4)
        if (jprint(1) .eq. 1) write(33,515)im,emu(im),
     *                      eizero(im),eizrol(im),eizror(im)
 515    format('evalitpol..im,emu,ei0,eil,eir',1x,i2,1x,4f8.5)
c     endif
c
c*****now calculate azimuth dependent intensities eiaz1 and eiaz2 at all
c*****azimuth angles (1 to iazmth) for current value of polar angle.
c*****           see eqs.(3.7) and (4.8)
c
      if(ipol.eq.0)then
         as=ematx(1,im)-1.0d0
         bs=admatx(1,jmuz)-1.0d0
      else
         as=emu(im)*emu(im)-1.0d0
         bs=emuz(jmuz)*emuz(jmuz)-1.0d0
      endif
c
      cs=emu(im)*emu(im)
      atb=as*bs
      absq=dsqrt(atb)*emuz(jmuz)
      atmu=as*emuz(jmuz)
c
c     write(33,134) ipol,im,jmuz,as,bs,atb,absq,atmu
 134  format('ipol,as,bs,atb,absq,atmu'/3i2,1x,1p5e12.4)
c
      do 225 jaz=1,iazmth
c     new statements added here
      if(ipol.eq.0)then
            eiaz1(im,jaz)=-0.375d0*absq*caz(jaz)*ztp1(4)
            eiaz2(im,jaz)=0.09375d0*(atb*caz2(jaz))*ztp2(4)/emu(im)
c     new code for polarization 
            eiaz1l(im,jaz)=eiaz1(im,jaz)
            eiaz1r(im,jaz)=0.0d0
            eiaz1u(im,jaz)=0.375d0*absq*dsqrt(1.0d0-caz(jaz)**2)*
     1                    ztp1(4)/emu(im)
            eiaz2l(im,jaz)=0.09375d0*cs*bs*caz2(jaz)*ztp2(4)/emu(im)
            eiaz2r(im,jaz)=-0.09375d0*bs*caz2(jaz)*ztp2(4)/emu(im)
            eiaz2u(im,jaz)=-0.09375d0*bs*
     1                     dsqrt(1.0d0-caz2(jaz)**2)*ztp2(4)
      else
            eiaz1(im,jaz)=-0.375d0*q2*absq*caz(jaz)*ztp1(4)
            eiaz2(im,jaz)=0.09375d0*q2*(atb*caz2(jaz))*ztp2(4)/emu(im)
c     new code for polarization 
            eiaz1l(im,jaz)=eiaz1(im,jaz)
            eiaz1r(im,jaz)=0.0d0
            eiaz1u(im,jaz)=0.375d0*q2*absq*dsqrt(1.0d0-caz(jaz)**2)*
     1                    ztp1(4)/emu(im)
            eiaz2l(im,jaz)=0.09375d0*q2*cs*bs*caz2(jaz)*ztp2(4)/emu(im)
            eiaz2r(im,jaz)=-0.09375d0*q2*bs*caz2(jaz)*ztp2(4)/emu(im)
            eiaz2u(im,jaz)=-0.09375d0*q2*bs*
     1                     dsqrt(1.0d0-caz2(jaz)**2)*ztp2(4)
      endif
c
c*****compute total intensity at top of atmosphere
c
      totint(im,jaz)=eizero(im)+eiaz1(im,jaz)+eiaz2(im,jaz)
c
c     compute degree of polarization
      eitl(im,jaz)=eizrol(im)+eiaz1l(im,jaz)+eiaz2l(im,jaz)
      eitr(im,jaz)=eizror(im)+eiaz1r(im,jaz)+eiaz2r(im,jaz)
      eitu(im,jaz)=eiaz1u(im,jaz)+eiaz2u(im,jaz)
      pol(im,jaz)=dsqrt((eitl(im,jaz)-eitr(im,jaz))**2
     1            +eitu(im,jaz)**2)/totint(im,jaz)
c
      if (jprint(1) .eq. 1) write(33,133)im,jaz,eizero(im),
     *      eiaz1(im,jaz),eiaz2(im,jaz),totint(im,jaz),pol(im,jaz)
  133 format('evalitpol...im,jaz,ei0,ei1,ei2,totint,pol'/
     1        2i4,1p5e11.3)
      if (jprint(1) .eq. 1) write(33,136)eitl(im,jaz),eitr(im,jaz),
     *        eitu(im,jaz)
 136  format(8x,1p3e11.3)
  225 continue
      if(jprint(1).eq.1) write(33,1100)im,thta(im),
     1 (ztop(4,j),j=1,4),ztp1(4),ztp2(4)
  105 continue
      if(jprint(1).eq.1) write(33,2000) gl(4),gr(4)
c
c*****compute the downward diffuse intensity gg
c
      if(ipol .eq. 0)then
c     gg=c316*(gr(4)*(admatx(1,jmuz)+1.0d0)+gl(4)*admatx(3,jmuz))
      f1=c316*(admatx(1,jmuz)+1.0d0)
      f2=c316*admatx(3,jmuz)
      gg=f1*gr(4)+f2*gl(4)
      else
      f1=0.25d0*q1*(admatx(1,jmuz)+1.0d0+q)
      f2=0.25d0*q1*admatx(3,jmuz)
      gg=(f1*gr(4)+f2*gl(4))/sdp
      endif
c 
c     for debugging
c     write(33,212)gl(4),gr(4),f1,f2,gg,q
 212  format('gl,gr,f1,f2,gg,q',6f8.4)
c     write(33,213)(admatx(iz,jmuz),iz=1,3)
 213  format(1x,'admatx 1-3',3f8.4)
c  
c**** store data for archive tape
      raray(1)=fs
      raray(2)=gg
      k1=3
      do 400 i=1,ncp1
      k5=810+i
      k6=1012+i
      raray(k1)=zma1(4,i)
      raray(k1+1)=zma2(4,i)
      raray(k1+2)=zma3(4,i)
      raray(k1+3)=zma4(4,i)
      raray(k5)=ej1(4,i)
      raray(k6)=ej2(4,i)
      k1=k1+4
  400 continue
C      write (34) raray,iarray
c      out48(1)=thnot(jmuz)
c      j=2
c      do 420 i=1,15
c      out48(j)=eizero(i)
c      out48(j+15)=eiaz1(i,1)
c      out48(j+30)=eiaz2(i,1)
c      j=j+1
c  420 continue
c      out48(47)=fs
c      out48(48)=gg
c      write(35)out48
  500 continue
c**** call intsum to compute total intensity
      call intsum
      if (jprint(10).eq.1) call sumry(jmuz)
      return
 1000 format(1h1,t5,'debug print out from evalit',10x,
     1       'source functions(extrapolated) integrated over',
     2       1x,'optical thickness',///,
     2       1x,'imu',2x,'theta',5x,'ztop1',9x,'ztop2',9x,'ztop3',9x,
     3       'ztop4',9x,' ztp1',9x,' ztp2',//)
 1100 format(1x,i3,2x,f5.2,2x,d12.6,2x,d12.6,2x,d12.6,2x,d12.6,
     1       2x,d12.6,2x,d12.6)
 2000 format(///,1x,'gl= ',d12.6,/,1x,'gr= ',d12.6)
      end