c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/evalrf.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      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          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.
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. Feb. 10, 2000...def
c purpose: add switch to turn-off call to geopro and also prevent call to geopro
c if itmax less than 3.
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
cccc
c***********************************************************************
c
c
      implicit integer*4(i-n),real*8(a-h,o-z)
      include 'parameter.inc'
      real*8  zst1(4,202),zst2(4,202)
      include 'emm.cmn'
      include 'eks.cmn'
      include 'contrl.cmn'
      include 'out.cmn'
      include 'switches.cmn'
c
c     -- common block in added on 11/19/92
      include 'in.cmn'
c     common block depolt added on 9/22/93
      include 'depolt.cmn'
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
      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 160 i=1,ncp1
           call geopro(zst1(1,i))
           call geopro(zst2(1,i))
160     continue
      endif
c
c
c        find sb
      sum2=0.d0
      sbp=0.0d0
c
      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
      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(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(2).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)
            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
  120 continue
c
  100 continue
c
c
      do 145 j=1,nalb
         qr(j)=alb(j)/(1.-sb*alb(j))
  145 continue
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
      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