!*****************************************************************!
!*                                                               *!
!*           SEA for 4C (FORESEE) Simulation Model               *!
!*                                                               *!
!*                                                               *!
!*						Subroutines:		                     *!
!*                                                               *!
!*  sea: control subroutine for sea                              *!
!*                                                               *!
!*  sort_mansort: first sorting of mansort                       *!
!*  sort_standsort: first sorting of standsort                   *!
!*  sort_industrial                                              *!
!*  calculate_harvest_costs                                      *!
!*  calculate_assets                                             *!
!*  calculate_costs                                              *!
!*  calculate_npv                                                *!
!*  read_sea_prices                                              *!
!*                                                               *!
!*                  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 sea

use data_simul
use data_wpm

implicit none
	

character(150) pricesFile

	
	! begin program
	call setFlags

	call allocate_in_output
	
	call ini_input_sea	

	
	! read prices
	pricesFile = trim(dirin)//'sea_prices.wpm'
    ! call fullPath(pricesFile, dirin)
	call read_sea_prices(pricesFile)

	! simulation
	if ( associated(first_mansort) ) then
		! first sorting
		call sort_mansort
	end if

	if ( associated(first_standsort) ) then
		call sort_standsort
	end if

	! sort 0.4 to the industrial wood
	call sort_industrial

	! harvest costs calculation
	call calculate_harvest_costs
	
	! timber selling assets calculation
	call calculate_assets
	
	! calculate rest costs
	call calculate_costs
	
	! calculate npv
	call calculate_npv



end subroutine sea


!***************************************************************
! calculate timber grades from the mansort and standsort input
! input:  data_wpm
! output: data_wpm
!***************************************************************

subroutine sort_mansort

use data_wpm
use data_simul

implicit none

integer i, index
real volume,  pi, diam
character(4) act_typus
integer act_spec, act_year, set_year
	

	pi = 3.1415926536  ! PI

    i = nr_years

	! set the first year, set an ima
	act_year = first_mansort%mansort%year
	set_year = first_manrec%manrec%year

	! set the run pointer to the begin of the list
	act_mansort => first_mansort
	act_manrec	=> first_manrec
	
	! check if actuelles management is not brushing or tending
	if (trim(act_manrec%manrec%management) .eq. 'brushing') then
		if( associated(act_manrec%next) ) then
			act_manrec => act_manrec%next
		endif
	endif
			
	if (trim(act_manrec%manrec%management) .eq. 'tending') then
		if( associated(act_manrec%next) ) then
			act_manrec => act_manrec%next
		endif
	endif



	! check if last management year was some years ealier
	do while (set_year < act_mansort%mansort%year)
		act_mansort => act_mansort%next	
	end do			

	act_year = act_mansort%mansort%year
	set_year = act_manrec%manrec%year


	if( associated(act_manrec%next) ) then
		act_manrec => act_manrec%next
	end if

	do while (associated(act_mansort))
		
		if ( act_year <= act_manrec%manrec%year .and.		&
				act_manrec%manrec%year /= set_year ) then
			! check the management in actual year
			if (trim(act_manrec%manrec%management) .eq. 'brushing') then
				if( associated(act_manrec%next) ) then
					act_manrec => act_manrec%next
				endif
			endif
			
			if (trim(act_manrec%manrec%management) .eq. 'tending') then
				if( associated(act_manrec%next) ) then
					act_manrec => act_manrec%next
				endif
			endif
			
			if (trim(act_manrec%manrec%management) .ne. 'tending'	.and.	&
				trim(act_manrec%manrec%management) .ne. 'brushing' ) then
					! set next value for actual manrec year
					set_year = act_manrec%manrec%year
			endif
			
			if( associated(act_manrec%next) ) then
				act_manrec => act_manrec%next
			end if
		endif 		               
		
	    	! set species index
		act_spec	= act_mansort%mansort%spec
		act_typus	= act_mansort%mansort%typus

		! calculate carbon for the actual line in mansort
		volume = act_mansort%mansort%volume * act_mansort%mansort%number 	
		

		select case (trim(act_typus))
			case ('ste1', 'ste2')				
				diam = act_mansort%mansort%diam_wob
				if (diam >=25 .and. diam < 30)	index = 8
				if (diam >=30 .and. diam < 35)	index = 9
				if (diam >=35)			index = 10					

			case ('sg1', 'sg2')
				diam = act_mansort%mansort%diam_wob
				if (diam < 15)					index = 3
				if (diam >=15 .and. diam < 20)	index = 4
				if (diam >=20 .and. diam < 25)	index = 5
				if (diam >=25 .and. diam < 30)	index = 6
				if (diam >=30)					index = 7					

			case ('in1', 'in2')
				index = 2
				case ('fue')
				index = 1
						
			end select
	
			mansort_tg(act_spec, index, set_year)	= mansort_tg(act_spec, index, set_year)	+ volume

		! after using the mansort list item, go to the next
		act_mansort => act_mansort%next
	
		if (associated(act_mansort)) then
			act_year = act_mansort%mansort%year		
		end if
	end do


end subroutine sort_mansort

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

subroutine sort_industrial

use data_wpm
use data_simul

implicit none

integer i, j
real ind
		
	! sort value*ind into the industrial wood
	ind = 0.4

	do i = 1, nr_spec
		do j = 1, nr_timb_grades
			select case (j>2)
			case (.TRUE.)
				mansort_tg(i, 2, :)   = mansort_tg(i, j, :)   * ind + mansort_tg(i, 2, :)			
				mansort_tg(i, j, :)   = mansort_tg(i, j, :)   * (1 - ind)
				standsort_tg(i, 2, :) = standsort_tg(i, j, :) * ind + standsort_tg(i, 2, :)			
				standsort_tg(i, j, :) = standsort_tg(i, j, :) * (1 - ind)
			end select
		end do
	end do


end subroutine sort_industrial

!****************************************************************
subroutine sort_standsort

use data_wpm

implicit none

integer j, index
real volume,  pi, diam
character(4) act_typus
integer act_spec, act_year
	

	pi = 3.1415926536  ! PI

	j = 1
	! set the first year
	act_year = first_standsort%mansort%year
	! set the run pointer to the begin of the list
	act_standsort => first_standsort
	
	! check if last management year was some years ealier
	do while (act_year < act_standsort%mansort%year)
		act_standsort => act_standsort%next	
	end do			

	do while (associated(act_standsort))
		
		! check the management year
		if ( act_year <= act_standsort%mansort%year )then
			! set next value for actual standsort year
			j = j+1
		endif
		
		! set species index
		act_spec	= act_standsort%mansort%spec

		act_typus	= act_standsort%mansort%typus
		! calculate carbon for the actual line in standsort
		volume = act_standsort%mansort%volume * act_standsort%mansort%number
		
		
		select case (trim(act_typus))
			
			case ('ste1', 'ste2')				
				diam = act_standsort%mansort%diam_wob
				if (diam >=25 .and. diam < 30)	index = 8
				if (diam >=30 .and. diam < 35)	index = 9
				if (diam >=35)					index = 10					


			case ('sg1', 'sg2')
				diam = act_standsort%mansort%diam_wob
				if (diam < 15)					index = 3
				if (diam >=15 .and. diam < 20)	index = 4
				if (diam >=20 .and. diam < 25)	index = 5
				if (diam >=25 .and. diam < 30)	index = 6
				if (diam >=30)					index = 7					

			case ('in1', 'in2')
				index = 2

			case ('fue')
				index = 1
							
		end select
	
		
		standsort_tg(act_spec, index, act_standsort%mansort%year)	= standsort_tg(act_spec, index, act_standsort%mansort%year)	+ volume
		

		! after using the standsort list item, go to the next
		act_standsort => act_standsort%next
		
		if (associated(act_standsort)) then
			act_year = act_standsort%mansort%year		
		end if

	end do


end subroutine sort_standsort


!*****************************************************************************
subroutine calculate_harvest_costs()

use data_wpm
use data_simul

implicit none

integer i, j

! calcultation of costs only implemented for monoculture stands 
! differentiation between coniferous and deciduous trees
do i = 1, nr_spec
	if (nr_spec.eq.2 .or. nr_spec.eq.3) then
		do j = 1, nr_timb_grades
			ms_costs(i,:) = hsystem(2) * mansort_tg(i,j,:) * chainsaw_prices(i,j) +		&
					hsystem(1) * mansort_tg(i,j,:) * harvester_prices(i,j) +	&
					ms_costs(i,:)
			st_costs(i,:) = hsystem(2) * standsort_tg(i,j,:) * chainsaw_prices(i,j) +	&
					hsystem(1) * standsort_tg(i,j,:) * harvester_prices(i,j) +	&
					st_costs(i,:)
		end do
	else
		do j = 1, nr_timb_grades
			ms_costs(i,:) = mansort_tg(i,j,:) * chainsaw_prices(i,j) + ms_costs(i,:)
			st_costs(i,:) = standsort_tg(i,j,:) * chainsaw_prices(i,j) + st_costs(i,:)
		end do
	endif
end do

end subroutine calculate_harvest_costs



!*****************************************************************************
subroutine calculate_assets()

use data_wpm
use data_simul

implicit none

integer i, j

	do i = 1, nr_spec
		do j = 1, nr_timb_grades
			ms_assets(i,:) = mansort_tg(i,j,:) * net_prices(i,j) +		&
							ms_assets(i,:)
			st_assets(i,:) = standsort_tg(i,j,:) * net_prices(i,j) +	&
							st_assets(i,:)
		end do
	end do

end subroutine calculate_assets



!*****************************************************************************
subroutine calculate_costs()

use data_wpm
use data_simul
use data_plant

implicit none

character(30) manag
integer i, act_year, spec

	
	! sum of standsort
	do i = 1, nr_spec
		sum_costs(2,:) = st_assets(i,:) - st_costs(i,:) + sum_costs(2,:)
	end do

	! sum of mansort
	do i = 1, nr_spec
		sum_costs(3,:) = ms_assets(i,:) - ms_costs(i,:) + sum_costs(3,:)
	end do

	! silvicultural costs like tending etc.
	act_manrec	=> first_manrec
	do while (associated(act_manrec))
		act_year = act_manrec%manrec%year
		manag  = trim(act_manrec%manrec%management)
		
		! sum up silvicultural costs and subsidies
		select case (trim(manag))
			
			case ('tending')
				! tending costs and subsidies
				sum_costs(4, act_year) = - tending_prices(1) + sum_costs(4, act_year)
				subsidy(2, act_year) = tending_prices(2) + subsidy(2, act_year)
						
			case ('brushing')
				! brushing costs and subsidies
				sum_costs(4, act_year) = - brushing(1) + sum_costs(4, act_year)
				subsidy(2, act_year) = brushing(2) + subsidy(2, act_year)

			case ('felling')
				! forest maintenance costs and subsidies
				sum_costs(4,act_year) = - ext_for(2,1) + sum_costs(4,act_year)
				sum_costs(5,act_year) = ext_for(2,2) + sum_costs(5,act_year)
				
			case ('shelterwood system1')
				! forest maintenance costs and subsidies
				sum_costs(4,act_year) = - ext_for(1,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(1,2) + sum_costs(4,act_year)
				sum_costs(4,act_year) = - ext_for(2,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(2,2) + sum_costs(4,act_year)

			case ('shelterwood system2')
				! brushing, forest maintenance, timber selling costs and subsidies
				sum_costs(4,act_year) = - ext_for(1,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(1,2) + sum_costs(4,act_year)
				sum_costs(4,act_year) = - ext_for(2,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(2,2) + sum_costs(4,act_year)
				sum_costs(4, act_year) = - brushing(1) + sum_costs(4, act_year)
				subsidy(2, act_year) = brushing(2) + subsidy(2, act_year)

			case ('felling after shelterwood s.')
				sum_costs(4,act_year) = - ext_for(2,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(2,2) + sum_costs(4,act_year)

			case ('thinning')
				! forest maintenance, timber selling
				sum_costs(4,act_year) = - ext_for(1,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(1,2) + sum_costs(4,act_year)
				sum_costs(4,act_year) = - ext_for(2,1) + sum_costs(4,act_year)
				sum_costs(4,act_year) = ext_for(2,2) + sum_costs(4,act_year)
			
		
		end select
		act_manrec => act_manrec%next
	end do

	! planting
	if (plant_year /= 0) then
		select case (flag_plant)
			case(8,7,6,5,4,33)
				sum_costs(4,plant_year)	= - sum(planting_prices * (npl_mix / 1000)) + sum_costs(4,plant_year)
				subsidy(2,plant_year)	= sum(planting_sub(1,plant_year)*npl_mix/npl_mix) + subsidy(2,plant_year)
							
			! pine, beech, oak, spruce, birch
			case (10)
				spec = 3
			case (11)
				spec = 1
			case (12)
				spec = 4
			case (13)
				spec = 2
			case (14)
				spec = 5
		end select
		sum_costs(4,plant_year)	= - planting_prices(spec) * numplant(spec)/1000 + sum_costs(4,plant_year)
		subsidy(2,plant_year)	= planting_sub(1,spec) + subsidy(2,plant_year)

		! fence
		sum_costs(4,plant_year) = - fence(1,spec) + sum_costs(4,plant_year)
		sum_costs(5,plant_year) = fence(2,spec) + sum_costs(5,plant_year)
	end if


	! sum up subsidies
	sum_costs(5,:) = subsidy(1,:) + subsidy(2,:) + fix(2) + sum_costs(5,:)

	! sum up all except standsort
	sum_costs(1,:) = sum_costs(3,:) + sum_costs(4,:) - fix(1) + sum_costs(5,:)

end subroutine calculate_costs


!*****************************************************************************
subroutine calculate_npv()

use data_wpm
use data_simul
use data_plant

implicit none

real, dimension(4, nr_years) :: rate
integer i, j

	rate(:,:) = 0.
	do i = 1, nr_years		
		do j = 1, 4
			rate(j, i) = (1+int_rate(j))**i
			npv(j,i) = (sum_costs(2,i) + sum_costs(3,i)) / rate(j, i)
		npv(j+4,i) = sum(sum_costs(1,1:i)/rate(j,1:i))
			npv(j+8,i) = npv(j+4,i) - npv(1,1) + npv(j,i)
		end do
	end do

end subroutine calculate_npv


!*****************************************************************************
subroutine read_sea_prices(pricesFile)

use data_wpm
use data_simul

implicit none

character(70) pricesFile
integer i,  unit, ios

	unit = getunit()

	open(unit,	FILE=pricesFile, STATUS='OLD', action='read')
	
	! Headers
	do i=1,5
      read (unit, *)
    enddo

	read(unit,	'(F6.2)',iostat=ios) fix(1)
    read (unit, *)
	read(unit, *)
	read(unit,	'(5(F6.2))',iostat=ios) &
		planting_prices(1),			&
		planting_prices(2),			&
		planting_prices(3),			&
		planting_prices(4),			&
		planting_prices(5)
	read(unit, *)
	read(unit,	'(5(F6.2))',iostat=ios)    &
		fence(1,1),						&
		fence(1,2),						&
		fence(1,3),						&
		fence(1,4),						&
		fence(1,5)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) brushing(1)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) tending_prices(1)

	read(unit, *)
	read(unit, *)
	read(unit,	'(2(F6.2))',iostat=ios) &
		hsystem(1), &
		hsystem(2)
	
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) dec_per

	do i=1,4
      read (unit, *)
    enddo
    do i=1, nr_timb_grades
		read(unit,	'(5(F6.2) )',iostat=ios)		&
				chainsaw_prices(1, i),			&
				chainsaw_prices(2, i),			&
				chainsaw_prices(3, i),			&
				chainsaw_prices(4, i),			&
				chainsaw_prices(5, i)
	end do

	read(unit, *)
	read(unit, *)	
	do i=1, nr_timb_grades
		read(unit, *)	&
				harvester_prices(1, i),	&
				harvester_prices(2, i), &
				harvester_prices(3, i), &
				harvester_prices(4, i), &
				harvester_prices(5, i)
	end do

	read(unit, *)
	read(unit, *)
	do i=1, nr_timb_grades
		read(unit,	'(5(F6.2) )',iostat=ios)	&
			net_prices(1, i),					&
			net_prices(2, i),					&
			net_prices(3, i),					&
			net_prices(4, i),					&
			net_prices(5, i)
	end do

	read(unit, *)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) ext_for(1,1)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) ext_for(1,2)

	read(unit, *)
	read(unit, *)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) fix(2)
    read (unit, *)
	read(unit, *)
	read(unit,	'(5(F6.2))',iostat=ios) &
		planting_sub(1,1),			&
		planting_sub(1,2),			&
		planting_sub(1,3),			&
		planting_sub(1,4),			&
		planting_sub(1,5)
	read(unit, *)
	read(unit,	'(5(F6.2))',iostat=ios) &
		planting_sub(2,1),			&
		planting_sub(2,2),			&
		planting_sub(2,3),			&
		planting_sub(2,4),			&
		planting_sub(2,5)
	read(unit, *)
	read(unit,	'(5(F6.2))',iostat=ios)    &
		fence(1,1),						&
		fence(1,2),						&
		fence(1,3),						&
		fence(1,4),						&
		fence(1,5)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) brushing(2)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) tending_prices(2)

	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) ext_for(2,2)
	read(unit, *)
	read(unit,	'(F6.2)',iostat=ios) ext_for(2,2)

	read(unit, *)
	read(unit, *) &
		int_rate(2),					&
		int_rate(3),					&
		int_rate(4)

	close(unit, status="keep")

end subroutine read_sea_prices
