!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*          Subroutine                                           *!
!* target thinning -                                             *!
!*          thinning routine with given values of biomass per    *!
!*          thinning year as target values                       *!
!*          rtargetm i given in kg DW/ha                         *!
!*                                                               *!
!*                  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 target_thinning_OC(i)

use data_stand
use data_manag
use data_simul
use data_species
use data_par

implicit none

real     :: rtargetm=0.             ! target value of stem biomass
integer     :: target_help  = 0.

real     ::  dbhmin=0,     &
             dbhmin_us = 0, &
             wpa=0,        &      ! Weibull parameter
             wpb=0,        &      !    -"-
             wpc=0,         &     !    -"-
             d63=0,         &
             help=0,        &
             pequal,      &
             tdbh=0,        &
             bas_area=0,    &
             bas_help=0.,   &
             target_help1=0,&
             dbh_h =0,      &
             db_l = 0.,     &
             db_u = 0.,     &
             d_est=0.,         &
             w_kb=0.,       &
             stembiom,      &
			 stembiom_us,   &
			 stembiom_re,    &
             stembiom_all,  &
             diff,          &
 			 mdiam,    &
			 mdiam_us


integer  ::  nrmin,      &
             lowtree,    &
			 undertree,   &
             flagth,     &
             taxnr,      &
             counth,     &
			 min_id,     &
			 max_id,     &
			 ih1,ih2,ncoh,   &
             coun1
! auxilarity for thinning routine 4: selective thinning
integer  :: count, i,  &
            idum,third, ipot, isel, ih
integer,dimension(0:anz_coh) :: cohl
integer, dimension(anz_coh) :: id_pot

real     :: h1, h2 , tar_h
real,external   :: gasdev
real:: ran0


! reacalculation of target to kg DW/patch
  h1 = 0.
  h2 = 0.
  count = 0
  cohl = -1
 flagth = 0
 coun1 = 0
 help=0.
 lowtree=0
 undertree = 0
 anz_tree_dbh = 0
 bas_area = 0.
! stem biomass of overstorey
 stembiom = 0.
! stem biomass of understorey
 stembiom_us = 0.
 stembiom_all = 0.
 tar_h = 300.

  if (time.eq.73.and. ip .eq.87) then

      stembiom = 0
  end if

taxnr = thin_spec(i)
mdiam = 0.
mdiam_us = 0.
! calculation of mean diameter (correspondung to med_diam) and basal area of stand
    zeig => pt%first
    DO
        IF (.NOT. ASSOCIATED(zeig)) EXIT

! Modification for V Kint: no test for diameter
             IF((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.0) THEN
! overstorey
				  stembiom = stembiom + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
                  help = help + zeig%coh%ntreea*(zeig%coh%diam**2)
                  bas_area = bas_area +  zeig%coh%ntreea*(zeig%coh%diam**2)*pi/4
                  if( zeig%coh%diam>0)  then
				        anz_tree_dbh = anz_tree_dbh +  zeig%coh%ntreea
                        mdiam = mdiam + zeig%coh%ntreea * (zeig%coh%diam**2)
				  end if

               ! Trees with DBH=0 for populations and per species; Baeume mit DBH =0 fuer Bestand und pro Spezie
             ELSE IF( (zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.1) THEN
! seedings/regeneration
			         stembiom_re = stembiom_re + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
                     lowtree = lowtree + zeig%coh%ntreea
             ELSE if((zeig%coh%ntreea>0).and.zeig%coh%species.eq.taxnr.and.zeig%coh%underst.eq.2) THEN
! understorey
			        stembiom_us = stembiom_us + (zeig%coh%x_sap + zeig%coh%x_hrt)*zeig%coh%ntreea
                     mdiam_us = mdiam_us + zeig%coh%ntreea * (zeig%coh%diam**2)
                     undertree = undertree + zeig%coh%ntreea

             ENDIF
        zeig => zeig%next
    ENDDO

! mean diameter for over and understorey
stembiom_all = stembiom + stembiom_us
if(anz_tree_dbh.ne.0) mdiam = sqrt(mdiam/real(anz_tree_dbh))
if(undertree.ne.0) mdiam_us = sqrt(mdiam_us/undertree)

third =  nint(anz_tree_dbh*0.333333)
anz_tree_ha =   nint(anz_tree_dbh*10000./kpatchsize)
anz_tree = anz_tree_dbh + undertree

 IF(anz_tree>0)THEN
    if(lowtree<anz_tree) help = sqrt(help/(anz_tree-lowtree))

 ENDIF

if(thin_stor(i).eq.0) then
    target_help = anz_tree_dbh
else
    target_help = undertree
end if

! tending
  if(thin_tysp(i).eq.4) target_help = stembiom_re

   if(target_mass(i).gt.1) then
      rtargetm = target_mass(i)*kpatchsize/10000.
   else
      rtargetm = target_mass(i)
   end if

! target value of biomass
   if(thin_tysp(i).eq.4) then
      rtargetm = stembiom_re - rtargetm*stembiom_re
    else
    end if

 ! cuttting
 if (rtargetm.eq.0.)then
  zeig => pt%first
      DO
        IF (.NOT. ASSOCIATED(zeig)) EXIT
        if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
               zeig%coh%ntreem =  zeig%coh%ntreea
               zeig%coh%ntreea = 0
               zeig%coh%nta = 0
        end if
        zeig=> zeig%next
      END DO

!tending of regeneration

else  if(thin_tysp(i).eq.4) then

  min_id = 1000
  max_id  = 0.
  zeig=>pt%first
   do
     if(.not.associated(zeig)) exit
     if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1) then
            ih1 = zeig%coh%ident
			if(ih1.lt.min_id) min_id = ih1
            ih2 = zeig%coh%ident
            if (ih2.gt.max_id) max_id = ih2
     end if
	 zeig=> zeig%next
   end do
   target_help1 = 0.
   do
       call random_number(pequal)
	   ncoh = min_id +(max_id-min_id)*pequal
	   zeig=>pt%first
       do
          if(.not.associated(zeig)) exit
          if(zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.1.and. zeig%coh%ident.eq.ncoh ) then
              zeig%coh%ntreea = zeig%coh%ntreea - 1
		      zeig%coh%nta = zeig%coh%nta-1
			  zeig%coh%ntreem = zeig%coh%ntreem  +1
              target_help = target_help - zeig%coh%ntreea
			  exit
          end   if
	      zeig=>zeig%next
        end do

	    diff = rtargetm -target_help
	    if(diff.lt.0.01) exit
	 end do

! different thinnings from below and above
else IF ( rtargetm .ne. 0.) then

    if(target_mass(i).lt.1.) then
     rtargetm = target_mass(i) * target_help
! No management if rtargetm=1
  else if (rtargetm.eq.1) then
     return
endif

     select case(thin_tysp(i))
          case(1)
! medium lower thinning
                d_est = 1.02
				w_kb = 2.5
          case(2)
! strong lower thinning
                d_est = 1.03
                w_kb = 1.5
          case(3)
! High thinning
                d_est = 1.04
                w_kb = 1.2
     end select



! calculation of Weibull-Parameter
    call min_dbh_overs(nrmin,dbhmin,taxnr)
    call min_dbh_unders(nrmin,dbhmin_us, taxnr)

    bas_help = bas_area
if (thin_stor(i).eq.0) then
     wpa = dbhmin
else
    wpa = dbhmin_us

end if 
if (thin_stor(i).eq.0) then
    d63 = mdiam*d_est
else
    d63 = mdiam_us * d_est
end if

    wpb = (d63 - wpa)/ w_kb
    wpc = 2
    wpc = 0.8
     
  if ((thin_tysp(i).ne.4) .and. rtargetm.lt.target_help) then
  
!selection of trees for thinning
  do
      call random_number(pequal)

        tdbh = wpa + wpb*(-log(1.-pequal))**(1./wpc)
        flagth = 0
! list of potential thinned tree chorts
   counth = 0
   id_pot = 0
   ipot = 1
    zeig => pt%first

    DO
        IF (.NOT. ASSOCIATED(zeig)) EXIT

          if(zeig%coh%notviable.eqv. .TRUE.) then
		   if(flag_mort.eq.0) then
		      id_pot(ipot)=zeig%coh%ident
		      ipot=ipot + 1
           endif
            else if(zeig%coh%ntreea.gt.0.and.zeig%coh%species.eq.taxnr.and. zeig%coh%underst.eq.thin_stor(i)) then
              dbh_h = zeig%coh%diam
              db_l =  dbh_h - 0.1*dbh_h
              db_u =  dbh_h + 0.1*dbh_h
			  counth = counth +1
              if (tdbh.ge.db_l.and.tdbh.le.db_u.and. zeig%coh%ntreea.ne. 0) then
                  id_pot(ipot) = zeig%coh%ident
                  ipot = ipot + 1
             end if           
         if(counth.gt. 100000) exit
         end if

         zeig=> zeig%next
    END DO  !  list of potential thinned tree cohorts

! selecting one equal distributed tree from the list of cohorts

    if ((ipot-1).ge.1) then
      if((ipot-1).eq.1) then

        isel = 1
      else
        call random_number(pequal)
		pequal = ran0(idum)
        isel = int(pequal*(ipot-1)) +1
      end if
      ih = id_pot(isel)
      zeig => pt%first
      DO
        IF (.NOT. ASSOCIATED(zeig)) EXIT

        if(zeig%coh%ident.eq. ih.and.zeig%coh%ntreea.ne. 0 ) then
                if(zeig%coh%notviable.eqv..TRUE.) then
				 if(flag_mort.eq.0) then
                   zeig%coh%ntreem = zeig%coh%ntreea
				   zeig%coh%ntreea=0
				   zeig%coh%nta=0
				   coun1=coun1+1
				   target_help = target_help - zeig%coh%ntreem
                  endif
		         else	  
			  zeig%coh%ntreea = zeig%coh%ntreea -1
              zeig%coh%nta =  zeig%coh%nta -1
              zeig%coh%ntreem = zeig%coh%ntreem +1
               coun1 = coun1 + 1
                target_help = target_help - 1
              exit
              end if
	     end if
         zeig =>zeig%next

      END DO ! thinning of one tree
    end if
    
    diff = target_help - rtargetm
    if(diff.le.0.1) exit
	if(coun1.gt.100000) exit
   end do  ! total thinning

 end if   ! thintype 1,2,3

 END IF ! all thinnings and tending

! adding biomasses to litter pools depending on stage of stand
 zeig=>pt%first

 do
   if(.not.associated(zeig)) exit
   if(zeig%coh%ntreem>0.and.zeig%coh%species.eq.taxnr) then
! all parts without stems of trees are input for litter
        h1 = zeig%coh%ntreem*zeig%coh%x_fol
        h2 = h2 + h1

         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
   endif
 zeig=>zeig%next

  enddo
call class_man
END SUBROUTINE target_thinning_OC
