!*****************************************************************!
!*                                                               *!
!*                    4C  Simulation Model                       *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*       Simulation of processes at subannual resolution         *!
!*                                                               *!
!*                                                               *!
!* Contains subroutines:                                         *!
!*                                                               *!
!* STAND_DAILY                                                   *!
!* SET_PS                                                        *!
!* DROUGHT     : Calculation of drought stress indices           *!
!* FIRE_RISK                                                     *!
!* calc_frost_index  : calculation of indices for frost damage   *!
!* calc_endbb   : calculation of end of the vegetation period    *!
!*                                                               *!
!*                  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_daily

  !*** Declaration part ***!

  USE data_stand
  USE data_simul
  USE data_species
  USE data_climate
  USE data_site
  USE data_soil_cn
  USE data_out
  USE data_par
  USE data_evapo
  USE data_soil
  use data_manag

  IMPLICIT NONE

  REAL    :: aveT,   &  ! average of temperature for PS/NPP models
             avDL,   &  ! average of daylength for PS/NPP model
             avRD,   &  ! average of radiation
             avPR,   &  ! average of pressure (hPa)
             PAR        ! average of PAR for PS/NPP model [mol quanta d-1]
  
  REAL    :: hdfr, hdt, hprs
  INTEGER :: i, jd, k, d, week, monthday, ns_pro_help
  

  real                  :: p_help, t_help

  REAL     :: photoper

  p_help=0.
  t_help=0.

  irelpool_ll=0.
  bgpool_ll=0.

  !*** Calculation part ***!

  week     = 0
  monthday = 0 
  monat    = 1
  woche    = 1

  ! daily loop
  DO jd = 1, recs(time)

    iday = jd
    monthday=monthday+1
   
    ! input of daily climate data
    CALL day_ini
    
    if(anz_coh .gt. 0) then     ! if no cohort, then no phaenology necessary
        IF(all_leaves_on==0) CALL pheno_begin
        CALL pheno_count
        IF(leaves_on) CALL pheno_shed
    endif
    IF(phen_flag==1 .OR. (.not.flag_tree .and. leaves_on)) THEN

       ! Calculate this year's crown geometry for each cohort, followed by
       ! leaf area and light profiles across the canopy
       CALL CANOPY
       if (anz_coh.eq.0) then
           irelpool_ll = 1.
       end if
       if(all_leaves_on.eq.1) then
           irelpool_ll = irelpool(0)
           bgpool_ll = bgpool(2)
       end if
       IF(flag_end.EQ.3) RETURN

       ! update of stand variables (LAI, cover)
       CALL standup
       phen_flag=0;
    END IF

   !call distubance after start day
    select case(flag_dis)
     case(1,2)
      if (dis_control(1,1) .eq. 1) then
        if(all_leaves_on .eq. 1 .and. dis_start(dis_control(1,2)) .eq. iday) then
         CALL disturbance_defoliator
         CALL CANOPY
         CALL stand_balance
         CALL standup
        endif
      endif
      if (dis_control(2,1) .eq. 1) then
        if(all_leaves_on .eq. 1 .and. dis_start(dis_control(2,2)) .eq. iday) CALL disturbance_xylem
      endif
      if (dis_control(3,1) .eq. 1) then
        if(dis_start(dis_control(3,2)) .eq. iday) CALL disturbance_phloem
      endif
      if (dis_control(4,1) .eq. 1) then
        if(dis_start(dis_control(4,2)) .eq. iday) then
         CALL disturbance_root
         CALL stand_balance
         CALL standup
        endif
      endif
      if (dis_control(5,1) .eq. 1) then
        if(dis_start(dis_control(5,2)) .eq. iday) CALL disturbance_stem
      endif
    end select
    ns_pro_help = ns_pro
    ! set ns_pro_help to length of last photosynthesis period at end of year
    IF(iday >int(recs(time)/ns_pro)*ns_pro .and. (MOD( iday, ns_pro )==1)) THEN
       ns_pro_help = recs(time) - int(recs(time)/ns_pro)*ns_pro
    END IF
    
    ! optimum photosynthesis submodel
    IF (ns_pro==1.OR.(MOD( iday, ns_pro )==1) .or. iday.eq.1) THEN
      ! assign averaged input variables for PS model
      aveT = 0.
      avDL = 0.
      avRD = 0.
      avPR = 0.
      hdfr = 0.
      ns_day = 1
      DO k = 1, ns_pro_help     ! this calculates 365 or 366, but is not included as a wwek value
                                ! ==> last week of the year is recieving this amount
        d = iday-1+k
        hdt  = Q10_T**((tp(d,time) - 15.) / 10.)
        hdfr = hdfr + hdt
        dayfract(k) = hdt
        aveT = aveT + tp(d,time) + deltaT
        avRD = avRD + rd(d,time)
            hprs = prs(d,time)
            if (hprs .lt. 800.) then
                  hprs = 1013
            endif
        avPR = avPR + hprs
        avDL = avDL + photoper( FLOAT(d), xLat )
      END DO
      aveT = aveT / ns_pro_help
      avDL = avDL / ns_pro_help
      avRD = avRD / ns_pro_help
      avPR = avPR / ns_pro_help
      ! PAR that is coming in stand reflection is substracted
      PAR    = (1.-pfref)* GR_in_PAR * avRD
      if (iday .gt. 364) then
        dayfract = 1.   ! at the last days of the year no temperature depending daily fraction of flux
      else
        dayfract = ns_pro * dayfract / hdfr   ! temperature depending daily fraction of flux, calc. from sum of ns_pro days
      endif
      CALL OPT_PS( aveT, avDL, PAR, avPR )
    ENDIF

    ! aggregation of stomatal conductance of the canopy
    gp_can_mean = gp_can_mean + gp_can
    gp_can_min  = min(gp_can_min, gp_can)
    gp_can_max  = max(gp_can_max, gp_can)

    ! soil submodel
    CALL SOIL
    CALL  drought

    ! NPP submodel
    IF (ns_pro==1.OR.(MOD( (iday-1), ns_pro )==0) .or. iday .eq. recs(time) .or. iday.eq.1) THEN
       CALL NPP( aveT, avDL, PAR, ns_pro_help )
       IF(.not.flag_tree .and. leaves_on.and.flag_sprout.eq.1)  CALL growth_seed_week (ns_pro_help)
       ! daily output every ns_pro days of dips- and gsdps-files
       IF (flag_dayout .ge. 1) CALL coh_out_d(2)
    ENDIF
    
    CALL calc_fire_risk
! calculation of the start of vegetation period 
    if(flag_vegper.eq.0) then
      if(airtemp.le.5. .and. flag_tveg .ne.0)  then
        flag_tveg=0
      else if(airtemp.gt.5. .and. flag_tveg.eq.0) then
        flag_tveg =1
      else if(airtemp.gt.5. .and. flag_tveg.eq.1) then
        flag_tveg =2
      else if(airtemp.gt.5. .and. flag_tveg.eq.2) then
        flag_tveg =3
      else  if(airtemp.gt.5. .and. flag_tveg.eq.3)then
         flag_tveg =4
      else if(airtemp.gt.5. .and. flag_tveg.eq.4) then
         flag_tveg =5
     end if

     if(flag_tveg .eq.5) then
       flag_vegper=1
       iday_vegper = iday
     end if
   endif
    
 ! call of SR for calculation of various indices for the frost index
   if(airtemp_min .gt. -90.) call calc_frost_index
 ! Calculation of maximal radiation (for information only)
    call glob_rad(dlength, iday, lat, rad_max)
    Cout%NEE(iday)      = respsoil - dailyNPP_C    ! g C/m²
    Cout%Resp_aut(iday) = dailyautresp_C * dayfract(ns_day)
    NPP_day          = dailyNPP_C * dayfract(ns_day)
    GPP_day          = (dailyNPP_C + dailyautresp_C) * dayfract(ns_day)
    TER_day          = dailyautresp_C * dayfract(ns_day) + respsoil
    
    IF (flag_dayout .ge. 1) CALL outday(1)
    IF (ns_pro==1.OR.(MOD( iday, ns_pro )==0) .or. iday .eq. recs(time) )  CALL SET_PS
  ! Wochen- und Monatswerte berechnen
    aet_mon(monat)  = aet_mon(monat) + aet
    aet_week(woche) = aet_week(woche) + aet
    pet_mon(monat)  = pet_mon(monat) + pet
    pet_week(woche) = pet_week(woche) + pet
    temp_mon(monat)  = temp_mon(monat) + airtemp
    temp_week(woche) = temp_week(woche) + airtemp
    prec_mon(monat)  = prec_mon(monat) + prec
    prec_week(woche) = prec_week(woche) + prec
    rad_mon(monat)   = rad_mon(monat) + rad
    hum_mon(monat)   = hum_mon(monat) + hum
    perc_mon(monat)  = perc_mon(monat) + perc(nlay)
    perc_week(woche) = perc_week(woche) + perc(nlay)
    resps_mon(monat) = resps_mon(monat) + respsoil
    resps_week(woche)= resps_week(woche) + respsoil
    GPP_mon(monat)   = GPP_mon(monat) + dailyNPP_C + dailyautresp_C
    GPP_week(woche)  = GPP_week(woche) + dailyNPP_C + dailyautresp_C
    NEE_mon(monat)   = NEE_mon(monat) + Cout%NEE(iday)                      ! g C/m²
    NPP_mon(monat)   = NPP_mon(monat) + dailyNPP_C
    NPP_week(woche)  = NPP_week(woche) + dailyNPP_C
    TER_mon(monat)   = TER_mon(monat) + dailyautresp_C + respsoil
    TER_week(woche)  = TER_week(woche) + dailyautresp_C + respsoil
    tempmean_mo(monat)   = tempmean_mo(monat) + airtemp   ! long-term monthly means
    
   ! summation output with variabel time steps
    photsum   = photsum + phot_C
    npppotsum = npppotsum + dailypotNPP_C
    nppsum    = nppsum + dailyNPP_C
    resosum   = resosum + respsoil
    nee       = nee + respsoil - dailyNPP_C
	gppsum    = gppsum + GPP_day
	sumGPP    = sumGPP + dailyNPP_C + dailyautresp_C
	sumTER    = sumTER + dailyautresp_C + respsoil 
	resautsum = resautsum + dailyautresp_C
    precsum   = precsum + prec
    tempmean  = tempmean + airtemp
	tempmeanh = tempmeanh +airtemp
	aet_sum   = aet_sum + aet
	pet_sum   = pet_sum + pet
    perc_sum  = perc_sum + perc(nlay)

    if(monthday==monrec(monat)) then
       tempmeanh = tempmeanh/monrec(monat)
	   if(monat.eq.1) med_air_cm = tempmeanh
	   if(tempmeanh.lt.med_air_cm) med_air_cm = tempmeanh
	   if(tempmeanh.gt.med_air_wm) med_air_wm = tempmeanh
	   tempmeanh = 0.

	   temp_mon(monat) = temp_mon(monat) / monrec(monat)
	   rad_mon(monat)  = rad_mon(monat) / monrec(monat)
	   hum_mon(monat)  = hum_mon(monat) / monrec(monat)
	   if(temp_mon(monat).lt.med_air_cm) med_air_cm = temp_mon(monat)
	   if(temp_mon(monat).gt.med_air_wm) med_air_wm = temp_mon(monat)       
	end if

	   if(airtemp.ge.10.) then
	       t_help= t_help + airtemp
	       p_help= p_help + prec
	   end if

    ns_day = ns_day + 1
 
  ! daily output
    IF(flag_sum .eq. 1) THEN
        write(unit_sum,'(2I5,13F10.3)') iday,time_cur,photsum,npppotsum,nppsum,resosum,     &
                                        lightsum,nee,abslightsum,precsum,tp(iday,time),     &
                                        exp(0.069*(tp(iday,time)-15.)), sumGPP, sumTER, resautsum
        photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.; precsum=0.
		sumGPP = 0.
		sumTER = 0.
        resautsum = 0.
    ENDIF
   ! output with time step of photosynthesis
    IF(flag_sum .eq. 2 .and. mod(iday,ns_pro)==0) THEN
        week = week + 1
        write(unit_sum,'(2I6,17F10.3)') week,time_cur,time_cur+(week-0.5)/52.,photsum,npppotsum,nppsum,resosum,    &
                                        lightsum,nee,abslightsum,precsum,aveT,exp(0.069*(aveT-15.)), &
                                        aet_sum, pet_sum, perc_sum, sumGPP, sumTER, resautsum
        photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.; precsum=0.
		aet_sum  = 0.; pet_sum = 0.
        perc_sum = 0.
		sumGPP = 0.
		sumTER = 0.
        resautsum = 0.
    ENDIF

    if(mod(iday,7) .eq. 0) then
        woche = woche + 1
    endif
     
    if(monthday .eq. monrec(monat)) then 
        IF(flag_sum .eq. 3 ) THEN
            tempmean = tempmean/monrec(monat)
	        if( temp_mon(monat) .le. 0.) then
		         ind_cout_mo = 12.* prec_mon(monat)
		         ind_cout_mo = 12*precsum
		    else
			     ind_cout_mo = 12.* prec_mon(monat) /(temp_mon(monat) + 10.)
			     ind_cout_mo = 12*precsum/(tempmean+10)
		    end if
		    if(temp_mon(monat) .le. 0.) then
		          ind_wiss_mo = 12.* prec_mon(monat)
		          ind_wiss_mo = 12*precsum
		    else
	             ind_wiss_mo = 12.* prec_mon(monat) /(temp_mon(monat) + 7.)
	             ind_wiss_mo = 12*precsum/(tempmean+7)
	        end if
		    if(ind_arid_mo.ne.0.) then
		       
		        ind_arid_mo = prec_mon(monat)/pet_sum
		    else
		        ind_arid_mo=0.
		    end if
    	    cwb_mo = prec_mon(monat) - pet_sum
     	    ind_cout_an = ind_cout_an + ind_cout_mo
    	    ind_wiss_an = ind_wiss_an + ind_wiss_mo
            write(unit_sum,'(I7,I5,20F10.3)') monat,time_cur,time_cur+(monat-0.5)/12.,photsum,npppotsum,nppsum,resosum,   &
                                            lightsum,nee,abslightsum, precsum, tempmean, aet_sum, pet_sum, ind_cout_mo, ind_wiss_mo, &
										    ind_arid_mo, cwb_mo, perc_sum, sumGPP, sumTER, resautsum
            photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.; precsum=0.; tempmean = 0.
		    aet_sum  = 0.; pet_sum = 0.; ind_cout_mo = 0.; ind_wiss_mo=0.; ind_arid_mo=0.; cwb_mo = 0.
            perc_sum = 0.
		    sumGPP = 0.
		    sumTER = 0.
            resautsum = 0.
        ENDIF   ! flag_sum

        monat    = monat+1
        monthday = 0
 
    endif  ! monthday
  END DO ! iday daily loop

!calculate the mean stress factor for root growth
if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
    do i=1,nlay
         do k=1,nspecies
	        svar(k)%Smean(i)=svar(k)%Smean(i)/recs(time)
         enddo
    enddo
endif
ind_shc = p_help/(t_help/10)

END SUBROUTINE stand_daily

!***************************************************************

SUBROUTINE SET_PS

  USE data_stand
  TYPE(coh_obj), POINTER :: p
  p => pt%first
  DO WHILE (ASSOCIATED(p))
! reset drought index & day counter to zero for next time step
    p%coh%drIndPS = 0.
    p%coh%nDaysPS = 0.
    p => p%next
  END DO

END SUBROUTINE SET_PS

!**************************************************************

 SUBROUTINE drought

! Calculation of drought stress indices
! Sum up of RedN
USE data_simul
USE data_stand
USE data_par
USE data_species

implicit none
integer i, ii
real, dimension(1:nspecies):: rhelp

rhelp = 0.
! drought index of trees
zeig => pt%first
do while (associated(zeig))
   ns = zeig%coh%species
  ! calculation of daily drought index
  if (zeig%coh%demand .gt. 10E-6) then
      if (ns.eq.nspec_tree+2) then                ! set drought index to 1 for mistletoe (no drought)
      zeig%coh%drIndD = 1
    else
      zeig%coh%drIndD = zeig%coh%supply / zeig%coh%demand
    endif
  else
    zeig%coh%drIndD = 1.
  endif

  select case (flag_limi)
  case (4, 5, 6, 7, 8, 9)
    rhelp(ns) = rhelp(ns) + zeig%coh%ntreeA * zeig%coh%RedNc  ! mean annual RedN 
  end select
   
  IF ((iday .ge. zeig%coh%day_bb) .AND. (iday .le. spar(zeig%coh%species)%end_bb)) THEN
     zeig%coh%drIndPS = zeig%coh%drIndPS + zeig%coh%drIndD
     zeig%coh%drIndAl = zeig%coh%drIndAl + zeig%coh%drIndD

    drIndD = drIndD + zeig%coh%ntreeA * zeig%coh%drIndD
  ENDIF

  zeig => zeig%next
enddo  ! zeig (cohorts)
if (flag_limi .ge. 4 .and. flag_limi .le. 9) then
    do i=1,anrspec 
        ii = nrspec(i)
        svar(ii)%RedN = rhelp(ii) * 10000. / (svar(ii)%sum_nTreeA * kpatchsize)    ! durch Anz. Tree pro patchsize teilen
    enddo
endif
do i=1,anrspec 
    ii = nrspec(i)
    svar(ii)%RedNm = svar(ii)%RedNm + svar(ii)%RedN
enddo
if(anz_tree.ne.0) then
    drIndD = drIndD / anz_tree
endif
END   subroutine drought

!***************************************************************

SUBROUTINE calc_fire_risk

!calculation of fire risk index
  USE data_biodiv
  USE data_climate
  USE data_simul
  USE data_soil
  USE data_species
  USE data_stand

implicit none
integer i, ii, nshelp
real hsum, hday, Tcrit_bi, cdays
real svp_13, vp_13, vpd_13, relhum_13
real k_prec      ! constant depending on precipitation
real k_phen
real  hh

if (iday.eq.1) then
    prec_flag1 = 0
    prec_flag2 = 0
    tsumrob    = 0.
    day_bb_rob = 0
    tsumbi     = 0.
    day_bb_bi  = -999.
    cdays      = 0.
    Tcrit_bi   = 0.
end if

!  calculation of day_bb for 'Robinie'
if(day_bb_rob.lt.1) then
    if(airtemp.gt.9.3) tsumrob = tsumrob + airtemp
    if(tsumrob.gt.537.) then
         day_bb_rob = iday
    end if
end if

!  calculation of day_bb for birch
nshelp = 5
! Temperature sum model Schaber 2002
if(day_bb_bi.lt.-99) then
   if(airtemp > spar(nshelp)%LTbT.and. iday.gt.47) then     
      tsumbi = tsumbi + airtemp - spar(nshelp)%LTbT
   end if
    if(tsumbi > spar(nshelp)%LTcrit) then
      day_bb_bi = iday
   end if
end if
! if birch is simulated
zeig=>pt%first
DO

   IF (.not.ASSOCIATED(zeig)) exit
      if(zeig%coh%species.eq.5) day_bb_bi  = zeig%coh%day_bb
       zeig=>zeig%next
END DO

!  fire index west
if (iday .ge. 60 .and. iday .lt. 270) then
   hday = iday/30.
   ii   = int(hday) - 1  ! month index
   hsum = SUM(clim_waterb)
   i = 1
   do i=1,4
      if (hsum .gt. risk_class(i,ii)) then
         fire_indw  = i
         fire(1)%index = i
         exit
      endif
      fire_indw  = 5
      fire(1)%index = 5
   enddo
    fd_fire_indw(fire_indw)=fd_fire_indw(fire_indw)+1
   fire(1)%frequ(fire(1)%index) = fire(1)%frequ(fire(1)%index) + 1
else
   fire(1)%index = 0
endif

if(airtemp_max .gt. -90.) then
    !  fire index east
    if (iday .ge. 46 .and. iday .lt. 275) then
       svp_13    = 6.1078 * exp(17.62 * airtemp_max / (243.12+airtemp_max)) ! saturated vapour pressure at 13.00
    ! estimation actual vapour pressure derived from mean air humidity 
       vp_13  = svp_13*hum/100
       vpd_13    = svp_13 - vp_13      ! vapour pressure deficit at 13.00 
       relhum_13 = 100. * vp_13 / svp_13

       if ((prec .ge. 1.0 .and. prec .lt. 5.0) .or. (snow_day .eq. 1)) then
          k_prec = 0.5
       else if ((prec .ge. 5.0 .and. prec .lt. 10.0) .or. (snow_day .eq. 2)) then
          k_prec = 0.25
       else if ((prec .ge. 10.0) .or. (snow_day .gt. 2)) then
          k_prec = 0.0
       else
          k_prec = 1.0
       endif

        if (iday .lt. day_bb_bi .or. day_bb_bi.eq.-999) then
          k_phen = 3.
        else if (prec.lt. 5 .and. iday .le. 227 .and. day_bb_rob.ne.0 .and. prec_flag1.eq.0) then
          k_phen = 2.
        else if (prec.ge. 5 .and. day_bb_rob.ne.0 .and. iday .gt. day_bb_rob .and. iday .lt. 227 .or. (prec_flag1.eq.1.and.iday.le.227)) then
          k_phen = 1.
          prec_flag1 = 1
       else if( day_bb_rob.eq.0) then
          k_phen = 2
       else if (iday.ge. 227.and. prec.ge. 5) then
          k_phen = 0.5
          prec_flag2 = 1
        else if(prec_flag2 .eq.1 .or. iday .gt. 243) then
          k_phen = 0.5
       else
          k_phen = 1.  ! no modification of forest fire index
       endif

      hh = (airtemp_max + 10)*vpd_13 
      fire_indi = k_prec * fire_indi + k_phen*(airtemp_max + 10)*vpd_13      
       if (fire_indi .gt. 4000) fire_indi_day = fire_indi_day + 1
       fire_indi_max = max(fire_indi, fire_indi_max)

       ! fire hazard level east
       if (fire_indi .le. 500.) then
          fire(2)%index = 1      ! no alarm level
       else if (fire_indi .le. 2000.) then
          fire(2)%index = 2      ! alarm level 1
       else if (fire_indi .le. 4000.) then
          fire(2)%index = 3      ! alarm level 2
       else if (fire_indi .le. 7000.) then
          fire(2)%index = 4      ! alarm level 3
       else
          fire(2)%index = 5      ! alarm level 4
       endif
       fire(2)%frequ(fire(2)%index) = fire(2)%frequ(fire(2)%index) + 1
    else
       fire_indi = 0.
       fire(2)%index = 0
    endif

    ! fire index Bruschek
    if (iday > 90 .AND. iday < 275) then
       if(airtemp_max .ge. 25.) Ndayshot = Ndayshot + 1  
       Psum_FP = Psum_FP + prec
    endif

    ! fire index Nesterov
    ! only calulated for vegetation and snow free period
    if (iday .ge. 60 .and. iday .lt. 275 .and. snow .lt. 0.01 .and. airtemp_max .gt. 0.) then
        if (prec .lt. 3.) then
            day_nest = day_nest + 1
            p_nest = p_nest + (airtemp_max -  dptemp) * airtemp_max 
        else
            day_nest = 0
            p_nest = 0.
        endif
       if (p_nest .le. 300.) then
          fire(3)%index = 1      ! minimal
       else if (p_nest .le. 1000.) then
          fire(3)%index = 2      ! moderate
       else if (p_nest .le. 4000.) then
          fire(3)%index = 3      ! high
       else 
          fire(3)%index = 4      ! extreme
       endif
       fire(3)%frequ(fire(3)%index) = fire(3)%frequ(fire(3)%index) + 1
    else
       p_nest = 0.
       fire(3)%index = 0
    endif
else
    fire(2)%index = -99.0
    fire(3)%index = -99.0
endif ! airtemp_max          

END   subroutine calc_fire_risk

!*******************************************************************************

subroutine calc_frost_index

USE data_frost
USE data_climate
USE data_simul
USE data_stand

implicit none
integer     :: day_bb, j, t, m, ii

!  absolute and annual last frost day during spring/ summer
if(airtemp_min .lt. temp_frost .and. iday .lt. 200 ) then
    if(iday.gt.dlfabs ) dlfabs = iday 
    if(iday.gt.date_lftot(time)) date_lftot(time)=iday
end if    

! annual number of frost days after start of the vegetation period and  annual last frost day
if(flag_vegper.eq.1. .and. iday.lt.200) then
    if(airtemp_min .lt. temp_frost) then
        dnlf(time) = dnlf(time) +1
! calculation  of last frost day after beginning of vegetation period due to 5°C threshold for the case of needle trees
       if( waldtyp.eq.10 .or. waldtyp.eq.40.or.waldtyp.eq.90 .and. iday.gt. date_lf(time)) date_lf(time)= iday
    end if
 end if
! calculation of the number of the actual month
  j= time_cur
  ii = iday
  call tzinda(t,m,j,ii)
  iday = ii
 
 if(m.eq.4 .or. m.eq.5 .or. m.eq.6) then
    if(airtemp_min .lt.0) then
        anzdlf(time)=anzdlf(time)+1
        sumtlf(time) = sumtlf(time) + airtemp_min
    end if
 endif
 ! annual minimum temperature may for year time
 if(airtemp_min.lt.tminmay_ann(time).and. m.eq.5) tminmay_ann(time) = airtemp_min
 ! absolute minimum temperature May

 if( airtemp_min .lt. tminmay .and. m.eq.5) tminmay = airtemp_min
 ! assuming mono species stand !!!
 zeig=>pt%first
DO
   IF (.not.ASSOCIATED(zeig)) exit
      taxnum= zeig%coh%species
      day_bb  = zeig%coh%day_bb
      exit
      zeig=>zeig%next
END DO

! caculation not for conifer stands (pine, spruce, douglas fir)
if(waldtyp .ne. 10 .and. waldtyp .ne. 40 .and. waldtyp .ne.90)then
   if(all_leaves_on.eq.1) then  
     if (iday.ge.day_bb .and. iday.lt.200) then
! calculation of number of frost day during vegetation period (bud burst) for year time
       if(airtemp_min .lt. temp_frost ) then
         dnlf_sp(time) = dnlf_sp(time) +1
   ! calculagtion of last frost day after beginning of vegetation period by bud burst
        if(iday .gt. date_lf(time)) date_lf(time)= iday

       end if
    end if
  end if   ! all_leaves_on
end if  ! waldtyp
  
END subroutine calc_frost_index

!*******************************************************************************

Subroutine calc_endbb

use data_climate
use data_stand
use data_species
use data_simul

implicit none
integer    :: tax,fl

if(iday.gt.180) then
zeig => pt%first
do while (associated(zeig))
   tax = zeig%coh%species
   fl = zeig%coh%flag_vegend
   if(spar(tax)%end_bb.ne.366) then
    if(spar(ns)%flag_endbb.eq.0) then
      if(airtemp.ge.5. .and. fl .ne.0)  then
       fl=0
      else if(airtemp.lt.5. .and. fl.eq.0) then
        fl =1
      else if(airtemp.lt.5. .and. fl.eq.1) then
        fl =2
      else if(airtemp.lt.5. .and. fl.eq.2) then
        fl =3
      else  if(airtemp.lt.5. .and. fl.eq.3)then
         fl =4
      else if(airtemp.lt.5. .and. fl.eq.4) then
         fl =5
      end if
      zeig%coh%flag_vegend = fl
     if(fl .eq.5) then
       spar(tax)%flag_endbb=1
       spar(tax)%end_bb = iday
       write(666,*) time, iday
     end if
   end if
   zeig => zeig%next
  end if
end do
end if
end subroutine calc_endbb
