Skip to content
Snippets Groups Projects
canopy.f 63.8 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 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutine canopy for:                     *!
!*    Calculation of canopy geometry & light absorption          *!
!*      with                                                     *!
!*     CALC_LA                                                   *!
!*     LIGHT_GROWTH                                              *!
!*     COV_AREA                                                  *!
!*     Light_1                                                   *!
!*     Light_2                                                   *!
!*     Light_3                                                   *!
!*     Light_4                                                   *!
!*     L_3_COH_LOOP                                              *!
!*     L_4_COH_LOOP                                              *!
!*    LIGHT_OUT_2                                                *!  
!*    CROWN_PROJ                                                 *! 
!*                                                               *!
!*                  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 CANOPY            *!
!**********************************!

SUBROUTINE CANOPY

  !*** Declaration part ***!

  USE data_out
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  integer i

  ! If no Cohorts on the patch, initialize properly
  IF( anz_coh == 0 ) THEN

      lowest_layer=0
      highest_layer=0
      vStruct%cumLAI= 0.
      vStruct%Irel  = 0.
      vStruct%sumBG = 0.
      Irelpool      = 0.
      BGpool        = 0.
      LAI           = 0.

      ! full light on the ground (layer = 0)
      ! Lightroutine 1,2
      vStruct(highest_layer)%Irel=1
      ! Lightroutine 3,4
      Irelpool(highest_layer)=1

      ! the whole patch is availabe for recruitment
      BGpool(highest_layer+1)=1
      BGpool(highest_layer+2)=1
      all_leaves_on=0

      ! Calculation of leaf area, lowest and highest layer, etc.
      ! for all cohorts in all respective layers
      CALL CALC_LA   ! leaf area etc. always calculate
      
      RETURN
  END IF

  ! Calculation of leaf area, lowest and highest layer, etc.
  ! for all cohorts in all respective layers
  CALL CALC_LA

  IF(flag_end.EQ.3) RETURN

   IF( flag_light == 1 )THEN

       CALL LIGHT_1

   ELSE IF ( flag_light == 2 ) THEN

       CALL LIGHT_2

   ELSE IF ( flag_light == 3 ) THEN

       CALL LIGHT_3

   ELSE IF ( flag_light == 4 ) THEN

       CALL LIGHT_4

   END IF
   DO i=1,anrspec
      ns = nrspec(i)
      IF(svar(ns)%act_sum_lai > svar(ns)%sum_lai) svar(ns)%sum_lai = svar(ns)%act_sum_lai
   ENDDO

  ! Determine relative light in the middle of each cohort canopy, the sla
  ! and the totFPAR per square meter patch and the total FPAR on the patch
  CALL LIGHT_GROWTH
  ! print relevant light parameters for the canopy for each layer and cohort
  if (time_out.gt.0 .and. out_flag_light.ne.0) CALL LIGHT_OUT_2

!------------------------------------------------
!------------------- SUBROUTINES ----------------
!------------------------------------------------
CONTAINS

SUBROUTINE CALC_LA

  ! Calculation of leaf area, lowest and highest layer, etc.
  ! for all cohorts in all respective layers

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: nl, i
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  ! auxiliary variable
  REAL    :: x    ! leaf area per crown unit [m**2/cm]

  vStruct%LA = 0.

  ! structure of the canopy is determined once at the start of the year
  ! initialisation
  IF(iday==1) THEN
      lowest_layer=250
      highest_layer=0
  END IF

  do i = 1, anrspec
     svar(nrspec(i))%act_sum_lai = 0.
  enddo

  p => pt%first

  DO WHILE (ASSOCIATED(p))
    ns = p%coh%species

    ! cohort loop for determination of lowest and highest canopy layer of the tree crown
    ! structure of the canopy must only be determined once at the start of the year
    IF(iday==1) THEN

       ! determine bottom of the crown in terms of  number of layers
       p%coh%botLayer = INT( p%coh%x_hbole / dz ) + 1

       ! determine top of the crown in terms of  number of layers
       IF (MODULO(p%coh%height,dz)==0.) THEN
           p%coh%topLayer = INT( p%coh%height /  dz )
       ELSE
           p%coh%topLayer = INT( p%coh%height /  dz ) + 1
       END IF

       ! remember the highest layer
       IF(p%coh%topLayer > highest_layer .AND. p%coh%toplayer < 250) THEN
          highest_layer=p%coh%topLayer

       ELSE IF(p%coh%toplayer >= 250) THEN
           if (.not.flag_mult8910) then
               CALL stop_mess(time,'FATAL EXCEPTION RAISED IN CANOPY CALC_LA')
               CALL error_mess(time,'maximal tree height of 125 m reached by cohort No.',REAL(p%coh%ident))
           endif  
           flag_end=3
           RETURN
       END IF

       !remember the lowest layer of the stand
       IF(p%coh%botLayer < lowest_layer) THEN
           lowest_layer=p%coh%botLayer
       END IF
    END IF

    p%coh%leafarea = 0.

    ! total leaf area of a tree in this cohort [m**2]
    IF((iday >= p%coh%day_bb) .AND. (iday <= spar(ns)%end_bb)) THEN
       p%coh%t_leaf = p%coh%med_sla * p%coh%x_fol

       ! amount of leaf area per tree in layers
       IF (p%coh%topLayer-p%coh%botLayer.GE.1) THEN
          ! now calculate leaf area per crown unit of this tree [m**2/cm]
          x = p%coh%t_leaf / ( p%coh%height - p%coh%x_hbole )
          p%coh%leafArea( p%coh%botLayer ) = ( dz - MODULO( p%coh%x_hbole, dz ) ) * x

          IF (MODULO(p%coh%height,dz)==0.) THEN
             p%coh%leafArea( p%coh%topLayer ) = dz * x
          ELSE
             p%coh%leafArea( p%coh%topLayer ) = MODULO( p%coh%height, dz ) * x
          END IF

          DO nl = p%coh%botLayer+1, p%coh%topLayer-1
             p%coh%leafArea(nl) = x * dz
          END DO
       ELSE
          p%coh%leafArea(p%coh%botLayer) = p%coh%t_leaf
       END IF
       ! Update vertical patch leaf area profile of the canopy
       DO nl = p%coh%botLayer, p%coh%topLayer
         vStruct(nl)%LA = vStruct(nl)%LA + p%coh%leafArea(nl) * p%coh%nTreeA
       END DO
    ELSE
        p%coh%leafArea=0.
    ENDIF
      IF(iday<=spar(ns)%end_bb) svar(ns)%act_sum_lai = svar(ns)%act_sum_lai + p%coh%ntreea*p%coh%t_leaf/kpatchsize
    p => p%next
  END DO
END SUBROUTINE CALC_LA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_GROWTH

  ! Determine  relative light in the middle of each cohort canopy, the sla,
  ! the total FPAR on the patch

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

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

  totFPARsum=0    ! sum of all totFPAR's
  totFPARcan=0    ! sum of all totFPAR's for the canopy
  p => pt%first

  DO WHILE (ASSOCIATED(p))

      ns=p%coh%species

      ! the new average specific leaf area per cohort depends
      ! on the light regime in the middle of the canopy
      ! this is the SLA which is used for the leaf area distr. in the next year
      ! the new average specific leaf area per cohort depends on the
      ! mean light regime in the middle in the canopy

      ! IrelCan modifies the growthfunction
      IF(all_leaves_on==1) THEN
         
         select case (flag_light)
         case (1,2)
            p%coh%med_sla = spar(ns)%psla_min+spar(ns)%psla_a*&
                    (1-(vStruct(p%coh%toplayer)%Irel+vStruct(p%coh%botlayer)%Irel)/2.)
            p%coh%IrelCan = vStruct(p%coh%toplayer)%Irel
         case default
            p%coh%med_sla = spar(ns)%psla_min+spar(ns)%psla_a*&
                           (1-(p%coh%Irel(p%coh%topLayer)+p%coh%Irel(p%coh%botLayer))/2.)
            select case (ns)
            case (10)   ! Douglas fir
                help = p%coh%botLayer+2*(p%coh%toplayer - p%coh%botLayer) / 3
                p%coh%IrelCan = p%coh%Irel(help)
            case default
                help = vStruct(p%coh%toplayer)%SumBG
                if (help .gt. 0.) then
                    p%coh%IrelCan = p%coh%Irel(p%coh%toplayer)*MIN(kpatchsize/help, 1.)
                else
                    p%coh%IrelCan = p%coh%Irel(p%coh%toplayer)
               endif
            end select    ! ns 
         end select    ! flag_light
      END IF

      totFPARsum = totFPARsum + p%coh%totFPAR*p%coh%nTreeA
      IF (p%coh%species .le. nspec_tree .or. p%coh%species.eq.nspec_tree+2) totFPARcan = totFPARcan + p%coh%totFPAR*p%coh%nTreeA
      p => p%next
  END DO

END SUBROUTINE LIGHT_GROWTH

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE COV_AREA

  ! calculate coverage-area as fraction of the patchsize per tree and layer

  !*** Declaration part ***!
  USE data_climate
  USE data_par
  USE data_stand
  USE data_site

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  ! Variables to test restriction in light model 4
  REAL  :: y         ! potential shadow cast of the cohort [m]
  REAL  :: w         ! effective shadow cast of the cohort [m]
  REAL  :: l         ! side length of a coort layer [m]
  REAL  :: reqarea   ! area of the patch required for the shadow cast for all cohorts per layer
  INTEGER :: layer_flag ! remember the highest layer where first LM4 restriction occurs
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  y = dz/100/TAN(beta)
  lm3layer=0
  layer_flag=0

  DO i = highest_layer, lowest_layer, -1
      reqarea=0.
      p => pt%first
      DO WHILE (ASSOCIATED(p))

          p%coh%BG(i) = 0.
          ! only those trees that have leaves
          IF((iday >= p%coh%day_bb) .AND. (iday <= spar(p%coh%species)%end_bb) .AND. &
             i <= p%coh%topLayer .AND. i >= p%coh%botLayer) THEN

             IF (vStruct(i)%sumBG > kpatchsize) THEN

                 p%coh%BG(i)=p%coh%crown_area/vStruct(i)%sumBG

             ELSE

                 p%coh%BG(i)=p%coh%crown_area/kpatchsize

             END IF

             l = SQRT(p%coh%BG(i)*kpatchsize)

             reqarea = reqarea + l*y*p%coh%nTreeA

          END IF

         p => p%next

     END DO  ! cohorts

     IF( kpatchsize > vStruct(i)%sumBG .AND. reqarea /= 0) THEN

       w = y*(kpatchsize-vStruct(i)%sumBG)/reqarea

     ELSE

       w = 0

     END IF

     p => pt%first

     DO WHILE (ASSOCIATED(p) .AND. layer_flag.EQ.0)
         ! only those trees that have leaves
         IF((iday >= p%coh%day_bb) .AND. (iday <= spar(p%coh%species)%end_bb) .AND. &
             i <= p%coh%topLayer .AND. i >= p%coh%botLayer) THEN

             l = SQRT(p%coh%BG(i)*kpatchsize)
             ! layer from that on light model 3 is used instead of light model 4
             ! because of LM4 restrictions
             IF( y-w > w+l ) THEN
                 layer_flag=1
                 lm3layer = i
                 EXIT ! do loop
             END IF
         END IF
         p => p%next
     END DO

  END DO     ! layers

END SUBROUTINE COV_AREA

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_1

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i, nl

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

  ! auxiliary variables
  REAL    :: radSum       ! sum of absorbed radiation (help variable)
  REAL    :: pfext=0.6    ! extinction coefficient. Only for one specie.


  !*** Calculation part ***!
  ! Intialization radiation summator
  radSum         = 0.
  vStruct%cumLAI = 0.
  vStruct%Irel   = 0.
  ! Calculate cumulative leaf area index and absorbed radiation per layer
  ! using Lambert-Beer
  vStruct(highest_layer)%Irel=1

  DO i = highest_layer, lowest_layer, -1
    vStruct(i)%cumLAI  = vStruct(i)%LA/kPatchsize + vStruct(i+1)%cumLAI
    vStruct( i )%radFrac =  1. - Exp(-pfext * vStruct(i)%cumLAI) - radSum
    radSum = radSum + vStruct(i)%radFrac
    vStruct(i-1)%Irel=vStruct(i)%Irel-vStruct(i)%radFrac

  END DO

  ! Light intensitiy unto the ground
  DO i = lowest_layer - 2, 0, -1
      vStruct(i)%Irel=vStruct(i+1)%Irel
  END DO

  ! total LAI is simply the value of cumLAI at the forest floor
  LAI = vStruct(lowest_layer)%cumLAI
  IF(lai>laimax) laimax=lai

  ! Determine layer-specific & total fraction of PAR absorbed by this tree
  p => pt%first
  DO WHILE (ASSOCIATED(p))

    p%coh%totFPAR = 0.
    p%coh%FPAR = 0.

    DO nl = p%coh%botLayer, p%coh%topLayer

        p%coh%FPAR(nl) = p%coh%leafArea(nl) / vStruct(nl)%LA * vStruct(nl)%radFrac
        p%coh%totFPAR  = p%coh%totFPAR + p%coh%FPAR(nl)

    END DO
    p => p%next
  END DO

   IF(all_leaves_on==1) THEN

      p => pt%first

      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=p%coh%FPAR(i)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

END SUBROUTINE LIGHT_1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_2

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  real    :: help

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

  !*** Calculation part ***!
  vStruct%cumLAI = 0.
  vStruct%Irel   = 0.

   ! cohort loop
   p => pt%first
   DO WHILE (ASSOCIATED(p))

       p%coh%FPAR = 0.
       p%coh%totFPAR = 0.

       p => p%next
   END DO        ! cohort loop

  ! Now calculate crown projection per tree and layer and
  ! the coverage sum over all layers
  CALL CROWN_PROJ

  ! now calculate coverage-area as fraction of the patchsize per tree and layer
  CALL COV_AREA

  vStruct(highest_layer)%Irel=1

  DO i = highest_layer, lowest_layer, -1

    p => pt%first

    help=0.

    vStruct(i)%cumLAI  = vStruct(i)%LA/kpatchsize + vStruct(i+1)%cumLAI

    DO WHILE (ASSOCIATED(p))

        ns=p%coh%species

         IF (p%coh%BG(i).ne.0.) THEN

             ! faction of absorbed light rel. to the light at the top of this layer
             ! the reference area is the whole patch (weighted by BG(i))!
             p%coh%FPAR(i)=(1-exp(-spar(ns)%pfext*p%coh%leafArea(i)/&
                                          kpatchsize/p%coh%BG(i)))*p%coh%BG(i)

             ! sum up the total absorbed fraction of this cohort,
             ! the total fraction of absorbed light in this layer
             ! is the fraction absorbed* fraction of light*BG
             ! the reference area is the whole patch!
             p%coh%totFPAR=p%coh%totFPAR+vStruct(i)%Irel*p%coh%FPAR(i)*&
                                (1+(0.5-vStruct(i)%Irel)*spar(ns)%fpar_mod/0.5)

             ! at first sum all the absorbed light fractions over the cohorts
             help=help+p%coh%FPAR(i)*p%coh%nTreeA

         ELSE

             p%coh%FPAR(i)=0.

         END IF

        p => p%next

    END DO

    ! then calculate the fraction of light which is available for the next layer
    vStruct(i-1)%Irel=vStruct(i)%Irel*(1-help)

  END DO

  ! Light intensitiy unto the ground
  DO i = lowest_layer - 2, 0, -1
      vStruct(i)%Irel=vStruct(i+1)%Irel
  END DO

   IF(all_leaves_on==1) THEN

      p => pt%first

      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=vStruct(i)%Irel*p%coh%FPAR(i)*(1+(0.5-vStruct(i)%Irel)*spar(ns)%fpar_mod/0.5)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

  ! total LAI is simply the value of cumLAI at the forest floor
  LAI = vStruct(lowest_layer)%cumLAI
  IF(lai>laimax) laimax=lai

END SUBROUTINE LIGHT_2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE L_3_COH_LOOP(i,j)

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

   INTEGER :: i, j   ! i= Schicht, j= Variante
   REAL    :: help

   p => pt%first

   ! cohort loop in layer i
   DO WHILE (ASSOCIATED(p))

       ns=p%coh%species

       IF((iday < p%coh%day_bb) .OR. (iday > spar(ns)%end_bb))  GOTO 1313

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           p%coh%FPAR(i)=1-exp(-spar(ns)%pfext*p%coh%leafArea(i)/&
                               kpatchsize/p%coh%BG(i))

           ! FPAR is related to the projection area and has to be modified
           ! by the same factor by that the projection area is being modified
           ! in case sumBG > patchsize
           p%coh%FPAR(i)=p%coh%FPAR(i)*MIN(kpatchsize/vStruct(i)%sumBG,1.)

           ! test wether the cohort is new, was there before or will not be
           ! represented in the next layer
           IF (i == p%coh%toplayer) THEN

               p%coh%Irel(i)=Irelpool(i)

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

           ELSE IF (i == p%coh%botlayer) THEN

               IF( j == 2 ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

               END IF

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

               ! The light which leaves the cohort is fed into the pool
               ! the light intensitiy is weighted by the overall BG of this cohort
               Irelpool(i-1)=(1/(p%coh%BG(i)*p%coh%nTreeA+BGpool(i)))*&
                             (p%coh%BG(i)*p%coh%nTreeA*p%coh%Irel(i-1)+BGpool(i)*Irelpool(i-1))

               ! BG of the pool available for the next layer increases
               BGpool(i)=BGpool(i)+p%coh%BG(i)*p%coh%nTreeA

           ELSE

               IF( j == 2 ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

               END IF

               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               p%coh%totFPAR=p%coh%totFPAR+p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)

               ! light available for this cohort in the next layer
               p%coh%Irel(i-1)=p%coh%Irel(i)*(1-p%coh%FPAR(i))

           END IF

       END IF   ! Layer test

1313   CONTINUE

       p => p%next

   END DO        ! cohort loop

END SUBROUTINE L_3_COH_LOOP

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE LIGHT_3

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  INTEGER :: i
  REAL    :: help

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

  !*** Calculation part ***!
  vStruct%cumLAI = 0.
  Irelpool      = 0.
  BGpool        = 0.
  vStruct%Irel = 0.       ! test variable for the light balance in layers
  vStruct%radFrac = 0.    ! test variable for the light balance in layers

  ! cohort loop
  p => pt%first
  DO WHILE (ASSOCIATED(p))

      p%coh%FPAR = 0.
      p%coh%totFPAR = 0.
      p%coh%Irel = 0.

      p => p%next
  END DO        ! cohort loop

    ! Now calculate crown projection per tree and layer and
    ! the coverage sum over all layers
    CALL CROWN_PROJ

    ! now calculate coverage-area as fraction of the patchsize per tree and layer
    CALL COV_AREA

    ! -----------------------------------------------------------
    ! now calculate per tree and layer the effective LAI
    ! this gives the absorbed light per tree and layer
    ! this gives the total fraction absorbes light per tree
    ! further each tree and each layer has an individual light regime. The area
    ! which is not covered by trees is treated as a pool
    !
    ! reference area for the total fracation absorbed is the patch area

    ! above the canopy there is 100 % rel. light
    Irelpool(highest_layer)=1.
    ! the size of the pool is defined as the fraction of the patch
    ! which can potentially be used by new cohorts in the next layer.
    ! Therefore is is the patch-fraction which is free anyway plus the
    ! fraction coverd by cohorts that will not be present in the next layer
    ! this means, the light intensity Irelpool(i) is available on the
    ! area BGpool(i+1)
    BGpool(highest_layer+1)=1.

    DO i = highest_layer, lowest_layer, -1

        vStruct(i)%cumLAI  = vStruct(i)%LA/kpatchsize + vStruct(i+1)%cumLAI

        ! two cases:
        ! first case: sumBG increases in this layer or remains the same
        IF (vStruct(i+1)%sumBG<=vStruct(i)%sumBG) THEN

            ! three subcases:
            ! first subcase of 'sumBG increases': sumBG stays below patchsize
            ! ( no BG modification) or does not change
            IF ((vStruct(i+1)%sumBG.LT.kpatchsize.AND.vStruct(i)%sumBG.LE.kpatchsize).OR.&
                vStruct(i+1)%sumBG == vStruct(i)%sumBG) THEN

                ! At the beginning the light intensity of the pool remains the same
                ! but it will be updated when cohorts drop out
                Irelpool(i-1)=Irelpool(i)

                ! until there are cohorts dropping out
                BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)

                CALL L_3_COH_LOOP(i,1)

            ! second and third subcase of 'sumBG increases or remains the same'
            ! the BG's of the cohorts change because sumBG exceeds patchsize.
            ! second subcase: sumBG was < patchsize before
            ! third subcase: sumBG was > patchsize before
            ELSE

                ! BG and light intensitiy of the pool for the next(!) layer
                ! is 0 as long as there are no cohorts dropping out
                Irelpool(i-1)=0.
                BGpool(i)=0.

                p => pt%first

                ! cohort loop 1
                DO WHILE (ASSOCIATED(p))

                    ! calculate the new fraction covered by the pool
                    ! which is the old pool plus the fractions which are lost
                    ! by the old cohorts due to new BG's
                    ! this also changes the light intensity of the pool
                    ! This pool will all be used by the new cohorts
                    ! consider only cohorts that have been there before (i<toplayer)
                   IF (i<p%coh%toplayer.AND.i>=p%coh%botlayer .AND.&
                       iday >= p%coh%day_bb .AND. iday <= spar(p%coh%species)%end_bb) THEN

                        help=BGpool(i+1)+(p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA

                        Irelpool(i)=(1/help)*(Irelpool(i)*BGpool(i+1)+p%coh%Irel(i)*&
                                              (p%coh%BG(i+1)-p%coh%BG(i))*p%coh%nTreeA)

                        BGpool(i+1)=help

                    END IF       ! layer test

                    p => p%next

                END DO   ! cohort loop1

                CALL L_3_COH_LOOP(i,1)

            END IF    ! subcases of 'sumBG increases

        ! second case: sumBG decreases
        ELSE

            ! two subcases
            ! first subcase of 'sumBG decrease': sumBG < patchsize before and after
            ! i.e. BG's do not change
            ! i.e. all projection area requirements can be fulfilled in the next layer
            IF (vStruct(i+1)%sumBG.LT.kpatchsize) THEN

                ! At the beginning the light intensity of the pool remains the same
                ! but it will be updated when cohorts drop out
                Irelpool(i-1)=Irelpool(i)

                ! until there are cohorts dropping out
                BGpool(i)=(kpatchsize-vStruct(i)%sumBG)/kpatchsize

                CALL L_3_COH_LOOP(i,1)

            ! second subcase of 'sumBG decrease': sumBG remains > patchsize or
            ! sumBG was > patchsize, i.e. BG's do change
            ELSE

                ! BG of the pool for the next layer as long as there are
                ! no cohorts dropping out
                BGpool(i)=MAX((kpatchsize-vStruct(i)%sumBG)/kpatchsize,0.)
                Irelpool(i-1)=Irelpool(i)

                CALL L_3_COH_LOOP(i,2)

            END IF            ! subcases

        END IF        ! three main cases

    END DO    ! end layer loop
    ! -----------------------------------------------------------
   IF(all_leaves_on==1) THEN

      p => pt%first
      DO WHILE (ASSOCIATED(p))

          DO i = highest_layer, lowest_layer, -1
             p%coh%antFPAR(i)=p%coh%Irel(i)*p%coh%FPAR(i)*p%coh%BG(i)/p%coh%totFPAR
             p%coh%sleafarea(i)=p%coh%leafarea(i)
          END DO    ! end layer loop

      p => p%next
      END DO ! cohort loop
   ENDIF

    ! total LAI is simply the value of cumLAI at the lowest layer
    LAI = vStruct(lowest_layer)%cumLAI
    IF(lai>laimax) laimax=lai

    ! light intensitiy and free patch space unto the ground
    DO i = lowest_layer - 2, 0, -1
        Irelpool(i)=Irelpool(i+1)
        BGpool(i+1)=BGpool(i+2)
    END DO

END SUBROUTINE LIGHT_3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE L_4_COH_LOOP(i,j,beta,y)

  !*** Declaration part ***!
  USE data_species
  USE data_simul
  USE data_stand

  IMPLICIT NONE

  ! variables required for technical reasons
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

   INTEGER :: i, j     ! i= layer, j= type
   REAL    :: y        ! potential shadow cast of a cohort layer [m]
   REAL    :: l        ! side length of a cohort layer [m]
   REAL    :: w        ! effective shadow cast of a cohort layer [m]
   REAL    :: helplai  ! LAI per layer and cohort
   REAL    :: help
   REAL    :: beta     ! sun inclination
   REAL    :: dropoutpool  ! relative area covered by cohort dropping out
   REAL    :: f1,f2,f3,f4,f5,f6,f7,f8 ! average fraction of absorbed radiation in different
                                      ! regions of the tree layer according to the 4C description paper
   REAL    :: k        ! extintion coefficient
   REAL    :: reqarea  ! area of the patch required for the shadow cast for all cohorts per layer

   reqarea=0.

   ! cohort loop
   p => pt%first
   DO WHILE (ASSOCIATED(p))

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           l = SQRT(p%coh%BG(i)*kpatchsize)

           reqarea = reqarea + l*y*p%coh%nTreeA

       END IF

       p => p%next
   END DO        ! cohort loop

   ! the size of the pool is defined as the fraction of
   ! the patch which is not covered by cohorts. This is the
   ! area covered by the sum of the 'shadows' of the cohorts,
   ! i.e. y's or rather w's + the area of cohorts dropping out in the next layer +
   ! the are that exeeds the maximal required area by the shadow-cast.
   ! This is updated in each layer
   ! w is the width of the shadow-cast of the cohorts that is maximal y.
   ! This maximal y also defines the maximal required area for all shadows
   ! 'reqarea' = required area
   ! When the maximal y cannot be satisfied, then this area is reduced by the
   ! relative share of the available space not covered by cohorts to the
   ! maximal required area for shadow cast

   IF( kpatchsize > vStruct(i)%sumBG ) THEN
     if (reqarea .gt. 1E-08) then
       w = y*(kpatchsize-vStruct(i)%sumBG)/reqarea
     else
       w = y*kpatchsize
     endif
   ELSE

       w = 0

   END IF

   BGpool(i)=0.
   dropoutpool=0

   p => pt%first

   ! cohort loop in layer i
   DO WHILE (ASSOCIATED(p))

       ns=p%coh%species

       IF((iday < p%coh%day_bb) .OR. (iday > spar(ns)%end_bb))  GOTO 1313

       k = spar(ns)%pfext

       IF (i<=p%coh%toplayer.AND.i>=p%coh%botlayer) THEN

           l = SQRT(p%coh%BG(i)*kpatchsize)

             if( p%coh%BG(i).ne.0) then

		         helplai=p%coh%leafArea(i)/kpatchsize/p%coh%BG(i)
		         if (helplai .le. 0.) then
		         continue
		         endif
		     else
			     helplai = 0.
			 end if

           IF (i == p%coh%toplayer) THEN

               p%coh%Irel(i)=Irelpool(i)

           ELSE IF( j == 2 .AND. i /= p%coh%toplayer ) THEN

                   help=p%coh%BG(i)-p%coh%BG(i+1)

                   p%coh%Irel(i)=(1/(p%coh%BG(i)))*&
                                 (p%coh%Irel(i)*p%coh%BG(i+1)+Irelpool(i)*help)

           END IF

           ! two main cases:
           ! first case : all light from the side comes from the pool
           ! second case : light from the side comes partially from the cohort itself

           IF( w >= y ) THEN

               ! subcases : 1.: light from the side of the layer
               ! does only leave at the bottom of the layer
               ! 2: light from the side does also leave on the other side
               ! totFPAR per patch! Since the projection area changes totFPAR has to
               ! be related to the patch in each layer
               IF( y <= l ) THEN

                   f1 = 1-exp(-k*helplai/SIN(beta))
                   if (helplai .lt. 1.E-6) then
                      f2 = 0.  
                   else
                      f2 = 1-SIN(beta)/(k*helplai)*f1
                      if (f2 .lt. 0.) then
                           continue
                           f2 = 0.   
                      endif
                   endif
                   p%coh%totFPAR=p%coh%totFPAR+(1/kpatchsize)*&
                                 ((l-y)*l*p%coh%Irel(i)*f1+& ! max. LAI
                                  ! exits layer at the side
                                  y*l*f2*p%coh%Irel(i)+&
                                  ! from the side to the next layer
                                  y*l*f2*Irelpool(i))

                   p%coh%FPAR(i)=p%coh%totFPAR