!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*                  finishing simulation                         *!
!*                                                               *!
!*   contains                                                    *!
!*   FINISH_SIMUL: deallocation of variables,                    *!
!*                 closing files for each simulation             *!
!*   FINISH_ALL  : Finish all processes after all simulations    *!
!*   DEALLOC_SOIL: deallocation of soil variables                *!
!*                 (also used in other routines)                 *!
!*                                                               *!
!*                  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 finish_simul

use data_climate
use data_depo
!use data_effect
use data_evapo
use data_init
use data_manag
use data_out
use data_simul
use data_soil
use data_soil_cn
use data_species
use data_stand
use data_site
use data_tsort
use data_frost

! test lambda_ts
use data_par

implicit none

integer i ,unitout
character(150)   :: filename, infile
REAL       ::   rsap, cform
CHARACTER  ::   source

 rsap = 0.
 cform=0.
 source='U'
 infile='planting'

  if(time_out.gt.0) then

   if (flag_dis .eq. 2) then
    ! output of new tree.ini at the end of the simulation
    ! hier neue version rausschreiben zum einlesen in 4C mit NSC Speicher
    ! x_nsc_tb    x_nsc_crt   x_nsc_sap als kg C per tree
     unitout=getunit()
!     filename = trim(site_name(site_nr))//'_tree.ini'//trim(anh)
     filename = trim(site_name(ip))//'_tree.ini'//trim(anh)
     open(unitout,file=trim(dirout)//filename,status='replace')
     write(unitout,'(A1,2X,I1,1F12.0,A55)')'C', flag_volfunc, kpatchsize,' ! = ini version flag_dis, volume function, patch size'
     write(unitout,'(A)')'!     x_fol       x_frt       x_sap       x_hrt       x_Ahb     height   x_hbole  x_age      n     sp       DC          DBH &
          &      x_nsc_tb    x_nsc_crt   x_nsc_sap'

     zeig => pt%first
     do while (associated(zeig))
        write(unitout,'(5f12.5,2f10.0,i7,f7.0,i7, 5f12.5)') zeig%coh%x_fol, zeig%coh%x_frt, zeig%coh%x_sap, zeig%coh%x_hrt,  &
                                            zeig%coh%x_Ahb, zeig%coh%height, zeig%coh%x_hbole, zeig%coh%x_age, zeig%coh%ntreea,  &
                                            zeig%coh%species, zeig%coh%dcrb, zeig%coh%diam, zeig%coh%x_nsc_tb, zeig%coh%x_nsc_crt, zeig%coh%x_nsc_sap 

        zeig => zeig%next
     end do
    close(unitout)
  
    else

  ! output of new tree.ini at the end of the simulation
     unitout=getunit()
     filename = trim(site_name(ip))//'_tree.ini'//trim(anh)
     open(unitout,file=trim(dirout)//filename,status='replace')
     write(unitout,'(I1,1F12.0,A32)')flag_volfunc,kpatchsize,' ! = volume function, patch size'
     write(unitout,'(A)')'!     x_fol       x_frt       x_sap       x_hrt       x_Ahb     height   x_hbole  x_age      n     sp       DC          DBH'

     zeig => pt%first
     do while (associated(zeig))
        write(unitout,'(5f12.5,2f10.0,i7,f7.0,i7, 2f12.5)') zeig%coh%x_fol, zeig%coh%x_frt, zeig%coh%x_sap, zeig%coh%x_hrt,  &
                                            zeig%coh%x_Ahb, zeig%coh%height, zeig%coh%x_hbole, zeig%coh%x_age, zeig%coh%ntreea,  &
                                            zeig%coh%species, zeig%coh%dcrb, zeig%coh%diam
        zeig => zeig%next
     end do
    close(unitout)

    endif

  ! output of new .lit-file at the end of the simulation
     if (flag_end .eq. 0) then
         unitout=getunit()
         filename = trim(site_name(ip))//'.lit'//trim(anh)
         open(unitout,file=trim(dirout)//filename,status='replace')
         write(unitout,'(A,A)')'! litter initialisation ', site_name(ip)
         write(unitout,'(A)')'! fraction     Fagus sylvatica       Picea abies  Pinus sylvestris     Quercus robur    Betula pendula    Pinus contorta     Pinus ponderosa   Populus tremula   ground cover'
         write(unitout,'(A12, 9F18.1)') ' C_opm_fol  ', (slit(i)%C_opm_fol, i=1,nspecies)
         write(unitout,'(A12, 9F18.1)') ' C_opm_tb   ', (slit(i)%C_opm_tb, i=1,nspecies)
         write(unitout,'(A12, 9F18.1)') ' C_opm_frt  ', (slit(i)%C_opm_frt(1), i=1,nspecies)
         write(unitout,'(A12, 9F18.1)') ' C_opm_crt  ', (slit(i)%C_opm_crt(1), i=1,nspecies)
         write(unitout,'(A12, 9F18.1)') ' C_opm_stem ', (slit(i)%C_opm_stem,i=1,nspecies)
         close(unitout)
     endif

  end if   ! time_out

! deallocate cohorts
if(flag_end.ne.1 .and. associated(pt%first)) then
  zeig => pt%first
  do while (associated(zeig))
     pt%first => zeig%next
     deallocate (zeig%coh%frtrel)


     deallocate(zeig%coh%frtrelc)
     deallocate (zeig%coh%rooteff)

     if (flag_wred .eq. 9) deallocate (zeig%coh%rld)
     deallocate(zeig)
     zeig => pt%first
  end do
end if

if(associated(pt%first)) deallocate (pt%first)

if (flag_eva .gt.10) close (unit_eva)

if (allocated(dayfract))deallocate(dayfract)

!  fields for frost index
if(allocated(dnlf)) deallocate(dnlf)
if(allocated(tminmay_ann))deallocate(tminmay_ann)
if(allocated(date_lf)) deallocate(date_lf)
if(allocated(date_lftot)) deallocate(date_lftot)
if(allocated(dnlf_sp)) deallocate(dnlf_sp)
if(allocated(anzdlf)) deallocate(anzdlf)
if (allocated(sumtlf)) deallocate(sumtlf)


if (flag_clim==1) then
   if (allocated(recs))deallocate(recs)
   if (allocated(dd))deallocate(dd)
   if (allocated(mm))deallocate(mm);
   if (allocated(yy))deallocate(yy)
   if (allocated(tp))deallocate(tp);
   if (allocated(hm))deallocate(hm)
   if (allocated(prc))deallocate(prc);
   if (allocated(prs))deallocate(prs)
   if (allocated(rd))deallocate(rd)
   if (allocated(wd))deallocate(wd)
   if (allocated(tx))deallocate(tx)
   if (allocated(tn))deallocate(tn)
   if (allocated(vp))deallocate(vp)
   if (allocated(sdu))deallocate(sdu)
   if (allocated(sde))deallocate(sde)
   if (allocated(bw))deallocate(bw)

   if (allocated(tempfield))deallocate(tempfield)
   if (allocated(globfield))deallocate(globfield)
   if (allocated(dayfield))deallocate(dayfield)
endif

if (.not.flag_mult910) then
    if (allocated(NHd))deallocate(NHd)
    if (allocated(NOd))deallocate(NOd)
endif

if (allocated(diam_class))deallocate(diam_class)
if (allocated(diam_class_t))deallocate(diam_class_t)
if (allocated(diam_class_h))deallocate(diam_class_h)
if (allocated(diam_class_age))deallocate(diam_class_age)
if (allocated(diam_class_mvol))deallocate(diam_class_mvol)
if (allocated(diam_classm))deallocate(diam_classm)
if (allocated(diam_classm_h))deallocate(diam_classm_h)
if (allocated(height_class))deallocate(height_class)

if (allocated(ngroups))deallocate(ngroups)

if (allocated(dead_wood)) then
    do i = 1, nspec_tree
        deallocate(dead_wood(i)%C_tb)
        deallocate(dead_wood(i)%N_tb)
        deallocate(dead_wood(i)%C_stem)
        deallocate(dead_wood(i)%N_stem)
    enddo
    deallocate(dead_wood)
endif

svar%sumvsdead = 0.
svar%sumvsdead_m3 = 0.   
svar%daybb = 0.

if (flag_multi .eq. 1 .or. flag_multi .eq. 6 .or. flag_multi .eq. 0) then
   if(allocated(spar)) deallocate(spar)
   if(allocated(nrspec)) deallocate(nrspec)
   
   ! clear subfields for stress variables of svar
   if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
   do i=1,nspecies
   deallocate(svar(i)%tstress)
   deallocate(svar(i)%sstr)
   deallocate(svar(i)%BDstr)
   deallocate(svar(i)%BDmax)
   deallocate(svar(i)%porcrit)
   deallocate(svar(i)%airstr)
   deallocate(svar(i)%phstr)
   deallocate(svar(i)%Rstress)
   deallocate(svar(i)%Smean)
   enddo
   endif
   
   if(allocated(svar)) deallocate(svar)
endif

 if(flag_multi .eq. 4 .or. flag_mult8910) then
   do i=1,nspecies
      svar(i)%RedN = -99.9
   enddo
 end if

call dealloc_soil    ! soil-files immer deallok.

do i = 1,outy_n
   if (outy(i)%out_flag .ne. 0) then
      close (outy(i)%unit_nr)
   endif
enddo
do i = 1,outd_n
   if (outd(i)%out_flag .ne. 0) then
      close (outd(i)%unit_nr)
   endif
enddo

C_bc_tot = 0.
N_bc_tot = 0.
if (flag_bc .gt. 0) then
    deallocate(C_bc)
    deallocate(N_bc)
    deallocate (C_bc_appl)
    deallocate (N_bc_appl)
    deallocate (bc_appl_lay)
    deallocate (cnv_bc)
    deallocate (dens_bc)
    deallocate (cpart_bc)
    deallocate (y_bc)
    flag_decomp = flag_decomp + 100   ! flag_decomp zurücksetzen
endif

if (flag_cohout .ge. 1) then
   do i = 1,outcy_n
      if (outcy(i)%out_flag .ne. 0) then
         close (outcy(i)%unit_nr)
      endif
   enddo
endif

if (flag_dayout .ge. 1) then
   do i = 1,outcd_n
      if (outcd(i)%out_flag .ne. 0) then
         close (outcd(i)%unit_nr)
      endif
   enddo

endif

if(time_out .gt. 0) then
   if (out_flag_light .ne. 0) close(unit_light)
   if (flag_cohout .eq. 2) then
      close(unit_prod)
      close(unit_allo)
   endif
 end if

if (flag_dayout .gt. 1) then
   close(unit_wat)
   close(unit_soicnd);close(unit_soicna)
endif

if (.not.flag_mult910) close (unit_soil)
if (flag_sum > 0) close(unit_sum)
if (flag_mg==1) then
 deallocate(thin_year);deallocate(thin_tree)
endif
if (flag_mg==3.or. flag_mg==33) then
    if (allocated(thin_year)) deallocate(thin_year)
    if( allocated(target_mass)) deallocate(target_mass)
    if (allocated(thin_tysp))deallocate(thin_tysp)
    if (allocated(thin_spec))deallocate(thin_spec)
    if (allocated(rot))deallocate(rot)
    if (allocated(thin_flag1))deallocate(thin_flag1)
    if (allocated(thinyear))deallocate(thinyear)
    if (allocated(thin_stor))deallocate(thin_stor)
endif
if (flag_mg==2.and. flag_end==0) then
 if (allocated(zbnr))deallocate(zbnr)
 if (allocated(tend))deallocate(tend)
 if (allocated(rot))deallocate(rot)
 if (allocated(regage))deallocate(regage)
 if (allocated(thin_flag1))deallocate(thin_flag1)
 if (allocated(thin_flag2))deallocate(thin_flag2)
 if (allocated(thin_flag3))deallocate(thin_flag3)
 if (allocated(thin_flag4))deallocate(thin_flag4)
 if (allocated(np_mod))deallocate(np_mod)
 if (allocated(specnr))deallocate(specnr)
 if (allocated(age_spec))deallocate(age_spec)
 if (allocated(anz_tree_spec))deallocate (anz_tree_spec)
 if (allocated(thinyear))deallocate(thinyear)
end if
if (flag_mg==4. .or. flag_mg == 5) then
    if (allocated(thin_flag1))  deallocate(thin_flag1)
end if
if(flag_mg == 10) then
  if (allocated(thin_flag1))deallocate(thin_flag1)
  if (allocated(dis_id))deallocate(dis_id)
  if (allocated(dis_type))deallocate(dis_type)
  if (allocated(fortype))deallocate(fortype)
  if (allocated(dis_year))deallocate(dis_year)
  if (allocated(dis_rel))deallocate(dis_rel)
  if (allocated(sum_dis))deallocate(sum_dis)
 end if
if(flag_dis == 1 .or. flag_dis == 2) then
  if (allocated(dis_year))deallocate(dis_year)
  if (allocated(dis_spec))deallocate(dis_spec)
  if (allocated(dis_start))deallocate(dis_start)
  if (allocated(dis_rel))deallocate(dis_rel)
  if (allocated(dis_type))deallocate(dis_type)
 end if

if(flag_mg == 9) then
  if (allocated(thin_flag1))deallocate(thin_flag1)
  if (allocated(yman))deallocate(yman)
  if (allocated(dbh_clm))deallocate(dbh_clm)
  if (allocated(rem_clm))deallocate(rem_clm)
  if (allocated(spec_man))deallocate(spec_man)
  if (allocated(act))deallocate(act)
   if (allocated(rel_part))deallocate(rel_part)
end if
if(flag_mg == 8) then
  if (allocated(thin_flag1))deallocate(thin_flag1)
  if (allocated(yman))deallocate(yman)
   if (allocated(rel_part))deallocate(rel_part)
end if

if(flag_wpm.ne.0) then
	! free the resources
	call deallocate_wpm

IF ( associated(st%first)) then
   ztim => st%first
   do while (associated(ztim))
      st%first => ztim%next
      deallocate(ztim)
      ztim => st%first
   end do
endif

 IF ( associated(st%first)) deallocate(st%first)
 if ( associated(ztim)) deallocate(ztim)
end if

! test lambda_ts
if (flag_lambda.eq.1) then
deallocate (lambda_ts)
end if
    
! compressed output for each simulation run
lcomp1 = .TRUE.
end subroutine finish_simul

!-----------------------------------------

SUBROUTINE finish_all

use data_simul
use data_climate
use data_depo
use data_mess
use data_out
use data_site
use data_soil
use data_soil_cn
use data_species
use data_stand

if (allocated(site_name))deallocate(site_name)
if (allocated(climfile))deallocate(climfile);
if (allocated(sitefile))deallocate(sitefile)
if (allocated(valfile))deallocate(valfile)
if (allocated(treefile))deallocate(treefile)
if (allocated(wpmfile))deallocate(wpmfile)
if (allocated(depofile))deallocate(depofile)
if (allocated(redfile))deallocate(redfile)
if (allocated(litfile))deallocate(litfile)
if (allocated(standid))deallocate(standid)

IF(ALLOCATED(thick)) CALL dealloc_soil

if(flag_multi .eq. 1 .or. flag_multi .ge. 3) then
    if ( allocated(sitenum))deallocate(sitenum)
    if ( allocated(clim_id))deallocate(clim_id)
    if ( allocated(soilid))deallocate(soilid)
    if ( allocated(gwtable))deallocate(gwtable)
    if ( allocated(NOdep))deallocate(NOdep)
    if ( allocated(NHdep))deallocate(NHdep)
endif

if(allocated(diam_class)) deallocate(diam_class)
if(allocated(diam_class_t)) deallocate(diam_class_t)
if(allocated(diam_class_h)) deallocate(diam_class_h)
if(allocated(diam_classm)) deallocate(diam_classm)
if(allocated(diam_classm_h)) deallocate(diam_classm_h)
if(allocated(height_class)) deallocate(height_class)

if (allocated(NHd))deallocate(NHd)
if (allocated(NOd))deallocate(NOd)

if(allocated(recs))then
 deallocate(recs)
 deallocate(dd);deallocate(mm);deallocate(yy)
 deallocate(tp);deallocate(hm);deallocate(prc);deallocate(prs)
 deallocate(rd)
 if (allocated(tempfield))deallocate(tempfield)
 if (allocated(globfield))deallocate(globfield)
 if (allocated(dayfield))deallocate(dayfield)
endif

if(time_out .ne. -2) then
   close(unit_comp1)
   close(unit_comp2)
endif

if (flag_stat .gt. 0) then
    close(unit_cons)
    close(unit_mess)
    close(unit_stat)
endif

if (flag_multi .gt.8)  close (output_unit_all)

if (flag_multi .eq. 2) close(unit_ctr)
if(flag_multi.eq.7) deallocate(fl_co2)

if(flag_multi .eq. 4 .or. flag_mult8910) then
    if (allocated(output_var))deallocate(output_var)
    if (allocated(output_varm))deallocate(output_varm)
    if (allocated(output_varw))deallocate(output_varw)
    if (allocated(climszenres))deallocate(climszenres)
    if (allocated(climszenyear))deallocate(climszenyear)
    if (allocated(climszenmon))deallocate(climszenmon)
    if (allocated(climszenweek))deallocate(climszenweek)
endif

if ((ip .eq. 1 .or. flag_multi .eq. 1 .or. flag_multi .eq. 6) .and. (time_out .ne. -2) ) close(unit_err)

end subroutine finish_all

!-----------------------------------------

SUBROUTINE  dealloc_soil

use data_soil
use data_soil_cn
use data_soil_t
use data_simul

implicit none

if (allocated(thick)) deallocate(thick)
if (allocated(mid)) deallocate(mid)
if (allocated(depth)) deallocate(depth)
if (allocated(pv)) deallocate(pv)
if (allocated(pv_v)) deallocate(pv_v)
if (allocated(dens)) deallocate(dens)
if (allocated(f_cap_v)) deallocate(f_cap_v)
if (allocated(wilt_p_v)) deallocate(wilt_p_v)
if (allocated(field_cap)) deallocate(field_cap)
if (allocated(wilt_p)) deallocate(wilt_p)
if (allocated(vol)) deallocate(vol)
if (allocated(quarzv)) deallocate(quarzv)
if (allocated(sandv)) deallocate(sandv)
if (allocated(clayv)) deallocate(clayv)
if (allocated(siltv)) deallocate(siltv)
if (allocated(humusv)) deallocate(humusv)
if (allocated(dmass)) deallocate(dmass)
if (allocated(fcaph)) deallocate(fcaph)
if (allocated(wiltph)) deallocate(wiltph)
if (allocated(pvh)) deallocate(pvh)
if (allocated(skelv)) deallocate(skelv)
if (allocated(skelfact)) deallocate(skelfact)
if (allocated(spheat)) deallocate(spheat)
if (allocated(phv)) deallocate(phv)
if (allocated(wlam)) deallocate(wlam)
if (allocated(wats)) deallocate(wats)
if (allocated(watvol)) deallocate(watvol)
if (allocated(wat_res)) deallocate(wat_res)
if (allocated(perc)) deallocate(perc)
if (allocated(wupt_r)) deallocate(wupt_r)
if (allocated(wupt_ev)) deallocate(wupt_ev)
if (allocated(s_drought)) deallocate(s_drought)
if (allocated(root_fr)) deallocate(root_fr)
if (allocated(temps)) deallocate(temps)
if (allocated(BDopt)) deallocate(BDopt)
if (allocated(fr_loss)) deallocate(fr_loss)
if (allocated(redis)) deallocate(redis)
if (allocated(C_opm)) deallocate(C_opm)
if (allocated(C_hum)) deallocate(C_hum)
if (allocated(C_opmfrt)) deallocate(C_opmfrt)
if (allocated(C_opmcrt)) deallocate(C_opmcrt)
if (allocated(N_opm)) deallocate(N_opm)
if (allocated(N_hum)) deallocate(N_hum)
if (allocated(N_opmfrt)) deallocate(N_opmfrt)
if (allocated(N_opmcrt)) deallocate(N_opmcrt)
if (allocated(NH4)) deallocate(NH4)
if (allocated(NO3)) deallocate(NO3)
if (allocated(Nupt)) deallocate(Nupt)
if (allocated(Nmin)) deallocate(Nmin)
if (allocated(rmin_phv)) deallocate(rmin_phv)
if (allocated(rnit_phv)) deallocate(rnit_phv)
if (allocated(cnv_opm)) deallocate(cnv_opm)
if (allocated(cnv_hum)) deallocate(cnv_hum)
if (allocated(slit)) deallocate(slit)
if (allocated(slit_1)) deallocate(slit_1)
if (allocated(sh)) deallocate(sh)
if (allocated(sv)) deallocate(sv)
if (allocated(sb)) deallocate(sb)
if (allocated(sbt)) deallocate(sbt)
if (allocated(t_cond)) deallocate(t_cond)
if (allocated(t_cb)) deallocate(t_cb)
if (allocated(h_cap)) deallocate(h_cap)
if (allocated(sxx)) deallocate(sxx)
if (allocated(svv)) deallocate(svv)
if (allocated(svva)) deallocate(svva)
if (allocated(soh)) deallocate(soh)
if (allocated(son)) deallocate(son)
if (allocated(wat_root)) deallocate(wat_root)
if (allocated(root_lay)) deallocate(root_lay)
if (allocated(gr_depth)) deallocate(gr_depth)

if (allocated(xwatupt)) deallocate (xwatupt)
if (allocated(xNupt)) deallocate (xNupt)
if (allocated(wat_left)) deallocate (wat_left)

end subroutine dealloc_soil
!-----------------------------------------------------------------


