//_____________________________________________________________
//  Simulation of stand structure change    ©  Takashi Kohyama
//   since 90/02/05 ~ for Sumatra Plots  
//
//  ver. 04/04 CONC.BAS
//  ver. 90/12/01   for 14 spp  S-MODEL   for Yakushima WTRF
//_____________________________________________________________
//  <<GAP MODEL>>: incorporating patch mosaic   29/Jul/1991 ~
//  GapM.bas for multi-species system           16/Aug/1991 ~ 
//  Publication: Kohyama T (1993) Journal of Ecology 81:131-143
//_____________________________________________________________
//  Version for SERIMBU DATA analysis           25/Jan/1999 ~
//  Translated to C++                           19/Mar/1999 ~
//  Completed by Matthew Potts                  22/Mar/1999 ~
//_____________________________________________________________
//  SAL model: "SAL" for (tree-)Size (patch-)Age (geographic-)Location 
//  Geographic-gradient/pach-mosaic combined   25/May/1999 ~
//  Kohyama('93) + Kohyama & Shigesada('95) hybrid model
//           Kohyama T & Shigesada N (1995) Vegetatio 121:117-126
//
//  KEY REFERENCE: 
//  Kohyama, T. (2005) Scaling up from shifting gap mosaic to 
//  geographic distribution in the modeling of forest dynamics.
//  Ecological Research 20 (3), in press.
//  
//  CONDITION FOR CODE USE
//  SAL code is for scientific use only, no commercial use;
//  References of Kohyama (2005), together with Kohyama (1993) and
//  Kohyama & Shigesada (1995) are requested.


#include <string.h>
#include <stdlib.h> 
#include <math.h>
#include <fstream.h> 
#include <iostream.h> 
#include <ctype.h> 

const int dpatchage = 16;
const int dsize_class = 60;
const int nspnn = 4;
const int g_site = 50;  // site resolution; the same as TEMA model (Kohyama&Shigesada '95)

const double PI_quart = 3.1416/4;

class patch_tree 
{	public:
	float tree_size[g_site][dpatchage][dsize_class];
};

class tree 
{	public:
	float size[dsize_class];
};

class patch
{	public:
	double patch_age[dpatchage];
};
	
class region
{	public:
	float site[g_site];
};

class Species
{	public:
	char name[10];
};

//class Data
//{	public:
//	float time;
//	float area;
//	float density;
//	float BA;
//     float biomass;
//};

// function decralations


int main()
{
// variable; parameter setting


int i, j, k, l, t,spn; //,numout;
int init, maxt, maxx, maxa, maxs, mins, time;
float dx, ds, maximt, maximx, minimx, maxima, mx, x, recr;
float Bt_sp_site;  // species BA at each region __ for print out
double dt, da, dadt;

char datainfile[100], dataoutfile[100];

spn = 3;        // three "species" (forest zone) system examined

dt = 1;                                  // time interval
dx = 4;                                  // size interval
da = 20;                                 // stand age interval
ds = 1;                                     // site interval
init = 0;                                // initial time
maximt = 30000;                           // maximum time/dt
maximx = 220;                            // maximum size
minimx = 0;                              // minimum size
maxima = 200;                              // maximum stand age/da
mx = 2;                                  // for mean x of class
maxs = 45;                            // max site location
mins = 6;                               // min site location

maxx=int(maximx/dx);
maxa=int(maxima/da);
maxt=int(maximt/dt);
dadt = da/dt;

Species plname[nspnn];

// dynamic equations variables, with initialization

// Data *Output;
// Output = (Data *) calloc(numout,sizeof(Data));



patch_tree *F;
F = (patch_tree *) calloc(nspnn, sizeof(patch_tree));

tree *F0, *Fa1;
F0 = (tree *) calloc(nspnn,sizeof(tree));
Fa1 = (tree *) calloc(nspnn,sizeof(tree));

patch *B0, *S, *S1;
B0 = (patch *) calloc(nspnn, sizeof(patch));
S = (patch *) calloc(g_site, sizeof(patch));
S1 = (patch *) calloc(g_site, sizeof(patch));

region *R, *sitef, *Blocal;
R =  (region *) calloc(nspnn, sizeof(region));
sitef = (region *) calloc(nspnn, sizeof(region));
Blocal = (region *) calloc(dpatchage, sizeof(region));

float *B, *Btot, *dens;
//double *S, *S1;
B = (float *) calloc(dsize_class,sizeof(float));
Btot = (float *) calloc(nspnn,sizeof(float));
dens = (float *) calloc(nspnn,sizeof(float));
//S = (double *) calloc(dpatchage,sizeof(double));
//S1 = (double *) calloc(dpatchage,sizeof(double));

double Sa = 0, Sa1 = 0, Srecr = 0;                // .. require a fine precision

// dynamics temporary variables (derivative abbreviation from above)

float freq, Fa, G, GFx, GFx1, dFdt, dFda, dGFdx, mc, surv;
double ms, dSdt;

// demographic parameters definition

float a[nspnn], a1[nspnn], a2[nspnn];            // growth parameters
float c[nspnn]; //, c1[nspnn],                    // mortality parameters
float d[nspnn], d1[nspnn];                    // recruitment parameters

// geographic dimensions________________________
float sopt[nspnn], spv[nspnn];

strcpy(plname[1].name, "TR");  //forest types
strcpy(plname[2].name, "WR");
strcpy(plname[3].name, "CR");

sopt[1]=40;                             //optimal site location
sopt[2]=30;
sopt[3]=20;

spv[1]=1.;                                //species vigor in optimal location
spv[2]=0.75;
spv[3]=0.5;

// dispersal parameter p
float dispf;
//dispf=.000025;                       //p = 1 km^2
dispf=.01;                                // p =100 km^2

//  Input demographic parameters from external file

strcpy(datainfile, "Basic_param.dat");       
ifstream infile("Basic_param.dat", ios::in);

//while (!infile.eof())
//{
j = 1;
//	for (j=1; j<=spn; j++) // (j=1; j<=spn; j++)
//	{
		//infile >> plname[j].name;
		//infile >> AA >> HH;
		infile >> a[j] >> a1[j] >> a2[j];
		infile >> c[j];
		infile >> d[j] >> d1[j];
//	}
//}
infile.close();

float rs;
ms = c[j];                                 // gap formation rate as asymptotic mortality

cout << "Parameter rs for adv. regen (0, .0005, 1000) ";
cin >> rs;
//rs = 0.0005;                                // coefficient of advance regeneration in new gap
                                               // rs>100 (say 1000); no adv reg. allowed, 
						// rs=0; no patch mosaic structure
// cout << "Coefficient of adv.regeneration?";
// cin >> rs

char ofilname[20];
cout << "Name of output file:  ";
cin >> ofilname;
strcpy(dataoutfile, ofilname);  //"Gmodel.out");

ofstream outfile(dataoutfile);
//      cout << AA << "\t" << HH << "\t" << a[j] << "\t" << a1[j] << "\t" << a2[j]<< "\t"  << c[j] << "\t" << c1[j] << "\t" << d[j] << "\t" << d1[j] << endl;

// Output simulation conditions prior to dynamics_loop

outfile << "running G_geogr.cpp:  Multispecies S model with gap mosaic along geograhic gradient " << endl;
outfile << "dt " << dt <<" yr" << endl;
outfile << "dl " << ds << "x 100 km" << endl;
outfile << "da " << da <<" yr" << endl;
outfile << "dx " << dx <<" cm" << endl;
outfile << " minus G not allowed" << endl;
//outfile << " D-H allom parameters:" << AA << HH << endl;
outfile << " species basic parameters (b, b1, b2, c, d): " << endl;

j=1; // for (j=1; j<=spn; j++)
//{
	outfile << a[j] << "\t" << a1[j] << "\t" << a2[j] << "\t" << c[j] << "\t"<< d[j] << "\t" << d1[j] << endl;
	cout << a[j] << "\t" << a1[j] << "\t" << a2[j] << "\t" << c[j] << "\t"<< d[j] << "\t" << d1[j] << endl;
//};
outfile << "dispersal parameter-p  " << dispf << endl;
outfile << "patch mosaic parameter-rs " << rs << endl;
cout << dispf << "\t" << rs << endl;

outfile << "   Start of run               " << endl;
outfile << endl;

//
// Fj(t,a,x) = s(t,a)*pj(t,a,x)   pj(t,a,x) is local density per unit area
//
// Initial condition: all at minimum patch age, .0001 tree per m^2 or 1 tree per ha

for (l=mins; l<=maxs; l++)
{
	S[l].patch_age[0] = 1;  
	for (j=1; j<=spn; j++)
	{
	        if (j==2)
		F[j].tree_size[l][0][1] = .0001; 
	        else
	        	F[j].tree_size[l][0][1] = 0.;
	}
	
	for (k=1; k<=maxa; k++)
	{
		S[l].patch_age[k] = 0;
		for (j=1; j<=spn; j++)
   		{
			for (i=2; i<=maxx; i++)
			{
     				F[j].tree_size[l][k][i] = 0.;
     			}
		}
	}
}

for (j=1; j<=spn; j++)
{
	for (i=1; i<=maxx; i++)
	{
		F0[j].size[i] = 0;
	}
}

//double SS = 0;          // to check whether area conserved at 1
float den = 0;               // for cumulation of tree number; per sq-m
float Bt = 0;                // for basal area cumulation

srand(100);

float seedy[nspnn];
float seed_year[nspnn];
float seed_prob = 0.25;        // seeding interval set at 4 years every spp.

for  (j=1; j<=spn; j++)
        {
	seedy[j] = 0.;
	seed_year[j] = 0.;
        }
//double h, W, Ws, Wb, Wl;     // these two lines for biomass calculation ...
//double biomass = 0.;          // ... based on Yamakura Seburu allometry

// time-course printout header line
outfile << "Time" <<"\t" << "Total area" << "\t" << "Total_density" <<"\t" <<"total_BA" << endl;

//  Dynamics main t loop
//_______________________________________________________________________


for (t=init; t<=maxt; t++)
{

	// seed year routine 
	for (j=1; j<=spn; j++){
		if (rand()/(1.+RAND_MAX) < seed_prob)
		{
			seed_year[j] = seedy[j];
			seedy[j] = 1.;
		}
		else 
		{
			seed_year[j] = 0.;
			seedy[j] += 1.;
		}
	}


/*
// no mast routine
for (j=1; j<=spn; j++){
seed_year[j]=1.;
}
*/

	//

	//cout << t << endl;
	
	// a century-long global warming definition
	if (t >= 10000/dt && t < 10100/dt)
	{
		for (j=1; j<=spn; j++)
      			sopt[j] -= 0.07*dt;
	}

	// time course printout
	time = (t*dt)/100;
	if ((t*dt == time*100) && (t*dt >= 9000) && (t*dt <= 30000))   // return change at t*10*dt yr interval
	{
		Bt = 0.;
		den = 0.;
      	  	//for (j=1; j<=spn; j++)
		for (j=2; j<=2; j++)
      	  	{
      	  		outfile << t*dt << "\t" ; //  << plname[j].name << "\t";
      	  		for (l=mins; l<=maxs; l++)
			{
				Bt_sp_site = 0;
				for (k=0; k<=maxa; k++)
				{
    	    				for (i=1; i<=maxx; i++)
      	    				{
          	    				x = (i - 1)*dx + minimx + mx;
          	       				den += F[j].tree_size[l][k][i];
          	        			Bt += x*x*F[j].tree_size[l][k][i];
          	        			Bt_sp_site +=  x*x*F[j].tree_size[l][k][i];
          	        	
  					};
				};
				outfile << Bt_sp_site*PI_quart << "\t";
 			};
 			outfile << endl;
 		} 
    		Bt = Bt*PI_quart;
		// outfile << endl;
		cout  << t*dt << "\t" << Bt/(maxs-mins+1) << endl;
  	}; // end of time-course printout section
  
	for (l=mins; l<=maxs; l++)
	{
		for (j=1; j<=spn; j++)
    		{
    			Btot[j] = 0;
    			sitef[j].site[l] = spv[j] - 0.004*(l-sopt[j])*(l-sopt[j]); // parabolic site-responce of vigor, K&S 95
    			if (sitef[j].site[l] < 0) 
    				sitef[j].site[l] = 0;
		};
  		Sa1 = Srecr = 0;
  	
		for (k=0; k<=maxa; k++)
		{
			for (j=1; j<=spn; j++)
			{
				B0[j].patch_age[k] = 0;
			}
    	
			//   Cumulative basal area routine
   	
   	 		for (i=maxx; i>=1; i--)
    			{
				B[i] = 0;
      				x = (i - 1)*dx + minimx + mx;
				for (j=1; j<=spn; j++)
				{
					freq = F[j].tree_size[l][k][i];
					
					if (S[l].patch_age[k]==0)
          					B0[j].patch_age[k]=0;
					else
						B0[j].patch_age[k] += x*x*freq*PI_quart/S[l].patch_age[k];
  					B[i] += B0[j].patch_age[k];               //  cumulative basal area at each k
				};
			};
    
			Blocal[k].site[l] = B[1];
    
			for (j=1; j<=spn; j++)
			{
				Btot[j] += B0[j].patch_age[k]*S[l].patch_age[k];
			}
	
			//    Stand age distribution routine
  	
			Sa = S[l].patch_age[k]/da;
			if (k==maxa)
				Sa = 0;
			ms = c[1];
			Srecr += ms*S[l].patch_age[k];
    
			//     Tree size distribution routine
   	
			for (j=1; j<=spn; j++)
			{
				for (i=1; i<=maxx; i++)
				{
					if (i==1)
        					GFx1 = 0;
					freq=F[j].tree_size[l][k][i];
                        	                       //  (1) size x dynamics
					x = (i - 1)*dx + minimx + mx;
        				G = a[1]*x*(1 - a1[1]*log(x) - a2[1]*B[i])*sitef[j].site[l];
        				if (G<0)                                    // negative G not allowed
        					G = 0;
        				GFx = G*freq/dx;
					dGFdx = GFx - GFx1;
                             	                  //   (2) age a dynamics
					Fa = freq/da;
        				if (k==maxa)
        					Fa = 0;
        				if (k==0)
        					Fa1[j].size[i] = 0;
					dFda = Fa - Fa1[j].size[i];
                         	                      //   (3) mortality by competition
					// fruction of survival with gap formation
         	
					surv = c[1] - rs*c[1]*x/(c[1] + rs*x);   // survived fruction
					mc = surv; // =c1[j]/x + surv;                   // non-gap mortality
				
					//        mc=cc[j]-(1-rs)*ms       // alternate Mc (constant as Koh '91b)
					//        mc=c+c1*a2[j]*B[i]/a2tot // alternate Mc (a2-proportional-Mc); The same as Thin-M model

							//   (4) gap formation
					F0[j].size[i] += surv*freq;   // summation of adv. regen. with respect to k
							//   (5) main dynamics equation
					dFdt = -dFda - dGFdx - ms*freq - mc*freq;
					F[j].tree_size[l][k][i] = freq + dFdt*dt;
        		
					Fa1[j].size[i] = Fa;
					GFx1 = GFx;
				};
			};
     
			//  Main equation for dynamics of patch age distribution
    
			dSdt = -Sa + Sa1 - ms*S[l].patch_age[k];
			S1[l].patch_age[k] = S[l].patch_age[k];
			S[l].patch_age[k] = S1[l].patch_age[k] + dSdt*dt;
			Sa1 = Sa;
		}; // end of k-loop
		
		// Boundary condition of patch age distribution
		S[l].patch_age[0] += Srecr*dt;

		for (j=1; j<=spn; j++)
		{    
			// definition of site-specific recruitment rate
			R[j].site[l] = seed_year[j]*d[1]*sitef[j].site[l]*Btot[j];//    *dt;  // repr. event-basis; *dt will be in later F summation
			if (R[j].site[l]<0)
				R[j].site[l]=0;

			//    F-oundary condition (2):  advance regeneration at a=0, supplied by gap formation
			for (i=1; i<=maxx; i++)
			{
				F[j].tree_size[l][0][i] += F0[j].size[i]*dt;
				F0[j].size[i] = 0;                         //  resetting of adv. regen. for next t
			};
		};
	} // l-loop end
	 
	//     F-boundary condition (1):  recruitment at x=minx in each stand age class, geographic dispersal involved
  			
	for (j=1; j<=spn; j++)
	{   	
		R[j].site[mins-1] = 0;          // boundary conditions for geographic gradient edges
		R[j].site[maxs+1] = 0;         // do
		
		for (l=mins; l<=maxs; l++)
		{
			for (k=0; k<=maxa; k++)
			{
				recr=R[j].site[l]+dispf*(R[j].site[l+1]-2*R[j].site[l]+R[j].site[l-1])/(ds*ds); // diffusion across sites
      				F[j].tree_size[l][k][1] += recr*exp(-d1[1]*Blocal[k].site[l])*S1[l].patch_age[k]*dt;       // dispersal process included
				//F[j].tree_size[l][k][1] += recr*S1[l].patch_age[k]*dt;    
			};
		};
	} //j
  
  
}  // end of t-loop

cout << "____END OF RUN____" << endl;

return 0;
};