c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/reflexpol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      subroutine reflex(itmax,nek,ncp1)
c
c***********************************************************************
cccc
c     subroutine reflex
c
c     version  aug 22,1977
c
c     purpose
c
c        reflex is a fortran iv routine to calculate the fract. of surf.
c        reflected radiation which is back scattered by the atmosphere.
c        it also computes the ratio of the intensity of light emerging
c        from the top of the atmosphere (which was surface reflected )
c        to the intensity of light reflected by the lambertian surface.
c        this quantity is tenstr=i-star /i-g,see eq 6.7
c
c     method
c
c        the scattering of radiations scattered by the ground is
c        computed.the z-star matrix elements are calculated.
c        by iteration.
c
c     calling sequence
c
c        call reflex (nek,itmax,ncp1,qr)
c
c        variable   type    i/o   description
c        --------   ----    ---   -----------
c
c        nek        i*4      i    fortran data set unit # for
c                                 file containing exponential integrals
c                                 created by expone
c        itmax      i*4      i    number of iteration steps to be used
c        ncp1       i*4      i    # layers
c        qr(10)     r*8      o    factor to be used in computing i-sb-g
c                                 see eqn 6.7
c
c     external references
c        evalrf
c
c     common areas referenced
c
c         eks
c         thkns
c      analysis and programming
c         k.f. klenk, p.m. smith sasc,aug 21 1977
c
c      modifications (date name purpose)
c
c      last modified 11/19/92....zia ahmad
c      modified for single iteration
c      last modified 10/19/94....zia ahmad
c      modified to print results after each itration
c      last modified 03/07/95....dave flittner
c      modified to initialize itp to zero, so that one iteration 
c      may be done.
c      last modified 03/08/95 by dave flittner
c      modified to call eva1pol after each iteration if switch 
c      write_iter_file is TRUE
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 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.
c
ccc
c**********************************************************************
c
      implicit none
      include 'parameter.inc'
c passed
      integer*4 itmax,ncp1
      real*8 nek(202,5,202)
c common
      include 'switches.cmn'
      include 'eks.cmn'
      include 'thkns.cmn'
c local
      integer itp,itmm2,index,i,it,k,ii,iii,j 
      real*8 ekary(202,5),b1(202),b2(202),z1(202),z2(202),zst1(4,202)
      real*8 zst2(4,202),temp1(4,202),temp2(4,202)
      real*8 sum1,sum2
c
c
c        construct zeroth z-star vector
c
       itp=1
       itmm2=itmax-2
       index=1
c
      do 100 i=1,ncp1
         z1(i)=ek4(i)/dtsp(i)
         z2(i)=ek5(i)/dtsp(i)
         zst1(1,i) = z1(i)
         zst2(1,i) = z2(i)
         temp1(1,i)=z1(i)
         temp2(1,i)=z2(i)
100   continue
      if(write_iter_file) call eva2pol(temp1,temp2,ncp1,itp)
      if(itmax.gt.0)then
c
c        do first and higher order iterations
c
        do 140 it=1,itmax
c
c        set b1,b2 equal to new z1,z2
c
           do 150 k=1,ncp1
              b1(k)=z1(k)
              b2(k)=z2(k)
  150      continue
c
c
           do 280 i=1,ncp1
              sum1=ek4(i)/dtsp(i)
              sum2=ek5(i)/dtsp(i)
c
              do 5100 iii = 1,3
                 do 5101 ii = 1,ncp1
                   ekary(ii,iii) = nek(ii,iii,i)
5101             continue
5100          continue
c
              do 190 j=1,ncp1
                 sum1=sum1+ekary(j,1)*b1(j)+ekary(j,2)*b2(j)
                 sum2=sum2+ekary(j,2)*b1(j)+ekary(j,3)*b2(j)
  190         continue
c
c update value of z1 and z2
              z1(i)=sum1
              z2(i)=sum2
c
              temp1(1,i)=z1(i)
              temp2(1,i)=z2(i)
c
              if(it.ge.itmm2)then
c        save zstar matrix of last three iterations
c
                zst1(index,i)=z1(i)
                zst2(index,i)=z2(i)
              endif
280        continue
           if(it.ge.itmm2)index=index+1
c
           itp=itp+1
           if(write_iter_file) call eva2pol(temp1,temp2,ncp1,itp)
c
140     continue
      endif
c
c        call evalrf to calculate istar/ig of eq (6.7)
c        for all polar look angles
c
      itp=itp+1
      if(write_iter_file) call eva2pol(zst1,zst2,ncp1,itp)
      call evalrf(zst1,zst2,ncp1)
c
      return
      end