subroutine eva2pol(zst1,zst2,ncp1,itz)
c
c**********************************************************************
c
c     subroutine eva2pol
c
c     this routine is adopted from subroutine evalrf to extract sb,tnstr 
c     after each iteration......zia ahmad 10/21/94
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  tenstr(4,15),zst1(4,202),zst2(4,202)
      include "emm.inc"
      include "eks.inc"
      include "contrl.inc"
      include "out.inc"
      include "consts.inc"
      include "prints.inc"
c     -- common block in added on 11/19/92
      include "in.inc"
c     common block depolt added on 9/22/93
      include "depolt.inc"
      include "buff1.inc"
c
      data c316/0.1875d0/
c
c     write(23,500)rhon,gama,q,q1,q2,delp,sdp,q12s,ipol
 500  format('rhon,gama,q,q1,q2,delp,sdp,q12s,ipol'/
     1        1p6e12.4/1p2e12.4,i5)
c
c        compute gometrical progression of zstar layer by layer
c
      itmaxp=itmax+1
      do i=1,ncp1
         if (itz .le. itmaxp) then
            zst1(4,i) = zst1(1,i)
            zst2(4,i) = zst2(1,i)
         else
            call geopro(zst1(1,i))
            call geopro(zst2(1,i))
         endif

      enddo
c
      it=4
c
c        find sb(it)
      sum2=0.d0
c
      sum2l=0.0d0
      sum2r=0.0d0
      do 110 i=1,ncp1
      sum2=sum2+ek4(i)*zst1(it,i)+ek5(i)*zst2(it,i)
  110 continue
c
c     new statements added here
      if(ipol.eq.0)then
            sb(it)=0.375d0*sum2
      else
            sb(it)=q12s*sum2
      endif
c
      sbz(itz)=sb(it)
      if (jprint(4) .eq. 1) write(33,501)itz,itmax,ipol,sb(it)
 501  format('eva2pol...itz,itmax,ipol,sb',3i5,1pe12.4)
c        find tnstr(it)
c
      do 120 j=1,imu
      suma=0.
      sumb=0.
c
      do 130 i=1,ncp1
      suma=suma+zst1(it,i)*extmu(i,j)
      sumb=sumb+zst2(it,i)*extmu(i,j)
  130 continue
c
c        multiply by m-matrix see eq(6.7)
c
      tenr=ematx(1,j)*suma+ematx(2,j)*sumb
      tenl=ematx(3,j)*suma
c
      tensum=tenl+tenr
c     write(33,565)suma,sumb,ematx(1,j),ematx(2,j),ematx(3,j),emu(j),
c    1             tenl,tenr,tensum
 565  format('evalrf',9f8.5)
c
c        compute istar/ig of eq (6.7)
c
c     new statements added here
      if(ipol.eq.0)then
            tnstr(j)=c316*(tenr+tenl)/emu(j)
            tnstrl(j)=c316*tenl/emu(j)
            tnstrr(j)=c316*tenr/emu(j)
      else
            tnstr(j)=0.25d0*(q1/sdp)*(tenr+tenl)/emu(j)
            tnstrl(j)=0.25d0*(q1/sdp)*tenl/emu(j)
            tnstrr(j)=0.25d0*(q1/sdp)*tenr/emu(j)
      endif
      tnstrz(itz,j)=tnstr(j) 
  120 continue
c
      return
      end