CalcCurrentDepFunctsInCell.f90

Go to the documentation of this file.

!==============================================================================
!
!  Name   :  CalcCurrentDepFunctsInCell
!
!  Purpose:  CalcCurrentDepFunctsInCell calculates the spatial forest dynamics
!             on a grid with or without seed dispersal. ContLightDep calculates the
!            continuous light dependence of growth and death.
!
!  Remarks:  contains the SUBROUTINEs
!
!            CalcCurrentDepFunctsInCell (<LocalDynOneTimeStep<SpatialDynOneTimeStep<
!                                          TimeLoop<TreeMig)
!            ContLightDep (<GetSpc<TreeMig AND/OR <Regen<LocalDynOneTimeStep<
!                            LocalDynOneTimeStep<SpatialDynOneTimeStep<TimeLoop<TreeMig)
!
!==============================================================================
! design          : H. Lischke, N. Zimmermann
! author(s)       : H. Lischke, N. Zimmermann
! implementation  : H. Lischke, N. Zimmermann
! cleaner         : T.J. Loeffler
! copyright       : (c) 1999, 03 by H. Lischke
!==============================================================================
!===================================================================== 
!===============================================================
SUBROUTINE calccurrentdepfunctsincell(nspc, ddegs, wtemp, drstr, avnit, brwpr, germDrought)
 
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    use all_par, only: maxspc, maxhc, maxlc, birht, includeenv, kn1, kn2, &
                       speccl, spec, estabeqhc0
 
    IMPLICIT NONE
 
    ! <Passed variables>-------------------------------------------------------------------------------
    INTEGER, INTENT(in) :: nspc           ! Passed number of species                             [=spec]
    REAL, INTENT(in) :: ddegs, & ! day-degrees of cell                                         [=dayDegs] !ENV
                        drstr, & ! drough stress of cell                                [=rDroughtStress] !ENV
                        avnit, & ! av. Nitrogen of cell                                     [=avNitrogen] !ENV
                        wtemp, & ! winter temperature of cell                               [=winterTemp] !ENV
                        brwpr, &           ! browsing pressure of cell                                [=browsPress] !ENV
                        germDrought ! specific drought values affecting germination
 
    ! <Local variables>---------------------------------------------------------------------------------
    INTEGER :: i               ! looping variable for species
    REAL    :: truht, & ! tmp var to calculate heights                             [=trueHeight]
               diams, & ! tmp var to calculate spec. diameters       [in DetermineProc..  =diam]
               alloc1, & ! tmp allocation var to calculate growth rate                   [=allo1]
               alloc2, & ! tmp allocation var to calculate growth rate                   [=allo2]
               grofac, & ! diameter increment growth factor                         [=diamIncFak]
               hcgrow, & ! height class dep. growth rate                             [=growthInH]
               lidep, & ! light dependent growth                                        [=lidep]
               tdep(maxspc), & ! temperature dependent growth factor                   [=growthTempDep] !ENV
               drdep(maxspc), & ! drought dependent growth factor                    [=growthDroughtDep] !ENV
               ndep(maxspc), & ! nitrogen dependent growth factor                         [=growthNDep] !ENV
               germDrDep, & ! germ suitability dependent germination factor                         [=growthNDep] !ENV
               envdpg, & ! local env. dependent growth factor                      [=allOtherDep] !ENV
               envetb, & ! temp var. for env. sensitivity to establishment         [=..otheretab] !ENV
               brwetb           ! temp var. for browsing sensit. to establ.               [=..otheretab] !ENV
    INTEGER::           htcl, ltcl
    ! <Here the SUBROUTINE starts>*********************************************************************
    !----- Check whether site data drive dynamics ----------------------------------------------------
    include_env: if (includeenv) then
        !----- Environment dependence of growth and survival                                                                                        !
        do i = 1, nspc
            !----- day degree dependence   (saturation function)                                                                                         !
            tdep(i) = max(0.0, 1.0 - exp((1.386/spec(i)%kddmax)*(spec(i)%kddmin - ddegs)))
            !----- drought dependence
            drdep(i) = (max(1.0 - drstr/spec(i)%kdrt, 0.0)**0.5)
                 !----- make drought dependence per cell an array later(!), very important!!
            !----- nitrogen dependence
            ndep(i) = max(1.0 - exp(kn1(spec(i)%kntol)*(avnit - kn2(spec(i)%kntol))), 0.0)
            spec(i)%envdep = tdep(i)*drdep(i)*ndep(i)
        end do
        !----- establishment dependence on day degree sum and winter temperature
        do_etab: do i = 1, nspc
            env_etab: if ((wtemp < real(spec(i)%kwit)) .OR. (ddegs < real(spec(i)%kddmin))) then
                envetb = 0.0
            else
                envetb = 1.0
            end if env_etab
            !----- establishment dependence on browsing
            brwetb = 1.0 - (real(spec(i)%kbrow - 1)*brwpr/30.0)
             !----- establishment dependence on drought 
            germdrdep  = max( 1.0  -   germdrought /  spec(i)%kGermDrSens, 0.0) ** 0.5 
           
            spec(i)%enveta = envetb * brwetb * germdrdep 
!             print *, "germDrought, drstr, spec, enveta,germDrDep , germDrought,kGermDrSens, drdep, kdrt",&
!                     germDrought, drstr, i, spec(i)%enveta,  &
!                     germDrDep,spec(i)%kGermDrSens, drdep(i), spec(i)%kdrt 
        end do do_etab
 
    else  !----- Include_ENV
        spec%envdep = 1.0
        spec%enveta = 1.0
    end if include_env
    !-----<Initialize env. dependent tree arrays [=DetermineProcRates]; this has moved from InitAndSetConstHeightClProps>-------
 ! >    environment dependent vitality, height growth and diameter growth for species, height class and light class
   byspcs2: do i = 1, nspc
        byhtcls2: do htcl = spec(i)%maxhtc, 1, -1
            envdpg = spec(i)%envdep**0.333333
            diams  = speccl(i)%diain(htcl)       ! diameter in height class
            truht  = speccl(i)%htmin(htcl)       ! lower bound of height class
            alloc1 = (2.0*(real(spec(i)%khmax) - birht)/spec(i)%maxdhm)
            alloc2 = -alloc1/2.0/spec(i)%maxdhm
            if (htcl <= spec(i)%maxhtc) then ! diameter growth for optimal climate, depending on current diameter
                grofac = (envdpg*spec(i)%maxgro*diams*(1.0 - (truht/real(spec(i)%khmax))))/ &
                         (2.74 + (((3.0*alloc1) + (4.0*alloc2*diams))*diams))
            else
                grofac = 0.0
            end if
 
            hcgrow = (grofac*(alloc1 + (2.0*alloc2*diams))) ! translate diameter growth to height growth
            !----- light dependence values in light classes of growth and mortality
            byltcls2: do ltcl = 1, maxlc  !get the environment dependent vitality, height growth and diameter growth for this species, height class and light class  
  !              lih = ltval(ltcl)
 
                !----- THIS STEP PRODUCES MOST OF THE COMPUTING TIME
                !----->OLD        call ContLightDep (lih, spec(i)%kligha, lidep)
                lidep = speccl(i)%lightDepA(ltcl)  ! light dependence of this species in this light class;! lightDepA has been calculated in GetSpc  
                speccl(i)%envVitality(htcl, ltcl)    = lidep*envdpg ! vitality        in this height and light class for this species   
                speccl(i)%hlgro(htcl, ltcl)          = lidep*hcgrow ! height growth   in this height and light class for this species  
                speccl(i)%growthVitality(htcl, ltcl) = lidep*grofac ! diameter growth in this height and light class for this species  
            end do byltcls2
        end do byhtcls2
 
        !----- function values of 0th height class (seedlings to 1.37 m)
 
        byltcls4: do ltcl = 1, maxlc
            if (estabeqhc0) then
                !------ if wished, light dependence of establishment (from FORCLIM) is contributed to 0th height class
!                lih = ltval(ltcl)
 
                !----- THIS STEP PRODUCES MOST OF THE COMPUTING TIME
                !----->OLD        call ContLightDep (lih, spec(i)%kligha, lidep)
                lidep = speccl(i)%lightDepE(ltcl) ! lightDepE has been calculated in GetSpc
                speccl(i)%envVitality(0, ltcl) = lidep*envdpg
 
                speccl(i)%hlgro(0, ltcl) = lidep*hcgrow
                speccl(i)%growthVitality(0, ltcl) = lidep*grofac
            else
                !----- otherwise the values for the first height class are also used for the 0th one (seedlings)
                speccl(i)%hlgro(0, ltcl)          = speccl(i)%hlgro(1, ltcl)
                speccl(i)%envVitality(0, ltcl)    = speccl(i)%envVitality(1, ltcl)
                speccl(i)%growthVitality(0, ltcl) = speccl(i)%growthVitality(1, ltcl)
            end if
        end do byltcls4
! set to zero growth, vitality and growthVitality in height classes higher than highest one 
        if (spec(i)%maxhtc < maxhc) then
            resthtcls2: do htcl = spec(i)%maxhtc + 1, maxhc
                byltcls3: do ltcl = 1, maxlc
                    speccl(i)%hlgro(htcl, ltcl) = 0.0
                    speccl(i)%envVitality(htcl, ltcl) = 0.0
                    speccl(i)%growthVitality(htcl, ltcl) = 0.0
                end do byltcls3
            end do resthtcls2
        end if
    end do byspcs2
end SUBROUTINE calccurrentdepfunctsincell
!****<Here the SUBROUTINE ends>********************************************************************
 
!****<Here the SUBROUTINE ContLightDep starts>*****************************************************
!===================================================================== 
!===============================================================
SUBROUTINE contlightdep(light, shade_sens_class, lidep)
    use loggermodule, only: logwarning, logmessage
    ! <CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    IMPLICIT NONE
 
    ! <Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(in) :: shade_sens_class    ! Spec. specific light threshold
    REAL, INTENT(in)    :: light        ! Light intensity
    REAL, INTENT(out)   :: lidep        ! Light dependence
 
    ! <Local variables>---------------------------------------------------------------------------------
    REAL :: gl1, gl1old, & ! light response for shade-tolerant species                  [=gl1]
            gl9, gl9old, & ! light response for shade-intolerant species                [=gl9]
            eps, steepness, &
            contstep1, contstep9
 
! <Here the SUBROUTINE starts>*********************************************************************
    eps = 0.00000001
    steepness = 5
    contstep1 = eps*(1 - exp(-((light/0.05)**steepness)))
    contstep9 = eps*(1 - exp(-((light/0.08)**steepness)))
    gl1old = 1.0 - exp(-4.64*(light - 0.05))
    gl1 = max(gl1old + eps, 0.0) + contstep1
    gl9old = 2.24*(1.0 - exp(-1.136*(light - 0.08)))
    gl9 = max(gl9old + eps, 0.0) + contstep9
    lidep = (gl1 + ((shade_sens_class - 1.0)*(gl9 - gl1)/8.0))
    if ((lidep < 0.0) .OR. (lidep > 2)) then
        write (logmessage, *) 'Negative LightDep:'
        call logwarning(logmessage)
        write (logmessage, "(9A8)") "light", "contstep1", "contstep9", "gl1old", "gl9old", "gl1", "gl9", "shade_sens_class", "lidep"
        call logwarning(logmessage)
        write (logmessage, "(7F8.3,I8,F8.3)") light, contstep1, contstep9, gl1old, gl9old, gl1, gl9, shade_sens_class, lidep
        call logwarning(logmessage)
    end if
 
    lidep = max(lidep, 0.0)**0.333333
end SUBROUTINE contlightdep
!done