!*****************************************************************!
!*                                                               *!
!*                     4C (FORESEE)                              *!
!*                                                               *!
!*                                                               *!
!*                Subroutines for:                               *!
!*                Austrian  management                           *!
!*                 contains:                                     *!
!*                SR aust_ini                                    *!
!*                SR aust_manag                                  *!
!*                SR plant_aust                                  *!
!*                SR calc_rel_class                              *!
!*                                                               *!
!*                  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 aust_ini
 use data_manag
 use data_species
 use data_simul
 use data_stand
 implicit none

 integer :: manag_unit,i, ih1,ih2,ios,ih4, flp , flag_help
 character(len=150) :: filename
 logical :: ex
 character ::text
 real    :: hp, ih3

 manag_unit=getunit()
 filename = manfile(ip)

 allocate(thin_flag1(nspec_tree))
 flag_help = 0
 thin_flag1=-1
 thin_dead = 1
 allocate(yman(1000))
 allocate(dbh_clm(1000))
 allocate(rem_clm(1000))
 allocate(spec_man(1000))
 allocate(act(1000))
 allocate(rel_part(1000))
yman = 0
 dbh_clm = 0
 rem_clm = 0.
 spec_man = 0
 act = 0
 rel_part = 0
 flp = 0
  call testfile(filename,ex)
  open(manag_unit,file=trim(filename))
! read head of data-file
 do
    read(manag_unit,*) text
    if(text .ne. 's')then

       backspace(manag_unit);exit
    endif
 enddo

 i=1
 do
    read(manag_unit,*,iostat=ios) ih1,ih2, ih3, hp,ih4
    if(ios<0) exit
      yman(i) = ih1            ! year of treatment
      if(ih2.eq.1) then
! Fichte/ Spruce
          spec_man(i) = 2
      else if(ih2.eq.2) then
! Kiefer/ Pine
          spec_man(i) = 3
      else if(ih2.eq.3) then
! Eiche/ oak
          spec_man(i) = 4
      else if(ih2.eq.4) then
          spec_man(i) = 1
      end if

        ! species  number
      act(i) = ih4
      if(ih1.ne.-999 ) then
         if(flp.eq.0) then
            dbh_clm(i) = int(ih3)    ! dbh-cluss number for treatment
            rem_clm(i) = hp          ! removal of biomass
            i = i+1
         else
            act(i) = ih4
            rel_part(i) = ih3
            rem_clm(i) = 0
            i = i+1
         end if
      else
       if(i.eq.1) thin_dead = 0
       flp = 1
       backspace(manag_unit)
      end if
 end do
 num_man = i-1

 close(manag_unit)
END SUBROUTINE aust_ini
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE aust_manag

 use data_manag
 use data_stand
 use data_simul
 use data_species
 use data_par

 implicit none
 integer    :: i,j,hcl, ha, taxnr,k,l, help_fl,helpz, helps
 real, dimension(5)   :: rel_biom, harv_biom
 real, dimension(5)   :: num_ccl, contr
 real,dimension(5)    :: help_rem_clm
 integer,dimension(5) :: help_rel_dbh, hrd
 real                 :: stump_dw, stump_v, hrb
 ha =0
 rel_biom = 0.
 num_ccl =0
 harv_biom = 0.
 help_rel_dbh = 0
 help_rem_clm = 0.
 hrd = 0
 helpz = 0
 helps = 0
  call calc_rel_class

! calculation of stem biomass of relative dbh-class
  zeig=>pt%first
  do
    if(.not.associated(zeig)) exit
    if(zeig%coh%diam.ne.0) then
      hcl = zeig%coh%rel_dbh_cl
      if(hcl.ne.0) then
          num_ccl(hcl)= num_ccl(hcl) +1
          rel_biom(hcl)= rel_biom(hcl) + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
      end if
     end if
       zeig=>zeig%next
  end do

 do l=1,nspecies
  help_rel_dbh = 0
  help_rem_clm = 0.
  helpz = 0
  helps = 0
! calculation of stem biomass of relative dbh-class
  zeig=>pt%first
  do
    if(.not.associated(zeig)) exit
    if(zeig%coh%diam.ne.0.and.zeig%coh%species.eq.l) then
      hcl = zeig%coh%rel_dbh_cl
      if(hcl.ne.0) then
          num_ccl(hcl)= num_ccl(hcl) +1
          rel_biom(hcl)= rel_biom(hcl) + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
      end if
     end if
       zeig=>zeig%next
  end do

  hrd=0
   do i=1,num_man

    if(yman(i).eq.time) then

      if(act(i) .eq.1.and.spec_man(i).eq.l) then
         zeig=>pt%first
         do
           if(.not.associated(zeig)) exit
           if(zeig%coh%diam.ne.0) then
            if(zeig%coh%species.eq.l) then
               hrd(zeig%coh%rel_dbh_cl)= 1
            end if
           end if
           zeig=>zeig%next
         end do
         help_rel_dbh(dbh_clm(i)) = 1
         help_rem_clm(dbh_clm(i)) = rem_clm(i)
       end if ! act(i)
    end if !yman(i)
   end do ! num_man
      do j=1,5
         if(help_rel_dbh(j).eq.1.and.hrd(j).eq.0) then
                if(j.eq.1.) then
                     do k=2,5
                       if(hrd(k).ne.0) then
                         help_rem_clm(k) = help_rem_clm(k) + help_rem_clm(j)
                         help_rel_dbh(k)=1
                         exit
                       end if
                     end do
                else if (j.eq.5.) then
                       do k= 4,1,-1
                           if(hrd(k).eq.1) then
                             help_rem_clm(k) = help_rem_clm(k) + help_rem_clm(j)
                             help_rel_dbh(k) = 1
                             exit
                           endif
                        end do
                 else
                    do k=j,5
                        if(hrd(k).eq.1) then
                             help_rem_clm(k) = help_rem_clm(k) + help_rem_clm(j)*0.5
                              help_rel_dbh(k) = 1
                              exit
                        end if
                    end do
                    do k=j,1,-1
                        if(hrd(k).eq.1) then
                             help_rem_clm(k) = help_rem_clm(k) + help_rem_clm(j)*0.5
                              help_rel_dbh(k) = 1
                              exit
                        end if
                    end do
                 end if
          help_rel_dbh(j) = 0
          help_rem_clm(j) = 0.
          end if

      end do
! thinning
  help_fl = 0
   do i=1,num_man

    if(yman(i).eq.time.and.help_fl.eq.0) then
                do k=1,5
                 helps = helps +  help_rel_dbh(k)
               end do

       help_fl=1
       zeig=>pt%first
       do
           if(.not.associated(zeig)) exit
           if(zeig%coh%diam.ne.0.and.zeig%coh%species.eq.l) then


               do k=1,5
                  if(zeig%coh%rel_dbh_cl.eq.k.and.help_rel_dbh(k).eq.1) then
                     if(help_rem_clm(k).gt.1.) help_rem_clm(k) = 1.
                     if( help_rem_clm(k) .eq. 1.)then

                         if(zeig%coh%underst.eq.0.and.zeig%coh%x_age.gt. 20) ha=int(help_rem_clm(k)* zeig%coh%ntreea+0.5)
                         helpz = helpz +1
                     else
                         ha=int(help_rem_clm(k)* zeig%coh%ntreea+0.5)

                     end if
                     if(ha.lt.1) ha = 1
                     if(help_rem_clm(k) .ne.1) then
                          harv_biom(k) = harv_biom(k) + ha* (zeig%coh%x_sap + zeig%coh%x_hrt)
                          hrb =  help_rem_clm(k)* rel_biom(k)
                          if(harv_biom(k).eq.rel_biom(k)) then
                               ha = ha -1
                          end if
                      end if

                     zeig%coh%ntreea = zeig%coh%ntreea - ha
                     zeig%coh%nta = zeig%coh%ntreea
                     zeig%coh%ntreem = zeig%coh%ntreem + ha

                  end if
               end do  ! k loop
            end if
            zeig=>zeig%next
      end do   ! zeig loop

    end if
  end do ! num_man

 if(helps.gt.0.and.helpz.ge.helps) then
    zeig=>pt%first
    do
           if(.not.associated(zeig)) exit
           if(zeig%coh%species.eq.l.and.zeig%coh%underst.eq.1) then
                   zeig%coh%underst = 0
           end if
           zeig => zeig%next
    end do
 end if

 write(9898,*) time, 'totbio', rel_biom
 write(9898,*) time, 'harvbio', harv_biom
  do i=1,5
  if(rel_biom(i).ne.0.) then
      contr(i) = harv_biom(i)/rel_biom(i)
  else
      contr(i) = 0.
  end if
  end do
 write(9898,*) time,l, contr

  rel_biom = 0.
  harv_biom = 0.
 end do ! nspecies

 ! planting

 do i=1,num_man
   if(yman(i).eq.time.and.act(i).ne.1) then
      call plant_aust(i)

    end if  ! act
 end do

stump_sum = 0
 zeig=>pt%first

 do
   if(.not.associated(zeig)) exit
   taxnr=zeig%coh%species

   if(zeig%coh%ntreem>0)then
! all parts without stems of trees are input for litter
         zeig%coh%litC_fol = zeig%coh%litC_fol + zeig%coh%ntreem*(1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart
         zeig%coh%litN_fol = zeig%coh%litN_fol + zeig%coh%ntreem*((1.-spar(taxnr)%psf)*zeig%coh%x_fol*cpart)/spar(taxnr)%cnr_fol
         zeig%coh%litC_frt = zeig%coh%litC_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart
         zeig%coh%litN_frt = zeig%coh%litN_frt + zeig%coh%ntreem*zeig%coh%x_frt*cpart/spar(taxnr)%cnr_frt
         zeig%coh%litC_tb = zeig%coh%litC_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart
         zeig%coh%litN_tb = zeig%coh%litN_tb + zeig%coh%ntreem*zeig%coh%x_tb*cpart/spar(taxnr)%cnr_tbc
         zeig%coh%litC_crt = zeig%coh%litC_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart
         zeig%coh%litN_crt = zeig%coh%litN_crt + zeig%coh%ntreem*zeig%coh%x_crt*cpart/spar(taxnr)%cnr_crt

! stumps into stem litter
          call stump( zeig%coh%x_ahb, zeig%coh%asapw,zeig%coh%dcrb,zeig%coh%x_hbole,  &
                      zeig%coh%height, taxnr,stump_v, stump_dw)
          zeig%coh%litC_stem = zeig%coh%litC_stem +  zeig%coh%ntreem*stump_dw*cpart
          zeig%coh%litN_stem = zeig%coh%litC_stem/spar(taxnr)%cnr_stem
          stump_sum = stump_sum + zeig%coh%ntreem*stump_dw
   endif
 zeig=>zeig%next
 enddo

 sumvsab = 0.
 sumvsab_m3 = 0.
 svar%sumvsab = 0.

 zeig=>pt%first
   do while (associated(zeig)) 

     ns = zeig%coh%species
     sumvsab          = sumvsab + zeig%coh%ntreem*(zeig%coh%x_sap + zeig%coh%x_hrt)
     sumvsab_m3       = sumvsab_m3 +  zeig%coh%ntreem*(zeig%coh%x_sap + zeig%coh%x_hrt)/(spar(ns)%prhos*1000000)
     svar(ns)%sumvsab = svar(ns)%sumvsab +  zeig%coh%ntreem*(zeig%coh%x_sap + zeig%coh%x_hrt)

     zeig=>zeig%next

   end do
  sumvsab = sumvsab *  10000./kpatchsize                 ! kg/ha
  sumvsab_m3 = sumvsab_m3 *  10000./kpatchsize           ! kg/ha


  do k = 1, nspec_tree
    svar(k)%sumvsab = svar(k)%sumvsab * 10000./kpatchsize           ! kg/ha
  end do
! cumulative harvested stem mass
  cumsumvsab = cumsumvsab + sumvsab

 if(thin_dead.ne.0) then
     call class_man
 end if

END SUBROUTINE aust_manag

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
SUBROUTINE plant_aust(mp)
use data_manag
use data_plant
use data_species
use data_stand

implicit none
integer     :: fl_plant, i, nplant,taxid,mp
 real       :: age,          &
               pl_height,    &
               sdev,         &
               plhmin

  infspec = 0
  npl_mix = 0
 fl_plant = act(mp)
 select case(fl_plant)
    case(2)
      infspec(2) = 1
      npl_mix(2) = 2500
      
    case(3)
      infspec(2) = 1
      npl_mix(2) = 10000
      
    case(4)
      infspec(3) = 1
      npl_mix(3) = 5000

    case(5)
      infspec(3) = 1
      npl_mix(3) = 2000

    case(6)
      infspec(1) = 1
      npl_mix(1) = 500

    case(7)
      infspec(1) = 1
      npl_mix(1) = 5000

    case(8)
      infspec(4) = 1
      npl_mix(4) = 5000

    case(9)
      infspec(1) = 1
      npl_mix(1) = 1000

      infspec(4) = 1
      npl_mix(4) = 3500

    case(10)
      infspec(3) = 1
      npl_mix(3) = 2500

      infspec(4) = 1
      npl_mix(4) = 2500

    case(11)
      infspec(3) = 1
      npl_mix(3) = 2500

      infspec(1) = 1
      npl_mix(1) = 2500

    case(12)
      infspec(3) = 1
      npl_mix(3) = 7000

    case(13)
      infspec(4) = 1
      npl_mix(4) = 2500

 end select
       do i = 1,nspec_tree
         if (infspec(i).eq.1) then

               taxid = i
! data for Austria
               age = pl_age(taxid)
               pl_height = plant_height(taxid)
               plhmin = plant_hmin(taxid)
               nplant = rel_part(mp)*nint(npl_mix(taxid)*kpatchsize/10000)
               sdev = hsdev(taxid)
               call gener_coh(taxid, age, pl_height, plhmin, nplant,sdev)

          end if
      end do  ! i

END SUBROUTINE plant_aust

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
    
SUBROUTINE calc_rel_class
 use data_manag
 use data_stand
 use data_species

 implicit none
 integer    :: nrmax, i, j, k, adm
 real,dimension(10)       :: maxdbh, mindbh,class_wd

 integer  :: nrmin
 real     :: help, help_h1
 class_wd =0.
 maxdbh = 0.
 mindbh = 0.

do  j= 1,nspecies
  call max_dbh(nrmax,help,adm, j)
  call min_dbh(nrmin,help_h1,adm, j)
  zeig=>pt%first

  do
    if(.not.associated(zeig)) exit

       if(zeig%coh%ident.eq.nrmax.and.zeig%coh%species.eq.j) then
          maxdbh(j) = help

       else if(zeig%coh%ident.eq.nrmin.and.zeig%coh%species.eq.j) then
          mindbh(j) = help_h1

       end if

     zeig=>zeig%next
  end do
end do

 do j=1,nspecies

     class_wd(j) =  (maxdbh(j)-mindbh(j))/5
     k = 5
    zeig=>pt%first

    do
       if(.not.associated(zeig)) exit
          if(zeig%coh%species.eq.j.and. zeig%coh%diam.gt.0) then

               do i=1,k
                 if(zeig%coh%diam.ge.(mindbh(j)+class_wd(j)*(i-1)).and.zeig%coh%diam.lt.(mindbh(j)+class_wd(j)*i)) then
                      zeig%coh%rel_dbh_cl = i
                      exit
                 else if (zeig%coh%diam.eq.maxdbh(j)) then
                       zeig%coh%rel_dbh_cl = 5
                 end if
               end do
           end if

         zeig=>zeig%next
    end do
 end do
END SUBROUTINE calc_rel_class
