### Rämö Janne ###
# Bollandsås et al. (2008) growth model, age structured #

### Species, set in .dat ###
param sp;

## Amount of age classes, default 12 ###
param ac;

### Recruit-function, probability ###
param rp {1..4, 1..4};

### Recruit-function parameters, conditional ###
param rc {1..4, 1..4};

### Diameter increment -function parameters ###
param di {1..8, 1..4};

### Mortality-function parameters ###
param mp {1..4, 1..4};

param maxt;							#Max time
set T := 0..maxt by 5;
set pl := sp..sp;
param q := 0.5;

param r;
param R := 1/(1+r);	
param s;
param SI;									#Site index
param pi = 4 * atan(1); 					#Set the value of pi to parameter
param LAT;									#Latitude
param y0 {i in 1..ac, j in 1..4} default 0;			#Initial stand, number of trees
var DBH {t in T, i in 1..ac, j in pl} >= 0;
param DBH0 {i in 1..ac} default 0;

param ps {j in 1..4};						#Tukin hinta, kantohinta
param pp {j in 1..4};						#Kuidun hinta, kantohinta


var y {t in T, i in 1..ac, j in pl} >= 0;	#Stand development parameter
var H {t in T, i in 1..ac, j in pl} >= 0;
param Hbool {t in T} default 0;

var BAsc {t in T, i in 1..ac, j in pl} = pi*(DBH[t,i,j]/2/1000)^2*y[t,i,j];		#Basal area of individual size classes, array, calculated from y & DBH
var BAsp {t in T, i in pl} = sum {e in 1..ac} BAsc[t,e,i];					#Basal area of species
var BA {t in T} = sum {u in pl} BAsp[t,u];						#Basal area, m^2/ha
var BAL {t in T, i in 1..ac} = sum {e in (i+1)..ac, u in pl} BAsc[t,e,u];		#Basal area of larger trees, m^2/ha
var PBA {t in T, i in pl} = BAsp[t,i]/BA[t]*100;								#Percentage of Basal area for the subject species

#Uncomment the species-sitetype combination that is optimized, estimated from the data used in size-structured model.
# Spruce volume functions, SI=15 
var Vs {t in T, i in 1..ac, j in pl} = (116.0906-31.1854*(DBH[t,i,j]/10)+1.9407*(DBH[t,i,j]/10)^2-0.0121*(DBH[t,i,j]/10)^3)/1000;
var Vp {t in T, i in 1..ac, j in pl} = (0.0068176*(DBH[t,i,j]/10)^3-0.660699*(DBH[t,i,j]/10)^2+18.2853*(DBH[t,i,j]/10)-72.8905)/1000;

# Pine volume functions, SI=15 
# var Vs {t in T, i in 1..ac} = ((-32.777+3623.353/(1+(DBH[t,i,3]/10/48.547)^-3.256))/1000;  
# var Vp {t in T, i in 1..ac} = (24.954+110.575/(1+((DBH[t,i,3]/10-15.797)/12.562)^2))/1000; 

# Pine  volume functions, SI=11
# var Vs {t in T, i in 1..ac} = (-31.764+3557.296/(1+(DBH[t,i,j]/10/49.073)^-3.233))/1000;
# var Vp {t in T, i in 1..ac} = (24.6+107.658/(1+((DBH[t,i,j]/10-15.771)/2.470)^2))/1000;


var standVol {t in T, j in pl} = sum {i in 1..ac} (y[t,i,j]*Vp[t,i,j]+y[t,i,j]*Vs[t,i,j]);

var Hs {t in T, i in 1..ac, j in pl} = H[t,i,j]*Vs[t+5,i+1,j]; #Harvests occur at the end of the period
var Hp {t in T, i in 1..ac, j in pl} = H[t,i,j]*Vp[t+5,i+1,j];


var m {t in T, i in 1..ac, j in pl} =  (1+exp(-(mp[1,j]+mp[2,j]*DBH[t,i,j]+mp[3,j]*10^-5*DBH[t,i,j]^2+mp[4,j]*BA[t])))^-1;		

var I {t in T, i in 1..ac, j in pl} = di[1,j]+di[2,j]*DBH[t,i,j]+di[3,j]*10^-5*DBH[t,i,j]^2+di[4,j]*10^-8*DBH[t,i,j]^3+di[5,j]*BAL[t,i]+di[6,j]*SI+di[7,j]*BA[t]+di[8,j]*LAT;  # ) <= 0 then 0 else (di[1,j]+di[2,j]*DBH[t,i,j]+di[3,j]*10^-5*DBH[t,i,j]^2+di[4,j]*10^-8*DBH[t,i,j]^3+di[5,j]*BAL[t,i]+di[6,j]*SI+di[7,j]*BA[t]+di[8,j]*LAT);

var CR {t in T, i in pl} = rc[1,i]*BA[t]^rc[2,i]*SI^rc[3,i]*PBA[t,i]^rc[4,i];						#Recruitment, combability
var Rprob {t in T, i in pl} = (1+exp(-(rp[1,i]+rp[2,i]*BA[t]+rp[3,i]*SI+rp[4,i]*PBA[t,i])))^-1;		#Recruitment, probability
var d {t in T, i in pl} = CR[t,i]*Rprob[t,i];														#Amount recruited

### Total cubic meters of wood harvested per size class ###
var Vtotal {t in T, j in pl} = sum {i in 1..ac} H[t,i,j]*(Vs[t+5,i+1,j]+Vp[t+5,i+1,j]); #Harvests occur at the end of the period
### Total number of trees harvested per size class ###
var Htotal {t in T, j in pl} = sum {i in 1..ac} H[t,i,j];
### Total amount of trees per size class ###
var ytotal {t in T, j in pl} = sum {i in 1..ac} y[t,i,j];


### Revenues ###
var c {t in T} = sum {i in 1..ac, j in pl} (Hs[t,i,j]*ps[j]+Hp[t,i,j]*pp[j]);


### Set the initial state ###
subject to initial_state {i in 1..ac, j in pl}:
	y[0,i,j] = y0[i,j];
	
subject to initialDBH {i in 2..ac, j in pl}:
  DBH[0,i,j] = DBH0[i];
  
subject to diameter1 {t in T, j in pl}:
  DBH[t,1,j] = 50;
 
subject to diameter {t in 0..maxt-5 by 5, i in 1..ac-1, j in pl}:
  DBH[t+5,i+1,j] = DBH[t,i,j]+I[t,i,j];


subject to harvests3 {t in T, i in 1..ac, j in pl}:
  H[t,i,j] = Hbool[t]*H[t,i,j];
  
subject to standstate1 {t in 0..maxt-5 by 5, j in pl}:
  y[t+5,1,j] = d[t,j]-H[t,1,j];

subject to standstate {t in 0..maxt-5 by 5, i in 1..ac-2, j in pl}:
  y[t+5,i+1,j] = (1-m[t,i,j])*y[t,i,j]-H[t,i+1,j];
  
subject to standstateF {t in 0..maxt-5 by 5, j in pl}:
  y[t+5,ac,j] = (1-m[t,ac-1,j])*y[t,ac-1,j]+(1-m[t,ac,j])*y[t,ac,j]-H[t,ac,j];
  
subject to ylimiter {t in T, i in 1..ac, j in pl}:
	y[t,i,j] <= 300;
	
	
  
maximize objective:
	(sum {t in 0..maxt-5 by 5} (c[t]*R^t/1000));
  
param i;