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


      subroutine canopy_extinction 


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

      integer true
      parameter(true=1)
      
      integer nl,ns
      parameter(nl=250)
      parameter(ns=5)
c     consider replacing the above with an include!
      
      complex tauv_v(nl),tauh_v(nl)
      
      integer specie,layer
      
      complex mv_leaf,mh_leaf
      common/leaf_prop/ mv_leaf,mh_leaf

      complex mv_grain,mh_grain
      common/grain_prop/ mv_grain,mh_grain

      complex tauv(nl),tauh(nl)      
      common /extinction_data/ tauv,tauh

      integer N_species
      common /canopy_data/ N_species
      
      integer N_layer      
      common /layer_segments/ N_layer
      
      integer leaf_flag,stalk_flag,grain_flag
      integer gl_flag,gs_flag,sl_flag,ll_flag
      common /scat_term_flags/leaf_flag,stalk_flag,grain_flag,
     &     gl_flag,gs_flag,sl_flag,ll_flag
      
      include 'error.incld'
      include 'status.incld'

      if (N_layer .gt. nl) then

        write(error,*)'N_layer = ',nl,' too large in sub extinction'
        code=28
        call the_death_code(die,error,code)

      else if (N_species .gt. ns) then   
      
        write(error,*)'N_specie = ',ns,' too large in sub extinction'
        code=29
        call the_death_code(die,error,code)

      else
         
        continue

      endif

      call status_report('canopy_extinction',null,null,null,start)

      do  layer=1,N_layer

        tauv(layer)=cmplx(0.0,0.0)
        tauh(layer)=cmplx(0.0,0.0)

      enddo

      do  specie=1,N_species

        call load_specie_data(specie)

        call plant_extinction(N_layer, tauv_v,tauh_v)

c        print*,'layer,tauv,tauh'
        
        do  layer=1,N_layer

          tauv(layer)=tauv(layer)+tauv_v(layer)
          tauh(layer)=tauh(layer)+tauh_v(layer)

c          print*,layer,tauv(layer),tauh(layer)

        enddo

      enddo

      call leaf_center_prop(mv_leaf,mh_leaf)
      
      if(grain_flag .eq. true)then
        
        call grain_center_prop(mv_grain,mh_grain)
      
      endif
      
      call status_report('canopy_extinction',null,null,null,stop)
      
      return
      end

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

      function tau_v(z)

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

      integer nl
      parameter(nl=250)

      complex tauv(nl),tauh(nl)      
      common /extinction_data/ tauv,tauh

      integer N_layer      
      common /layer_segments/ N_layer

      real d,del_leaf,z_lc,hole
      common /layer_struct/ d,del_leaf,z_lc,hole
      
      complex tau_v_lc,tau_h_lc
      common /tau_lc/ tau_v_lc,tau_h_lc

      complex tau_v_gc,tau_h_gc
      common /tau_gc/ tau_v_gc,tau_h_gc
      
      real del_grain,z_gc
      common /grain_layer/del_grain,z_gc
      
      include 'error.incld'

      complex tau_v,m,b 

      real z,z_delta

      integer layer

      if(z .eq. -d)then

        tau_v=tauv(N_layer)

      else if(z .eq. z_lc)then

        tau_v=tau_v_lc

      else if(z .eq. z_gc)then
                
        tau_v=tau_v_gc
        
      else if(z .gt. 0.0) then
        
        write(error, *)'z = ',z,' too large in func tau_v'
        code=30
        call the_death_code(die,error,code)

      else if(z .lt. -d) then

        write(error,*)'N_specie = ',z,' too large in sub extinction'
        code=31
        call the_death_code(die,error,code)

      else

        z_delta=-d/N_layer
        layer=int(z/z_delta)+1

        if (layer .eq. N_layer+1) then

          layer=N_layer

        endif

        if(layer .eq. 1)then
        
          m=tauv(layer)/z_delta

        else

          m=(tauv(layer)-tauv(layer-1))/z_delta

        endif

        b=tauv(layer)-m*(layer*z_delta)

        tau_v=m*z+b

      endif

      return
      end
      
      
c***********************************************************************

      function tau_h(z)

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

      integer nl
      parameter(nl=250)

      complex tauv(nl),tauh(nl)      
      common /extinction_data/ tauv,tauh

      integer N_layer      
      common /layer_segments/ N_layer

      real d,del_leaf,z_lc,hole
      common /layer_struct/ d,del_leaf,z_lc,hole
      
      complex tau_v_lc,tau_h_lc
      common /tau_lc/ tau_v_lc,tau_h_lc

      complex tau_v_gc,tau_h_gc
      common /tau_gc/ tau_v_gc,tau_h_gc
      
      real del_grain,z_gc
      common /grain_layer/del_grain,z_gc
      
      include 'error.incld'

      complex tau_h,m,b 

      real z,z_delta

      integer layer


      if(z .eq. -d)then

        tau_h=tauh(N_layer)
        
      else if(z .eq. z_lc)then

        tau_h=tau_h_lc

      else if(z .eq. z_gc)then

        tau_h=tau_h_gc
              
      else if(z .gt. 0.0) then
        
        write(error,*)'z = ',z,' too large in func tau_h'
        code=32
        call the_death_code(die,error,code)

      else if(z .lt. -d) then

        write(error,*)'N_specie = ',z,' too large in sub extinction'
        code=33
        call the_death_code(die,error,code)

      else

        z_delta=-d/N_layer
        layer=int(z/z_delta)+1

        if (layer .eq. N_layer+1) then
        
          layer=N_layer
        
        endif

        if(layer .eq. 1)then
        
          m=tauh(layer)/z_delta

        else

          m=(tauh(layer)-tauh(layer-1))/z_delta

        endif

        b=tauh(layer)-m*(layer*z_delta)

        tau_h=m*z+b

      endif
        
      return
      end


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

      subroutine plant_extinction(N_layer, tauv,tauh)

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

      integer true
      parameter(true=1)
      
      integer nl
      parameter(nl=250)

      complex tauv(nl),tauh(nl)
      complex tauh_g(nl),tauv_g(nl)
      complex tauh_s(nl),tauv_s(nl)
      complex tauh_l(nl),tauv_l(nl)

      common/tau_save/ tauh_l,tauv_l,tauh_s,tauv_s,tauh_g,tauv_g
      
      integer N_layer,layer

      integer leaf_flag,stalk_flag,grain_flag
      integer gl_flag,gs_flag,sl_flag,ll_flag
      common /scat_term_flags/leaf_flag,stalk_flag,grain_flag,
     &     gl_flag,gs_flag,sl_flag,ll_flag

      include 'status.incld'

      if(grain_flag .eq. true)then

        call grain_extinction(N_layer, tauv_g,tauh_g)

        call status_report('grain_extinction',null,null,null,elapsed)

      endif

      if(stalk_flag .eq. true)then

        call stalk_extinction(N_layer, tauv_s,tauh_s)
      
        call status_report('stalk_extinction',null,null,null,elapsed)
      
      endif

      if(leaf_flag .eq. true)then

        call leaf_extinction(N_layer, tauv_l,tauh_l)

        call status_report('leaf_extinction',null,null,null,elapsed)
      
      endif
      
c        print*,'layer,tauh,tauh_s,tauh_l,
c     &       tauv_s,tauv_l'
        
      do  layer=1,N_layer
        
        tauh(layer)=tauh_g(layer)+tauh_s(layer)+tauh_l(layer)
        tauv(layer)=tauv_g(layer)+tauv_s(layer)+tauv_l(layer)

c        print*,layer,tauh(layer),tauh_s(layer),tauh_l(layer),
c     &       layer,tauv(layer),tauv_s(layer),tauv_l(layer)

      enddo
      
      return
      end
      
c***********************************************************************
      
      subroutine leaf_center_prop(mv_leaf,mh_leaf)
      
c***********************************************************************

      integer nl
      parameter(nl=250)

      complex tauv(nl),tauh(nl)      
      common /extinction_data/ tauv,tauh

      integer N_layer      
      common /layer_segments/ N_layer

      real d,del_leaf,z_lc,hole
      common /layer_struct/ d,del_leaf,z_lc,hole

      complex tau_v_lc,tau_h_lc
      common /tau_lc/ tau_v_lc,tau_h_lc

      complex mv_leaf,mh_leaf
      complex m_v,m_h,b_v,b_h


      real z_delta,z_plus,z_minus

      integer z_lc_layer
      
      include 'error.incld'
      
      z_delta=-d/N_layer
      z_lc_layer=int((z_lc)/z_delta)

      z_plus=int((z_lc+0.5*del_leaf)/z_delta)
      z_minus=int((z_lc-0.5*del_leaf)/z_delta)

      if(z_plus .eq. z_lc_layer)then

        z_plus=z_lc_layer-1

      else if(z_minus .eq. z_lc_layer)then

        z_plus=z_lc_layer+1

      endif
      
C      z_plus=((z_lc_layer-1)-1)*z_delta
C      z_minus=((z_lc_layer+1)-1)*z_delta
C      
C      if((z_plus .gt. z_lc+del_leaf/2.0)
C     &     .or.(z_minus .lt. z_lc-del_leaf/2.0))then
C        
C        write(error,*)'z_delta = ',z_delta,' apparently .gt. ',
C     &       'del_leaf/3 = ',del_leaf/3.0,'in sub leaf_center_prop'
C        code=78
C        call the_death_code(die,error,code)
C
C      endif
      
      mv_leaf=(tauv(z_minus)-tauv(z_plus))
     &     /((z_minus-z_plus)*z_delta)
      
      mh_leaf=(tauh(z_minus)-tauh(z_plus))
     &     /((z_minus-z_plus)*z_delta)

C      m_v=(tauv(z_lc_layer)-tauv(z_lc_layer-1))/z_delta
C      m_h=(tauh(z_lc_layer)-tauh(z_lc_layer-1))/z_delta
C      
C      b_v=tauv(z_lc_layer)-m_v*(z_lc_layer*z_delta)
C      b_h=tauh(z_lc_layer)-m_h*(z_lc_layer*z_delta)
      
      tau_v_lc=tauv(z_plus)-mv_leaf*(z_plus*z_delta-z_lc)
      tau_h_lc=tauh(z_plus)-mh_leaf*(z_plus*z_delta-z_lc)
      
C      mv_leaf=(tauv(z_lc_layer+1)-tauv(z_lc_layer-1))
C     &     /(z_minus-z_plus)
C      
C      mh_leaf=(tauh(z_lc_layer+1)-tauh(z_lc_layer-1))
C     &     /(z_minus-z_plus)
C
C      m_v=(tauv(z_lc_layer)-tauv(z_lc_layer-1))/z_delta
C      m_h=(tauh(z_lc_layer)-tauh(z_lc_layer-1))/z_delta
C      
C      b_v=tauv(z_lc_layer)-m_v*(z_lc_layer*z_delta)
C      b_h=tauh(z_lc_layer)-m_h*(z_lc_layer*z_delta)
C      
C      tau_v_lc=m_v*z_lc+b_v
C      tau_h_lc=m_h*z_lc+b_h

      print*,'tau_h_lc,mh_leaf,z_lc_layer,z_plus,z_minus'
      print*,tau_h_lc,mh_leaf,z_lc_layer,z_plus,z_minus

      print*,'tau_v_lc,mv_leaf,z_lc_layer,z_plus,z_minus'
      print*,tau_v_lc,mv_leaf,z_lc_layer,z_plus,z_minus

      
      return
      end

      
c***********************************************************************
      
      subroutine grain_center_prop(mv_grain,mh_grain)
      
c***********************************************************************

      integer nl
      parameter(nl=250)

      complex tauv(nl),tauh(nl)      
      common /extinction_data/ tauv,tauh

      integer N_layer      
      common /layer_segments/ N_layer

      real d,del_leaf,z_lc,hole
      common /layer_struct/ d,del_leaf,z_lc,hole

      complex tau_v_gc,tau_h_gc
      common /tau_gc/ tau_v_gc,tau_h_gc

      complex mv_grain,mh_grain
      complex m_v,m_h,b_v,b_h
      
      real del_grain,z_gc
      common /grain_layer/del_grain,z_gc
      
      real z_delta,z_plus,z_minus

      integer z_gc_layer
      
      include 'error.incld'
      
      z_delta=-d/N_layer
      z_gc_layer=int(z_gc/z_delta)+1
      z_plus=((z_gc_layer-1)-1)*z_delta
      z_minus=((z_gc_layer+1)-1)*z_delta
      
      if((z_plus .gt. z_gc+del_grain/2.0)
     &     .or.(z_minus .lt. z_gc-del_grain/2.0))then
        
        write(error,*)'z_delta = ',-z_delta,' apparently .gt. ',
     &       'del_grain/3 = ',del_grain/3.0,'in sub grain_center_prop'
        code=78
        call the_death_code(die,error,code)

      endif

      mv_grain=(tauv(z_gc_layer+1)-tauv(z_gc_layer-1))
     &     /(z_minus-z_plus)
      
      mh_grain=(tauh(z_gc_layer+1)-tauh(z_gc_layer-1))
     &     /(z_minus-z_plus)

      print*,'mv_grain,mh_grain',mv_grain,mh_grain
      print*,'tauv(z_gc_layer+1),tauv(z_gc_layer-1)'
      print*,tauv(z_gc_layer+1),tauv(z_gc_layer-1)
      print*,'z_gc_layer+1,z_gc_layer-1'
      print*,(z_gc_layer+1),(z_gc_layer-1)
      
      m_v=(tauv(z_gc_layer)-tauv(z_gc_layer-1))/z_delta
      m_h=(tauh(z_gc_layer)-tauh(z_gc_layer-1))/z_delta
      
      b_v=tauv(z_gc_layer)-m_v*(z_gc_layer*z_delta)
      b_h=tauh(z_gc_layer)-m_h*(z_gc_layer*z_delta)
      
      tau_v_gc=m_v*z_gc+b_v
      tau_h_gc=m_h*z_gc+b_h
      
      return
      end

Cc***********************************************************************
C      
C      subroutine  grain_extinction(N_layer, tauv_g,tauh_g)
C      
Cc***********************************************************************
C
C      integer nl
C      parameter(nl=250)
C      
C      integer N_layer,i
C
C      complex tauv_g(nl),tauh_g(nl)
C
C      do i=1,N_layer
C
C        tauv_g(i)=cmplx(0.0,0.0)
C
C        tauh_g(i)=cmplx(0.0,0.0)
C
C      enddo
C
C      return
C      end