#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>


/* NOTE: Now totalGrids is 70, from 144e-78w  */
   
/* the model of forced Rossby waves and their reflected Kelvin waves. Please see
documentation in Kelvin.c */


#define Boolean {True = 1; False = 0}
#define dx 220000.0
#define dt 86400
#define climSteps 365*86400/dt
#define totalSteps 5020*86400/dt   /* the number of days in this run */
#define totalGrids 70
#define const 0.0000000003675
#define drag0 1./(12.*2635200.)   /*  month number of drag friction */
#define initStep -9999
#define maxVerModes 8

float C[maxVerModes]= {2.73, 1.74, 1.06, 0.76, 0.59, 0.50, 0.42, 0.37};
float Chi[maxVerModes] = {4.29, 3.90, 1.70, 1.25, 1.83, 2.20, 1.57, 1.38};

float phaseSpeed()
{}
  
int effectMomentClim(float end, float origin, int currentStep, float cc)
{
  int steps, effectStep;

  steps = (int) (end - origin)/(dt*cc);
  /*
  printf("%10.0f%10.0f%4d \n", end, origin, steps); 
  */
  effectStep = currentStep - steps;
  
  while (effectStep < 0)
    effectStep = climSteps + effectStep;
  
  return effectStep;
}

int effectMoment(float end, float origin, int currentStep, float cc)
{
  int steps, effectStep;

  steps = (int) (end - origin)/(dt*cc);
  /*
  printf("%10.0f%10.0f%4d \n", end, origin, steps); 
  */
  effectStep = currentStep - steps;
  
  /* if the wave message has not got the current location, effect should be none, so
     set effectStep = -9999
  */

  if (effectStep < 0)
    effectStep = initStep;
  
  return effectStep;
}

void integral(FILE *ofp, int meriMode, int ii, float tauxc[totalGrids][climSteps], float taux[totalGrids][totalSteps], float westbound[totalSteps])
{
  int i, k, l, effectStep;
  float drag, response, damp;

  drag = drag0*(C[0]/C[ii])*(C[0]/C[ii]);

  k = 0;
  while (k < totalSteps)
    {
      for (i=0; i<totalGrids; i++)
	{
	  response = 0;
	  for (l=i; l<totalGrids; l++)
	    {
	      effectStep = effectMoment(i*dx, l*dx, k, C[ii]/(-(2*meriMode+1)));
	      damp = exp(-drag*(l-i)*dx*(2*meriMode+1)/C[ii]);
	      
	      if (effectStep == initStep)
		{
		  effectStep = effectMomentClim(i*dx, l*dx, k, C[ii]/(-(2*meriMode+1)));
		  response = response + tauxc[l][effectStep]*damp*dx*const*Chi[ii]/C[ii];
		}
	      else
		response = response + taux[l][effectStep]*damp*dx*const*Chi[ii]/C[ii];
	    }
	  
	  /* write to output file once every 5 days */
	  if ((k-2)%5 == 0)
	    fprintf(ofp, "%8.3f", response*100);

	  /* write western boundary for calculation of reflection */
	  if (i == 0)
	    westbound[k] = response;
	}

      if ((k-2)%5 == 0)    
	fprintf(ofp, "\n");
      
      ++k ;
    }
}

void integralClim(FILE *ofpc, int meriMode, int ii, float tauxc[totalGrids][climSteps], float westbound[climSteps])
{
  int i, k, l, effectStep;
  float drag, response, damp;

  drag = drag0*(C[0]/C[ii]);

  k = 0;
  while (k < climSteps)
    {
      for (i=0; i<totalGrids; i++)
	{
	  response = 0;
	  for (l=i; l<totalGrids; l++)
	    {
	      damp = exp(-drag*(l-i)*dx*(2*meriMode+1)/C[ii]);
	      effectStep = effectMomentClim(i*dx, l*dx, k, C[ii]/(-(2*meriMode+1)));
	      response = response + tauxc[l][effectStep]*damp*dx*const*Chi[ii]/C[ii];
	    }
	  
	  /* write to output file once every 5 days */
	  if ((k-2)%5 == 0)
	    fprintf(ofpc, "%8.3f", response*100);

	  /* write western boundary for calculation of reflection */
	  if (i == 0)
	    westbound[k] = response;
	}

      if ((k-2)%5 == 0)      
	fprintf(ofpc, "\n");
      
      ++k;
    }
}

void reflectClim(FILE *ofprc, float coeff, int ii, float westbound[climSteps])
{
  int i, k, effectStep;
  float drag, response, damp;

  drag = drag0*(C[0]/C[ii]);

  /* coeff is reflection factor */

  k = 2;
  while (k < climSteps)
    {
      for (i=0; i<totalGrids; i++)
	{
	  damp = exp(-drag*i*dx/C[ii]);
	  effectStep = effectMomentClim(i*dx, 0, k, C[ii]);
	  response = westbound[effectStep]*damp*coeff;
	  
	  fprintf(ofprc, "%8.3f", response*100);
	}
      
      fprintf(ofprc, "\n");
      
      k = k + 5;
    }
}

void reflect(FILE *ofpr, float coeff, int ii, float westbound[totalSteps], float westboundc[climSteps])
{
  int i, k, effectStep;
  float drag, response, damp;

  drag = drag0*(C[0]/C[ii]);

  /* coeff is reflection factor */

  k = 2;
  while (k < totalSteps)
    {
      for (i=0; i<totalGrids; i++)
	{
	  damp = exp(-drag*i*dx/C[ii]);
	  effectStep = effectMoment(i*dx, 0, k, C[ii]);
	  if (effectStep == initStep)
	    {
	      effectStep = effectMomentClim(i*dx, 0, k, C[ii]);	  
	      response = westboundc[effectStep]*damp*coeff;
	    }
	  else
	    response = westbound[effectStep]*damp*coeff;

	  fprintf(ofpr, "%8.3f", response*100);
	}
      
      fprintf(ofpr, "\n");
      
      k = k + 5;
    }
}


void main(int argc, char* argv[])
{
  int i,k,l, ii, totalVerModes, meriMode;
  FILE *ifpc, *ifp, *ofp, *ofpc, *ofpr, *ofprc, *input;
  float tauxc[totalGrids][climSteps], taux[totalGrids][totalSteps];
  float westb[totalSteps], westbc[climSteps];
  char *climWindFile, *windFile,*climOutFile, *outFile, *refClimOutFile, *refOutFile;

  /* open input data file */
  if(!(input = fopen(argv[1], "r")))
    {
      printf("open inputfile error: %s!! \n", argv[1]);
      exit(1);
    }

  /* read in input data */
  fscanf(input, "%d,%d \n", &totalVerModes,&meriMode);

  climWindFile = calloc(40, sizeof(char));
  windFile = calloc(40, sizeof(char));
  climOutFile = calloc(40, sizeof(char));
  outFile = calloc(40, sizeof(char));
  refClimOutFile = calloc(40, sizeof(char));
  refOutFile = calloc(40, sizeof(char));

  for (ii=0; ii<totalVerModes; ii++)
    {
      /* input file names in the input data file */
      fscanf(input, "%s \n", climWindFile);
      fscanf(input, "%s \n", windFile);
      fscanf(input, "%s \n", climOutFile);
      fscanf(input, "%s \n", outFile);
      fscanf(input, "%s \n", refClimOutFile);
      fscanf(input, "%s \n", refOutFile);

      /*
	Check if input files can be accessed
      */
      if(!(ifpc = fopen(climWindFile, "r")))
	{
	  printf("open file error: %s!! \n", climWindFile);
	  exit(1);
	}
      if(!(ifp = fopen(windFile, "r")))
	{
	  printf("open file error: %s!! \n", windFile);
	  exit(1);
	}
      if(!(ofpc = fopen(climOutFile, "w")))
	{
	  printf("open file error: %s!! \n", climOutFile);
	  exit(1);
	}
      if(!(ofp = fopen(outFile, "w")))
	{
	  printf("open file error: %s!! \n", outFile);
	  exit(1);
	}
      if(!(ofprc = fopen(refClimOutFile, "w")))
	{
	  printf("open file error: %s!! \n", refClimOutFile);
	  exit(1);
	}
      if(!(ofpr = fopen(refOutFile, "w")))
	{
	  printf("open file error: %s!! \n", refOutFile);
	  exit(1);
	}

      /* Read in Taux_clim and time series */
      for(k=0; k<climSteps; k++)
	for (i=0; i<totalGrids; i++)
	  fscanf(ifpc, "%f", &tauxc[i][k]);

      for(k=0; k<totalSteps; k++)
	for (i=0; i<totalGrids; i++)
	  fscanf(ifp, "%f", &taux[i][k]);

      integralClim(ofpc, meriMode, ii, tauxc, westbc);      
      integral(ofp, meriMode, ii, tauxc, taux, westb);

      /* calculate reflection */
      if (meriMode == 1)
	{
	  reflectClim(ofprc, 0.5, ii, westbc);
	  reflect(ofpr, 0.5, ii, westb, westbc);
	}
      if (meriMode == 3)
	{
	  reflectClim(ofprc, 0.1, ii, westbc);
	  reflect(ofpr, 0.1, ii, westb, westbc);
	}
      if (meriMode == 5)
	{
	  reflectClim(ofprc, 0.05, ii, westbc);
	  reflect(ofpr, 0.05, ii, westb, westbc);
	}


      /* close all the file for a mode */
      fclose(ifpc);
      fclose(ifp);
      fclose(ofpc);
      fclose(ofp);
      fclose(ofprc);
      fclose(ofpr);
    }    
}