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

