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


/*********************Description******************** 

The integration of forced Kelvin wave and reflected Rossby waves along the wave
characteristic paths. The reflection is based on Cane and Sarachik (1981) that
meridionally integrated zonal transport is zero. The code is not well written, 
but I try to document it well-:)

******************************************************/

/* NOTE: Now totalGrids is 70, from 144e-78w  */
   

#define Boolean {True = 1; False = 0}
#define dt 86400                   /* seconds a day */
#define climSteps 365*86400/dt     /* always assume 365 days a year */
#define totalSteps 5020*86400/dt   /* the number of days in this run */
#define const 0.0000000003675      /* RHOair*drag/D */
#define drag0 1./(12.*2635200.)    /*  month number of drag friction. 
				      NOTE: this actually is damping */
#define initStep -9999
#define maxVerModes 8
#define totalGrids 70              /* total zonal grid points */
#define dx 220000.0                /* grid size */ 

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};

/*never designed this. My original purpose is to adjust Cg with latitudes.*/
float phaseSpeed()
{}

/*find effective timestep along the characteric path. Since seasonal cycle is
cyclic, adding 12 months if the effective timestep becomes negative */  
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;
}

/*find effective timestep along the characteric path for whole time series. seasonal cycle is
used to initialize a run. That is, add 12 months to the effective timestep if it becomes
negative */  
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;
}

/* Actual integration along the path for whole time series*/
void integral(FILE *ofp, int ii, float tauxc[totalGrids][climSteps], float taux[totalGrids][totalSteps], float eastbound[totalSteps])
{
  int i, k, l, effectStep;
  float drag, response, damp;

  /* dumping depends on vertical modes based on McCreary (1981)*/
  drag = drag0*(C[0]/C[ii])*(C[0]/C[ii]);

  k = 0;
  while (k < totalSteps)
    {
      /*look at every zonal grid */
      for (i=0; i<totalGrids; i++)
	{
	  response = 0;
	  /* for the grid i, only grids west of it matters */
	  for (l=0; l<=i; l++)
	    {
	      effectStep = effectMoment(i*dx, l*dx, k, C[ii]);
	      damp = exp(-drag*(i-l)*dx/C[ii]);
	      
	      /*if effective step is negative, use mean clim to initialize */
	      if (effectStep == initStep)
		{
		  effectStep = effectMomentClim(i*dx, l*dx, k, C[ii]);
		  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 == totalGrids-1)
	    eastbound[k] = response;
	}

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

/* Actual integration along the path for seasonal cycle*/
void integralClim(FILE *ofpc, int ii, float tauxc[totalGrids][climSteps], float eastbound[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=0; l<=i; l++)
	    {
	      damp = exp(-drag*(i-l)*dx/C[ii]);
	      effectStep = effectMomentClim(i*dx, l*dx, k, C[ii]);
	      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 == totalGrids-1)
	    eastbound[k] = response;
	}

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

/* this should be easy to understand. parameters ii and coeff are the vertical mode number
and reflection coefficient */
void reflectClim(FILE *ofprc, int meriMode, int ii, float coeff, float eastbound[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*(totalGrids-1-i)*dx*(2*meriMode+1)/C[ii]);
	  effectStep = effectMomentClim(i*dx, (totalGrids-1)*dx, k, C[ii]/(-(2*meriMode+1)));
	  response = eastbound[effectStep]*damp*coeff;
	  
	  fprintf(ofprc, "%8.3f", response*100);
	}
      
      fprintf(ofprc, "\n");
      
      k = k + 5;
    }
}

/* this should be easy to understand. parameters ii and coeff are the vertical mode number
and reflection coefficient. Again reflected seasonal cycle Rossby waves are for initializing
whole time series of reflected Rossby waves */
void reflect(FILE *ofpr, int meriMode, int ii, float coeff, float eastbound[totalSteps], float eastboundc[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*(totalGrids-1-i)*dx*(2*meriMode+1)/C[ii]);
	  effectStep = effectMoment(i*dx, (totalGrids-1)*dx, k, C[ii]/(-(2*meriMode+1)));
	  if (effectStep == initStep)
	    {
	      effectStep = effectMomentClim(i*dx, (totalGrids-1)*dx, k, C[ii]/(-(2*meriMode+1)));	  
	      response = eastboundc[effectStep]*damp*coeff;
	    }
	  else
	    response = eastbound[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;
  FILE *ifpc, *ifp, *ofp, *ofpc, *ofprM1, *ofprcM1,*ofprM3, *ofprcM3, *input;
  float tauxc[totalGrids][climSteps], taux[totalGrids][totalSteps];
  float eastb[totalSteps], eastbc[climSteps];
  char *climWindFile, *windFile,*climOutFile, *outFile, *refClimM1OutFile, *refM1OutFile;
  char *refClimM3OutFile, *refM3OutFile;

  /* 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 \n", &totalVerModes);

  climWindFile = calloc(40, sizeof(char));
  windFile = calloc(40, sizeof(char));
  climOutFile = calloc(40, sizeof(char));
  outFile = calloc(40, sizeof(char));
  refClimM1OutFile = calloc(40, sizeof(char));
  refM1OutFile = calloc(40, sizeof(char));
  refClimM3OutFile = calloc(40, sizeof(char));
  refM3OutFile = 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", refClimM1OutFile);
      fscanf(input, "%s \n", refM1OutFile);
      fscanf(input, "%s \n", refClimM3OutFile);
      fscanf(input, "%s \n", refM3OutFile);

      /*
	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(!(ofprcM1 = fopen(refClimM1OutFile, "w")))
	{
	  printf("open file error: %s!! \n", refClimM1OutFile);
	  exit(1);
	}
      if(!(ofprM1 = fopen(refM1OutFile, "w")))
	{
	  printf("open file error: %s!! \n", refM1OutFile);
	  exit(1);
	}

      if(!(ofprcM3 = fopen(refClimM3OutFile, "w")))
	{
	  printf("open file error: %s!! \n", refClimM3OutFile);
	  exit(1);
	}
      if(!(ofprM3 = fopen(refM3OutFile, "w")))
	{
	  printf("open file error: %s!! \n", refM3OutFile);
	  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, ii, tauxc, eastbc);      
      integral(ofp, ii, tauxc, taux, eastb);

      /* calculate reflection of 1st and 3rd meri models */
      reflectClim(ofprcM1, 1, ii, 1., eastbc);
      reflect(ofprM1, 1, ii, 1., eastb, eastbc);
      reflectClim(ofprcM3, 3, ii, 1.225, eastbc);
      reflect(ofprM3, 3, ii, 1.225, eastb, eastbc);

      /* close all the file for a mode */
      fclose(ifpc);
      fclose(ifp);
      fclose(ofpc);
      fclose(ofp);
      fclose(ofprcM1);
      fclose(ofprM1);
      fclose(ofprcM3);
      fclose(ofprM3);
    }    
}