c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/evalitpol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
         subroutine evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
c***********************************************************************
c
c     subroutine evalit
c
c     purpose
c        evalit is used to perform the integration of the reduced
c        source functions over the total optical thickness to each
c        level in the model atmosphere.
c        evalit is also used to transfer solar angle dependent data
c        to the archive and tables tapes.
c
c    method
c       the iterated reduced source functions as calculated by itrate
c       are integrated over the total cumulative optical thickness.the
c       resultant integrated quantities (ztop,ztp1,ztp2,gg) are
c       extrapolated using a geometric progression. finally certain
c       matrix multiplications are carried out to transform the
c       integrated quantities into intensity terms.
c
c       calling sequence
c        call evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
c
c        variable     type    i/o    description
c        --------     ----    ---    -----------
c
c        zma1(4,202)  r*8      i     first element of z matrix
c                                    for each source function
c                                    level for final three iters.
c                                    fourth row of this matrix is
c                                    reserved for extrapolated values.
c        zma2,3,4     r*8     i/o    second,third,fourth elements
c                                    of z matrix.
c        ej1(4,202)     r*8      o   source functions for cos(phi)
c                                    dependent terms for last three
c                                    iterations of itrate.
c        ej2(4,202)   r*8      i/o   source functions for cos(2*phi)
c                                    dependent terms.
c        jmuz         i*4       i    current index counter for zenith
c                                    angle.
c        ncp1         i*4       i    number of levels in model atmos.
c     common /out/ contains intensity terms computed by
c     evalrf,evalit and intsum routines.
c
c     common /out/
c
c       variables    type    i/o  description
c       ---------    ----    ---  -----------
c
c       sb           r*8      i   fractional atmos. back scatter
c                                 factor.sb is extrapolated value.
c       tnstr(15)     r*8     i   extrapolated values of istar/ig for
c                                 each polar look angle.
c       qr(10)       r*8      i   combined reflectivity -sb term for
c                                 all albedo values.
c       eizero(15)   r*8      o   azimuth independent intensity at each
c                                 polar look angle.
c       eiaz1(15,8)  r*8      o   cos(phi) azimuth dependent term for
c                                 each polar look angle and scan plane
c                                 angle.
c       eiaz2(15,8)  r*8      o   cos(2*phi) dependent term for each
c                                 polar look angle and scan plane angle.
c       gg           r*8      o   integrated downward diffuse flux at
c                                 ground for current zenith angle.
c       fs           r*8      i   direct solar flux at ground for curren
c                                 current solar zenith angle.
c       totint(15,8) r*8      o   total intensity at top of atmosphere
c                                 for each polar and azimuth angle.
c
c       analysis and programming   k. f. klenk, p. m. smith sasc aug. 77
c
c       modifications (date name purpose)
c
c       last modified 11/18/92 by zia ahmad
c       modified for single iteration
c
c       last modified 9/22/93 ... zia ahmad
c       purpose: to add the effect of molecular depolarization
c
c       last modified 10/12/94... zia ahmad
c       purpose: modified to compute and print polarization of the 
c       diffused radiation.  
c
c       last modified 03/08/95... dave flittner
c       purpose: store Z1func=I1/(-3/8*muo*sqrt(1-muo^2)*sqrt(1-mu^2)) and
c       Z2func=I2/(3/32*(1-muo^2)*(1-mu^2)/mu) so can be used as input into
c       interpolation tables.
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 05/01/96...dave flittner
c       purpose: delete lines using matmul incorrectly to multiply ztop
c       by b to give c, and then multiplying a by c to get d.
c
c       last modified 09/03/96...dave flittner
c       purpose: allow computation of the transmitted radiance at the surface.
c       This is activated by the logical switch ldown set to TRUE.
c     
c     last modified: 09/20/96...dave flittner
c     purpose: Add spherical single scatter to radiance if ipsudo=2.
c     Will also rename some local variables as follows:
c     zma1top(*)=ztop(*,1)
c     zma2top(*)=ztop(*,2)
c     zma3top(*)=ztop(*,3)
c     zma4top(*)=ztop(*,4)
c     ej1top(*)=ztp1(*)
c     ej2top(*)=ztp2(*)
c
c     last modified: 8/30/97...dave flittner
c     purpose: Compute the variables needed for calculation of the
c     actinic flux at the surface, ggp (this is the contribution from 
c     diffuse light with zero surface reflectance).
c
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn.
c 
c last mod: May 10, 1999...def
c purpose: change the name to matmlt to avoid conflict with linex f77 compiler.
c
c last mod: May 12, 1999...def
c purpose: change .ne. to .neqv. to compile with linex f77 compiler.
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***********************************************************************
c
      implicit none
      include 'parameter.inc'
c passed in
      integer jmuz,ncp1
      real*8 zma1(4,202),zma2(4,202),zma3(4,202),zma4(4,202)
      real*8 ej1(4,202),ej2(4,202)
c from commons
      include 'czsing.cmn'
      include 'switches.cmn'
      include 'czfunc.cmn'
      include 'uptoa.cmn'
      include 'energy.cmn'
      include 'eks.cmn'
      include 'out.cmn'
      include 'emm.cmn'
      include 'contrl.cmn'
      include 'in.cmn'
      include 'depolt.cmn'
      include 'buff3.cmn'
c local
      integer i,im,k,index
      real*8 zma1top(max_az),zma2top,zma3top,zma4top(max_az)
      real*8 ej1top(max_az),ej2top(max_az),a(4),b(4),c(4),d(4)
      real*8 gr,gl,hold(4),xhold(4),gruptoa,gluptoa,tmu,grp,glp
      real*8 eiaz1l(max_scan,max_az),eiaz1r(max_scan,max_az)
      real*8 eiaz1u(max_scan,max_az),eiaz2l(max_scan,max_az)
      real*8 eiaz2r(max_scan,max_az),eiaz2u(max_scan,max_az)
      real*8 pol(max_scan,max_az),eizrol(max_scan,max_az)
      real*8 eizror(max_scan,max_az),as,bs,atmu,cs,atb,c316,absq,f1,f2
      logical lprtflux
c
      data c316/0.1875d0/
c
      lprtflux=.FALSE.
c
c     extrapolate source functions 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 300 i=1,ncp1
            zma1(4,i) = zma1(index,i)
            zma2(4,i) = zma2(index,i)
            zma3(4,i) = zma3(index,i)
            zma4(4,i) = zma4(index,i)
            ej1(4,i)  = ej1(index,i)
            ej2(4,i)  = ej2(index,i)
  300    continue
      else
         do i=1,ncp1
            call geopro(zma1(1,i))
            call geopro(zma2(1,i))
            call geopro(zma3(1,i))
            call geopro(zma4(1,i))
            call geopro(ej1(1,i))
            call geopro(ej2(1,i))
         enddo
      endif
c
c*****put adjoint matrix for current solar zenith angle into b
c
      b(1)=admatx(1,jmuz)
      b(2)=admatx(2,jmuz)
      b(3)=admatx(3,jmuz)
      b(4)=0.0d0
c
      if(jprint(1).eq.1) write(33,1000)
c
c compute the flux integration for downward at surface and upward at top
c
      gr=0.0d0
      gl=0.0d0
      grp=0.0d0
      glp=0.0d0
      gruptoa=0.0d0
      gluptoa=0.0d0
      do i=1,ncp1
         gr=gr+ek4(i)*zma1(4,i)+ek5(i)*zma3(4,i)
         gl=gl+ek4(i)*zma2(4,i)+ek5(i)*zma4(4,i)
         grp=grp+ek8(i)*zma1(4,i)+ek9(i)*zma3(4,i)
         glp=glp+ek8(i)*zma2(4,i)+ek9(i)*zma4(4,i)
         gruptoa=gruptoa+ek6(i)*zma1(4,i)+ek7(i)*zma3(4,i)
         gluptoa=gluptoa+ek6(i)*zma2(4,i)+ek7(i)*zma4(4,i)
      enddo
c
c now for the reflected radiance
c
c if ipsudo=2, then eizero will actually depend upon the azimuth angle
c due to the single scattering along the line of sight.  So, will need
c to loop iazmth times to do the integration. Otherwise, only integrate
c once and then duplicate for the other azimuth directions since
c there really is no azimuth dependance in that case.
c
      do 105 im=1,imu
         zma2top=0.0d0
         zma3top=0.0d0
         do k=1,iazmth
            zma1top(k)=0.0d0
            zma4top(k)=0.0d0
            ej1top(k)=0.0d0
            ej2top(k)=0.0d0
         enddo
c
c*****do integration over optical thickness
c
         do 110 i=1,ncp1
            tmu=extmu(i,im)
            zma2top=zma2top+zma2(4,i)*tmu
            zma3top=zma3top+zma3(4,i)*tmu
            if(ipsudo.lt.2)then
              zma1top(1)=zma1top(1)+zma1(4,i)*tmu
              zma4top(1)=zma4top(1)+zma4(4,i)*tmu
              ej1top(1)=ej1top(1)+ej1(4,i)*tmu
              ej2top(1)=ej2top(1)+ej2(4,i)*tmu
            else
              do k=1,iazmth
                 zma1top(k)=zma1top(k)+(zma1(4,i)+zsp(i,k,im))*tmu
                 zma4top(k)=zma4top(k)+(zma4(4,i)+zsp(i,k,im))*tmu
                 ej1top(k)=ej1top(k)+(ej1(4,i)+zsp(i,k,im))*tmu
                 ej2top(k)=ej2top(k)+(ej2(4,i)+zsp(i,k,im))*tmu
              enddo
            endif
110      continue
         if(ipsudo.lt.2)then
           do k=2,iazmth
              ej1top(k)=ej1top(1)
              ej2top(k)=ej2top(1)
           enddo
         endif
         do k=1,iazmth
            if(ipol.eq.0)then
              Z1func(im,k)=ej1top(k)
              Z2func(im,k)=ej2top(k)
            else
              Z1func(im,k)=q2*ej1top(k)
              Z2func(im,k)=q2*ej2top(k)
            endif
         enddo
c
         if(jprint(4).eq.1)write(33,132)ipol,ej1top(1),ej2top(1)
 132  format('debug evalit..ipol,ej1top,ej1top',i2,1x,1p2e12.4)
c
c*****calculate the azimuth -independent intensity eizero at current
c*****value of the polar angle. see eqs(3.7) and (4.16)
c
c load ematx for the current scan angle
c
         a(1)=ematx(1,im)
         a(2)=ematx(2,im)
         a(3)=ematx(3,im)
         a(4)=0.0d0
c
c loop over the azimuth angles
c
         do k=1,iazmth
c
c**** store extrapolated z-matrix in hold(4)
c
            if((k.eq.1).or.(ipsudo.eq.2))then
              hold(1)=zma1top(k)
              hold(2)=zma2top
              hold(3)=zma3top
              hold(4)=zma4top(k)
c*****multiply matrix z(1-4) times the adjoint matrix admatx(=b). cal
c*****call the product c
c
              call matmlt(hold,b,c)
c
c*****multiply c by ematx
c
              call matmlt(a,c,d)
              xhold(4)=d(1)+d(2)+d(3)+d(4)
            endif
c     new statements added here
            if(ipol.eq.0) then
               eizero(im,k)=0.09375d0*xhold(4)/emu(im)
               eizrol(im,k)=0.09375d0*(d(1)+d(2))/emu(im)
               eizror(im,k)=0.09375d0*(d(3)+d(4))/emu(im)
            else
               eizero(im,k)=0.125d0*(q1/sdp)*xhold(4)/emu(im)
               eizrol(im,k)=0.125d0*(q1/sdp)*(d(1)+d(2))/emu(im)
               eizror(im,k)=0.125d0*(q1/sdp)*(d(3)+d(4))/emu(im)
            endif
c
            if (jprint(1) .eq. 1)then
              write(33,514)im,emu(im),d(1),d(2),d(3),d(4)
514     format('im,emu,d1,d2,d3,d4',i3,1x,5f8.4)
              write(33,515)im,emu(im),eizero(im,k),eizrol(im,k),
     &                     eizror(im,k)
 515    format('evalitpol..im,emu,ei0,eil,eir',1x,i2,1x,4f8.5)
            endif
         enddo
c
c*****now calculate azimuth dependent intensities eiaz1 and eiaz2 at all
c*****azimuth angles (1 to iazmth) for current value of polar angle.
c*****           see eqs.(3.7) and (4.8)
c
         if(ipol.eq.0)then
            as=ematx(1,im)-1.0d0
            bs=admatx(1,jmuz)-1.0d0
         else
            as=emu(im)*emu(im)-1.0d0
            bs=emuz(jmuz)*emuz(jmuz)-1.0d0
         endif
c
         cs=emu(im)*emu(im)
         atb=as*bs
         absq=dsqrt(atb)*emuz(jmuz)
         atmu=as*emuz(jmuz)
c
         if(jprint(1).eq.1)write(33,134)ipol,im,jmuz,as,bs,atb,absq,atmu
 134  format('ipol,as,bs,atb,absq,atmu'/3i2,1x,1p5e12.4)
c
         do 225 k=1,iazmth
c     new statements added here
            if(ipol.eq.0)then
              eiaz1(im,k)=-0.375d0*absq*caz(k)*ej1top(k)
              eiaz2(im,k)=0.09375d0*(atb*caz2(k))*ej2top(k)/emu(im)
c     new code for polarization 
              eiaz1l(im,k)=eiaz1(im,k)
              eiaz1r(im,k)=0.0d0
              eiaz1u(im,k)=0.375d0*absq*saz(k)*ej1top(k)/emu(im)
              eiaz2l(im,k)=0.09375d0*cs*bs*caz2(k)*ej2top(k)/emu(im)
              eiaz2r(im,k)=-0.09375d0*bs*caz2(k)*ej2top(k)/emu(im)
              eiaz2u(im,k)=-0.09375d0*bs*saz2(k)*ej2top(k)
            else
              eiaz1(im,k)=-0.375d0*q2*absq*caz(k)*ej1top(k)
              eiaz2(im,k)=0.09375d0*q2*(atb*caz2(k))*ej2top(k)/emu(im)
c     new code for polarization 
              eiaz1l(im,k)=eiaz1(im,k)
              eiaz1r(im,k)=0.0d0
              eiaz1u(im,k)=0.375d0*q2*absq*saz(k)*ej1top(k)/emu(im)
              eiaz2l(im,k)=.09375d0*q2*cs*bs*caz2(k)*ej2top(k)/emu(im)
              eiaz2r(im,k)=-0.09375d0*q2*bs*caz2(k)*ej2top(k)/emu(im)
              eiaz2u(im,k)=-0.09375d0*q2*bs*saz2(k)*ej2top(k)
            endif
            if(ldown)then
              eiaz2u(im,k)=-eiaz2u(im,k)
            endif
            if(lbelowhorizon)then
              eiaz1u(im,k)=-eiaz1u(im,k)
            endif
            if(ldown.neqv.lbelowhorizon)then
              eiaz1l(im,k)=-eiaz1l(im,k)
              eiaz1(im,k)=-eiaz1(im,k)
            endif
c
c*****compute total intensity at top of atmosphere
c
            totint(im,k)=eizero(im,k)+eiaz1(im,k)+eiaz2(im,k)
c
c     compute degree of polarization
            eitl(im,k)=eizrol(im,k)+eiaz1l(im,k)+eiaz2l(im,k)
            eitr(im,k)=eizror(im,k)+eiaz1r(im,k)+eiaz2r(im,k)
            eitu(im,k)=eiaz1u(im,k)+eiaz2u(im,k)
            pol(im,k)=dsqrt((eitl(im,k)-eitr(im,k))**2
     1            +eitu(im,k)**2)/totint(im,k)
c
            if (jprint(1) .eq. 1) write(33,133)im,k,eizero(im,k),
     *      eiaz1(im,k),eiaz2(im,k),totint(im,k),pol(im,k)
  133 format('evalitpol...im,k,ei0,ei1,ei2,totint,pol'/
     1        2i4,1p5e11.3)
            if (jprint(1) .eq. 1) write(33,136)eitl(im,k),
     *        eitr(im,k),eitu(im,k)
 136  format(8x,1p3e11.3)
  225    continue
         if(jprint(1).eq.1) write(33,1100)im,thta(im),
     1 zma1top(1),zma2top,zma3top,zma4top(1),
     2 ej1top(k),ej2top(k)
  105 continue
      if(jprint(1).eq.1) write(33,2000) gl,gr
c
c*****compute the downward diffuse intensity gg
c
      if(ipol .eq. 0)then
c     gg=c316*(gr*(admatx(1,jmuz)+1.0d0)+gl*admatx(3,jmuz))
        f1=c316*(admatx(1,jmuz)+1.0d0)
        f2=c316*admatx(3,jmuz)
        gg=f1*gr+f2*gl
        ggp=f1*grp+f2*glp
        gguptoa=f1*gruptoa+f2*gluptoa
      else
        f1=0.25d0*q1*(admatx(1,jmuz)+1.0d0+q)
        f2=0.25d0*q1*admatx(3,jmuz)
        gg=(f1*gr+f2*gl)/sdp
        ggp=(f1*grp+f2*glp)/sdp
        gguptoa=(f1*gruptoa+f2*gluptoa)/sdp
      endif
c 
      if(jprint(1).eq.1)then
        write(33,212)gl,gr,f1,f2,gg,q
 212  format('gl,gr,f1,f2,gg,q',6f8.4)
        write(33,213)(admatx(k,jmuz),k=1,3)
 213  format(1x,'admatx 1-3',3f8.4)
      endif
c  
      if(jprint(9).eq.1.or.write_iter_file)then
c**** call intsum to compute total intensity
        call intsum
        if(jprint(9).eq.1)call sumry(jmuz)
      endif
      if(lprtflux)then
        write(6,*)'*****************************'
        write(6,'(2(a4,1PE12.4),a9,1Pe12.4,a15,1Pe12.4)')
     &'fs=',fs,'gg=',gg,'gguptoa=',gguptoa,
     &'fs+gg+gguptoa=',fs+gg+gguptoa
        write(6,'(a5,1Pe12.4,a8,2(1Pe12.4),/,a18,1Pe12.4)')
     &'mu0=',emuz(jmuz),'engy_in=',engyin,1./engyin,
     &'engy_out/engy_in=',(fs+gg+gguptoa)/engyin
c emuz(jmuz)
        write(6,*)'*****************************'
      endif
c
      if(lcontrb)then
c compute the contribution function for this solar zenith angle
         call contrb(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
      endif
      return
 1000 format(1h1,t5,'debug print out from evalit',10x,
     1       'source functions(extrapolated) integrated over',
     2       1x,'optical thickness',///,
     2       1x,'imu',2x,'theta',5x,'zma1 ',9x,'zma2 ',9x,'zma3 ',9x,
     3       'zma4 ',9x,' ej1 ',9x,' ej2 ',//)
 1100 format(1x,i3,2x,f5.2,2x,d12.6,2x,d12.6,2x,d12.6,2x,d12.6,
     1       2x,d12.6,2x,d12.6)
 2000 format(///,1x,'gl= ',d12.6,/,1x,'gr= ',d12.6)
2010   format(1P6e12.4)
      end