c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/eva1pol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
         subroutine eva1pol(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1,itz)
c***********************************************************************
c
c     subroutine eva1pol
c
c     purpose
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 04/15/96...dave flittner
c     purpose: print upward flux at top of the atmosphere
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 
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. 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.
c
c***********************************************************************
c
      implicit none
      include 'parameter.inc'
c passed in
      integer jmuz,ncp1,itz
      real*8 zma1(4,202),zma2(4,202),zma3(4,202),zma4(4,202)
      real*8 ej1(4,202),ej2(4,202)
c commons
      include 'uptoa.cmn'
      include 'eks.cmn'
      include 'contrl.cmn'
      include 'out.cmn'
      include 'emm.cmn'
      include 'in.cmn'
      include 'depolt.cmn'
      include 'buff2.cmn'
      include 'czsing.cmn'
      include 'switches.cmn'
c local
      real*8 zma1top(max_az),zma2top,zma3top,zma4top(max_az),f1,f2
      real*8 ej1top(max_az),ej2top(max_az),a(4),b(4),c(4),d(4),atmu,tmu
      real*8 gr,gl,hold(4),xhold(4),c316,gruptoa,gluptoa,as,bs,atb,absq
      real*8 glp,grp
      integer i,itmaxp,im,k,index
c
      data c316/0.1875d0/
c
c     extrapolate source functions layer by layer
c
      itmaxp=itmax+1
      if (itz .le. itmaxp) then
        do i=1,ncp1
           zma1(4,i) = zma1(1,i)
           zma2(4,i) = zma2(1,i)
           zma3(4,i) = zma3(1,i)
           zma4(4,i) = zma4(1,i)
           ej1(4,i)  = ej1(1,i)
           ej2(4,i)  = ej2(1,i)
        enddo
      else
        if ((itmax.lt.3).or.lnoextrap) then
          if(itmax.lt.1)then
            index=1
          else
            index=min(itmax,3)
          endif
          do 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)
          enddo
        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
      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
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
         if(jprint(4).eq.1)write(33,130)
130      format('debug eva1pol. level zma1top zma1 tmu')
         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
            if(jprint(4).eq.1)write(33,131)i,zma1top(1),zma1(4,i),tmu
131         format(i4,1P3e12.4)
 110     continue
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
c*****multiply matrix ztop(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)
            else
              eizero(im,k)=0.125d0*(q1/sdp)*xhold(4)/emu(im)
            endif
c
            e0za(itz,im,k)=eizero(im,k)
c
            if(im.eq.1.and.jprint(4).eq.1)then
              write(33,515)itz,im,emu(im),eizero(im,k),e0za(itz,im,k)
 515          format('sub eva1pol.itz,im,emu,eizero,e0za',2i5,1x,3f8.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
         atb=as*bs
         absq=dsqrt(atb)*emuz(jmuz)
         atmu=as*emuz(jmuz)
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)
            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)
            endif
            if(ldown.neqv.lbelowhorizon)eiaz1(im,k)=-eiaz1(im,k)
c
            e1za(itz,im,k)=eiaz1(im,k)
            e2za(itz,im,k)=eiaz2(im,k)
c
c*****compute total intensity at top of atmosphere
c
            totint(im,k)=eizero(im,k)+eiaz1(im,k)+eiaz2(im,k)
c
            if (jprint(4) .eq. 1)
     *    write(33,133)itz,im,k,eizero(im,k),eiaz1(im,k),
     *             eiaz2(im,k),totint(im,k)
  133 format('eva1.itz,im,k,ez0,ez1,ez2,tot',3i4,4f9.4)
  225    continue
  105 continue
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 
      ggz(itz)=gg
      ggzuptoa(itz)=gguptoa
c     for debugging
      if (jprint(4) .eq. 1) write(33,212)gl,gr,f1,f2,gg
      if (jprint(4) .eq. 1) write(33,312)(e0za(itz,im,1),im=1,imu)
      if (jprint(4) .eq. 1) write(33,313)ggz(itz)
 312  format('e0za',1p6e11.3)
 313  format('ggz',1pe11.3)
 212  format('eva1pol...gl,gr,f1,f2,gg',5f8.4)
c  
      return
      end