c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/contrb.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
         subroutine contrb(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1)
c***********************************************************************
c
c     subroutine contrb
c
c     purpose
c        compute the contribution function.  This is defined as the
c        integrand of the intensity integral, such that the integral of
c        the contribution function will equal the intensity.  So, it
c        is the source function times the attenuation out of the
c        atmosphere.  Much of this routine is taken from evalitpol.
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 contrb(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
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 last mod.: Mar. 14, 2000...def
c purpose: Change sb and sbp to single elements.
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
c      include 'czfunc.cmn'
c      include 'uptoa.cmn'
c      include 'energy.cmn'
c      include 'consts.cmn'
c      include 'buff3.cmn'
      include 'in.cmn'
      include 'switches.cmn'
      include 'contrl.cmn'
      include 'eks.cmn'
      include 'emm.cmn'
      include 'czsing.cmn'
      include 'depolt.cmn'
      include 'out.cmn'
      include 'thkns.cmn'
      include 'contrb.cmn'
c local
      integer i,im,k
      real*8 a(4),b(4),c(4),d(4)
      real*8 hold(4),xhold(4),trans
      real*8 as,bs,atb,absq
      real*8 dJ1,dJ2,dJ0,sum,const1,const2
c
      if(ipol.eq.0)then
        const1=0.09375d0
        const2=1.0d0
      else
        const1=0.125d0*(q1/sdp)
        const2=q2
      endif
c
c extrapolation of source functions has already been done in evalitpol
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*****calculate the azimuth -independent intensity eizero at current
c*****value of the polar angle. see eqs(3.7) and (4.16)
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
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*****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)
c
c
c*****do integration over optical thickness
c
         do 110 i=1,ncp1
c transmission term from each level to the top of the atmosphere
            trans=dexp(-tt(i)/emu(imu))
c loop over the azimuth angles
c
            do k=1,iazmth
c
c**** store extrapolated z-matrix in hold(4)
c
               hold(1)=zma1(4,i)
               hold(2)=zma2(4,i)
               hold(3)=zma3(4,i)
               hold(4)=zma4(4,i)
               if(ipsudo.eq.2)then
                 hold(1)=hold(1)+zsp(i,k,im)
                 hold(4)=hold(4)+zsp(i,k,im)
               endif
c*****multiply matrix z(1-4) times the adjoint matrix admatx(=b).
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)
               dJ0=const1*xhold(4)/emu(im)
               if(ipsudo.lt.2)then
                 dJ1=-0.375d0*const2*absq*caz(k)*ej1(4,i)
                 dJ2=0.09375d0*const2*(atb*caz2(k))*ej2(4,i)/emu(im)
               else
                 dJ1=-0.375d0*const2*absq*caz(k)*(ej1(4,i)+zsp(i,k,im))
                 dJ2=0.09375d0*const2*(atb*caz2(k))*(ej2(4,i)+
     &zsp(i,k,im))/emu(im)
               endif
               if(ldown.ne.lbelowhorizon)then
                 dJ1=-dJ1
               endif
c
c*****compute total intensity at top of atmosphere
c
               dJ_trans(i,im,k)=trans*(dJ0+dJ1+dJ2)
            enddo
110      continue
  105 continue
c
c check the results
      do im=1,imu
         do k=1,iazmth
      sum=0.0d0
      do i=1,ncp1
c         dJ_trans(i,im,k)=dJ_trans(i,im,k)*dtsp(i)
c         sum=sum+dJ_trans(i,im,k)
         dJ_trans(i,im,k)=dJ_trans(i,im,k)
         sum=sum+dJ_trans(i,im,k)*dtsp(i)
      enddo
      write(6,*)'.contrb sum totint ',im,k,sum,totint(im,k),
     &sum/totint(im,k)-1.
      enddo
      enddo
c  
      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