


c***********************************************************************

      
      subroutine stalk_leaf(ss_sl,ss_ls)


c***********************************************************************


      complex ss_sl(4,4,4),ss_ls(4,4,4)
      complex s(5,4,4),t(2,4)
      complex ss,tt

      integer i,j,k


      do 10 i=1,4

        do 20 j=1,4

          do 30 k=1,5
            
           if(i .eq. 1) then
            
            if (j .eq. 1) then

              if (k .eq. 1) then
 
                call stalk_lphi_int(i,j,k, ss)

                call leaf_field_int(i,j,k, tt)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)

                call leaf_field_int(i,j,k, tt)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)
                
              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)


                

                call leaf_field_int(i,j,k, tt)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else

                print*, 'error in stalk_leaf a'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))
                
              else if(k .eq. 3) then

                 ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else

                print*, 'error in stalk_leaf b'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

              else if(k .eq. 2) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 4) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)
                
              else

                print*, 'error in  stalk_leaf c'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))

              else if (k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, ss)

                ss_sl(i,j,k)=s(i,j,k)*t(i,k)

                ss_ls(i,k,j)=conjg(s(i,j,k))*conjg(t(i,k))


              else

                print*, 'error in  stalk_leaf d'
                
              endif

            else
              
              print*, 'error in  stalk_leaf e'
                
            endif
            
            
           else if(i .eq. 2) then
            
            if (j .eq. 1) then

              if (k .eq. 1) then
 
                call stalk_lphi_int(i,j,k, ss)

                call leaf_field_int(i,j,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)

                call leaf_field_int(i,j,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)
                
              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)

                call leaf_field_int(i,j,k, tt)

                s(i,j,k)=ss

                t(i,k)=tt

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)
                
              else

                print*, 'error in stalk_leaf f'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)
                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else

                print*, 'error in stalk_leaf g'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

              else if(k .eq. 2) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)

                
              else if(k .eq. 4) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
                ss_ls(i,k,j)=ss_ls(i,k,2)
                
              else

                print*, 'error in  stalk_leaf h'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else if (k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

                ss_ls(i,k,j)=ss_ls(i,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(i,k))

                ss_ls(i,k,j)=conjg(s(i,j,k))*t(i,k)

              else

                print*, 'error in  stalk_leaf i'
                
              endif

            else
              
              print*, 'error in  stalk_leaf j'
                
            endif

          else if(i .eq. 3)
            
            ss_sl(i,j,k)=cmplx(0.0,0.0)

            ss_ls(i,k,j)=cmplx(0.0,0.0)
                
            endif
            

          else if(i .eq. 4) then
            
            if (j .eq. 1) then

              if (k .eq. 1) then
 
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))


              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))
                
              else

                print*, 'error in stalk_leaf k'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else

                print*, 'error in stalk_leaf l'

              endif

            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else if(k .eq. 2) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else if(k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else if(k .eq. 4) then

                ss_sl(i,j,k)=ss_sl(i,2,k)
                
              else

                print*, 'error in  stalk_leaf m'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else if (k .eq. 3) then

                ss_sl(i,j,k)=ss_sl(i,j,2)

              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_sl(i,j,k)=s(i,j,k)*conjg(t(1,k))

              else

                print*, 'error in  stalk_leaf n'
                
              endif

            else
              
              print*, 'error in  stalk_leaf o'
                
            endif

          else if(i .eq. 5) then
            
            if (j .eq. 1) then

              if (k .eq. 1) then
 
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)

              else if(k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)

              else if(k .eq. 3) then

                ss_ls(4,k,j)=ss_ls(4,2,j)
                
              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)
                
              else

                print*, 'error in stalk_leaf p'

              endif

            else if(j .eq. 2) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)
                
              else if(k .eq. 2) then
                
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)
                
              else if(k .eq. 3) then

                ss_ls(4,k,j)=ss_ls(4,2,j)

              else if(k .eq. 4) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)

              else

                print*, 'error in stalk_leaf q'

              endif
              
            else if(j .eq. 3) then

              if(k .eq. 1) then

                ss_ls(4,k,j)=ss_ls(4,k,2)

              else if(k .eq. 2) then

                ss_ls(4,k,j)=ss_ls(4,k,2)

              else if(k .eq. 3) then

                ss_ls(4,k,j)=ss_ls(4,k,2)

              else if(k .eq. 4) then

                ss_ls(4,k,j)=ss_ls(4,k,2)
                
              else

                print*, 'error in  stalk_leaf r'

              endif

            else if (j .eq. 4) then

              if (k .eq. 1) then

                call stalk_lphi_int(i,j,k, ss)

                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)
                
              else if (k .eq. 2) then

                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)

              else if (k .eq. 3) then

                ss_ls(4,k,j)=ss_ls(4,2,j)
               
              else if (k .eq. 4) then
                
                call stalk_lphi_int(i,j,k, ss)
                
                s(i,j,k)=ss

                ss_ls(4,k,j)=conjg(s(i,j,k))*t(2,k)

              else

                print*, 'error in  stalk_leaf s'
                
              endif

            else
              
              print*, 'error in  stalk_leaf t'
                
            endif

          else
            
              print*, 'error in  stalk_leaf u'

          endif

          
   30     continue  

   20   continue

   10 continue

      return

      end






c************************************************************************

      subroutine stalk_lphi_int(i,j,k, p)

c************************************************************************

      complex h_sl,p
      real pi,mpi
      integer i,j,k

      common /index/ i,j,k

      pi=4.0*atan(1.0)
      mpi=-pi
      

      external h_sl
      call qgausx(h_sl,mpi,pi,p)

      return

      end
      

c************************************************************************

      function h_sl(phii)

c************************************************************************      

      complex h_sl,s,g_sl

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real phii
      
      external g_sl
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      phi_s=phii

      call qgausy(g_sl,amin,amax,s)

      h_sl=s

      return

      end

c************************************************************************

      function g_sl(aa)

c************************************************************************


      complex g_sl,i_sl,s

      real amin,amax,z0min,z0max
      real phi_s,a_s,z0
      real aa

      external i_sl
      
      common /stalk_limits/ amin,amax,z0min,z0max
      common /stalk_state/ phi_s,a_s,z0

      a_s=aa

      call qgausy(i_sl,z0min,z0max,s)

      g_sl=s

      return

      end

c************************************************************************

      function i_sl(zz0)

c************************************************************************


      complex i_sl,stalk_phil_ave

      real phi_s,a_s,z0
      real zz0
      
      common /stalk_state/ phi_s,a_s,z0
      common /index/ i,j,k
      
      z0=zz0
      
      i_sl=stalk_phil_ave(phi_s,a_s,z0,i,j,k)

      return

      end

c************************************************************************

      function stalk_phil_ave(alpha_s,a,z0,i,j,k)

c************************************************************************
      complex p11vva_s,p33vva_s,p13vva_s
      complex p11vvb_s,p33vvb_s,p13vvb_s
      complex p11vvc_s,p33vvc_s,p13vvc_s
      complex p11vvc_s,p33vvd_s,p13vvd_s

      complex p11hha_s,p33hha_s,p13hha_s
      complex p11hhb_s,p33hhb_s,p13hhb_s
      complex p11hhc_s,p33hhc_s,p13hhc_s
      complex p11hhc_s,p33hhd_s,p13hhd_s

      complex p11hva_s,p33hva_s,p13hva_s
      complex p11hvb_s,p33hvb_s,p13hvb_s
      complex p11hvc_s,p33hvc_s,p13hvc_s
      complex p11hvc_s,p33hvd_s,p13hvd_s

      complex p11a,p12a,p13a,p22a,p23a,p33a
      complex p11b,p12b,p13b,p22b,p23b,p33b
      complex p11c,p12c,p13c,p22c,p23c,p33c
      complex p11d,p12d,p13d,p22d,p23d,p33d

      complex stalk_state_ave
      
      real dm,z0
      real a,alpha_s
      real theta_i,phi_i,k_0
      real zero
      real stalk_alpha_pdf,stalk_a_pdf,stalk_z0_pdf
      real dm,del_leaf,zlc
      
      integer i,j,k

      external p11vva_s,p33vva_s,p13vva_s
      external p11vvb_s,p33vvb_s,p13vvb_s
      external p11vvc_s,p33vvc_s,p13vvc_s
      external p11vvc_s,p33vvd_s,p13vvd_s

      external p11hha_s,p33hha_s,p13hha_s
      external p11hhb_s,p33hhb_s,p13hhb_s
      external p11hhc_s,p33hhc_s,p13hhc_s
      external p11hhc_s,p33hhd_s,p13hhd_s

      external p11hva_s,p33hva_s,p13hva_s
      external p11hvb_s,p33hvb_s,p13hvb_s
      external p11hvc_s,p33hvc_s,p13hvc_s
      external p11hvc_s,p33hvd_s,p13hvd_s      
      
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      zero=0.0
      

      if (i .eq. 1) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)


          else if (k .eq. 2) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &        *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
             
          else if (k .eq. 3) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)

          else if (k .eq. 4) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
                      
          else

            print*, 'error in stalk_phil_ave a'

          endif

        else if (j .eq. 2) then
          
          if (k .eq. 1) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)
            
          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else if (k .eq. 3)


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else if (k .eq. 4)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else

            print*, 'error in stalk_phil_ave b'

          endif

        else if (j .eq. 3) then
          
          if (k .eq. 1) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)
            
          else if (k .eq. 2) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else if (k .eq. 3)


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else if (k .eq. 4)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else

            print*, 'error in stalk_phil_ave c'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivva_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)


          else if (k .eq. 2) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
             
          else if (k .eq. 3) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)

          else if (k .eq. 4) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phivvd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
                      
          else

            print*, 'error in stalk_phil_ave d'

          endif
          
        else

          print*, 'error in stalk_phil_ave e'

        endif

      else if (i .eq. 2) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave f'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihhc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave g'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave h'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihha_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phihhd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)            
        
           else

             print*, 'error in stalk_phil_ave i'

           endif

         else

           print*, 'error in stalk_phil_ave j'

         endif
         

       else if (i .eq. 3) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hva_s,dm,z0,p21a)

            call qgaus(p23hva_s,dm,z0,p23a)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21a-sin(theta_i)*p23a)*conjg(phi_l)

         else

           print*, 'error in stalk_phil_ave k'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvb_s,dm,z0,p21b)

            call qgaus(p23hvb_s,dm,z0,p23b)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21b-sin(theta_i)*p23b)*conjg(phi_l)

         else

           print*, 'error in stalk_phil_ave l'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvc_s,dm,z0,p21c)

            call qgaus(p23hvc_s,dm,z0,p23c)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (cos(theta_i)*p21c-sin(theta_i)*p23c)*conjg(phi_l)

         else
           

           print*, 'error in stalk_phil_ave m'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvb_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

         
         else if (k .eq. 3) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvc_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p21hvd_s,dm,z0,p21d)

            call qgaus(p23hvd_s,dm,z0,p23d)

            call phihvd_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*
     &       (-cos(theta_i)*p21d-sin(theta_i)*p23d)*conjg(phi_l)

           else

             print*, 'error in stalk_phil_ave n'

           endif

         else

           print*, 'error in stalk_phil_ave o'

         endif


      else if (i .eq. 4) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hha_s,dm,z0,p22a)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22a*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave p'

         endif

       else if (j .eq. 2)

          if (k .eq. 1) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivvc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhb_s,dm,z0,p22b)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22b*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave q'

         endif

       else if (j .eq. 3)

          if (k .eq. 1) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhc_s,dm,z0,p22c)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22c*conjg(phi_l)            
        
         else

           print*, 'error in stalk_phil_ave r'

         endif

       else if (j .eq. 4)

          if (k .eq. 1) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivva_l(phi_l)            

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)

          else if (k .eq. 2) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvb_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
         
         else if (k .eq. 3) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvc_l(phi_l)  

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)
                     
         else if (k .eq. 4) then

            call qgaus(p22hhd_s,dm,z0,p22d)

            call phivvd_l(phi_l)
            
            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*p22d*conjg(phi_l)            
        
           else

             print*, 'error in stalk_phil_ave s'

           endif

         else

           print*, 'error in stalk_phil_ave t'

         endif
         

      else if (i .eq. 5) then

        if (j .eq. 1) then

          if (k .eq. 1) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)


          else if (k .eq. 2) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &        *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
             
          else if (k .eq. 3) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
             *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)

          else if (k .eq. 4) then

            call qgaus(p11vva_s,dm,z0,p11a)

            call qgaus(p33vva_s,dm,z0,p33a)

            call qgaus(p13vva_s,dm,z0,p13a)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11a
     &       -2.0*cos(theta_i)*sin(theta_i)*p13a
     &       +sin(theta_i)**2*p33a)
                      
          else

            print*, 'error in stalk_phil_ave u'

          endif

        else if (j .eq. 2) then
          
          if (k .eq. 1) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)
            
          else if (k .eq. 2) then


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else if (k .eq. 3)


            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else if (k .eq. 4)

            call qgaus(p11vvb_s,dm,z0,p11b)

            call qgaus(p33vvb_s,dm,z0,p33b)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11b+sin(theta_i)**2*p33b)

          else

            print*, 'error in stalk_phil_ave v'

          endif

        else if (j .eq. 3) then
          
          if (k .eq. 1) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)
            
          else if (k .eq. 2) then


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else if (k .eq. 3)


            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else if (k .eq. 4)

            call qgaus(p11vvc_s,dm,z0,p11c)

            call qgaus(p33vvc_s,dm,z0,p33c)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha_s)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (-cos(theta_i)**2*p11c+sin(theta_i)**2*p33c)

          else

            print*, 'error in stalk_phil_ave w'

          endif

        else if (j .eq. 4) then

          if (k .eq. 1) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihha_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)


          else if (k .eq. 2) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhb_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
             
          else if (k .eq. 3) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhc_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)

          else if (k .eq. 4) then

            call qgaus(p11vvd_s,dm,z0,p11d)

            call qgaus(p33vvd_s,dm,z0,p33d)

            call qgaus(p13vvd_s,dm,z0,p13d)

            call phihhd_l(phi_l)

            stalk_phil_ave=stalk_alpha_pdf(alpha)
     &       *stalk_a_pdf(a)*stalk_z0_pdf(z0)*conjg(phi_l)*
     &       (cos(theta_i)**2*p11d
     &       +2.0*cos(theta_i)*sin(theta_i)*p13d
     &       +sin(theta_i)**2*p33d)
                      
          else

            print*, 'error in stalk_phil_ave x'

          endif
          
        else

          print*, 'error in stalk_phil_ave y'

        endif

      else

        print*, 'error in stalk_phil_ave z'

      endif

      return

      end

      
c************************************************************************

      subroutine phivva_l(r)

c************************************************************************

      complex ii,phia
      complex b,r
      complex mv_leaf,mh_leaf

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z_lc)
     &+ 2.0*tau_v(z_lc)/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))
     &+2.0*mv_leaf/cos(theta_i)

      r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end

      

c************************************************************************

      subroutine phihha_l(r)

c************************************************************************

      complex ii,phia
      complex b,r
      complex mv_leaf,mh_leaf

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z_lc)
     &+ 2.0*tau_h(z_lc)/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))
     &+2.0*mh_leaf/cos(theta_i)

      r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end
      
      
c************************************************************************

      subroutine phihva_l(r)

c************************************************************************

      complex ii,phia
      complex b,r
      complex mv_leaf,mh_leaf

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phia=cexp(ii*2.0*k_0*(sin(theta_i)*x+cos(theta_i)*z_lc)
     &+(tau_v(z_lc)+tau_h(z_lc))/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a+cos(theta_i))
     &+(mv_leaf+mh_leaf)/cos(theta_i)

      r=phia*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end




c************************************************************************

      subroutine phivvb_l(r)

c************************************************************************

      complex ii,phib
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_v

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phib=cexp(ii*2.0*k_0*sin(theta_i)*x)
     &*cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_v(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a

      r=r_v(theta_i)*
     & phib*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end


c************************************************************************

      subroutine phihhb_l(r)

c************************************************************************

      complex ii,phib
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_h

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phib=cexp(ii*2.0*k_0*sin(theta_i)*x)
     &*cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_h(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a

      r=r_h(theta_i)*
     & phib*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end

c************************************************************************

      subroutine phihvb_l(r)

c************************************************************************

      complex ii,phib
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_v

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phib=cexp(ii*2.0*k_0*sin(theta_i)*x)
     & +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     & *cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_v(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
      +(mh_leaf-mv_leaf)/cos(theta_i)
      
      r=r_v(theta_i)*
     & phib*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end


c************************************************************************

      subroutine phivvc_l(r)

c************************************************************************

      complex ii,phic
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_v

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phic=cexp(ii*2.0*k_0*sin(theta_i)*x)
     &*cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_v(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a

      r=r_v(theta_i)*
     & phic*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end

      
c************************************************************************

      subroutine phihhc_l(r)

c************************************************************************

      complex ii,phic
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_h

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phic=cexp(ii*2.0*k_0*sin(theta_i)*x)
     &*cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_h(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a

      r=r_h(theta_i)*
     & phic*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end


      
c************************************************************************

      subroutine phihvc_l(r)

c************************************************************************

      complex ii,phic
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_v

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phic=cexp(ii*2.0*k_0*sin(theta_i)*x)
     & +(tau_h(z_lc)-tau_v(z_lc))/cos(theta_i)
     & *cexp(ii*2.0*k_0*cos(theta_i)*dm+ 2.0*tau_v(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*sin(theta_i)*cos(alpha_s)*a
      +(mh_leaf-mv_leaf)/cos(theta_i)
      
      r=r_v(theta_i)*
     & phic*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b

      return

      end

c************************************************************************

      subroutine phivvd_l(r)

c************************************************************************

      complex ii,phid
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_v
      
      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc

      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & -2.0*tau_v(z_lc)/cos(theta_i))
     & *cexp(ii*4.0*k_0*cos(theta_i)*dm+4.0*tau_v(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))
     &-2.0*mv_leaf/cos(theta_i)

      r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
     & *r_v(theta_i)**2
      
      return

      end
      
c************************************************************************

      subroutine phihhd_l(r)

c************************************************************************

      complex ii,phid
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_h
      
      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc
      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & -2.0*tau_h(z_lc)/cos(theta_i))
     & *cexp(ii*4.0*k_0*cos(theta_i)*dm+4.0*tau_h(dm)/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))
     & -2.0*mh_leaf/cos(theta_i)

      r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
     & *r_h(theta_i)**2
      
      return

      end

      

c************************************************************************

      subroutine phihvd_l(r)

c************************************************************************

      complex ii,phid
      complex b,r
      complex mv_leaf,mh_leaf
      complex r_h,r_v

      real a,alpha_s
      real theta_i,phi_i,k_0
      real dm,del_leaf,z_lc
      
      common /leaf_prop/ mv_leaf,mh_leaf
      common /stalk_struct/ a,alpha_s
      common /incident/ theta_i,phi_i,k_0
      common /layer_struct/ dm,del_leaf,z_lc

      ii=cmplx(0.0,1.0)

      rho=a*(z_lc-dm)

      x=rho*cos(alpha_s)

      phid=cexp(ii*2.0*k_0*(sin(theta_i)*x-cos(theta_i)*z_lc)
     & -(tau_v(z_lc)+tau_h(z_lc))/cos(theta_i))
     & *cexp(ii*4.0*k_0*cos(theta_i)*dm
     & +2.0*(tau_v(dm)+tau_h(dm))/cos(theta_i))
      
      b=ii*2.0*k_0*(sin(theta_i)*cos(alpha_s)*a-cos(theta_i))
     & -(mv_leaf+mh_leaf)/cos(theta_i)

      r=phid*(cexp(b*(del_leaf/2))-cexp(-b*(del_leaf/2)))/b
     & *r_h(theta_i)*r_v(theta_i)
      
      return

      end
      
