Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • foresee/4C
  • gutsch/4C
2 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 6355 additions and 0 deletions
source_code/version_2.3_windows/sorting.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for: *!
!* - dimsort: sorting of cohorts according to a *!
!* charcteristic variable *!
!* - sort2: subroutine from Numerical recipes *!
!* *!
!* 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 dimsort(n,var,ranktable)
USE data_species
USE data_stand ! state variables of stand, cohort and cohort element
IMPLICIT NONE
INTEGER :: isort,n
INTEGER :: ranktable(n)
CHARACTER(3) :: var
REAL :: sortarray(n)
ranktable=0
sortarray=0
isort=1
zeig=>pt%first
DO
IF (.not.ASSOCIATED(zeig)) exit
IF (zeig%coh%species .le. nspec_tree) THEN ! for trees only
ranktable(isort) = zeig%coh%ident
IF(var=='hei') sortarray(isort) = zeig%coh%height
IF(var=='dbh') sortarray(isort) = zeig%coh%diam
isort=isort+1
ENDIF
zeig=>zeig%next
END DO
CALL sort2(n,sortarray,ranktable)
END SUBROUTINE dimsort
!******************************************************************************
SUBROUTINE sort2(n,arr,brr)
! sorts array arr(1:n) into an ascending order and
! makes the corresponding rearrangement of the array brr(1:n)
INTEGER n,M,NSTACK
REAL arr(n)
INTEGER brr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,b,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
b=brr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
brr(i+1)=brr(i)
11 continue
i=0
2 arr(i+1)=a
brr(i+1)=b
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
temp=brr(k)
brr(k)=brr(l+1)
brr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
temp=brr(l+1)
brr(l+1)=brr(ir)
brr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
temp=brr(l)
brr(l)=brr(ir)
brr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
temp=brr(l+1)
brr(l+1)=brr(l)
brr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
b=brr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
temp=brr(i)
brr(i)=brr(j)
brr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
brr(l)=brr(j)
brr(j)=b
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END Subroutine
! (C) Copr. 1986-92 Numerical Recipes Software "!D#+.
source_code/version_2.3_windows/stand_bal.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for: *!
!* - STAND_BALANCE: Recalculation of stand variables *!
!* contains: *!
!* UPDATE_AGE *!
!* - STAND_BAL_SPEC *!
!* - CLASS *!
!* - CLASST *!
!* - CLASS_MAN *!
!* - CALC_HEIDOM *!
!* - MAX_HEIGHT(nrmax,anz,cohl) *!
!* - STANDUP: Update of cover and ceppot *!
!* - LITTER: Summation variables of litter fractions *!
!* - calc_ind_rep: calculation of representation index *!
!* - overstorey *!
!* *!
!* 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_balance
use data_species
use data_stand
use data_climate
use data_simul
use data_manag
use data_out
use data_par
implicit none
integer i, ntr, nd, hanz
integer, dimension(nspecies) :: helpin, helpout
real conv ! conversion factor
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' stand_balance'
if(time>0. .and. flag_standup.ne.2) call update_age
! calculation of total dead biomass per cohort and total biomass of allcohorts
! calc. of ceppot
anz_sveg = 0
anz_tree = 0.
anz_tree_in = 0.
anz_tree_out = 0.
anz_spec_in = 0.
anz_spec_out = 0.
anz_coh_in = 0.
anz_coh_out = 0.
anz_coh_act = 0.
lai_in = 0.
lai_out = 0.
totfol_in = 0.
totfol_out = 0.
med_diam_in = 0.
med_diam_out = 0.
hmean_in = 0.
hmean_out = 0.
mean_height = 0.
sumbio = 0.
sumbio_in = 0.
sumbio_out = 0.
sumNPP = 0.
drIndAl = 0.
Ndem = 0.
helpin = 0
helpout = 0
basal_area = 0.
totstem_m3 = 0.
totsteminc_m3 = 0.
totsteminc = 0
autresp = 0.
totfol = 0.
totsap = 0.
totfrt = 0.
totfrt_p = 0.
totcrt = 0.
tottb = 0.
tothrt = 0.
sumbio_sv = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
ns = zeig%coh%species
ntr = zeig%coh%ntreeA
svar(ns)%daybb = zeig%coh%day_bb
if(ns.le.nspec_tree) then
if(zeig%coh%ident .le. coh_ident_max) then
anz_coh_act = anz_coh_act + 1
anz_tree = anz_tree + ntr
zeig%coh%totBio = zeig%coh%x_fol + zeig%coh%x_sap + zeig%coh%x_hrt + zeig%coh%x_tb + zeig%coh%x_frt +zeig%coh%x_crt
zeig%coh%Dbio = zeig%coh%nTreeD * zeig%coh%totBio
sumbio = sumbio + ntr * zeig%coh%totBio
sumNPP = sumNPP + ntr * zeig%coh%NPP
Ndem = Ndem + ntr * zeig%coh%Ndemc_c
select case (flag_dis)
case (0,1)
autresp = autresp + ntr * zeig%coh%maintres
case (2)
autresp = autresp + ntr * (zeig%coh%maintres+zeig%coh%biocost_all)
end select
totfol = totfol + ntr * zeig%coh%x_fol
totsap = totsap + ntr * zeig%coh%x_sap
totfrt = totfrt + ntr * zeig%coh%x_frt
totcrt = totcrt + ntr * zeig%coh%x_crt
tottb = tottb + ntr * zeig%coh%x_tb
tothrt = tothrt + ntr * zeig%coh%x_hrt
if (zeig%coh%height.le.thr_height) then
seedlfrt = seedlfrt + zeig%coh%x_frt * ntr
endif
totstem_m3 = totstem_m3 + (ntr*zeig%coh%x_sap + ntr*zeig%coh%x_hrt) &
/(spar(ns)%prhos*1000000) ! conversion kg/patch ---m³/ha
nd = zeig%coh%nDaysGr
if (nd .gt. 0) drIndAl = drIndAl + ntr * zeig%coh%drIndAl * zeig%coh%NPP / nd
endif
if(zeig%coh%ident > coh_ident_max) then
anz_tree_in = anz_tree_in + ntr
sumbio_in = sumbio_in + ntr * zeig%coh%totBio
anz_coh_in = anz_coh_in + 1
helpin(ns) = ns
lai_in = lai_in + ntr * zeig%coh%t_leaf/kpatchsize
totfol_in = totfol_in + ntr * zeig%coh%x_fol
med_diam_in = med_diam_in + ntr * (zeig%coh%diam**2)
hmean_in = hmean_in + ntr * zeig%coh%height
totfrt = totfrt + ntr * zeig%coh%x_frt
endif
if((zeig%coh%nTreeD > 0.1) .or. (zeig%coh%nTreeM > 0.1) .or. (zeig%coh%nTreet > 0.1)) then
hanz = zeig%coh%nTreeD + zeig%coh%nTreeM + zeig%coh%nTreet
anz_tree_out = anz_tree_out + hanz
sumbio_out = sumbio_out + hanz * zeig%coh%totBio
sumNPP = sumNPP + hanz * zeig%coh%NPP ! eliminated (died or harvested) trees produce during the year as well;
autresp = autresp + hanz * zeig%coh%maintres
anz_coh_out = anz_coh_out + 1
helpout(ns) = ns
lai_out = lai_out + hanz * zeig%coh%t_leaf/kpatchsize
totfol_out = totfol_out + hanz * zeig%coh%x_fol
med_diam_out = med_diam_out + hanz * (zeig%coh%diam**2)
hmean_out = hmean_out + hanz * zeig%coh%height
endif
else
ntr = zeig%coh%ntreeA
anz_sveg = anz_sveg +1
zeig%coh%totBio = zeig%coh%x_fol + (1.+spar(ns)%alphac)*(zeig%coh%x_sap + zeig%coh%x_hrt) + zeig%coh%x_frt
sumbio_sv = sumbio_sv + ntr * zeig%coh%totBio
totfrt_p = totfrt_p + ntr * zeig%coh%x_frt
end if !only trees
zeig=>zeig%next
end do
if (flag_cumNPP .eq. 1) then
cum_sumNPP = cum_sumNPP + sumNPP
flag_cumNPP = 0
endif
if (sumNPP .gt. 1E-06) drIndAl = drIndAl / sumNPP
! conversion kg/patch ---> kg/ha; N/patch ---> N/ha
conv = 10000./kpatchsize
totfrt_p = totfrt_p + totfrt ! Rootmass f. patch (trees and soil veg.) save before conversion; Wurzelmenge vor Umrechnung sichern
if (totfrt_p .gt. zero) then
totfrt_1 = 1./totfrt_p ! reciprocal for later calculationshKehrwert f. spaetere Berechnungen
else
totfrt_1 = 0.
endif
totfrt = totfrt * conv
totfol = totfol * conv
totfol_in = totfol_in * conv
totfol_out = totfol_out * conv
tottb = tottb * conv
totsap = totsap * conv
tothrt = tothrt * conv
totcrt = totcrt * conv
sumbio = sumbio * conv
sumbio_in = sumbio_in * conv
sumbio_out = sumbio_out * conv
sumbio_sv = sumbio_sv * conv
Ndem = Ndem / kpatchsize ! g per tree --> g/m2
totstem_m3 = totstem_m3* conv ! m3/ha
anz_tree_ha = anz_tree * conv
anz_tree_in = anz_tree_in * conv
anz_tree_out = anz_tree_out * conv
do i=1, nspec_tree+1 ! for all but mistletoe
if (helpin(i) > 0) anz_spec_in = anz_spec_in + 1
if (helpout(i) > 0) anz_spec_out = anz_spec_out + 1
enddo
if(anz_tree_in > 0.) then
med_diam_in = sqrt(med_diam_in/anz_tree_in)
hmean_in = hmean_in / anz_tree_in
endif
if(anz_tree_out > 0.) then
med_diam_out = sqrt(med_diam_out/anz_tree_out)
hmean_out = hmean_out / anz_tree_out
endif
! call species values
call classt
call stand_bal_spec
call calc_ind_rep
!call classification of stand diameter and height
call class
! moving of understorey tree cohorts to overstorey tree cohorts
if(flag_mg.ne.33) call overstorey
contains
!**************************************************************
subroutine update_age
if(flag_standup.ne. 2) then
zeig=>pt%first
do
if(.not.associated(zeig)) exit
zeig%coh%x_age=zeig%coh%x_age + 1
zeig=>zeig%next
end do
end if
end subroutine update_age
end subroutine stand_balance
!**************************************************************
subroutine stand_bal_spec
use data_climate
use data_out
use data_simul
use data_site
use data_stand
use data_species
use data_par
use data_manag
implicit none
integer :: i, j, k, ntr, nd, lowtree, hntr, spec_new
real,dimension(nspec_tree):: vgldom1, vgldom2, vgldom_spec1, vgldom_spec2
integer,dimension(nspec_tree):: anzdom1, anzdom2, anzdom_spec1, anzdom_spec2, &
helpdiam
integer,dimension(nspecies):: helpanz
integer helpntr
integer help_nr_inf_trees
logical lhelp
INTEGER leapyear
real atemp, hh, help_height_top
real triangle
real, external :: daylength
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' stand_bal_spec'
! Initialisation
vgldom1 = 0.
vgldom2 = 0.
anzdom1 = 0
anzdom2 = 0
med_diam = 0.
mean_diam = 0.
mean_height = 0.
hdom = 0. ! dominante height (highest or two highest cohorst); Hoehe (hoehste oder die beiden hoechsten Kohorten)
anz_spec = 0 ! currently existing species
lowtree = 0 ! amount of trees with DBH=0 for the whole population; Anzahl Baeume mit DBH = 0 fuer gesamten Bestand
hntr = 0
helpanz = 0 ! auxiliary variable to count species; Hilfsvariable um Spezies zu zaehlen
helpdiam = 0 ! amount of trees with DBH=0 per species; Anzahl Baeume mit DBH = 0 pro Spezies
vgldom_spec1 = 0.
vgldom_spec2 = 0.
anzdom_spec1 = 0
anzdom_spec2 = 0
svar%med_diam = 0.
svar%mean_diam = 0.
svar%mean_jrb = 0.
svar%mean_height= 0.
svar%basal_area = 0.
svar%sumNPP = 0.
svar%Ndem = 0.
svar%Nupt = 0.
svar%sum_ntreea = 0.
svar%sum_ntreed = 0.
svar%sum_bio = 0.
svar%sum_lai = 0.
svar%anz_coh = 0.
svar%drIndAl = 0.
svar%totsteminc = 0.
svar%totsteminc_m3 = 0.
svar%fol = 0.
svar%sap = 0.
svar%hrt = 0.
svar%frt = 0.
!update height of mistletoe to height of mistletoe-infected cohort
if (flag_mistle.ne.0) then
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%mistletoe.eq.1) then
help_height_top=zeig%coh%height
end if
zeig => zeig%next
ENDDO
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%height = help_height_top !upper crown
zeig%coh%x_hbole = zeig%coh%height-50. !lower crown
end if
zeig => zeig%next
ENDDO
end if ! end update of height of Mistletoe
! update of #of mistletoe upon dist_manag
if (flag_mistle.ne.0) then
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%mistletoe.eq.1) then
help_nr_inf_trees=zeig%coh%nTreeA
end if
zeig => zeig%next
ENDDO
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%nTreeA= help_nr_inf_trees*AMAX1(1.,dis_rel(time))
zeig%coh%nta=zeig%coh%nTreeA
end if
zeig => zeig%next
ENDDO
end if ! end update #of mistletoe
zeig => pt%first
DO WHILE (ASSOCIATED(zeig))
ns = zeig%coh%species
helpanz(ns) = helpanz(ns) + 1 ! all species incl. ground vegetation;
IF(zeig%coh%ident .le. coh_ident_max) THEN
ntr = zeig%coh%ntreea
IF(ns .le. nspec_tree) THEN
!
IF((ns .le. nspec_tree) .and. (ntr > 0.) .and. (zeig%coh%diam > 0.)) THEN
svar(ns)%med_diam = svar(ns)%med_diam + ntr * (zeig%coh%diam**2)
med_diam = med_diam + ntr * (zeig%coh%diam**2)
mean_diam = mean_diam + ntr*zeig%coh%diam
svar(ns)%mean_diam = svar(ns)%mean_diam + ntr*zeig%coh%diam
svar(ns)%mean_height = svar(ns)%mean_height + ntr*zeig%coh%height
svar(ns)%mean_jrb = svar(ns)%mean_jrb + ntr*zeig%coh%jrb
mean_height = mean_height + ntr*zeig%coh%height
hntr = hntr + ntr
ELSE
! Trees with DBH=0 for population and per species; Baeume mit DBH =0 fuer Bestand und pro Spezies
lowtree = lowtree + ntr
helpdiam(ns) = helpdiam(ns) + ntr
ENDIF ! ns
IF(zeig%coh%height > vgldom1(ns)) THEN
vgldom2(ns) = vgldom1(ns)
anzdom2(ns) = anzdom1(ns)
vgldom1(ns) = zeig%coh%height
anzdom1(ns) = ntr
ELSE
if(zeig%coh%height > vgldom2(ns))then
vgldom2(ns) = zeig%coh%height
anzdom2(ns) = ntr
endif
ENDIF ! vgldom1
IF(zeig%coh%height > vgldom_spec1(ns)) THEN
vgldom_spec2(ns) = vgldom_spec1(ns)
anzdom_spec2(ns) = anzdom_spec1(ns)
vgldom_spec1(ns) = zeig%coh%height
anzdom_spec1(ns) = ntr
ELSE
if(zeig%coh%height > vgldom_spec2(ns))then
vgldom_spec2(ns) = zeig%coh%height
anzdom_spec2(ns) = ntr
endif
ENDIF ! vgldom_spec2
ELSE
svar(ns)%dom_height = zeig%coh%height
ENDIF ! end loop across trees
svar(ns)%sumNPP = svar(ns)%sumNPP + ntr * zeig%coh%NPP
svar(ns)%sum_ntreea = svar(ns)%sum_ntreea + ntr
svar(ns)%sum_ntreed = svar(ns)%sum_ntreed + zeig%coh%nTreeD + zeig%coh%nTreeM ! died or harvested trees of current year; ausgeschiedene Bäume des akt. Jahres
svar(ns)%Ndem = svar(ns)%Ndem + ntr * zeig%coh%Ndemc_c
svar(ns)%Nupt = svar(ns)%Nupt + ntr * zeig%coh%Nuptc_c
svar(ns)%sum_bio = svar(ns)%sum_bio + ntr * zeig%coh%totBio
svar(ns)%sum_lai = svar(ns)%sum_lai + ntr * zeig%coh%t_leaf/kpatchsize
svar(ns)%anz_coh = svar(ns)%anz_coh + 1
svar(ns)%totsteminc = svar(ns)%totsteminc + ntr * zeig%coh%stem_inc
if (zeig%coh%species.ne.nspec_tree+2) then !no stem increment for mistletoe
svar(ns)%totsteminc_m3 = svar(ns)%totsteminc_m3 + ntr * zeig%coh%stem_inc /(spar(ns)%prhos*1000000)
endif
svar(ns)%fol = svar(ns)%fol + ntr * zeig%coh%x_fol
svar(ns)%sap = svar(ns)%sap + ntr * zeig%coh%x_sap
svar(ns)%hrt = svar(ns)%hrt + ntr * zeig%coh%x_hrt
svar(ns)%frt = svar(ns)%frt + ntr * zeig%coh%x_frt
nd = zeig%coh%nDaysGr
if (nd .gt. 0) svar(ns)%drIndAl = svar(ns)%drIndAl + ntr * zeig%coh%drIndAl * zeig%coh%NPP / nd
ENDIF ! coh%ident
zeig%coh%ntreed = 0.
zeig%coh%ntreem = 0.
zeig => zeig%next
ENDDO ! cohort loop
! neue Spezies feststellen und belegen
if (time .gt. 1) then
do i=1,nspecies
if (helpanz(i) > 0) then
spec_new = 0
lhelp = .True.
do j=1,anrspec
if (nrspec(j) .eq. i) lhelp = .False.
enddo
if (lhelp) then
spec_new = i
if(spec_new.le.nspec_tree) then
IF(spar(spec_new)%Phmodel==1) THEN
svar(spec_new)%Pro = 0.
svar(spec_new)%Inh = 1.
ELSE
svar(spec_new)%Pro = 0.
svar(spec_new)%Inh = 0.
svar(spec_new)%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, time)
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
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
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 + 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
IF(spar(spec_new)%phmodel==4) THEN
svar(spec_new)%daybb = svar(spec_new)%ext_daybb
ELSE
svar(spec_new)%daybb = 181 + leapyear(time_cur)
ENDIF
endif
endif
endif
enddo
endif ! time
k = 0
do i=1,nspecies
if (helpanz(i) > 0) then
k = k + 1
anrspec = k
nrspec(k) = i
endif
ntr = svar(i)%sum_ntreea
if (svar(i)%sumNPP .gt. 1E-06) svar(i)%drIndAl = svar(i)%drIndAl / svar(i)%sumNPP
if (i .le. nspec_tree) then
IF(helpanz(i) > 0) THEN
anz_spec = anz_spec + 1
IF(helpdiam(i) < ntr) THEN
svar(i)%med_diam = SQRT(svar(i)%med_diam / (ntr - helpdiam(i)))
ENDIF
svar(i)%Ndem = svar(i)%Ndem / kpatchsize ! g per tree --> g/m2
svar(i)%Nupt = svar(i)%Nupt / kpatchsize ! g per tree --> g/m2
if (ntr .ne. 0) then
svar(i)%mean_height = svar(i)%mean_height / ntr
svar(i)%mean_diam = svar(i)%mean_diam / ntr
svar(i)%mean_jrb = svar(i)%mean_jrb / ntr
svar(i)%basal_area = pi*(ntr-helpdiam(i))*(svar(i)%med_diam*svar(i)%med_diam/40000)*10000/kpatchsize
else
svar(i)%sum_ntreea = 0.
endif
call calc_heidom_spec(i)
ENDIF
end if ! nspec_tree
! conversion kg/patch ---> kg/ha; N/patch ---> N/ha
helpntr = svar(i)%sum_nTreeA* 10000./kpatchsize
if(helpntr.eq.0 .and. svar(i)%sum_nTreeA.eq.1) then
svar(i)%sum_nTreeA = 1
else
svar(i)%sum_nTreeA = helpntr
end if
svar(i)%sum_bio = svar(i)%sum_bio * 10000./kpatchsize
svar(i)%fol = svar(i)%fol * 10000./kpatchsize
svar(i)%sap = svar(i)%sap* 10000./kpatchsize
svar(i)%hrt= svar(i)%hrt* 10000./kpatchsize
svar(i)%frt= svar(i)%frt* 10000./kpatchsize
svar(i)%totstem_m3= ( svar(i)%sap + svar(i)%hrt)/ (spar(i)%prhos*1000000) ! m3/ha
svar(i)%totsteminc = svar(i)%totsteminc * 10000./kpatchsize ! kg/ha
svar(i)%totsteminc_m3 = svar(i)%totsteminc_m3 * 10000./kpatchsize ! kg/ha
totsteminc_m3 = totsteminc_m3 + svar(i)%totsteminc_m3
totsteminc = totsteminc + svar(i)%totsteminc
end do
! new calculation of dominant height
call calc_heidom
if(anz_tree>0)then
if(lowtree<anz_tree) then
med_diam = sqrt(med_diam /(anz_tree-lowtree))
mean_diam = mean_diam /(anz_tree-lowtree)
mean_height = mean_height /(anz_tree-lowtree)
basal_area = pi*(anz_tree-lowtree)*(med_diam*med_diam/40000)*10000/kpatchsize
endif
else
if (hntr .ne. 0) then
med_diam = sqrt(med_diam /hntr)
mean_diam = mean_diam / hntr
mean_height = mean_height / hntr
else
med_diam = 0.
mean_diam = 0.
mean_height = 0.
endif
endif
end subroutine stand_bal_spec
!**************************************************************
subroutine class
use data_stand
use data_simul
use data_species
use data_par
implicit none
integer i,k
diam_class=0;height_class=0
diam_class_age=0.
diam_class_h = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
k = zeig%coh%species
if (k.ne.nspec_tree+2) then !exclusion of mistletoe
if(zeig%coh%diam<=dclass_w) then
diam_class(1,k)=diam_class(1,k)+zeig%coh%ntreea
diam_class_h(1,k) = diam_class_h(1,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(1,k) = diam_class_age(1,k)+zeig%coh%x_age*zeig%coh%ntreea
end if
do i=2,num_class
if(zeig%coh%diam.le.(dclass_w + dclass_w*(i-1)) .and. zeig%coh%diam>(dclass_w + dclass_w*(i-2))) then
diam_class(i,k)=diam_class(i,k) + zeig%coh%ntreea
diam_class_h(i,k) = diam_class_h(i,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(i,k) = diam_class_age(i,k)+zeig%coh%x_age*zeig%coh%ntreea
else if(zeig%coh%diam.gt. (dclass_w + dclass_w*(num_class-2))) then
diam_class(num_class,k)=diam_class(num_class,k) + zeig%coh%ntreea
diam_class_h(num_class,k) = diam_class_h(num_class,k) + zeig%coh%ntreea*zeig%coh%height
diam_class_age(num_class,k) = diam_class_age(num_class,k) + zeig%coh%x_age+zeig%coh%ntreea
exit
end if
enddo
if(zeig%coh%height<=100) height_class(1) = height_class(1)+zeig%coh%ntreea
if(zeig%coh%height>100.and.zeig%coh%height<500) height_class(2) = height_class(2)+zeig%coh%ntreea
do i=3,num_class-2
if(zeig%coh%height>(i+2)*100.and.zeig%coh%height<=(i+3)*100) height_class(i) = height_class(i)+zeig%coh%ntreea
enddo
if(zeig%coh%height>5000.and.zeig%coh%height<5500) height_class(num_class-1) = height_class(num_class-1)+zeig%coh%ntreea
if(zeig%coh%height>5500) height_class(num_class) = height_class(num_class)+zeig%coh%ntreea
endif!exclusion of mistletoe
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspec_tree
if(diam_class(i,k).ne.0) diam_class_h(i,k) = (diam_class_h(i,k)/diam_class(i,k))*10000./kpatchsize
if(diam_class_age(i,k).ne.0.and.diam_class(i,k).ne.0 ) diam_class_age(i,k) =diam_class_age(i,k)/diam_class(i,k)
diam_class(i,k) = diam_class(i,k)*10000./kpatchsize
end do
end do
end subroutine class
!**************************************************************
subroutine classt
use data_stand
use data_simul
use data_species
implicit none
integer i,k
diam_class_t=0;height_class=0
diam_class_h = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
k = zeig%coh%species
if (k.ne.nspec_tree+2) then ! exclusion mistletoe
if(zeig%coh%diam<=dclass_w) then
diam_class_t(1,k)=diam_class_t(1,k)+zeig%coh%ntreed
end if
do i=2,num_class
if(zeig%coh%diam.le.(dclass_w + dclass_w*(i-1)) .and. zeig%coh%diam>(dclass_w + dclass_w*(i-2))) then
diam_class_t(i,k)=diam_class_t(i,k) + zeig%coh%ntreed
else if(zeig%coh%diam.gt. (dclass_w + dclass_w*(num_class-2))) then
diam_class_t(num_class,k)=diam_class_t(num_class,k) + zeig%coh%ntreed
exit
end if
enddo
endif !exclusion of mistletoe
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspec_tree
diam_class_t(i,k)=diam_class_t(i,k)*10000./kpatchsize
end do
end do
end subroutine classt
!**************************************************************
subroutine class_man
use data_stand
use data_simul
use data_species
use data_manag
implicit none
integer i , k
real anz
diam_classm=0
diam_classm_h=0.
diam_class_mvol = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem.ne.0.or.(zeig%coh%ntreed.gt.0 .and. zeig%coh%diam.gt.tardiam_dstem)) then
if(zeig%coh%diam.le.tardiam_dstem) then
anz = zeig%coh%ntreem
else
anz = zeig%coh%ntreem + zeig%coh%ntreed
end if
k = zeig%coh%species
if(zeig%coh%diam<=dclass_w) then
diam_classm(1,k)=diam_classm(1,k)+anz
diam_classm_h(1,k) = diam_classm_h(1,k) + anz*zeig%coh%height
diam_class_mvol(1,k) = diam_class_mvol(1,k) +anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
if(zeig%coh%diam<=dclass_w*2.and.zeig%coh%diam.gt.dclass_w) then
diam_classm(2,k)=diam_classm(2,k)+anz
diam_classm_h(2,k) = diam_classm_h(2,k) + anz*zeig%coh%height
diam_class_mvol(2,k) = diam_class_mvol(2,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
do i=3,num_class
if(zeig%coh%diam.le.(dclass_w*2 + dclass_w*(i-2)) .and. zeig%coh%diam>(dclass_w*2 + dclass_w*(i-3))) then
diam_classm(i,k) = diam_classm(i,k) + anz
diam_classm_h(i,k) = diam_classm_h(i,k) + anz*zeig%coh%height
diam_class_mvol(i,k) = diam_class_mvol(i,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
else if(zeig%coh%diam.gt. (dclass_w*2 + dclass_w*(num_class-3))) then
diam_classm(num_class,k)=diam_classm(num_class,k) + anz
diam_classm_h(num_class,k) = diam_classm_h(num_class,k) + anz*zeig%coh%height
diam_class_mvol(num_class,k) = diam_class_mvol(num_class,k) + anz*(zeig%coh%x_sap + zeig%coh%x_hrt)
end if
enddo
end if
zeig=>zeig%next
enddo
do i=1,num_class
do k=1,nspecies
if(diam_classm(i,k).ne.0) diam_classm_h(i,k) = diam_classm_h(i,k)/diam_classm(i,k)
if(diam_class_mvol(i,k).ne.0.) then
diam_class_mvol(i,k) = diam_class_mvol(i,k)/(spar(k)%prhos*1000000)*10000/kpatchsize
end if
diam_classm(i,k) = diam_classm(i,k)*10000./kpatchsize
end do
end do
end subroutine class_man
!**************************************************************
subroutine calc_heidom
use data_out
use data_simul
use data_stand
implicit none
real :: mh
integer :: nhelp, &
nh1,nh2, &
testflag=0, &
j
allocate (height_rank(anz_coh))
nh1=0
nh2=0
mh = 0
testflag = 0
nhelp = nint(kpatchsize/100)
if(anz_tree.le.nhelp) nhelp = anz_tree
! sorting by height of cohorts
call dimsort(anz_coh, 'height',height_rank)
if(anz_tree>1) then
do j=anz_coh, 1,-1
call dimsort(anz_coh, 'height',height_rank)
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ident.eq.height_rank(j)) then
nh2 = nh1
nh1 = nh1 + zeig%coh%ntreea
if(nh1.le. nhelp) then
mh = mh + zeig%coh%ntreea*zeig%coh%height
else
mh = mh + zeig%coh%height*( nhelp - nh2)
testflag=1
exit
end if
endif
zeig=>zeig%next
if(testflag.eq.1) exit
end do
if(testflag.eq.1) exit
if(nh1.eq.nhelp) exit
end do
if (nhelp.lt. nh1) then
hdom = mh/nhelp
else
hdom = mh/nh1
end if
end if
deallocate(height_rank)
end subroutine calc_heidom
!**************************************************************
subroutine calc_heidom_spec(ispec)
!*****************************************************
! species specific dominant height calculation
!*****************************************************
use data_out
use data_simul
use data_stand
implicit none
real :: mh
integer :: nhelp, &
nh1,nh2, &
testflag=0, &
j, &
ispec
allocate (height_rank(anz_coh))
hdom = 0
nh1=0
nh2=0
mh = 0
testflag = 0
! calculation of number of trees used for H100 ( 100/ ha = nhelp/ kpachtsize)
nhelp = nint(kpatchsize/100)
if(anz_tree.le.nhelp) nhelp = anz_tree
! sorting by height of cohorts
call dimsort(anz_coh, 'height',height_rank)
if(anz_tree>1) then
do j=anz_coh, 1,-1
call dimsort(anz_coh, 'height',height_rank)
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ident.eq.height_rank(j).and. zeig%coh%species.eq.ispec) then
nh2 = nh1
nh1 = nh1 + zeig%coh%ntreea
if(nh1.le. nhelp) then
mh = mh + zeig%coh%ntreea*zeig%coh%height
else
mh = mh + zeig%coh%height*( nhelp - nh2)
testflag=1
exit
end if
endif
zeig=>zeig%next
if(testflag.eq.1) exit
end do
if(testflag.eq.1) exit
if(nh1.eq.nhelp) exit
end do
if (nhelp.lt. nh1.and. nhelp.ne.0) then
hdom = mh/nhelp
else if(nh1.ne.0) then
hdom = mh/nh1
end if
else if(anz_tree.eq.1) then
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.ispec) hdom=zeig%coh%height
zeig=>zeig%next
end do
end if
deallocate(height_rank)
svar(ispec)%dom_height = hdom
end subroutine calc_heidom_spec
!**************************************************************
subroutine max_height(nrmax,anz,cohl)
use data_out
use data_simul
use data_stand
implicit none
integer :: nrmax
integer :: nrmax_h
integer :: anz, testflag,i
real :: help_h1, help_h2
integer,dimension(0:anz_coh) :: cohl
testflag=0
nrmax = -1
nrmax_h = -1
help_h2=0.
help_h1=0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
do i=0,anz-1
if(cohl(i).eq.zeig%coh%ident) then
testflag=1
endif
end do
if (testflag.eq.0) then
help_h2= zeig%coh%height
nrmax_h = zeig%coh%ident
if(help_h2.gt. help_h1) then
help_h1 = help_h2
nrmax = nrmax_h
end if
end if
zeig=>zeig%next
testflag = 0
end do
anz = anz +1
cohl(anz-1) = nrmax
end subroutine max_height
!**************************************************************
SUBROUTINE standup
! update of stand variables (LAI, cover, waldtyp)
USE data_par
USE data_stand
USE data_soil
USE data_species
use data_out
use data_simul
implicit none
integer i
REAL :: sumfol_can = 0.
REAL :: sumfol_sveg= 0.
REAL :: ntr, cover3
! estimating degree of covering
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time, ' standup'
cover3 = 0.
sumfol_can = 0.
sumfol_sveg= 0.
crown_area = 0.
do i = 1, anrspec
svar(nrspec(i))%crown_area = 0.
enddo
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
if (zeig%coh%crown_area .ge. 0) then
ntr = zeig%coh%nTreeA
ns = zeig%coh%species
cover3 = cover3 + ntr * zeig%coh%crown_area
if (ns .le. nspec_tree) then
sumfol_can = sumfol_can + ntr * zeig%coh%x_fol
crown_area = crown_area + ntr * zeig%coh%crown_area
else
sumfol_sveg = sumfol_sveg + ntr * zeig%coh%x_fol
endif
svar(ns)%crown_area = svar(ns)%crown_area + ntr * zeig%coh%crown_area
endif
zeig=>zeig%next
end do
cover3 = cover3 / kpatchsize
anz_tree = 0
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
ns=zeig%coh%species
if (ns .le. nspec_tree) then
zeig%coh%rel_fol = zeig%coh%ntreea * zeig%coh%x_fol/sumfol_can
ceppot_can = ceppot_can + zeig%coh%rel_fol * spar(ns)%ceppot_spec
anz_tree = anz_tree + zeig%coh%ntreea
else if (ns.eq.nspec_tree+1) then
zeig%coh%rel_fol = zeig%coh%ntreea * zeig%coh%x_fol/sumfol_sveg
ceppot_sveg = ceppot_sveg + zeig%coh%rel_fol * spar(ns)%ceppot_spec
endif
zeig=>zeig%next
end do
!Berechnung LAI und ceppot
ceppot_can = 0.
ceppot_sveg = 0.
LAI_can = 0.
LAI_sveg = 0.
DO i=1,anrspec
ns = nrspec(i)
IF (ns .le. nspec_tree) THEN
LAI_can = LAI_can + svar(ns)%act_sum_lai
ELSE
LAI_sveg = LAI_sveg + svar(ns)%act_sum_lai
ENDIF
ENDDO
DO i=1,anrspec
ns = nrspec(i)
IF (ns .le. nspec_tree) THEN
IF(LAI_can .gt. 0.) THEN
ceppot_can = ceppot_can + svar(ns)%act_sum_lai/LAI_can * spar(ns)%ceppot_spec
ELSE
ceppot_can = 0.
ENDIF
ELSE
IF(LAI_sveg .gt. 0.) THEN
ceppot_sveg = ceppot_sveg + svar(ns)%act_sum_lai/LAI_sveg * spar(ns)%ceppot_spec
ELSE
ceppot_sveg= 0.
ENDIF
END IF
ENDDO
if (LAI .gt. 1.) then
cover = cover3
else if (LAI .le. zero) then
cover = 0.1 * cover3
else
cover = LAI * cover3 ! to combine with leave surface; an Blattflaeche koppeln
endif
call wclas(waldtyp) ! forest type
END SUBROUTINE standup
!******************************************************************************
SUBROUTINE senescence
! update of senescence rates
USE data_stand
USE data_species
USE data_simul
IMPLICIT NONE
! senescence rates
zeig=>pt%first
DO
IF(.not.associated(zeig)) exit
if (zeig%coh%species.ne.nspec_tree+2) then ! exclude mistletoe from senescence
select case (flag_dis)
! case (1,2)
! zeig%coh%sfol = spar(zeig%coh%species)%psf * zeig%coh%x_fol + zeig%coh%x_fol_loss
! zeig%coh%sfrt = spar(zeig%coh%species)%psr * zeig%coh%x_frt + zeig%coh%x_frt_loss
case (0,1,2)
zeig%coh%sfol = spar(zeig%coh%species)%psf * zeig%coh%x_fol
zeig%coh%sfrt = spar(zeig%coh%species)%psr * zeig%coh%x_frt
end select
IF (zeig%coh%height.ge.thr_height .and.zeig%coh%species.LE. nspec_tree) THEN
zeig%coh%ssap = spar(zeig%coh%species)%pss * zeig%coh%x_sap
ELSE
zeig%coh%ssap = 0
if(zeig%coh%species.GT.nspec_tree) zeig%coh%ssap = spar(zeig%coh%species)%pss*zeig%coh%x_sap
ENDIF
end if !exclusion of mistletoe
zeig=>zeig%next
END DO
END SUBROUTINE senescence
!**************************************************************
SUBROUTINE litter
! Calculation of summation variables of litter fractions
use data_par
use data_out
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
implicit none
real hconvd
integer taxnr, i
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' litter'
zeig => pt%first
do while (associated(zeig))
taxnr = zeig%coh%species
if(taxnr.le.nspec_tree) then
totfol_lit_tree = totfol_lit_tree + zeig%coh%litC_fol
totfrt_lit_tree = totfrt_lit_tree + zeig%coh%litC_frt
end if
totfol_lit = totfol_lit + zeig%coh%litC_fol
totfrt_lit = totfrt_lit + zeig%coh%litC_frt
tottb_lit = tottb_lit + zeig%coh%litC_tb
totcrt_lit = totcrt_lit + zeig%coh%litC_crt
totstem_lit = totstem_lit + zeig%coh%litC_stem
zeig => zeig%next
enddo ! zeig (cohorts)
! litter biomass: x kg C/tree to kg/ha (n*x*1000g/(kPatchSize m2)/cpart==> kg/ha)
hconvd = (1000*gm2_in_kgha) / (kpatchsize * cpart) !
totfrt_lit = totfrt_lit * hconvd
totfol_lit = totfol_lit * hconvd
tottb_lit = tottb_lit * hconvd
totcrt_lit = totcrt_lit * hconvd
totstem_lit = totstem_lit * hconvd
totfol_lit_tree = totfol_lit_tree * hconvd
totfrt_lit_tree = totfrt_lit_tree * hconvd
do i = 1,nspec_tree
tottb_lit = tottb_lit + dead_wood(i)%C_tb(1)*gm2_in_kgha
totstem_lit = totstem_lit + dead_wood(i)%C_stem(1)*gm2_in_kgha
enddo
END subroutine litter
!**************************************************************
SUBROUTINE calc_ind_rep
USE data_stand
USE data_species
USE data_simul
implicit none
integer :: i
real :: hi
real, dimension(nspecies) :: rindex_spec
rindex1 = 0.
rindex2 = 0.
if(anz_spec.ne.0) then
hi = 1/real(anz_spec)
rindex_spec = 0.
do i = 1, nspecies
if(sumbio.ne.0) then
rindex_spec(i) = svar(i)%sum_bio/sumbio
end if
end do
rindex1 = 0.
rindex2 = 1.
do i = 1, nspecies
if(rindex_spec(i).ne.0) then
rindex1 = rindex1 + abs(hi -rindex_spec(i))
rindex2 = rindex2 * abs(hi -rindex_spec(i))
end if
end do
if(hi.ne.1) then
rindex1 = 1. - rindex1/(2*(1.-hi))
else
rindex1 = 0.
end if
rindex2 = rindex2**anz_spec
end if
END subroutine calc_ind_rep
!**************************************************************
SUBROUTINE overstorey
use data_out
USE data_stand
USE data_species
USE data_simul
implicit none
real,dimension(nspec_tree) :: mindbh, maxdbh, dminage, dmaxage
integer :: i, nrmin, taxnr, agedmin, agedmax
real :: dbhmin, dbhmax
integer :: anzoverst, nrmax
anzoverst = 0
mindbh=0.
do i =1,nspec_tree
call min_dbh(nrmin,dbhmin,agedmin, i)
mindbh(i) = dbhmin
dminage(i) = agedmin
call max_dbh(nrmax, dbhmax, agedmax, i)
maxdbh(i) = dbhmax
dmaxage(i) = agedmax
end do
if (time.eq.0) then
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
taxnr = zeig%coh%species
if(taxnr .le.nspec_tree) then
if(zeig%coh%x_age.lt. (dminage(taxnr) +20) .and. dminage(taxnr).lt. dmaxage(taxnr)) then
zeig%coh%underst =2
end if
end if
zeig=>zeig%next
end do
else
zeig=>pt%first
do
IF(.not.associated(zeig)) exit
taxnr = zeig%coh%species
if(zeig%coh%height.gt. 130..and. zeig%coh%underst.eq.4) then
zeig%coh%underst = 2
end if
zeig=>zeig%next
end do
end if ! time
END SUBROUTINE overstorey
source_code/version_2.3_windows/stand_mort.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - stand_mort *!
!* - int_mort intrinsic mortality rate *!
!* - stress_mort stress mortality rate *!
!* - int_mort_weib *!
!* *!
!* - Calculation of dead trees per cohort and species *!
!* deterministic approach *!
!* - relative mortality rate is determined by intr. mortality *!
!* and stress mortality *!
!* - stress mortality is calculated depending on *!
!* npp, ystress, yhealth *!
!* - intrinsic probability is optionally calculated on *!
!* age of cohort *!
!* - for each tree of the cohort mortality is decided *!
!* by means of the Mortality probability and *!
!* a uniformly distributed variable *!
!* *!
!* 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 *!
!* *!
!*****************************************************************!
! input variables:
! pro cohort NPP
! state variables:
! pro cohort nTreeA nTreeD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE stand_mort
USE data_stand
USE data_species
USE data_simul
USE data_manag
use data_out
use data_par
IMPLICIT NONE
!local variables
INTEGER :: flag_hgrowth
INTEGER :: taxnr
REAL :: intmort
REAL :: strmort
REAL :: totmort
REAL :: totmort_m
REAL :: besmort
REAL :: ntdead
REAL :: ntdead_m
REAL :: nhelp
REAL :: survpfunct
INTEGER :: hage
REAL :: besp1,besp2 ! parameters for besetting mortality
real :: help1, help2
real :: intmorth
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' stand_mort'
if (flag_standup .eq. 0) flag_standup = 1 ! call stand_balance later
ntdead=0.
nhelp=0.
ntdead_m=0
flag_hgrowth=0
sumvsdead = 0.
sumvsdead_m3 = 0.
besmort = 0
strmort = 0.
totmort = 0.
svar%sumvsdead = 0
svar%sumvsdead_m3 = 0.
zeig=>pt%first
DO
IF(.not.associated(zeig)) exit
IF (zeig%coh%height.ge.thr_height.and. zeig%coh%species.le.nspec_tree) then
taxnr=zeig%coh%species
IF(time.eq.1) then
zeig%coh%nta=zeig%coh%nTreeA
ELSE
IF (flag_mg.eq.1) then
IF (thin_year(act_thin_year-1).eq.(time-1)) zeig%coh%nta=zeig%coh%nTreeA
ELSE IF (flag_mg.eq.2) then
if(flag_adapm .eq. 1) zeig%coh%nta=zeig%coh%nTreeA
ENDIF
ENDIF
IF(zeig%coh%notViable) then
print*,time, zeig%coh%notViable
zeig%coh%nTreeD = zeig%coh%ntreeA
zeig%coh%nta = 0.
zeig%coh%ntreeA = 0
goto 1000
ENDIF
! calculation of stress and health indicator based on foliage biomass increment
hage = zeig%coh%x_age
IF (flag_hgrowth==0) THEN
IF(zeig%coh%fol_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = zeig%coh%x_stress - 1
ENDIF
ENDIF
IF (flag_hgrowth==1) THEN
IF(zeig%coh%bio_inc.le.0.0) then
zeig%coh%x_stress = zeig%coh%x_stress + 1
zeig%coh%x_health= 0
ELSE
zeig%coh%x_health = zeig%coh%x_health + 1
IF(zeig%coh%x_stress.eq.1.and. zeig%coh%x_health.gt.0) zeig%coh%nta = zeig%coh%ntreeA
IF(zeig%coh%x_stress .gt.0) zeig%coh%x_stress = MAX(zeig%coh%x_stress - 3,0)
ENDIF
ENDIF
! calculation of relative mortality rates
! intrinsic mortality
! constant
call int_mort(taxnr,intmorth)
! age-depending using Weibull function
call int_mort_weib(taxnr,intmort,hage)
! stress mortality
IF(zeig%coh%x_stress.gt.0) then
IF(flag_hgrowth==0) strmort = weibal*spar(taxnr)%weibla*zeig%coh%x_stress**(weibal-1)
IF(flag_hgrowth==1) strmort = weibal*spar(taxnr)%weibla*(zeig%coh%x_stress/3.)**(weibal-1)
ELSE
strmort = 0.
ENDIF
!mortality depending on gross growth rate foliage
IF(strmort==0.0 .AND. flag_hgrowth==2) THEN
IF(zeig%coh%sfol/zeig%coh%gfol.GT.0.9) THEN
strmort=((zeig%coh%sfol/zeig%coh%gfol-0.9)*10.)**2
ELSE
ENDIF
ENDIF
if(strmort==0. .and. flag_hgrowth==3) then
help1 = 10**((log10(4.5)-log10((zeig%coh%x_sap + zeig%coh%x_hrt))*zeig%coh%ntreea)/1.5)
help2 = help1/zeig%coh%ntreea
end if
! mortality caused by besetting for oak
besp1= 0.018
besp2= 0.0216
if(zeig%coh%species.eq.4) then
if( zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = besp1*zeig%coh%bes- besp2
end if
else if (zeig%coh%species.eq.1.) then
if(zeig%coh%bes.le. 1.2) then
besmort = 0.
else
besmort = 0.04*zeig%coh%bes- 0.04
end if
end if
!total mortality rate depending on flag_mort
IF(flag_mort.eq.1) THEN
totmort = strmort
ELSE IF(zeig%coh%x_age.le.30) then
totmort=intmort+(1-intmort)*strmort
if(taxnr.eq.8) totmort = strmort
ELSE
totmort=intmort+(1-intmort)*strmort
ENDIF
! if species type oak then combination of stress mortality and besetting mortality
if(zeig%coh%species.eq.4.or.zeig%coh%species.eq.1) then
totmort = besmort + (1-besmort)* strmort
end if
survpfunct = exp(- spar(taxnr)%weibla * zeig%coh%x_stress**weibal)
ntdead = totmort*zeig%coh%nTreeA
IF(totmort.GT.1) CALL error_mess(time,"totmort greater 1: ",totmort)
! calculation of real stem number
zeig%coh%nta = zeig%coh%nta - ntdead
! calculation of integer stem number
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) zeig%coh%nTreeA=0.
if (zeig%coh%mistletoe.eq.1) then ! in case Mist.infect. tree dies, number of mistletoes dies, too
totmort_m = zeig%coh%nTreeD/(zeig%coh%nTreeD+zeig%coh%nTreeA) ! share of trees removed of total trees. used for the share of mistletoe that dies
ntdead_m = 1. !flag
end if
! calculation of the litter pool of a died tree of a cohort
1000 IF (zeig%coh%ntreeD.ne.0) then
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreeD*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreeD*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreeD*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreeD*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreeD*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
if(flag_mg.ne.0) then
! if(thin_dead.eq.0.and.thin_flag1(1).lt.0.) then
if(thin_dead.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
end if
else if(zeig%coh%diam.le.tardiam_dstem.or. flag_mg.eq.0) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(zeig%coh%diam.gt.tardiam_dstem.and.flag_mg.ne.0.or.flag_mg.eq.5) then
sumvsdead = sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
svar(taxnr)%sumvsdead= svar(taxnr)%sumvsdead + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)
sumvsdead_m3 = sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt) /(spar(taxnr)%prhos*1000000)
svar(taxnr)%sumvsdead_m3 = svar(taxnr)%sumvsdead_m3 + zeig%coh%ntreeD*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(taxnr)%prhos*1000000)
else if(thin_dead.eq.1.and.zeig%coh%diam.le.tardiam_dstem) then
zeig%coh%litC_stem =zeig%coh%litC_stem + zeig%coh%ntreeD*(zeig%coh%x_sap+zeig%coh%x_hrt)*cpart
zeig%coh%litN_stem =zeig%coh%litC_stem/spar(taxnr)%cnr_stem
end if
ENDIF
ENDIF
zeig=>zeig%next
ENDDO
! if tree cohort with mistletoe changed, change number of mistletoes too
if (ntdead_m.eq.1.) then
zeig => pt%first
do while (associated(zeig))
if (zeig%coh%species.eq.nspec_tree+2) then
zeig%coh%nta=zeig%coh%nTreeA
ntdead_m = totmort_m*zeig%coh%nTreeA
zeig%coh%nta = zeig%coh%nta - ntdead_m
zeig%coh%nTreeD = zeig%coh%nTreeA-NINT(zeig%coh%nta)
zeig%coh%nTreeA = NINT(zeig%coh%nta)
IF(zeig%coh%nTreeA.lt.1.) then
zeig%coh%nTreeA=0.
flag_mistle=0 !set flag mistletoe back to zero
ENDIF
endif
zeig=>zeig%next
enddo ! zeig (cohorts)
end if
ntdead_m=0.
! recalculation sumvsdead
sumvsdead = sumvsdead * 10000./kpatchsize ! kg/patch ---> ! kg/ha
sumvsdead_m3 = sumvsdead_m3 * 10000./kpatchsize ! kg/patch ---> ! kg/ha
cumsumvsdead = cumsumvsdead + sumvsdead
END SUBROUTINE stand_mort
SUBROUTINE int_mort(taxnr,intmort)
USE data_species
IMPLICIT NONE
REAL :: intmort
INTEGER :: taxnr
intmort=1.-exp(-spar(taxnr)%intr)
END SUBROUTINE int_mort
SUBROUTINE int_mort_weib(taxnr,intmort,hage)
USE data_species
USE data_stand
USE data_simul
IMPLICIT NONE
REAL :: intmort, weibla_int
INTEGER :: taxnr
INTEGER :: hage
! Weibull functions depending on age
weibla_int = -log(0.01)/(spar(taxnr)%max_age**weibal_int)
intmort = weibal_int*weibla_int*(hage)**(weibal_int-1.)
END SUBROUTINE int_mort_weib
source_code/version_2.3_windows/stand_regen.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* subroutine for regeneration *!
!* including the SR: *!
!* - gener_seed *!
!* - seed_ini *!
!* - simseed *!
!* - growth_seed *!
!* - mort_seed *!
!* *!
!* 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_regen
USE data_simul
USE data_stand
IMPLICIT NONE
flag_standup = 2
CALL mort_seed
CALL gener_seed
END SUBROUTINE stand_regen
SUBROUTINE gener_seed
USE data_stand
USE data_species
USE data_simul
use data_out
USE data_plant
USE data_soil
IMPLICIT NONE
real :: seedla ! leaf area of all seedling cohorts
real :: laiseed ! lai ----"------
integer :: nseed ! number of generated seeds
real :: redseed
integer :: i
integer, dimension(5) :: agemin, seedpot
real, dimension(5,3) :: latg
real :: help,help1, help2
real :: pequal
integer :: hlayer
integer :: flag_reg_help
TYPE(coh_obj), POINTER :: p
DATA latg /1.,0.3,0.1,1.,0.1,1.,0.9,0.5,1.,0.5,1.,1.,0.9,1.,0.9/
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' gener_seed'
flag_reg_help = 0
agemin = 0
seedpot = 0
seedla = 0.
help1 = 0.
! calculation of leafarea of all seedling cohorts
SELECT CASE (flag_reg)
! according to FORGRA (Vincent Kint)
CASE (30)
call random_number(pequal)
DO i= 1, nspecies
nseed = 0
p =>pt%first
DO WHILE (ASSOCIATED(p))
if(p%coh%species.eq.i) then
if (i.eq.1) then
agemin(i) = 50 + (30* pequal)
call random_number(pequal)
seedpot(i) = 810*pequal
else if(i.eq.3) then
agemin(i) = 15 + (45* pequal)
call random_number(pequal)
seedpot(i) = 1000*pequal
else if(i.eq.4) then
agemin(i) = 15 +(35* pequal)
call random_number(pequal)
seedpot(i) = 1125*pequal
else if(i.eq.5) then
agemin(i) = 10 +(10* pequal)
call random_number(pequal)
seedpot(i) = 8750*pequal
end if
if(p%coh%x_age.ge. agemin(i).and.p%coh%diam.gt.0.) then
nseed = nseed + seedpot(i)*(p%coh%ntreem + p%coh%ntreea)
end if
end if
p => p%next
END DO ! cohort
help2 = irelpool_ll
if(help2.lt.0) help2 =0
if(help2.eq. 0.) then
redseed = 0.
else if( help2.gt.0. .and. help2.le.latg(i,1)) then
redseed = help2*latg(i,1)/0.4
else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
redseed = help2*latg(i,2)/0.6
else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
redseed = help2*latg(i,3)/0.8
else if(help2.gt.latg(i,3)) then
redseed = help2* latg(i,3)
end if
nseed = nseed * redseed
! for birch 1 year old saplings are used
if (i.eq.5) then
numplant(i) = nseed
flag_reg = 14
if(nseed.ne.0) call planting
flag_reg= 0
else
call seed_multi(nseed,i)
end if
END DO ! species
CASE(1,2,3)
p =>pt%first
DO WHILE (ASSOCIATED(p))
if(p%coh%height.lt.thr_height) then
seedla = seedla + p%coh%t_leaf*p%coh%ntreea
help1 = help1 + p%coh%x_fol*p%coh%ntreea
end if
p => p%next
END DO
! calculation LAI of lowest_layer
laiseed=seedla/kpatchsize
DO i=1,nspecies
IF (spar(i)%regflag.eq.1) THEN
CALL simseed(i,nseed)
IF(laiseed.lt.1) THEN
! reduction of seedling number nseed depending on light module and free space in the lowest_layer
SELECT CASE (flag_light)
CASE(1)
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
CASE (2)
if (anz_coh.eq. 0) then
if(time.eq.1) then
hlayer = 0
else
hlayer = 1
end if
else
hlayer = lowest_layer -1
end if
help = vstruct(hlayer)%Irel
if (help.lt.0.05) help = 0
IF(help.eq.0) THEN
nseed = 0
ELSE
nseed = nseed*help
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
END IF
CASE(3)
redseed= bgpool_ll
nseed = nseed*redseed
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
CALL seed_multi(nseed,i)
END IF
CASE(4)
if(i.gt.5) then
redseed= irelpool_ll
else
! according to FORGRA, not for all species (i=1,5)
help2 = irelpool_ll
if(help2.lt. 0.01) then
redseed = 0.
else if( help2.gt.0.01 .and. help2.le.latg(i,1)) then
redseed = help2*latg(i,1)/0.4
else if ( help2.gt.latg(i,1).and. help2.le.latg(i,2)) then
redseed = help2*latg(i,2)/0.6
else if ( help2.gt.latg(i,2).and. help2.le.latg(i,3)) then
redseed = help2*latg(i,3)/0.8
else if(help2.gt.latg(i,3)) then
redseed = help2* latg(i,3)
end if
end if
nseed = redseed * nseed
IF(flag_reg.ne.3) THEN
CALL seed_ini(nseed,i)
ELSE
if (i.eq.5) then
numplant(i) = nseed
flag_reg_help = flag_reg
flag_reg = 14
if(nseed.ne.0) call planting
flag_reg = flag_reg_help
else
CALL seed_multi(nseed,i)
end if
END IF
END SELECT
ELSE
nseed = 0.
END IF
ELSE
nseed=0.
END IF
END DO
END SELECT ! flag_reg
END subroutine gener_seed
SUBROUTINE simseed(specnum,nseed)
USE data_species
use data_simul
use data_stand
IMPLICIT NONE
REAL :: pequal
INTEGER :: nseed,specnum
REAL :: seedmax
! calculation of max. seedrate of patch from max. seedrate per m2
seedmax=spar(specnum)%seedrate*kpatchsize
CALL random_number(pequal)
CALL random_number(pequal)
nseed=-seedmax*alog(1.-pequal)
IF (flag_mg ==4 .and. time.eq.1) THEN
nseed = NINT(spar(specnum)%seedrate*kpatchsize)
ELSE IF(flag_mg ==4.and. time.gt.1)THEN
nseed = 0
END IF
end
SUBROUTINE seed_ini(nseed,nsp)
USE data_species
use data_stand
use data_help
use data_out
use data_simul
use data_soil
IMPLICIT NONE
integer :: nseed, nr, j
integer :: nsp
REAL :: shoot
REAL :: x1,x2,xacc,root
REAL :: rtflsp
REAL :: troot2
TYPE(cohort) ::tree_ini
external weight
external rtflsp
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' seed_ini'
IF(nseed.eq.0) RETURN
hnspec = nsp
max_coh = max_coh + 1
! nullify of all elements
call coh_initial (tree_ini)
tree_ini%ident = max_coh
tree_ini%species = nsp
tree_ini%ntreea = nseed
tree_ini%nta = nseed
tree_ini%x_age = 1
mschelp = spar(nsp)%seedmass/1000. ! g ---> kg
x1 = 0.
x2 = 0.1
xacc = (1.0e-10) * (x1+x2)/2
root = rtflsp(weight,x1,x2,xacc)
tree_ini%x_sap = root
shoot = root*1000. ! [kg]
tree_ini%x_fol= (spar(nsp)%seeda*(tree_ini%x_sap** spar(nsp)%seedb)) ![kg]
tree_ini%x_frt = tree_ini%x_fol ! [kg]
tree_ini%med_sla = spar(nsp)%psla_min + spar(nsp)%psla_a*0.5
tree_ini%t_leaf = tree_ini%med_sla* tree_ini%x_fol ! [m-2]
tree_ini%ca_ini = tree_ini%t_leaf
! initialize pheno state variables
IF(spar(tree_ini%species)%Phmodel==1) THEN
tree_ini%P=0
tree_ini%I=1
ELSE
tree_ini%P=0
tree_ini%I=0
tree_ini%Tcrit=0
END IF
! tranformation of shoot biomass kg --> mg
if(nsp.ne.2)tree_ini%height = spar(nsp)%pheight1*(shoot*1000.)**spar(nsp)%pheight2 ! [cm] calculated from shoot biomass (mg); berechnet aus shoot biomass (mg)
if(nsp.eq.2) tree_ini%height = 10**(spar(nsp)%pheight1+ spar(nsp)%pheight2*LOG10(shoot*1000.)+ &
spar(nsp)%pheight3*(LOG10(shoot*1000.))**2)
IF(nseed.ne.0.) then
IF (.not. associated(pt%first)) THEN
ALLOCATE (pt%first)
pt%first%coh = tree_ini
NULLIFY(pt%first%next)
! root distribution
call root_depth (1, pt%first%coh%species, pt%first%coh%x_age, pt%first%coh%height, pt%first%coh%x_frt, pt%first%coh%x_crt, nr, troot2, pt%first%coh%x_rdpt, pt%first%coh%nroot)
pt%first%coh%nroot = nr
do j=1,nr
pt%first%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
pt%first%coh%rooteff = 0. ! layers with no roots
enddo
ELSE
ALLOCATE(zeig)
zeig%coh = tree_ini
zeig%next => pt%first
pt%first => zeig
! root distribution
call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
zeig%coh%nroot = nr
do j=1,nr
zeig%coh%rooteff = 1. ! assumption for the first use
enddo
do j=nr+1, nlay
zeig%coh%rooteff = 0. ! layers with no roots
enddo
END IF
anz_coh=anz_coh+1
END IF
END SUBROUTINE seed_ini
SUBROUTINE growth_seed
USE data_stand
USE data_species
USE data_simul
USE data_par
use data_out
IMPLICIT NONE
REAL :: lambdaf = 0., & ! partitioning functions
lambdas = 0., &
lambdar = 0., &
NPP = 0., & ! annual NPP
F = 0., & ! state variables: foliage,
S = 0., & ! shoot biomass,
R = 0., & ! fine roots,
H = 0., & ! total tree height
FNew, SNew, & ! new states
RNew, &
sigmaf = 0., & ! current leaf activity rate
sigman = 0., & ! current root activity rate
betar = 0., &
ar = 0
REAL :: Sf, & ! senescence rates
Sr, &
Gf, & ! growth rates
Gs, &
Gr
real :: pab,helpdr,helpsum
TYPE(coh_obj), POINTER :: p
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' growth_seed'
flag_standup = 2 ! call stand_balance and root_distribution later
p=>pt%first
DO
if(.not.associated(p)) exit
if( p%coh%height.lt.thr_height.and. p%coh%species.le.nspec_tree) then
ns = p%coh%species
F = p%coh%x_fol
S = p%coh%x_sap
R = p%coh%x_frt
NPP = p%coh%NPP
IF (flag_reg .eq. 2) NPP = p%coh%NPPpool ! [kg]
H = p%coh%height
Sf = p%coh%sfol
Sr = p%coh%sfrt
! only allocate if enough NPP is available
1 IF (NPP>1.0E-9.or. NPP.ge.(Sf+Sr).and.(sr+Sf)>1.0E-9) THEN
! calculate leaf activity based on net PS and leaf mass
sigmaf = NPP/F
! calculate root activity based on drought index
helpdr= p%coh%drIndAl / p%coh%nDaysGr
IF (flag_sign.eq.1) THEN
sigman = amax1(spar(ns)%sigman*10*(((5.-spar(ns)%stol)*1.-p%coh%crown_area)/(5.-spar(ns)%stol)*1.),spar(ns)%sigman) * p%coh%drIndAl / p%coh%nDaysGr
ELSE
sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
END IF
! auxiliary variables for fine roots
if(helpdr.lt.0.001) ar = 0.
ar = spar(ns)%pcnr * sigmaf / sigman
betar = (Sr - R + ar*(F-Sf)) / NPP
! calculate coefficients for roots and foliage and shoot
select case (ns)
case (1)
pab = 0.487
case(2)
pab = 0.826
case(3)
pab=1.9
case(4)
pab=1.002
! Pinus contorta
case(6)
! Gholz
pab = 0.236
! Populus tremula
case(8)
pab = 0.3233
case(9)
! Pinus halepensis
pab = 1.42335
case(10)
! pseudotsuga menziesii
pab = 0.8515
case(11)
! Robinia
pab = 0.8594
end select
lambdaf = (pab*(1-betar)+ (Sf/NPP))/(1 + pab*(1. + ar))
lambdar = ar * lambdaf + betar
lambdas = 1.- lambdaf - lambdar
! consistency
ELSE
lambdaf = 0.
lambdas = 0.
lambdar = 0.
END IF
if ( lambdas.lt.0.) then
lambdas = 0.
lambdaf = (1.-betar)/(ar+1)
lambdar = 1.-lambdaf
if( lambdar.lt.0) then
lambdar=0
lambdaf=1
end if
if(lambdaf.lt.0) then
lambdaf =1
lambdar = 0.
end if
endif
helpsum = lambdaf + lambdar+ lambdas
Gf = lambdaf*NPP
Gr = lambdar*NPP
Gs = lambdas*NPP
p%coh%gfol = Gf
p%coh%gfrt = Gr
p%coh%gsap = Gs
! update of state vector
FNew = F + Gf - Sf
SNew = S + Gs
RNew = R + Gr - Sr
p%coh%x_fol = FNew
p%coh%x_sap = SNew
p%coh%x_frt = RNew
p%coh%fol_inc_old = p%coh%fol_inc
p%coh%fol_inc = Gf - Sf
p%coh%stem_inc = Gs
! update height and shoot base diameter (regression functions from Schall 1998)
IF(ns.ne.2) p%coh%height = spar(ns)%pheight1* (snew*1000000.) **spar(ns)%pheight2
IF(ns.eq.2) p%coh%height = 10**(spar(ns)%pheight1+ spar(ns)%pheight2*LOG10(snew*1000000.)+ &
spar(ns)%pheight3*(LOG10(snew*1000000.))**2)
p%coh%height_ini = p%coh%height
! update foliage area, parameter med_sla
SELECT CASE (flag_light)
CASE (1:2)
p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.- vstruct(lowest_layer)%irel)
CASE(3,4)
p%coh%med_sla = spar(ns)%psla_min + spar(ns)%psla_a*(1.-irelpool(lowest_layer))
END SELECT
! total leaf area of a tree in this cohort [m**2]
p%coh%ca_ini = p%coh%med_sla * p%coh%x_fol
! update age -now not necessary this is done in stand_bal
p%coh%notViable = (FNew <= 0.) .OR. (SNew <= 0.) .OR. &
(RNew <= 0.)
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreea * Sf * cpart
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreea * Sr * cpart
! with species specific N content and reallocation factor (see species.par)
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * spar(ns)%reallo_fol / spar(ns)%cnr_fol
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreea * Sr * cpart * spar(ns)%reallo_frt / spar(ns)%cnr_frt
end if ! seedling cohort test
p=> p%next
END DO
END SUBROUTINE growth_seed
SUBROUTINE mort_seed
USE data_species
USE data_simul
use data_stand
use data_par
use data_out
IMPLICIT NONE
integer :: taxnr
integer :: hage
real :: intmort
real :: strmort
real :: totmort
real :: ntdead
real :: ntahelp
TYPE(coh_obj), POINTER :: p
if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' mort_seed'
p=>pt%first
DO
IF(.not.associated(p)) EXIT
IF(p%coh%height.lt.thr_height) THEN
IF(p%coh%notViable) then
PRINT*,time, p%coh%notViable
p%coh%ntreed = p%coh%ntreea
p%coh%ntreea = 0
ENDIF
END IF
p => p%next
END DO
p => pt%first
DO
IF(.not.associated(p)) EXIT
IF(p%coh%height.lt.thr_height .and. p%coh%species.le.nspec_tree) THEN
taxnr = p%coh%species
if(p%coh%ntreea .eq.0) goto 1000
hage = p%coh%x_age
IF( p%coh%fol_inc.le.0.) THEN
p%coh%x_stress = p%coh%x_stress +1
p%coh%x_health = 0
ELSE
p%coh%x_health = p%coh%x_health +1
IF(p%coh%x_stress .gt. 0.) p%coh%x_stress = p%coh%x_stress -1
ENDIF
! intrinsic mortality
CALL int_mort_weib(taxnr, intmort, hage)
! stress mortality
IF(p%coh%x_stress.gt.0) THEN
strmort = weibal*spar(taxnr)%weibla*p%coh%x_stress**(weibal-1)
ELSE
strmort = 0.
ENDIF
totmort=intmort+(1-intmort)*strmort
! calculation of real number of dying stems
ntdead = totmort*p%coh%ntreeA
! update of real stem number nta and number of dead stems
p%coh%nta = p%coh%nta -ntdead
p%coh%nTreeD = p%coh%nTreeA-NINT(p%coh%nta)
! help variable for comparison
ntahelp = p%coh%nTreeA
! update of integer stem number
p%coh%nTreeA = NINT(p%coh%nta)
! update of integer stem number
if(p%coh%nta.lt.1.) p%coh%nTreeA=0.
! update of real stem number if integer stem number was changed
if (ntahelp .ne. p%coh%nTreeA ) p%coh%nta = p%coh%nTreeA
1000 if (p%coh%ntreeD.ne.0) then
p%coh%litC_fol = p%coh%litC_fol + p%coh%ntreeD*(1.-spar(taxnr)%psf)*p%coh%x_fol*cpart
p%coh%litN_fol = p%coh%litN_fol + p%coh%ntreeD*((1.-spar(taxnr)%psf)*p%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
p%coh%litC_frt = p%coh%litC_frt + p%coh%ntreeD*p%coh%x_frt*cpart
p%coh%litN_frt = p%coh%litN_frt + p%coh%ntreeD*p%coh%x_frt*cpart/spar(taxnr)%cnr_frt
p%coh%litC_tb = p%coh%litC_tb + p%coh%ntreeD*p%coh%x_tb*cpart
p%coh%litN_tb = p%coh%litN_tb + p%coh%ntreeD*p%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
p%coh%litC_crt = p%coh%litC_crt + p%coh%ntreeD*p%coh%x_crt*cpart
p%coh%litN_crt = p%coh%litN_crt + p%coh%ntreeD*p%coh%x_crt*cpart/spar(taxnr)%cnr_tbc
p%coh%litC_stem = p%coh%litC_stem + p%coh%ntreeD*(p%coh%x_sap)*cpart
p%coh%litN_stem = p%coh%litC_stem/spar(taxnr)%cnr_stem
endif
END IF
p => p%next
ENDDO
END SUBROUTINE mort_seed
source_code/version_2.3_windows/statistik.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* *!
!* statistical analysis of model quality *!
!* *!
!* Author: F. Suckow *!
!* *!
!* contains: *!
!* residuen *!
!* statistik *!
!* mean (n, arr) *!
!* variance (n, meanv, arr) *!
!* correl (n, meanv1, arr1, meanv2, arr2) *!
!* sumsq (n, arr) *!
!* stat_mon *!
!* *!
!* 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 residuen (ip)
use data_mess
implicit none
integer i,j, ires, ip
! calculate and save residues, with date, simulation as well as measurement value
! Residuen berechnen, mit Datum, Sim.- und Messwert speichern
if (ip .eq. 1) then
allocate (val(imkind))
do i=1,imkind
ires = 0
val(i)%tkind = tkind
allocate (val(i)%day(1:anz_val))
allocate (val(i)%year(1:anz_val))
allocate (val(i)%resid(1:anz_val))
allocate (val(i)%sim(1:anz_val))
allocate (val(i)%mess(1:anz_val))
val(i)%day = -99
val(i)%year = -99
val(i)%resid = -9999.0
val(i)%mkind = sim_kind(i)
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%day(ires) = stz(1,j)
val(i)%year(ires) = stz(2,j)
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
val(i)%mess(ires) = mess2(j,i)
else
endif
enddo
val(i)%imes = ires
enddo
else
do i=1,imkind
ires = 0
val(i)%resid = -9999.9
val(i)%sim = -9999.9
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
else
endif
enddo
enddo
endif
END SUBROUTINE residuen
!**************************************************************
SUBROUTINE statistik
use data_mess
use data_simul
implicit none
integer imt ! aktueller Messwert-Typ
real, external :: mean, variance, correl, sumsq
integer i, n, nhelp
real help, h1, h2
real, allocatable, dimension(:):: arr, arrs, arrm, harr
real:: avs, & ! mean value simulation; Mittelwert Simulation
mins, & ! Minimum Simulation
maxs, & ! Maximum Simulation
stdevs, & ! standard deviation simulation; Standardabweichung Simulation
varis, & ! scattering of simulation; Streuung Simulation
varcos, & ! coefficient of variation for simulation; Variationskoeffizient Simulation
avm, & ! mean value measurements; Mittelwert Messwerte
minm, & ! minimum value measurements; Minimum Messwerte
maxm, & ! maximum value measurements; Maximum Messwerte
stdevm, & ! standard deviation measurements; Standardabweichung Messwerte
varim, & ! scattering of measurements; Streuung Messwerte
varcom, & ! coefficient of variation of measurements; Variationskoeffizient Messwerte
corrco, & ! coefficient of correlation; Korrelationskoeffizient
rsq, & ! coefficient of determination; Bestimmtheitsmass
avr, & ! mean error residues; Mittlerer Fehler Residuen
minr, & ! minimum residues; Minimum Residuen
maxr, & ! maximum residues; Maximum Residuen
stdevr, & ! standard deviation residues; Standardabweichung Residuen
varir, & ! scattering of residues; Streuung Residuen
varcor, & ! coefficient of variation residues; Variationskoeffizient Residuen
nme, & ! normalised mean error; Normalisierter mittlerer Fehler
mae, & ! mean absolute error of residues; Mittlerer absoluter Fehler Residuen
nmae, & ! normalised mean absolute error; Normalisierter mittlerer absoluter Fehler
sse , & ! sum of squared errors; Fehlerquadratsumme
rmse, & ! Root mean square error
nrmse, & ! Normalised root mean square error
pme, & ! mean procental error; Mittlerer prozentualer Fehler
prmse, & ! mean squared procentual error; Mittlerer quadratischer prozentualer Fehler
tic, & ! Theilsch imbalance coefficient; Theilscher Ungleichheitskoeffizient
meff ! Model efficiency (Medlyn et al. 2005)
do imt = 1, imkind
n = val(imt)%imes
if (n .gt. 0) then
allocate (arr(n))
allocate (arrs(n))
allocate (arrm(n))
allocate (harr(n))
! simulation
arrs = val(imt)%sim(1:n)
avs = mean(n, arrs)
mins = minval(arrs)
maxs = maxval(arrs)
varis = variance(n, avs, arrs)
if (varis .ge. 0.) then
stdevs = sqrt(varis)
else
stdevs = 0.
endif
if (avs .ne. 0.) then
varcos = stdevs / avs
else
varcos = -9999.0
endif
! observed
arrm = val(imt)%mess(1:n)
avm = mean(n, arrm)
minm = minval(arrm)
maxm = maxval(arrm)
varim = variance(n, avm, arrm)
if (varim .ge. 0.) then
stdevm = sqrt(varim)
else
stdevm = 0.
endif
! residuals
arr = val(imt)%resid(1:n)
avr = mean(n, arr)
minr = minval(arr)
maxr = maxval(arr)
varir = variance(n, avr, arr)
if (varir .ge. 0.) then
stdevr = sqrt(varir)
else
stdevr = 0.
endif
if (avr .ne. 0.) then
varcor = stdevr / avr
else
varcor = -9999.0
endif
corrco = correl(n, avs, arrs, avm, arrm)
if (corrco .ge. -1.) then
rsq = corrco * corrco
rsq_av = rsq_av + rsq
else
imk_rsq = imk_rsq - 1
rsq = -9999.0
endif
if (avs .ne. 0.) then
nme = (avm - avs) / avs
nme_av = nme_av + nme
else
imk_nme = imk_nme - 1
nme = -9999.0
endif
mae = mean(n, abs(arr))
sse = sumsq(n, arr)
rmse = sqrt(sse / n)
if (avm .ne. 0.) then
varcom = stdevm / avm
nrmse = rmse / abs(avm)
nrmse_av = nrmse_av + nrmse
nmae = mae / abs(avm)
nmae_av = nmae_av + nmae
else
imk_nrmse = imk_nrmse - 1
imk_nmae = imk_nmae - 1
varcom = -9999.0
nrmse = -9999.0
nmae = -9999.0
endif
nhelp = n
do i = 1, n
if (arrm(i) .ne. 0.) then
harr(i) = abs(arr(i)/arrm(i))
else
nhelp = nhelp -1
harr(i) = 0
endif
enddo
pme = mean(nhelp, harr)
prmse = sumsq(nhelp, harr)
prmse = sqrt(prmse / nhelp)
tic = sse / (sumsq(n, arrs) + sumsq(n, arrm))
h1=sumsq(n, arr)
harr = arrm-avm
h2=sumsq(n, harr)
if (h2.ne.0) then
meff = 1. - (h1 / h2)
else
meff=1
end if
! Mittelwert
pme_av = pme_av + pme
prmse_av = prmse_av + prmse
tic_av = tic_av + tic
meff_av = meff_av + meff
deallocate (arr)
deallocate (arrm)
deallocate (arrs)
deallocate (harr)
write (unit_stat, '(I5,2X, A20,1X,A10,I8,1X,30E13.5)') ip, site_name(ip), val(imt)%mkind, val(imt)%imes, &
avr, minr, maxr, stdevr, varir, varcor, nme, mae, nmae, sse, rmse, nrmse, pme, prmse, tic,meff, corrco, rsq, &
avs, mins, maxs, stdevs, varis, varcos, avm, minm, maxm, stdevm, varim, varcom
endif
enddo
END SUBROUTINE statistik
!**************************************************************
REAL FUNCTION mean (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i)
enddo
mean = help / n
END FUNCTION mean
!**************************************************************
REAL FUNCTION variance (n, meanv, arr)
integer n, i
real, dimension(n):: arr
real meanv, help, xx
help = 0.
if (n .gt. 1) then
do i = 1,n
xx = arr(i) - meanv
help = help + xx * xx
enddo
variance = help / (n -1)
else
variance = -9999.0
endif
END FUNCTION variance
!**************************************************************
REAL FUNCTION correl (n, meanv1, arr1, meanv2, arr2)
integer n, i
real, dimension(n):: arr1, arr2
real meanv1, meanv2, help1, help2, help3, xx1, xx2
help1 = 0.
help2 = 0.
help3 = 0.
do i = 1,n
xx1 = arr1(i) - meanv1
xx2 = arr2(i) - meanv2
help1 = help1 + xx1 * xx2
help2 = help2 + xx1 * xx1
help3 = help3 + xx2 * xx2
enddo
if ((help2 .gt. 1.E-06) .and. (help3 .gt. 1.E-06)) then
correl = help1 / sqrt(help2*help3)
else
correl = -9999.0
endif
END FUNCTION correl
!**************************************************************
REAL FUNCTION sumsq (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i) * arr(i)
enddo
sumsq = help
END FUNCTION sumsq
!**************************************************************
Subroutine stat_mon
! Statistics of monthly values, derived from daily observed values
use data_mess
use data_out
use data_simul
implicit none
integer i, j, k, l
integer dd, mm, yy, doy, yanz, arranz
integer :: outunit ! output unit
character(250) text, filename
character(20) idtext, datei, vunit, obskind
character(150) htext
real, allocatable, dimension(:):: helparr ! help array with montly values of one month for all years
real, allocatable, dimension(:,:):: help_mon ! array with monthly values for each year, year
real, allocatable, dimension(:,:):: help_day ! array with mean daily values per month for each year, year
integer, allocatable, dimension(:,:):: help_num ! array with number of measurement values for each year, year
yanz = mtz(2,imess) - mtz(2,1) + 1
if (.not. allocated(unit_mon)) then
allocate(unit_mon(imkind))
allocate(unit_mon_stat(imkind))
allocate(helparr(yanz))
endif
if (.not. allocated(help_mon)) then
allocate(help_mon(12,yanz))
allocate(help_day(12,yanz))
allocate(help_num(12,yanz))
endif
do k = 1, imkind
help_mon = 0.0
help_num = 0
obskind = sim_kind(k)
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs'//'.out'
unit_mon(k) = getunit()
open(unit_mon(k),file=filename,status='replace')
! Calculate mmonthly sums
do j = 1, imess
doy = mtz(1,j)
yy = mtz(2,j)
call TZINDA(dd,mm,yy,doy)
yy = mtz(2,j) - mtz(2,1) + 1
if (mess1(j,k) .Gt. -9990.) then
if (sim_kind(k) .eq. 'AET') then
if (mess1(j,k) .lt. 0.) then
mess_info = '# negative AET set to zero'
mess1(j,k) = 0. ! avoid negative AET
endif
endif
help_mon(mm,yy) = help_mon(mm,yy) + mess1(j,k)
help_num(mm,yy) = help_num(mm,yy) + 1
endif
enddo ! j
do j = 1, yanz
do i = 1,12
if (help_num(i,j) .gt. 0) then
help_day(i,j) = help_mon(i,j) / help_num(i,j)
else
help_mon(i,j) = -9999.
help_day(i,j) = -9999.
endif
enddo
enddo
! Output of monthly sums
select case (trim(obskind))
case ('AET')
vunit = 'mm'
case ('GPP', 'NPP', 'TER')
vunit = 'g C/m'
case ('Snow')
return
case default
vunit = ' '
end select
write (unit_mon(k), '(A)') '# Monthly sum, daily mean of month and number of values per month of observed '//trim(obskind)
write (unit_mon(k), '(A)') '# '//trim(vunit)
write (unit_mon(k), '(A)', advance='no') '# Year'
do i = 1,12
write (unit_mon(k), '(A8,I2)', advance='no') trim(obskind)//'_',i
enddo
write (unit_mon(k), '(A)')
l = 0
do j = mtz(2,1), mtz(2,imess)
l = l + 1
write (unit_mon(k), '(A,I6,12F10.2)') 'sum ', j, (help_mon(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12F10.2)') 'daily mean ', j, (help_day(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12I10)') 'number ', j, (help_num(i,l), i=1,12)
enddo
! Statistics
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs_stat'//'.res'
outunit = getunit()
open(outunit,file=filename,status='replace')
write (outunit, '(A)') '# Statistics over all years for each monthly sum and daily mean per month of '//trim(obskind)
write (outunit, '(A)') '# '//trim(vunit)
write (outunit, '(A, I6)') '# Simulation period (years): ', year
write (outunit, '(A)') '# site_id Month number Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
write (outunit, '(A)') 'monthly sum'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_mon(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_mon(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
write (outunit, '(A)') ' '
write (outunit, '(A)') 'daily mean per month'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_day(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_day(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
enddo ! k
continue
End Subroutine stat_mon
\ No newline at end of file
source_code/version_2.3_windows/target_thin.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of biomass per *!
!* thinning year as target values *!
!* targetm i given in kg DW/ha *!
!* *!
!* 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 target_thinning(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: targetm ! target value of stem biomass
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpa_us , &
wpb=0, & ! -"-
wpb_us, &
wpc=0, & ! -"-
d63=0, &
d63_us, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
rtarget_help=0, &
target_help1=0,&
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us = 0. , &
stembiom_re = 0. , &
stembiom_all = 0. , &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
! auxilary for thinning routine 4: selective thinning
integer :: count,i, &
idum , third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for population and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
ENDIF
zeig => zeig%next
ENDDO
! mean diamteer for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
! setting of aux. variable target_help
rtarget_help = stembiom_all
! tending
if(thin_tysp(i).eq.4.or.(stembiom_re.ne.0. .and. stembiom_all.eq.0)) then
rtarget_help = stembiom_re
end if
! Umrechnung in Biomasse pro patch in kg
targetm = target_mass(i)*1000*kpatchsize/10000.
! target value of biomass
if(thin_tysp(i).eq.4 .or.(stembiom_re.ne.0. .and. stembiom_all.eq.0) ) then
! tending
targetm = stembiom_re - targetm*stembiom_re
else
end if
if( targetm.eq.1) targetm = 0.
! targetm (kg DW/patch)
! cuttting
if (targetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig=>zeig%next
end do
diff = targetm - rtarget_help
if(diff.lt.0.01) exit
end do
else IF ( targetm .ne. 0.) then
if(target_mass(i).lt.1.) then
targetm = target_mass(i) * rtarget_help
end if
! different thinnings from below and above
select case(thin_tysp(i))
case(1)
! moderate lower thinning;
d_est = 1.02
w_kb = 1.8
case(2)
! intensive lower thinning;
d_est = 1.03
w_kb = 1.5
case(3)
! high thinning;
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin,dbhmin_us, taxnr)
bas_help = bas_area
wpa = dbhmin
wpa_us = dbhmin_us
d63 = mdiam*d_est
d63_us = mdiam_us * d_est
wpb = (d63 - wpa)/ w_kb
wpb_us = (d63_us-wpa_us)/w_kb
wpc = 2
if (thin_tysp(i).eq. 3) then
! starting with overstorey!, continuing with the understorey
if(targetm.lt.(stembiom_all-stembiom)) then
! total removing of overstorey
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! defining the remaining part of stem biomass which has to be remove
rtarget_help = stembiom_us
! understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
! start understorey thinning
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt.10000) exit
end if
zeig=> zeig%next
END DO ! list of potential thinned tees cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
count = count +1
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(count.gt.100000) exit
end do ! understorey thinning
end if ! mediam_us.ne.0
else ! targetm.lt.(stembiom_all-stembiom)
!selection of trees for thinning
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.2*dbh_h
db_u = dbh_h + 0.2*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
! write(1234,*) time, ipot, pequal, isel
! if(isel.eq.0) isel =1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
end do ! total thinning
end if !targetm.lt.(stembiom_all-stembiom)
end if ! thintype 3
if(thin_tysp(i).eq.1.or.thin_tysp(i).eq.2) then
if(targetm.lt.(stembiom_all-stembiom_us)) then
! total removing of understorey
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.2) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
! definging the remaining part of stem biomass which has to remove
rtarget_help = stembiom
if(mdiam.ne.0) then
! additional thinning from the overstorey
counth = 0
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.0) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.3*dbh_h
db_u = dbh_h + 0.3*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1) exit
if(counth.gt.100000) exit
end do ! overstorey thinning
end if ! mdiam.ne.0
else ! targtem.lt.(stembiom_all-stembiom_us)
! first thinning from the understorey
!selection of trees for thinning
if(mdiam_us.ne.0) then
do
call random_number(pequal)
tdbh = wpa_us + wpb_us*(-log(1.-pequal))**(1./wpc)
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth .gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
rtarget_help = rtarget_help - (zeig%coh%x_sap+zeig%coh%x_hrt)
exit
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = rtarget_help - targetm
if(diff.le.0.1 .or. (stembiom_all-stembiom_us).eq.rtarget_help) exit
if(coun1.gt.100000) exit
end do ! understorey thinning
end if ! mdiam_us
end if
end if !! thin_tysp.eq.1 or. thin-tysp.eq.2
END IF ! all thinnings and tending
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
END SUBROUTINE target_thinning
source_code/version_2.3_windows/target_thin_bas.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of stem number per *!
!* thinning year as target values *!
!* rtargetm i given inN/ha *!
! or target value 0 < rtargetm <1 for *!
!* basal area reduction *!
!* *!
!* 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 target_thinning_bas(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: rtargetm=0. ! target value of stem biomass
integer :: target_help = 0.
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpb=0, & ! -"-
wpc=0, & ! -"-
d63=0, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
bas_area_us = 0., &
bas_area_test=0., &
target_help1=0, &
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us, &
stembiom_re, &
stembiom_all, &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
nrmin_us, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
integer :: flag_bas
! auxilarity for thinning routine 4: selective thinning
integer :: count, i, &
idum,third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2 , tar_h
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
flag_bas=0
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
bas_area_us = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
tar_h = 300.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for populations and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
bas_area_us = bas_area_us + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
ENDIF
zeig => zeig%next
ENDDO
! mean diameter for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
anz_tree = anz_tree_dbh + undertree
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
! new basal area management
if(target_mass(i).le.1) then
flag_bas=1
if(thin_stor(i).eq.0) then
target_help = bas_area
else
target_help = bas_area_us
end if
else
if(thin_stor(i).eq.0) then
target_help = anz_tree_dbh
else
target_help = undertree
end if
end if
! tending
if(thin_tysp(i).eq.4) then
target_help = bas_area_us
end if
if(target_mass(i).gt.1) then
rtargetm = target_mass(i)*kpatchsize/10000.
else
rtargetm = target_mass(i)
end if
! cuttting
if (rtargetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
rtargetm = target_mass(i) * target_help
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1 .or. zeig%coh%underst.eq.2) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1 .or.zeig%coh%underst.eq.2 .and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
if (flag_bas.eq.1) then
target_help=target_help - zeig%coh%diam*zeig%coh%diam*pi/4
else
target_help = target_help - zeig%coh%ntreea
end if
exit
end if
zeig=>zeig%next
end do
diff = target_help - rtargetm
if(diff.lt.0.01) exit
end do
! different thinnings from below and above
else IF ( rtargetm .ne. 0.) then
if(target_mass(i).lt.1.) then
rtargetm = target_mass(i) * target_help
! No management if rtargetm=1
else if (rtargetm.eq.1) then
return
endif
select case(thin_tysp(i))
case(1)
! medium lower thinning
d_est = 1.02
w_kb = 2.5
case(2)
! strong lower thinning
d_est = 1.03
w_kb = 1.5
case(3)
! High thinning
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin_us,dbhmin_us, taxnr)
! write (7777,*) time, nrmin,nrmin_us, dbhmin, dbhmin_us
bas_help = bas_area
if (thin_stor(i).eq.0) then
wpa = dbhmin
else
wpa = dbhmin_us
end if
if (thin_stor(i).eq.0) then
d63 = mdiam*d_est
else
d63 = mdiam_us * d_est
end if
wpb = (d63 - wpa)/ w_kb
wpc = 2
wpc = 0.8
if ((thin_tysp(i).ne.4) .and. rtargetm.lt.target_help) then
!selection of trees for thinning
DO ! 11 begin selecting
IF (thin_stor(i) .eq. 2 .and. nrmin_us .eq. -1) EXIT ! exit thinning loop if there are no trees with dbh in understory an understory shall be thinned
coun1=coun1+1
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO ! 2 list preparing
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%notviable.eqv. .TRUE.) then
if(flag_mort.eq.0) then
id_pot(ipot)=zeig%coh%ident
ipot=ipot + 1
endif
else if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! 2 list of potential thinned tree cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO ! 3 removing one tree
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
if(zeig%coh%notviable.eqv..TRUE.) then
if(flag_mort.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea=0
zeig%coh%nta=0
if (flag_bas.eq.1) then
target_help = target_help - (zeig%coh%diam**2)*pi/4
else
target_help = target_help -1
endif
end if
else
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
if(flag_bas.eq.1) then
! write (8888,*) 'target_help', time, target_help, rtargetm
target_help = target_help - (zeig%coh%diam**2)*pi/4
else
target_help = target_help -1
end if
exit
! end if
end if
end if
zeig =>zeig%next
ENDDO ! 3 thinning of one tree
end if
diff = target_help - rtargetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
ENDDO ! 11 total thinning
end if ! thintype 1,2,3
END IF ! all thinnings and tending
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreea>0.and.zeig%coh%species.eq.taxnr) then
bas_area_test = bas_area_test + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
end if
zeig=>zeig%next
end do
! adding biomasses to litter pools depending on stage of stand
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
call class_man
END SUBROUTINE target_thinning_bas
source_code/version_2.3_windows/target_thinstemnum.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutine *!
!* target thinning - *!
!* thinning routine with given values of biomass per *!
!* thinning year as target values *!
!* rtargetm i given in kg DW/ha *!
!* *!
!* 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 target_thinning_OC(i)
use data_stand
use data_manag
use data_simul
use data_species
use data_par
implicit none
real :: rtargetm=0. ! target value of stem biomass
integer :: target_help = 0.
real :: dbhmin=0, &
dbhmin_us = 0, &
wpa=0, & ! Weibull parameter
wpb=0, & ! -"-
wpc=0, & ! -"-
d63=0, &
help=0, &
pequal, &
tdbh=0, &
bas_area=0, &
bas_help=0., &
target_help1=0,&
dbh_h =0, &
db_l = 0., &
db_u = 0., &
d_est=0., &
w_kb=0., &
stembiom, &
stembiom_us, &
stembiom_re, &
stembiom_all, &
diff, &
mdiam, &
mdiam_us
integer :: nrmin, &
nrmin_us, &
lowtree, &
undertree, &
flagth, &
taxnr, &
counth, &
min_id, &
max_id, &
ih1,ih2,ncoh, &
coun1
! auxilarity for thinning routine 4: selective thinning
integer :: count, i, &
idum,third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot
real :: h1, h2 , tar_h
real,external :: gasdev
real:: ran0
! reacalculation of target to kg DW/patch
h1 = 0.
h2 = 0.
count = 0
cohl = -1
flagth = 0
coun1 = 0
help=0.
lowtree=0
undertree = 0
anz_tree_dbh = 0
bas_area = 0.
! stem biomass of overstorey
stembiom = 0.
! stem biomass of understorey
stembiom_us = 0.
stembiom_all = 0.
tar_h = 300.
if (time.eq.73.and. ip .eq.87) then
stembiom = 0
end if
taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Modification for V Kint: no test for diameter
IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
bas_area = bas_area + zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
if( zeig%coh%diam>0) then
anz_tree_dbh = anz_tree_dbh + zeig%coh%ntreea
mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
end if
! Trees with DBH=0 for populations and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
lowtree = lowtree + zeig%coh%ntreea
ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
undertree = undertree + zeig%coh%ntreea
ENDIF
zeig => zeig%next
ENDDO
! mean diameter for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)
third = nint(anz_tree_dbh*0.333333)
anz_tree_ha = nint(anz_tree_dbh*10000./kpatchsize)
anz_tree = anz_tree_dbh + undertree
IF(anz_tree>0)THEN
if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))
ENDIF
if(thin_stor(i).eq.0) then
target_help = anz_tree_dbh
else
target_help = undertree
end if
! tending
if(thin_tysp(i).eq.4) target_help = stembiom_re
if(target_mass(i).gt.1) then
rtargetm = target_mass(i)*kpatchsize/10000.
else
rtargetm = target_mass(i)
end if
! target value of biomass
if(thin_tysp(i).eq.4) then
rtargetm = stembiom_re - rtargetm*stembiom_re
else
end if
! cuttting
if (rtargetm.eq.0.)then
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea = 0
zeig%coh%nta = 0
end if
zeig=> zeig%next
END DO
!tending of regeneration
else if(thin_tysp(i).eq.4) then
min_id = 1000
max_id = 0.
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
ih1 = zeig%coh%ident
if(ih1.lt.min_id) min_id = ih1
ih2 = zeig%coh%ident
if (ih2.gt.max_id) max_id = ih2
end if
zeig=> zeig%next
end do
target_help1 = 0.
do
call random_number(pequal)
ncoh = min_id +(max_id-min_id)*pequal
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
zeig%coh%ntreea = zeig%coh%ntreea - 1
zeig%coh%nta = zeig%coh%nta-1
zeig%coh%ntreem = zeig%coh%ntreem +1
target_help = target_help - zeig%coh%ntreea
exit
end if
zeig=>zeig%next
end do
diff = rtargetm -target_help
if(diff.lt.0.01) exit
end do
! different thinnings from below and above
else IF ( rtargetm .ne. 0.) then
if(target_mass(i).lt.1.) then
rtargetm = target_mass(i) * target_help
! No management if rtargetm=1
else if (rtargetm.eq.1) then
return
endif
select case(thin_tysp(i))
case(1)
! medium lower thinning
d_est = 1.02
w_kb = 2.5
case(2)
! strong lower thinning
d_est = 1.03
w_kb = 1.5
case(3)
! High thinning
d_est = 1.04
w_kb = 1.2
end select
! calculation of Weibull-Parameter
call min_dbh_overs(nrmin,dbhmin,taxnr)
call min_dbh_unders(nrmin_us,dbhmin_us, taxnr)
write (7777,*) time, nrmin, nrmin_us, dbhmin, dbhmin_us
bas_help = bas_area
if (thin_stor(i).eq.0) then
wpa = dbhmin
else
wpa = dbhmin_us
end if
if (thin_stor(i).eq.0) then
d63 = mdiam*d_est
else
d63 = mdiam_us * d_est
end if
wpb = (d63 - wpa)/ w_kb
wpc = 2
wpc = 0.8
if ((thin_tysp(i).ne.4) .and. rtargetm.lt.target_help) then
!selection of trees for thinning
do
call random_number(pequal)
tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
flagth = 0
! list of potential thinned tree chorts
counth = 0
id_pot = 0
ipot = 1
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%notviable.eqv. .TRUE.) then
if(flag_mort.eq.0) then
id_pot(ipot)=zeig%coh%ident
ipot=ipot + 1
endif
else if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
dbh_h = zeig%coh%diam
db_l = dbh_h - 0.1*dbh_h
db_u = dbh_h + 0.1*dbh_h
counth = counth +1
if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
id_pot(ipot) = zeig%coh%ident
ipot = ipot + 1
end if
if(counth.gt. 100000) exit
end if
zeig=> zeig%next
END DO ! list of potential thinned tree cohorts
! selecting one equal distributed tree from the list of cohorts
if ((ipot-1).ge.1) then
if((ipot-1).eq.1) then
isel = 1
else
call random_number(pequal)
pequal = ran0(idum)
isel = int(pequal*(ipot-1)) +1
end if
ih = id_pot(isel)
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
if(zeig%coh%notviable.eqv..TRUE.) then
if(flag_mort.eq.0) then
zeig%coh%ntreem = zeig%coh%ntreea
zeig%coh%ntreea=0
zeig%coh%nta=0
coun1=coun1+1
target_help = target_help - zeig%coh%ntreem
endif
else
zeig%coh%ntreea = zeig%coh%ntreea -1
zeig%coh%nta = zeig%coh%nta -1
zeig%coh%ntreem = zeig%coh%ntreem +1
coun1 = coun1 + 1
target_help = target_help - 1
exit
end if
end if
zeig =>zeig%next
END DO ! thinning of one tree
end if
diff = target_help - rtargetm
if(diff.le.0.1) exit
if(coun1.gt.100000) exit
end do ! total thinning
end if ! thintype 1,2,3
END IF ! all thinnings and tending
! adding biomasses to litter pools depending on stage of stand
zeig=>pt%first
do
if(.not.associated(zeig)) exit
if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
h1 = zeig%coh%ntreem*zeig%coh%x_fol
h2 = h2 + h1
zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt
endif
zeig=>zeig%next
enddo
call class_man
END SUBROUTINE target_thinning_OC
source_code/version_2.3_windows/timsort.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* SUBROUTINE *!
!* timsort - for sorting of harvested timber to *!
!* different timber qualities *!
!* definition: *!
!* ste - stems *!
!* sg1/sg2 - stem segments *!
!* in1/in2 - industrial wood *!
!* fue - fuelwood *!
!* Subroutine: *!
!* out_tim - generating field sort *!
!* out_timlist *!
!* fuction rabf *!
!* *!
!* 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 timsort
use data_stand
use data_species
use data_tsort
use data_simul
use data_par
use data_manag
!***************** wpm ************************
use data_wpm
use data_stand
use data_species
!***************** wpm ************************
implicit none
integer i
real h,dbh,db, dcrb, hbo, llazmin, lx,lxz,lldmin ,llasdmin,liszmin,help,lisdmin,llzmin
real h1,h2,dcrb_org,db_org,suml, h_org, sumbio_help,sumvol, h3
real (KIND = dg) :: calcvol
character(4) K
character(2) standt
integer count,count_old,taxid
real, external :: rabf
real :: diam_base=0. ! diameter at basis
llazmin=0.;lx=0.;lldmin=0.;llasdmin=0.; liszmin=0.;lisdmin=0.;llzmin=0.
count = 1
sumbio_help = 0
DO i=1,nspec_tree
zeig => pt%first
DO
IF (.NOT. ASSOCIATED(zeig)) EXIT
! Douglasie --- > Kiefer
taxid = zeig%coh%species
IF(taxid .EQ. 10) taxid = 3
IF(taxid .EQ. i) THEN
calcvol = 0.
sumvol = 0.
h = zeig%coh%height -stoh(i) ! stump correction
hbo = zeig%coh%x_hbole
dbh = zeig%coh%diam
dcrb = zeig%coh%dcrb
dcrb_org = dcrb
suml = stoh(i) ! hight of stump; stockhhe
h_org = zeig%coh%height
! calculation of stump biomass for harvesting
! selection of small trees with out dbh
! IF (dbh.eq.0.) THEN
!litter
diam_base= sqrt((zeig%coh%x_ahb+zeig%coh%asapw)*4/pi)
if(hbo.ne.0) then
db = dcrb + (hbo-stoh(i))*(diam_base-dcrb)/hbo
else if (hbo.eq.0)then
db = diam_base*h/(h + stoh(i))
end if
db_org = db
if( hbo.eq.0) then
llzmin=0
lx=0
end if
! stems
help = rabf(i,lzmin(i))
if(db.ge.(lzmin(i) + rabf(i,lzmin(i)))) then
! calculation of lenght at diameter = lzmin
! llzmin > h can occur
if ( dcrb .gt. lzmin(i)+rabf(i,lzmin(i))) then
llzmin = h-(h-hbo)*(lzmin(i)+rabf(i,lzmin(i)))/dcrb
else if (dcrb.le.lzmin(i)+rabf(i,lzmin(i))) then
llzmin = hbo -hbo*(lzmin(i)+rabf(i,lzmin(i))-dcrb)/(db-dcrb)
end if
! calculation of diameter at llzmin/2
if (llzmin/2.lt. hbo) then
lx = dcrb + (db-dcrb)*(hbo-llzmin/2)/hbo
else
lx = dcrb*(h- llzmin/2)/(h-hbo)
end if
! begin of sorting stem lumbre
if( flag_sort.eq.0) then
if( llzmin.ge. lmin(i).and. lx.ge. ldmin(i)+rabf(i,ldmin(i))) then
k = 'ste'
if(zeig%coh%ntreem.ne.0.) then
standt='ab'
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreem,&
standt,h,hbo,db,dcrb,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt='vb'
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreea,&
standt,h,hbo,db,dcrb,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
standt='tb'
flag_deadsort = 1
call out_tim(count,i,k,llzmin + zug ,lx,lzmin(i), zeig%coh%ntreed,&
standt,h,hbo,db,dcrb,calcvol)
end if
end if
suml = suml +llzmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h-llzmin-zug
if(h.lt.0.) h = 0.
db = lzmin(i)
else
! Hhe Hx berechnen, wo lx = ldmin(i)+rabf(i,ldmin(i)) , dann testen, ob 2*Hx >= lmin ist,
!wenn ja, abspeichern in sto
if (dcrb .gt.(ldmin(i)+rabf(i,ldmin(i)))) then
lldmin= h-(h-hbo)*(ldmin(i)+rabf(i,ldmin(i)))/dcrb
else if (dcrb.le.(ldmin(i)+rabf(i,ldmin(i)))) then
lldmin = hbo - hbo*(ldmin(i)+rabf(i,ldmin(i))-dcrb)/(db_org-dcrb)
end if
! calculation of diameter at 2*lldmin and test against lzmin
!(Durchmesser an der Spitze des Stammstckes
if (2*lldmin.lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-lldmin*2)/hbo
else
lx = dcrb*(h- lldmin*2)/(h-hbo)
end if
if (2*lldmin .ge. lmin(i) .and. lx.ge.lzmin(i)+rabf(i,lx)) then
!second test Stammholz!
k = 'ste'
if(zeig%coh%ntreem.ne.0.) then
standt='ab'
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt='vb'
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
standt='tb'
flag_deadsort = 1
call out_tim(count,i,k,2*lldmin + zug,ldmin(i),lx, zeig%coh%ntreed,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lldmin + zug
sumvol=sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lldmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
end if ! db.gt. lzmin(i)
end if ! flag_sortst
! ende test Stammholz
! begin test auf Stammstcke
! calculation of length at diamter lazmin(i)
Do while((h.ge.lasfixl1(i).or.h.ge.lasfixl2(i)).and. hbo.ne.0.)
count_old = count
IF(db .ge.laszmin(i)+rabf(i,laszmin(i)).and. db.gt.lasdmin(i)+ rabf(i,lasdmin(i))) THEN
if (dcrb .eq.0.) then
llazmin = (h_org-suml) -(h_org-stoh(i))*(laszmin(i)+rabf(i,laszmin(i)))/db
else if ( dcrb .gt. laszmin(i)+rabf(i,laszmin(i))) then
llazmin =(h_org-suml)- (h_org-hbo)*(laszmin(i)+rabf(i,laszmin(i)))/dcrb
else if (dcrb.le.laszmin(i)+rabf(i,laszmin(i))) then
llazmin = (hbo-suml) -hbo*(laszmin(i)+rabf(i,laszmin(i))-dcrb)/(db_org-dcrb)
end if
if(llazmin.ge.lasfixl1(i)) then
if(flag_sort.eq.2) then
llazmin = lasfixl2(i)
else
llazmin = lasfixl1(i)
end if
else if(llazmin.ge.lasfixl2(i)) then
llazmin = lasfixl2(i)
end if
! calculation of diameter lx at llazmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-(suml+llazmin/2))/(h_org-stoh(i))
else if ((suml+llazmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin/2))/hbo
else
lx = dcrb*(h_org-(suml+ llazmin/2))/(h_org-hbo)
end if
! calculation of diameter at llazmin
if (dcrb .eq.0.) then
lxz = db*(h_org-(suml+llazmin))/(h_org-stoh(i))
else if ((suml+llazmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin))/hbo
else
lxz = dcrb*(h_org-(suml+ llazmin))/(h_org-hbo)
end if
! test
help = lasdmin(i)+rabf(i,lasdmin(i))
! if flag_sort = 2 only lasfixl2 is used
h3 = lasdmin(i)+rabf(i,lasdmin(i))
if (llazmin.ge. lasfixl1(i).and. lx.ge. lasdmin(i)+rabf(i,lasdmin(i)).and. flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreed,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol +calcvol
calcvol =0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
! test
! if flag_sort = 3 only lasfixl1 is unsed
else if (llazmin.ge.lasfixl2(i).and.llazmin.lt.lasfixl1(i).and. lx.ge. lasdmin(i)+rabf(i,lasdmin(i)).and. flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreem,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreea,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreed,standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
llasdmin = (h_org-suml)- (h_org-stoh(i))*(lasdmin(i)+rabf(i,lasdmin(i)))/db
else if (dcrb .gt. lasdmin(i)+rabf(i,lasdmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(lasdmin(i)+rabf(i,lasdmin(i)))/dcrb
else if (dcrb.le.lasdmin(i)+rabf(i,lasdmin(i))) then
llasdmin = (hbo-suml)-hbo*(lasdmin(i)+rabf(i,lasdmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.lasfixl1(i)) then
llasdmin = lasfixl1(i)/2.
else if(2*llasdmin.ge.lasfixl2(i)) then
llasdmin = lasfixl2(i)/2.
end if
!calculation lx diameter at 2*llasdmin
if (dcrb .eq.0.) then
lx = db*(h_org-suml-llasdmin*2)/(h_org-stoh(i))
else if ((suml+2*llasdmin).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+2*llasdmin))/hbo
else
lx = dcrb*(h_org- suml-2*llasdmin)/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
if(2*llasdmin.ge.lasfixl1(i).and.lx.ge.laszmin(i)+rabf(i,laszmin(i)).and. flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
! if flag_sort = 3 only lasfixl1 is unsed
else if (2*llasdmin.ge.lasfixl2(i).and.2*llasdmin.lt.lasfixl2(i) .and.lx.ge.laszmin(i)+rabf(i,laszmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,lasdmin(i),lx, zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db.gt. laszmin(i)
if(count.eq.count_old) exit
END DO
! end test
! assortment LAS1a for pine
Do while((h.ge.lasfixl1(i).or.h.ge.lasfixl2(i)).and.i.eq.3)
count_old = count
IF(db .ge.las1zmin(i)+rabf(i,las1zmin(i)).and. db.gt.las1dmin(i)+ rabf(i,las1dmin(i))) THEN
if (dcrb .eq.0.) then
llazmin = (h_org-suml) -(h_org-stoh(i))*(las1zmin(i)+rabf(i,las1zmin(i)))/db
else if ( dcrb .gt. las1zmin(i)+rabf(i,las1zmin(i))) then
llazmin =(h_org-suml)- (h_org-hbo)*(las1zmin(i)+rabf(i,las1zmin(i)))/dcrb
else if (dcrb.le.las1zmin(i)+rabf(i,las1zmin(i))) then
llazmin = (hbo-suml) -hbo*(las1zmin(i)+rabf(i,las1zmin(i))-dcrb)/(db_org-dcrb)
end if
if(llazmin.ge.lasfixl1(i)) then
if(flag_sort.eq.2) then
llazmin = lasfixl2(i)
else
llazmin = lasfixl1(i)
end if
else if(llazmin.ge.lasfixl2(i)) then
llazmin = lasfixl2(i)
end if
! calculation of diameter lx at llazmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-(suml+llazmin/2))/(h_org-stoh(i))
else if ((suml+llazmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin/2))/hbo
else
lx = dcrb*(h_org-(suml+ llazmin/2))/(h_org-hbo)
end if
! calculation of diameter at llazmin
if (dcrb .eq.0.) then
lxz = db*(h_org-(suml+llazmin))/(h_org-stoh(i))
else if ((suml+llazmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+llazmin))/hbo
else
lxz = dcrb*(h_org-(suml+ llazmin))/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
help = las1dmin(i)+rabf(i,las1dmin(i))
if (llazmin.ge. lasfixl1(i).and. lx.ge. las1dmin(i)+rabf(i,las1dmin(i)).and.flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreem,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin+zug,lx,lxz, zeig%coh%ntreea,&
standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol +calcvol
calcvol =0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
! if flag_sort = 3 only lasfixl1 is used
else if (llazmin.ge.lasfixl2(i).and.llazmin.lt.lasfixl1(i).and. lx.ge. las1dmin(i)+rabf(i,las1dmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreem,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,llazmin + zug,lx, lxz,zeig%coh%ntreea,standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,llazmin+zug,lx,lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + llazmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - llazmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
llasdmin = (h_org-suml)- (h_org-stoh(i))*(las1dmin(i)+rabf(i,las1dmin(i)))/db
else if (dcrb .gt. las1dmin(i)+rabf(i,las1dmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(las1dmin(i)+rabf(i,las1dmin(i)))/dcrb
else if (dcrb.le.las1dmin(i)+rabf(i,las1dmin(i))) then
llasdmin = (hbo-suml)-hbo*(las1dmin(i)+rabf(i,las1dmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.lasfixl1(i)) then
llasdmin = lasfixl1(i)/2.
else if(2*llasdmin.ge.lasfixl2(i)) then
llasdmin = lasfixl2(i)/2.
end if
!calculation lx diameter at 2*llasdmin
if (dcrb .eq.0.) then
lx = db*(h_org-suml-llasdmin*2)/(h_org-stoh(i))
else if ((suml+2*llasdmin).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+2*llasdmin))/hbo
else
lx = dcrb*(h_org- suml-2*llasdmin)/(h_org-hbo)
end if
! if flag_sort = 2 only lasfixl2 is used
if(2*llasdmin.ge.lasfixl1(i).and.lx.ge.las1zmin(i)+rabf(i,las1zmin(i)).and.flag_sort.ne.2) then
k = 'sg1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
! if flag_sort = 3 only lasfixl1 is used
else if (2*llasdmin.ge.lasfixl2(i).and.2*llasdmin.lt.lasfixl2(i) .and.lx.ge.las1zmin(i)+rabf(i,las1zmin(i)).and.flag_sort.ne.3) then
k = 'sg2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx, zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*llasdmin + zug,las1dmin(i),lx,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
suml = suml + 2*llasdmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
h = h - 2*llasdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db.gt. laszmin(i)
if(count.eq.count_old) exit
END DO
! end test LAS1a for pine
! begin test industrial wood
Do while((h.ge.isfixl1(i).or.h.ge.isfixl2(i)).and.hbo.ne.0)
count_old = count
IF(db.gt.iszmin(i)+rabf(i,iszmin(i)).and.db.gt. isdmin(i)+rabf(i,isdmin(i))) THEN
help = iszmin(i)+rabf(i,iszmin(i))
! calculation of length at diameter iszmin(i)
if (dcrb .eq.0.) then
liszmin = h_org -suml -(h_org-stoh(i))*(iszmin(i)+rabf(i,iszmin(i)))/db
else if ( dcrb .gt. iszmin(i)+rabf(i,iszmin(i))) then
liszmin = (h_org-suml)- (h_org-hbo)*(iszmin(i)+rabf(i,iszmin(i)))/dcrb
else if (dcrb.le.(i)+rabf(i,iszmin(i))) then
liszmin = (hbo-suml) -hbo*(iszmin(i)+rabf(i,iszmin(i))-dcrb)/(db_org-dcrb)
end if
if(liszmin.ge.isfixl1(i)) then
liszmin = isfixl1(i)
else if (liszmin.ge.isfixl2(i)) then
liszmin = isfixl2(i)
end if
! calculation of diameter lx at liszmin/2
if (dcrb .eq.0.) then
lx = db*(h_org-suml-liszmin/2)/(h_org-stoh(i))
else if ((suml+liszmin/2).lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-(suml+liszmin/2))/hbo
else
lx = dcrb*(h_org-suml- liszmin/2)/(h_org-hbo)
end if
! calculation of diameter at liszmin
if (dcrb .eq.0.) then
lxz = db*(h_org-suml-liszmin)/(h_org-stoh(i))
else if ((suml+liszmin).lt. hbo) then
lxz = dcrb + (db_org-dcrb)*(hbo-(suml+liszmin))/hbo
else
lxz = dcrb*(h_org-(suml+ liszmin))/(h_org-hbo)
end if
! test industrial wood Fix length 1
if (liszmin.ge. isfixl1(i).and. lx.ge. isdmin(i)+rabf(i,isdmin(i))) then
k = 'in1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + liszmin + zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - liszmin-zug
if(h.lt.0.) h = 0.
db = lxz
! test industrial wood fix length 2
else if (liszmin.ge.isfixl2(i).and.liszmin.lt.isfixl1(i).and. lx.ge. isdmin(i)+rabf(i,isdmin(i))) then
k = 'in2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,liszmin + zug,lx, lxz,zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + liszmin + zug
sumvol = sumvol + calcvol
calcvol =0.
count = count + 1
h = h - liszmin-zug
if(h.lt.0.) h = 0.
db = lxz
else
if (dcrb.eq.0) then
h1 = isdmin(i)
h2 = rabf(i,isdmin(i))
llasdmin = h_org - suml-(h_org-stoh(i))*(isdmin(i)+rabf(i,isdmin(i)))/db
else if (dcrb .gt. isdmin(i)+rabf(i,isdmin(i))) then
llasdmin = (h_org-suml)-(h_org-hbo)*(isdmin(i)+rabf(i,isdmin(i)))/dcrb
else if (dcrb.le.isdmin(i)+rabf(i,isdmin(i))) then
llasdmin = hbo-suml -hbo*(isdmin(i)+rabf(i,isdmin(i))-dcrb)/(db_org-dcrb)
end if
if(2*llasdmin.ge.isfixl1(i)) then
llasdmin = isfixl1(i)/2.
else if (2*llasdmin.ge.isfixl2(i)) then
llasdmin = isfixl2(i)/2.
end if
!calculation lx diameter at 2*lisdmin
if (dcrb .eq.0.) then
lx = db*(h_org -suml -llasdmin*2)/(h_org -stoh(i))
else if (2*llasdmin.lt. hbo) then
lx = dcrb + (db_org-dcrb)*(hbo-suml-llasdmin)/hbo
else
lx = dcrb*(h_org- llasdmin)/(h_org-hbo)
end if
! test isfixl1
if(2*lisdmin.ge.isfixl1(i).and.lx.ge.iszmin(i)+rabf(i,iszmin(i))) then
k = 'in1'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i), zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lisdmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lisdmin-zug
if(h.lt.0.) h = 0.
db = lx
! test isfixl2
else if (2*lisdmin.ge.isfixl2(i).and.2*lisdmin.lt.isfixl2(i) .and.lx.ge.iszmin(i)+rabf(i,iszmin(i))) then
k = 'in2'
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreem, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
standt = 'vb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreea, standt,h,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,2*lisdmin+zug,lx, isdmin(i),zeig%coh%ntreed, standt,h,hbo,db,dcrb_org,calcvol)
end if
end if
suml = suml + 2*lisdmin+zug
sumvol = sumvol + calcvol
calcvol = 0.
count = count + 1
h = h - 2*lisdmin-zug
if(h.lt.0.) h = 0.
db = lx
end if
end if
END IF ! db .ge. iszmin
if(count.eq.count_old) exit
END DO
! ende test industrial wood
! begin fuelwood
if (h.ne.0.and. db .ne.0) then
k = 'fue'
lx=0.
if(zeig%coh%ntreem.ne.0.) then
standt = 'ab'
if (suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
else
! ! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreem, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreea.ne.0.) then
if(suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
help = zeig%coh%x_sap + zeig%coh%x_hrt
else
! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
standt = 'vb'
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreea, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
if(zeig%coh%ntreed.ne.0..and.flag_mg.ne.0.and.zeig%coh%diam.gt.tardiam_dstem) then
if(thin_flag1(1).ge.0) then
if(suml.eq.stoh(i)) then
! calculation of fuel wood in the case of total use of stem for fuel wood
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. ! kg DW/tree ---> m/tree
help = zeig%coh%x_sap + zeig%coh%x_hrt
else
! calculation of fuelwood volume from all stem segments and total volume of stem, error because stump is not considered
calcvol = (zeig%coh%x_sap + zeig%coh%x_hrt)/spar(i)%prhos/1000000. - sumvol ! m/tree
end if
flag_deadsort = 1
standt = 'tb'
call out_tim(count,i,k,h,db, lx, zeig%coh%ntreed, standt,h_org,hbo,db,dcrb_org,calcvol)
end if
end if
count = count + 1
end if
end if
zeig => zeig%next
end do
end do
end subroutine timsort
subroutine out_tim(cou,nr,k, len, d,zapf, anz,standt,h,hbo,db,dcrb,calcvol)
use data_tsort
use data_simul
use data_par
!***************** wpm ************************
use data_wpm
use data_stand
use data_species
use data_manag
!***************** wpm ************************
type(mansort_type) :: mansort_ini
integer nr,cou
real len, d, anz,zapf, volume, v1,v2,r,r1,rc,vhelp
real (KIND = dg) :: calcvol
character(4) k
character(2) standt
type(timber) ::tim_ini
tim_ini%year = time
tim_ini%count= cou
tim_ini%ttype = k
tim_ini%specnr = nr
tim_ini%length = len
tim_ini%dia = d
tim_ini%diaor = d -rabf(nr,d)
if(tim_ini%diaor.lt.0) tim_ini%diaor=0
!calculaiton of volume for stem segment, depending on the charcteristics (cone, 2 cones, or frustum of a cone)
! cone: vol=1./3.(pi*h*r)
! frustum: vol = pi*h(r1+r1*r2+r2)/3
r = db*0.5
r1 = zapf*0.5
rc = dcrb*0.5
if(k.eq. 'ste') then
if((len + 10.).lt.hbo) then
volume =( pi*len*(r*r + r*r1 + r1*r1)/3.)/1000000. ! frustum
else
v1 = pi*(hbo-stoh(nr))*(r*r +r*rc + rc*rc)/3.
v2 = pi*(len-stoh(nr)-hbo)*(rc*rc+ rc*r1 + r1*r1)/3.
volume = (v1+v2)/1000000.
end if
else if (k.eq.'fue'.and.hbo.ne.0)then
if( db.gt.dcrb) then
if(len.lt.(h-hbo)) then
volume = ( pi * len* r*r/3 )/1000000.
else
v1 = pi* (len-h+hbo)*(r*r + r*rc + rc*rc)/3. ! frustum
vhelp = pi*hbo*(r*r + r*rc + rc*rc)/3.
v2 = pi* (h-hbo)* rc*rc/3. ! cone
volume = (v1+v2)/1000000.
end if
else
volume = (pi*len*r*r/3.)/1000000.
end if
else if (k.eq.'fue'.and.hbo.eq.0)then
volume = (pi*len*r*r/3.)/1000000.
else
! stem timber or industrial timber
if(hbo .eq.0.) then
volume = (pi*len*r*r/3.)/1000000.
else
volume = ( pi*len*(r*r +r*r1 + r1*r1)/3.)/1000000.
end if
end if
if( volume.lt.0) then
volume = volume
end if
if(calcvol.eq.0.) then
tim_ini%vol = volume
calcvol = volume
else
tim_ini%vol = volume
end if
tim_ini%zapfd = zapf
tim_ini%zapfdor = zapf - rabf(nr,zapf)
if ( tim_ini%zapfdor.lt.0) tim_ini%zapfdor=0.
tim_ini%tnum = anz
tim_ini%stype = standt
tim_ini%hei_tree = h
tim_ini%hbo_tree = hbo
tim_ini%dcrb = dcrb
IF (.not. associated(st%first)) THEN
ALLOCATE (st%first)
st%first%tim = tim_ini
NULLIFY(st%first%next)
anz_list = 1
ELSE
ALLOCATE(ztim)
ztim%tim = tim_ini
ztim%next => st%first
st%first => ztim
anz_list = anz_list +1
END IF
!***************** wpm ************************
! information needed for wpm
if ( flag_wpm > 0 .and. (tim_ini%stype .eq. 'ab'.or.tim_ini%stype .eq. 'tb')) then
if (flag_manreal.eq.1.and.maninf.ne.'tending'.and.maninf.ne.'brushing') then
mansort_ini%year = tim_ini%year
mansort_ini%count = tim_ini%count
mansort_ini%spec = tim_ini%specnr
mansort_ini%typus = tim_ini%ttype
mansort_ini%diam = tim_ini%dia
mansort_ini%diam_wob = tim_ini%diaor
mansort_ini%volume = (tim_ini%vol/kpatchsize)*10000. ! m/patchsize ---> m3/ha
mansort_ini%dw = (tim_ini%vol/kpatchsize)*10000*spar(tim_ini%specnr)%prhos*1000000.*cpart ! m/patchsize ---> kg C/ha
mansort_ini%number = tim_ini%tnum
if (.not. associated(first_mansort)) then
allocate (first_mansort)
first_mansort%mansort = mansort_ini
nullify(first_mansort%next)
else
! build new mansort object
allocate(act_mansort)
act_mansort%mansort = mansort_ini
! chain into the list
act_mansort%next => first_mansort
! set the first pointer to the new object
first_mansort => act_mansort
end if
end if
end if
! information needed for sea or wpm+sea
if ( (flag_wpm == 2 .or. flag_wpm == 3) .and. tim_ini%stype .eq. 'vb') then
mansort_ini%year = tim_ini%year
mansort_ini%count = tim_ini%count
mansort_ini%spec = tim_ini%specnr
mansort_ini%typus = tim_ini%ttype
mansort_ini%diam = tim_ini%dia
mansort_ini%diam_wob = tim_ini%diaor
mansort_ini%volume = (tim_ini%vol/kpatchsize)*10000. ! m/patchsize ---> m3/ha
mansort_ini%dw = (tim_ini%vol/kpatchsize)*10000*spar(tim_ini%specnr)%prhos*1000000.*cpart ! m/patchsize ---> kg C/ha
mansort_ini%number = tim_ini%tnum
if (.not. associated(first_standsort)) then
allocate (first_standsort)
first_standsort%mansort = mansort_ini
nullify(first_standsort%next)
else
! build new mansort object
allocate(act_standsort)
act_standsort%mansort = mansort_ini
! chain into the list
act_standsort%next => first_standsort
! set the first pointer to the new object
first_standsort => act_standsort
end if
end if
!***************** wpm ************************
end subroutine out_tim
subroutine out_timlist
use data_tsort
use data_simul
integer timunit
timunit = getunit()
open (timunit,file = 'timlist.dat', status='unknown')
write( timunit,*) ' year ','count ',' spec','type ',' length',' diameter', 'diam wo bark', 'top diam. ',' top d. wo bark','Volume(m)',' number '
ztim=>st%first
do
IF (.not.ASSOCIATED(ztim)) exit
write(timunit,'(3I6,1x,A5,1x,F8.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f7.3,1x,f10.2)') ztim%tim%year, ztim%tim%count, &
ztim%tim%specnr,ztim%tim%ttype,ztim%tim%length,ztim%tim%dia,ztim%tim%diaor,ztim%tim%zapfd,ztim%tim%zapfdor, ztim%tim%vol,ztim%tim%tnum
ztim=>ztim%next
end do
close(timunit)
end subroutine out_timlist
real function rabf(spec, db)
! calculation of rabz i.A.
use data_tsort
use data_species
integer iz, spec
do iz = 1,nspec_tree
if(iz.eq.spec) then
if(db.lt.rabth(spec,1)) then
rabf = rabz(spec,1)
else if (db.ge.rabth(spec,1).and. db.lt.rabth(spec,2)) then
rabf = rabz(spec,2)
else
rabf = rabz(spec,3)
end if
end if
end do
end function rabf
\ No newline at end of file
source_code/version_2.3_windows/tool.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* Subroutines for standard tasks *!
!* *!
!* contains: *!
!* SOLV_QUADR solving quadratic equation, real*4 *!
!* DSOLV_QUADR solving quadratic equation, real*8 *!
!* NEWT Newton method *!
!* TRICOF Harmonic Analysis *!
!* SORT_INDEX Sorts two arrays *!
!* SORT sort an array by quicksort method *!
!* MOMENT Descriptive statistics of a data set *!
!* *!
!* 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 solv_quadr (meth, p, q, x1, x2, res1, res2, rnum)
! Solution of quadratic equation in normal form
! x*x + p*x + q = 0
IMPLICIT NONE
! Input
integer meth ! solver method
real p, q ! parameter of the quadratic equation
! Output
real x1, x2 ! solutions; initial value of Newton method
real res1, res2 ! residua
integer rnum ! return code
real discr
! Variables of Newton method
real df ! quotient of quadratic function and its first derivative
real :: precision = 1E-5
integer:: maxloop, iloop
! Variables of solver program ZPORC of ISML Library
!external ZPLRC
discr = (p*p/4.)-q
if (discr .lt. 0.) then
rnum = -1 ! no real solution
return
else
select case (meth)
case (1)
! standard solution
discr = SQRT(discr)
x1 = -p/2. + discr
x2 = -p/2. - discr
rnum = 0
case (2)
! Vieta's formulae (root theorem)
discr = SQRT(discr)
x2 = -p/2. - discr
x1 = q / x2
rnum = 0
case (3)
! Newton method
! initial value x2
maxloop = 100
iloop = 1
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
do while (abs(df) .gt. precision .and. iloop .le. maxloop)
x2 = x2 - df
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
iloop = iloop + 1
enddo
if (iloop .lt. maxloop) then
rnum = 0
else
rnum = 1
endif
end select
endif ! discr
res1 = x1*x1 + p*x1 + q
res2 = x2*x2 + p*x2 + q
END SUBROUTINE solv_quadr
!**************************************************************
SUBROUTINE dsolv_quadr (meth, p, q, x1, x2, res1, res2, rnum)
! Solution of quadratic equation in normal form
! x*x + p*x + q = 0
! with double precision
IMPLICIT NONE
! Input
integer meth ! solver method
real (kind (0.0D0)) :: p, q ! parameter of the quadratic equation
! Output
real (kind (0.0D0)) :: x1, x2 ! solutions
real (kind (0.0D0)) :: res1, res2 ! residua
integer rnum ! return code
real (kind (0.0D0)) :: discr
! Variables of Newton method
real (kind (0.0D0)) :: df ! quotient of quadratic function and its first derivative
real (kind (0.0D0)) :: precision = 1E-5
integer:: maxloop, iloop
! Variables for solver program ZPORC of ISML Library
real (kind (0.0D0)) :: coeff(3)
real (kind (0.0D0)) :: zero(2)
discr = (p*p/4.)-q
if (discr .lt. 0.) then
rnum = -1 ! no real solution
return
else
select case (meth)
case (1)
! standard solution
discr = DSQRT(discr)
x1 = -p/2. + discr
x2 = -p/2. - discr
rnum = 0
case (2)
! Vieta's formulae (root theorem)
discr = DSQRT(discr)
x2 = -p/2. - discr
x1 = q / x2
rnum = 0
case (3)
! Newton method
! initial value x2
maxloop = 100
iloop = 1
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
do while (abs(df) .gt. precision .and. iloop .le. maxloop)
x2 = x2 - df
df = (x2*x2 + p*x2 + q) / (2.* x2 + p)
iloop = iloop + 1
enddo
if (iloop .lt. maxloop) then
rnum = 0
else
rnum = 1
endif
end select
endif ! discr
res1 = x1*x1 + p*x1 + q
res2 = x2*x2 + p*x2 + q
END SUBROUTINE dsolv_quadr
!**************************************************************
SUBROUTINE newt (x, f, df, ddf, prec, maxit, rnum)
! Newton method
implicit none
integer :: maxit ! maximum number of iteration
real :: x ! initial value and result
real :: prec ! precision
real :: dx ! quotient of function and its first derivative
real :: hf, hdf, hddf
integer :: rnum ! options: 0 - change of sign allowed,
! 1 - no change of sign
! return code
integer :: i
real, external :: f, df, ddf ! function and its first derivative
hf = f(x)
hdf = df(x)
hddf = ddf(x)
dx = hddf * hf
if (abs(dx) .lt. abs(hdf*hdf)) then ! Test of convergence
! Iteration
i = 1
if (abs(dx) .gt. 0.) then
dx = hf / hdf
endif
do while (abs(hf) .gt. prec .and. i .le. maxit)
if (dx .gt. x .and. rnum .gt. 0) dx = x/2.
x = x - dx
hdf = df(x)
if (abs(hdf) .gt. 0.) then
hf = f(x)
dx = hf/hdf
endif
i = i + 1
enddo
if (i .lt. maxit) then
rnum = 0
else
rnum = 1 ! not enough iteration steps
endif
else
rnum = -1 ! no convergence
endif
END SUBROUTINE newt
!**************************************************************
SUBROUTINE TRICOF(F,NF,A,NE,B,NO,IOP)
! PURPOSE = TO COMPUTE THE COEFFICIENTS IN A TRIGONOMETRIC EXPANSION
! FOR A FUNCTION GIVEN IN EQUIDISTANT POINTS
! PARAMETERS
! F = AN ARRAY USED FOR STORING THE FUNCTION VALUES F(X)
! NF = THE NUMBER OF FUNCTION VALUES IN THE ARRAY F.NF MUST
! HAVE THE STRUCTURE , NF = 2*N+1 , THE GENERAL CASE
! NF = N+1 , THE EVEN CASE
! NF = N-1 , THE ODD CASE
! A = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE COS-
! INE TERMS
! NE = THE NUMBER OF COEFFICIENTS IN THE ARRAY A , NE = N+1
! B = AN ARRAY USED FOR RETURNING THE COEFFICIENTS OF THE SINE
! TERMS
! NO = THE NUMBER OF COEFFICIENTS IN THE ARRAY B , NO = N-1
! IOP = OPTION NUMBER , IOP = 1 , THE GENERAL CASE
! IOP = 2 , THE EVEN CASE
! IOP = 3 , THE ODD CASE
DIMENSION F(NF) , A(NE) , B(NO)
REAL KSI0 , KSI1 , KSIK
DATA ZERO , FOURTH , HALF , ONE , TWO , PI / 0. , .25 , .5 , 1. , 2. , 3.14159265358979 /
! COMPUTE THE NUMBER N (SEE EXPLANATION OF PARAMETERS)
1000 N=0
IF (IOP.EQ.1) N=(NF-1)/2
IF (IOP.EQ.2) N=NF-1
IF (IOP.EQ.3) N=NF+1
IF (N.EQ.0) STOP
! STOP IF IOP DOES NOT HAVE A CORRECT VALUE
IF (IOP.GT.1) GO TO 1030
IF ((2*N-NF+1).NE.0) STOP
! STOP IF NF DOES NOT HAVE THE CORRECT STRUCTURE IN THE GENERAL CASE
! SPLIT THE FUNCTION F(X) IN AN EVEN AND ODD PART
M=N+1
DO 1020 J=1,N
COF1=HALF*(F(M+J)+F(M-J))
COF2=HALF*(F(M+J)-F(M-J))
F(M+J)=COF2
F(M-J)=COF1
1020 CONTINUE
! REWRITE N IN POWERS OF 2 I.E. N=NBASE*2**NEXP
1030 NBASE=N
NEXP =0
1040 NINT =NBASE/2
IF ((NBASE-2*NINT).NE.0) GO TO 1050
NBASE=NINT
NEXP =NEXP+1
GO TO 1040
! DO SOME INITIAL CALCULATIONS
1050 REALN=NBASE
ARG =HALF*PI/REALN
KSI0 =COS(ARG)
ETA0 =SIN(ARG)
! START CALCULATION OF COEFFICIENTS
IF (IOP.EQ.3) GO TO 1160
! ********** EVEN COEFFICIENT CALCULATION **********
! COMPUTE THE BASIC COEFFICIENTS A(K) , K=1(1)(NBASE+1)
! START CALCULATION OF A(1)
NN =NBASE-1
NPOINT=1
NINCRE=2**NEXP
NLOCAL=NINCRE+1
BASEIN=ONE/REALN
A(1) =HALF*(F(1)+F(N+1))
IF (NN.EQ.0) GO TO 1065
DO 1060 J=1,NN
A(1) =A(1)+F(NLOCAL)
NLOCAL=NLOCAL+NINCRE
1060 CONTINUE
1065 A(1) =TWO*BASEIN*A(1)
! START CALCULATION OF A(K) , K=2(1)(NBASE+1)
KSI1=KSI0
KSIK=KSI1
ETA1=ETA0
ETAK=ETA1
CONST=HALF*F(N+1)
DO 1090 K=1,NBASE
COF1=TWO*(TWO*KSIK**2-ONE)
A2 =ZERO
A1 =A2
A0 =CONST
NLOCAL=N+1-NINCRE
DO 1070 J=1,NBASE
A2=A1
A1=A0
A0=F(NLOCAL)+COF1*A1-A2
NLOCAL=NLOCAL-NINCRE
1070 CONTINUE
1080 A(K+1)=BASEIN*(A0-A2)
COF1 =KSIK
COF2 =ETAK
KSIK =KSI1*COF1-ETA1*COF2
ETAK =ETA1*COF1+KSI1*COF2
1090 CONTINUE
! CALCULATION OF THE BASIC EVEN COEFFICIENTS FINISHED
IF (NEXP.EQ.0) GO TO 1145
! CONTINUE CALCULATION OF EVEN COEFFICIENTS
NUMCOF=NBASE
DO 1140 NSTEP=1,NEXP
NINCRE=2**(NEXP-NSTEP)
NPOINT=NINCRE+1
NINCRE=2*NINCRE
NLOCAL=NPOINT
NUMBER=2*NUMCOF+1
! COMPUTE CONSTANT TERM IN MID-POINT APPROXIMATION I.E. K=1
SUM=ZERO
DO 1100 J=1,NUMCOF
SUM=SUM+F(NLOCAL)
NLOCAL=NLOCAL+NINCRE
1100 CONTINUE
SUM =TWO*BASEIN*SUM
COF1=A(1)
A(1)=HALF*(COF1+SUM)
A(NUMBER)=HALF*(COF1-SUM)
IF (NUMCOF.EQ.1) GO TO 1135
! COMPUTE MID-POINT APPROXIMATION FOR K=2(1)NUMCOF
1105 NN =NUMCOF-1
KSIK=KSI1
ETAK=ETA1
DO 1130 K=1,NN
COF1=TWO*(TWO*KSIK**2-ONE)
A2=ZERO
A1=A2
NLOCAL=N+2-NPOINT
A0=F(NLOCAL)
DO 1110 J=1,NN
A2=A1
A1=A0
NLOCAL=NLOCAL-NINCRE
A0=F(NLOCAL)+COF1*A1-A2
1110 CONTINUE
1120 SUM=TWO*BASEIN*(A0-A1)*KSIK
COF1=A(K+1)
A(K+1)=HALF*(COF1+SUM)
A(NUMBER-K)=HALF*(COF1-SUM)
COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1130 CONTINUE
1135 A(NUMCOF+1)=HALF*A(NUMCOF+1)
! CALCULATIONS OF MID-POINT APPROXIMATIONS FINISHED
! DO CHANGES RELATED TO HALVING OF THE INTERVAL
ARG =HALF*ARG
COF1=ETA1
ETA1=SIN(ARG)
KSI1=HALF*COF1/ETA1
BASEIN=HALF*BASEIN
NUMCOF=2*NUMCOF
1140 CONTINUE
1145 IF (NEXP.EQ.0) NUMBER=NBASE+1
A(NUMBER)=HALF*A(NUMBER)
! CALULATION OF EVEN COEFFICIENTS FINISHED
1150 IF (IOP.EQ.2) RETURN
! RETURN TO CALLING PROGRAM IF F(X) WAS AN EVEN FUNCTION
! IF IOP=1 CHANGE SIGN OF EACH SECOND COEFFICIENTS
NINT=(N+1)/2
IF (NINT.EQ.0) GO TO 1166
DO 1164 K=1,NINT
A(2*K)=-A(2*K)
1164 CONTINUE
! ********** ODD COEFFICIENT CALCULATION **********
! COMPUTE THE BASIC COEFFICIENTS B(K) , K=1(1)NBASE
1166 ARG=HALF*PI/REALN
1160 IF (IOP.EQ.1) NMAX=2*N+1
IF (IOP.EQ.3) NMAX=N
NINCRE=2**NEXP
NPOINT=NMAX-NINCRE
NLOCAL=NPOINT
BASEIN=ONE/REALN
B(1)=ZERO
IF (NBASE.EQ.1) GO TO 1200
KSI1=TWO*KSI0**2-ONE
KSIK=KSI1
ETA1=TWO*KSI0*ETA0
ETAK=ETA1
NN =NBASE-1
NNN=NN-1
DO 1190 K=1,NN
COF1=TWO*KSIK
A2 =ZERO
A1 =A2
A0 =F(NPOINT)
NLOCAL=NPOINT-NINCRE
IF (NNN.EQ.0) GO TO 1180
DO 1170 J=1,NNN
A2=A1
A1=A0
A0=F(NLOCAL)+COF1*A1-A2
NLOCAL=NLOCAL-NINCRE
1170 CONTINUE
1180 B(K)=TWO*BASEIN*A0*ETAK
COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1190 CONTINUE
! CALCULATION OF THE BASIC ODD COEFFICIENTS FINISHED
1200 IF (NEXP.EQ.0) GO TO 1260
! CONTINUE CALCULATION OF ODD COEFFICIENTS
KSI1=KSI0
ETA1=ETA0
NUMCOF=NBASE
DO 1250 NSTEP=1,NEXP
KSIK=KSI1
ETAK=ETA1
NINCRE=2**(NEXP-NSTEP)
NPOINT=NMAX-NINCRE
NINCRE=2*NINCRE
NUMBER=2*NUMCOF
B(NUMCOF)=ZERO
! COMPUTE MID-POINT APPROXIMATIONS FOR K=1(1)NUMCOF
NN =NUMCOF-1
DO 1240 K=1,NUMCOF
COF1=TWO*(TWO*KSIK**2-ONE)
A2 =ZERO
A1 =A2
NLOCAL=NPOINT
A0 =F(NLOCAL)
IF (NN.EQ.0) GO TO 1220
DO 1210 J=1,NN
A2=A1
A1=A0
NLOCAL=NLOCAL-NINCRE
A0=F(NLOCAL)+COF1*A1-A2
1210 CONTINUE
1220 SUM=TWO*BASEIN*(A0+A1)*ETAK
COF1=B(K)
B(K)=HALF*(COF1+SUM)
IF (K.EQ.NUMCOF) GO TO 1230
B(NUMBER-K)=-HALF*(COF1-SUM)
1230 COF1=KSIK
COF2=ETAK
KSIK=KSI1*COF1-ETA1*COF2
ETAK=ETA1*COF1+KSI1*COF2
1240 CONTINUE
! CALCULATION OF MID-POINT APPROXIMATION FINISHED
! DO CHANGES RELATED TO HALVING OF INTERVAL
ARG =HALF*ARG
COF1=ETA1
ETA1=SIN(ARG)
KSI1=HALF*COF1/ETA1
BASEIN=HALF*BASEIN
NUMCOF=2*NUMCOF
1250 CONTINUE
! CALCULATION OF ODD COEFFICIENTS FINISHED
1260 IF (IOP.EQ.3) RETURN
! IF IOP=1 RECOMPUTE FUNCTION VALUES
DO 1270 J=1,N
COF2=F(M+J)
COF1=F(M-J)
F(M+J)=COF1+COF2
F(M-J)=COF1-COF2
1270 CONTINUE
RETURN
END SUBROUTINE TRICOF
!**************************************************************
SUBROUTINE sort_index(n,arr,brr)
! variation of sort2 for integer array
! sorts array arr(1:n) into an ascending order and
! makes the corresponding rearrangement of the array brr(1:n)
INTEGER n,M,NSTACK
Integer arr(n)
INTEGER brr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,b,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
b=brr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
brr(i+1)=brr(i)
11 continue
i=0
2 arr(i+1)=a
brr(i+1)=b
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
temp=brr(k)
brr(k)=brr(l+1)
brr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
temp=brr(l+1)
brr(l+1)=brr(ir)
brr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
temp=brr(l)
brr(l)=brr(ir)
brr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
temp=brr(l+1)
brr(l+1)=brr(l)
brr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
b=brr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
temp=brr(i)
brr(i)=brr(j)
brr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
brr(l)=brr(j)
brr(j)=b
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort2'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END
! (C) Copr. 1986-92 Numerical Recipes Software "!D#+.
!**************************************************************
SUBROUTINE sort(n,arr)
! sort a n-dimensional array arr(1:n) by quicksort method
INTEGER n,M,NSTACK
REAL arr(n)
PARAMETER (M=7,NSTACK=50)
INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
REAL a,temp
jstack=0
l=1
ir=n
1 if(ir-l.lt.M)then
do 12 j=l+1,ir
a=arr(j)
do 11 i=j-1,1,-1
if(arr(i).le.a)goto 2
arr(i+1)=arr(i)
11 continue
i=0
2 arr(i+1)=a
12 continue
if(jstack.eq.0)return
ir=istack(jstack)
l=istack(jstack-1)
jstack=jstack-2
else
k=(l+ir)/2
temp=arr(k)
arr(k)=arr(l+1)
arr(l+1)=temp
if(arr(l+1).gt.arr(ir))then
temp=arr(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
endif
if(arr(l).gt.arr(ir))then
temp=arr(l)
arr(l)=arr(ir)
arr(ir)=temp
endif
if(arr(l+1).gt.arr(l))then
temp=arr(l+1)
arr(l+1)=arr(l)
arr(l)=temp
endif
i=l+1
j=ir
a=arr(l)
3 continue
i=i+1
if(arr(i).lt.a)goto 3
4 continue
j=j-1
if(arr(j).gt.a)goto 4
if(j.lt.i)goto 5
temp=arr(i)
arr(i)=arr(j)
arr(j)=temp
goto 3
5 arr(l)=arr(j)
arr(j)=a
jstack=jstack+2
if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
if(ir-i+1.ge.j-l)then
istack(jstack)=ir
istack(jstack-1)=i
ir=j-1
else
istack(jstack)=j-1
istack(jstack-1)=l
l=i
endif
endif
goto 1
END
! C (C) Copr. 1986-92 Numerical Recipes Software )$!.
!**************************************************************
SUBROUTINE moment(array,n,ave,adev,sdev,var,skew,curt)
! Calculates statistics of array (n-dimensional array of data)
! n - number of observations
! adev - average deviation
! ave - average
! curt - curtosis
! sdev - standard deviation
! skew - skewness
! var - variance
INTEGER n
REAL adev,ave,curt,sdev,skew,var,array(n)
INTEGER j
REAL p,s,ep
if(n.le.1)pause 'n must be at least 2 in moment'
s=0.
do 11 j=1,n
s=s+array(j)
11 continue
ave=s/n
adev=0.
var=0.
skew=0.
curt=0.
ep=0.
do 12 j=1,n
s=array(j)-ave
ep=ep+s
adev=adev+abs(s)
p=s*s
var=var+p
p=p*s
skew=skew+p
p=p*s
curt=curt+p
12 continue
adev=adev/n
var=(var-ep**2/n)/(n-1)
sdev=sqrt(var)
if(var.ne.0.)then
skew=skew/(n*sdev**3)
curt=curt/(n*var**2)-3.
else
skew = -99.
curt = -99.
endif
return
END
! (C) Copr. 1986-92 Numerical Recipes Software )$!.
FUNCTION rtbis(func,x1,x2,xacc)
INTEGER JMAX
REAL rtbis,x1,x2,xacc,func
EXTERNAL func
PARAMETER (JMAX=40)
INTEGER j
REAL dx,f,fmid,xmid
fmid=func(x2)
f=func(x1)
if(f.lt.0.)then
rtbis=x1
dx=x2-x1
else
rtbis=x2
dx=x1-x2
endif
do j=1,JMAX
dx=dx*.5
xmid=rtbis+dx
fmid=func(xmid)
if(fmid.le.0.)rtbis=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
end do
END function
\ No newline at end of file
source_code/version_2.3_windows/tool1.f 0 → 100644
View file @ 82df1085
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for standard tasks *!
!* *!
!* contains: *!
!* DAINTZ Date to day of the year *!
!* TZINDA Day of the year to date *!
!* TAB_INT Table function *!
!* CHARACTER_IN_INTEGER Conversion of character in integer *!
!* INTEGER_IN_CHARACTER Conversion of integer in character *!
!* QUANTILE calculates the 0.05 and 0.95 quantile *!
!* QUANT_CALC calculates a quantile of a sorted array *!
!* *!
!* 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 DAINTZ(IT,IM,IJ,TZ)
! Umrechnen von Datum in Tageszaehler
implicit none
INTEGER IT, IM, IJ, TZ
INTEGER I, ME
REAL, DIMENSION(12):: MNL
! COMMON /MONTH/ MMM(12),JAHR,JS,IC
DATA MNL /31,28,31,30,31,30,31,31,30,31,30,31/
TZ=IT
IF (IM.EQ.1) RETURN
ME=IM-1
if (mod(IJ,4).EQ.0) MNL(2)=29
if ((IJ .eq. 1900) .or. (IJ .eq. 1800) .or. (IJ .eq. 1700)) MNL(2)=28
DO I=1,ME
TZ=TZ+MNL(I)
enddo
MNL(2)=28
END SUBROUTINE DAINTZ
!***********************************************************************
SUBROUTINE TZINDA(T,M,J,TZ)
! Umrechnen von Tageszaehler in Datum
implicit none
INTEGER MNL(12)
INTEGER T, M, J, TZ
DATA MNL /31,28,31,30,31,30,31,31,30,31,30,31/
if (mod(J,4).EQ.0) MNL(2)=29
if ((J .eq. 1900) .or. (J .eq. 1800) .or. (J .eq. 1700)) MNL(2)=28
T=TZ
M=1
do while (T .gt. MNL(M))
T=T-MNL(M)
M=M+1
if (M .gt. 12) return
enddo
MNL(2)=28
END SUBROUTINE TZINDA
!***********************************************************************
SUBROUTINE tab_int(x,y,idim,arg,val)
! Read a table function with ordered pairs x,y (sortet)
! linear interpolation between
implicit none
! input
integer idim ! dimension of array x, y
real, dimension(idim) :: x, y ! table values
real arg ! argument of function
! output
real val ! result
integer i
if (arg .le. x(1)) then
val = y(1)
else if (arg .ge. x(idim)) then
val = y(idim)
else
i = 2
do while ((i .lt. idim) .and. (arg .gt. x(i)))
i = i+1
enddo
if (arg .eq. x(i)) then
val = y(i)
else
val = y(i) + (y(i)-y(i-1)) * (arg-x(i)) / (x(i)-x(i-1))
endif
endif
END subroutine tab_int
!***********************************************************************
SUBROUTINE character_in_integer(string, vint)
! Conversion of character variable in integer variable
implicit none
integer vint
character (100) string
character (10) help
write(help,'(A)') string
read(help,*) vint
END subroutine character_in_integer
!**************************************************************
SUBROUTINE integer_in_character(vint, string)
! Conversion of integer variable in character variable
implicit none
integer vint
character (10) string
character (10) help
write(help,'(I10)') vint
read(help,*) string
END subroutine integer_in_character
!**************************************************************
SUBROUTINE quantile(idim, arr, quant05, quant95, median)
! sorts and calculates the 0.05 and 0.95 quantile of an array with dimension idim
implicit none
! input
integer idim ! dimension of array arr
real, dimension(idim) :: arr ! array
! output
real quant05, quant95, median ! 0.05 and 0.95 quantile
call sort(idim,arr)
call quant_calc(idim, arr, 0.05, quant05) ! 0.05 quantile
call quant_calc(idim, arr, 0.95, quant95) ! 0.95 quantile
call quant_calc(idim, arr, 0.5, median) ! 0.95 quantile
END SUBROUTINE quantile
!**************************************************************
SUBROUTINE quant_calc(idim, arr, pord, quant)
! calculates a quantile of a sorted array with dimension idim
implicit none
integer idim ! dimension of array arr
real, dimension(idim) :: arr ! array
real quant ! quantile
real pord, help ! order
integer ihelp
help = idim * pord
ihelp = int(help)
if (ihelp*1.0 .lt. help) then
quant = arr(ihelp+1)
else
quant = (arr(ihelp+1) + arr(ihelp)) / 2.
endif
END SUBROUTINE quant_calc
!**************************************************************