c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/eva2pol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      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     last modified 8/30/97...dave flittner
c     purpose: compute downward diffuse from surface reflectance
c     contribution to the actinic flux, sbp
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 the # of multiple scattering events.  So,
c modify to allow itmax=0.  Also change sb and sbp to single elements.
c
c
c***********************************************************************
c
c
      implicit none
      include 'parameter.inc'
c passed
      real*8  zst1(4,202),zst2(4,202)
      integer ncp1,itz
c common
      include 'eks.cmn'
      include 'emm.cmn'
      include 'contrl.cmn'
      include 'out.cmn'
      include 'switches.cmn'
c     -- common block in added on 11/19/92
      include 'in.cmn'
c     common block depolt added on 9/22/93
      include 'depolt.cmn'
      include 'buff1.cmn'
c local
      real*8 c316,sum2,suma,sumb,tenr,tenl,tensum
      integer itmaxp,i,j,index
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        compute gometrical progression of zstar layer by layer
c
      itmaxp=itmax+1
      if (itz .le. itmaxp) then
        do i=1,ncp1
            zst1(4,i) = zst1(1,i)
            zst2(4,i) = zst2(1,i)
        enddo
      else
        if ((itmax.lt.3).or.lnoextrap) then
          if(itmax.lt.1)then
            index=1
          else
            index=min(itmax,3)
          endif
          do i=1,ncp1
             zst1(4,i) = zst1(index,i)
             zst2(4,i) = zst2(index,i)
          enddo
        else
          do i=1,ncp1
             call geopro(zst1(1,i))
             call geopro(zst2(1,i))
          enddo
        endif
      endif
c
c        find sb
c
      sum2=0.d0
      sbp=0.0d0
      do 110 i=1,ncp1
         sum2=sum2+ek4(i)*zst1(4,i)+ek5(i)*zst2(4,i)
         sbp=sbp+ek8(i)*zst1(4,i)+ek9(i)*zst2(4,i)
  110 continue
c
c     new statements added here
      if(ipol.eq.0)then
        sb=0.375d0*sum2
        sbp=0.375d0*sbp
      else
        sb=q12s*sum2
        sbp=q12s*sbp
      endif
c
      sbz(itz)=sb
      if (jprint(4) .eq. 1) write(33,501)itz,itmax,ipol,sb
 501  format('eva2pol...itz,itmax,ipol,sb',3i5,1pe12.4)
c        find tnstr(4)
c
      do 120 j=1,imu
         suma=0.
         sumb=0.
c
         do 130 i=1,ncp1
            suma=suma+zst1(4,i)*extmu(i,j)
            sumb=sumb+zst2(4,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
         if(jprint(4).eq.1)write(33,565)suma,sumb,ematx(1,j),
     & ematx(2,j),ematx(3,j),emu(j),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)
         else
           tnstr(j)=0.25d0*(q1/sdp)*(tenr+tenl)/emu(j)
         endif
         tnstrz(itz,j)=tnstr(j) 
  120 continue
c
      return
      end