!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*                      Interception                             *!
!*                                                               *!
!*   contains:                                                   *!
!*   INTERCEP                                                    *!
!*   INTERCEP_SVEG                                               *!
!*   INT_LAYER                                                   *!
!*   INT_COH_LOOP1                                               *!
!*   INT_COH_LOOP2                                               *!
!*   INT_COH_LOOP3                                               *!
!*                                                               *!
!*                  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 intercep

! Interception of the whole stand
! Stand variables are calculated in stand_balance

use data_climate
use data_inter
use data_evapo
use data_par
use data_simul
use data_species
use data_soil
use data_stand

implicit none

type(Coh_Obj), pointer :: p  ! pointer to cohort list
real aev_c, helplai, hxx, hsum, harea,  &
     cepmax, cepmax_can, cepmax_sveg,   &
     prec_eff, &        ! effective crown precipitation
     R_crown, &
     interc_c, &        ! interception per cohort
     pet_c              ! pet per cohort

! effective crown precipitation like Anders et al., 2002, S. 95
  prec_eff = prec * (1 + 0.13 * wind * (1-crown_area/kpatchsize))
  aev_i = 0.
  select case (flag_inth)

  case (0)   ! nach Jansson (SOIL)
!  Evaporation calculated at the start (==> interception is possible to be higher)
!    evaporation of intercepted water aev_i is limited by potential evaporation
   aev_c      = max(min(int_st_can, pet), 0.)               		
   int_st_can = max(int_st_can - aev_c, 0.)  ! interception storage from actual day

! Canopy interception 
    if (lai_can .gt. 0.) then
      cepmax_can = ceppot_can * lai_can     ! max. int. cap. of the whole stand
      if (airtemp .ge. temp_snow) then 			! frost conditions 
         lint_snow = .false.
         hxx = 0.
         if (cepmax_can .ge. int_st_can) hxx = cepmax_can-int_st_can
         interc_can  = min(hxx, prec)
      else
         lint_snow = .true.
         hxx = 0.
         if (2.*cepmax_can .ge. int_st_can) hxx = 2.*cepmax_can-int_st_can
         interc_can  = min(hxx, prec)
      endif
    else
      cepmax_can = ceppot_can * LAI_can     ! max. int. cap. of the whole stand, only canopy
      hxx = 0.
      if (cepmax_can .ge. int_st_can) hxx = cepmax_can-int_st_can
      interc_can = crown_area/kpatchsize * 0.15 * prec
      interc_can  = min(hxx, interc_can)
      aev_c      = 0.
    endif
   int_st_can = int_st_can + interc_can
 ! interception of ground vegetation
    if (flag_sveg .gt. 0) call intercep_sveg (aev_c)
 ! interception and interc.-evaporation of cohorts 
    call interc_coh (aev_c)
    aev_i = aev_i + aev_c 
	    
!......................................
  
  case (1)   ! interception for each cohort
             ! with distribution of precipit. over all canopy layers
       int_st_can  = 0.
       int_st_sveg = 0.
	   interc_can  = 0.
	   interc_sveg = 0.
	   aev_i       = 0.
       hsum        = 0.

     if (prec .gt. 0. .and. highest_layer .gt. 0) then
        call Int_layer
     else
       p => pt%first
       do while (associated(p))
          p%coh%interc = 0.
          p => p%next
       enddo  ! p (cohorts)          
	 endif

       p => pt%first
       do while (associated(p))
	     ns = p%coh%species
         if (all_leaves_on .eq. 0) then
            if((anz_tree.ne.0) .and. (pet .gt. 0.)) then
			    pet_c = pet  * p%coh%ntreea / anz_tree 
			else
			     pet_c = 0.
		    end if      
         else 
             if (flag_eva .eq. 2 .or. flag_eva .eq. 4) then
                pet_c = p%coh%demand
             else
                if((anz_tree.ne.0) .and. (pet .gt. 0.)) then
			        pet_c = pet * p%coh%rel_fol 
			    else
			         pet_c = 0.
		        end if      
                p%coh%demand = pet_c
             endif
         endif

         interc_c  = p%coh%interc
         select case (ns)  ! species

         case (1,12,13)   ! Fagus sylvatica 
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st
         
         case (2,10,15)   ! Picea abies ... Mistletoe
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, 2.*pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

         case (3,6,7,9)   ! Pinus sylvestris
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

         case (4,5,8,11) ! Quercus robur, Betula pendula 
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, 2.*pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

         case (14) ! Ground vegetation 
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_sveg = interc_sveg + interc_c
           int_st_sveg = int_st_sveg + p%coh%interc_st
                    
         end select

	       p%coh%aev_i= aev_c
		   aev_i      = aev_i + aev_c
		   p => p%next
       enddo  ! p (cohorts)  
	    
!......................................

  case (2)   ! interception for each cohort 
             ! with relativ part of precipit. accord. to foliage
       int_st_can  = 0.
       int_st_sveg = 0.
	   interc_can  = 0.
	   interc_sveg = 0.
	   aev_i       = 0.
       hsum        = 0.
       stem_flow   = 0.
       
       p => pt%first
       do while (associated(p))
	     ns = p%coh%species
         if (flag_eva .eq. 2 .or. flag_eva .eq. 4) then
            pet_c = p%coh%demand
         else
            pet_c = pet * p%coh%rel_fol
         endif

         select case (ns)  ! species

         case (1)   ! Fagus sylvatica 
           if ((iday .ge. p%coh%day_bb) .and. (iday .le. spar(ns)%end_bb)) then
             helplai   = p%coh%t_leaf/p%coh%crown_area
             cepmax    = spar(ns)%ceppot_spec * p%coh%rel_fol * helplai 
             if (airtemp .ge. temp_snow) then 			! frost conditions 
                hxx = 0.
                if (cepmax .ge. p%coh%interc_st) hxx = cepmax - p%coh%interc_st
                interc_c  = min(hxx, prec * p%coh%rel_fol)
                stem_flow  = stem_flow + 0.2 * (prec * p%coh%rel_fol - interc_c)
             else
                interc_c  = 0.35 * prec * p%coh%rel_fol
             endif
           else 
             interc_c  = 0.1 * prec * p%coh%rel_fol
             stem_flow = stem_flow + 0.16 * prec * p%coh%rel_fol
           endif
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, 2.*pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st
           stem_flow  = stem_flow + 0.16 * prec * p%coh%rel_fol
         
         case (2,10,15)   ! Picea abies   ... Mistletoe
           helplai   = p%coh%t_leaf/p%coh%crown_area
           cepmax    = spar(ns)%ceppot_spec * p%coh%rel_fol * helplai 
           if (airtemp .ge. temp_snow) then 			! frost conditions 
              hxx = 0.
              if (cepmax .ge. p%coh%interc_st) hxx = p%coh%interc_st
              interc_c  = min(cepmax-hxx, prec * p%coh%rel_fol)
           else
              interc_c  = 0.35 * prec * p%coh%rel_fol
           endif
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, 2.*pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

         case (3,6,7,9)   ! Pinus sylvestris
           helplai   = p%coh%t_leaf/p%coh%crown_area
           cepmax    = spar(ns)%ceppot_spec * p%coh%rel_fol * helplai 
           if (airtemp .ge. temp_snow) then 			! frost conditions 
              hxx = 0.
              if (cepmax .ge. p%coh%interc_st) hxx = p%coh%interc_st
              interc_c  = min(cepmax-hxx, prec * p%coh%rel_fol)
           else
              interc_c  = 0.35 * prec * p%coh%rel_fol
           endif
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

         case (4,5,8,11) ! Quercus robur, Betula pendula 
           if ((iday .ge. p%coh%day_bb) .and. (iday .le. spar(ns)%end_bb)) then
             helplai   = p%coh%t_leaf/p%coh%crown_area
             cepmax    = spar(ns)%ceppot_spec * p%coh%rel_fol * helplai 
             if (airtemp .ge. temp_snow) then 			! frost conditions 
                hxx = 0.
                if (cepmax .ge. p%coh%interc_st) hxx = p%coh%interc_st
                interc_c  = min(cepmax-hxx, prec * p%coh%rel_fol)
             else
                interc_c  = 0.35 * prec * p%coh%rel_fol
             endif
           else 
              interc_c = 0.05 * prec * p%coh%rel_fol
           endif
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, 2.*pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_can = interc_can + interc_c
           int_st_can = int_st_can + p%coh%interc_st

           case (14) ! Ground vegetation 
           if ((iday .ge. p%coh%day_bb) .and. (iday .le. spar(ns)%end_bb)) then
             helplai   = p%coh%t_leaf/p%coh%crown_area
             cepmax    = spar(ns)%ceppot_spec * p%coh%rel_fol * helplai 
             if (airtemp .ge. temp_snow) then 			! frost conditions 
                hxx = 0.
                if (cepmax .ge. p%coh%interc_st) hxx = p%coh%interc_st
                interc_c  = min(cepmax-hxx, prec * p%coh%rel_fol)
             else
                interc_c  = 0.35 * prec * p%coh%rel_fol
             endif
           else
              if (iday .eq. spar(ns)%end_bb+1) then
                interc_c = p%coh%interc_st
              else
                interc_c = 0.
              endif 
           endif
	       p%coh%interc_st = p%coh%interc_st + interc_c
		   aev_c           = min(p%coh%interc_st, pet_c)
	       p%coh%interc_st = max(p%coh%interc_st - aev_c, 0.)
           interc_sveg = interc_sveg + interc_c
           int_st_sveg = int_st_sveg + p%coh%interc_st
         end select
	       p%coh%aev_i= aev_c
		   aev_i      = aev_i + aev_c
		   p => p%next
       enddo  ! p (cohorts)  
	    
!......................................

  case (3)   ! interception pine like Anders et al., 2002, S. 95
      cepmax_can = ceppot_can * lai_can     ! max. int. cap. of the whole stand
      cepmax_can = 2.9   ! effect. crown storage capacity of pine according to Anders 
      R_crown    = 0.083 ! s/m  aerodyn. resistance of the crown of pine (Anders)
         if (cepmax_can .gt. prec_eff) then
            interc_can  = (crown_area/kpatchsize) * prec_eff 
         else
            interc_can  = cepmax_can + (prec_eff - cepmax_can) * wind * R_crown
            interc_can  = (crown_area/kpatchsize) * interc_can
         endif
      int_st_can = int_st_can + interc_can
      aev_c      = int_st_can    ! imediate total evaporation
      int_st_can = 0.            ! interception storage from actual day
	    
!......................................

   case (4)   ! from Refr.-Bez. (reference notation) (polynom.) for Level II, Brandenburg
      interc_can = 0.2 * prec
      int_st_can = int_st_can + interc_can
      ! evaporation of intercepted water aev_i is limited by potential evaporation
      aev_c      = min(int_st_can, pet)
      int_st_can = max(int_st_can - aev_c, 0.)  ! interception storage from actual day

   ! Interception of ground vegetation
      if (flag_sveg .gt. 0) call intercep_sveg (aev_c)

   ! interception and interc.-evaporation of cohorts 
      call interc_coh (aev_c)
      aev_i = aev_i + aev_c 
                     		
   case (5)   ! 35% of precipitation (for spruce)
      interc_can = 0.3 * prec
      int_st_can = int_st_can + interc_can
      ! evaporation of intercepted water aev_i is limited by potential evaporation
      aev_c      = min(int_st_can, pet)
      int_st_can = max(int_st_can - aev_c, 0.)  ! interception storage from actual day

   ! interception of ground vegetation
      if (flag_sveg .gt. 0) call intercep_sveg (aev_c)

   ! interception and interc.-evaporation of cohorts 
      call interc_coh (aev_c)
      aev_i = aev_i + aev_c 
        
   case (6)   ! no interception
      interc_can  = 0.
      aev_c       = 0.
      interc_sveg = 0.       
   end select

if (flag_dayout .eq. 3) then
 write(666,*) 'day, prec, prec_eff: ', iday, prec, prec_eff
endif

! cumul. interc.
 int_cum_can  = int_cum_can + interc_can
 int_cum_sveg = int_cum_sveg + interc_sveg
 if(flag_dayout.eq.3) write(1414,*) iday, aev_i
END	subroutine intercep

!**************************************************************

SUBROUTINE Int_layer

! Interception per canopy layer
! calculation for each cohort in subroutine int_coh_loop1 (rain)
! and int_coh_loop3 (int_coh_loop2 old) for snow

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

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  REAL    :: intlay, itest  ! interception per layer
  REAL    :: help
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  !*** Calculation part ***!
  precpool      = 0.
  itest = 0.
  intlay = 0.

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

      p%coh%intcap = 0.
      p%coh%interc = 0.
      p%coh%prel = 0.

      p => p%next
  END DO        ! cohort loop

   ! above the canopy there is 100 % precipitation
    precpool(highest_layer) = prec

    if (airtemp .ge. temp_snow) then 			! frost conditions 
      
      lint_snow = .false.
      do i = highest_layer, lowest_layer, -1

                intlay = 0.
				CALL int_coh_loop1(i,intlay)
         ! Assum.: all layers are above eachother, that means precip. is reduc. layer by layer due to interception.
                precpool(i-1) = precpool(i) - intlay 
             itest = itest + intlay  
      enddo    ! end layer loop

    else
      
      lint_snow = .true.
      CALL int_coh_loop3(intlay)
      endif    ! airtemp

    ! stand precipitation unto the ground
    DO i = lowest_layer - 2, 0, -1
        precpool(i)=precpool(i+1)
    END DO
 itest = 0.

END SUBROUTINE Int_layer

!**************************************************************

SUBROUTINE int_coh_loop1(i,intlay)

! interception for each canopy layer of each cohort 

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

  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list
  INTEGER :: i, itop ! layer
  REAL    :: intlay  ! interception per layer
  REAL    :: interc_c, &  ! interception per cohort
             cepcap       ! Int.-Kapaz. fuer diese Variante reduzieren
  REAL    :: help, hxx

  interc_c = 0.
  p => pt%first

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

    ns=p%coh%species

	IF ((iday >= p%coh%day_bb) .AND. (iday <= spar(p%coh%species)%end_bb))  then	    
       IF (i <= p%coh%toplayer .AND. i >= p%coh%botlayer) THEN
         p%coh%prel(i) = precpool(i) * p%coh%BG(i) * p%coh%nTreeA

         select case (ns)  ! species

         case (1)   ! Fagus sylvatica 
           if (p%coh%t_leaf .gt. 0.) then
              cepcap          = spar(ns)%ceppot_spec * 0.5
              p%coh%intcap(i) = cepcap * p%coh%leafArea(i) * p%coh%rel_fol / &
                              (kpatchsize * p%coh%BG(i))
               ! intcap 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%intcap(i)=p%coh%intcap(i) * MIN(kpatchsize/vStruct(i)%sumBG, 1.)
               ! interc per patch! Since the projection area changes interc has to
               ! be related to the patch in each layer
              hxx = 0.
              if (p%coh%intcap(i) .ge. p%coh%interc_st/dz) hxx = p%coh%interc_st/dz   ! interc storage spead across all layers
              interc_c = min(p%coh%prel(i), p%coh%intcap(i)-hxx) 
            else
              interc_c = 0.1 * p%coh%prel(i)
			endif
         
         case (2,10,15)   ! Picea abies ... mistletoe
            cepcap          = spar(ns)%ceppot_spec * 0.5
            p%coh%intcap(i) = cepcap * p%coh%leafArea(i) * p%coh%rel_fol / &
                             (kpatchsize * p%coh%BG(i))
             ! intcap 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%intcap(i)=p%coh%intcap(i) * MIN(kpatchsize/vStruct(i)%sumBG, 1.)
             ! interc per patch! Since the projection area changes interc has to
             ! be related to the patch in each layer
            hxx = 0.
            if (p%coh%intcap(i) .ge. p%coh%interc_st/dz) hxx = p%coh%interc_st/dz   ! interc storage spead across all layers
            interc_c  = min(p%coh%prel(i), p%coh%intcap(i)-hxx) 
         
         case (3,6,7,9)   ! Pinus sylvestris
            cepcap          = spar(ns)%ceppot_spec * 0.5
            p%coh%intcap(i) = cepcap * p%coh%leafArea(i) * p%coh%rel_fol / &
                             (kpatchsize * p%coh%BG(i))
             ! intcap 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%intcap(i)=p%coh%intcap(i) * MIN(kpatchsize/vStruct(i)%sumBG, 1.)
             ! interc per patch! Since the projection area changes interc has to
             ! be related to the patch in each layer
            hxx = 0.
            if (p%coh%intcap(i) .ge. p%coh%interc_st/dz) hxx = p%coh%interc_st/dz   ! interc storage spead across all layers
            interc_c = min(p%coh%prel(i), p%coh%intcap(i)-hxx) 

         case (4,5,8,11) ! Quercus robur, Betula pendula 
           if (p%coh%t_leaf .gt. 0.) then
              cepcap          = spar(ns)%ceppot_spec * 0.5
              p%coh%intcap(i) = cepcap * p%coh%leafArea(i) * p%coh%rel_fol / &
                              (kpatchsize * p%coh%BG(i))
               ! intcap 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%intcap(i)=p%coh%intcap(i) * MIN(kpatchsize/vStruct(i)%sumBG, 1.)
               ! interc per patch! Since the projection area changes interc has to
               ! be related to the patch in each layer
              hxx = 0.
              if (p%coh%intcap(i) .ge. p%coh%interc_st/dz) hxx = p%coh%interc_st/dz   ! interc storage spead across all layers
              interc_c = min(p%coh%prel(i), p%coh%intcap(i)-hxx) 
            else
              interc_c = 0.1 * p%coh%prel(i)
			endif
			        
         case (14) ! Ground vegetation 
           if (p%coh%t_leaf .gt. 0.) then
              cepcap          = spar(ns)%ceppot_spec * 0.5
              p%coh%intcap(i) = cepcap * p%coh%leafArea(i) * p%coh%rel_fol / &
                              (kpatchsize * p%coh%BG(i))
               ! intcap 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%intcap(i)=p%coh%intcap(i) * MIN(kpatchsize/vStruct(i)%sumBG, 1.)
               ! interc per patch! Since the projection area changes interc has to
               ! be related to the patch in each layer
              hxx = 0.
              if (p%coh%intcap(i) .ge. p%coh%interc_st/dz) hxx = p%coh%interc_st/dz   ! interc storage spead across all layers
              interc_c = min(p%coh%prel(i), p%coh%intcap(i)-hxx)
            else
              interc_c = 0.0
			endif
			        
         end select

	   ENDIF  ! i - layer
    ELSE
       IF (i == p%coh%toplayer) THEN
         itop = i
         if(cover.ne.0) p%coh%prel(itop) = precpool(i) * p%coh%nTreeA *p%coh%crown_area/crown_area

         select case (ns)  ! species

         case (1)   ! Fagus sylvatica    p%coh%x_tb
           interc_c  = 0.2 * p%coh%prel(itop)

         case (2,10,15)   ! Picea abies ... Mistletoe
           interc_c      = 0.1 * p%coh%prel(itop) 

         case (3,6,7,9)   ! Pinus sylvestris
           interc_c      = 0.1 * p%coh%prel(itop) 

         case (4,5,8,11) ! Quercus robur, Betula pendula 
              interc_c  = 0.1 * p%coh%prel(itop)

         case (14) ! Ground vegetation 
              interc_c  = 0.
                    
         end select
       ENDIF  ! i - layer
    END IF    ! iday
    
    if (interc_c .le. 1E-15) interc_c = 0. 
    p%coh%interc = p%coh%interc + interc_c
    intlay = intlay + interc_c
	interc_c = 0.
	p => p%next
  END DO        ! cohort loop

END SUBROUTINE int_coh_loop1

!**************************************************************

SUBROUTINE int_coh_loop2(i,intlay)

! snow interception for each canopy layer of each cohort 
  !*** Declaration part ***!
  USE data_simul
  USE data_soil
  USE data_species
  USE data_stand
  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list
  INTEGER :: i       ! layer
  REAL    :: intlay  ! interception per layer
  REAL    :: interc_c, &  ! interception per cohort
             cepcap       ! Int.-Kapaz. fuer diese Variante reduzieren
  REAL    :: help, hxx

  interc_c = 0.
  p => pt%first

  ! cohort loop in layer i
  DO WHILE (ASSOCIATED(p))
    ns=p%coh%species
       IF (i <= p%coh%toplayer .AND. i >= p%coh%botlayer) THEN

         select case (ns)  ! species

         case (1)   ! Fagus sylvatica    
          if(cover.ne.0)  p%coh%prel(i) = precpool(i) * p%coh%nTreeA *p%coh%crown_area/(kpatchsize*cover)
           if (p%coh%t_leaf .gt. 0.) then
              interc_c  = 0.35 * p%coh%prel(i) 
           else 
              interc_c  = 0.1 * p%coh%prel(i)
           endif

         case (2,10,15)   ! Picea abies... Mistletoe
           p%coh%prel(i) = precpool(i) * p%coh%BG(i) * p%coh%nTreeA
           interc_c      = 0.35 * p%coh%prel(i) 

         case (3,6,7,9)   ! Pinus sylvestris
           p%coh%prel(i) = precpool(i) * p%coh%BG(i) * p%coh%nTreeA
           interc_c      = 0.35 * p%coh%prel(i) 

         case (4,5,8,11) ! Quercus robur, Betula pendula 
           p%coh%prel(i) = precpool(i) * p%coh%nTreeA *p%coh%crown_area/kpatchsize
           if (p%coh%t_leaf .gt. 0.) then
              interc_c  = 0.35 * p%coh%prel(i) 
           else 
              interc_c  = 0.1 * p%coh%prel(i)
           endif

         case (14) ! Ground vegetation 
           if (p%coh%t_leaf .gt. 0.) then
              p%coh%prel(i) = precpool(i) * p%coh%BG(i) * p%coh%nTreeA
              interc_c      = 0.35 * p%coh%prel(i) 
           else 
              interc_c  = 0.
           endif
         end select
         if (interc_c .le. 1E-15) interc_c = 0. 
         p%coh%interc = p%coh%interc + interc_c
	END IF

1313   CONTINUE
    intlay = intlay + interc_c
	interc_c = 0.
	p => p%next
  END DO        ! cohort loop

END SUBROUTINE int_coh_loop2

!**************************************************************

SUBROUTINE int_coh_loop3(intlay)

! snow interception for each cohort 
  !*** Declaration part ***!
  USE data_climate
  USE data_simul
  USE data_soil
  USE data_species
  USE data_stand
  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list
  INTEGER :: itop    ! toplayer
  REAL    :: intlay  ! canopy interception 
  REAL    :: interc_c, &  ! interception per cohort
             cepcap       ! Int.-Kapaz. fuer diese Variante reduzieren
  REAL    :: help, hxx
  real test_prel

  test_prel = 0.
  interc_c = 0.
  p => pt%first

  ! cohort loop 
  DO WHILE (ASSOCIATED(p))
    ns=p%coh%species
       itop    = p%coh%toplayer
        if(cover.ne.0)  p%coh%prel(itop) = prec * p%coh%nTreeA *p%coh%crown_area/crown_area
        test_prel = test_prel + p%coh%prel(itop)

         select case (ns)  ! species

         case (1)   ! Fagus sylvatica    p%coh%x_tb
           if (p%coh%t_leaf .gt. 0.) then
              interc_c  = 0.35 * p%coh%prel(itop) 
           else 
              interc_c  = 0.1 * p%coh%prel(itop)
           endif

         case (2,10,15)   ! Picea abies, Douglas Fir,  Mistletoe (better: nspec_tree+2)
           interc_c      = 0.35 * p%coh%prel(itop) 

         case (3,6,7,9)   ! Pinus sylvestris, P. contorta, P. ponder. P. halep.
           interc_c      = 0.6 * p%coh%prel(itop) 

         case (4,5,8,11) ! Quercus robur, Betula pendula, Populus, Robinia
           if (p%coh%t_leaf .gt. 0.) then
              interc_c  = 0.35 * p%coh%prel(itop) 
           else 
              interc_c  = 0.1 * p%coh%prel(itop)
           endif

         case (14) ! Ground vegetation 
           if (p%coh%t_leaf .gt. 0.) then
              interc_c      = 0.35 * p%coh%prel(itop) 
           else 
              interc_c  = 0.
           endif                    
         end select
         if (interc_c .le. 1E-15) interc_c = 0. 
         p%coh%interc = p%coh%interc + interc_c

1313   CONTINUE
    intlay = intlay + interc_c
	interc_c = 0.
	p => p%next
  END DO        ! cohort loop
continue

END SUBROUTINE int_coh_loop3

!**************************************************************

SUBROUTINE intercep_sveg (aev_c)

! Interception of ground vegetation

use data_climate
use data_inter
use data_evapo
use data_par
use data_species
use data_stand

implicit none

real aev_c, &    ! canopy interception evaporation      
     hxx,   &
     cepmax_sveg

        cepmax_sveg = ceppot_sveg * lai_sveg ! max. int. cap. of the whole stand
        if (airtemp .ge. temp_snow) then 			 ! frost conditions 
           hxx = 0.
           if (cepmax_sveg .ge. int_st_sveg) hxx = cepmax_sveg-int_st_sveg
           interc_sveg  = min(hxx, prec-interc_can)
        else
           interc_sveg  = 0.35 * (prec-interc_can)
        endif
        int_st_sveg = int_st_sveg + interc_sveg
       ! evaporation of intercepted water aev_i is limited by potential evaporation
        aev_i  = min(int_st_sveg, pet-aev_c)               		
        int_st_sveg = max(int_st_sveg - aev_i, 0.)  ! interception storage from actual day

END SUBROUTINE intercep_sveg

!**************************************************************

SUBROUTINE interc_coh (aev_c)

! Interception of ground vegetation

use data_climate
use data_inter
use data_evapo
use data_species
use data_stand

implicit none

type(Coh_Obj), pointer :: p  ! pointer to cohort list
integer ns
real aev_c, &    ! canopy interception evaporation      
     cepmax_sveg
    p => pt%first
    do while (associated(p))
      ns = p%coh%species
      if (ns .le. nspec_tree .OR. ns .eq. nspec_tree+2) then
       ! trees and mistletoe
          p%coh%interc_st = int_st_can * p%coh%rel_fol        
          p%coh%aev_i= aev_c * p%coh%rel_fol
      else
       ! ground vegetation
          p%coh%interc_st = int_st_sveg * p%coh%rel_fol        
          p%coh%aev_i= aev_i * p%coh%rel_fol
      endif

      p => p%next
    enddo  ! p (cohorts)  

END SUBROUTINE interc_coh
