!*****************************************************************!
!* *!
!* 4C Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* Simulation of processes at subannual resolution *!
!* *!
!* Contains subroutines: *!
!* *!
!* - pheno_ini *!
!* - pheno_begin *!
!* - pheno_count *!
!* - pheno_shed *!
!* *!
!* functions: *!
!* triangle *!
!* *!
!* 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 pheno_ini
USE data_climate
USE data_simul
USE data_site
USE data_species
USE data_stand
IMPLICIT NONE
integer i, j
integer leapyear
real atemp, hh, htemp
real triangle
real, external :: daylength
leaves_on = .false.
all_leaves_on = 0
phen_flag=1 ! CANOPY is calculated once at the beginning of each year
! Initialising of all species is done at the beginning, since if species information wouldnt be initialised
IF(time==1) THEN
do i=1,nspec_tree
ns = i
IF(spar(ns)%Phmodel==1) THEN
svar(ns)%Pro = 0.
svar(ns)%Inh = 1.
ELSE
svar(ns)%Pro = 0.
svar(ns)%Inh = 0.
svar(ns)%Tcrit = 0.
END IF
! initialize pheno state variables with climate from the actual year
do j = spar(ns)%end_bb+1, 365
atemp = tp(j, 1)
hh = DAYLENGTH(j,lat)
SELECT CASE(ns)
CASE(1,8)
!Fagus
! Promotor-Inhibitor model 11
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)* &
(1-svar(ns)%Inh)*hh/24 - &
spar(ns)%PPb*svar(ns)%Pro*(24-hh)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa* &
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)* &
svar(ns)%Inh*hh/24
CASE(4)
! Quercus
! Promotor-Inhibitor model 12
htemp = triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa * htemp * &
(1-svar(ns)%Inh) * hh/24
htemp = triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa * htemp * &
svar(ns)%Inh * hh/24 + spar(ns)%PPb*(24-hh)/24
CASE(5, 11)
! Betula, Robinia
IF(spar(ns)%Phmodel==1) THEN
! Promotor-Inhibitor model 2
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,atemp)* &
(1-svar(ns)%Inh) - spar(ns)%PPb*svar(ns)%Pro*(24-hh)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa* &
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,atemp)*svar(ns)%Inh
END IF
END SELECT
enddo ! j
Enddo ! nspec_tree
END IF
! latest day of bud burst 30. of June (DOY 181+leapyear(time_cur))
do i=1, anrspec
ns = nrspec(i)
if(ns.le.nspec_tree) then
IF(spar(ns)%phmodel==4) THEN
svar(ns)%daybb = svar(ns)%ext_daybb
ELSE
svar(ns)%daybb = 181 + leapyear(time_cur)
ENDIF
end if
END DO ! anrspec
end SUBROUTINE pheno_ini
!*******************************************************************
SUBROUTINE pheno_begin
! calculation of day_bb, latest day of bud burst 30. june (DOY 181)
USE data_simul
USE data_species
USE data_stand
USE data_climate
USE data_site
IMPLICIT NONE
REAL triangle
INTEGER leapyear
real hh, htemp
integer i
hh = dlength
do i=1, anrspec
ns = nrspec(i)
if (iday .ge.364) then
continue
endif
if(ns.le.nspec_tree .OR. ns.eq.nspec_tree+2) then !either tree or mistletoe
! Pheno model
select Case (spar(ns)%Phmodel)
Case(0) ! no model
!Picea, Pinus, Mistletoe
IF(iday.EQ.1) THEN
svar(ns)%daybb = iday
phen_flag = 1
leaves_on = .TRUE.
ENDIF
Case(1)
! Phenology starts after leaf coloring/shedding and ends not later than 30. June
IF (iday > spar(ns)%end_bb+1 .OR. iday <= svar(ns)%daybb) THEN
SELECT CASE(ns)
CASE(1,8)
!Fagus
! Promotor-Inhibitor model 11
htemp = triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,airtemp)
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa * htemp * &
(1-svar(ns)%Inh) * dlength/24 - &
spar(ns)%PPb*svar(ns)%Pro * (24-dlength)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa*&
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,airtemp)*&
svar(ns)%Inh*dlength/24
IF (svar(ns)%Pro >= 1) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < 1 .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
CASE(4)
! Quercus
! Promotor-Inhibitor model 12
all_leaves_on=0
if (svar(ns)%Inh .gt. 1.) then
continue
svar(ns)%Inh = 1.
endif
if (svar(ns)%Pro .lt. 0.) then
continue
svar(ns)%Pro = 0.
endif
htemp = triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,airtemp)
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa * htemp * &
(1-svar(ns)%Inh) * dlength/24
htemp = triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,airtemp)
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa * htemp * &
svar(ns)%Inh * dlength/24 + spar(ns)%PPb*(24-dlength)/24
IF (svar(ns)%Pro >= 1) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < 1 .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
CASE(5, 11)
! Betula, Robinia
all_leaves_on=0
IF(spar(ns)%Phmodel==1) THEN
! Promotor-Inhibitor model 2
svar(ns)%Pro = svar(ns)%Pro + spar(ns)%PPa* &
triangle(spar(ns)%PPtmin,spar(ns)%PPtopt,spar(ns)%PPtmax,airtemp)* &
(1-svar(ns)%Inh) - spar(ns)%PPb*svar(ns)%Pro*(24-dlength)/24
svar(ns)%Inh = svar(ns)%Inh - spar(ns)%PIa* &
triangle(spar(ns)%PItmin,spar(ns)%PItopt,spar(ns)%PItmax,airtemp)*svar(ns)%Inh
IF (svar(ns)%Pro >= 1) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < 1 .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
END IF
END SELECT
Endif
Case(2)
! Cannel-Smith model
IF(iday >= 305 + leapyear(time_cur) .OR. iday <= svar(ns)%daybb) THEN
IF(airtemp < spar(ns)%CSTbC) THEN
svar(ns)%Inh = svar(ns)%Inh + 1
svar(ns)%Tcrit = spar(ns)%CSa + spar(ns)%CSb*LOG(svar(ns)%Inh)
END IF
IF(airtemp > spar(ns)%CSTbT .AND. iday >= 32 .AND. iday <= svar(ns)%daybb) THEN
svar(ns)%Pro = svar(ns)%Pro + airtemp - spar(ns)%CSTbT;
END IF
IF(svar(ns)%Pro > svar(ns)%Tcrit) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < svar(ns)%Tcrit .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
END IF
Case(3)
! Temperature sum model
SELECT CASE(ns)
CASE(11)
! Robinia
IF(iday >= spar(ns)%Lstart .AND. iday <= svar(ns)%daybb) THEN
IF(airtemp > spar(ns)%LTbT) THEN
svar(ns)%Pro = svar(ns)%Pro + airtemp
END IF
IF(svar(ns)%Pro > spar(ns)%LTcrit) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < spar(ns)%LTcrit .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
END IF
CASE default
IF(iday >= spar(ns)%Lstart .AND. iday <= svar(ns)%daybb) THEN
IF(airtemp > spar(ns)%LTbT) THEN
svar(ns)%Pro = svar(ns)%Pro + airtemp - spar(ns)%LTbT
END IF
IF(svar(ns)%Pro > spar(ns)%LTcrit) THEN
svar(ns)%daybb=iday
phen_flag = 1
leaves_on=.TRUE.
ELSE IF (svar(ns)%Pro < spar(ns)%LTcrit .AND. iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
END IF
END SELECT
Case(4)
! externally prescribed day of budburst
IF(iday==svar(ns)%daybb) THEN
phen_flag = 1
leaves_on=.TRUE.
END IF
Case default
IF(iday.EQ.1) THEN
svar(ns)%daybb=iday
phen_flag=1
leaves_on=.TRUE.
ENDIF
end select
else if(iday==svar(ns)%daybb) then
phen_flag = 1
leaves_on=.TRUE.
end if
END DO
zeig=>pt%first
do while (associated(zeig))
ns = zeig%coh%species
zeig%coh%day_bb = svar(ns)%daybb
zeig=>zeig%next
enddo
END SUBROUTINE pheno_begin
!*******************************************************************
SUBROUTINE pheno_count
USE data_simul
USE data_species
USE data_stand
IMPLICIT NONE
zeig=>pt%first
DO
if(.not. associated(zeig)) exit
! vegetation period per PS-time step and per season
IF((iday >= zeig%coh%day_bb) .AND. (iday <= spar(zeig%coh%species)%end_bb)) THEN
zeig%coh%nDaysPS = zeig%coh%nDaysPS + 1. ! set to 0 in npp
zeig%coh%nDaysGr = zeig%coh%nDaysGr + 1. ! set to 0 year_ini
END IF
zeig=>zeig%next
END DO
END SUBROUTINE pheno_count
!*******************************************************************
SUBROUTINE pheno_shed
USE data_simul
USE data_species
USE data_stand
IMPLICIT NONE
integer i
leaves_on=.FALSE.
all_leaves_on=1
DO i=1, anrspec
ns = nrspec(i)
IF(iday == spar(ns)%end_bb +1) THEN
phen_flag=1
all_leaves_on=0
! reset pheno state variable
IF(spar(ns)%Phmodel==1) THEN
svar(ns)%Pro = 0.
svar(ns)%Inh = 1.
ELSE
svar(ns)%Pro = 0.
svar(ns)%Inh = 0.
svar(ns)%Tcrit = 0.
END IF
ELSE IF((iday < svar(ns)%daybb) .OR. (iday > spar(ns)%end_bb)) THEN
all_leaves_on=0
ELSE IF((iday >= svar(ns)%daybb) .AND. (iday <= spar(ns)%end_bb)) THEN
leaves_on=.TRUE.
END IF
END DO
END SUBROUTINE pheno_shed
!*******************************************************************
FUNCTION triangle(min,opt,max,x)
REAL :: min,opt,max,x,triangle
IF( min <= x .AND. x <= opt) THEN
triangle = (x - min)/(opt - min)
ELSE IF( opt < x .AND. x <= max) THEN
triangle = (max - x)/(max - opt)
ELSE
triangle = 0
END IF
END FUNCTION triangle
FUNCTION leapyear(year)
INTEGER :: year,leapyear
IF( MOD(year,400)==0 .OR. ( MOD(year,100)/=0 .AND. MOD(year,4)==0 )) THEN
leapyear = 1
ELSE
leapyear = 0
END IF
END FUNCTION leapyear