!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*    - Simulation control: SIM_CONTROL                          *!
!*                          SIMULATION_4C                        *!
!*                                                               *!
!*                  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 sim_control

use data_climate
use data_simul
use data_site
use data_out

implicit none

integer run_nr, ipp, irl, icl, i
character a
character(8) actdate
character(10) acttime, helpsim, text1, text2
real time1, time2, time3
logical lhelp

unit_err=getunit()
if(flag_multi.eq.5) dirout = './'
open(unit_err,file=trim(dirout)//trim(site_name(1))//'_error.log',status='replace', position='append')
unit_trace=getunit()
open(unit_trace,file=trim(dirout)//trim(site_name(1))//'_trace.log',status='replace', position='append')
write (unit_trace, *) ' Trace of calls  -  subroutines of 4C '
write (unit_trace, *) 
write (unit_trace, *) 'iday  time_cur    subroutine '
write (unit_trace, '(I4,I10,A)') iday, time_cur, ' sim_control'

! check daily output
if (year > 5 .and. flag_dayout .ge. 1) then

    lhelp = .true.
    do i = 1,outd_n
       if (outd(i)%out_flag .eq. flag_dayout) then

          select CASE (outd(i)%kind_name)

          CASE ('day_short')
            lhelp = .false.
          end select
       endif
    enddo

   if (lhelp) then
       write(*,*) '  Warning: Your choice of daily output is ON with a simulation time of'
       write(*,'(I6,A,I8,A)') year,'  years. This option will create ',365*year,' data records per file!'
       write(*,'(A)',advance='no') '  Do you really want do use daily output (y/n)? '
       read *,a
       IF (a .eq. 'n' .or. a .eq. 'N') then
          flag_dayout = 0
       ENDIF
   endif  ! lhelp
ENDIF

! open file ycomp (yearly compressed output (multi run))
IF (time_out .ne. -2) call prep_out_comp

!call epsilon
	IF(flag_multi.eq.1) THEN
	   run_nr = site_nr
	ELSE IF (flag_multi.eq.5) THEN
	   run_nr = 1
	ELSE
	   run_nr = repeat_number
	ENDIF

call date_and_time(actdate, acttime)
write (unit_err, *)
time3 = 0.

  if (.not.flag_mult910) then
    nrreal = 1
    nrclim = 1
  endif
   
  do icl = 1, nrclim   ! climate scenarios
    iclim = icl
	DO ipp = 1, run_nr ! sites
     ip = ipp
     if (flag_trace) write (unit_trace, '(I4,I10,A,I3)') iday, time_cur, ' sim_control ip=',ip  
      do irl = 1, nrreal  ! realization f climate scenarios
           if (flag_mult910) then
               climfile(ip) = climszenfile(ip, icl, irl)
               site_name (ip) = trim(site_name1)//'_'//trim(sitenum(ip))
               write (helpsim,'(I10)') icl
               read (helpsim,*) text1
               write (helpsim,'(I10)') irl
               read (helpsim,*) text2
               site_name (ip) = trim(site_name (ip))//'_'//trim(text1)//'_'//trim(text2)
               write (unit_err, *)
               write (unit_err, '(A,3I5,2X, A50)')'*  ip, cli-scenario, real., site:  ', ip, icl, irl, site_name(ip)
               write (unit_err, '(A,A)') 'Climate file:  ', trim(climfile(ip))
  	       else
               write (unit_err, *)
               write (unit_err, '(A10,I5,2X, A50)') ' ip/site ', ip, trim(site_name(ip))
               site_name1 = trim(site_name(ip))
           endif
           call CPU_time (time1)
           if(ip.ne.0) then
	           CALL sim_ini
	           CALL prepare_site
	           if (flag_multi.eq.5) then
!                  call m4c_simenv_in
                  unit_comp2 = 6    ! standard output 
               end if

	           if(flag_end.gt.0) then
                  select case (flag_end)
                  case (1)
     		        print*,ip, '  stop in prepare_stand (see error.log)'
                  case (2)
     		        print*, ip, 'stop in prepare_stand, stand ',   &
     		        'identificator not found in prepare_stand'
                  case (3)
                    print*,ip, 'stop in canopy'
                  case (4)
                    print*,ip, 'stop in readsim, climate ID not found'
                  case (5)
                    print*,ip, ' stop in readsoil, soil ID not found ', adjustl(soilid(ip))
                  case (6)
                    write(*,'(A,I5)') '  >>>foresee message: stop in read_cli - no climate data for year ',time_b
                    call finish_simul
                    stop
                  case default
                    print*,ip, 'flag_end = ', flag_end
                  end select

	              call finish_simul
	              flag_end = 0
	           else
	              IF(flag_multi==2) CALL fixclimscen
                    if (.not.flag_mult910) then
                        write (*,*)
                        write (*,*)  '>>> Start FORESEE-Simulation site ', ipp
                        write (*,*)
                    endif

	              CALL simulation_4c
	              CALL finish_simul
	           endif
               if (flag_mult910)  then
                    call out_var_stat(1, irl)
               else
                    if ((flag_multi .ne. 8) .and. (nvar .gt. 1)) call out_var_stat(3, 1)
               endif 
           endif  ! flag_end
 	       call CPU_time (time2)
           if (.not.flag_mult910)  then
               print *, ' run time for simulation ',ip, time2-time1, ' sec'
           endif
           write (unit_err, *) ' run time for simulation ',ip, time2-time1, ' sec'
           time3 = time3 + (time2-time1)
	  enddo  ! irl 
      if (flag_mult910)  call out_var_stat(2, -99) 
      write (unit_err, *)
      write (unit_err, *)
      write (unit_err, *) '* * * * *  New ip/site  * * * * *'
    ENDDO   ! ip until site_nr (page number)
    write (unit_err, *)
    write (unit_err, *) '************  New climate scenario  **********'
  enddo   ! icl

	if (nvar .gt. 1) then 
        select case (flag_multi)
        case (5, 9, 10)
            continue
        case (1)
            continue
        case default    
            call out_var_file
        end select
    endif

!  comparison with measurements
    if (flag_stat .gt. 0) CALL mess	
    
 	call CPU_time (time1)
    time3 = time3 + (time1-time2)
    write (unit_err, *)
    write (unit_err, *) ' total run time ', time3, ' sec'
    CALL finish_all
	PRINT *,'All simulations finished!'

END SUBROUTINE sim_control

!**************************************************************

SUBROUTINE simulation_4C

  !*** Declaration part ***!

  USE data_simul
  USE data_species    ! species specific parameters
  USE data_site       ! site specific data
  USE data_climate    ! climate data
  USE data_soil
  USE data_soil_cn
  USE data_stand           ! state variables of stand, cohort and cohort element
  USE data_out
  USE data_manag
  USE data_plant
  USE data_par
  IMPLICIT NONE

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' simulation_4C'

! allocation of environmental variable fields
 if(flag_wpm.ne.4 .and. flag_wpm.ne.5.and.flag_wpm.ne.6) then
  ! time loop
  DO time = 1, year
     iday = 1
  ! Update population variable for new year if population is changed through interventions 
     if (flag_standup .gt. 0 .or. flag_dis==1) then
        call stand_balance
        call standup
        flag_standup = 0
     endif

      CALL year_ini
      ! Calculate RedN from soil C/N
      ! read or create Redn for areal application 

     IF (time.EQ.1 .and. flag_redn) CALL RedN_ini
      IF (flag_dis .eq. 1) CALL dis_manag
  
      ! simulation of processes with subannual resolution (fluxes and soil)
      CALL stand_daily
      if(flag_end.ge.1) exit ! exit do loop time

      ! compressed output of start values
      IF (lcomp1) THEN                           
          CALL out_comp(unit_comp1)
          lcomp1 = .FALSE.
      ENDIF

      ! cohort litter production
      CALL senescence

      ! calculation of stand variables over all patches
       CALL stand_balance

      ! calculation of soil variables for yearly output
      CALL s_year

      ! calculation of fire variables for yearly output
      CALL fire_year

      ! calculation of indices for yearly output
      CALL t_indices(temp_mon)

      ! summation output
      IF(flag_sum.eq.4) THEN
         write(unit_sum,'(I5,9F11.3)') time_cur,photsum,npppotsum,nppsum,resosum,lightsum,nee,abslightsum,precsum, tempmean
         photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.;precsum=0.
      ENDIF

       totsteminc = 0.
       totsteminc_m3 = 0.
      ! cohort loop for change in crown dimensions, allocation and tree dimension calculations
      zeig=>pt%first
      DO
          IF (.not.ASSOCIATED(zeig)) exit
 
          IF (zeig%coh%height.ge.thr_height .and. zeig%coh%species.le.nspec_tree) then

              ! determine crown movement
              CALL CROWN( zeig )

              ! allocate NPP to the various tree compartments
              CALL PARTITION( zeig )
             if(flag_end.ge.1) exit  ! do loop
          ENDIF
          IF (zeig%coh%species.EQ.nspec_tree+1) then            ! Ground vegetation
              ! allocate NPP to the various ground vegetation compartments
              CALL PARTITION_SV( zeig )
          ENDIF
          IF (zeig%coh%species.eq.nspec_tree+2)  then             ! Mistletoe
              CALL PARTITION_MI( zeig )
              if(flag_end.ge.1) exit  ! do loop
          ENDIF
        zeig=>zeig%next
      END DO
      if(flag_end.ge.1) exit ! exit do loop time

 ! calculation of annual mortality
      IF(flag_mort.ge.1) CALL stand_mort

! annual growth of trees below thr_height, which are initialized by planting (not seeded!)
! at the beginning of the simulation or during management (shelter-wood)
        if(flag_reg.ne.2.and.flag_sprout.eq.0) CALL growth_seed
         CALL mort_seed
      if(flag_sprout.eq.1) flag_sprout=0
      IF(flag_mg==1) then
	     if(thin_year(act_thin_year)==time_cur) then
            CALL management
             act_thin_year = act_thin_year+1
		 end if
       ELSE IF((flag_mg.ge.2 .or. flag_mg.eq.3 .or. flag_mg.eq.33.or. flag_mg.eq.9 .or. flag_mg.eq.10).and.anz_spec.ne.0) THEN
 
         CALL management
         if(flag_wpm.ne.0) CALL timsort
       ENDIF

! no assortment if wpm is not called
      if(flag_mg.eq.0.and.anz_spec.ne.0) then
         if(flag_wpm.ne.0) call timsort
      end if
      CALL litter

! input of dead biomass into soil compartments
      CALL cn_inp

 !     if(flag_multi.eq.5) call m4c_simenv_out
! annual establishment for all species
      IF (flag_reg.eq.1.or.flag_reg.eq.2.or.flag_reg.eq.3.or.flag_reg.eq.30)  CALL stand_regen

! cumsteminc = cumsteminc + totsteminc   
! planting of seedlings/saplings at the beginning of simulation
      if(flag_reg.ge.9..and. flag_reg.lt.100. .and. time.eq.1) call planting
      if(flag_reg.ge.9..and. flag_reg.lt.100. .and. flag_mg .eq.44) call planting

! Update stand variables if stand changed
     if (flag_standup.gt.0 .or. anz_spec.eq.0) then
        call stand_balance
!        if (flag_standup .gt. 1) call root_distr   ! wird generell in year_ini berechnet
     endif
    cumsteminc = cumsteminc + totsteminc   

	

      ! yearly output
      IF (time_out .gt. 0) THEN
         IF (mod(time,time_out) .eq. 0) then
             CALL outyear (1)
             CALL outyear (2)
         endif
     ENDIF

      ! store of output variables (multi run 4, 8, 9)
      IF (nvar .gt. 1) CALL outstore

    ! RedN calculation
     if ((flag_limi .eq. 10) .or. (flag_limi .eq. 15)) call RedN_calc

!     CALL list_cohort
     CALL del_cohort

     if (.not.flag_mult910) PRINT *, ' * Year ', time, time_cur,' finished... '

   END DO ! time
    ! calculation of stand variables over all patches at the end!
      CALL stand_balance

         !*****    wpm    ******
	! check if management
	if(flag_mg == 0) then
		flag_wpm = 0
	endif	

	if (flag_wpm == 1 .or. flag_wpm == 21 .or. flag_wpm == 11) call wpm
	if (flag_wpm == 2) call sea
	if (flag_wpm == 3) then
		call wpm
		call sea
	end if
         !*** * * * * * * * ****
else
    call wpm
end if
   if (flag_wpm .gt. 0) call out_wpm(1)   

   CALL out_comp(unit_comp2)

   if(flag_end.eq.1) print*,ip, 'stop in partitio'
   if(flag_end.eq.3) print*,ip, 'stop in calc_la in canopy: toplayer = 125 m'
   flag_end = 0
   if (.not.flag_mult910) PRINT *, ' * Simulation ',ip,' finished.'

END SUBROUTINE simulation_4C


!**************************************************************
