Skip to content
Snippets Groups Projects
prepstand.f 31.1 KiB
Newer Older
Petra Lasch-Born's avatar
Petra Lasch-Born committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886
!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                Stand initialisation                           *!
!*                                                               *!
!*          CONTAINS SUBROUTINES :                               *!
!*              PREPARE_STAND                                    *!
!*              internal subroutines:                            *!
!*                  SLA_INI                                      *!
!*                                                               *!
!*              CALC_INT                                         *!
!*              CALC_WEIBLA                                      *!
!*              READ_STAND (treeunit)                            *!
!*              COH_INITIAL (coh)                                *!
!*              CREATE_MISTLETOE                                 *!
!*              CREATE_SOILVEG                                   *!
!*                                                               *!
!*              CONTAiNS FUNCTIONS :                             *!
!*              SURVAGE                                          *!
!*                                                               *!
!*                  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 prepare_stand

  USE data_simul
  USE data_site
  USE data_stand
  USE data_species
  use data_climate
  use data_par
  USE data_manag

  IMPLICIT NONE

  CHARACTER      :: a, with_storage
  CHARACTER(30)  :: text
  CHARACTER(50)  :: test_stand_id
  INTEGER        :: ios,treeunit
  LOGICAL        :: exs, lstin
  INTEGER        :: help_ip, test_vf
  REAL           :: test_patchsize, xx


  REAL help_height_top  ! auxiliary var. for setting mistletoe height at uppermost crown layer
  INTEGER which_cohort
  INTEGER nr_infect_trees
  INTEGER nr_mist_per_tree
  INTEGER i
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  IF(site_nr==1) THEN
      help_ip=site_nr
  ELSE
      help_ip=ip
  END IF

  pt = neu()
  anz_coh=0
  max_coh=0
  ios = -1
  nr_mist_per_tree=0

  IF(flag_stand>0) then

    exs = .false.
    stand_id = standid(help_ip)
    ! reading stand information from treefile
    inquire (File = treefile(help_ip), exist = exs)
    IF((exs .eqv. .false.) .or. (flag_stand==2)) then
      IF(exs .eqv. .false.) write(*,*) '  Stand initialization file not exists!'
      IF(flag_stand==2) write(*,*)'  Stand initialization with new file'
      write(*,'(A)',advance='no') '  Creating new file (y/n): '
      READ *, a
      IF(a.eq.'y'.or. a.eq.'Y') CALL initia

    ! planting of small trees
      if(flag_reg.eq.20) then
         call planting
         flag_reg=100
      end if
      flag_stand=1
      exs=.true.
    ENDIF
    ! read values from treefile
    IF (exs.eqv. .true.)  then
       treeunit=getunit()
       OPEN(treeunit,file=treefile(help_ip),action='read', pad='YES')  !!!! NSC implementation for ini-file
       READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage
       backspace treeunit
       if(with_storage .eq. 'C') then  !flag_dis .eq. 2 .and. 
          READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage,test_vf, test_patchsize
          call read_stand_with_nsc (treeunit)
          CLOSE(treeunit)
          kpatchsize = test_patchsize
          anz_coh       = max_coh
          coh_ident_max = anz_coh
        else
          READ(treeunit,'(I1,F12.0)',iostat=ios) test_vf, test_patchsize
  !       write(8888,*) ip, test_vf, flag_volfunc
          if(flag_multi.ne.4 .or. (flag_multi.eq.4.and.ip.eq.1) .or. (flag_multi.eq.8.and.ip.eq.1)) then
            IF(test_vf.NE.flag_volfunc) THEN
              if (.not.flag_mult8910) then
                CALL error_mess(time,"volume function in sim-file and the one used for initialisation do not match",REAL(flag_volfunc))
                CALL error_mess(time,"volume function (flag_volfunc) is set to",REAL(test_vf))
              endif
             flag_volfunc = test_vf
            ENDIF
          endif
         IF(test_patchsize .GT. 0.) THEN
           lmulti = .FALSE.
           IF(test_patchsize.NE.kpatchsize) THEN
             if (.not.flag_mult8910) then
               CALL error_mess(time,"patch size in sim-file and the one used for initialisation do not match",kpatchsize)
               CALL error_mess(time,"value in ini-file",test_patchsize)
               CALL error_mess(time,"value in sim-file",kpatchsize)
             endif
           kpatchsize = test_patchsize
           ENDIF
          ELSE
            lmulti = .TRUE.
         ENDIF
          do
           READ(treeunit,'(A)',iostat=ios) a
            IF (a .ne. '!') exit
          end do
          backspace treeunit
          ! generation of mistletoe cohort; mistletoe cohort need to be generated BEFORE tree cohorts as otherwise the light model becomes messy
          if (flag_dis.eq.1) then
            do i= 1, dis_row_nr
             if (dis_type(i) .eq. 'M') then
               if (flag_mistle.eq.0) then        !set #of mist. only once
                print *,"!! Note, implementation of mistletoe is restricted to trees of Pinus sylvestris"
                nr_mist_per_tree = dis_rel(i)
                flag_mistle=1                          ! flag indicating mistletoes
                call create_mistletoe  ! initialisation of Mistletoe
               endif
              anz_coh = max_coh
             endif
            enddo
          endif
       
          lstin = .TRUE.
          if(flag_multi.eq.4 .or. flag_multi.eq.8) stand_id = standid(help_ip)
            do while (lstin)
              IF (lmulti) THEN
               read(treeunit,*,iostat=ios) test_stand_id,  test_patchsize,text
               IF (ios .lt. 0) then
                if (.not.flag_mult8910) then
                   CALL error_mess(time,"stand identificator not found"//adjustl(stand_id)//"ip No.",real(help_ip))
                   write (*,*) '*** PREPSTAND:  program aborted'
                   write (*,*) '                stand identificator',stand_id,'  not found'
                   write (*,'(A, 2x,A)') '                 in initialisation file',treefile(help_ip)
                endif
               flag_end = 2
               return
               ENDIF
               IF (test_stand_id .ne. stand_id) THEN
                 read (treeunit,*) xx
                   do while (xx .gt. -90.0)
                    read (treeunit,*) xx
                   enddo  ! xx
                ELSE
                 lstin = .FALSE.
                 kpatchsize = test_patchsize
                 call read_stand (treeunit)
               END IF ! stand_id
              ELSE
               lstin = .FALSE.
               call read_stand (treeunit)
              END IF    ! lmulti
            end do      ! lstin
       CLOSE(treeunit)
       anz_coh       = max_coh
       coh_ident_max = anz_coh
       endif !w/o storage
    ENDIF    !exs.eqv. .true.
  END IF     !if stand >0

!if treefile not exists and not created:
IF(ios .ne. 0 .or. exs .eqv. .false.)THEN
  if (.not.flag_mult8910) PRINT *,' >>> No Stand Initialization possible '
  flag_stand=0
END IF

! Setting of height and number of mistletoe
if (flag_mistle.ne.0) then
    help_height_top=1.
    p=>pt%first
    DO WHILE (ASSOCIATED(p))
        if (p%coh%species.eq.3 .AND. p%coh%height.gt.help_height_top) then  !only on Pinus
            help_height_top=p%coh%height
            which_cohort=p%coh%ident
            nr_infect_trees=p%coh%nTreeA
        end if
        p=>p%next
    end do

    p=>pt%first
    DO WHILE (ASSOCIATED(p))
        if (p%coh%species.eq.nspec_tree+2) then
            p%coh%height  = help_height_top               !upper crown
            p%coh%x_hbole = p%coh%height-50.              !lower crown
            p%coh%nTreeA = nr_infect_trees*nr_mist_per_tree             !number of mistletoes
        end if
        if (p%coh%ident.eq.which_cohort) then             !mark uppermost tree cohort with flag mistletoe
            p%coh%mistletoe=1
        end if
        p=>p%next
    end do
end if ! end set height/number of mistletoe

! Soil Vegetation
if (flag_sveg .gt. 0) then
   call create_soilveg  ! initialisation of ground vegetation
   anz_coh = max_coh
endif

IF(flag_stand>0) CALL sla_ini
IF(flag_stand>0) CALL stand_bal_spec
CALL calc_int
CALL calc_weibla
if(flag_mg.ne.33) call overstorey

contains

SUBROUTINE sla_ini

 USE data_stand
 USE data_species

 IMPLICIT NONE
 TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

 p => pt%first

 DO WHILE (ASSOCIATED(p))
     ns=p%coh%species
     p%coh%med_sla=spar(ns)%psla_min+spar(ns)%psla_a*0.5
     p%coh%t_leaf = p%coh%med_sla * p%coh%x_fol
     p =>p%next
 END DO
end subroutine sla_ini

end subroutine prepare_stand

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

subroutine calc_int  !  calculation of intrinsic mortality rate

use data_species
implicit none
INTEGER j

do j=1,nspecies
spar(j)%intr = -log(0.01)/spar(j)%max_age
end do
end subroutine calc_int

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

subroutine calc_weibla
!   calculation of parameter lamda for Weibull-distribution of sress mortality

use data_species
implicit none
INTEGER j
REAL survage

do j=1,nspecies
spar(j)%weibla = -log(0.01)/(survage(j)**weibal)
end do

end subroutine calc_weibla

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

REAL function survage(ispec)
!  calculation of survival time per species depending on shade tolerance class stol

use data_species
implicit none
INTEGER :: ispec

IF(spar(ispec)%stol.eq.1) survage=20.
IF (spar(ispec)%stol.eq.2) survage=40.
IF (spar(ispec)%stol.eq.3) survage=60.
IF (spar(ispec)%stol.eq.4) survage=80.
IF (spar(ispec)%stol.eq.5) survage=100.
end function

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

SUBROUTINE read_stand (treeunit)

!  Read of stand initialisation

  USE data_par
  USE data_simul
  USE data_species
  USE data_stand

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini
  REAL           :: hdquo   ! help variable for stress initilization
  INTEGER        :: ios,treeunit, loc, i
  logical        :: treegroup_decid
  integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)

              do
                call coh_initial (coh_ini)
                READ(treeunit,'(5f12.5,2f10.0,i7, f10.0,i7, f9.5, f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                    coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
                    coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam
                IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit

                coh_ini%nTreeD = 0.
                coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
                coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)

				!NSC Speicher initialisieren
				!hier neue version auslesen
                treegroup_decid = .False.
                do i = 1, 5
                  if (decidous(i) .eq. coh_ini%species) then
                  treegroup_decid = .True.
                  exit
                  endif
                end do
                If (treegroup_decid .eq. .True.) then
                coh_ini%x_nsc_sap = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
                coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
                coh_ini%x_nsc_tb = coh_ini%x_tb * decid_tb_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
                coh_ini%x_nsc_crt = coh_ini%x_crt * decid_crt_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
                endif
                If (treegroup_decid .eq. .False.) then
                coh_ini%x_nsc_sap = coh_ini%x_sap * conif_sap_allo * cpart
                coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
                coh_ini%x_nsc_tb = coh_ini%x_tb * conif_tb_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
                coh_ini%x_nsc_crt = coh_ini%x_crt * conif_crt_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
                endif
				
                coh_ini%ident = max_coh + 1
                coh_ini%Fmax = coh_ini%x_fol
                coh_ini%x_health = 0
                coh_ini%x_hsap = 0.
                ns  = coh_ini%species
                coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol     ! kg * mg/g --> g
                 if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
                  CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
                else
                  coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
                  coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
                end if

 ! Stress calculation
                IF (coh_ini%diam.ne. 0.) THEN
                  hdquo = coh_ini%height/ (coh_ini%diam*100)
                  IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
                    coh_ini%x_stress =  coh_ini%x_age/2
                  ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
                    coh_ini%x_stress =  coh_ini%x_age*3./7.
                  ELSE
                    coh_ini%x_stress = 0.
                  END IF
                ELSE
                  coh_ini%x_stress = 0.
                END IF  ! coh_ini

                coh_ini%x_stress = 0.
                coh_ini%nta = coh_ini%nTreeA

                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
                ELSE
                  allocate(zeig)
                  zeig%coh = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
                END IF
                max_coh = max_coh + 1
              enddo

END SUBROUTINE read_stand

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

SUBROUTINE read_stand_with_nsc (treeunit)

!  Read of stand initialisation

  USE data_par
  USE data_simul
  USE data_species
  USE data_stand


  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini
  REAL           :: hdquo   ! help variable for stress initilization
  INTEGER        :: ios,treeunit, loc, i
  logical        :: treegroup_decid
  integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)
  character      :: a

              do
               READ(treeunit,'(A)',iostat=ios) a
               IF (a .ne. '!') exit
              end do
              backspace treeunit
              do
                call coh_initial (coh_ini)
!                READ(treeunit,'(5f12.5,2f10.0,i7, f7.0,i7, 2f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                READ(treeunit,'(5f12.5,2f10.0,i7,f7.0,i7, 5f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                    coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
                    coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam, coh_ini%x_nsc_tb, coh_ini%x_nsc_crt, coh_ini%x_nsc_sap 

                IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit

                coh_ini%nTreeD = 0.
                coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
                coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)

				!hier nur den NSC Max-Speicher initialisieren
				!hier neue version auslesen
                treegroup_decid = .False.
                do i = 1, 5
                  if (decidous(i) .eq. coh_ini%species) then
                  treegroup_decid = .True.
                  exit
                  endif
                end do
                If (treegroup_decid .eq. .True.) then
                coh_ini%x_nsc_sap_max = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
                coh_ini%x_nsc_tb_max = coh_ini%x_tb * decid_tb_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_crt * decid_crt_allo * cpart
                endif
                If (treegroup_decid .eq. .False.) then
                coh_ini%x_nsc_sap_max = coh_ini%x_sap * conif_sap_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_tb * conif_tb_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_crt * conif_crt_allo * cpart
                endif
!      IF(coh_ini%species==3.and.coh_ini%height.le.900) then
!      coh_ini%x_stress  =6
!      ELSE IF(coh_ini%species==1.and.coh_ini%height.le.2500) then
!      coh_ini%x_stress  =25
!      ELSE
!      coh_ini%x_stress = 0
!      ENDIF
                coh_ini%ident = max_coh + 1
                coh_ini%Fmax = coh_ini%x_fol
                coh_ini%x_health = 0
                coh_ini%x_hsap = 0.
                ns  = coh_ini%species
                coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol     ! kg * mg/g --> g
            ! calculate diameter at breast height and pipe length by call of subroutine in partitio
                 if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
!                if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.137.) then
                  CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
                else
                  coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
                  coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
                end if

 ! Stress calculation
                IF (coh_ini%diam.ne. 0.) THEN
                  hdquo = coh_ini%height/ (coh_ini%diam*100)
                  IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
                    coh_ini%x_stress =  coh_ini%x_age/2
                  ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
                    coh_ini%x_stress =  coh_ini%x_age*3./7.
                  ELSE
                    coh_ini%x_stress = 0.
                  END IF
                ELSE
                  coh_ini%x_stress = 0.
                END IF  ! coh_ini
! provisorisch stress auf Null setzen
                coh_ini%x_stress = 0.

                coh_ini%nta = coh_ini%nTreeA

                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
                ELSE
                  allocate(zeig)
                  zeig%coh = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
                END IF
                max_coh = max_coh + 1
              enddo

END SUBROUTINE read_stand_with_nsc

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



SUBROUTINE coh_initial (coh_ini)

  USE data_simul
  USE data_soil
  USE data_stand
  USE data_species

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini

        coh_ini%nTreeA = 0.
        coh_ini%nTreeD = 0.
        coh_ini%nTreeM = 0.
        coh_ini%nTreet = 0.
        coh_ini%nta    = 0.
        coh_ini%mistletoe = 0

        coh_ini%x_age  = 0.
        coh_ini%x_fol  = 0.
        coh_ini%x_sap  = 0.
        coh_ini%x_frt  = 0.
        coh_ini%x_hrt  = 0.
        coh_ini%x_crt  = 0.
        coh_ini%x_tb   = 0.
        coh_ini%x_hsap = 0.
        coh_ini%x_hbole= 0.
        coh_ini%x_Ahb  = 0.

        coh_ini%x_stress = 0
        coh_ini%x_health = 0

        coh_ini%bes     = 0.
        coh_ini%med_sla = 0.
        coh_ini%Fmax    = 0
        coh_ini%totBio  = 0.
        coh_ini%Dbio    = 0.
        coh_ini%height  = 0.
        coh_ini%deltaB  = 0.
        coh_ini%dcrb    = 0.
        coh_ini%diam    = 0.
        coh_ini%assi    = 0.
        coh_ini%LUE     = 0.
        coh_ini%resp    = 0.
        coh_ini%netAss  = 0.
        coh_ini%NPP     = 0.
        coh_ini%weekNPP = 0.
        coh_ini%NPPpool = 0.
        coh_ini%t_Leaf  = 0.
        coh_ini%geff    = 0.
        coh_ini%Asapw   = 0.
        coh_ini%crown_area = 0.

        coh_ini%BG        = 0.
        coh_ini%leafArea  = 0.
        coh_ini%sleafArea = 0.
        coh_ini%FPAR      = 0.
        coh_ini%antFPAR   = 0.
        coh_ini%Irel      = 0.

        coh_ini%totFPAR  = 0
        coh_ini%IrelCan  = 0
        coh_ini%botLayer = 0
        coh_ini%topLayer = 0
        coh_ini%survp    = 0.
        coh_ini%rel_fol  = 0.
        coh_ini%gfol     = 0.
        coh_ini%gfrt     = 0.
        coh_ini%gsap     = 0.
        coh_ini%sfol     = 0.
        coh_ini%sfrt     = 0.
        coh_ini%ssap     = 0.
        coh_ini%grossass = 0.
        coh_ini%maintres = 0.
        coh_ini%respsap  = 0.
        coh_ini%respfrt  = 0.
        coh_ini%respbr   = 0.

        coh_ini%height_ini = 0.
        coh_ini%ca_ini     = 0.

        coh_ini%rel_dbh_cl = 0
        coh_ini%underst    = 0

        coh_ini%fol_inc     = 0.
        coh_ini%fol_inc_old = 0.
        coh_ini%bio_inc     = 0.
        coh_ini%stem_inc    = 0.
        coh_ini%frt_inc     = 0.
        coh_ini%notViable   = .FALSE.

        coh_ini%intcap   = 0.
        coh_ini%prel     = 0.
        coh_ini%interc   = 0.
        coh_ini%prelCan  = 0.
        coh_ini%interc_st= 0.
        coh_ini%aev_i    = 0.
        coh_ini%demand   = 0.
        coh_ini%supply   = 0.
        coh_ini%watuptc  = 0.
        coh_ini%gp       = 0.
        coh_ini%drIndd   = 0.
        coh_ini%drIndPS  = 0.
        coh_ini%drIndAl  = 0.
        coh_ini%nDaysGr  = 0
        coh_ini%isGrSDay = .false.

        coh_ini%litC_fol  = 0.
        coh_ini%litC_fold = 0.
        coh_ini%litN_fol  = 0.
        coh_ini%litN_fold = 0.
        coh_ini%litC_frt  = 0.
        coh_ini%litC_frtd = 0.
        coh_ini%litN_frt  = 0.
        coh_ini%litN_frtd = 0.
        coh_ini%litC_stem = 0.
        coh_ini%litN_stem = 0.
        coh_ini%litC_tb   = 0.
        coh_ini%litC_crt  = 0.
        coh_ini%litC_tbcd = 0.
        coh_ini%litN_tb   = 0.
        coh_ini%litN_crt  = 0.
        coh_ini%litN_tbcd = 0.
        coh_ini%Nuptc_c   = 0.
        coh_ini%Nuptc_d   = 0.
        coh_ini%Ndemc_d   = 0.
        coh_ini%RedNc     = 1.
        coh_ini%N_pool    = 0.
        coh_ini%N_fol     = 0.
        coh_ini%wat_mg    = 0.   ! soley forflag_wred=9

        coh_ini%nroot   = 0
        coh_ini%shelter = 0
		coh_ini%day_bb  = 0
        coh_ini%x_nsc_sap = 0.
        coh_ini%x_nsc_tb = 0.
        coh_ini%x_nsc_crt = 0.
        coh_ini%x_nsc_sap_max = 0.
        coh_ini%x_nsc_tb_max = 0.
        coh_ini%x_nsc_crt_max = 0.
		
      if (coh_ini%species .ne. nspec_tree+2) then  ! no root allocation for mistletoe
        allocate (coh_ini%frtrel(nlay))
        allocate (coh_ini%frtrelc(nlay))
        if (flag_wred .eq. 9) then
            allocate (coh_ini%rld(nlay))
            coh_ini%rld  = 0.
		endif
		allocate (coh_ini%rooteff(nlay))
	    coh_ini%frtrel  = 0.
		coh_ini%rooteff = 0.
      end if   ! end exclude mistletoe
END SUBROUTINE coh_initial
!*************************************************************************
SUBROUTINE create_mistletoe
  USE data_plant
  USE data_simul
  USE data_species
  USE data_stand
  USE data_climate
  USE data_soil
  USE data_species
  USE data_par
  IMPLICIT NONE
  TYPE(cohort)   :: coh_ini
  real            :: help_height_top, help_height_bot
  REAL, EXTERNAL  ::  fi_lf, dfi_lf, ddfi_lf

 ! initialising of cohort of mistletoe
 call coh_initial (coh_ini)
 ! set mistletoe here to 20 m height, will be changed after, when cohorts of trees will be initialised
 help_height_top=2000 
 help_height_bot=help_height_top-50 
 ! following values are from sample calcul. of 10 year old V.austr. from Pfiz 2010
 coh_ini%ident      = max_coh + 1
 coh_ini%species    = nspec_tree+2    ! Species = species after all tree species and ground veg.
 coh_ini%nTreeA     = 1               ! #of mistletoes, to be read-in in management file
 coh_ini%nTreeD     = 0               ! dead trees
 coh_ini%nta        = coh_ini%nTreeA  ! alive trees internal calc.
 coh_ini%x_age      = 10
 coh_ini%x_fol      = mistletoe_x_fol ! fol biomass per tree [kg DW/tree], 1 Viscum (10years) see Pfiz 2010
 coh_ini%x_sap      = 0.              ! set near-zero for partitioning
 coh_ini%x_frt      = 0.              ! set near-zero for partitioning
 coh_ini%height     = help_height_top ! highest_layer   ! highest_layer of all cohorts
 coh_ini%x_hbole    = help_height_bot !
 coh_ini%med_sla    = 0.              ! average cohort specific leaf area [m2/kg] is being calculated internal
 coh_ini%Fmax       = 0               ! anual change of leaf biomass, for now: now change
 coh_ini%crown_area = 0.0189          ! max. projected crown area (m2) per individuum, calculated from Pfiz 2010
 coh_ini%t_leaf     = coh_ini%med_sla* coh_ini%x_fol      !leaf area per tree [m2]  !
 coh_ini%day_bb     = 1               ! evergreen
! no partitioning of NPP into stem/leaf etc.
! no root allocation
 allocate(zeig)
  zeig%coh = coh_ini
  zeig%next => pt%first
  pt%first => zeig
 max_coh = max_coh + 1
END SUBROUTINE create_mistletoe

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

SUBROUTINE create_soilveg

!  Read of stand initialisation

  USE data_plant
  USE data_simul
  USE data_species
  USE data_stand
  USE data_climate
  USE data_soil

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini

  real            :: lai_help, irel_help, FRsum
  integer         :: age_stand, nr, j
  integer         :: flag_SV_allo, rnum
  real			  :: troot2

  REAL, EXTERNAL  ::  fi_lf, dfi_lf, ddfi_lf
  
  age_stand = 0
  lai_help = 0.
  irel_help = 0.
  call wclas(waldtyp)

	   zeig=>pt%first
       
	   DO WHILE (ASSOCIATED(zeig))
	     ns = zeig%coh%species
	     lai_help = lai_help + zeig%coh%ntreea*zeig%coh%x_fol* spar(ns)%psla_min
		  age_stand = MAX(zeig%coh%x_age,age_stand)
		 zeig=>zeig%next
       end do

   IF((flag_stand==0 .or. age_stand .le. 5) .AND. flag_sveg ==2) THEN
      NPP_est = 10.  
   ELSE if(age_stand.le.5) then
      if(ns.eq.4) then
	     NPP_est = 5
	   else
	      NPP_est = 10.
	   end if
   ELSE if(flag_reg.ne.0) then
         NPP_est = 10
   ELSE
      lai_help = lai_help/kpatchsize
      irel_help = exp(-0.5*lai_help)
      if( svar(nspec_tree+1)%RedN .lt.0.) then
	         NPP_est = irel_help * med_rad1 * 365./100. *0.5
	   else
             NPP_est = irel_help * med_rad1 * 365./100. *0.5 * svar(nspec_tree+1)%RedN
      end if
  ENDIF

     call coh_initial (coh_ini)

       coh_ini%species = nspec_tree+1    ! numbre of species determined automatically
         ns = coh_ini%species
         flag_SV_allo=1
  IF(flag_SV_allo==0) THEN
  ! the parameters pdiam in the species.par file are used for allocation fractions
         coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
		   coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
		   coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize    
  ELSE  
     FRsum=0.8*NPP_est/1000.    ! start value as fraction of NPP in kg DM m-2   
     CALL newt (FRsum, fi_lf, dfi_lf, ddfi_lf, 0.001, 100, rnum)
     IF(rnum==-1) THEN
        if (.not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort: ',real(ns))
         coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
		   coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
		   coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize    
     ELSE
         coh_ini%x_sap = (ksi*FRsum**kappa)*kpatchsize
		   coh_ini%x_fol = (FRsum/2.)*kpatchsize
		   coh_ini%x_frt = (FRsum/2.)*kpatchsize    
     ENDIF
  ENDIF 

         coh_ini%height  = 60.
         coh_ini%x_age   = 1
         coh_ini%nTreeA  = 1
         coh_ini%ident   = max_coh + 1
         coh_ini%Fmax    = coh_ini%x_fol
         coh_ini%med_sla = spar(coh_ini%species)%psla_min + spar(coh_ini%species)%psla_a*irel_help
         coh_ini%t_leaf  = coh_ini%med_sla* coh_ini%x_fol      ! [m2]

         coh_ini%nta     = coh_ini%nTreeA
         coh_ini%ca_ini  = kpatchsize
         coh_ini%day_bb  = 100            ! assumption budding on 8.April 

! root allocation
                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
				  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 = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
				  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
                max_coh = max_coh + 1

END SUBROUTINE create_soilveg

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

!***************************!
! FUNCTION fi_lf           *!
!***************************!

REAL FUNCTION fi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  fi_lf = spar(nspec_tree+1)%pss*ksi*x**kappa + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2.*x - NPP_est/1000.
END ! FUNCTION fi_lf

!***************************!
! FUNCTION dfi_lf          *!
!***************************!

REAL FUNCTION dfi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  dfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*x**(kappa-1.) + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2. 
END ! FUNCTION dfi_lf

!***************************!
! FUNCTION ddfi_lf         *!
!***************************!

REAL FUNCTION ddfi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  ddfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*(kappa-1.)*x**(kappa-2.) 
END ! FUNCTION ddfi_lf