-
Petra Lasch-Born authored
version 2.3
Petra Lasch-Born authoredversion 2.3
initia.f 61.31 KiB
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* - Initialisation of cohorts = *!
!* reads cohort information and calculates missing values *!
!* which are needed for stand initialisation *!
!* initia *!
!* treeini *!
!* sapini *!
!* header *!
!* crown_base *!
!* crown_base_eg *!
!* fdfahc: function *!
!* ini_gener_sap *!
!* NEWTON: function 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 INITIA !
!***********************!
SUBROUTINE INITIA
! begin declaration section
USE data_init
USE data_par
USE data_simul
USE data_species
USE data_stand
use data_help
IMPLICIT none
REAL :: area !area of database in m^2 (10000=1ha)
INTEGER :: area_factor !factor for calculation per patch (=area/kpatchsize)
REAL :: hlp_lai,share, ager
INTEGER :: taxid, & ! species number
age, & ! tree age
n, & ! number of trees
n_koh, & !
k, & ! number of tree classes
ng_locid ! ID stand
INTEGER :: inunit, parunit,outunit,tmpunit,ctrlunit,listunit !units
CHARACTER*85 zeile
CHARACTER*80 :: infile
CHARACTER :: source
INTEGER :: nlines, nlines_comp, istart, fl_num, nhelp, numstand, ihelp
INTEGER :: tax_of_BRA_id
INTEGER,DIMENSION(:),ALLOCATABLE :: locid_comp
REAL rsap, cform, dummy
REAL aux
LOGICAL :: select_lines
real standsz(10000)
CHARACTER*40, allocatable, dimension(:) ::helptmp
INTEGER :: helpz
! Stand data (model initialisation)
INTEGER baum(10),alt(10),klimid,gwa,lbanr,wgeb,lein,zei
REAL mhoe(10),dm(10),gf(10),bon(10),en(10),psi(10)
! Parameters for missing data algorithms
REAL p0(nspec_tree),p1(nspec_tree),p2(nspec_tree),p3(nspec_tree),p4(nspec_tree), &
c1(nspec_tree),c2(nspec_tree),ku_a0(nspec_tree),ku_a1(nspec_tree),ku_a2(nspec_tree),&
ku_b0(nspec_tree),ku_b1(nspec_tree),ku_b2(nspec_tree),ku_c0(nspec_tree),&
ku_c1(nspec_tree),ku_c2(nspec_tree),wei_k1(nspec_tree),wei_k2(nspec_tree)
! ------------------------------------------------------------------
! INTEGER ncl !Number of classes after classification
integer ncl1
REAL dg,dmin,dmax,g,gpatch,b,c,bhd,height,hbc,hg
REAL tot_crown_area, mixed_tot_ca, corr_la
INTEGER pass
REAL saquad, genDg, nbhd,x,gx,bhdmax,bhdmin,clwdth,Fint(0:100)
REAL ku_a,ku_b,ku_c,wei_f,thdmax,p1n,p4n
REAL, allocatable, dimension(:) :: nz
REAL, allocatable, dimension(:) :: zheigh,zbhd,zhbc
REAL xxr,xyr, &
kd, &
h_para,h_parb !parameter of the height function of level II sites
INTEGER idum,anzahl, data_flag,start,baumid,dir_flag,inwahl,bz,imax
INTEGER i,j,anzit,iz,id,icl,ios,xid,xnr,xxi,xyi, &
bhdcl, & !diameter classes level II
dclmin, & !smallest diameter class level II
ndcl, & !amount of diameter classes of level II
dcwdth, & !class wideness diameter classes of level II
n_dc(30) !stem figure of level II diameter class
LOGICAL ehkwei, wfirst, kfirst, optimi
LOGICAL, allocatable, dimension(:) :: smaldc, bigdc
CHARACTER*20 fnam2
CHARACTER*5 datasets
CHARACTER status
real nzsum
! ------------------------------------------------------------------
! ----Function----
REAL ran0
REAL crown_base
real crown_base_eg
! ------------------------------------------------------------------
REAL T
DATA T/7.0/
! ------------------------------------------------------------------
!
! end of declaration section
!******************************************************************************
!ncl1 = 60
ncl1=40
allocate (zheigh(ncl1), zbhd(ncl1), zhbc(ncl1), nz(ncl1))
allocate (smaldc(ncl1), bigdc(ncl1))
print *,' '
print *, ' *** Choice of forest stand data set: '
print *, ' 1 - Datenspeicher Waldfond'
print *, ' 2 - single tree data; classification must be performed (e.g. SILVA data)'
print *, ' 3 - Level2-data'
print *, ' 4 - already existing class file'
print *, ' 5 - FORGRA data'
print *, ' 6 - Bavarian inventory data'
WRITE(*,'(A)',advance='no') ' ***Make your choice: '
READ *, data_flag
print *,' '
clwdth=15 !set diameter class-class width
corr_la=1. !standard value for leaf area correction in stands of high sum of crown projection areas
mixed_tot_ca=0. !sum of crown projection area for mixed stands
pass = 1 !counter for number of passes through calculation loop for mixed standsrsap=0.3 !standard value of rsap for cases where rsap is not determined dynamically
! get unit number and open units used in all of the above cases
ctrlunit=GETUNIT()
WRITE(*,*)site_name(ip)
OPEN (ctrlunit,FILE=TRIM(site_name(ip))//'.initctrl',STATUS='replace')
WRITE(ctrlunit,*)'# number of trees in cohort = n trees'
WRITE(ctrlunit,*)'# age = age'
WRITE(ctrlunit,*)'# height = H'
WRITE(ctrlunit,*)'# height to the base of crown = Hbc'
WRITE(ctrlunit,*)'# breast height diameter = bhd'
WRITE(ctrlunit,*)'# sapwood fraction of trunc cross sectional area at breast height = rsap'
WRITE(ctrlunit,*)'# trunc diameter at tree base = D'
WRITE(ctrlunit,*)'# trunc diameter at crown base = Dc'
WRITE(ctrlunit,*)'# sapwood cross sectional area inside bole = Asap'
WRITE(ctrlunit,*)'# heartwood cross sectional area at crown base = Ahc'
WRITE(ctrlunit,*)'# heartwood cross sectional area at tree base = Ahb'
WRITE(ctrlunit,*)'# Vol for no heartwood in crown space = Vmin'
WRITE(ctrlunit,*)'# Vol prescribed according to empiracal volume function = Vpre'
WRITE(ctrlunit,*)'# stem vol inherent in initialisation = Veff'
WRITE(ctrlunit,'(A150)')'# n trees age H Hbc bhd rsap D Dc Asap Ahc Ahb Vmin Vpre Veff'
outunit=GETUNIT()
OPEN (outunit,FILE=TRIM(treefile(ip)),STATUS='replace')
! ------------------------------------------------------------------
! read in parameter for the missing-data-generator:
! bhd-distribution from Nagel & Biging (1995),
! crown starting height from Nagel (1995), uni-height curve according to Weimann (1980) bzw. Kuleschis (1981)
parunit=GETUNIT()
OPEN (parunit, FILE='input/generreg.par', STATUS='old')
do i=1,nspec_tree
READ (parunit,*) p0(i),p1(i),p2(i),p3(i),p4(i),c1(i),c2(i),ku_a0(i),ku_a1(i),ku_a2(i), &
ku_b0(i),ku_b1(i),ku_b2(i),ku_c0(i),ku_c1(i),ku_c2(i),wei_k1(i),wei_k2(i)
ENDDO
CLOSE(parunit)
! ---------------------------------------------------------------------
inunit=GETUNIT()
SELECT CASE(data_flag)
! ****************************************************************************
! case(1) stand generation if data source is Datenspeicher Waldfond
CASE(1)
print *, ' Forest stand data set: Datenspeicher Waldfond'
! preliminary: here make a choice and compile
! datasets='singl' sets the choice of the old version which uses one single
! set (i.e. the first one in an input file) which contains
! the complete imformation for the stand in one single line
! datasets='multi' sets the choice of a version reading a file with line by
! line information as in the original Datenspeicher and then
! writes a *.ini file for many stands with individual stand
! information separated by lines with stand identifiers
print*, 'choose data set (multi/singl):'
read(*,*) datasets
print*, ' file name (with directory):'
read(*,'(A)') infile
source='D'
standsz = 0.
OPEN (inunit, FILE=TRIM(infile), STATUS='old')
! ------------------------------------------------------------------
! generating standard value out of data from the data storage unit
! based on estimation routine from Nagel und Biging (1995),
! Nagel (1995) und Gerold (1990).
! ------------------------------------------------------------------
!
! The following variables are read from forest inventory data:
! Species(baum),Age(alt),Quadratic Mean Diameter(dm),Height of tree with dm(mhoe),
! Basal area(gf),Yield Class(bon),"Ertragsniveau"(en)
! Additional Site variables:
! Climate station(klimid),distance of groundwater table(gwa),soil type(lbanr),! forest region 'Wuchsgebiet'(wgeb),last management operation(lein), number of tree layers(zei)
! currently not used for initialisation: xid, klimid, gwa, lbanr, wgeb, lein, bon(i), en(i)
! lbanr (check difference to declaration!),
! check if alt and baum can be skipped as variable names and age and species directly used
! check idendity of hg and mhoe, dg and dm, gf and g
! ------------------------------------------------------------------
! input of data from a dataset, first row
IF (datasets=='singl') THEN
READ (inunit,*)xid,klimid,lbanr,gwa,wgeb,lein, &
zei,(baum(i), alt(i),mhoe(i),dm(i),gf(i),bon(i),en(i),i=1,zei)
ALLOCATE(ngroups(zei))
DO i=1,zei
IF(baum(i).EQ.8) ngroups(i)%taxid=1
IF(baum(i).EQ.10) ngroups(i)%taxid=2
IF(baum(i).EQ.11) ngroups(i)%taxid=3
IF(baum(i).EQ.15) ngroups(i)%taxid=4
if(baum(i).eq.12) ngroups(i)%taxid = 10
! Eucalyptus
IF(baum(i).EQ.30) ngroups(i)%taxid=12
IF(baum(i).EQ.31) ngroups(i)%taxid=13
IF (dm(i).eq.0) dm(i) = 0.5
IF (mhoe(i).eq.0) mhoe(i) = 1.0
IF (gf(i).eq.0) gf(i) = 0.25
ngroups(i)%locid=xid
ngroups(i)%alter=alt(i)
ngroups(i)%mhoe=mhoe(i)
ngroups(i)%gf=gf(i)
ngroups(i)%dm=dm(i)
ngroups(i)%patchsize=10000
ENDDO
CLOSE(inunit)
nlines=zei
cform=1;hlp_lai=0
! Initialisastion of stand data: area = 1ha
area=10000
area_factor=int(area/kpatchsize)
! read file head for description, write in ini-file
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
ENDIF !block for reading of input data DSW 'singl' = specially prepared for FORSKA
! read in stand dataEinlesen out of data storage for many stands
IF (datasets=='multi') THEN
select_lines=.false.
fl_num=0
ALLOCATE(ngroups(10000))
numstand= 0
nlines=1
ngroups%taxid=0
ngroups%schicht=-99
DO
READ (inunit,*,END=3333)xid,klimid,lbanr,gwa,wgeb,lein, &
zei,(baum(i),alt(i), psi(i), mhoe(i),dm(i),gf(i),bon(i),en(i),i=1,zei)
numstand = numstand +1
ngroups(nlines)%standsize= 0
DO i=1,zei
IF(baum(i).EQ.5) ngroups(nlines)%taxid=5
IF(baum(i).EQ.8) ngroups(nlines)%taxid=1
IF(baum(i).EQ.10) ngroups(nlines)%taxid=2
IF(baum(i).EQ.11) ngroups(nlines)%taxid=3
IF(baum(i).EQ.15) ngroups(nlines)%taxid=4
! the following species are preliminarily assigned
IF(baum(i).EQ.1) ngroups(nlines)%taxid=2 ! Abies alba
IF(baum(i).EQ.2) ngroups(nlines)%taxid=1 ! Acer platanoides IF(baum(i).EQ.3) ngroups(nlines)%taxid=1 ! Acer pseudoplatanus
IF(baum(i).EQ.4) ngroups(nlines)%taxid=5 ! Alnus glutinosa
IF(baum(i).EQ.6) ngroups(nlines)%taxid=1 ! Carpinus betulus
IF(baum(i).EQ.7) ngroups(nlines)%taxid=4 ! Castanea sativa
IF(baum(i).EQ.9) ngroups(nlines)%taxid=4 ! Fraxinus excelsior
IF(baum(i).EQ.12) ngroups(nlines)%taxid=5 ! Populus tremula
IF(baum(i).EQ.13) ngroups(nlines)%taxid=4 ! Quercus petraea
IF(baum(i).EQ.14) ngroups(nlines)%taxid=4 ! Quercus pubescencs
IF(baum(i).EQ.16) ngroups(nlines)%taxid=1 ! Tilia cordata
IF(baum(i).EQ.17) ngroups(nlines)%taxid=4 ! Ulmus glabra
iF(baum(i).EQ.21) ngroups(nlines)%taxid=10 ! Douglasie
iF(baum(i).EQ.22) ngroups(nlines)%taxid=6 ! Larix
iF(baum(i).EQ.23) ngroups(nlines)%taxid=7 ! Pinus strobus
iF(baum(i).EQ.24) ngroups(nlines)%taxid=10 ! Douglasie
IF (dm(i).eq.0) dm(i) = 0.5
IF (mhoe(i).eq.0) mhoe(i) = 1.0
IF (gf(i).eq.0) gf(i) = 0.25
ngroups(nlines)%locid=xid
ngroups(nlines)%alter=alt(i)
ngroups(nlines)%mhoe=mhoe(i)
ngroups(nlines)%gf=gf(i)
ngroups(nlines)%dm=dm(i)
ngroups(nlines)%patchsize=psi(i)*10000
ngroups(nlines)%standsize=psi(i)*10000
nlines=nlines+1
standsz(numstand) = standsz(numstand) + psi(i)*10000
ENDDO
ENDDO ! read loop
3333 CONTINUE
nlines=nlines-1
WRITE(*,*) nlines,'sets of data', numstand, 'sets of stands'
CLOSE(inunit)
! read in file headder for description, write into ini-file
cform=1;hlp_lai=0
! initilisation for stand data: area = stand area based on fractions of areas
area_factor=1
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,-99.)
WRITE(*,*) 'number of data lines: ', nlines
write(*,*)'number of plots for calculations: ', fl_num
ENDIF ! block for reading input data DSW, many lines = 'multi'
id=1
tmpunit=getunit()
ihelp = 1
istart=-99
DO iz=1,nlines
IF(select_lines) THEN
DO i=1,nlines_comp
IF(locid_comp(i)==ngroups(iz)%locid) GOTO 2233
ENDDO ! comparison of site id to list of sites to be selected
CYCLE
ENDIF ! end of site selection
2233 CONTINUE
WRITE(*,*) iz, nlines, ngroups(iz)%locid,ngroups(iz)%schicht
IF(datasets=='multi'.AND.(istart.NE.ngroups(iz)%locid)) THEN
WRITE(outunit,*) ngroups(iz)%locid,ngroups(iz)%standsize,'stand identifier, stand area'
ihelp = ihelp +1
istart=ngroups(iz)%locid
ENDIF
IF(datasets=='multi'.AND.ngroups(iz)%taxid==0.) THEN
WRITE(*,*) 'not the right species'
GOTO 2222
ENDIF
IF(datasets=='multi'.AND.ngroups(iz)%schicht==20) THEN
! retention trees
age=ngroups(iz)%alter
taxid=ngroups(iz)%taxid
height=ngroups(iz)%mhoe
bhd=ngroups(iz)%dm n_koh=ngroups(iz)%baumzahl
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
GOTO 2222
ENDIF ! end special treatment of retention trees
IF(datasets=='multi'.AND.ngroups(iz)%dm==0.) THEN
WRITE(4444,*)'data insufficient for: ',ngroups(iz)%locid,' line: ',iz
GOTO 2222
ENDIF
IF(datasets=='multi'.AND.ngroups(iz)%mhoe<h_sapini*0.01 .or. ngroups(iz)%gf.eq.0.) THEN
aux = ngroups(iz)%standsize/10000.
height=ngroups(iz)%mhoe
n_koh=ngroups(iz)%baumzahl* aux
age=ngroups(iz)%alter
taxid = ngroups(iz)%taxid
WRITE(4444,*)'sapling init needed for: ',ng_locid,' line: ',iz
call ini_gener_sap(outunit, taxid,age,height,n_koh)
GOTO 2222
ENDIF
optimi=.false.
anzahl= 0;start=1
allocate(helptmp(10000000))
helptmp = ' '
! generation of single trees out of population mean values
DO
helptmp = ' '
IF((start==1).or.(.not.optimi))THEN
T =7
anzahl=0
start=0
wfirst=.true.
kfirst=.true.
WRITE(*,*)ngroups(iz)%locid,ngroups(iz)%patchsize
age=ngroups(iz)%alter
dg=ngroups(iz)%dm !quadratic mean diameter
hg=ngroups(iz)%mhoe !corresponding height to dg
taxid=ngroups(iz)%taxid !species
g=ngroups(iz)%gf !basal area/ha
gpatch=g/area_factor !basal area/patch
IF (datasets=='multi') gpatch=g*ngroups(iz)%standsize/10000.
! selection of uni-height curve: Beech, Spruce, Oak calculated according to Weimann,
! other species of tree according to Kuleschis (vergl. Gerold 1990)
IF (taxid==3.OR.taxid==5) THEN
ehkwei=.false.
ELSE
ehkwei=.true.
ENDIF
IF ((dg-T).lt. 3.0) THEN
T=dg-4.0
IF (T.lt.0.3) T=0.3
ENDIF
! Estimation of Dmax out of dg (Gerold 1990)
Dmax=8.2+1.8*dg-0.01*dg**2
IF (dg.le.2) Dmax=dg+2
! calculation for the Weibull-distribution function
! in case b or c are calcuted too small, p1 and p4 respectively have to be modified
p1n=p1(taxid)
IF (p1n.lt.((1.0001-p0(taxid))/Dg)) p1n=(1.0001-p0(taxid))/Dg
b=p0(taxid)+p1n*Dg
p4n=p4(taxid)
IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax
anzit=0
thdmax=5.0
ENDIF ! end of introductory calculation and repetitions without optimisation
genDg=0 nbhd=0
saquad=0
bhdmax=0
bhdmin=100
clwdth=0
gx=0
idum=1
x=0
!----------------------------
! generation of single trees
DO
IF (gx.gt.gpatch) exit
x = ran0(idum)
bhd=b*((T/b)**c-log(1.-x))**(1./c)
if ( bhd.ge. 0.5*Dg) then
IF (bhd.gt.bhdmax) bhdmax=bhd
IF (bhd.lt.bhdmin) bhdmin=bhd
IF ((.not. optimi) .and. (bhd.gt.(1.5*dmax))) bhd=1.5*dmax
!***height calculation according to uni-height curve
IF (ehkwei) THEN
! uni-height curve of Weimann (1980)
IF (wfirst) THEN
wei_f=wei_k1(taxid)+wei_k2(taxid)*hg
wfirst=.false.
ENDIF
IF (bhd.ge.(dg-hg/2.)) THEN
height=hg+wei_f*(log(hg-dg+bhd)-log(hg))
ELSE
height=(hg+wei_f*(log(hg/2.)-log(hg))-1.3)*(bhd/(dg-hg/2.))**0.5+1.3
ENDIF
ELSE
! uni-height curve of Kuleschis (1981)
IF (kfirst) THEN
ku_a=1-(ku_a0(taxid)+ku_a1(taxid)*dg+ku_a2(taxid)*dg**2)
ku_b=ku_b0(taxid)+ku_b1(taxid)*dg+ku_b2(taxid)*dg**2
ku_c=ku_c0(taxid)+ku_c1(taxid)*dg+ku_c2(taxid)*dg**2
kfirst=.false.
ENDIF
height=hg*(ku_a+(ku_b/(bhd+dg/2.))*dg+(ku_c/(bhd+dg/2.)**2)*dg**2)
ENDIF
if(taxid.eq.10) then
! height curve after Bwinpro Douglas fir
height = 1.3 +(hg-1.3)*exp(-(0.199651*dg+4.63277655)*((1/bhd) - (1/dg)))
end if
if(taxid.eq.12.or. taxid.eq.13) then
! Medhurst et al. 1999
height = 3.665629*bhd**0.541
end if
! solution for small stands; tree dimensions below 3 m = rubbish
IF (height.gt.(bhd*3.)) height=bhd*3.
IF (height.lt.1.35) height=1.35+bhd
if(taxid.eq.12.or. taxid.eq.13) then
! Eucalyptus
hbc = crown_base_eg(height, bhd)
else
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
end if
IF ((height-hbc).lt. 0.5) hbc= height - 0.5
write(helptmp(nbhd+1), '(3f7.1,2i7)') bhd,height,hbc,age,taxid
gx=gx+1E-4*pi*(bhd/2.)**2
nbhd=nbhd+1
anzahl=anzahl+1
saquad=saquad+bhd**2
end if ! BHD test
ENDDO ! single tree calculation
!---calculates the generated Dg and test deviations of Dg and Dmax of the population value.! if deviation greater 20% a fittinf of the parameters acording to the Weibull-distribution is done
! the standard generation is repeated in several iterations.
!---the optimisation can be shut off with optimi=.false.
genDg=SQRT(saquad/nbhd)
IF((.not. optimi) .or. (Dg .lt. 7)) exit
IF(ABS(genDg-Dg).gt.(Dg/10.).or.(bhdmax-Dmax).gt. (Dmax/thdmax)) THEN
IF (ABS(genDg-Dg).gt.(Dg/10.))THEN
p1n=p1n*Dg/genDg
IF (p1n.lt.((1.0001-p0(taxid))/Dg)) p1n=(1.0001-p0(taxid))/Dg
b=p0(taxid)+p1n*Dg
ELSE
p4n=p4n*Dmax/bhdmax
IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) &
p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax
ENDIF
anzahl=anzahl-Int(nbhd)
anzit=anzit+1
IF (anzit.ge.50) THEN
IF (thdmax.eq.2) THEN
print *,'id/zei: ',id,iz,' Optimization not successful. Biased STAND.INI will be generated'
optimi=.false.
ELSE
anzit=0
thdmax=2.0
b=p0(taxid)+p1(taxid)*Dg
c=p2(taxid)+p3(taxid)*Dg+p4(taxid)*Dmax
ENDIF
ENDIF
ELSE
exit
ENDIF
ENDDO
! end of generation of single trees
! classification of single values in diameter cohorts
clwdth=1+AINT((bhdmax-bhdmin)/ncl1) !calculation of class widths
! write(4444,*) 'clwdth', clwdth, bhdmax, bhdmin, ncl1
DO i=1,ncl1
nz(i)=0
zbhd(i)=0
zheigh(i)=0
zhbc(i)=0
ENDDO
DO j=1,nbhd
read(helptmp(j), *) bhd,height,hbc,age,taxid
IF(height<1.3) WRITE(4444,*)'bhd ',bhd,'height ',height,'art ',taxid
icl=INT(bhd/clwdth)+1
IF(icl.gt.ncl1) icl=ncl1
nz(icl)=nz(icl)+1 !addition stem numbre of diameter classes
zbhd(icl)=zbhd(icl)+bhd !sum of diametes of diameter calsses
zheigh(icl)=zheigh(icl)+height !sum of height value of classes
zhbc(icl)=zhbc(icl)+hbc !sum of crown starting height of classes
ENDDO
deallocate(helptmp)
tot_crown_area=0.
DO i=1,ncl1
IF (nz(i).ne.0) THEN
bhd=zbhd(i)/nz(i)
height=zheigh(i)/nz(i)
hbc=zhbc(i)/nz(i)
n_koh=NINT(nz(i)/area_factor)
tot_crown_area=tot_crown_area+n_koh*PI*(MIN(spar(taxid)%crown_a*bhd+spar(taxid)%crown_b,spar(taxid)%crown_c))**2
ENDIF
ENDDO
IF(tot_crown_area>1.1*kpatchsize) THEN
corr_la=kpatchsize/tot_crown_area
ELSE
corr_la=1.
ENDIF
DO i=1,ncl1
IF (nz(i).ne.0) THEN
bhd=zbhd(i)/nz(i)
height=zheigh(i)/nz(i)
hbc=zhbc(i)/nz(i)
n_koh=NINT(nz(i)/area_factor)
! --- 4C-specific calculations:
IF(height<1.3) WRITE(4444,*)ngroups(iz)%locid,'bhd ',bhd,'height ',height,'art ',taxid
IF(height*100<h_sapini) THEN
CALL sapini(outunit,taxid, height,hbc, n_koh,age)
WRITE(4444,*)ngroups(iz)%locid,bhd,taxid
ELSE
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
ENDIF
ENDDO !cohort loop
2222 CONTINUE
if(datasets=='multi') then
IF (iz.ne.nlines.AND.datasets=='multi'.AND.(istart.NE.ngroups(iz+1)%locid)) WRITE(outunit,*) '-99.9'
end if
2244 CONTINUE
ENDDO !line loop
CLOSE(outunit)
CLOSE(ctrlunit)
RETURN
! ****************************************************************************
! case(6) stand generation if data source is from Bavarian inventories
CASE(6)
print *, ' Forest stand data set: Bavarian inventories'
infile='/data/safe/4C/4C_input/stand/Bayernw.dat'
source='B'
OPEN (inunit, FILE=TRIM(infile), STATUS='old')
listunit=GETUNIT()
OPEN (listunit, FILE='/home/lasch/4c/v0.99e1/input/koord.txt', STATUS='old')
! ------------------------------------------------------------------
! generated standard values of data from data storage based on
! estimation routines of Nagel and Biging (1995), Nagel (1995) and
! Gerold (1990).
! ------------------------------------------------------------------
!
! The following variables are read from forest inventory data:
! Species(baum),Age(alt),Quadratic Mean Diameter(dm),Height of tree with dm(mhoe),
! Basal area(gf),Yield Class(bon),"Ertragsniveau"(en)
!
! ------------------------------------------------------------------
! read in stad data of multiple stands out of records
select_lines=.true.
datasets='multi'
fl_num=0
IF(select_lines) THEN
READ(listunit,*)nlines_comp
ALLOCATE(locid_comp(nlines_comp))
DO i=1,nlines_comp ! reading list of sites to be initialised
READ(listunit,*) locid_comp(i)
ENDDO ! end reading list of sites to be initialised
ENDIF ! end of reading file with sites to be selected
IF(select_lines) CLOSE(listunit)
CALL assign_BAY
CALL init_plenter_param
READ (inunit,*) READ (inunit,*)nlines
ALLOCATE(ngroups(nlines))
istart=1
READ(inunit,*) dummy, dummy, dummy, ngroups(1)%locid, dummy, &
ngroups(1)%schicht, ngroups(1)%BRAid, dummy, dummy, ngroups(1)%alter, &
dummy, dummy, ngroups(1)%dm, ngroups(1)%mhoe, ngroups(1)%baumzahl, &
ngroups(1)%gf, ngroups(1)%volume, dummy
ngroups(1)%taxid=tax_of_BRA_id(ngroups(1)%BRAid)
ngroups(1)%standsize=40000
IF(ngroups(1)%alter==0.OR.ngroups(1)%mhoe==0.OR.ngroups(1)%dm==0.OR.ngroups(1)%volume==0.OR.ngroups(1)%gf==0) CALL data_gap_fill_DSW(1)
DO i=2,nlines
READ(inunit,*) dummy, dummy, dummy, ngroups(i)%locid, dummy, &
ngroups(i)%schicht, ngroups(i)%BRAid, dummy, dummy, ngroups(i)%alter, &
dummy, dummy, ngroups(i)%dm, ngroups(i)%mhoe, ngroups(i)%baumzahl, &
ngroups(i)%gf, ngroups(i)%volume, dummy
WRITE(*,*) 'set no', i, 'read'
ngroups(i)%taxid=tax_of_BRA_id(ngroups(i)%BRAid)
ngroups(i)%standsize=40000
! preliminary solution: larches mapped to pine
IF(ngroups(i)%taxid==6) ngroups(i)%taxid=3
IF(ngroups(i)%taxid==0) THEN
ELSE
IF(ngroups(i)%alter==0.OR.ngroups(i)%mhoe==0.OR.ngroups(i)%dm==0.OR.ngroups(i)%gf==0) THEN
WRITE(7333,*)'set ',i,'not enough data or below 1.3 m height'
! CALL data_gap_fill_DSW(i)
ENDIF
ENDIF
IF(ngroups(i)%locid.NE.ngroups(istart)%locid) THEN
istart=i
fl_num=fl_num+1
ENDIF
ENDDO ! readin loop for multi data-set
CLOSE(inunit)
! read file headder for description, write in ini-file
cform=1;hlp_lai=0
! initialisation of stand records: area =
! stand area calculated according to partial areas.
area_factor=1
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,-99.)
id=1
WRITE (fnam2,'(a,i1,a)') 'schicht',id,'.tmp'
tmpunit=getunit()
istart=-99
DO iz=1,nlines
ng_locid = ngroups(iz)%locid
taxid=ngroups(iz)%taxid
IF(select_lines) THEN
DO i=1,nlines_comp
IF(locid_comp(i)==ng_locid) GOTO 2255
ENDDO ! comparison of site id to list of sites to be selected
CYCLE
ENDIF ! end of site selection
2255 CONTINUE
IF(datasets=='multi'.AND.(istart.NE.ng_locid)) THEN
WRITE(outunit,*) ng_locid,ngroups(iz)%standsize,'stand identifier, stand area'
istart=ng_locid
aux = ngroups(iz)%standsize/10000.
ENDIF
IF(datasets=='multi'.AND.taxid==0.) THEN
! solution for bushes must be found
WRITE(*,*) 'not the right species'
GOTO 2277
ENDIF
IF(ngroups(iz)%baumzahl<30.AND.ngroups(iz)%baumzahl>0) ngroups(iz)%schicht=5
IF(datasets=='multi'.AND.ngroups(iz)%schicht==5) THEN ! retention trees can be directly initialized since they are not distributed onto different height cohorts
WRITE(4444,*) 'single type ',ngroups(iz)%schicht
age=ngroups(iz)%alter
height=ngroups(iz)%mhoe
bhd=ngroups(iz)%dm
n_koh=ngroups(iz)%baumzahl*aux
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
GOTO 2277
ENDIF ! end special treatment of retention trees
IF(datasets=='multi'.AND.ngroups(iz)%dm==0.and.ngroups(iz)%mhoe>h_sapini*0.01) THEN
WRITE(4444,*)'data insufficient for: ',ng_locid,' line: ',iz
GOTO 2277
ENDIF
IF(datasets=='multi'.AND.ngroups(iz)%mhoe<h_sapini*0.01) THEN
height=ngroups(iz)%mhoe
n_koh=ngroups(iz)%baumzahl* aux
age=ngroups(iz)%alter
call ini_gener_sap(outunit, taxid,age,height,n_koh)
GOTO 2277
ENDIF
T=7
age=ngroups(iz)%alter
dg=ngroups(iz)%dm !quadratic mean diameter
hg=ngroups(iz)%mhoe !corresponding height to dg
g=ngroups(iz)%gf !basal area/ha
gpatch=g*4. !basal area/patch
bz=ngroups(iz)%baumzahl*4. !tree numbre/patch
! clwdth=dg/20.
clwdth=dg/5
! selection of uni-height curve: beech, spruce, oak calculation according to Weimann,
! other species of trees after Kuleschis (vergl. Gerold 1990)
IF (taxid==3.OR.taxid==5) THEN
ehkwei=.false.
ELSE
ehkwei=.true.
ENDIF
! zuweisen der PArameterwerte fr Einheitshhenkurve
IF (ehkwei) THEN
! uni-height curve from Weimann (1980)
wei_f=wei_k1(taxid)+wei_k2(taxid)*hg
ELSE
! uni-height curve from Kuleschis (1981)
ku_a=1-(ku_a0(taxid)+ku_a1(taxid)*dg+ku_a2(taxid)*dg**2)
ku_b=ku_b0(taxid)+ku_b1(taxid)*dg+ku_b2(taxid)*dg**2
ku_c=ku_c0(taxid)+ku_c1(taxid)*dg+ku_c2(taxid)*dg**2
ENDIF
IF ((dg-T).lt. 3.0) THEN
T=dg-4.0
IF (T.lt.0.3) T=0.3
ENDIF
! Estimation of Dmax from dg (Gerold 1990)
Dmax=8.2+1.8*dg-0.01*dg**2
IF (dg.le.2) Dmax=dg+2
! Calculation of parameter for Weibull-distribution
! in case b or c is calculated too small,
! p1 and p4 respectively have to be modified
p1n=p1(taxid)
IF (p1n.lt.((1.0001-p0(taxid))/Dg)) p1n=(1.0001-p0(taxid))/Dg
b=p0(taxid)+p1n*Dg
Dmin = 0.1*Dg
IF(Dg>70) Dmin = 2.*Dg - Dmax
p4n=p4(taxid)
IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax
anzit=0 thdmax=5.0
helpz=0
DO
imax=INT((Dmax-Dmin)/clwdth)
if(imax.gt.30) then
imax= 30
clwdth= (Dmax-Dmin)/30.
end if
if(helpz.gt.50) goto 2277
helpz= helpz + 1
Fint(0)=0.
gx=0.
bhd=Dmin+0.5*clwdth
DO i = 1,imax
Fint(i)=1-exp(-((bhd-Dmin)/b)**c)
gx=gx+(Fint(i)-Fint(i-1))*bhd**2
bhd=bhd+clwdth
END DO
gx=gx*PI/4*1e-4*bz
IF(ABS(gx-gpatch)>0.02*gpatch) THEN
IF(gx>gpatch) THEN
c=c*gpatch/gx
ELSE
IF(Dmin<0.8*Dg) THEN
Dmin=Dmin*1.05
ELSE
c=c*gx/gpatch
ENDIF
ENDIF
ELSE
EXIT
ENDIF
END DO
bhd=Dmin+0.5*clwdth
DO i = 1,imax
n_koh=NINT((Fint(i)-Fint(i-1))*bz)
!***calculate height according to uni-height curve
IF (ehkwei) THEN
! uni-height curve from Weimann (1980)
IF (bhd.ge.(dg-hg/2.)) THEN
height=hg+wei_f*(log(hg-dg+bhd)-log(hg))
ELSE
height=(hg+wei_f*(log(hg/2.)-log(hg))-1.3)*(bhd/(dg-hg/2.))**0.5+1.3
ENDIF
ELSE
! uni-height curve from Kuleschis (1981)
height=hg*(ku_a+(ku_b/(bhd+dg/2.))*dg+(ku_c/(bhd+dg/2.)**2)*dg**2)
ENDIF
! solution for small stands; tree dimensions below 3 m = rubbish
IF (height.gt.(bhd*3.)) height=bhd*3.
IF (height.lt.1.35) height=1.35+bhd
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
IF ((height-hbc).lt. 0.5) hbc= height - 0.5
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
if(fail.eq.1) write(4444,*) 'negative root in newton', ng_locid,iz
bhd=bhd+clwdth
END DO
2277 CONTINUE
IF (iz.ne.nlines.AND. datasets=='multi'.AND.(istart.NE.ngroups(iz+1)%locid)) WRITE(outunit,*) '-99.9'
2266 CONTINUE
ENDDO !sign loop
CLOSE(outunit)
CLOSE(ctrlunit)
RETURN
CASE(2)
334 CONTINUE CALL assign_DSW
inwahl=0
source='S'
PRINT *, 'If you want to use SILVA data, type: 1'
PRINT *, 'If you want to use levelII data from Sachsen, type: 2'
PRINT *, 'If you want to use single tree data with tree class information, type: 3'
PRINT *, ' if you want to use data like level II single tree data and generate one tree cohorts, type: 4'
READ(*,*) inwahl
IF (inwahl<1.OR.inwahl>4) THEN
WRITE(*,*) 'You should use integer 1, 2,3 or 4 for the choice of data source'
GOTO 334
ENDIF
333 CONTINUE
IF (inwahl==1) PRINT *, ' Forest stand data set: SILVA (classification must be performed)'
IF (inwahl==2) PRINT *, ' Forest stand data set: levelII Sachsen (classification must be performed)'
IF (inwahl==3) PRINT *, ' Forest stand data set: single tree data with tree type information (classification must be performed)'
IF (inwahl==4) PRINT *, ' Forest stand data set: single tree data without clissification'
WRITE(*,'(A)')
WRITE(*,'(A)')' Do you want to read the input file from '
WRITE(*,'(A)')' 1 - the Standard 4C stand directory on data/safe/4C/4C_input/stand'
WRITE(*,'(A)')' 2 - or do you want to specify another directory?'
WRITE(*,'(A)',advance='no') ' ***Make your choice: '
READ(*,*) dir_flag
IF(dir_flag.EQ.1) THEN
WRITE(*,'(A)',advance='no')' Input file: '
READ (*,'(A)') infile
ELSEIF(dir_flag.EQ.2) THEN
WRITE(*,'(A)',advance='no')' Input directory and file: '
READ (*,'(A)') infile
ELSE
WRITE(*,*) 'You should use integer 1 or 2 for the choice of the input mode. Please try again!'
GOTO 333
ENDIF
337 CONTINUE
cform=1;hlp_lai=0
IF(dir_flag.EQ.1) OPEN (inunit,FILE='/data/safe/4C/4C_input/stand/'//trim(infile),STATUS='old')
IF(dir_flag.EQ.2) OPEN (inunit,FILE=trim(infile),STATUS='old')
! initialising for stand records: area = 1ha
area=10000
IF(inwahl==2.OR.inwahl==3.OR.inwahl==4) THEN
! class width
clwdth=1 !set diameter of classes width
READ(inunit,'(a85)')zeile
READ(inunit,*) area
READ(inunit,'(a85)')zeile
ENDIF
area_factor = 1.
kpatchsize = area
! read in file headder for descriptions, write in ini-file
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
! classification of single values into diameter cohorts
IF(inwahl==1) THEN
READ(inunit,'(a85)')zeile
READ(inunit,'(a85)')zeile
ENDIF
335 CONTINUE
DO i=1,ncl1
nz(i)=0
zbhd(i)=0
zheigh(i)=0
zhbc(i)=0
ENDDO
nhelp=0
DO
IF(inwahl==1) READ(inunit,*,IOSTAT=ios)xnr,baumid,bhd,height,hbc,kd,xxr,xyr,xxi,xyi
IF(inwahl==2.or.inwahl.eq.4) THEN READ(inunit,*,IOSTAT=ios)xnr,taxid,bhd,height,hbc,age
nhelp = nhelp+1
if(bhd.le.10) bhd=11.
bhd=bhd/10.
IF(hbc>-99.99.AND.hbc<-99.8) THEN
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
IF(height-hbc<0.5) CALL error_mess(time,"crown to shallow in tree",REAL(xnr))
ENDIF
ENDIF
IF(inwahl==3) THEN
READ(inunit,*,IOSTAT=ios)xnr,taxid,bhd,height,hbc,ager,status
IF(taxid>=100) taxid=tax_of_BRA_id(taxid)
age = INT(ager)
bhd=bhd/10.
IF(hbc>-99.99.AND.hbc<-99.8) THEN
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
IF(height-hbc<0.5) CALL error_mess(time,"crown to shallow in tree",REAL(xnr))
IF((height-hbc)/height<0.5) hbc=0.5*height
IF(bhd<=3.) hbc=0.
ENDIF
ENDIF
IF (ios<0) exit
IF (xnr==-9999) exit
IF (inwahl==4) exit
icl=INT(bhd/clwdth)+1
IF(inwahl.eq.4.or.(inwahl==3.AND.status.NE.'F'.AND.status.NE.'Z'.AND.status.NE.'V'.and.status.NE.'H'.and.status.NE.'U'.and. status.NE.'B'))THEN
ELSE
IF(icl.gt.ncl1) icl=ncl1
nz(icl)=nz(icl)+1 !sum stem numbre of diameter class
zbhd(icl)=zbhd(icl)+bhd !sum up the diameters of a class
zheigh(icl)=zheigh(icl)+height !sum up height value of a class
zhbc(icl)=zhbc(icl)+hbc !sum up crown startin height of a class
ENDIF
ENDDO
nzsum=sum(nz)
IF(inwahl.ne.4) THEN
tot_crown_area=0.
DO i=1,ncl1
IF (nz(i).ne.0) THEN
bhd=zbhd(i)/nz(i)
height=zheigh(i)/nz(i)
hbc=zhbc(i)/nz(i)
if(hbc<0.025) hbc = 0.
if(hbc>=0.025.and.hbc<0.05) hbc =0.05
n_koh=NINT(nz(i)/area_factor)
IF(inwahl==1) THEN
SELECT CASE(baumid)
CASE(5)
taxid=1
CASE(1)
taxid=2
CASE(3)
taxid=3
CASE default
taxid=99
END select
ENDIF
tot_crown_area=tot_crown_area+n_koh*PI*(MIN(spar(taxid)%crown_a*bhd+spar(taxid)%crown_b,spar(taxid)%crown_c))**2
ENDIF
ENDDO
IF(tot_crown_area>1.1*kpatchsize) THEN
corr_la=kpatchsize/tot_crown_area
ELSE
corr_la=1.
ENDIF
IF(pass==1) THEN
mixed_tot_ca = mixed_tot_ca + tot_crown_area
ELSE
corr_la=kpatchsize/mixed_tot_ca ENDIF
DO i=1,ncl1
IF (nz(i).ne.0) THEN
bhd=zbhd(i)/nz(i)
height=zheigh(i)/nz(i)
hbc=zhbc(i)/nz(i)
if(hbc<0.025) hbc = 0.
if(hbc>=0.025.and.hbc<0.05) hbc =0.05
n_koh=NINT(nz(i)/area_factor)
IF(inwahl==1) THEN
SELECT CASE(baumid)
CASE(5)
taxid=1
CASE(1)
taxid=2
CASE(3)
taxid=3
CASE default
taxid=99
END select
ENDIF
! --- 4C-specific calculation:
WRITE(*,*) 'call :', taxid,bhd,height,hbc,nz(i),n_koh
IF( height<(h_sapini/100.)) then
call sapini(outunit,taxid, height, hbc, n_koh,age)
ELSE
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
ENDIF
ENDDO
else if(xnr.ne.-999) then
n_koh = 1
print*, 'xnr:', xnr
IF( height<(h_sapini/100.)) then
call sapini(outunit,taxid, height, hbc, n_koh,age)
ELSE
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
end if
IF (xnr==-9999) GOTO 335
if(inwahl==4.and.xnr==-999) then
CLOSE(inunit)
CLOSE(outunit)
CLOSE(ctrlunit)
RETURN
end if
if(inwahl==4) goto 335
CLOSE(inunit)
CLOSE(outunit)
IF(mixed_tot_ca>1.1*kpatchsize .AND. pass == 1) THEN
OPEN (outunit,FILE=TRIM(treefile(ip)),STATUS='replace')
pass = 2
GOTO 337
ENDIF
CLOSE(ctrlunit)
RETURN
CASE(3)
444 print *, ' Forest stand data set: Level2-Daten'
source='L'
WRITE(*,'(A)')
WRITE(*,'(A)')' Do you want to read the input file from '
WRITE(*,'(A)')' 1 - the Standard 4C stand directory on data/safe/4C/4C_input/stand'
WRITE(*,'(A)')' 2 - or do you want to specify another directory?'
WRITE(*,'(A)',advance='no') ' ***Make your choice: '
READ(*,*) dir_flag
IF(dir_flag.EQ.1) THEN
WRITE(*,'(A)',advance='no')' Input file: ' READ (*,'(A)') infile
ELSEIF(dir_flag.EQ.2) THEN
WRITE(*,'(A)',advance='no')' Input directory and file: '
READ (*,'(A)') infile
ELSE
WRITE(*,*) 'You should use integer 1 or 2 for the choice of the input mode. Please try again!'
GOTO 444
ENDIF
cform=1;hlp_lai=0
IF(dir_flag.EQ.1) OPEN (inunit,FILE='/data/safe/4C/4C_input/stand/'//trim(infile),STATUS='old')
IF(dir_flag.EQ.2) OPEN (inunit,FILE=trim(infile),STATUS='old')
!------------------------------------------------------------------
! Read in level II data according to diamter classes
READ(inunit,'(a85)')zeile
READ(inunit,'(a85)')zeile
READ(inunit,'(a85)')zeile
READ(inunit,*)age,taxid,area, rsap, &
dclmin, & !smallest diameter of experimentation patches
ndcl, & !amount diameter class
dcwdth !class width
READ(inunit,*)h_para,h_parb, & !parameter of height function after Lockow
(n_dc(i),i=1,ndcl) !stem numbre per diameter class
close(inunit)
clwdth=dcwdth
! ---------------------------------------------------------------------
! current patch size = value specified by kpatchsize
area_factor=int(area/kpatchsize)
! read in file headder for desciption, write into ini-file
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
DO i=1,ncl1
nz(i)=0
zbhd(i)=0
zheigh(i)=0
zhbc(i)=0
ENDDO
bhdcl=dclmin
DO i=1,ndcl
bhd=bhdcl
height=h_para*(0.01*bhd)**h_parb !height function after regression from Lockow
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
IF ((height-hbc).lt. 0.5) hbc= height - 0.5
icl=INT(bhd/clwdth)+1
IF(icl.gt.ncl1) icl=ncl1
nz(icl)=nz(icl)+n_dc(i) !sum stem numbre of diameter class
zbhd(icl)=zbhd(icl)+bhd*n_dc(i) !sum up diameters of a class
zheigh(icl)=zheigh(icl)+height*n_dc(i) !sum up height values of a class
zhbc(icl)=zhbc(icl)+hbc*n_dc(i) !sum up crown starting height of a class
bhdcl=bhdcl+dcwdth
ENDDO
smaldc(1)=.false.
DO i=1,ncl1
IF (smaldc(i)) THEN
IF (i<ncl1) smaldc(i+1)=.true.
ELSE
IF (i<ncl1) smaldc(i+1)=.false.
n_koh=NINT(nz(i)/area_factor)
IF (n_koh>0) THEN
IF (i<ncl1) smaldc(i+1)=.true.
ENDIF
ENDIF
ENDDO
bigdc(ncl1)=.false.
DO i=ncl1,1,-1
IF (bigdc(i)) THEN IF (i>1) bigdc(i-1)=.true.
ELSE
IF (i>1) bigdc(i-1)=.false.
n_koh=NINT(nz(i)/area_factor)
IF (n_koh>0) THEN
IF (i>1) bigdc(i-1)=.true.
ENDIF
ENDIF
ENDDO
DO i=1,ncl1
IF (nz(i).ne.0) THEN
n_koh=NINT(nz(i)/area_factor)
IF (n_koh==0) THEN !if no trees in cohorte, shift trees to next class
zbhd(i+1)=zbhd(i+1)+zbhd(i) !add diameter to sum of next class
zheigh(i+1)=zheigh(i+1)+zheigh(i) !add height to sum of next class
zhbc(i+1)=zhbc(i+1)+zhbc(i) !add hbc to sum of next class
nz(i+1)=nz(i+1)+nz(i) !add trees to next class
nz(i)=0 !empty class
ELSE
bhd=zbhd(i)/nz(i)
height=zheigh(i)/nz(i)
hbc=zhbc(i)/nz(i)
! --- 4C-specific calculations:
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
ENDIF
IF (.not.bigdc(i+1)) exit
ENDDO
DO j=ncl1,(i+1),-1
IF (nz(j).ne.0) THEN
n_koh=NINT(nz(j)/area_factor)
IF (n_koh==0) THEN !if no trees in cohorte, shift trees to next class
zbhd(j-1)=zbhd(j-1)+zbhd(j) !add diameter to sum of next class
zheigh(j-1)=zheigh(j-1)+zheigh(j) !add height to sum of next class
zhbc(j-1)=zhbc(j-1)+zhbc(j) !add hbc to sum of next class
nz(j-1)=nz(j-1)+nz(j) !add trees to next class
nz(j)=0 !empty class
ELSE
bhd=zbhd(j)/nz(j)
height=zheigh(j)/nz(j)
hbc=zhbc(j)/nz(j)
! --- 4C-specific calculation:
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
ENDIF
IF (.not. smaldc(i)) exit
ENDDO
CLOSE(outunit)
CLOSE(ctrlunit)
RETURN
CASE(4)
WRITE(*,*) 'Do you want to use the standard procedure - type: S'
WRITE(*,*) 'or Manfred Lexers input format - type: L'
READ(*,*) source
WRITE(*,'(A)',advance='no')' Input file: '
READ(*,'(A)') infile
cform=1;hlp_lai=0
IF(flag_volfunc.EQ.0) THEN
WRITE(*,'(A)',advance='no')' Input form factor (Default in 4C = 1): '
READ *, cform
ENDIF
OPEN (inunit,FILE=TRIM(infile),STATUS='old')
! read in data from input-file
IF (source=='S') THEN
READ(inunit,*)source, taxid, rsap READ(inunit,*) area
READ(inunit,*,END=10)n,k,age
area_factor = 1.
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
!read in data
DO i=1,k
READ(inunit,*,END=10)bhd,height,share,hbc
IF(hbc>-99.99.AND.hbc<-99.8) THEN
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
END IF
n_koh = NINT(share*n)
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDDO
ELSE
READ(inunit,*) area
kpatchsize= area
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
!read in data
DO
READ(inunit,*,iostat=ios)bhd,taxid,height,n_koh,age
if(ios < 0) exit
IF(height.ne.0 .AND. n_koh.ne.0) then
IF(height<h_sapini*0.01) then
CALL ini_gener_sap(outunit,taxid,age,height,n_koh)
else
hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
end if
ENDIF
ENDDO
ENDIF
10 continue
PRINT*, 'Bestandesblattflche (pro ha): ', hlp_lai*area_factor
CLOSE(inunit)
CLOSE(outunit)
CLOSE(ctrlunit)
! FORGRA data input
CASE(5)
WRITE(*,'(A)',advance='no')' Input file: '
READ(*,'(A)') infile
cform=1;hlp_lai=0
IF(flag_volfunc.EQ.0) THEN
WRITE(*,'(A)',advance='no')' Input form factor (Default in 4C = 1): '
READ *, cform
ENDIF
OPEN (inunit,FILE=TRIM(infile),STATUS='old')
! read in data from input file
READ(inunit,*)source, rsap
READ(inunit,*) area
READ(inunit,*,END=20)n,k
area_factor=int(area/kpatchsize)
CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
!read in data
DO i=1,k
READ(inunit,*,END=20)bhd,height,share,hbc,age,taxid
n_koh=NINT(share*n/area_factor)
IF(height<h_sapini) THEN
CALL sapini(outunit,taxid, height,hbc, n_koh,age)
ELSE
CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
ENDIF
ENDDO
20 CONTINUE
CLOSE(outunit)
CLOSE(ctrlunit)
CASE default
PRINT *,' False number'
RETURN
END select
WRITE(*,*) 'initialisation terminated'
deallocate (zheigh, zbhd, zhbc, nz)
deallocate (smaldc, bigdc)
if (allocated(locid_comp))deallocate(locid_comp)
END subroutine initia
!****************************!
!* SUBROUTINE TREEINI *!
!****************************!
SUBROUTINE treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
! Species (taxid) must be handed over (Beech 1, Spruce 2, Pine 3, Oak 4)
! Source is specifying data source
! height and hbc are read in meter and is converted later to cm
! n_koh numbre of trees in a cohort
! -------------------------------------------------------------------------
USE data_init
USE data_par
USE data_simul
USE data_species
USE data_stand
USE data_help
IMPLICIT none
! ----VARIABLEN---
REAL :: bhd,height,hbc,hlp_lai,hfd,vd,VS,Vg,k1,k2,k3,hm,Ahc,Veff,dbc,corr_la
REAL :: swheight,stembio,afol,asap,dbase, dcbase,volratio,d1,d2,h1,h2,a1,b0, x_ges
INTEGER :: taxid, & ! species number
age, & ! tree age
n_koh
INTEGER :: outunit,ctrlunit !units
CHARACTER*85 zeile
CHARACTER(75):: infile
CHARACTER :: source
REAL rsap, cform, sicrsap, lifrac, rsapfit
INTEGER taumax, ring
! function
REAL newton
sicrsap=rsap
! since the fraction of wood which is sapwood generally is not measured at the
! plots for which the model is initialized, it needs to be approximated
! the following rsap initialisation has been fitted to a pine run at Kienhorst
rsapfit=1.-1.544e-8*age**4+4.343e-6*age**3-3.359e-4*age**2-4.557e-4*age
! estimation of rsap from average diameter increase
! attention: age of tree when first ring has been grown at 1.3 m must be estimated
! for the time being this is set to 5
! If hbc < h_breast, rsap and Asap (below) have to be calculated at lower height
hm=height
height=height*100
hbc=hbc*100
lifrac=1.-spar(taxid)%pss
IF(age>6) THEN
IF(hbc<h_breast) THEN
taumax=age-INT(hbc/h_breast*5.)
ELSE
taumax=age-5
ENDIF rsap=0.
DO ring = 0,taumax-1
rsap=rsap+exp(ring*log(lifrac))*(2.*(taumax-ring)-1.)
END DO
rsap=rsap/taumax**2
ELSE
rsap=1.
ENDIF
rsap=rsap*corr_la
! --- calculate height of Sapwood-Pipes and stem-mass
swheight=2.*hbc/3.+height/3.
if(taxid.ne.12. .and. taxid.ne.13) then
if(taxid.eq.10) then
! after BWINpro , Bergel 1974
hfd = (-200.31914/(height*bhd*bhd))+(0.8734/bhd) - 0.0052*log(bhd*bhd) + 7.3594/(height*bhd) + 0.46155
else
k1=par_S(taxid,1)+par_S(taxid,2)*log(bhd)+par_S(taxid,3)*log(bhd)**2
k2=par_S(taxid,4)+par_S(taxid,5)*log(bhd)+par_S(taxid,6)*log(bhd)**2
k3=par_S(taxid,7)+par_S(taxid,8)*log(bhd)+par_S(taxid,9)*log(bhd)**2
hfd=exp(k1+k2*log(hm)+k3*log(hm)**2)
end if
! vd volume with SILVA equations
vd=(hfd*pi*bhd**2)/40000
else
! Eucalyptus, Binkley et al 2002
vd = 0.00005447*bhd**1.921157*(height/100)**0.950581
! Stape et. al 2010 Fkt. VER
vd = (0.027*bhd**2.221*(height/100)**0.625)/500
! Stape et al 2010 Fkt ARA
vd = (0.004*bhd**1.959*(height/100)**1.512)/500
end if
! vs volume with Eberswalde equations
if(taxid.eq.3) vs = exp(parEBW(10,1)+parEBW(10,2)*log(bhd)+parEBW(10,3)*log(hm))
IF(taxid==3) vd = vs
IF(flag_volfunc.EQ.0) THEN
IF(source.ne.'S') stembio= swheight*spar(taxid)%prhos*cform*pi*(bhd/2.)**2
IF(source.eq.'S') THEN
stembio=vd*spar(taxid)%prhos*1000000
bhd= SQRT(stembio*4/(swheight*spar(taxid)%prhos*cform*pi))
ENDIF
! --- seperation of sap wood and heartwood and sap wood cross section
x_Ahb= 0.
x_sap=rsap*stembio
x_hrt=(1-rsap)*stembio
asap=rsap*pi*(bhd/2.)**2
! --- estimation of leafe matter and leave area
x_fol=asap*spar(taxid)%pnus
afol=x_fol*(spar(taxid)%psla_min+0.5*spar(taxid)%psla_a)
hlp_lai=hlp_lai+afol*n_koh
! --- fine root matter roughly estimated
x_frt=x_fol
IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,height,hbc,age,n_koh,taxid
ELSEIF(flag_volfunc.EQ.1) THEN
IF (hbc>h_breast.AND.hbc<h_breast+h_bo_br_diff) hbc=h_breast
IF (hbc==h_breast) dbc=bhd
IF (hbc<h_breast) THEN
dbc=bhd/height*(h_breast-hbc)+bhd ! dbc = diameter at base of the crown
asap=PI/4.*dbc**2.*rsap
ELSE
asap=PI/4.*bhd**2.*rsap !change Martin bhd>>dbc as written ins description and rsap weg
ENDIF
rsap = asap/((pi*bhd*bhd)/4)
x_sap=spar(taxid)%prhos*asap*swheight
! first guess for start values of Ahc
IF (hbc<=h_breast) THEN
Ahc=PI/4.*dbc**2.-asap
x_Ahb=PI/4.*(dbc*age/taumax)**2.-asap ELSE
Ahc=PI/4.*bhd**2.*(1.-rsap)*0.04
Ahc=Newton(Ahc,asap,bhd,hbc,height,Vd)
if(fail.eq.1) return
x_Ahb=PI/4.*((bhd-(4./PI*(asap+Ahc))**0.5*h_breast/hbc)/(1.-h_breast/hbc))**2-asap
ENDIF
! Vg for test purposes = volume if no heartwood in crown space
Vg=1./3.*height*asap+2./3.*hbc*asap+1./3.*hbc*x_Ahb
! --- seperation of sap wood and heartwood and splitting of sap wood cross section
stembio=spar(taxid)%prhos*(1./3.*height*(asap+Ahc)+1./3.*hbc*(2.*asap+x_ahb+(x_ahb*ahc)**0.5))
volratio=1.0
if(infile=='input/bwi2_blmwert1.prn') then
!Spruce
if(taxid.eq.2)then
!after Wirth et al. 2002 Tree physiology
b0=-2.83958
d1=2.55203
d2=-0.14991
h1=-0.19172
h2=0.25739
a1=-0.08278
volratio=(exp(b0+d1*log(bhd)+d2*(log(bhd))**2+h1*log(height/100)+h2*(log(height/100))**2+a1*log(age+0.01)))/stembio
endif
!Pine
if(taxid.eq.3)then
!after Zianis et al. 2005 Silva Fennica EFI BEFs Europe
volratio=exp(-2.6768+7.5939*(bhd/(bhd+13))+0.0151*height/100+0.8799*log(height/100))/stembio
endif
!for douglas fir (correction after bartelink 1996, forest ecol. manag.)
if(taxid.eq.10)then
volratio=exp(-3.229+1.901*log(bhd)+0.807*log(height/100))/stembio
endif
end if
x_sap=x_sap*volratio
x_hrt=stembio*volratio-x_sap
x_ges=x_hrt+x_sap
x_Ahb=x_Ahb*volratio
asap=asap*volratio
if (x_hrt/x_ges .gt. 0.5 .and. taxid .eq. 2 .and. age .gt. 100) then !query too heigh heart wood percentage
x_hrt=0.5*stembio*volratio
x_sap=0.5*stembio*volratio
endif
if (x_hrt/x_ges .gt. 0.35 .and. taxid .eq. 3 .or. taxid .eq. 10) then !query too heigh heart wood percentage
x_hrt=0.35*stembio*volratio
x_sap=0.65*stembio*volratio
endif
Veff=(1./3.*height*(asap+Ahc)+1./3.*hbc*(2.*asap+x_ahb+(x_ahb*ahc)**0.5))*0.000001
dbase = ((x_Ahb+asap)*4./PI)**0.5
dcbase = ((Ahc+asap)*4./PI)**0.5
WRITE(ctrlunit,'(2I5, 12F12.5)') n_koh,age,height,hbc,bhd,rsap,dbase,dcbase,asap,ahc,x_ahb,Vg/1000000,Vd,Veff
! --- estimation leaf matter and leaf area
x_fol=asap*spar(taxid)%pnus*volratio
afol=x_fol*(spar(taxid)%psla_min+0.5*spar(taxid)%psla_a)
hlp_lai=hlp_lai+afol*n_koh
! --- fine root matter rough estimate
x_frt=x_fol
IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7, 2f12.5)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,height,hbc,age,n_koh,taxid, dcbase,bhd
ENDIF
END subroutine treeini
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! SUBROUTINE SAPINI !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! initilization of seedling cohorts with given height according to relations used in growth_seed
SUBROUTINE sapini(outunit,taxid, height, hbc, n_koh,iage)
USE data_species
USE data_stand
use data_help
IMPLICIT none
REAL :: height,hbc,hhelp
INTEGER :: outunit,n_koh ,taxid,iage
REAL :: x1,x2,xacc,shelp
real :: rtflsp, sapwood
external sapwood
external rtflsp
! Shootbiomass kg from height (cm), originally x_sap [mg]
hhelp = height * 100.
IF (taxid.ne.2) x_sap = exp(( LOG(hhelp)-LOG(spar(taxid)%pheight1))/spar(taxid)%pheight2)/1000000.
IF (taxid.eq.2) THEN
x1 = 1.
x2 = 2.
xacc=(1.0e-10)*(x1+x2)/2
! solve equation for calculation of sapwood from height; determine root
heihelp = hhelp
hnspec = taxid
shelp=rtflsp(sapwood,x1,x2,xacc)
x_sap = (10**shelp)/1000000 ! transformation mg ---> kg
ENDIF
! leaf matter
x_fol = (spar(taxid)%seeda*(x_sap** spar(taxid)%seedb)) ![kg]
! fine root matter rough estimate
x_frt = x_fol
! cross sectional area of heartwood
x_ahb = 0.
x_hrt = 0.
IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,hhelp,hbc,iage,n_koh,taxid
END subroutine sapini
FUNCTION ran0(idum)
INTEGER idum,IA,IM,IQ,IR,MASK
REAL ran0,AM
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,MASK=123459876)
INTEGER kran
idum=ieor(idum,MASK)
kran=idum/IQ
idum=IA*(idum-kran*IQ)-IR*kran
IF (idum.lt.0) idum=idum+IM
ran0=AM*idum
idum=ieor(idum,MASK)
RETURN
END
! (C) Copr. 1986-92 Numerical Recipes Software 0)+0143$!-.
SUBROUTINE header(outunit,infile,source,cform,rsap,flag_volfunc,patchsize)
! write file headder into ini-file
INTEGER :: outunit, flag_volfunc
REAL :: rsap, cform, patchsize
CHARACTER(75) :: infile
CHARACTER :: source
WRITE(outunit,'(I1,1F12.0,A32)')flag_volfunc,patchsize,' ! = volume function, patch size'
WRITE(outunit,'(A15,A1,A13,A80)') '! data source= ',source,' source file= ',infile
WRITE(outunit,'(A57)') '! sapwood fraction and form factor now dynamic per cohort '
WRITE(outunit,'(a37)')'! 4C Tree Initialization File (Stand)'
WRITE(outunit,'(a1)')'!'
WRITE(outunit,'(a51)')'! contains the following data (single tree values):'
WRITE(outunit,'(a1)')'!'
WRITE(outunit,'(a31)')'! x_fol: foliage biomass (kg)' WRITE(outunit,'(a33)')'! x_frt: fine root biomass (kg)'
WRITE(outunit,'(a31)')'! x_sap: sapwood biomass (kg)'
WRITE(outunit,'(a33)')'! x_hrt: heartwood biomass (kg)'
WRITE(outunit,'(a65)')'! x_Ahb: cross sectional area of heartwood at stem base (cm**2)'
WRITE(outunit,'(a27)')'! height: tree height (cm)'
WRITE(outunit,'(a27)')'! x_hbole: bole height (cm)'
WRITE(outunit,'(a27)')'! x_age: tree age (years)'
WRITE(outunit,'(a26)')'! n: number of trees'
WRITE(outunit,'(a35)')'! sp: species (integer number)'
WRITE(outunit,'(a33)')'! DC: diameter at crown base'
WRITE(outunit,'(a37)')'! DBH: diameter at breast height'
WRITE(outunit,'(a1)')'!'
WRITE(outunit,'(a120)')'! x_fol x_frt x_sap x_hrt x_Ahb height x_hbole x_age n sp DC DBH'
END subroutine header
FUNCTION crown_base(height,c1,c2,bhd)
IMPLICIT NONE
REAL crown_base
REAL height,bhd,c1,c2
!--- estimate crown starting height according to Nagel (1995)
crown_base=height*(1.-exp(-1.*(c1+c2*height/bhd)**2))
END function crown_base
Function crown_base_eg(height,bhd)
IMPLICIT NONE
real crown_base_eg
real height, bhd
! after Nutto etal. 2006
crown_base_eg= -5.12 -0.407*bhd + 1.193*height
if ( crown_base_eg.lt. 0.) crown_base_eg = 0.
END function crown_base_eg
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE fdfahc(X,F,DF,asap,bhd,hbc,height,Vd,J)
USE data_par
USE data_simul
use data_help
IMPLICIT none
REAL X,F,DF,asap,bhd,hbc,height,Vd,C,dCdX
INTEGER J
fail=0
IF (asap+X.LE.0) THEN
WRITE(*,*) 'negative root at calculation C in fdfahc, program will stop'
STOP
ENDIF
C=(bhd-(4./PI*(asap+X))**0.5*h_breast/hbc)/(1.-h_breast/hbc)
dCdX=(-h_breast)/hbc/(1.-h_breast/hbc)/(4./PI*(asap+X))**0.5*2./PI
IF (X*(PI/4.*C**2.-asap).LE.0) THEN
fail=1
return
ENDIF
F=1./3.*height*(asap+X)+1./3.*hbc*(asap+PI/4.*C**2.+(X*(PI/4.*C**2.-asap))**0.5)-Vd*1000000.
DF=1./3.*(height+hbc*PI/2.*C*dCdX+hbc*0.5/(X*(PI/4.*C**2.-asap))**0.5*(PI/4.*C**2+X*PI/2.*C*dCdX-asap))
END subroutine fdfahc
FUNCTION NEWTON(X,asap,bhd,hbc,height,Vd)
use data_help
IMPLICIT NONE
REAL newton
REAL F,DF,X,DX,asap,bhd,hbc,height,Vd
INTEGER J,stepmax
! Newton is to be called with a start value for X
! a subroutine NEWFDF is to be included in the main program which
! calculates the value of the function and its derivative at X and
! returns them in the variables F and DF
PARAMETER (stepmax=5000)
DO 7 J=1,stepmax
CALL fdfAhc(X,F,DF,asap,bhd,hbc,height,Vd,J) if(fail.eq.1) return
IF(DF.EQ.0.0) THEN
DX=0.01*X
ELSE
DX=F/DF
ENDIF
Newton=X
IF(DX.GT.X) DX=X/2.
X=X-DX
IF(ABS(DX).LT.0.0005) RETURN
7 END DO
END
SUBROUTINE ini_gener_sap(outunit,taxid,age,pl_height, nplant)
USE data_stand
USE data_par
USE data_species
USE data_soil
USE data_help
USE data_plant
USE data_manag
IMPLICIT NONE
integer :: nplant, &
taxid, &
nclass, &
i,nr, &
age, &
outunit
real :: pl_height, &
height, &
hhelp, &
hbc, &
sdev, &
help, &
nstot
real :: rtflsp, sapwood
real :: hmin_est ! empirical estimated minimum height
real, dimension(:), allocatable :: hei, &
nschelp
integer,dimension(:),allocatable :: nsc
external sapwood
external rtflsp
sdev = hsdev(taxid)
if (nplant.eq.0) nplant= numplant(taxid)
height = pl_height*100
if(height .lt. 100) then
hmin_est = height - height*0.2
else
hmin_est = height - height*0.1
end if
if(nplant.eq.1) hmin_est = height
nclass= nint((height+2*sdev) - hmin_est) + 1
if(nplant.eq.1) nclass =1
if(nplant.lt.200) nclass=1
allocate(hei(nclass))
allocate(nschelp(nclass))
allocate(nsc(nclass))
nstot = 0
help = (1/(sqrt(2*pi)*sdev))
do i = 1, nclass
! height per class
hei(i) = hmin_est + (i-1)
nschelp(i) = help*exp(-((hei(i)-height)**2)/(2*(sdev)**2))
nstot = nstot + nschelp (i)
end do
! scaling of plant number per cohort do i = 1,nclass
nsc(i) = nint((nschelp(i)*nplant/nstot) + 0.5)
end do
if(nplant.eq.1) nsc(1) = nplant
do i = 1,nclass
hhelp = hei(i)*0.01
hbc=0
call sapini(outunit,taxid, hhelp, hbc,nsc(i),age)
end do
END SUBROUTINE ini_gener_sap