!*****************************************************************!
!* *!
!* 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