!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Stand initialisation *!
!* *!
!* CONTAINS SUBROUTINES : *!
!* PREPARE_STAND *!
!* internal subroutines: *!
!* SLA_INI *!
!* *!
!* CALC_INT *!
!* CALC_WEIBLA *!
!* READ_STAND (treeunit) *!
!* COH_INITIAL (coh) *!
!* CREATE_MISTLETOE *!
!* CREATE_SOILVEG *!
!* *!
!* CONTAiNS FUNCTIONS : *!
!* SURVAGE *!
!* *!
!* 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 prepare_stand
USE data_simul
USE data_site
USE data_stand
USE data_species
use data_climate
use data_par
USE data_manag
IMPLICIT NONE
CHARACTER :: a, with_storage
CHARACTER(30) :: text
CHARACTER(50) :: test_stand_id
INTEGER :: ios,treeunit
LOGICAL :: exs, lstin
INTEGER :: help_ip, test_vf
REAL :: test_patchsize, xx
REAL help_height_top ! auxiliary var. for setting mistletoe height at uppermost crown layer
INTEGER which_cohort
INTEGER nr_infect_trees
INTEGER nr_mist_per_tree
INTEGER i
TYPE(Coh_Obj), Pointer :: p ! pointer to cohort list
IF(site_nr==1) THEN
help_ip=site_nr
ELSE
help_ip=ip
END IF
pt = neu()
anz_coh=0
max_coh=0
ios = -1
nr_mist_per_tree=0
IF(flag_stand>0) then
exs = .false.
stand_id = standid(help_ip)
! reading stand information from treefile
inquire (File = treefile(help_ip), exist = exs)
IF((exs .eqv. .false.) .or. (flag_stand==2)) then
IF(exs .eqv. .false.) write(*,*) ' Stand initialization file not exists!'
IF(flag_stand==2) write(*,*)' Stand initialization with new file'
write(*,'(A)',advance='no') ' Creating new file (y/n): '
READ *, a
IF(a.eq.'y'.or. a.eq.'Y') CALL initia
! planting of small trees
if(flag_reg.eq.20) then
call planting
flag_reg=100
end if
flag_stand=1
exs=.true.
ENDIF
! read values from treefile
IF (exs.eqv. .true.) then
treeunit=getunit()
OPEN(treeunit,file=treefile(help_ip),action='read', pad='YES') !!!! NSC implementation for ini-file
READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage
backspace treeunit
if(with_storage .eq. 'C') then !flag_dis .eq. 2 .and.
READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage,test_vf, test_patchsize
call read_stand_with_nsc (treeunit)
CLOSE(treeunit)
kpatchsize = test_patchsize
anz_coh = max_coh
coh_ident_max = anz_coh
else
READ(treeunit,'(I1,F12.0)',iostat=ios) test_vf, test_patchsize
! write(8888,*) ip, test_vf, flag_volfunc
if(flag_multi.ne.4 .or. (flag_multi.eq.4.and.ip.eq.1) .or. (flag_multi.eq.8.and.ip.eq.1)) then
IF(test_vf.NE.flag_volfunc) THEN
if (.not.flag_mult8910) then
CALL error_mess(time,"volume function in sim-file and the one used for initialisation do not match",REAL(flag_volfunc))
CALL error_mess(time,"volume function (flag_volfunc) is set to",REAL(test_vf))
endif
flag_volfunc = test_vf
ENDIF
endif
IF(test_patchsize .GT. 0.) THEN
lmulti = .FALSE.
IF(test_patchsize.NE.kpatchsize) THEN
if (.not.flag_mult8910) then
CALL error_mess(time,"patch size in sim-file and the one used for initialisation do not match",kpatchsize)
CALL error_mess(time,"value in ini-file",test_patchsize)
CALL error_mess(time,"value in sim-file",kpatchsize)
endif
kpatchsize = test_patchsize
ENDIF
ELSE
lmulti = .TRUE.
ENDIF
do
READ(treeunit,'(A)',iostat=ios) a
IF (a .ne. '!') exit
end do
backspace treeunit
! generation of mistletoe cohort; mistletoe cohort need to be generated BEFORE tree cohorts as otherwise the light model becomes messy
if (flag_dis.eq.1) then
do i= 1, dis_row_nr
if (dis_type(i) .eq. 'M') then
if (flag_mistle.eq.0) then !set #of mist. only once
print *,"!! Note, implementation of mistletoe is restricted to trees of Pinus sylvestris"
nr_mist_per_tree = dis_rel(i)
flag_mistle=1 ! flag indicating mistletoes
call create_mistletoe ! initialisation of Mistletoe
endif
anz_coh = max_coh
endif
enddo
endif
lstin = .TRUE.
if(flag_multi.eq.4 .or. flag_multi.eq.8) stand_id = standid(help_ip)
do while (lstin)
IF (lmulti) THEN
read(treeunit,*,iostat=ios) test_stand_id, test_patchsize,text
IF (ios .lt. 0) then
if (.not.flag_mult8910) then
CALL error_mess(time,"stand identificator not found"//adjustl(stand_id)//"ip No.",real(help_ip))
write (*,*) '*** PREPSTAND: program aborted'
write (*,*) ' stand identificator',stand_id,' not found'
write (*,'(A, 2x,A)') ' in initialisation file',treefile(help_ip)
endif
flag_end = 2
return
ENDIF
IF (test_stand_id .ne. stand_id) THEN
read (treeunit,*) xx
do while (xx .gt. -90.0)
read (treeunit,*) xx
enddo ! xx
ELSE
lstin = .FALSE.
kpatchsize = test_patchsize
call read_stand (treeunit)
END IF ! stand_id
ELSE
lstin = .FALSE.
call read_stand (treeunit)
END IF ! lmulti
end do ! lstin
CLOSE(treeunit)
anz_coh = max_coh
coh_ident_max = anz_coh
endif !w/o storage
ENDIF !exs.eqv. .true.
END IF !if stand >0
!if treefile not exists and not created:
IF(ios .ne. 0 .or. exs .eqv. .false.)THEN
if (.not.flag_mult8910) PRINT *,' >>> No Stand Initialization possible '
flag_stand=0
END IF
! Setting of height and number of mistletoe
if (flag_mistle.ne.0) then
help_height_top=1.
p=>pt%first
DO WHILE (ASSOCIATED(p))
if (p%coh%species.eq.3 .AND. p%coh%height.gt.help_height_top) then !only on Pinus
help_height_top=p%coh%height
which_cohort=p%coh%ident
nr_infect_trees=p%coh%nTreeA
end if
p=>p%next
end do
p=>pt%first
DO WHILE (ASSOCIATED(p))
if (p%coh%species.eq.nspec_tree+2) then
p%coh%height = help_height_top !upper crown
p%coh%x_hbole = p%coh%height-50. !lower crown
p%coh%nTreeA = nr_infect_trees*nr_mist_per_tree !number of mistletoes
end if
if (p%coh%ident.eq.which_cohort) then !mark uppermost tree cohort with flag mistletoe
p%coh%mistletoe=1
end if
p=>p%next
end do
end if ! end set height/number of mistletoe
! Soil Vegetation
if (flag_sveg .gt. 0) then
call create_soilveg ! initialisation of ground vegetation
anz_coh = max_coh
endif
IF(flag_stand>0) CALL sla_ini
IF(flag_stand>0) CALL stand_bal_spec
CALL calc_int
CALL calc_weibla
if(flag_mg.ne.33) call overstorey
contains
SUBROUTINE sla_ini
USE data_stand
USE data_species
IMPLICIT NONE
TYPE(Coh_Obj), Pointer :: p ! pointer to cohort list
p => pt%first
DO WHILE (ASSOCIATED(p))
ns=p%coh%species
p%coh%med_sla=spar(ns)%psla_min+spar(ns)%psla_a*0.5
p%coh%t_leaf = p%coh%med_sla * p%coh%x_fol
p =>p%next
END DO
end subroutine sla_ini
end subroutine prepare_stand
!*************************************************************************
subroutine calc_int ! calculation of intrinsic mortality rate
use data_species
implicit none
INTEGER j
do j=1,nspecies
spar(j)%intr = -log(0.01)/spar(j)%max_age
end do
end subroutine calc_int
!*************************************************************************
subroutine calc_weibla
! calculation of parameter lamda for Weibull-distribution of sress mortality
use data_species
implicit none
INTEGER j
REAL survage
do j=1,nspecies
spar(j)%weibla = -log(0.01)/(survage(j)**weibal)
end do
end subroutine calc_weibla
!*************************************************************************
REAL function survage(ispec)
! calculation of survival time per species depending on shade tolerance class stol
use data_species
implicit none
INTEGER :: ispec
IF(spar(ispec)%stol.eq.1) survage=20.
IF (spar(ispec)%stol.eq.2) survage=40.
IF (spar(ispec)%stol.eq.3) survage=60.
IF (spar(ispec)%stol.eq.4) survage=80.
IF (spar(ispec)%stol.eq.5) survage=100.
end function
!*************************************************************************
SUBROUTINE read_stand (treeunit)
! Read of stand initialisation
USE data_par
USE data_simul
USE data_species
USE data_stand
IMPLICIT NONE
TYPE(cohort) :: coh_ini
REAL :: hdquo ! help variable for stress initilization
INTEGER :: ios,treeunit, loc, i
logical :: treegroup_decid
integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)
do
call coh_initial (coh_ini)
READ(treeunit,'(5f12.5,2f10.0,i7, f10.0,i7, f9.5, f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam
IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit
coh_ini%nTreeD = 0.
coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)
!NSC Speicher initialisieren
!hier neue version auslesen
treegroup_decid = .False.
do i = 1, 5
if (decidous(i) .eq. coh_ini%species) then
treegroup_decid = .True.
exit
endif
end do
If (treegroup_decid .eq. .True.) then
coh_ini%x_nsc_sap = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
coh_ini%x_nsc_tb = coh_ini%x_tb * decid_tb_allo * cpart
coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
coh_ini%x_nsc_crt = coh_ini%x_crt * decid_crt_allo * cpart
coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
endif
If (treegroup_decid .eq. .False.) then
coh_ini%x_nsc_sap = coh_ini%x_sap * conif_sap_allo * cpart
coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
coh_ini%x_nsc_tb = coh_ini%x_tb * conif_tb_allo * cpart
coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
coh_ini%x_nsc_crt = coh_ini%x_crt * conif_crt_allo * cpart
coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
endif
coh_ini%ident = max_coh + 1
coh_ini%Fmax = coh_ini%x_fol
coh_ini%x_health = 0
coh_ini%x_hsap = 0.
ns = coh_ini%species
coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol ! kg * mg/g --> g
if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
else
coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
end if
! Stress calculation
IF (coh_ini%diam.ne. 0.) THEN
hdquo = coh_ini%height/ (coh_ini%diam*100)
IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
coh_ini%x_stress = coh_ini%x_age/2
ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
coh_ini%x_stress = coh_ini%x_age*3./7.
ELSE
coh_ini%x_stress = 0.
END IF
ELSE
coh_ini%x_stress = 0.
END IF ! coh_ini
coh_ini%x_stress = 0.
coh_ini%nta = coh_ini%nTreeA
IF (.not. associated(pt%first)) THEN
max_coh = 0
allocate(pt%first)
pt%first%coh = coh_ini
nullify(pt%first%next)
ELSE
allocate(zeig)
zeig%coh = coh_ini
zeig%next => pt%first
pt%first => zeig
END IF
max_coh = max_coh + 1
enddo
END SUBROUTINE read_stand
!*************************************************************************
SUBROUTINE read_stand_with_nsc (treeunit)
! Read of stand initialisation
USE data_par
USE data_simul
USE data_species
USE data_stand
IMPLICIT NONE
TYPE(cohort) :: coh_ini
REAL :: hdquo ! help variable for stress initilization
INTEGER :: ios,treeunit, loc, i
logical :: treegroup_decid
integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)
character :: a
do
READ(treeunit,'(A)',iostat=ios) a
IF (a .ne. '!') exit
end do
backspace treeunit
do
call coh_initial (coh_ini)
! READ(treeunit,'(5f12.5,2f10.0,i7, f7.0,i7, 2f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
READ(treeunit,'(5f12.5,2f10.0,i7,f7.0,i7, 5f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam, coh_ini%x_nsc_tb, coh_ini%x_nsc_crt, coh_ini%x_nsc_sap
IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit
coh_ini%nTreeD = 0.
coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)
!hier nur den NSC Max-Speicher initialisieren
!hier neue version auslesen
treegroup_decid = .False.
do i = 1, 5
if (decidous(i) .eq. coh_ini%species) then
treegroup_decid = .True.
exit
endif
end do
If (treegroup_decid .eq. .True.) then
coh_ini%x_nsc_sap_max = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
coh_ini%x_nsc_tb_max = coh_ini%x_tb * decid_tb_allo * cpart
coh_ini%x_nsc_crt_max = coh_ini%x_crt * decid_crt_allo * cpart
endif
If (treegroup_decid .eq. .False.) then
coh_ini%x_nsc_sap_max = coh_ini%x_sap * conif_sap_allo * cpart
coh_ini%x_nsc_tb_max = coh_ini%x_tb * conif_tb_allo * cpart
coh_ini%x_nsc_crt_max = coh_ini%x_crt * conif_crt_allo * cpart
endif
! IF(coh_ini%species==3.and.coh_ini%height.le.900) then
! coh_ini%x_stress =6
! ELSE IF(coh_ini%species==1.and.coh_ini%height.le.2500) then
! coh_ini%x_stress =25
! ELSE
! coh_ini%x_stress = 0
! ENDIF
coh_ini%ident = max_coh + 1
coh_ini%Fmax = coh_ini%x_fol
coh_ini%x_health = 0
coh_ini%x_hsap = 0.
ns = coh_ini%species
coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol ! kg * mg/g --> g
! calculate diameter at breast height and pipe length by call of subroutine in partitio
if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
! if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.137.) then
CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
else
coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
end if
! Stress calculation
IF (coh_ini%diam.ne. 0.) THEN
hdquo = coh_ini%height/ (coh_ini%diam*100)
IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
coh_ini%x_stress = coh_ini%x_age/2
ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
coh_ini%x_stress = coh_ini%x_age*3./7.
ELSE
coh_ini%x_stress = 0.
END IF
ELSE
coh_ini%x_stress = 0.
END IF ! coh_ini
! provisorisch stress auf Null setzen
coh_ini%x_stress = 0.
coh_ini%nta = coh_ini%nTreeA
IF (.not. associated(pt%first)) THEN
max_coh = 0
allocate(pt%first)
pt%first%coh = coh_ini
nullify(pt%first%next)
ELSE
allocate(zeig)
zeig%coh = coh_ini
zeig%next => pt%first
pt%first => zeig
END IF
max_coh = max_coh + 1
enddo
END SUBROUTINE read_stand_with_nsc
!*************************************************************************
SUBROUTINE coh_initial (coh_ini)
USE data_simul
USE data_soil
USE data_stand
USE data_species
IMPLICIT NONE
TYPE(cohort) :: coh_ini
coh_ini%nTreeA = 0.
coh_ini%nTreeD = 0.
coh_ini%nTreeM = 0.
coh_ini%nTreet = 0.
coh_ini%nta = 0.
coh_ini%mistletoe = 0
coh_ini%x_age = 0.
coh_ini%x_fol = 0.
coh_ini%x_sap = 0.
coh_ini%x_frt = 0.
coh_ini%x_hrt = 0.
coh_ini%x_crt = 0.
coh_ini%x_tb = 0.
coh_ini%x_hsap = 0.
coh_ini%x_hbole= 0.
coh_ini%x_Ahb = 0.
coh_ini%x_stress = 0
coh_ini%x_health = 0
coh_ini%bes = 0.
coh_ini%med_sla = 0.
coh_ini%Fmax = 0
coh_ini%totBio = 0.
coh_ini%Dbio = 0.
coh_ini%height = 0.
coh_ini%deltaB = 0.
coh_ini%dcrb = 0.
coh_ini%diam = 0.
coh_ini%assi = 0.
coh_ini%LUE = 0.
coh_ini%resp = 0.
coh_ini%netAss = 0.
coh_ini%NPP = 0.
coh_ini%weekNPP = 0.
coh_ini%NPPpool = 0.
coh_ini%t_Leaf = 0.
coh_ini%geff = 0.
coh_ini%Asapw = 0.
coh_ini%crown_area = 0.
coh_ini%BG = 0.
coh_ini%leafArea = 0.
coh_ini%sleafArea = 0.
coh_ini%FPAR = 0.
coh_ini%antFPAR = 0.
coh_ini%Irel = 0.
coh_ini%totFPAR = 0
coh_ini%IrelCan = 0
coh_ini%botLayer = 0
coh_ini%topLayer = 0
coh_ini%survp = 0.
coh_ini%rel_fol = 0.
coh_ini%gfol = 0.
coh_ini%gfrt = 0.
coh_ini%gsap = 0.
coh_ini%sfol = 0.
coh_ini%sfrt = 0.
coh_ini%ssap = 0.
coh_ini%grossass = 0.
coh_ini%maintres = 0.
coh_ini%respsap = 0.
coh_ini%respfrt = 0.
coh_ini%respbr = 0.
coh_ini%height_ini = 0.
coh_ini%ca_ini = 0.
coh_ini%rel_dbh_cl = 0
coh_ini%underst = 0
coh_ini%fol_inc = 0.
coh_ini%fol_inc_old = 0.
coh_ini%bio_inc = 0.
coh_ini%stem_inc = 0.
coh_ini%frt_inc = 0.
coh_ini%notViable = .FALSE.
coh_ini%intcap = 0.
coh_ini%prel = 0.
coh_ini%interc = 0.
coh_ini%prelCan = 0.
coh_ini%interc_st= 0.
coh_ini%aev_i = 0.
coh_ini%demand = 0.
coh_ini%supply = 0.
coh_ini%watuptc = 0.
coh_ini%gp = 0.
coh_ini%drIndd = 0.
coh_ini%drIndPS = 0.
coh_ini%drIndAl = 0.
coh_ini%nDaysGr = 0
coh_ini%isGrSDay = .false.
coh_ini%litC_fol = 0.
coh_ini%litC_fold = 0.
coh_ini%litN_fol = 0.
coh_ini%litN_fold = 0.
coh_ini%litC_frt = 0.
coh_ini%litC_frtd = 0.
coh_ini%litN_frt = 0.
coh_ini%litN_frtd = 0.
coh_ini%litC_stem = 0.
coh_ini%litN_stem = 0.
coh_ini%litC_tb = 0.
coh_ini%litC_crt = 0.
coh_ini%litC_tbcd = 0.
coh_ini%litN_tb = 0.
coh_ini%litN_crt = 0.
coh_ini%litN_tbcd = 0.
coh_ini%Nuptc_c = 0.
coh_ini%Nuptc_d = 0.
coh_ini%Ndemc_d = 0.
coh_ini%RedNc = 1.
coh_ini%N_pool = 0.
coh_ini%N_fol = 0.
coh_ini%wat_mg = 0. ! soley forflag_wred=9
coh_ini%nroot = 0
coh_ini%shelter = 0
coh_ini%day_bb = 0
coh_ini%x_nsc_sap = 0.
coh_ini%x_nsc_tb = 0.
coh_ini%x_nsc_crt = 0.
coh_ini%x_nsc_sap_max = 0.
coh_ini%x_nsc_tb_max = 0.
coh_ini%x_nsc_crt_max = 0.
if (coh_ini%species .ne. nspec_tree+2) then ! no root allocation for mistletoe
allocate (coh_ini%frtrel(nlay))
allocate (coh_ini%frtrelc(nlay))
if (flag_wred .eq. 9) then
allocate (coh_ini%rld(nlay))
coh_ini%rld = 0.
endif
allocate (coh_ini%rooteff(nlay))
coh_ini%frtrel = 0.
coh_ini%rooteff = 0.
end if ! end exclude mistletoe
END SUBROUTINE coh_initial
!*************************************************************************
SUBROUTINE create_mistletoe
USE data_plant
USE data_simul
USE data_species
USE data_stand
USE data_climate
USE data_soil
USE data_species
USE data_par
IMPLICIT NONE
TYPE(cohort) :: coh_ini
real :: help_height_top, help_height_bot
REAL, EXTERNAL :: fi_lf, dfi_lf, ddfi_lf
! initialising of cohort of mistletoe
call coh_initial (coh_ini)
! set mistletoe here to 20 m height, will be changed after, when cohorts of trees will be initialised
help_height_top=2000
help_height_bot=help_height_top-50
! following values are from sample calcul. of 10 year old V.austr. from Pfiz 2010
coh_ini%ident = max_coh + 1
coh_ini%species = nspec_tree+2 ! Species = species after all tree species and ground veg.
coh_ini%nTreeA = 1 ! #of mistletoes, to be read-in in management file
coh_ini%nTreeD = 0 ! dead trees
coh_ini%nta = coh_ini%nTreeA ! alive trees internal calc.
coh_ini%x_age = 10
coh_ini%x_fol = mistletoe_x_fol ! fol biomass per tree [kg DW/tree], 1 Viscum (10years) see Pfiz 2010
coh_ini%x_sap = 0. ! set near-zero for partitioning
coh_ini%x_frt = 0. ! set near-zero for partitioning
coh_ini%height = help_height_top ! highest_layer ! highest_layer of all cohorts
coh_ini%x_hbole = help_height_bot !
coh_ini%med_sla = 0. ! average cohort specific leaf area [m2/kg] is being calculated internal
coh_ini%Fmax = 0 ! anual change of leaf biomass, for now: now change
coh_ini%crown_area = 0.0189 ! max. projected crown area (m2) per individuum, calculated from Pfiz 2010
coh_ini%t_leaf = coh_ini%med_sla* coh_ini%x_fol !leaf area per tree [m2] !
coh_ini%day_bb = 1 ! evergreen
! no partitioning of NPP into stem/leaf etc.
! no root allocation
allocate(zeig)
zeig%coh = coh_ini
zeig%next => pt%first
pt%first => zeig
max_coh = max_coh + 1
END SUBROUTINE create_mistletoe
!*************************************************************************
SUBROUTINE create_soilveg
! Read of stand initialisation
USE data_plant
USE data_simul
USE data_species
USE data_stand
USE data_climate
USE data_soil
IMPLICIT NONE
TYPE(cohort) :: coh_ini
real :: lai_help, irel_help, FRsum
integer :: age_stand, nr, j
integer :: flag_SV_allo, rnum
real :: troot2
REAL, EXTERNAL :: fi_lf, dfi_lf, ddfi_lf
age_stand = 0
lai_help = 0.
irel_help = 0.
call wclas(waldtyp)
zeig=>pt%first
DO WHILE (ASSOCIATED(zeig))
ns = zeig%coh%species
lai_help = lai_help + zeig%coh%ntreea*zeig%coh%x_fol* spar(ns)%psla_min
age_stand = MAX(zeig%coh%x_age,age_stand)
zeig=>zeig%next
end do
IF((flag_stand==0 .or. age_stand .le. 5) .AND. flag_sveg ==2) THEN
NPP_est = 10.
ELSE if(age_stand.le.5) then
if(ns.eq.4) then
NPP_est = 5
else
NPP_est = 10.
end if
ELSE if(flag_reg.ne.0) then
NPP_est = 10
ELSE
lai_help = lai_help/kpatchsize
irel_help = exp(-0.5*lai_help)
if( svar(nspec_tree+1)%RedN .lt.0.) then
NPP_est = irel_help * med_rad1 * 365./100. *0.5
else
NPP_est = irel_help * med_rad1 * 365./100. *0.5 * svar(nspec_tree+1)%RedN
end if
ENDIF
call coh_initial (coh_ini)
coh_ini%species = nspec_tree+1 ! numbre of species determined automatically
ns = coh_ini%species
flag_SV_allo=1
IF(flag_SV_allo==0) THEN
! the parameters pdiam in the species.par file are used for allocation fractions
coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize
ELSE
FRsum=0.8*NPP_est/1000. ! start value as fraction of NPP in kg DM m-2
CALL newt (FRsum, fi_lf, dfi_lf, ddfi_lf, 0.001, 100, rnum)
IF(rnum==-1) THEN
if (.not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort: ',real(ns))
coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize
ELSE
coh_ini%x_sap = (ksi*FRsum**kappa)*kpatchsize
coh_ini%x_fol = (FRsum/2.)*kpatchsize
coh_ini%x_frt = (FRsum/2.)*kpatchsize
ENDIF
ENDIF
coh_ini%height = 60.
coh_ini%x_age = 1
coh_ini%nTreeA = 1
coh_ini%ident = max_coh + 1
coh_ini%Fmax = coh_ini%x_fol
coh_ini%med_sla = spar(coh_ini%species)%psla_min + spar(coh_ini%species)%psla_a*irel_help
coh_ini%t_leaf = coh_ini%med_sla* coh_ini%x_fol ! [m2]
coh_ini%nta = coh_ini%nTreeA
coh_ini%ca_ini = kpatchsize
coh_ini%day_bb = 100 ! assumption budding on 8.April
! root allocation
IF (.not. associated(pt%first)) THEN
max_coh = 0
allocate(pt%first)
pt%first%coh = coh_ini
nullify(pt%first%next)
call root_depth (1, pt%first%coh%species, pt%first%coh%x_age, pt%first%coh%height, pt%first%coh%x_frt, pt%first%coh%x_crt, nr, troot2, pt%first%coh%x_rdpt, pt%first%coh%nroot)
pt%first%coh%nroot = nr
do j=1,nr
pt%first%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
pt%first%coh%rooteff = 0. ! layers with no roots
enddo
ELSE
allocate(zeig)
zeig%coh = coh_ini
zeig%next => pt%first
pt%first => zeig
call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
zeig%coh%nroot = nr
do j=1,nr
zeig%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
zeig%coh%rooteff = 0. ! layers with no roots
enddo
END IF
max_coh = max_coh + 1
END SUBROUTINE create_soilveg
!*************************************************************************
!***************************!
! FUNCTION fi_lf *!
!***************************!
REAL FUNCTION fi_lf(x)
USE data_stand
USE data_plant
USE data_species
REAL :: x
fi_lf = spar(nspec_tree+1)%pss*ksi*x**kappa + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2.*x - NPP_est/1000.
END ! FUNCTION fi_lf
!***************************!
! FUNCTION dfi_lf *!
!***************************!
REAL FUNCTION dfi_lf(x)
USE data_stand
USE data_plant
USE data_species
REAL :: x
dfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*x**(kappa-1.) + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2.
END ! FUNCTION dfi_lf
!***************************!
! FUNCTION ddfi_lf *!
!***************************!
REAL FUNCTION ddfi_lf(x)
USE data_stand
USE data_plant
USE data_species
REAL :: x
ddfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*(kappa-1.)*x**(kappa-2.)
END ! FUNCTION ddfi_lf