!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* SR SEED_multi *!
!* *!
!* including SR/Function *!
!* function rtflsp (regula falsi solving equation) *!
!* function weight *!
!* function weight1 *!
!* *!
!* generation of a variety of seedling cohorts for *!
!* one seed number according to seedmass distribution *!
!* (for given mean value and standard deviation) *!
!* *!
!* 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 seed_multi(nseed,nsp)
USE data_species
use data_stand
use data_help
use data_par
use data_soil
use data_simul
IMPLICIT NONE
integer :: nseed, nseedha, nsclass , k, j, nr
integer,dimension(:),allocatable :: nsc
real, dimension(:), allocatable :: msc, &
shooth, &
nschelp
integer :: nsp
REAL :: shoot
REAL :: ms, msclass, x1,x2,xacc,shelp, nshelp,ntot,help
REAL :: troot2
real :: standdev
real :: rtflsp, weight
TYPE(cohort) ::tree_ini
external weight
external rtflsp
if(nseed.eq.0) return
standdev = spar(nsp)%seedsd*1000.
hnspec = nsp
ms = spar(nsp)%seedmass *1000. ! g ---> mg
nseedha = nseed
nshelp = nseedha/10000.
! calculation of seed class number
if(flag_reg.eq.3) then
nsclass = int(100.*nshelp**0.6)
else if(flag_reg.eq.30) then
nsclass = int(10.*nshelp**0.6)+1
end if
allocate(nsc(nsclass))
allocate(nschelp(nsclass))
allocate(msc(nsclass))
allocate(shooth(nsclass))
! seed weight and number of seeds per class
msclass = 6.*standdev/nsclass
ntot = 0
help = (1/(sqrt(2*pi)*standdev))
do k=1, nsclass
msc(k) = (ms - 3.*standdev) + msclass*(k-1)
nschelp(k) = help*exp(-((msc(k)-ms)**2)/(2*(standdev)**2))
ntot = ntot + nschelp(k)
end do
do k= 1,nsclass
nsc(k) = nint((nschelp(k)*nseedha/ntot) + 0.5)
end do
! calculation of shoot weight per seed class and initilization
do k = 1,nsclass
mschelp = msc(k)/1000000. ! mg ---> kg
x1 = 0.
x2 = 0.1
xacc=(1.0e-10)*(x1+x2)/2
! solve mass equation; determine root
shelp=rtflsp(weight,x1,x2,xacc)
shooth(k)= shelp
max_coh = max_coh + 1
call coh_initial (tree_ini)
tree_ini%ident = max_coh
tree_ini%species = nsp
tree_ini%ntreea = nsc(k)
tree_ini%nta = nsc(k)
shoot = shooth(k)
tree_ini%x_sap = shoot ! [kg]
shoot = shoot * 1000. ! [g]
tree_ini%x_fol= (spar(nsp)%seeda*(tree_ini%x_sap** spar(nsp)%seedb)) ![kg]
tree_ini%x_frt = tree_ini%x_fol ! [kg]
! Leder
tree_ini%x_hrt = 0.
tree_ini%med_sla = spar(nsp)%psla_min + spar(nsp)%psla_a*0.5
tree_ini%t_leaf = tree_ini%med_sla* tree_ini%x_fol ! [m-2]
tree_ini%ca_ini = tree_ini%t_leaf
tree_ini%crown_area = tree_ini%ca_ini
tree_ini%underst = 1
! tranformation of shoot biomass kg --> mg
if(nsp.ne.2)tree_ini%height = spar(nsp)%pheight1*(shoot*1000.)**spar(nsp)%pheight2 ! [cm] berechnet aus shoot biomass (mg)
! Leder
if(nsp.eq.2) tree_ini%height = 10**(spar(nsp)%pheight1+ spar(nsp)%pheight2*LOG10(shoot*1000.)+ &
spar(nsp)%pheight3*(LOG10(shoot*1000.))**2)
IF(nsc(k).ne.0.) then
IF (.not. associated(pt%first)) THEN
ALLOCATE (pt%first)
pt%first%coh = tree_ini
NULLIFY(pt%first%next)
! root distribution
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 = tree_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
anz_coh=anz_coh+1
END IF
end do
deallocate(nsc)
deallocate(nschelp)
deallocate(msc)
deallocate(shooth)
END SUBROUTINE seed_multi
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! weight: seed mass function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function weight (x)
use data_help
use data_species
implicit none
real :: x
real :: p1,p2, weight
p1 = spar(hnspec)%seeda
p2 = spar(hnspec)%seedb
weight = p1*2*(x**p2) + x - 0.7*mschelp
end function weight
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! weight1: coarse root mass function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function weight1 (x)
use data_help
use data_species
real :: x
real :: p1,p2
p1 = spar(hnspec)%seeda
p2 = spar(hnspec)%seedb
weight1 = p1*(x**p2) + x - mschelp
end function weight1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! rtflsp: regula falsi solving euation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
FUNCTION rtflsp(func,x1,x2,xacc)
INTEGER MAXIT
REAL rtflsp,x1,x2,xacc,func
EXTERNAL func
PARAMETER (MAXIT=30)
INTEGER j
REAL del,dx,f,fh,fl,swap,xh,xl
fl=func(x1)
fh=func(x2)
if(fl.lt.0.)then
xl=x1
xh=x2
else
xl=x2
xh=x1
swap=fl
fl=fh
fh=swap
endif
dx=xh-xl
do j=1,MAXIT
rtflsp=xl+dx*fl/(fl-fh)
f=func(rtflsp)
if(f.lt.0.) then
del=xl-rtflsp
xl=rtflsp
fl=f
else
del=xh-rtflsp
xh=rtflsp
fh=f
endif
dx=xh-xl
if(abs(del).lt.xacc.or.f.eq.0.)return
end do
END function rtflsp