!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*     SUBROUTINES                                               *!
!*     - assign_DSW                                              *!
!*     - assign_Bay                                              *!
!*     - parthe_param                                            *!
!*     - data_gap_fill_DSW                                       *!
!*     - init_plenter_param                                      *!
!*     - fdfk                                                    *!
!*     FUNCTIONS                                                 *!
!*     - tax_of_BRA_id                                           *!
!*     - wachsfunc                                               *!
!*     - inv_wachsfunc                                           *!
!*     - agefunc                                                 *!
!*     - newton_plenter                                          *!
!*     - n0ofvol                                                 *!
!*                                                               *!
!*                  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 assign_DSW
! Table of species numbers and names used in Datenspeicher Waldfonds (DSW)
! data structure of species data DSW
! this table is based on the Brandenburg version of the BRA (Betriebsregelungsanweisung)
! in other states (Länder) the numbering can be different. Therefore the table must be checked against
! their definitions, when new data sources are to be used
USE data_init
IMPLICIT NONE
INTEGER         :: i,imax
! <400 conifers
! >100 <200 pines and larches
i=1
spec_nrDSW(i)=111; spec_code(i)='GKI'; GER_name(i)='Gemeine Kiefer'; LAT_name(i)='Pinus sylvertris L.' 
spec_4c(i)=3
i=i+1
spec_nrDSW(i)=112; spec_code(i)='WKI'; GER_name(i)='Weymouthkiefer'; LAT_name(i)='Pinus strobus L.'
spec_4c(i)=7
i=i+1
spec_nrDSW(i)=113; spec_code(i)='SKI'; GER_name(i)='Schwarzkiefer'; LAT_name(i)='Pinus nigra ARN.'
spec_4c(i)=3
i=i+1
spec_nrDSW(i)=114; spec_code(i)='MKI'; GER_name(i)='Murraykiefer'; LAT_name(i)='Pinus contorta DOUGL. Ex LOUD.'
spec_4c(i)=6
i=i+1
spec_nrDSW(i)=115; spec_code(i)='RKI'; GER_name(i)='Rumelische Kiefer'; LAT_name(i)='Pinus peuce GRISEB.'
spec_4c(i)=3
i=i+1
spec_nrDSW(i)=116; spec_code(i)='BKI'; GER_name(i)='Bergkiefer'; LAT_name(i)='Pinus mugo TURRA'
spec_4c(i)=6
i=i+1
spec_nrDSW(i)=117; spec_code(i)='ZKI'; GER_name(i)='Zirbelkiefer'; LAT_name(i)='Pinus cembra L.'
spec_4c(i)=6
i=i
spec_nrDSW(i)=118; spec_code(i)='PKI'; GER_name(i)='Gelbkiefer'; LAT_name(i)='Pinus ponderosa DOUGL. Ex LAWS.'
spec_4c(i)=7
i=i+1
spec_nrDSW(i)=119; spec_code(i)='KIS'; GER_name(i)='Sonst. Kiefern'
spec_4c(i)=3
i=i+1
spec_nrDSW(i)=171; spec_code(i)='ELA'; GER_name(i)='Europ. Lärche'; LAT_name(i)='Larix decidua MILL.'
spec_4c(i)=6
i=i+1
spec_nrDSW(i)=172; spec_code(i)='JLA'; GER_name(i)='Japan. Lärche'; LAT_name(i)='Larix kaempferi (LAMB.) CARR.'
spec_4c(i)=6
i=i+1
spec_nrDSW(i)=173; spec_code(i)='HLA'; GER_name(i)='Hybridlärche'; LAT_name(i)='Larix x eurolepis HENRY'
spec_4c(i)=6
i=i+1
spec_nrDSW(i)=179; spec_code(i)='LAS'; GER_name(i)='Sonst. Lärchen'
spec_4c(i)=6
! >2oo <300 spruces
i=i+1
spec_nrDSW(i)=211; spec_code(i)='GFI'; GER_name(i)='Gemeine Fichte'; LAT_name(i)='Picea abies (L.) KARST.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=212; spec_code(i)='SFI'; GER_name(i)='Sitkafichte'; LAT_name(i)='Picea sitchensis (BONG.) CARR.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=213; spec_code(i)='WFI'; GER_name(i)='Weißfichte'; LAT_name(i)='Picea glauca (MOENCH) VOSS'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=214; spec_code(i)='OFI'; GER_name(i)='Omorikafichte'; LAT_name(i)='Picea omorika (PANC.) PURK.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=215; spec_code(i)='BFI'; GER_name(i)='Stechfichte, Blaufichte'; LAT_name(i)='Picea pungens ENGELM. + P.p. Glauca'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=216; spec_code(i)='EFI'; GER_name(i)='Engelmannfichte'; LAT_name(i)='Picea engelmannii ENGELM.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=217; spec_code(i)='MFI'; GER_name(i)='Schwarzfichte'; LAT_name(i)='Picea mariana (MILL.) B. S. P.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=218; spec_code(i)='RFI'; GER_name(i)='Rotfichte'; LAT_name(i)='Picea rubens SARG.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=219; spec_code(i)='FIS'; GER_name(i)='Sonst. Fichten'
spec_4c(i)=2
! >300 <400 firs, douglas fir, thuja, hemlock fir
i=i+1
spec_nrDSW(i)=311; spec_code(i)='WTA'; GER_name(i)='Weißtanne'; LAT_name(i)='Abies alba MILL.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=312; spec_code(i)='KTA'; GER_name(i)='Küstentanne'; LAT_name(i)='Abies grandis (D. DON) LINDL.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=313; spec_code(i)='CTA'; GER_name(i)='Coleradotanne'; LAT_name(i)='Abies concolor (GORD. et GLEND.) LINDL.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=314; spec_code(i)='NTA'; GER_name(i)='Nordmanntanne'; LAT_name(i)='Abies nordmanniana (STEV.) SPACH.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=315; spec_code(i)='ETA'; GER_name(i)='Amerikanische Edeltanne'; LAT_name(i)='Abies procera REHD.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=316; spec_code(i)='HTA'; GER_name(i)='Nikkotanne'; LAT_name(i)='Abies homolepis SIEB. et ZUCC.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=317; spec_code(i)='VTA'; GER_name(i)='Veitchtanne'; LAT_name(i)='Abies veitchii LINDL.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=319; spec_code(i)='TAS'; GER_name(i)='Sonst. Tannen'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=321; spec_code(i)='GDG'; GER_name(i)='Grüne Douglasie'; LAT_name(i)='Pseudotsuga menziesii (MIRBEL) FRANCO var. menziesii'
spec_4c(i)=10
i=i+1
spec_nrDSW(i)=322; spec_code(i)='BDG'; GER_name(i)='Blaue Douglasie'; LAT_name(i)='Pseudotsuga menziesii var. glauca (BEISSN.) FRANCO'
spec_4c(i)=10
i=i+1
spec_nrDSW(i)=323; spec_code(i)='CDG'; GER_name(i)='Graue Douglasie'; LAT_name(i)='Pseudotsuga menziesii var. caesia (SCHWERIN) FRANCO'
spec_4c(i)=10
i=i+1
spec_nrDSW(i)=329; spec_code(i)='DGS'; GER_name(i)='Sonst. Douglasien'
spec_4c(i)=10
i=i+1
spec_nrDSW(i)=331; spec_code(i)='RLB'; GER_name(i)='Riesenlebensbaum'; LAT_name(i)='Thuja plicata DONN ex D. DON'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=332; spec_code(i)='MLB'; GER_name(i)='Morgenländischer Lebensbaum'; LAT_name(i)='Thuja orientalis L.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=333; spec_code(i)='ALB'; GER_name(i)='Abendländischer Lebensbaum'; LAT_name(i)='Thuja occidentalis L.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=339; spec_code(i)='LBS'; GER_name(i)='Sonst. Lebensb.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=341; spec_code(i)='LLZ'; GER_name(i)='Lawson-Scheinzypresse'; LAT_name(i)='Chamaecyparis lawsoniana (A. MURR.) PARL.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=349; spec_code(i)='LZS'; GER_name(i)='Sonstige Scheinzypressen'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=351; spec_code(i)='KHT'; GER_name(i)='Kanadische Hemlockstanne'; LAT_name(i)='Tsuga canadensis (L.) CARR.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=352; spec_code(i)='WHT'; GER_name(i)='Westamerikanische Hemlockstanne'; LAT_name(i)='Tsuga heterophylla (RAF.) SARG.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=359; spec_code(i)='HTS'; GER_name(i)='Hemlockstannen'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=361; spec_code(i)='EIB'; GER_name(i)='(Beeren-) Eibe'; LAT_name(i)='Taxus baccata L.'
spec_4c(i)=2
i=i+1
spec_nrDSW(i)=371; spec_code(i)='GWA'; GER_name(i)='Gemeiner Wachholder'; LAT_name(i)='Juniperus communis L.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=379; spec_code(i)='WAS'; GER_name(i)='Sonstige Wacholder'; LAT_name(i)='Juniperus spec.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=381; spec_code(i)='MA '; GER_name(i)='Mammutbäume'; LAT_name(i)='Metasequoia spec., Sequia spec.'
spec_4c(i)=10
i=i+1
spec_nrDSW(i)=399; spec_code(i)='NDS'; GER_name(i)='Sonstige Nadelbaumarten'
spec_4c(i)=2
! >400  broad leaved trees
! >400 <500  oaks
i=i+1
spec_nrDSW(i)=410; spec_code(i)='EI '; GER_name(i)='Eichen-Bastarde (SEI-/TEI-Bastarde)'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=411; spec_code(i)='SEI'; GER_name(i)='Stieleiche'; LAT_name(i)='Quercus robur L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=412; spec_code(i)='TEI'; GER_name(i)='Traubeneiche'; LAT_name(i)='Quercus petraea (MAT.) LIEBL.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=413; spec_code(i)='ZEI'; GER_name(i)='Zerreiche Quercus cerris L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=414; spec_code(i)='PEI'; GER_name(i)='Sumpfeiche'; LAT_name(i)='Quercus palustris MUENCHH.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=415; spec_code(i)='REI'; GER_name(i)='Roteiche'; LAT_name(i)='Quercus rubra L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=419; spec_code(i)='EIS'; GER_name(i)='Sonst. Eichen'
spec_4c(i)=4
! >500 <600 Buchen, beeches
i=i+1
spec_nrDSW(i)=511; spec_code(i)='RBU'; GER_name(i)='Rotbuche'; LAT_name(i)='Fagus sylvatica L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=519; spec_code(i)='BUS'; GER_name(i)='Sonst. Buchen'
spec_4c(i)=1
i=i+1
! >600 <700  Hard wood specieces, except oaks and beeches
spec_nrDSW(i)=611; spec_code(i)='HBU'; GER_name(i)='Hainbuche'; LAT_name(i)='Carpinus betulus L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=621; spec_code(i)='GES'; GER_name(i)='Gemeine Esche'; LAT_name(i)='Fraxinus excelsior L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=622; spec_code(i)='WES'; GER_name(i)='Weißesche'; LAT_name(i)='Fraxinus americana L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=629; spec_code(i)='ESS'; GER_name(i)='Sonstige Eschen'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=631; spec_code(i)='BAH'; GER_name(i)='Bergahorn'; LAT_name(i)='Acer pseudoplatanus L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=632; spec_code(i)='SAH'; GER_name(i)='Spitzahorn'; LAT_name(i)='Acer platanoides L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=633; spec_code(i)='FAH'; GER_name(i)='Feldahorn'; LAT_name(i)='Acer campestre L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=634; spec_code(i)='IAH'; GER_name(i)='Silberahorn'; LAT_name(i)='Acer saccharinum L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=635; spec_code(i)='EAH'; GER_name(i)='Eschenblättriger Ahorn'; LAT_name(i)='Acer negundo L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=639; spec_code(i)='AHS'; GER_name(i)='Sonst. Ahornarten'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=641; spec_code(i)='BRU'; GER_name(i)='Bergrüster,Bergulme'; LAT_name(i)='Ulmus glabra HUDS.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=642; spec_code(i)='WRU'; GER_name(i)='Weißrüster, Flatterulme'; LAT_name(i)='Ulmus laevis PALL.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=643; spec_code(i)='FRU'; GER_name(i)='Feldrüster, Feldulme'; LAT_name(i)=''; LAT_name(i)='Ulmus minor MILL.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=644; spec_code(i)='HRU'; GER_name(i)='Hölländische Rüster, Bastardulme'; LAT_name(i)='Ulmus x hollandica MILL.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=649; spec_code(i)='RUS'; GER_name(i)='(UL)  Sonstige Rüstern, (Heimische) Rüstern - Ulmen'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=651; spec_code(i)='RO '; GER_name(i)='Gem. Robinie'; LAT_name(i)='Robinia pseudoacacia L.'
spec_4c(i)=11
i=i+1
spec_nrDSW(i)=654; spec_code(i)='GLE'; GER_name(i)='Amerikanische Gleditschie'; LAT_name(i)='Gleditsia triacanthos L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=661; spec_code(i)='EK '; GER_name(i)='Edelkastanie'; LAT_name(i)='Castanea sativa MILL.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=662; spec_code(i)='NB '; GER_name(i)='Nußbaumarten'; LAT_name(i)='Juglans spec.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=665; spec_code(i)='HI '; GER_name(i)='Hickory-Arten'; LAT_name(i)='Carya spec.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=667; spec_code(i)='EHA'; GER_name(i)='Europäische Hasel'; LAT_name(i)='Corylus avellana L.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=668; spec_code(i)='BHA'; GER_name(i)='Baumhasel'; LAT_name(i)='Corylus colurna L.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=671; spec_code(i)='VKB'; GER_name(i)='Vogelkirsche (-baum)'; LAT_name(i)='Cerasus avium (L.) MOENCH ssp. Avium (Prunus avium L.)'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=672; spec_code(i)='GTK'; GER_name(i)='Gewöhnliche Traubenkirsche'; LAT_name(i)='Padus avium MILL. (Prunus padus L.)'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=673; spec_code(i)='STK'; GER_name(i)='Spätbl. Traubenk.'; LAT_name(i)='Padus serotina (EHRH.) BORKH. (Prunus serotina EHRH.)'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=675; spec_code(i)='AB '; GER_name(i)='Wildapfel (-baum)'; LAT_name(i)='Malus sylvestris MILL.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=676; spec_code(i)='BB '; GER_name(i)='Wildbirne (-baum)'; LAT_name(i)='Pyrus spec.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=679; spec_code(i)='KBS'; GER_name(i)='Sonstige Obstbäume'; LAT_name(i)='Prunus spec.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=681; spec_code(i)='PLT'; GER_name(i)='Platanen'; LAT_name(i)='Platanus spec.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=699; spec_code(i)='HLS'; GER_name(i)='Sonst. Hartlaubbaumarten'
spec_4c(i)=4
! >700 <800 soft (deciduous) wood species
i=i+1
spec_nrDSW(i)=711; spec_code(i)='GBI'; GER_name(i)='Gemeine Birke, Sandbirke'; LAT_name(i)='Betula pendula ROTH.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=712; spec_code(i)='MBI'; GER_name(i)='Moorbirke'; LAT_name(i)='Betula pubescens EHRH.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=719; spec_code(i)='BIS'; GER_name(i)='Sonst. Birken'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=721; spec_code(i)='RER'; GER_name(i)='Roterle, Schwarzerle'; LAT_name(i)='Alnus glutinosa (L.) GAERTN.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=722; spec_code(i)='WER'; GER_name(i)='Weißerle, Grauerle'; LAT_name(i)='Alnus incana (L.) MOENCH'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=723; spec_code(i)='GER'; GER_name(i)='Grünerle'; LAT_name(i)='Alnus viridis (CHAIX) DC.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=731; spec_code(i)='WLI'; GER_name(i)='Winterlinde'; LAT_name(i)='Tilia cordata MILL.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=732; spec_code(i)='SLI'; GER_name(i)='Sommerlinde'; LAT_name(i)='Tilia platyphyllos SCOP.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=739; spec_code(i)='LIS'; GER_name(i)='Sonstige Linden'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=741; spec_code(i)='SPA'; GER_name(i)='Europäische Schwarzpappel'; LAT_name(i)='Populus nigra L.'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=742; spec_code(i)='HPA'; GER_name(i)='Schwarzpappel-Hybriden'; LAT_name(i)='Populus canadensis MOENCH.'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=743; spec_code(i)='TPA'; GER_name(i)='Trichocarpa-Pappel'; LAT_name(i)='Populus trichocarpa TORR. et A. GRAY ex HOOK'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=744; spec_code(i)='BPA'; GER_name(i)='Balsampappel-Hybriden'; LAT_name(i)='Populus trichocarpa x maximoviczii HENRY (Androscoggin)'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=745; spec_code(i)='GPA'; GER_name(i)='Graupappel + Hybriden'; LAT_name(i)='Populus x canescens SMITH + P. can. X grandidentata MICHX.'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=746; spec_code(i)='WPA'; GER_name(i)='Silberpappel (Weißpappel)'; LAT_name(i)='Populus Populus alba L.'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=747; spec_code(i)='AS '; GER_name(i)='Aspe'; LAT_name(i)='Populus tremula L.'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=748; spec_code(i)='HAS'; GER_name(i)='Aspen-Hybriden'; LAT_name(i)='Populus tremula l. x Populus tremuloides'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=749; spec_code(i)='PAS'; GER_name(i)='Sonst. Pappeln (z.B. Balsam-Schwarzpappel-Hybriden)'
spec_4c(i)=8
i=i+1
spec_nrDSW(i)=751; spec_code(i)='WWE'; GER_name(i)='Weißweide (Silberweide)'; LAT_name(i)='Salix alba L.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=752; spec_code(i)='BWE'; GER_name(i)='Bruchweide (Knackweide)'; LAT_name(i)='Salix fragilis L.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=753; spec_code(i)='FWE'; GER_name(i)='Fahlweide (Baumweiden-Hybriden)'; LAT_name(i)='Salix x rubens SCHRANK (= Salix alba x fragilis)'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=754; spec_code(i)='SWE'; GER_name(i)='Salweide'; LAT_name(i)='Salix caprea L.'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=759; spec_code(i)='WEB'; GER_name(i)='Baumweiden'
spec_4c(i)=5
i=i+1
spec_nrDSW(i)=761; spec_code(i)='RK '; GER_name(i)='Roßkastanie'; LAT_name(i)='Aesculus hippocastanum L.'
spec_4c(i)=1
i=i+1
spec_nrDSW(i)=771; spec_code(i)='EB '; GER_name(i)='Gemeine Eberesche'; LAT_name(i)='Sorbus aucuparia L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=772; spec_code(i)='EEB'; GER_name(i)='Edel-Eberesche'; LAT_name(i)='Sorbus a. var. Edulis DIECK'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=773; spec_code(i)='ME '; GER_name(i)='Echte Mehlbeere'; LAT_name(i)='Sorbus aria CRANTZ'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=774; spec_code(i)='EL '; GER_name(i)='Elsbeere'; LAT_name(i)='Sorbus torminalis CRANTZ'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=775; spec_code(i)='SG '; GER_name(i)='Speierling'; LAT_name(i)='Sorbus domestica L.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=779; spec_code(i)='MES'; GER_name(i)='Sonst. Mehlbeeren'; LAT_name(i)='Sorbus spec.'
spec_4c(i)=4
i=i+1
spec_nrDSW(i)=781; spec_code(i)='GO '; GER_name(i)='Gemeiner Götterbaum'; LAT_name(i)='Ailanthus altissima (MILL.) SWINGLE'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=786; spec_code(i)='MB '; GER_name(i)='Maulbeeren'; LAT_name(i)='Morus spec.'
spec_4c(i)=0
i=i+1
spec_nrDSW(i)=799; spec_code(i)='WLS'; GER_name(i)='Sonstige Weichlaubbaumarten'
spec_4c(i)=5
imax=i
spnum_for_DSW=0
DO i=1,imax
   spnum_for_DSW(spec_nrDSW(i))=i
ENDDO
END ! subroutine assign_DSW

! Baumartenkodierung fuer Bayern 2003
!
! BA-ID	BA_Gruppe	BA_Typ	Kurzname	Name	Ertragstafel-ID	X1	X2
! 10	1	1	Fi	Fichte	10	1	0.81
! 11	1	1	OFi	Omorikafichte	10	1	0.81
! 12	1	1	SFi	Sitkafichte	10	1	0.81
! 20	2	1	Kie	Kiefer	20	1	0.79
! 21	2	1	Stro	Strobe	10	1	0.79
! 22	2	1	SKie	Schwarzkiefer	20	1	0.79
! 23	2	1	Spir	Spirke	20	1	0.79
! 24	2	1	Zir	Zirbe	20	1	0.79
! 25	2	1	Lat	Latsche			
! 30	3	1	Ta	Tanne	30	6	0.81
! 35	3	1	Eib	Eibe	30	6	0.81
! 40	4	1	ELae	Lärche (europ.)	40	7	0.72
! 41	4	1	JLae	Jap.Lärche	41	7	0.72
! 50	5	1	Dgl	Douglasie	50	8	0.79
! 60	6	2	Bu	Buche	60	1	0.846
! 61	8	2	HBu	Hainbuche	60	1	0.81
! 62	9	2	WLi	Winterlinde	60	1	0.81
! 63	9	2	Es	Esche	63	12	0.81
! 64	9	2	BAh	Bergahorn	63	12	0.81
! 65	8	2	SAh	Spitzahorn	63	12	0.81
! 66	8	2	FAh	Feldahorn	75	13	0.81
! 67	8	2	Rob	Robinie	60	1	0.81
! 68	9	2	Kir	Kirsche	60	1	0.81
! 69	9	2	Wob	Wildobst	60	1	0.81
! 70	7	2	Ei	Eiche	70	9	0.79
! 71	9	2	REi	Roteiche	71	10	0.79
! 72	9	2	Ul	Ulme	70	9	0.81
! 73	9	2	Elsb	Elsbeere	70	9	0.81
! 74	8	2	Mebe	Mehlbeere	70	9	0.81
! 75	8	2	SBi	Sandbirke	75	13	0.81
! 76	8	2	Vobe	Vogelbeere	75	13	0.81
! 77	9	2	Kast	Edelkastanie	60	1	0.81
! 78	9	2	Nuss	Nußarten	60	1	0.81
! 79	9	2	Spei	Speierling	70	9	0.81
! 80	8	2	SLbh	Sonst.Laubholz	75	13	0.81
! 81	8	2	WErl	Weißerle	75	13	0.81
! 82	8	2	As	Aspe	75	13	0.81
! 83	8	2	Pa	Pappel	83	14	0.81
! 84	8	2	Wei	Weide	75	13	0.81
! 85	8	2	GErl	Grünerle	75	13	0.81
! 86	8	2	SErl	Schwarzerle	86	11	0.81
! 87	9	2	ELbh	Edellaubholz	63	12	0.81
! 88	9	2	SLi	Sommerlinde	60	1	0.81
! 89	8	2	MBi	Moorbirke	60	1	0.81
! 90	1	1	SNdh	Sonst.Nadelholz	10	1	0.81


SUBROUTINE assign_BAY
! Table of species numbers and names used in Bavaria (Bayern)
USE data_init
IMPLICIT NONE
INTEGER         :: i,imax

! <60 conifers and 90 =  other conifers
! >=10 <20  spruces
i=1
spec_nrBAY(i)=10; spec_code(i)='FI'; GER_name(i)='Fichte'; LAT_name(i)='Picea abies (L.) KARST.'; spec_4c(i)=2
i=i+1
spec_nrBAY(i)=11; spec_code(i)='OFI'; GER_name(i)='Omorikafichte'; LAT_name(i)='Picea omorika (PANC.) PURK.'; spec_4c(i)=2
i=i+1
spec_nrBAY(i)=12; spec_code(i)='SFI'; GER_name(i)='Sitkafichte'; LAT_name(i)='Picea sitchensis (BONG.) CARR.'; spec_4c(i)=2
! >= 20 <30 Scots pine
i=i+1
spec_nrBAY(i)=20; spec_code(i)='KIE'; GER_name(i)='Kiefer'; LAT_name(i)='Pinus sylvertris L.'; spec_4c(i)=3
i=i+1
spec_nrBAY(i)=21; spec_code(i)='STR'; GER_name(i)='Strobe'; LAT_name(i)='Pinus strobus L.'; spec_4c(i)=3
! Weymouth pine (Pinus strobus) is assumably classed as spruce in Bavaria, so it is adopted here  
i=i+1
spec_nrBAY(i)=22; spec_code(i)='SKI'; GER_name(i)='Schwarzkiefer'; LAT_name(i)='Pinus nigra ARN.'; spec_4c(i)=3
i=i+1
spec_nrBAY(i)=23; spec_code(i)='SPI'; GER_name(i)='Spirke'; LAT_name(i)='Pinus uncinata RAMOND'; spec_4c(i)=3
i=i+1
! Spirke/local mountain pine no distinction between mountain pine  mugo pine
spec_nrBAY(i)=24; spec_code(i)='ZKI'; GER_name(i)='Zirbelkiefer'; LAT_name(i)='Pinus cembra L.'; spec_4c(i)=3
i=i+1
spec_nrBAY(i)=25; spec_code(i)='BKI'; GER_name(i)='Latsche'; LAT_name(i)='Pinus mugo TURRA'; spec_4c(i)=3
! arolla pine
! >= 30 <40   firs
i=i+1
spec_nrBAY(i)=30; spec_code(i)='TA'; GER_name(i)='Tanne'; LAT_name(i)='Abies alba MILL.'; spec_4c(i)=2
i=i+1
spec_nrBAY(i)=35; spec_code(i)='EIB'; GER_name(i)='Eibe'; LAT_name(i)='Taxus baccata L.'; spec_4c(i)=2
! >= 40 <50  larches
i=i+1
spec_nrBAY(i)=40; spec_code(i)='ELA'; GER_name(i)='Europ. Lärche'; LAT_name(i)='Larix decidua MILL.'; spec_4c(i)=6
i=i+1
spec_nrBAY(i)=41; spec_code(i)='JLA'; GER_name(i)='Japan. Lärche'; LAT_name(i)='Larix kaempferi (LAMB.) CARR.'; spec_4c(i)=6
! >=50 <60  douglas firs
i=i+1
spec_nrBAY(i)=50; spec_code(i)='DGL'; GER_name(i)='Douglasie'; LAT_name(i)='Pseudotsuga menziesii (MIRBEL) FRANCO var. menziesii'; spec_4c(i)=2
i=i+1
! >= 60 < deciduous tree species
spec_nrBAY(i)=60; spec_code(i)='BU'; GER_name(i)='Buche'; LAT_name(i)='Fagus sylvatica'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=61; spec_code(i)='HBU'; GER_name(i)='Hainbuche'; LAT_name(i)='Carpinus betulus L.'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=62; spec_code(i)='WLi'; GER_name(i)='Winterlinde'; LAT_name(i)='Tilia cordata'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=63; spec_code(i)='Es'; GER_name(i)='Esche'; LAT_name(i)='Fraxinus excelsior'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=64; spec_code(i)='BAh'; GER_name(i)='Bergahorn'; LAT_name(i)='Acer pseudoplatanus'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=65; spec_code(i)='SAh'; GER_name(i)='Spitzahorn'; LAT_name(i)='Acer platanoides'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=66; spec_code(i)='FAh'; GER_name(i)='Feldahorn'; LAT_name(i)='Acer campestre'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=67; spec_code(i)='Rob'; GER_name(i)='Robinie'; LAT_name(i)='Robinia pseudoacacia L.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=68; spec_code(i)='Kir'; GER_name(i)='Kirsche'; LAT_name(i)='??? L.'; spec_4c(i)=0
i=i+1
spec_nrBAY(i)=69; spec_code(i)='Wob'; GER_name(i)='Wildobst'; LAT_name(i)='???'; spec_4c(i)=0
i=i+1
spec_nrBAY(i)=70; spec_code(i)='Ei'; GER_name(i)='Eiche'; LAT_name(i)='Quercus sp.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=71; spec_code(i)='REi'; GER_name(i)='Roteiche'; LAT_name(i)='Quercus rubra L.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=72; spec_code(i)='Ul'; GER_name(i)='Ulme'; LAT_name(i)='Ulmus sp.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=73; spec_code(i)='Elsb'; GER_name(i)='Elsbeere'; LAT_name(i)='Sorbus torminalis CRANTZ'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=74; spec_code(i)='Mebe'; GER_name(i)='Mehlbeere'; LAT_name(i)='Sorbus aria CRANTZ'; spec_4c(i)=0
i=i+1
spec_nrBAY(i)=75; spec_code(i)='SBi'; GER_name(i)='Sandbirke'; LAT_name(i)='Betula pendula ROTH'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=76; spec_code(i)='Vobe'; GER_name(i)='Vogelbeere'; LAT_name(i)='Sorbus aucuparia L.'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=77; spec_code(i)='Kast'; GER_name(i)='Edelkastanie'; LAT_name(i)='Castanea sativa MILL.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=78; spec_code(i)='Nuss'; GER_name(i)='Nußarten'; LAT_name(i)='Juglans spec.'; spec_4c(i)=4
i=i+1
spec_nrBAY(i)=79; spec_code(i)='Spei'; GER_name(i)='Speierling'; LAT_name(i)='Sorbus domestica L.'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=80; spec_code(i)='SLbh'; GER_name(i)='Sonst. Laubholz'; LAT_name(i)=''; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=81; spec_code(i)='WErl'; GER_name(i)='Weißerle'; LAT_name(i)='Alnus incana (L.) MOENCH'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=82; spec_code(i)='As'; GER_name(i)='Aspe'; LAT_name(i)='Populus tremula L.'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=83; spec_code(i)='Pa'; GER_name(i)='Pappel'; LAT_name(i)='Populus spec.'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=84; spec_code(i)='Wei'; GER_name(i)='Weide'; LAT_name(i)='Salix spec.'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=85; spec_code(i)='GErl'; GER_name(i)='Grünerle'; LAT_name(i)='Alnus viridis (CHAIX) DC.'; spec_4c(i)=0
i=i+1
spec_nrBAY(i)=86; spec_code(i)='SErl'; GER_name(i)='Schwarzerle'; LAT_name(i)='Alnus glutinosa (L.) GAERTN.'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=87; spec_code(i)='ELbh'; GER_name(i)='Edellaubholz'; LAT_name(i)=''; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=88; spec_code(i)='SLi'; GER_name(i)='Sommerlinde'; LAT_name(i)='Tilia platyphyllos SCOP.'; spec_4c(i)=1
i=i+1
spec_nrBAY(i)=89; spec_code(i)='SLi'; GER_name(i)='Moorbirke'; LAT_name(i)='Betula pubescens EHRH.'; spec_4c(i)=5
i=i+1
spec_nrBAY(i)=90; spec_code(i)='SNdh'; GER_name(i)='Sonst. Nadelholz'; LAT_name(i)=''; spec_4c(i)=2
i=i+1

imax = i-1
spnum_for_DSW=0
DO i=1,imax
   spnum_for_DSW(spec_nrBay(i))=i
ENDDO
END ! subroutine assign_BAY


FUNCTION tax_of_BRA_id(BRAid)
USE data_init
IMPLICIT NONE
INTEGER BRAid, tax_of_BRA_id
   tax_of_BRA_id=spec_4c(spnum_for_DSW(BRAid))
END

SUBROUTINE parthe_param(species,schichtin,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
USE data_init
IMPLICIT NONE
INTEGER spezies,schicht,species,schichtin
REAL    hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe
! assignment of parameter values for data gap filling on height and diameter
spezies=species
schicht=schichtin
IF(schicht==50) schicht=10
IF(schicht==20) GOTO 2222
1111 CONTINUE
IF(spezies==111) THEN                          ! Pinus sylvestris
   hymax_Parthe=23.74697;    hb_Parthe=0.003;    hT_Parthe=21.94225
   dymax_Parthe=34.83703;    db_Parthe=0.00146;  dT_Parthe=34.72167
ELSEIF (spezies==211) THEN                     ! Picea abies
   hymax_Parthe=25.93201;    hb_Parthe=0.00186;  hT_Parthe=26.76250
   dymax_Parthe=42.86844;    db_Parthe=0.00029;  dT_Parthe=15.53258
ELSEIF (spezies==171) THEN                     !  Larix decidua
   hymax_Parthe=25.65709;    hb_Parthe=0.00295;  hT_Parthe=18.05441
   dymax_Parthe=50.63337;    db_Parthe=0.00027;  dT_Parthe=9.03576
ELSEIF (spezies==711) THEN                     ! Betula pendula
   hymax_Parthe=24.63548;    hb_Parthe=0.00298;  hT_Parthe=18.02402
   dymax_Parthe=36.45272;    db_Parthe=0.00112;  dT_Parthe=36.2542
ELSEIF (spezies==411.AND.schicht==10) THEN     !  Quercus robur
   hymax_Parthe=22.22929;    hb_Parthe=0.00224;  hT_Parthe=24.73157
   dymax_Parthe=87.64567;    db_Parthe=0.00012;  dT_Parthe=89.0633
ELSEIF (spezies==411.AND.schicht==40) THEN     !  Quercus robur
   hymax_Parthe=14.34897;    hb_Parthe=0.00970;  hT_Parthe=20.76731
   dymax_Parthe=12.78134;    db_Parthe=0.02083;  dT_Parthe=25.9982
ELSEIF (spezies==412) THEN                     !  Quercus petraea
   hymax_Parthe=22.39128;    hb_Parthe=0.003;    hT_Parthe=25.4039
   dymax_Parthe=54.13989;    db_Parthe=0.00037;  dT_Parthe=62.1369
ELSEIF (spezies==511.AND.schicht==10) THEN     ! Fagus sylvatica
   hymax_Parthe=28.6865;     hb_Parthe=0.00172;  hT_Parthe=28.46973
   dymax_Parthe=68.5734;     db_Parthe=0.00032;  dT_Parthe=73.12856
ELSEIF (spezies==511.AND.schicht==40) THEN     ! Fagus sylvatica
   hymax_Parthe=31.28959;    hb_Parthe=0.00162;  hT_Parthe=39.51603
   dymax_Parthe=21.01226;    db_Parthe=0.00363;  dT_Parthe=32.94303
ELSEIF (spezies==631) THEN                     ! Acer pseudoplatanus
   hymax_Parthe=28.36913;    hb_Parthe=0.00123;  hT_Parthe=12.72464
   dymax_Parthe=63.8451;     db_Parthe=0.00016;  dT_Parthe=19.84293
ELSEIF (spezies==621) THEN                     ! Fraxinus excelsior
   hymax_Parthe=28.69626;    hb_Parthe=0.00138;  hT_Parthe=15.23287
   dymax_Parthe=76.37174;    db_Parthe=0.0001;   dT_Parthe=16.90759
ELSEIF (spezies==611.AND.schicht==10) THEN     ! Carpinus betulus
   hymax_Parthe=24.60247;     hb_Parthe=0.00132;  hT_Parthe=11.40522
   dymax_Parthe=45.57378;     db_Parthe=0.00047;  dT_Parthe=55.59576
ELSEIF (spezies==611.AND.schicht==40) THEN     !  Carpinus betulus
   hymax_Parthe=19.04968;     hb_Parthe=0.00174;  hT_Parthe=5.76216
   dymax_Parthe=38.45864;     db_Parthe=0.00042;  dT_Parthe=46.93101
ELSEIF (spezies==731.AND.schicht==10) THEN     !  Tilia cordata
   hymax_Parthe=27.69013;     hb_Parthe=0.00156;  hT_Parthe=23.76142
   dymax_Parthe=50.06284;     db_Parthe=0.00044;  dT_Parthe=53.24075
ELSEIF (spezies==731.AND.schicht==40) THEN     ! Tilia cordata
   hymax_Parthe=17.46179;     hb_Parthe=0.00371;  hT_Parthe=19.00039
   dymax_Parthe=13.19608;     db_Parthe=0.00586;  dT_Parthe=16.4324
ELSE
   ! if no parameters provided for the species then the parameters for the
   ! assigned 4c species will be used
   IF(spec_4c(spnum_for_DSW(spezies))==1)spezies=511
   IF(spec_4c(spnum_for_DSW(spezies))==2)spezies=211
   IF(spec_4c(spnum_for_DSW(spezies))==3)spezies=111
   IF(spec_4c(spnum_for_DSW(spezies))==4)spezies=411
   IF(spec_4c(spnum_for_DSW(spezies))==5)spezies=711
   IF(spec_4c(spnum_for_DSW(spezies))==6)spezies=171
   GOTO 1111
ENDIF
! assignment of parameter values for missing data generation on
spezies=species
2222 CONTINUE
IF(spezies==111) THEN          !  Pinus sylvestris
   uh_Parthe=24;    um_Parthe=3.462;   un_Parthe=110;   uxu_Parthe=30
ELSEIF (spezies==211) THEN     ! Picea abies; Larix parameters used
   uh_Parthe=25;    um_Parthe=0.417;   un_Parthe=90;    uxu_Parthe=30
ELSEIF (spezies==171) THEN     !  Larix decidua
   uh_Parthe=25;    um_Parthe=0.417;   un_Parthe=90;    uxu_Parthe=30
ELSEIF (spezies==711) THEN     !  Betula pendula; Larix parameters used
   uh_Parthe=25;    um_Parthe=0.417;   un_Parthe=90;    uxu_Parthe=30
ELSEIF (spezies==411) THEN     !  Quercus robur
   uh_Parthe=24;    um_Parthe=1.63;    un_Parthe=145;   uxu_Parthe=40
ELSEIF (spezies==412) THEN     ! Quercus petraea
   uh_Parthe=23;    um_Parthe=0.395;   un_Parthe=150;   uxu_Parthe=40
ELSEIF (spezies==511) THEN     ! Fagus sylvatica
   uh_Parthe=28;    um_Parthe=1.17;    un_Parthe=125;   uxu_Parthe=45
ELSEIF (spezies==621) THEN     ! Fraxinus excelsior; Fagus parameters used
   uh_Parthe=28;    um_Parthe=1.17;    un_Parthe=125;   uxu_Parthe=45
ELSEIF (spezies==731) THEN     !  Tilia cordata; Fagus parameters used
   uh_Parthe=28;    um_Parthe=1.17;    un_Parthe=125;   uxu_Parthe=45
ELSEIF (spezies==611) THEN     !  Carpinus betulus; Fagus parameters used except height
   uh_Parthe=23;    um_Parthe=1.17;    un_Parthe=125;   uxu_Parthe=45
ELSE
   ! if no parameters provided for the species then the parameters for the
   ! assigned 4c species will be used
   IF(spec_4c(spnum_for_DSW(spezies))==1)spezies=511
   IF(spec_4c(spnum_for_DSW(spezies))==2)spezies=211
   IF(spec_4c(spnum_for_DSW(spezies))==3)spezies=111
   IF(spec_4c(spnum_for_DSW(spezies))==4)spezies=411
   IF(spec_4c(spnum_for_DSW(spezies))==5)spezies=711
   IF(spec_4c(spnum_for_DSW(spezies))==6)spezies=171
   GOTO 2222
ENDIF
END ! subroutine parthe_param

FUNCTION wachsfunc(x,ymax,b,T)
! data gap filling parameters for height and diameter as a function of age
IMPLICIT NONE
INTEGER x
REAL ymax,b,T,wachsfunc
   wachsfunc=ymax*(1.-(1./(1.+(EXP(b*ymax))**(x-T)-(EXP(b*ymax))**(-T))))
END   ! function wachsfunc

FUNCTION inv_wachsfunc(x,ymax,b,T)
! inverse function of wachsfunc for retrieval of age corresponding to a given diameter
IMPLICIT NONE
REAL x,ymax,b,T
INTEGER inv_wachsfunc
   inv_wachsfunc=NINT(LOG(1./(1.-x/ymax)-1.+(EXP(b*ymax))**(-T))/(b*ymax)+T)
END   ! function inv_wachsfunc

FUNCTION agefunc(x,m,xu,n)
! data gap filling parameters for age as a function of diameter for seed trees
IMPLICIT NONE
REAL x,m,xu,n,agefunc
   agefunc=m*(x-xu)+n
END   ! function agefunc

SUBROUTINE init_plenter_param
! determines and sets ages at which dbh of 2 cm is reached and harvest age in plenter wald respectively
USE data_init
IMPLICIT NONE
INTEGER inv_wachsfunc
REAL    hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe
    high_age(1)=180    ! average estimated harvest age in plenter wald for Fagus sylvatica
    high_age(2)=140    ! average estimated harvest age in plenter wald for Picea abies
    high_age(3)=170    ! average estimated harvest age in plenter wald for Pinus silvestris
    high_age(4)=190    ! average estimated harvest age in plenter wald for Quercus sp.
    CALL parthe_param(511,10,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
    low_age(1)=inv_wachsfunc(2.,dymax_Parthe,db_Parthe,dT_Parthe)
    CALL parthe_param(211,10,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
    low_age(2)=inv_wachsfunc(2.,dymax_Parthe,db_Parthe,dT_Parthe)
    CALL parthe_param(111,10,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
    low_age(3)=inv_wachsfunc(2.,dymax_Parthe,db_Parthe,dT_Parthe)
    CALL parthe_param(411,10,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
    low_age(4)=inv_wachsfunc(2.,dymax_Parthe,db_Parthe,dT_Parthe)
    END ! subroutine init_plenter_param

SUBROUTINE data_gap_fill_DSW(i)
! fills gaps in input data
USE data_init
USE data_par
USE data_simul
USE data_species
IMPLICIT NONE
INTEGER i,n0ofvol,inv_wachsfunc
REAL    formfactor,wachsfunc,agefunc,k_age,newton_plenter
REAL    hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe
LOGICAL init_plent
  IF(ngroups(i)%taxid==2.OR.ngroups(i)%taxid==3) THEN
      formfactor=0.45
   ELSE
      formfactor=0.5
   ENDIF
   CALL parthe_param(ngroups(i)%BRAid,ngroups(i)%schicht,hymax_Parthe,hb_Parthe,hT_Parthe,dymax_Parthe,db_Parthe,dT_Parthe,uh_Parthe,um_Parthe,un_Parthe,uxu_Parthe)
   IF(ngroups(i)%schicht==20) THEN
      ! gap filling for Überhälter (seed or shelter trees)
      ngroups(i)%mhoe=uh_Parthe
      ngroups(i)%alter=agefunc(ngroups(i)%dm,um_Parthe,uxu_Parthe,un_Parthe)
      ngroups(i)%baumzahl=NINT(ngroups(i)%volume/(PI/4.*ngroups(i)%dm**2*ngroups(i)%mhoe*formfactor))
      IF(ngroups(i)%baumzahl==0.AND.ngroups(i)%volume/(PI/4.*ngroups(i)%dm**2*ngroups(i)%mhoe*formfactor)>=0.) ngroups(i)%baumzahl=1
   ELSEIF(ngroups(i)%schicht==10.OR.ngroups(i)%schicht==40) THEN
      ! gap filling for Oberstand and Unterstand (upper and lower canopy strata)
      ! with missing diameter and/or height information
      IF(ngroups(i)%alter==0.) CALL error_mess(ngroups(i)%locid,'no age information for stand: ',REAL(ngroups(i)%locid))
      IF(ngroups(i)%alter==0.) WRITE(8999,*) i,ngroups(i)%locid,ngroups(i)%BRAid,ngroups(i)%alter
      IF(ngroups(i)%patchsize==0.) CALL error_mess(ngroups(i)%locid,'no area information for stand: ',ngroups(i)%patchsize)
      IF(ngroups(i)%mhoe==0.) ngroups(i)%mhoe=wachsfunc(ngroups(i)%alter,hymax_Parthe,hb_Parthe,hT_Parthe)
      IF(ngroups(i)%dm==0.)   ngroups(i)%dm=wachsfunc(ngroups(i)%alter,dymax_Parthe,db_Parthe,dT_Parthe)
      IF(ngroups(i)%gf==0.AND.ngroups(i)%volume==0.AND.ngroups(i)%baumzahl==0.) ngroups(i)%gf=PI/4.*(ngroups(i)%dm/100.)**2*10000./(PI*(spar(ngroups(i)%taxid)%crown_a*ngroups(i)%dm+spar(ngroups(i)%taxid)%crown_b)**2)
   ELSEIF(ngroups(i)%schicht==50) THEN
      ! gap filling for plenterwald
      ! this routine is built on the use of the so called plenterwaldkurve (plenterwaldcurve)
      ! i.e. exponential decrease in number of trees in age classes
      init_plent=.false.
      IF(init_plent) THEN
   
        k_age=newton_plenter(0.15,low_age(ngroups(i)%taxid),high_age(ngroups(i)%taxid),dymax_Parthe,db_Parthe,dT_Parthe,ngroups(i)%dm)
        ngroups(i)%baumzahl=n0ofvol(k_age,low_age(ngroups(i)%taxid),high_age(ngroups(i)%taxid),dymax_Parthe,db_Parthe,dT_Parthe,hymax_Parthe,hb_Parthe,hT_Parthe,ngroups(i)%volume,formfactor)
        WRITE(8989,*) i,k_age,ngroups(i)%baumzahl,ngroups(i)%patchsize,ngroups(i)%baumzahl/ngroups(i)%patchsize,ngroups(i)%dm
      ELSE
        ngroups(i)%alter=inv_wachsfunc(ngroups(i)%dm,dymax_Parthe,db_Parthe,dT_Parthe)
        ngroups(i)%mhoe=wachsfunc(ngroups(i)%alter,hymax_Parthe,hb_Parthe,hT_Parthe)
        ngroups(i)%baumzahl=ngroups(i)%volume/(PI/4.*(ngroups(i)%dm/100.)**2*ngroups(i)%mhoe*formfactor)
        ngroups(i)%gf=PI/4.*(ngroups(i)%dm/100.)**2*ngroups(i)%baumzahl*10000./ngroups(i)%patchsize
        WRITE(8999,*) i,ngroups(i)%baumzahl,ngroups(i)%patchsize,ngroups(i)%gf,ngroups(i)%dm
      ENDIF
   ELSE
      CALL error_mess(ngroups(i)%locid,'unknown schicht_id occured: ',real(ngroups(i)%schicht))
   END IF  ! end of distinction according to layer (schicht)
END   ! subroutine data_gap_fill_DSW

FUNCTION newton_plenter(X,low_age,high_age,dmax,b,T,dg)
IMPLICIT NONE
REAL newton_plenter
REAL F,DF,X,DX,dmax,b,T,dg
INTEGER J,stepmax,low_age,high_age
! Newton-plenter is to be called with a start value for X 
! which is k_age here
! 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=50)
      DO 7 J=1,stepmax
        CALL fdfk(X,low_age,high_age,dmax,b,T,dg,F,DF)
      IF (J==stepmax) WRITE(8989,*) F, DF, X
!       IF(J==15) STOP
        IF(DF.EQ.0.0) THEN
          DX=0.01*X
        ELSE
          DX=F/DF
        ENDIF
        newton_plenter=X
        IF(DX.GT.X) DX=X/2.
        X=X-DX
        IF(ABS(DX).LT.0.0005) RETURN
7     END DO
END

SUBROUTINE fdfk(k_age,low_age,high_age,dmax,b,T,dg,F,DF)
! calculates function value and derivative for newton_plenter
USE data_par
USE data_simul
IMPLICIT NONE
INTEGER  :: low_age,high_age,age
REAL     :: term(1:4),sum(1:4),F,DF,k_age,dg,dmax,b,T,wachsfunc
   sum=0.
   DO age=low_age,high_age
      term(1)=exp(-k_age*age)
      term(2)=PI/4.*wachsfunc(age,dmax,b,T)**2*term(1)
      term(3)=-term(2)*age
      term(4)=-term(1)*age
      sum(1)=sum(1)+term(2)
      sum(2)=sum(2)+term(1)
      sum(3)=sum(3)+term(3)
      sum(4)=sum(4)+term(4)
   END DO
   F=(sum(1)/sum(2)*4/PI)**0.5-dg
   DF=((1./PI)**0.5*(sum(3)*sum(2)-sum(4)*sum(1)))/(sum(2)**2.*(sum(1)/sum(2))**0.5)
END    ! subroutine fdfk

FUNCTION n0ofvol(k_age,low_age,high_age,dmax,db,dT,hmax,hb,hT,vol,formfactor)
! calcualtes number of trees at dbh = 2cm for the plenter wald curve
! called by data_gap_fill_DSW if schicht=20 and init_plent=true; schicht is the word for layer
USE data_par
USE data_simul
IMPLICIT NONE
INTEGER  :: low_age,high_age,age,n0ofvol
REAL     :: sum,k_age,dmax,db,dT,hmax,hb,hT,vol,wachsfunc,formfactor
   sum=0.
   DO age=low_age,high_age
      sum=PI/4.*wachsfunc(age,dmax,db,dT)**2*exp(-k_age*age)*wachsfunc(age,hmax,hb,hT)
   END DO
   n0ofvol=NINT(vol/(sum/10000.*pi/4.*formfactor))
END    ! function n0ofvol
