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,tnstr,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        tnstr(15)  r*8      o    intensity ratio described above for
c                                 each polar back scatter angle.
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.
ccc
c**********************************************************************
c
      implicit integer*4(i-n),real*8(a-h,o-z)
c
      include "eks.inc"
      include "consts.inc"
      include "thkns.inc"
      include "cwrite_iter.inc"
c
      real*8 ekary(202,5),b1(202),b2(202),z1(202),z2(202),
     1       zst1(4,202),zst2(4,202),tnstr(15)
      real*8 nek(202,5,202)
      real*8 temp1(4,202),temp2(4,202)
c
c        construct zeroth z-star vector
c
       itp=0                                           !def
       itmm2=itmax-2
       index=1
       itmaxp=itmax+1
c
      do 100 i=1,ncp1
         z1(i)=ek4(i)/dtsp(i)
         z2(i)=ek5(i)/dtsp(i)
c        if (itmax.eq.1) then
            zst1(1,i) = z1(i)
            zst2(1,i) = z2(i)
            temp1(1,i)=z1(i)
            temp2(1,i)=z2(i)
c        endif
  100 continue
      if(write_iter_file) call eva2pol(temp1,temp2,ncp1,1)   !def
      if (itmax.eq.1) goto 101
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 ii = 1,202
         do 5101 iii = 1,5
            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
      z1(i)=sum1
      z2(i)=sum2
c
      temp1(1,i)=z1(i)
      temp2(1,i)=z2(i)

c        save zstar matrix of last three iterations
c
      if(it.lt.itmm2) go to 280
c
      zst1(index,i)=z1(i)
      zst2(index,i)=z2(i)
  280 continue
c
      itp=it+1
      if((itp .le. itmaxp).or.write_iter_file)then            !def
         call eva2pol(temp1,temp2,ncp1,itp)
      endif
c
      if(it.ge.itmm2) index=index+1
  140 continue
101   continue
c
c        call evalrf to calculate istar/ig of eq (6.7)
c        for all polar look angles
c
      itpp=itp+1
      if(write_iter_file) call eva2pol(zst1,zst2,ncp1,itpp)   !def
      call evalrf(zst1,zst2,ncp1)
c
      return
      end