         subroutine eva1pol(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1,itz)
c***********************************************************************
c
c     subroutine eva1pol
c
c     purpose
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***********************************************************************
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)
c
      include "out.inc"
      include "emm.inc"
      include "eks.inc"
      include "consts.inc"
      include "contrl.inc"
      include "prints.inc"
      include "in.inc"
      include "buff2.inc"
c
c     new statement 
      include "depolt.inc"
c     end of statement
c
      data c316/0.1875d0/
c
c     extrapolate source functions layer by layer
c
      itmaxp=itmax+1
      do 300 i=1,ncp1
         if (itz .le. itmaxp) 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
      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
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)
      else
            eizero(im)=0.125d0*(q1/sdp)*xhold(4)/emu(im)
      endif
c
      e0za(itz,im)=eizero(im)
c
      if(im.eq.1 .and. jprint(4) .eq. 1)then
        write(33,515)itz,im,emu(im),eizero(im),e0za(itz,im)
 515    format('sub eva1pol.itz,im,emu,eizero,e0za',2i5,1x,3f8.5)
      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
      atb=as*bs
      absq=dsqrt(atb)*emuz(jmuz)
      atmu=as*emuz(jmuz)
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)
      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)
      endif
c
      e1za(itz,im,jaz)=eiaz1(im,jaz)
      e2za(itz,im,jaz)=eiaz2(im,jaz)
c
c*****compute total intensity at top of atmosphere
c
      totint(im,jaz)=eizero(im)+eiaz1(im,jaz)+eiaz2(im,jaz)
c
      if (jprint(4) .eq. 1)
     *    write(33,133)itz,im,jaz,eizero(im),eiaz1(im,jaz),
     *             eiaz2(im,jaz),totint(im,jaz)
  133 format('eva1.itz,im,jaz,ez0,ez1,ez2,tot',3i4,4f9.4)
  225 continue
  105 continue
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 
      ggz(itz)=gg
c     for debugging
      if (jprint(4) .eq. 1) write(33,212)gl(4),gr(4),f1,f2,gg
      if (jprint(4) .eq. 1) write(33,312)(e0za(itz,ij),ij=1,imu)
      if (jprint(4) .eq. 1) write(33,313)ggz(itz)
 312  format('e0za',1p6e11.3)
 313  format('ggz',1pe11.3)
 212  format('eva1pol...gl,gr,f1,f2,gg',5f8.4)
c  
      return
      end
