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