Skip to content
Snippets Groups Projects
Select Git revision
  • dev-mahnken
  • cherry-pick-c40afffb
  • cherry-pick-f1d9f851
  • master default protected
4 results
Blame Permalink
Forked from 4C / FORESEE
206 commits behind the upstream repository.
stand_regen.f 18.87 KiB 
!*****************************************************************!
!*                                                               *!
!*             4C (FORESEE) Simulation Model                     *!
!*                                                               *!
!*                                                               *!
!*   subroutine for regeneration                                 *!
!*   including the SR:                                           *!
!*   - gener_seed                                                *!
!*   - seed_ini                                                  *!
!*   - simseed                                                   *!
!*   - growth_seed                                               *!
!*   - mort_seed                                                 *!
!*                                                               *!
!*                  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 stand_regen

  USE data_simul
  USE data_stand

  IMPLICIT NONE

 flag_standup = 2   

 CALL mort_seed

 CALL gener_seed

END SUBROUTINE stand_regen


SUBROUTINE gener_seed

USE data_stand
USE data_species
USE data_simul
use data_out
USE data_plant
USE data_soil

IMPLICIT NONE
real        :: seedla               ! leaf area of all seedling cohorts
real        :: laiseed              ! lai            ----"------
integer     :: nseed                ! number of generated seeds
real      :: redseed
integer        :: i
integer, dimension(5)  :: agemin, seedpot
real, dimension(5,3)   :: latg
real        :: help,help1, help2
real        :: pequal
integer     :: hlayer
integer     :: flag_reg_help
TYPE(coh_obj), POINTER :: p
DATA latg /1.,0.3,0.1,1.,0.1,1.,0.9,0.5,1.,0.5,1.,1.,0.9,1.,0.9/

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' gener_seed'

flag_reg_help = 0
agemin = 0
seedpot = 0
seedla = 0.
help1 = 0.
! calculation of leafarea of all seedling cohorts

SELECT CASE (flag_reg)
! according to FORGRA (Vincent Kint)
   CASE (30)
   call random_number(pequal)
   
   DO i= 1, nspecies
   nseed = 0
   p =>pt%first
    DO WHILE (ASSOCIATED(p))
        if(p%coh%species.eq.i) then
         if (i.eq.1) then
            agemin(i) = 50 + (30* pequal)
            call random_number(pequal)
            seedpot(i) = 810*pequal
         else if(i.eq.3) then
            agemin(i) = 15 + (45* pequal)
            call random_number(pequal)
            seedpot(i) = 1000*pequal
         else if(i.eq.4) then
            agemin(i) = 15 +(35* pequal)
            call random_number(pequal)
            seedpot(i) = 1125*pequal
         else if(i.eq.5) then
            agemin(i) = 10 +(10* pequal)
            call random_number(pequal)
            seedpot(i) = 8750*pequal
         end if
         if(p%coh%x_age.ge. agemin(i).and.p%coh%diam.gt.0.) then
              nseed = nseed + seedpot(i)*(p%coh%ntreem + p%coh%ntreea)
          end if
        end if
       p => p%next
    END DO  ! cohort

     help2 = irelpool_ll
	 if(help2.lt.0) help2 =0
     if(help2.eq. 0.) then
         redseed = 0.
     else if( help2.gt.0. .and. help2.le.latg(i,1)) then
         redseed = help2*latg(i,1)/0.4
     else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
         redseed = help2*latg(i,2)/0.6
     else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
         redseed = help2*latg(i,3)/0.8
     else if(help2.gt.latg(i,3)) then
         redseed = help2* latg(i,3)
     end if

     nseed = nseed * redseed


! for birch 1 year old saplings are used
     if (i.eq.5) then
       numplant(i) = nseed
       flag_reg = 14
       if(nseed.ne.0) call planting
       flag_reg= 0
     else
       call seed_multi(nseed,i)
     end if

   END DO   ! species

   CASE(1,2,3)

p =>pt%first
  DO WHILE (ASSOCIATED(p))
      if(p%coh%height.lt.thr_height) then

          seedla = seedla + p%coh%t_leaf*p%coh%ntreea
          help1 = help1 + p%coh%x_fol*p%coh%ntreea
      end if
      p => p%next

  END DO

! calculation LAI of lowest_layer
   laiseed=seedla/kpatchsize

DO i=1,nspecies


   IF (spar(i)%regflag.eq.1) THEN
       CALL simseed(i,nseed)
        IF(laiseed.lt.1) THEN
! reduction of seedling number nseed depending on light module and free space in the lowest_layer
          SELECT CASE (flag_light)

              CASE(1)
                        IF(flag_reg.ne.3) THEN
                           CALL seed_ini(nseed,i)
                        ELSE
                          CALL seed_multi(nseed,i)
                        END IF
              CASE (2)
                  

                if (anz_coh.eq. 0) then
                  if(time.eq.1) then
                    hlayer = 0
                  else
                    hlayer = 1
                  end if
                 else
                   hlayer = lowest_layer -1
                end if

                 help =  vstruct(hlayer)%Irel
                 if (help.lt.0.05) help = 0
                 IF(help.eq.0) THEN
                        nseed = 0
                 ELSE
                        nseed = nseed*help
                        IF(flag_reg.ne.3) THEN
                           CALL seed_ini(nseed,i)
                        ELSE
                          CALL seed_multi(nseed,i)
                        END IF
                 END IF

              CASE(3)
                   redseed= bgpool_ll
                   nseed = nseed*redseed
                    IF(flag_reg.ne.3) THEN
                       CALL seed_ini(nseed,i)
                    ELSE
                       CALL seed_multi(nseed,i)
                    END IF

             CASE(4)
                if(i.gt.5) then

                    redseed= irelpool_ll
                else
 ! according to FORGRA, not for all species (i=1,5)
                 help2 = irelpool_ll
                  if(help2.lt. 0.01) then
                      redseed = 0.
 
                  else if( help2.gt.0.01 .and. help2.le.latg(i,1)) then

                      redseed = help2*latg(i,1)/0.4
                 else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
                      redseed = help2*latg(i,2)/0.6
                 else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
                       redseed = help2*latg(i,3)/0.8
                 else if(help2.gt.latg(i,3)) then
                      redseed = help2* latg(i,3)

                 end if
                end if
				 nseed = redseed * nseed     

                    IF(flag_reg.ne.3) THEN
                       CALL seed_ini(nseed,i)
                    ELSE
                         if (i.eq.5) then
                           numplant(i) = nseed
						   flag_reg_help = flag_reg
                           flag_reg = 14
                           if(nseed.ne.0) call planting
					       flag_reg = flag_reg_help
                         else
					      CALL seed_multi(nseed,i)
						 end if
                    END IF
             END SELECT
        ELSE
           nseed = 0.
        END IF
   ELSE
       nseed=0.
   END IF
   
END DO

END SELECT   ! flag_reg


END subroutine gener_seed

SUBROUTINE simseed(specnum,nseed)
USE data_species
use data_simul
use data_stand
IMPLICIT NONE
REAL	        :: pequal
INTEGER	:: nseed,specnum
REAL	        :: seedmax

! calculation of max. seedrate of patch from max. seedrate per m2
  seedmax=spar(specnum)%seedrate*kpatchsize

   CALL  random_number(pequal)
   CALL  random_number(pequal)
   nseed=-seedmax*alog(1.-pequal)
   IF (flag_mg ==4 .and. time.eq.1) THEN
         nseed  = NINT(spar(specnum)%seedrate*kpatchsize)
   ELSE IF(flag_mg ==4.and. time.gt.1)THEN
          nseed = 0

   END IF
    end
    


SUBROUTINE seed_ini(nseed,nsp)

USE data_species
use data_stand
use data_help
use data_out
use data_simul
use data_soil

IMPLICIT NONE
integer    :: nseed, nr, j
integer    :: nsp
REAL       :: shoot
REAL       :: x1,x2,xacc,root
REAL       :: rtflsp
REAL	   :: troot2
TYPE(cohort)    ::tree_ini


external weight
external rtflsp

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' seed_ini'

IF(nseed.eq.0) RETURN


hnspec = nsp
max_coh = max_coh + 1

! nullify of all elements

 call coh_initial (tree_ini)

tree_ini%ident =  max_coh
tree_ini%species = nsp
tree_ini%ntreea = nseed
tree_ini%nta = nseed
 tree_ini%x_age = 1

mschelp = spar(nsp)%seedmass/1000.                              ! g ---> kg
x1 = 0.
x2 = 0.1
xacc = (1.0e-10) * (x1+x2)/2
root = rtflsp(weight,x1,x2,xacc)
tree_ini%x_sap = root
shoot = root*1000.                                           ! [kg]
tree_ini%x_fol= (spar(nsp)%seeda*(tree_ini%x_sap** spar(nsp)%seedb))   ![kg]            
tree_ini%x_frt = tree_ini%x_fol                                 ! [kg]
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

    ! initialize pheno state variables
    IF(spar(tree_ini%species)%Phmodel==1) THEN
       tree_ini%P=0
       tree_ini%I=1
    ELSE
       tree_ini%P=0
       tree_ini%I=0
       tree_ini%Tcrit=0
    END IF

! tranformation of shoot biomass kg --> mg

if(nsp.ne.2)tree_ini%height = spar(nsp)%pheight1*(shoot*1000.)**spar(nsp)%pheight2      ! [cm] calculated from shoot biomass (mg); berechnet aus shoot biomass (mg)
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(nseed.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

!     root distribution
        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 SUBROUTINE  seed_ini



SUBROUTINE growth_seed

USE data_stand
USE data_species
USE data_simul
USE data_par
use data_out

IMPLICIT NONE
  REAL   :: lambdaf = 0.,  &      ! partitioning functions
            lambdas = 0.,  &
            lambdar = 0.,  &
            NPP = 0.,      &      ! annual NPP
            F = 0.,        &      ! state variables: foliage,
            S = 0.,        &      ! shoot biomass,

            R = 0.,        &      ! fine roots,

            H = 0.,        &      ! total tree height

            FNew, SNew,    &      ! new states
            RNew,          &
            sigmaf = 0.,   &      ! current leaf activity rate
            sigman = 0.,   &      ! current root activity rate
            betar = 0.,    &
            ar = 0
 REAL     :: Sf,           &          ! senescence rates
            Sr,            &
            Gf,            &          ! growth rates
            Gs,            &
            Gr
 real     :: pab,helpdr,helpsum
 TYPE(coh_obj), POINTER :: p

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' growth_seed'

 flag_standup = 2   ! call stand_balance and root_distribution later 

  p=>pt%first
 DO
  if(.not.associated(p)) exit
  if( p%coh%height.lt.thr_height.and. p%coh%species.le.nspec_tree) then
  ns   = p%coh%species

  F    = p%coh%x_fol
  S    = p%coh%x_sap
  R    = p%coh%x_frt

  NPP = p%coh%NPP
  IF (flag_reg .eq. 2) NPP  = p%coh%NPPpool              ! [kg]
  H   = p%coh%height

  Sf = p%coh%sfol
  Sr = p%coh%sfrt

! only allocate if enough NPP is available

1  IF (NPP>1.0E-9.or. NPP.ge.(Sf+Sr).and.(sr+Sf)>1.0E-9) THEN

! calculate leaf activity based on net PS and leaf mass
  	sigmaf = NPP/F

! calculate root activity based on drought index
      helpdr= p%coh%drIndAl / p%coh%nDaysGr
      IF (flag_sign.eq.1) THEN
         sigman = amax1(spar(ns)%sigman*10*(((5.-spar(ns)%stol)*1.-p%coh%crown_area)/(5.-spar(ns)%stol)*1.),spar(ns)%sigman) * p%coh%drIndAl / p%coh%nDaysGr
      ELSE
         sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
      END IF
     ! auxiliary variables for fine roots
      if(helpdr.lt.0.001) ar = 0.

        ar = spar(ns)%pcnr * sigmaf / sigman
        betar = (Sr - R + ar*(F-Sf)) / NPP

! calculate coefficients for roots and foliage and shoot

      select case (ns)
         case (1)
           pab = 0.487
         case(2)
           pab = 0.826
         case(3)
           pab=1.9
         case(4)
          pab=1.002
! Pinus contorta
         case(6)
! Gholz
          pab = 0.236
! Populus tremula
          case(8)
             pab = 0.3233
      	  case(9)
! Pinus halepensis
             pab = 1.42335
         case(10)
! pseudotsuga menziesii
             pab = 0.8515
         case(11)
! Robinia
            pab = 0.8594
		end select

      lambdaf = (pab*(1-betar)+ (Sf/NPP))/(1 + pab*(1. + ar))
      lambdar = ar * lambdaf + betar
      lambdas = 1.- lambdaf - lambdar
      
! consistency
  ELSE

     lambdaf   = 0.
     lambdas   = 0.
     lambdar   = 0.

   END IF
   if ( lambdas.lt.0.) then
            lambdas = 0.
    		lambdaf = (1.-betar)/(ar+1)
    		lambdar = 1.-lambdaf
            if( lambdar.lt.0) then
			    lambdar=0
				lambdaf=1
			end if
			if(lambdaf.lt.0) then
			    lambdaf =1
				lambdar = 0.
			end if
	endif
helpsum = lambdaf + lambdar+ lambdas
   Gf = lambdaf*NPP
   Gr = lambdar*NPP
   Gs = lambdas*NPP
   p%coh%gfol = Gf
   p%coh%gfrt = Gr
   p%coh%gsap = Gs

  ! update of state vector 
  FNew  = F + Gf - Sf
  SNew  = S + Gs
  RNew  = R + Gr - Sr

  p%coh%x_fol  = FNew
  p%coh%x_sap  = SNew
  p%coh%x_frt  = RNew
  p%coh%fol_inc_old = p%coh%fol_inc
  p%coh%fol_inc = Gf - Sf
  p%coh%stem_inc = Gs
  
! update height and shoot base diameter (regression functions from Schall 1998) 

  IF(ns.ne.2) p%coh%height = spar(ns)%pheight1* (snew*1000000.) **spar(ns)%pheight2
  IF(ns.eq.2) p%coh%height = 10**(spar(ns)%pheight1+ spar(ns)%pheight2*LOG10(snew*1000000.)+ &
                              spar(ns)%pheight3*(LOG10(snew*1000000.))**2)
   p%coh%height_ini =   p%coh%height
   
! update foliage area, parameter med_sla 

      SELECT CASE (flag_light)
        CASE (1:2)
   p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.- vstruct(lowest_layer)%irel)
        CASE(3,4)
    p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.-irelpool(lowest_layer))
      END SELECT

! total leaf area of a tree in this cohort [m**2]
    p%coh%ca_ini = p%coh%med_sla * p%coh%x_fol

   ! update age -now not necessary this is done in stand_bal 
   p%coh%notViable = (FNew <= 0.) .OR. (SNew <= 0.) .OR.    &
                    (RNew <= 0.)
  p%coh%litC_fol  = p%coh%litC_fol +  p%coh%ntreea * Sf * cpart
  p%coh%litC_frt  = p%coh%litC_frt +  p%coh%ntreea * Sr * cpart

  ! with species specific N content and reallocation factor (see species.par)
  p%coh%litN_fol  = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * spar(ns)%reallo_fol / spar(ns)%cnr_fol
  p%coh%litN_frt  = p%coh%litN_frt +  p%coh%ntreea * Sr * cpart * spar(ns)%reallo_frt / spar(ns)%cnr_frt

  end if ! seedling cohort test

  p=> p%next

  END DO

END SUBROUTINE growth_seed


SUBROUTINE mort_seed

USE data_species
USE data_simul
use data_stand
use data_par
use data_out

IMPLICIT NONE

integer     ::  taxnr
integer     ::  hage
real     :: intmort
real     :: strmort
real     :: totmort
real     :: ntdead
real     :: ntahelp
TYPE(coh_obj), POINTER :: p

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' mort_seed'

p=>pt%first
DO
   IF(.not.associated(p)) EXIT
   IF(p%coh%height.lt.thr_height) THEN
       IF(p%coh%notViable) then
          PRINT*,time, p%coh%notViable
          p%coh%ntreed = p%coh%ntreea
          p%coh%ntreea = 0

        ENDIF
    END IF
   p => p%next

END DO

p => pt%first

DO

   IF(.not.associated(p)) EXIT
   IF(p%coh%height.lt.thr_height .and. p%coh%species.le.nspec_tree) THEN
   taxnr = p%coh%species
       if(p%coh%ntreea .eq.0) goto 1000

       hage =  p%coh%x_age
       IF( p%coh%fol_inc.le.0.) THEN
           p%coh%x_stress = p%coh%x_stress +1
           p%coh%x_health = 0

       ELSE
           p%coh%x_health = p%coh%x_health +1
           IF(p%coh%x_stress .gt. 0.) p%coh%x_stress = p%coh%x_stress -1

       ENDIF
! intrinsic mortality 
       CALL int_mort_weib(taxnr, intmort, hage)

! stress mortality 
       IF(p%coh%x_stress.gt.0) THEN
           strmort = weibal*spar(taxnr)%weibla*p%coh%x_stress**(weibal-1)
       ELSE
          strmort = 0.
       ENDIF

        totmort=intmort+(1-intmort)*strmort

! calculation of real number of dying stems
       ntdead = totmort*p%coh%ntreeA
! update of real stem number nta and number of dead stems 
       p%coh%nta =  p%coh%nta -ntdead
       p%coh%nTreeD = p%coh%nTreeA-NINT(p%coh%nta)
! help variable for comparison
        ntahelp = p%coh%nTreeA
! update of integer stem number 
        p%coh%nTreeA = NINT(p%coh%nta)
! update of  integer stem number 
       if(p%coh%nta.lt.1.) p%coh%nTreeA=0.
! update of real stem number if integer stem number was changed 
       if (ntahelp .ne.  p%coh%nTreeA ) p%coh%nta = p%coh%nTreeA
1000   if (p%coh%ntreeD.ne.0) then
         p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreeD*(1.-spar(taxnr)%psf)*p%coh%x_fol*cpart
         p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreeD*((1.-spar(taxnr)%psf)*p%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
         p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreeD*p%coh%x_frt*cpart
         p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreeD*p%coh%x_frt*cpart/spar(taxnr)%cnr_frt
         p%coh%litC_tb = p%coh%litC_tb + p%coh%ntreeD*p%coh%x_tb*cpart
         p%coh%litN_tb = p%coh%litN_tb + p%coh%ntreeD*p%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
         p%coh%litC_crt = p%coh%litC_crt + p%coh%ntreeD*p%coh%x_crt*cpart
         p%coh%litN_crt = p%coh%litN_crt + p%coh%ntreeD*p%coh%x_crt*cpart/spar(taxnr)%cnr_tbc

         p%coh%litC_stem = p%coh%litC_stem + p%coh%ntreeD*(p%coh%x_sap)*cpart
         p%coh%litN_stem = p%coh%litC_stem/spar(taxnr)%cnr_stem
       endif

   END IF
   p => p%next

ENDDO

END SUBROUTINE mort_seed