!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutine canopy for:                     *!
!*    Calculation of canopy geometry & light absorption          *!
!*      with                                                     *!
!*     CALC_LA                                                   *!
!*     LIGHT_GROWTH                                              *!
!*     COV_AREA                                                  *!
!*     Light_1                                                   *!
!*     Light_2                                                   *!
!*     Light_3                                                   *!
!*     Light_4                                                   *!
!*     L_3_COH_LOOP                                              *!
!*     L_4_COH_LOOP                                              *!
!*    LIGHT_OUT_2                                                *!  
!*    CROWN_PROJ                                                 *! 
!*                                                               *!
!*                  Copyright (C) 1996-2018                      *!
!*     Potsdam Institute for Climate Impact Reserach (PIK)       *!
!*          Authors and contributors see AUTHOR file             *!
!*  This file is part of 4C and is licensed under BSD-2-Clause   *!
!*                   See LICENSE file or under:                  *!
!*     http://www.https://opensource.org/licenses/BSD-2-Clause   *!
!*                           Contact:                            *!
!*       https://gitlab.pik-potsdam.de/foresee/4C                *!
!*                                                               *!
!*****************************************************************!

!**********************************!
!*   SUBROUTINE CANOPY            *!
!**********************************!

SUBROUTINE CANOPY

  !*** Declaration part ***!

  USE data_out
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  integer i

  ! If no Cohorts on the patch, initialize properly
  IF( anz_coh == 0 ) THEN

      lowest_layer=0
      highest_layer=0
      vStruct%cumLAI= 0.
      vStruct%Irel  = 0.
      vStruct%sumBG = 0.
      Irelpool      = 0.
      BGpool        = 0.
      LAI           = 0.

      ! full light on the ground (layer = 0)
      ! Lightroutine 1,2
      vStruct(highest_layer)%Irel=1
      ! Lightroutine 3,4
      Irelpool(highest_layer)=1

      ! the whole patch is availabe for recruitment
      BGpool(highest_layer+1)=1
      BGpool(highest_layer+2)=1
      all_leaves_on=0

      ! Calculation of leaf area, lowest and highest layer, etc.
      ! for all cohorts in all respective layers
      CALL CALC_LA   ! leaf area etc. always calculate
      
      RETURN
  END IF

  ! Calculation of leaf area, lowest and highest layer, etc.
  ! for all cohorts in all respective layers
  CALL CALC_LA

  IF(flag_end.EQ.3) RETURN

   IF( flag_light == 1 )THEN

       CALL LIGHT_1

   ELSE IF ( flag_light == 2 ) THEN

       CALL LIGHT_2

   ELSE IF ( flag_light == 3 ) THEN

       CALL LIGHT_3

   ELSE IF ( flag_light == 4 ) THEN

       CALL LIGHT_4

   END IF
   DO i=1,anrspec
      ns = nrspec(i)
      IF(svar(ns)%act_sum_lai > svar(ns)%sum_lai) svar(ns)%sum_lai = svar(ns)%act_sum_lai
   ENDDO

  ! Determine relative light in the middle of each cohort canopy, the sla
  ! and the totFPAR per square meter patch and the total FPAR on the patch
  CALL LIGHT_GROWTH
  ! print relevant light parameters for the canopy for each layer and cohort
  if (time_out.gt.0 .and. out_flag_light.ne.0) CALL LIGHT_OUT_2

!------------------------------------------------
!------------------- SUBROUTINES ----------------
!------------------------------------------------
CONTAINS

SUBROUTINE CALC_LA

  ! Calculation of leaf area, lowest and highest layer, etc.
  ! for all cohorts in all respective layers

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: nl, i
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  ! auxiliary variable
  REAL    :: x    ! leaf area per crown unit [m**2/cm]

  vStruct%LA = 0.

  ! structure of the canopy is determined once at the start of the year
  ! initialisation
  IF(iday==1) THEN
      lowest_layer=250
      highest_layer=0
  END IF

  do i = 1, anrspec
     svar(nrspec(i))%act_sum_lai = 0.
  enddo

  p => pt%first

  DO WHILE (ASSOCIATED(p))
    ns = p%coh%species

    ! cohort loop for determination of lowest and highest canopy layer of the tree crown
    ! structure of the canopy must only be determined once at the start of the year
    IF(iday==1) THEN

       ! determine bottom of the crown in terms of  number of layers
       p%coh%botLayer = INT( p%coh%x_hbole / dz ) + 1

       ! determine top of the crown in terms of  number of layers
       IF (MODULO(p%coh%height,dz)==0.) THEN
           p%coh%topLayer = INT( p%coh%height /  dz )
       ELSE
           p%coh%topLayer = INT( p%coh%height /  dz ) + 1
       END IF

       ! remember the highest layer
       IF(p%coh%topLayer > highest_layer .AND. p%coh%toplayer < 250) THEN
          highest_layer=p%coh%topLayer

       ELSE IF(p%coh%toplayer >= 250) THEN
           if (.not.flag_mult8910) then
               CALL stop_mess(time,'FATAL EXCEPTION RAISED IN CANOPY CALC_LA')
               CALL error_mess(time,'maximal tree height of 125 m reached by cohort No.',REAL(p%coh%ident))
           endif  
           flag_end=3
           RETURN
       END IF

       !remember the lowest layer of the stand
       IF(p%coh%botLayer < lowest_layer) THEN
           lowest_layer=p%coh%botLayer
       END IF
    END IF

    p%coh%leafarea = 0.

    ! total leaf area of a tree in this cohort [m**2]
    IF((iday >= p%coh%day_bb) .AND. (iday <= spar(ns)%end_bb)) THEN
       p%coh%t_leaf = p%coh%med_sla * p%coh%x_fol

       ! amount of leaf area per tree in layers
       IF (p%coh%topLayer-p%coh%botLayer.GE.1) THEN
          ! now calculate leaf area per crown unit of this tree [m**2/cm]
          x = p%coh%t_leaf / ( p%coh%height - p%coh%x_hbole )
          p%coh%leafArea( p%coh%botLayer ) = ( dz - MODULO( p%coh%x_hbole, dz ) ) * x

          IF (MODULO(p%coh%height,dz)==0.) THEN
             p%coh%leafArea( p%coh%topLayer ) = dz * x
          ELSE
             p%coh%leafArea( p%coh%topLayer ) = MODULO( p%coh%height, dz ) * x
          END IF

          DO nl = p%coh%botLayer+1, p%coh%topLayer-1
             p%coh%leafArea(nl) = x * dz
          END DO
       ELSE
          p%coh%leafArea(p%coh%botLayer) = p%coh%t_leaf
       END IF
       ! Update vertical patch leaf area profile of the canopy
       DO nl = p%coh%botLayer, p%coh%topLayer
         vStruct(nl)%LA = vStruct(nl)%LA + p%coh%leafArea(nl) * p%coh%nTreeA
       END DO
    ELSE
        p%coh%leafArea=0.
    ENDIF
      IF(iday<=spar(ns)%end_bb) svar(ns)%act_sum_lai = svar(ns)%act_sum_lai + p%coh%ntreea*p%coh%t_leaf/kpatchsize
    p => p%next
  END DO
END SUBROUTINE CALC_LA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_GROWTH

  ! Determine  relative light in the middle of each cohort canopy, the sla,
  ! the total FPAR on the patch

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE
  integer help
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  totFPARsum=0    ! sum of all totFPAR's
  totFPARcan=0    ! sum of all totFPAR's for the canopy
  p => pt%first

  DO WHILE (ASSOCIATED(p))

      ns=p%coh%species

      ! the new average specific leaf area per cohort depends
      ! on the light regime in the middle of the canopy
      ! this is the SLA which is used for the leaf area distr. in the next year
      ! the new average specific leaf area per cohort depends on the
      ! mean light regime in the middle in the canopy

      ! IrelCan modifies the growthfunction
      IF(all_leaves_on==1) THEN
         
         select case (flag_light)
         case (1,2)
            p%coh%med_sla = spar(ns)%psla_min+spar(ns)%psla_a*&
                    (1-(vStruct(p%coh%toplayer)%Irel+vStruct(p%coh%botlayer)%Irel)/2.)
            p%coh%IrelCan = vStruct(p%coh%toplayer)%Irel
         case default
            p%coh%med_sla = spar(ns)%psla_min+spar(ns)%psla_a*&
                           (1-(p%coh%Irel(p%coh%topLayer)+p%coh%Irel(p%coh%botLayer))/2.)
            select case (ns)
            case (10)   ! Douglas fir
                help = p%coh%botLayer+2*(p%coh%toplayer - p%coh%botLayer) / 3
                p%coh%IrelCan = p%coh%Irel(help)
            case default
                help = vStruct(p%coh%toplayer)%SumBG
                if (help .gt. 0.) then
                    p%coh%IrelCan = p%coh%Irel(p%coh%toplayer)*MIN(kpatchsize/help, 1.)
                else
                    p%coh%IrelCan = p%coh%Irel(p%coh%toplayer)
               endif
            end select    ! ns 
         end select    ! flag_light
      END IF

      totFPARsum = totFPARsum + p%coh%totFPAR*p%coh%nTreeA
      IF (p%coh%species .le. nspec_tree .or. p%coh%species.eq.nspec_tree+2) totFPARcan = totFPARcan + p%coh%totFPAR*p%coh%nTreeA
      p => p%next
  END DO

END SUBROUTINE LIGHT_GROWTH

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE COV_AREA

  ! calculate coverage-area as fraction of the patchsize per tree and layer

  !*** Declaration part ***!
  USE data_climate
  USE data_par
  USE data_stand
  USE data_site

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  ! Variables to test restriction in light model 4
  REAL  :: y         ! potential shadow cast of the cohort [m]
  REAL  :: w         ! effective shadow cast of the cohort [m]
  REAL  :: l         ! side length of a coort layer [m]
  REAL  :: reqarea   ! area of the patch required for the shadow cast for all cohorts per layer
  INTEGER :: layer_flag ! remember the highest layer where first LM4 restriction occurs
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  y = dz/100/TAN(beta)
  lm3layer=0
  layer_flag=0

  DO i = highest_layer, lowest_layer, -1
      reqarea=0.
      p => pt%first
      DO WHILE (ASSOCIATED(p))

          p%coh%BG(i) = 0.
          ! only those trees that have leaves
          IF((iday >= p%coh%day_bb) .AND. (iday <= spar(p%coh%species)%end_bb) .AND. &
             i <= p%coh%topLayer .AND. i >= p%coh%botLayer) THEN

             IF (vStruct(i)%sumBG > kpatchsize) THEN

                 p%coh%BG(i)=p%coh%crown_area/vStruct(i)%sumBG

             ELSE

                 p%coh%BG(i)=p%coh%crown_area/kpatchsize

             END IF

             l = SQRT(p%coh%BG(i)*kpatchsize)

             reqarea = reqarea + l*y*p%coh%nTreeA

          END IF

         p => p%next

     END DO  ! cohorts

     IF( kpatchsize > vStruct(i)%sumBG .AND. reqarea /= 0) THEN

       w = y*(kpatchsize-vStruct(i)%sumBG)/reqarea

     ELSE

       w = 0

     END IF

     p => pt%first

     DO WHILE (ASSOCIATED(p) .AND. layer_flag.EQ.0)
         ! only those trees that have leaves
         IF((iday >= p%coh%day_bb) .AND. (iday <= spar(p%coh%species)%end_bb) .AND. &
             i <= p%coh%topLayer .AND. i >= p%coh%botLayer) THEN

             l = SQRT(p%coh%BG(i)*kpatchsize)
             ! layer from that on light model 3 is used instead of light model 4
             ! because of LM4 restrictions
             IF( y-w > w+l ) THEN
                 layer_flag=1
                 lm3layer = i
                 EXIT ! do loop
             END IF
         END IF
         p => p%next
     END DO

  END DO     ! layers

END SUBROUTINE COV_AREA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_1

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i, nl

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  ! auxiliary variables
  REAL    :: radSum       ! sum of absorbed radiation (help variable)
  REAL    :: pfext=0.6    ! extinction coefficient. Only for one specie.


  !*** Calculation part ***!
  ! Intialization radiation summator
  radSum         = 0.
  vStruct%cumLAI = 0.
  vStruct%Irel   = 0.
  ! Calculate cumulative leaf area index and absorbed radiation per layer
  ! using Lambert-Beer
  vStruct(highest_layer)%Irel=1

  DO i = highest_layer, lowest_layer, -1
    vStruct(i)%cumLAI  = vStruct(i)%LA/kPatchsize + vStruct(i+1)%cumLAI
    vStruct( i )%radFrac =  1. - Exp(-pfext * vStruct(i)%cumLAI) - radSum
    radSum = radSum + vStruct(i)%radFrac
    vStruct(i-1)%Irel=vStruct(i)%Irel-vStruct(i)%radFrac

  END DO

  ! Light intensitiy unto the ground
  DO i = lowest_layer - 2, 0, -1
      vStruct(i)%Irel=vStruct(i+1)%Irel
  END DO

  ! total LAI is simply the value of cumLAI at the forest floor
  LAI = vStruct(lowest_layer)%cumLAI
  IF(lai>laimax) laimax=lai

  ! Determine layer-specific & total fraction of PAR absorbed by this tree
  p => pt%first
  DO WHILE (ASSOCIATED(p))

    p%coh%totFPAR = 0.
    p%coh%FPAR = 0.

    DO nl = p%coh%botLayer, p%coh%topLayer

        p%coh%FPAR(nl) = p%coh%leafArea(nl) / vStruct(nl)%LA * vStruct(nl)%radFrac
        p%coh%totFPAR  = p%coh%totFPAR + p%coh%FPAR(nl)

    END DO
    p => p%next
  END DO

   IF(all_leaves_on==1) THEN

      p => pt%first

      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=p%coh%FPAR(i)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

END SUBROUTINE LIGHT_1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_2

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  real    :: help

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  !*** Calculation part ***!
  vStruct%cumLAI = 0.
  vStruct%Irel   = 0.

   ! cohort loop
   p => pt%first
   DO WHILE (ASSOCIATED(p))

       p%coh%FPAR = 0.
       p%coh%totFPAR = 0.

       p => p%next
   END DO        ! cohort loop

  ! Now calculate crown projection per tree and layer and
  ! the coverage sum over all layers
  CALL CROWN_PROJ

  ! now calculate coverage-area as fraction of the patchsize per tree and layer
  CALL COV_AREA

  vStruct(highest_layer)%Irel=1

  DO i = highest_layer, lowest_layer, -1

    p => pt%first

    help=0.

    vStruct(i)%cumLAI  = vStruct(i)%LA/kpatchsize + vStruct(i+1)%cumLAI

    DO WHILE (ASSOCIATED(p))

        ns=p%coh%species

         IF (p%coh%BG(i).ne.0.) THEN

             ! faction of absorbed light rel. to the light at the top of this layer
             ! the reference area is the whole patch (weighted by BG(i))!
             p%coh%FPAR(i)=(1-exp(-spar(ns)%pfext*p%coh%leafArea(i)/&
                                          kpatchsize/p%coh%BG(i)))*p%coh%BG(i)

             ! sum up the total absorbed fraction of this cohort,
             ! the total fraction of absorbed light in this layer
             ! is the fraction absorbed* fraction of light*BG
             ! the reference area is the whole patch!
             p%coh%totFPAR=p%coh%totFPAR+vStruct(i)%Irel*p%coh%FPAR(i)*&
                                (1+(0.5-vStruct(i)%Irel)*spar(ns)%fpar_mod/0.5)

             ! at first sum all the absorbed light fractions over the cohorts
             help=help+p%coh%FPAR(i)*p%coh%nTreeA

         ELSE

             p%coh%FPAR(i)=0.

         END IF

        p => p%next

    END DO

    ! then calculate the fraction of light which is available for the next layer
    vStruct(i-1)%Irel=vStruct(i)%Irel*(1-help)

  END DO

  ! Light intensitiy unto the ground
  DO i = lowest_layer - 2, 0, -1
      vStruct(i)%Irel=vStruct(i+1)%Irel
  END DO

   IF(all_leaves_on==1) THEN

      p => pt%first

      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=vStruct(i)%Irel*p%coh%FPAR(i)*(1+(0.5-vStruct(i)%Irel)*spar(ns)%fpar_mod/0.5)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

  ! total LAI is simply the value of cumLAI at the forest floor
  LAI = vStruct(lowest_layer)%cumLAI
  IF(lai>laimax) laimax=lai

END SUBROUTINE LIGHT_2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE L_3_COH_LOOP(i,j)

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

   INTEGER :: i, j   ! i= Schicht, j= Variante
   REAL    :: help

   p => pt%first

   ! cohort loop in layer i
   DO WHILE (ASSOCIATED(p))

       ns=p%coh%species

       IF((iday < p%coh%day_bb) .OR. (iday > spar(ns)%end_bb))  GOTO 1313

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           p%coh%FPAR(i)=1-exp(-spar(ns)%pfext*p%coh%leafArea(i)/&
                               kpatchsize/p%coh%BG(i))

           ! FPAR is related to the projection area and has to be modified
           ! by the same factor by that the projection area is being modified
           ! in case sumBG > patchsize
           p%coh%FPAR(i)=p%coh%FPAR(i)*MIN(kpatchsize/vStruct(i)%sumBG,1.)

           ! test wether the cohort is new, was there before or will not be
           ! represented in the next layer
           IF (i == p%coh%toplayer) THEN

               p%coh%Irel(i)=Irelpool(i)

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

           ELSE IF (i == p%coh%botlayer) THEN

               IF( j == 2 ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

               END IF

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

               ! The light which leaves the cohort is fed into the pool
               ! the light intensitiy is weighted by the overall BG of this cohort
               Irelpool(i-1)=(1/(p%coh%BG(i)*p%coh%nTreeA+BGpool(i)))*&
                             (p%coh%BG(i)*p%coh%nTreeA*p%coh%Irel(i-1)+BGpool(i)*Irelpool(i-1))

               ! BG of the pool available for the next layer increases
               BGpool(i)=BGpool(i)+p%coh%BG(i)*p%coh%nTreeA

           ELSE

               IF( j == 2 ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

               END IF

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

           END IF

       END IF   ! Layer test

1313   CONTINUE

       p => p%next

   END DO        ! cohort loop

END SUBROUTINE L_3_COH_LOOP

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_3

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  REAL    :: help

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  !*** Calculation part ***!
  vStruct%cumLAI = 0.
  Irelpool      = 0.
  BGpool        = 0.
  vStruct%Irel = 0.       ! test variable for the light balance in layers
  vStruct%radFrac = 0.    ! test variable for the light balance in layers

  ! cohort loop
  p => pt%first
  DO WHILE (ASSOCIATED(p))

      p%coh%FPAR = 0.
      p%coh%totFPAR = 0.
      p%coh%Irel = 0.

      p => p%next
  END DO        ! cohort loop

    ! Now calculate crown projection per tree and layer and
    ! the coverage sum over all layers
    CALL CROWN_PROJ

    ! now calculate coverage-area as fraction of the patchsize per tree and layer
    CALL COV_AREA

    ! -----------------------------------------------------------
    ! now calculate per tree and layer the effective LAI
    ! this gives the absorbed light per tree and layer
    ! this gives the total fraction absorbes light per tree
    ! further each tree and each layer has an individual light regime. The area
    ! which is not covered by trees is treated as a pool
    !
    ! reference area for the total fracation absorbed is the patch area

    ! above the canopy there is 100 % rel. light
    Irelpool(highest_layer)=1.
    ! the size of the pool is defined as the fraction of the patch
    ! which can potentially be used by new cohorts in the next layer.
    ! Therefore is is the patch-fraction which is free anyway plus the
    ! fraction coverd by cohorts that will not be present in the next layer
    ! this means, the light intensity Irelpool(i) is available on the
    ! area BGpool(i+1)
    BGpool(highest_layer+1)=1.

    DO i = highest_layer, lowest_layer, -1

        vStruct(i)%cumLAI  = vStruct(i)%LA/kpatchsize + vStruct(i+1)%cumLAI

        ! two cases:
        ! first case: sumBG increases in this layer or remains the same
        IF (vStruct(i+1)%sumBG<=vStruct(i)%sumBG) THEN

            ! three subcases:
            ! first subcase of 'sumBG increases': sumBG stays below patchsize
            ! ( no BG modification) or does not change
            IF ((vStruct(i+1)%sumBG.LT.kpatchsize.AND.vStruct(i)%sumBG.LE.kpatchsize).OR.&
                vStruct(i+1)%sumBG == vStruct(i)%sumBG) THEN

                ! At the beginning the light intensity of the pool remains the same
                ! but it will be updated when cohorts drop out
                Irelpool(i-1)=Irelpool(i)

                ! until there are cohorts dropping out
                BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)

                CALL L_3_COH_LOOP(i,1)

            ! second and third subcase of 'sumBG increases or remains the same'
            ! the BG's of the cohorts change because sumBG exceeds patchsize.
            ! second subcase: sumBG was < patchsize before
            ! third subcase: sumBG was > patchsize before
            ELSE

                ! BG and light intensitiy of the pool for the next(!) layer
                ! is 0 as long as there are no cohorts dropping out
                Irelpool(i-1)=0.
                BGpool(i)=0.

                p => pt%first

                ! cohort loop 1
                DO WHILE (ASSOCIATED(p))

                    ! calculate the new fraction covered by the pool
                    ! which is the old pool plus the fractions which are lost
                    ! by the old cohorts due to new BG's
                    ! this also changes the light intensity of the pool
                    ! This pool will all be used by the new cohorts
                    ! consider only cohorts that have been there before (i<toplayer)
                   IF (i<p%coh%toplayer.AND.i>=p%coh%botlayer .AND.&
                       iday >= p%coh%day_bb .AND. iday <= spar(p%coh%species)%end_bb) THEN

                        help=BGpool(i+1)+(p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA

                        Irelpool(i)=(1/help)*(Irelpool(i)*BGpool(i+1)+p%coh%Irel(i)*&
                                              (p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA)

                        BGpool(i+1)=help

                    END IF       ! layer test

                    p => p%next

                END DO   ! cohort loop1

                CALL L_3_COH_LOOP(i,1)

            END IF    ! subcases of 'sumBG increases

        ! second case: sumBG decreases
        ELSE

            ! two subcases
            ! first subcase of 'sumBG decrease': sumBG < patchsize before and after
            ! i.e. BG's do not change
            ! i.e. all projection area requirements can be fulfilled in the next layer
            IF (vStruct(i+1)%sumBG.LT.kpatchsize) THEN

                ! At the beginning the light intensity of the pool remains the same
                ! but it will be updated when cohorts drop out
                Irelpool(i-1)=Irelpool(i)

                ! until there are cohorts dropping out
                BGpool(i)=(kpatchsize-vStruct(i)%sumBG)/kpatchsize

                CALL L_3_COH_LOOP(i,1)

            ! second subcase of 'sumBG decrease': sumBG remains > patchsize or
            ! sumBG was > patchsize, i.e. BG's do change
            ELSE

                ! BG of the pool for the next layer as long as there are
                ! no cohorts dropping out
                BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)
                Irelpool(i-1)=Irelpool(i)

                CALL L_3_COH_LOOP(i,2)

            END IF            ! subcases

        END IF        ! three main cases

    END DO    ! end layer loop
    ! -----------------------------------------------------------
   IF(all_leaves_on==1) THEN

      p => pt%first
      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

    ! total LAI is simply the value of cumLAI at the lowest layer
    LAI = vStruct(lowest_layer)%cumLAI
    IF(lai>laimax) laimax=lai

    ! light intensitiy and free patch space unto the ground
    DO i = lowest_layer - 2, 0, -1
        Irelpool(i)=Irelpool(i+1)
        BGpool(i+1)=BGpool(i+2)
    END DO

END SUBROUTINE LIGHT_3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE L_4_COH_LOOP(i,j,beta,y)

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

   INTEGER :: i, j     ! i= layer, j= type
   REAL    :: y        ! potential shadow cast of a cohort layer [m]
   REAL    :: l        ! side length of a cohort layer [m]
   REAL    :: w        ! effective shadow cast of a cohort layer [m]
   REAL    :: helplai  ! LAI per layer and cohort
   REAL    :: help
   REAL    :: beta     ! sun inclination
   REAL    :: dropoutpool  ! relative area covered by cohort dropping out
   REAL    :: f1,f2,f3,f4,f5,f6,f7,f8 ! average fraction of absorbed radiation in different
                                      ! regions of the tree layer according to the 4C description paper
   REAL    :: k        ! extintion coefficient
   REAL    :: reqarea  ! area of the patch required for the shadow cast for all cohorts per layer

   reqarea=0.

   ! cohort loop
   p => pt%first
   DO WHILE (ASSOCIATED(p))

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           l = SQRT(p%coh%BG(i)*kpatchsize)

           reqarea = reqarea + l*y*p%coh%nTreeA

       END IF

       p => p%next
   END DO        ! cohort loop

   ! the size of the pool is defined as the fraction of
   ! the patch which is not covered by cohorts. This is the
   ! area covered by the sum of the 'shadows' of the cohorts,
   ! i.e. y's or rather w's + the area of cohorts dropping out in the next layer +
   ! the are that exeeds the maximal required area by the shadow-cast.
   ! This is updated in each layer
   ! w is the width of the shadow-cast of the cohorts that is maximal y.
   ! This maximal y also defines the maximal required area for all shadows
   ! 'reqarea' = required area
   ! When the maximal y cannot be satisfied, then this area is reduced by the
   ! relative share of the available space not covered by cohorts to the
   ! maximal required area for shadow cast

   IF( kpatchsize > vStruct(i)%sumBG ) THEN
     if (reqarea .gt. 1E-08) then
       w = y*(kpatchsize-vStruct(i)%sumBG)/reqarea
     else
       w = y*kpatchsize
     endif
   ELSE

       w = 0

   END IF

   BGpool(i)=0.
   dropoutpool=0

   p => pt%first

   ! cohort loop in layer i
   DO WHILE (ASSOCIATED(p))

       ns=p%coh%species

       IF((iday < p%coh%day_bb) .OR. (iday > spar(ns)%end_bb))  GOTO 1313

       k = spar(ns)%pfext

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           l = SQRT(p%coh%BG(i)*kpatchsize)

             if( p%coh%BG(i).ne.0) then

		         helplai=p%coh%leafArea(i)/kpatchsize/p%coh%BG(i)
		         if (helplai .le. 0.) then
		         continue
		         endif
		     else
			     helplai = 0.
			 end if

           IF (i == p%coh%toplayer) THEN

               p%coh%Irel(i)=Irelpool(i)

           ELSE IF( j == 2 .AND. i /= p%coh%toplayer ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

           END IF

           ! two main cases:
           ! first case : all light from the side comes from the pool
           ! second case : light from the side comes partially from the cohort itself

           IF( w >= y ) THEN

               ! subcases : 1.: light from the side of the layer
               ! does only leave at the bottom of the layer
               ! 2: light from the side does also leave on the other side
               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               IF( y <= l ) THEN

                   f1 = 1-exp(-k*helplai/SIN(beta))
                   if (helplai .lt. 1.E-6) then
                      f2 = 0.  
                   else
                      f2 = 1-SIN(beta)/(k*helplai)*f1
                      if (f2 .lt. 0.) then
                           continue
                           f2 = 0.   
                      endif
                   endif
                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ((l-y)*l*p%coh%Irel(i)*f1+& ! max. LAI
                                  ! exits layer at the side
                                  y*l*f2*p%coh%Irel(i)+&
                                  ! from the side to the next layer
                                  y*l*f2*Irelpool(i))

                   p%coh%FPAR(i)=p%coh%totFPAR

                   ! average light leaving the bottom of the cohort
                   p%coh%Irel(i-1)=(1/l)*&
                                   ! max. LAI
                                   ((l-y)*p%coh%Irel(i)*(1-f1)+&
                                   ! from the side to the next layer
                                   y*(1-f2)*Irelpool(i))

                   ! Light in the pool.
                   IF(i /= p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+y*l*p%coh%nTreeA)*&
                                    ! amount present in the pool
                                    (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                     ! exits layer at the side
                                     y*l*p%coh%nTreeA*(1-f2)*p%coh%Irel(i))

                       BGpool(i)=BGpool(i)+y*l*p%coh%nTreeA/kpatchsize

                   ELSE

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+(y+l)*l*p%coh%nTreeA)*&
                                     ! amount present in the pool
                                     (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                     ! exits layer at the side
                                     y*l*p%coh%nTreeA*(1-f2)*p%coh%Irel(i)+&
                                     ! from layer onto next layer
                                     l*l*p%coh%nTreeA*p%coh%Irel(i-1))

                       ! BG of the pool available for the next layer increases
                       BGpool(i)=BGpool(i)+p%coh%nTreeA*(y*l/kpatchsize+p%coh%BG(i))
                       dropoutpool=dropoutpool+p%coh%nTreeA*p%coh%BG(i)

                   END IF

               ! y > l
               ELSE

                   f3 = 1-exp(-k*helplai*l/(SIN(beta)*y))
                   f4 = 1-SIN(beta)*y/(l*k*helplai)*f3

                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ((y-l)*l*f3*Irelpool(i)+& ! red. max. LAI
                                   ! exits layer at the side
                                   l*l*f4*p%coh%Irel(i)+&
                                   ! from the side to next layer
                                   l*l*f4*Irelpool(i))

                   p%coh%FPAR(i)=p%coh%totFPAR

                   ! average light leaving the cohort
                   p%coh%Irel(i-1)=(1-f4)*Irelpool(i)

                   ! Light in the pool. Even when the area of the pool is
                   ! equal to zero, there is virtual light in the pool
                   ! which is used as light coming from the side
                   ! the area weighted mean over all y is calculated
                   IF(i /= p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+y*l*p%coh%nTreeA)*&
                                      ! amount present in pool 
                                      (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                      ! red. max. LAI
                                      (y-l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i)+&
                                      ! exits layer at side
                                      l*l*p%coh%nTreeA*(1-f4)*p%coh%Irel(i))

                       BGpool(i)=BGpool(i)+y*l*p%coh%nTreeA/kpatchsize

                   ELSE

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+(l+y)*l*p%coh%nTreeA)*&
                                       ! amount present in the pool
                                       (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                       ! red. max. LAI
                                       (y-l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i)+&
                                       ! exits layer at side
                                       l*l*p%coh%nTreeA*(1-f4)*p%coh%Irel(i)+&
                                       ! from layer to next layer
                                       l*l*p%coh%nTreeA*p%coh%Irel(i-1))

                       ! BG of the pool available for the next layer increases
                       BGpool(i)=BGpool(i)+p%coh%nTreeA*(y*l/kpatchsize+p%coh%BG(i))
                       dropoutpool=dropoutpool+p%coh%nTreeA*p%coh%BG(i)

                   END IF ! bottom layer or not

               END IF ! light entering sideways also leaving sideways or not

           ! second main case : light from the side comes partially from the
           ! cohort itself
           ELSE

               ! Exit, when average light from the side needs itself as input
               ! should not happen because this is taken care for in COV_AREA
               IF( y-w > w+l ) THEN

                 if (.not.flag_mult8910) then
                   CALL stop_mess(time,'FATAL EXCEPTION RAISED IN CANOPY LIGHT ROUTINE 4')
                   CALL error_mess(time,'Light leaving the side of cohort needs itself as input. Cohort No.',REAL(p%coh%ident))
                   CALL error_mess(time,'Try decreasing layer height dz or increasing average sun inclination.',0.)
                 endif  
                   STOP
               END IF

               ! subcases : 1.: light from the side of the layer
               ! does only leave at the bottom of the layer
               ! 2: light from the side does also leave on the other side but light from the top
               ! still goes into the pool
               ! 3. light from the side does also leave on the other side and light from the top
               ! is all used as input again
               ! totFPAR per patch! because the projection area changes totFPAR has to
               ! be related to the patch in each layer
               IF( y <= l ) THEN

                   IF( w /= 0 ) THEN
                       ! max LAI
                       f1 = 1-exp(-k*helplai/SIN(beta))
                       ! edge piece
                       f5 = 1+SIN(beta)*y/((y-w)*k*helplai)*(exp(-k*helplai*(y-w)/(SIN(beta)*y))-1)
                       ! red. LAI
                       f6 = 1+SIN(beta)*y/(w*k*helplai)*(1-f1-exp(-k*helplai*(y-w)/(SIN(beta)*y)))

                   ELSE
                       ! max LAI
                       f1 = 1-exp(-k*helplai/SIN(beta))
                       f5 = 1+SIN(beta)*y/((y-w)*k*helplai)*(exp(-k*helplai*(y-w)/(SIN(beta)*y))-1)
                       f6 = 0
                   END IF

                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ! enters from above into the pool
                                 (w*l*f6*p%coh%Irel(i)+&
                                 ! from above on own side 
                                 (y-w)*l*f5*p%coh%Irel(i)+&
                                 ! max. LAI
                                 (l-y)*l*f1*p%coh%Irel(i)+&
                                 ! from pool to next layer
                                 w*l*f6*Irelpool(i)+&
                                 ! from the side to the next layer
                                 (y-w)*l*(1-f5)*f5*p%coh%Irel(i))

                   p%coh%FPAR(i)=p%coh%totFPAR

                   ! average light leaving the bottom of the cohort
                   p%coh%Irel(i-1)=(1/l)*&
                                   ! max. LAI
                                   ((l-y)*(1-f1)*p%coh%Irel(i)+&
                                   ! from pool to next layer
                                   w*(1-f6)*Irelpool(i)+&
                                   ! from the sides to the next layer
                                   (y-w)*(1-f5)*(1-f5)*p%coh%Irel(i))

                   ! Light in the pool.
                   IF(i /= p%coh%botlayer .AND. w/=0) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+w*l*p%coh%nTreeA)*&
                                    ! present in the pool
                                    (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                     ! exits layer at the side
                                     w*l*p%coh%nTreeA*(1-f6)*p%coh%Irel(i))

                       BGpool(i)=BGpool(i)+w*l*p%coh%nTreeA/kpatchsize

                   ELSE IF(i == p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+(w+l)*l*p%coh%nTreeA)*&
                                     ! present in pool
                                     (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                     ! exits layer to the side
                                     w*l*p%coh%nTreeA*(1-f6)*p%coh%Irel(i)+&
                                     ! from layer to next layer
                                     l*l*p%coh%nTreeA*p%coh%Irel(i-1))

                       ! BG of the pool available for the next layer increases
                       BGpool(i)=BGpool(i)+p%coh%nTreeA*(w*l/kpatchsize+p%coh%BG(i))
                       dropoutpool=dropoutpool+p%coh%nTreeA*p%coh%BG(i)

                   END IF
               ! light from the top still goes into the pool.
               ! The case w=0 is no longer permissible
               ELSE IF(y > l .AND. w >= y-l) THEN

                   IF( w /= y-l ) THEN
                       f3 = 1-exp(-k*helplai*l/(SIN(beta)*y))
                       f5 = 1+SIN(beta)*y/((y-w)*k*helplai)*(exp(-k*helplai*(y-w)/(SIN(beta)*y))-1)
                       f7 = 1+SIN(beta)*y/((l-y+w)*k*helplai)*(exp(-k*helplai*l/(SIN(beta)*y))-&
                                                               exp(-k*helplai*(y-w)/(SIN(beta)*y)))
                   ELSE
                      f3 = 1-exp(-k*helplai*l/(SIN(beta)*y))
                      f5 = 1+SIN(beta)*y/((y-w)*k*helplai)*(exp(-k*helplai*(y-w)/(SIN(beta)*y))-1)
                      f7 = 0
                   END IF

                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ! enters pool from above
                                 ((l-y+w)*l*f7*p%coh%Irel(i)+&
                                 ! from above into own side
                                 (y-w)*l*f5*p%coh%Irel(i)+&
                                 ! red. max. LAI
                                 (y-l)*l*f3*Irelpool(i)+&
                                 ! from the side into the next layer
                                 (l-y+w)*l*f7*Irelpool(i)+&
                                 ! from the side into the next layer
                                 (y-w)*l*f5*(1-f5)*p%coh%Irel(i))

                   p%coh%FPAR(i)=p%coh%totFPAR

                   ! average light leaving the cohort
                   p%coh%Irel(i-1)=(1/l)*((l-y+w)*((1-f7)*Irelpool(i)+&
                                     (y-w)*(1-f5)*(1-f5)*p%coh%Irel(i)))


                   ! Light in the pool.
                   IF(i /= p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+w*l*p%coh%nTreeA)*&
                                       ! present in the pool
                                       (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                       ! exits from top to the side
                                       (l-y+w)*l*p%coh%nTreeA*(1-f7)*p%coh%Irel(i)+&
                                       ! from the side into the pool
                                       (y-l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i))

                       BGpool(i)=BGpool(i)+w*l*p%coh%nTreeA/kpatchsize

                   ELSE IF (i == p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+(l+w)*l*p%coh%nTreeA)*&
                                       ! present in the pool
                                       (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                       ! exits from the sides
                                       (l-y+w)*l*p%coh%nTreeA*(1-f7)*p%coh%Irel(i)+&
                                       ! enters from the sied into the pool
                                       (y-l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i)+&
                                       ! from layer to next layer
                                       l*l*p%coh%nTreeA*p%coh%Irel(i-1))

                       ! BG of the pool available for the next layer increases
                       BGpool(i)=BGpool(i)+p%coh%nTreeA*(w*l/kpatchsize+p%coh%BG(i))
                       dropoutpool=dropoutpool+p%coh%nTreeA*p%coh%BG(i)

                   END IF ! bottom layer or not

               ! light from the top still goes into the pool
               ELSE IF(y > l .AND. w < y-l) THEN

                       f3 = 1-exp(-k*helplai*l/(SIN(beta)*y))
                       f4 = 1-SIN(beta)*y/(l*k*helplai)*f3
                       f8 = 1/(y-w)*(l*f4+(y-w-l)*f3)

                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ! from above to own side
                                 (l*l*f4*p%coh%Irel(i)+&
                                 ! from side to the own side and into the pool
                                 y*l*f3*Irelpool(i)+&
                                 ! from the side to the next layer and into the pool
                                 l*f8*(1-f8)*(l*p%coh%Irel(i)+(y-w-l)*Irelpool(i)))

                   p%coh%FPAR(i)=p%coh%totFPAR

                   ! average light leaving the cohort
                   p%coh%Irel(i-1)=(1-f4)*(1-f8)*(l*p%coh%Irel(i)+(y-w-l)*Irelpool(i))

                   ! Light in the pool.
                   IF(i /= p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+w*l*p%coh%nTreeA)*&
                                       ! present in the pool
                                       (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                       ! from the side into the pool
                                       (2*w-y+l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i)+&
                                       (y-w-l)*l*p%coh%nTreeA*(1-f3)*(1-f8)*&
                                                              (l*p%coh%Irel(i)+(y-w-l)*Irelpool(i)))

                       BGpool(i)=BGpool(i)+w*l*p%coh%nTreeA/kpatchsize

                   ELSE IF (i == p%coh%botlayer) THEN

                       Irelpool(i-1)=1/(BGpool(i)*kpatchsize+(l+w)*l*p%coh%nTreeA)*&
                                       ! present in the pool
                                       (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                                       ! from the side into the pool
                                       (2*w-y+l)*l*p%coh%nTreeA*(1-f3)*Irelpool(i)+&
                                       (y-w-l)*l*p%coh%nTreeA*(1-f3)*(1-f8)*&
                                                              (l*p%coh%Irel(i)+(y-w-l)*Irelpool(i))+&
                                       ! from layer to next layer 
                                       l*l*p%coh%nTreeA*(1-f4)*(1-f8)*&
                                                              (l*p%coh%Irel(i)+(y-w-l)*Irelpool(i)))

                       ! BG of the pool available for the next layer increases
                       BGpool(i)=BGpool(i)+p%coh%nTreeA*(w*l/kpatchsize+p%coh%BG(i))
                       dropoutpool=dropoutpool+p%coh%nTreeA*p%coh%BG(i)

                   END IF ! bottom layer or not

               END IF ! light entering sideways also leaving sideways or not


           END IF  ! two main cases
       END IF

1313   CONTINUE
if (p%coh%FPAR(i) .lt. 0. .or. p%coh%totFPAR .lt. 0.) then
    continue
    p%coh%FPAR(i) = 0.    ! intercept negative radiation 
    p%coh%totFPAR = 0.
endif
       p => p%next

   END DO        ! cohort loop

   ! Treelayers are distributed on the patch such that their y's
   ! cover the free space as good as possible
   IF( w > y ) THEN

       Irelpool(i-1)=1/(kpatchsize*(1+dropoutpool)-vStruct(i)%sumBG)*&
                      (BGpool(i)*kpatchsize*Irelpool(i-1)+&
                       (kpatchsize-vStruct(i)%sumBG-(BGpool(i)-dropoutpool)*kpatchsize)*Irelpool(i))

       BGpool(i)=(kpatchsize-vStruct(i)%sumBG)/kpatchsize + dropoutpool
   END IF

END SUBROUTINE L_4_COH_LOOP

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_4

  !*** Declaration part ***!
  USE data_climate
  USE data_par
  USE data_species
  USE data_stand
  use data_site

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  REAL    :: help
  REAL    :: y         ! potential shadow cast of the stand [m]

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  !*** Calculation part ***!
  vStruct%cumLAI = 0.
  Irelpool      = 0.
  BGpool        = 0.
  vStruct%Irel = 0.       ! test variable for the balance in layers
  vStruct%radFrac = 0.    ! test variable for the balance in layers

  y = dz/100/TAN(beta)
   ! cohort loop
   p => pt%first

   DO WHILE (ASSOCIATED(p))

       p%coh%FPAR = 0.
       p%coh%totFPAR = 0.
       p%coh%Irel = 0.
       p => p%next
   END DO        ! cohort loop

if (time .eq. 8 .and. iday .eq. 134) then
continue
endif 
   ! Now calculate crown projection per tree and layer and
   ! the coverage sum over all layers
   CALL CROWN_PROJ

   ! now calculate coverage-area as fraction of the patchsize per tree and layer
   CALL COV_AREA

   ! -----------------------------------------------------------
   ! now calculate per tree and layer the effective LAI
   ! this gives the absorbed light per tree and layer
   ! this gives the total fraction absorbes light per tree
   ! further each tree and each layer has an individual light regime. The area
   ! which is not covered by trees is treated as a pool
   ! whose light is available for all new cohorts.
   ! reference area for the total fraction absorbed is the patch area.

   ! GBpool is exactly defined in subroutine L_4_COH_LOOP
   BGpool(highest_layer+1)=1.

   ! above the canopy there is 100 % rel. light
   Irelpool(highest_layer)=1.

   DO i = highest_layer, lowest_layer, -1

       vStruct(i)%cumLAI  = vStruct(i)%LA/kpatchsize + vStruct(i+1)%cumLAI

       ! two cases:
       ! first case: sumBG increases in this layer or remains the same
       IF (vStruct(i+1)%sumBG<=vStruct(i)%sumBG) THEN

           ! three subcases:
           ! first subcase of 'sumBG increases': sumBG stays below patchsize
           ! ( no BG modification) or does not change
           IF ((vStruct(i+1)%sumBG.LT.kpatchsize.AND.vStruct(i)%sumBG.LE.kpatchsize).OR.&
               vStruct(i+1)%sumBG == vStruct(i)%sumBG) THEN

               !until light model 4 restriction apply
               IF ( i <= lm3layer ) THEN

                  ! At the beginning the light intensity of the pool remains the same
                  ! but it will be updated when cohorts drop out
                  Irelpool(i-1)=Irelpool(i)

                  ! until there are cohorts dropping out
                  BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)

                  CALL L_3_COH_LOOP(i,1)

                  ! FPAR in light model 3 defined differently has
                  ! to be redefined here to cause no conflict in crown.f
                  p => pt%first
                  DO WHILE (ASSOCIATED(p))
                      p%coh%FPAR(i)=p%coh%totFPAR
                      p => p%next
                  END DO   ! cohort loop1
               ELSE
                   CALL L_4_COH_LOOP(i,1,beta,y)
               END IF

           ! second and third subcase of 'sumBG increases or remains the same'
           ! the BG's of the cohorts change because sumBG exceeds patchsize.
           ! second subcase: sumBG was < patchsize before
           ! third subcase: sumBG was > patchsize before
           ELSE

               p => pt%first

               ! cohort loop 1
               DO WHILE (ASSOCIATED(p))

                   ! calculate the new fraction covered by the pool
                   ! which is the old pool plus the fractions which are lost
                   ! by the old cohorts due to new BG's
                   ! this also changes the light intensity of the pool
                   ! consider only cohorts that have been there before (i<toplayer)
                   ! consider only cohorts that have leafed out already, otherwise
                   ! it may happen that help=0
                   IF (i<p%coh%toplayer.AND.i>=p%coh%botlayer .AND.&
                       iday >= p%coh%day_bb .AND. iday <= spar(p%coh%species)%end_bb) THEN

                       help=BGpool(i+1)+(p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA

                       if( help.ne.0) then
					       Irelpool(i)=(1/help)*(Irelpool(i)*BGpool(i+1)+p%coh%Irel(i)*&
                                                 (p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA)
                           BGpool(i+1)=help
					   end if

                   END IF       ! layer test

                   p => p%next

               END DO   ! cohort loop1

               !until light model 4 restriction apply
               IF ( i <= lm3layer ) THEN
                    CALL L_3_COH_LOOP(i,1)

                  ! FPAR in light model 3 defined differently has
                  ! to be redefined here to cause no conflict in crown.f
                  p => pt%first
                  DO WHILE (ASSOCIATED(p))
                      p%coh%FPAR(i)=p%coh%totFPAR
                      p => p%next
                  END DO   ! cohort loop1
               ELSE
                   CALL L_4_COH_LOOP(i,1,beta,y)
               END IF

           END IF    ! subcases of 'sumBG increases

       ! second case: sumBG decreases
       ELSE

           ! two subcases
           ! first subcase of 'sumBG decrease': sumBG < patchsize before and after
           ! i.e. BG's do not change
           ! i.e. all projection area requirements can be fulfilled in the next layer
           IF (vStruct(i+1)%sumBG.LT.kpatchsize) THEN

               !until light model 4 restriction apply
               IF ( i <= lm3layer ) THEN

                  ! At the beginning the light intensity of the pool remains the same
                  ! but it will be updated when cohorts drop out
                  Irelpool(i-1)=Irelpool(i)

                  ! until there are cohorts dropping out
                  BGpool(i)=(kpatchsize-vStruct(i)%sumBG)/kpatchsize

                  CALL L_3_COH_LOOP(i,1)

                  ! FPAR in light model 3 defined differently has
                  ! to be redefined here to cause no conflict in crown.f
                  p => pt%first
                  DO WHILE (ASSOCIATED(p))
                      p%coh%FPAR(i)=p%coh%totFPAR
                      p => p%next
                  END DO   ! cohort loop1
               ELSE
                   CALL L_4_COH_LOOP(i,1,beta,y)
               END IF

           ! second subcase of 'sumBG decrease': sumBG remains > patchsize or
           ! sumBG was > patchsize, i.e. BG's do increase
           ELSE

               !until light model 4 restriction apply
               IF ( i <= lm3layer ) THEN

                  ! BG of the pool for the next layer as long as there are
                  ! no cohorts dropping out
                  BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)
                  Irelpool(i-1)=Irelpool(i)

                  CALL L_3_COH_LOOP(i,2)

                  ! FPAR in light model 3 defined differently has
                  ! to be redefined here to cause no conflict in crown.f
                  p => pt%first
                  DO WHILE (ASSOCIATED(p))
                      p%coh%FPAR(i)=p%coh%totFPAR
                      p => p%next
                  END DO   ! cohort loop1
               ELSE
                   CALL L_4_COH_LOOP(i,2,beta,y)
               END IF

           END IF            ! subcases

       END IF        ! three main cases

   END DO    ! end layer loop
   ! -----------------------------------------------------------
   IF(all_leaves_on==1) THEN

      p => pt%first

      DO WHILE (ASSOCIATED(p))

          p%coh%bes = 0.
          DO i = highest_layer, lowest_layer, -1
             if(p%coh%totFPAR.ne.0) p%coh%antFPAR(i)=(p%coh%FPAR(i)-p%coh%FPAR(i+1))/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
       ! besetting here weighted with relative leaf area in layer, could also be done with nimber of layers
             IF((vstruct(i)%sumBG > kpatchsize) .and. (p%coh%t_leaf .gt. zero)) p%coh%bes = p%coh%bes + p%coh%leafarea(i)/p%coh%t_leaf*(vstruct(i)%sumBG/kpatchsize)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

   ! total LAI is simply the value of cumLAI at the lowest canopy layer
    LAI = vStruct(lowest_layer)%cumLAI
    IF(lai>laimax) laimax=lai

    ! light intensitiy and free patch space unto the ground
    DO i = lowest_layer - 2, 0, -1
        Irelpool(i)=Irelpool(i+1)
        BGpool(i+1)=BGpool(i+2)
    END DO

END SUBROUTINE LIGHT_4

END SUBROUTINE CANOPY

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! writes essential light paramerter into light.res1
! seperated to cohorts and layers
SUBROUTINE LIGHT_OUT_2

  use data_simul
  USE data_out
  USE data_stand
  USE data_species

  INTEGER:: i=0,j=0

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  ! Header
  write(unit_light,'(2A5,5A9)') 'YEAR ','layer ','  Coh1   ', &
      '  Coh2   ','  Coh3   ','  Coh4   ','...'

  p => pt%first

  WRITE(unit_light,'(i3,A)',ADVANCE='NO') time,'      '

  ! the crown cover area for cohorts
  DO WHILE (ASSOCIATED(p))

      WRITE(unit_light,'(F8.2)',ADVANCE='NO') p%coh%crown_area

      p => p%next

  END DO

  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') '-----------------------------------------------------------------------'

  SELECT CASE (flag_light)

  CASE(3,4)

      DO i = highest_layer, lowest_layer, -1

          IF(i.EQ.lm3layer) WRITE(unit_light,'(A)',ADVANCE='NO') 'ab hier LM3!'
          WRITE(unit_light,'(A,i3)',ADVANCE='NO') 'IREL  ',i

       ! relativ light intensity that hits layers and cohorts
          p => pt%first

          DO j=1, anz_coh

              IF (p%coh%Irel(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%Irel(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A,A7)',ADVANCE='NO') 'BG',' '

       ! cover degree per cohort and layer
          p => pt%first

          DO j=1, anz_coh

              IF (p%coh%BG(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%BG(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A,A5)',ADVANCE='NO') 'FPAR',' '

       ! the fraction absorbed by corhort and layer
          p => pt%first

          DO j=1, anz_coh

              IF (p%coh%FPAR(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%FPAR(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '
          WRITE(unit_light,'(A,F8.4)') 'BGpool in dieser schicht :', BGpool(i)
          WRITE(unit_light,'(A,F8.4)') 'relative Ueberdeckung in dieser Schicht :', vStruct(i)%sumBG/kpatchsize
          WRITE(unit_light,'(A,F8.4)') 'Summer der Ueberdeckungen :', BGpool(i)+vStruct(i)%sumBG/kpatchsize
          WRITE(unit_light,'(A)') ' '
          WRITE(unit_light,'(A,F8.4)') 'Rel. Licht unter dieser schicht :', VStruct(i)%Irel
          WRITE(unit_light,'(A,F8.4)') 'totFparsum bis zu dieser schicht :', VStruct(i)%radFrac
          WRITE(unit_light,'(A,F8.4)') ' Lichtbilanz : ', vStruct(i)%Irel+VStruct(i)%radFrac
          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A)') '-----------------------------------------------------------------------'

      END DO ! layers loop


  CASE(2)

      DO i = highest_layer, lowest_layer, -1

          WRITE(unit_light,'(A,i3)',ADVANCE='NO') 'Irel  ',i

       ! relative light intensity that hits the layer and cohorts
          DO j=1, anz_coh

              WRITE(unit_light,'(F8.4)',ADVANCE='NO') vStruct(i)%Irel

          END DO

          WRITE(unit_light,'(A)') ' '

       ! cover degree per cohort and layers
          p => pt%first

          WRITE(unit_light,'(A,A7)',ADVANCE='NO') 'BG',' '


          DO j=1, anz_coh

              IF (p%coh%BG(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%BG(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A,A5)',ADVANCE='NO') 'FPAR',' '

       ! fraction absorbed by cohort and layer
          p => pt%first

          DO j=1, anz_coh

              IF (p%coh%FPAR(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%FPAR(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A)') '-----------------------------------------------------------------------'

      END DO

  CASE(1)

      DO i = highest_layer, lowest_layer, -1

          WRITE(unit_light,'(A,i3)',ADVANCE='NO') 'IREL  ',i

       ! relative light inensity that hits layers and cohorts
          DO j=1, anz_coh

              WRITE(unit_light,'(F8.4)',ADVANCE='NO') vStruct(i)%Irel

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A,A5)',ADVANCE='NO') 'FPAR',' '

       ! fraction absirbed by cohort and layer
          p => pt%first

          DO j=1, anz_coh

              IF (p%coh%FPAR(i) == 0.) THEN

                  WRITE(unit_light,'(F8.2)',ADVANCE='NO') -99.99

              ELSE

                  WRITE(unit_light,'(F8.4)',ADVANCE='NO') p%coh%FPAR(i)

              END IF

              p => p%next

          END DO

          WRITE(unit_light,'(A)') ' '

          WRITE(unit_light,'(A)') '-----------------------------------------------------------------------'

      END DO

  END SELECT

  WRITE(unit_light,'(A,A2)',ADVANCE='NO') 'totFPAR',' '

  p => pt%first

  DO j=1, anz_coh

      WRITE(unit_light,'(F8.5)',ADVANCE='NO') p%coh%totFPAR

      p => p%next

  END DO

      WRITE(unit_light,'(A)') ' '

      WRITE(unit_light,'(A,F8.4)') 'Summe totFPAR : ',totFPARsum

  SELECT CASE(flag_light)

  CASE(3,4)

      WRITE(unit_light,'(A,F8.4)') 'Irel(lowest-1) : ', Irelpool(lowest_layer-1)

      WRITE(unit_light,'(A,F8.4)') ' Lichtbilanz : ', Irelpool(lowest_layer-1)+totFPARsum

  CASE(1,2)

      WRITE(unit_light,'(A,F8.4)') 'Irel(lowest-1) : ', vStruct(lowest_layer-1)%Irel

      WRITE(unit_light,'(A,F8.4)') ' Lichtbilanz : ', vStruct(lowest_layer-1)%Irel+totFPARsum

  END SELECT

  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') '------------------------------------------------------------------------------------'
  WRITE(unit_light,'(A)') ' '
  WRITE(unit_light,'(A)') ' '

END SUBROUTINE LIGHT_OUT_2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE CROWN_PROJ

  ! Now calculate crown projection per tree and layer and
  ! the coverage sum over all layers

  !*** Declaration part ***!
  USE data_par
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  real    :: help, help1

  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  vStruct%sumBG=0.

  p => pt%first

  DO WHILE (ASSOCIATED(p))
     ns=p%coh%species

   ! SMALL TREES OR GROUND VEGETATION
   IF (p%coh%height.lt.thr_height .or. ns .eq. nspec_tree+1) THEN
       p%coh%crown_area = p%coh%t_leaf                              ! small trees or ground vegetation
   ELSEIF (p%coh%species.eq.nspec_tree+2) then                      ! Case mistletoe
       p%coh%crown_area=pi*(real(p%coh%nTreeA)*0.000475)**(0.6666)        ! 1 big ball: volume = sum of mistletoe standard balls (10 years, pfiz 2000)
                                                                    ! V=4/3*Pi*r^3 , r= (3*V/4*PI)^1/3, (set V=n*4/3*Pi*512, with r=0.08 standard ball), r=(n*5.12*10-4)^1/3,A=pi*(n*5.12*10-4)^2/3
   ELSE
       !            Formel nach Biber 1996 S. 121, Kronenradius [dm]= a*DBH [cm]+b
       help1 = MIN(spar(ns)%crown_c,spar(ns)%crown_a*(p%coh%diam)+spar(ns)%crown_b)
       help=PI*(help1)**2
       ! adaptation of seedling crown projected area
       IF(p%coh%ca_ini.GT.help) THEN
          p%coh%crown_area=p%coh%ca_ini
       ELSE IF (p%coh%ca_ini.LT.help.AND.p%coh%diam == 0) THEN
          if(p%coh%height_ini.eq.137. .or. p%coh%height.eq.p%coh%height_ini) then
             p%coh%crown_area=p%coh%ca_ini
          else
                p%coh%crown_area=(p%coh%height-p%coh%height_ini)/(137.-p%coh%height_ini)*&
               (PI*(spar(ns)%crown_b)**2-p%coh%ca_ini)+p%coh%ca_ini
          end if
       ELSE
          p%coh%crown_area=help
       END IF
   END IF

	   if(p%coh%crown_area.lt.0) then
	          p%coh%crown_area = p%coh%ca_ini
	   end if

        DO i=p%coh%topLayer,p%coh%botLayer,-1

            vStruct(i)%sumBG=vStruct(i)%sumBG+p%coh%crown_area*p%coh%nTreeA

        END DO

      p => p%next

  END DO

END SUBROUTINE CROWN_PROJ
