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