InitAndSetConstHeightClProps.f90

Go to the documentation of this file.

!==============================================================================
!===============================================================
SUBROUTINE initandsetconstheightclprops(nspc)
    
    ! CALL modules for TYPE definitions and +- fixed variables>---------------------------------------------------------------
    use all_par, only:   tabsepoutput,         &
            maxhc, birht, hcwidth, htmax,          &
            pltsiz,  maxlat, maxlon,               &
            maxlc, toplt, ltfrq, ltcum, ltval, la, &
            spec, speccl, stategrid, seedrain 
    
    use filelistmodule, only:  lightandleafareaperlightclass,  biompertreehcloutfile, laipertreehcloutfile, &
            seedprodfractpertreehcloutfile, dbhpertreehcloutfile,  heightperhcloutfile   
    IMPLICIT NONE
    
    ! Passed variables>--------------------------------------------------------------------------------
    INTEGER, INTENT(in) :: nspc    ! Passed number of species                                   [=spec]
    
    ! Local variables>---------------------------------------------------------------------------------
    INTEGER ::  latIni, lonIni,   &  ! tmp vars to initialize grid with trees            
            isp,      &  ! looping variable for species               
            ltcl,htcl    ! looping variable for light and height classes               
    REAL    ::  &
            truht,    &  ! tmp var to calculate heights                               [=trueHeight]
            nexht,    &  ! tmp var to calculate heights                               [=nextHeight]
            diamt,    &  ! tmp var to calculate diameters               [in DetermineTree..  =diam]
            nextdm,    &  ! tmp var to calculate diameters                               [=nextDiam]
            bio , &          ! tmp var to calculate biomass                      [=BioMassAtHeight(..)]
            DBHFromHeight,&  ! function to calculate diameter from height        <- Allometries.f90                    
            LAFromDBH,&      ! function to calculate leaf area from  diameter    <- Allometries.f90                    [ 
            BiomassFromDBH   ! function to calculate biomass from  diameter      <- Allometries.f90       
    character*1 :: tab
    ! Here the SUBROUTINE starts>**********************************************************************
    if (tabsepoutput) then
        tab = char(9)
    else
        tab = ','
    end if
    
    !----- initialize the light frequencies. First no trees, i.e. full light, highest lightclass frequency 1
    
    ltfrq        = 0.0
    ltfrq(maxlc) = 1.0
    ltcum        = 1.0
    
    
    !-----  The light values and the LAI values belonging to them in the light classes are calculated
    byltcls1: do ltcl = 1, maxlc
        ltval(ltcl) = real(ltcl*(ltcl + 1))/real(maxlc*(maxlc + 1))
        la(ltcl)    = - log( max( ltval(ltcl)/toplt, 0.0000001 ) )/.25 ! inverse Beer-Lambert with extinction factor 0.25
    end do byltcls1
    
    byspcs1: do isp = 1, nspc  ! species loop
 
    !----- <Initialize the tree species in grid cells and height classes array [=DetermineTreeSizeVariables]>     
        latloop1: do latini = 1, maxlat   ! spatial loops
            lonloop1: do lonini = 1, maxlon ! spatial loops
                !----- Here the  state variable values for all height classes are initialized to zero
                byhtcls11: do htcl = 0, maxhc
                    stategrid(latini, lonini)%sp(isp)%numin(htcl) = 0.0
                end do byhtcls11
                !----- Here the  state variable values for the antagonists, the seedbank and the seedrain are initialized to zero           
                stategrid(latini, lonini)%sp(isp)%antagonist = 10.0
                stategrid(latini, lonini)%sp(isp)%seedBank = 0.0
                seedrain(latini, lonini)%sp(isp)%newSeeds = 0.0
            end do lonloop1
        end do latloop1
        
        !----- Here height class dependent but constant values are calculated: biomass, diameter, LA, height 
        nexht  = birht
        nextdm = dbhfromheight(birht, birht, spec(isp)%khmax,spec(isp)%maxdhm)
        byhtcls1: do htcl = 1, spec(isp)%maxhtc  !---- from bottom to top
            truht = nexht
            nexht = min((real(htcl)*hcwidth) + birht, real(spec(isp)%khmax)) ! height belonging to height class
            speccl(isp)%htmin(htcl) = truht                                  ! store height
            diamt = nextdm
            nextdm= dbhfromheight(nexht, birht, spec(isp)%khmax,spec(isp)%maxdhm) ! get diameter
            speccl(isp)%diain(htcl) = diamt                                   ! store diameter
            speccl(isp)%lain(htcl) = lafromdbh(diamt, spec(isp)%tb,spec(isp)%stn) 
            if (nextdm > diamt) then
                bio =  biomassfromdbh(diamt,  spec(isp)%tb, spec(isp)%stn )
                speccl(isp)%bioin(htcl) = bio * 10 / pltsiz
                !NZ: use factor 10000 instead of 10 if biomass is stem volume
                ! factor 10 transforms kg/m2 to t/ha, if you use biomass (normal version)
                ! factor 10000 transforms into dm3 per ha (NOT m3!!) if stem volume is used instead of biomass
                ! (see SUBROUTINEs BiomassFromDBH, and Report)
                
            else ! if at top, then the mean between this biomass and the  one at the next height (=diameter) class is taken
                bio = biomassfromdbh(diamt, spec(isp)%tb, spec(isp)%stn)  
                speccl(isp)%bioin(htcl) = bio*10/pltsiz
                bio = biomassfromdbh(nextdm, spec(isp)%tb, spec(isp)%stn) 
                speccl(isp)%bioin(htcl) = (speccl(isp)%bioin(htcl) + (bio*10/pltsiz))/2.0                
            end if
        end do byhtcls1
        
        !----- Here the  values for the lowest height class (seedlings up to 1.37 m) are set... 
        speccl(isp)%htmin(0) = 0.0   
        speccl(isp)%diain(0) = 0.0 
        speccl(isp)%lain(0)  = speccl(isp)%lain(1) /10.0  ! rough guess
        speccl(isp)%bioin(0) = speccl(isp)%bioin(1)/10.0 ! rough guess
        !-----  ...and the states initialized  
        latloop2: do latini = 1, maxlat
            lonloop2: do lonini = 1, maxlon
                stategrid(latini, lonini)%sp(isp)%numin(0) = 0.0
            end do lonloop2
        end do latloop2
        
        !----- Here the seed production is determined proportional to LA 
        byhtcls3: do htcl = 0, spec(isp)%maxhtc
            if (speccl(isp)%htmin(htcl) >= spec(isp)%minmat) then
                speccl(isp)%seedprod(htcl) = speccl(isp)%lain(htcl) / &
                        speccl(isp)%lain(spec(isp)%maxhtc)
            else
                speccl(isp)%seedprod(htcl) = 0
            end if
        end do byhtcls3
        
        !----- Here the  values for all  height classes higher than the max. one are set
        if (spec(isp)%maxhtc < maxhc) then
            resthtcls2: do htcl = spec(isp)%maxhtc + 1, maxhc
                speccl(isp)%htmin(htcl) = real(spec(isp)%khmax) + 1.0
                speccl(isp)%diain(htcl) = spec(isp)%maxdhm
                speccl(isp)%lain(htcl) = 0.0
                speccl(isp)%bioin(htcl) = 0.0
                speccl(isp)%seedprod(htcl) = 0
                
                latloop3: do latini = 1, maxlat
                    lonloop3: do lonini = 1, maxlon
                        stategrid(latini, lonini)%sp(isp)%numin(htcl) = 0.0
                    end do lonloop3
                end do latloop3
            end do resthtcls2
            
            speccl(isp)%htmin(spec(isp)%maxhtc + 1) = spec(isp)%khmax
            speccl(isp)%htmin(maxhc) = htmax
        end if
                
    end do byspcs1
    
    !----- Write out height and species dependent fixed transformation factors to files
    byhtcls4: do htcl = 0, maxhc
        write (biompertreehcloutfile%unit,          101) 0, tab, 0, tab, 0, tab,&
                (speccl(isp)%bioin(htcl), tab,  isp=1, nspc) , 0.
        write (laipertreehcloutfile%unit,           101) 0, tab, 0, tab, 0, tab,&
                (speccl(isp)%lain(htcl), tab,   isp=1, nspc) , 0.
        write (dbhpertreehcloutfile%unit,           101) 0, tab, 0, tab, 0, tab,&
                (speccl(isp)%diain(htcl), tab,  isp=1, nspc) , 0.
        write (seedprodfractpertreehcloutfile%unit, 101) 0, tab, 0, tab, 0, tab,&
                (speccl(isp)%seedprod(htcl), tab, isp=1, nspc) , 0.
        write (heightperhcloutfile%unit,            101) 0, tab, 0, tab, 0, tab,&
                (speccl(isp)%htmin(htcl), tab,  isp=1, nspc) , 0.
    end do byhtcls4  
    !----- Write out light class dependent relativ light intensities and LAI's leading to these
    write (lightandleafareaperlightclass%unit,102)  0,tab, 0,tab,0,tab,0, tab, ( ltval(ltcl),tab, ltcl=1,maxlc)
    write (lightandleafareaperlightclass%unit,102)  0,tab, 0,tab,0,tab,0, tab, ( la(ltcl),tab,    ltcl=1,maxlc)
    
    
    
    101 FORMAT(3(i1, a1), 60(e13.7, a1))
    102 FORMAT(4(i1, a1), 20(e11.5, a1))
    
    
     
    
end SUBROUTINE initandsetconstheightclprops
!
!done