subroutine evalrf(zst1,zst2,ncp1)
c
c**********************************************************************
c
cccc
c     subroutine evalrf
c
c     version aug 22,1977
c
c     purpose
c
c        evalrf is a fortran iv routine which is used to perfform
c         the integration of eq(6.7) to find tnstr (istar/ig)
c
c     method
c
c        using a trapezoidal integration , evalrf computes tenstr which
c        represents the sum of all intensities upto and including the
c        intensity  due to itmax+1 order scattering of surface reflected
c        radiation.
c        the remaining higher order terms are calulated by gepro.
c
c     calling sequence
c
c        call evalrf (zst1,zst2,tnstr,qr,sb)
c
c        variable   type    i/o   description
c        --------   ----    ---   -----------
c
c        zst1(4,202) r*8     i    last three iterations of
c                                 zstar matrix and extrapolated value
c        zst2(4,202) r**     i    similar to zst1
c        tnstr(4,15) r*8     o    integrated intensities from last
c                                 three iterations of z-star matrix
c                                 and extrapolated intensity.
c                                 contains results for each polar
c                                 look angle
c        qr(11)      r*8     o    factor used in computing ig
c        sb(4)       r*8     o    fraction of surface recflected
c                                 radiation reaching top of atmosphere
c
c     external references
c        geopro
c
c     common areas referenced
c
c
c     analysis and programming
c        k.f. klenk,p.m. smith sasc, aug 22 1977
c
c     modifications (date name purpose)
c   
c     last modified 11/19/92 by zia ahmad
c     modified for single iteration
c     last modified 9/22/93 .... zia ahmad
c     purpose: to add the effect of molecular anisotropy 
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.
cccc
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)
c
      include "hedout.inc"
      include "tabbuf.inc"
      include "emm.inc"
      include "eks.inc"
      include "contrl.inc"
      include "consts.inc"
      include "out.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"
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
      do 160 i=1,ncp1
         if (itmax.eq.1) 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
160   continue
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
      if (jprint(10).eq.1) write(23,501)sb,ipol
      if (jprint(2) .eq. 1) write(33,501)sb,ipol
 501  format('evalrf..sb,ipol',1p4e12.4,i5)
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)
      else
            tnstr(j)=0.25d0*(q1/sdp)*(tenr+tenl)/emu(j)
      endif
  120 continue
c
  100 continue
c
c
      do 145 j=1,nalb
      qr(j)=alb(j)/(1.-sb(4)*alb(j))
  145 continue
c
c
c**** store data for archive tape
c**** check if tables data set to be created
c
      rarray(407)=sb(4)
      j=420
      do 200 i=1,ncp1
      k=j+1
      rarray(j)=zst1(4,i)
      rarray(k)=zst2(4,i)
      j=j+2
  200 continue
c**** write data to tape
c
c      write (34) rarray,psaray,iarray
c
c**** check if tables data set to be created
c
c**** load r*4 output buffer for later fwrite
c      out48(1)=float(lambda)
c      out48(2)=alfaef
c      out48(3)=beta
c      out48(4)=code
c      out48(5)=pnot
c      out48(6)=sb(4)
c      j=7
c      do 275 iload=1,15
c      out48(j)=tnstr(iload)
c      j=j+1
c  275 continue
c      out48(22)=float(imuz)
c      out48(23)=float(imu)
c      j=24
c      do 260 i=1,15
c      out48(j)=thta(i)
c      j=j+1
c  260 continue
c      write(35) out48
c
c*****optional print out of z-star array
c
      if (jprint(2) .eq. 1) then
         write(33,1000)
         write(33,2000) (k,(zst1(j,k),j=1,4),(zst2(j,k), j=1,4),
     1    k=1,ncp1)
         write(33,2500)
         write(33,3000) (k,thta(k),tnstr(k),k=1,imu)
         write(33,3500) sb(it)
      endif
c
      return
1000  format(1h1,t50,'debug printout for evalrf',///,25x,
     1 'zst1(i=1,4)',47x,'zst2(i=1,4)')
2000  format(2x,i3,4e13.5,5x,4e13.5)
2500  format(///,2x,'imu',20x,'tnstr(j)')
3000  format(5x,'theta.',5x,'tnstr(k)',//,(1x,i2,3x,f6.2,2x,d12.6))
3500  format(///,' sb=',e15.5)
      end