/* C Program example to compute mean and std. deviation for */
/* precipitation data in homework #3 */

/* ASCII file version */


/* External libraries */
/* These are need by the program to do things like IO */
/* and math functions                                 */


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

/* Define array limits */
/* For the HW #3 data these correspond to the data dimensions */
/* in space and time */

#define XPOINTS 81
#define YPOINTS 51
#define MONTHS 12
#define YEARS 55

/* Main program */
/* Note: you would put any function declarations to be called
   within the main program before this */

main(){

  /* Variable declarations: each variable must be specified as a certain */
  /* data type (i.e. int, float, double, etc.) */

  /* Counting dimensions for space and time (integers) */

  int i,j,p,q;

  /* Other variable declarations */

  int xmaxpoints, ymaxpoints;
  float missingval;


  /* Precipitation declared as three dimensional array */

  float precip[XPOINTS][YPOINTS][MONTHS][YEARS];
  float missingarray[XPOINTS][YPOINTS][MONTHS];
  float sumprecip[XPOINTS][YPOINTS][MONTHS];
  float sumprecipsq[XPOINTS][YPOINTS][MONTHS];
  float meanprecip[XPOINTS][YPOINTS][MONTHS];
  float stdevprecip[XPOINTS][YPOINTS][MONTHS];
  float years;
 

  /* Input and output files */

  FILE *inPtr, *outPtr1, *outPtr2;


  /* Open input and output files */
  /* Format: fopen statement, file name, read/write */

  inPtr=fopen("precipmonthly_usmx.asci", "r");


  outPtr1=fopen("Mon_means.asci", "w");
  outPtr2=fopen("Mon_std.asci", "w");


  /* Declare missing value */

  missingval=-999.9;

  years=55;


  /* Initialize missing array */
  
  for(p=0;p<MONTHS;p++){
    for(j=0;j<YPOINTS;j++){
      for(i=0;i<XPOINTS;i++){
        missingarray[i][j][p]=0;
      }
    }
  }


  /* Read in data.  Keep track of missing values in the process */

  for(q=0;q<YEARS;q++){
    printf("q=%d\n", q);
    for(p=0;p<MONTHS;p++){
      for(j=0;j<YPOINTS;j++){
        for(i=0;i<XPOINTS;i++){
          printf("i=%d\n", i);
          fscanf(inPtr, "%f", &precip[i][j][p][q]);
          if(precip[i][j][p][q]==missingval){
            missingarray[i][j][p]=1;
	  }
        }
      }
    }
  }

  /* Initialize summation arrays */

  for(p=0;p<MONTHS;p++){
    for(j=0;j<YPOINTS;j++){
      for(i=0;i<XPOINTS;i++){
        sumprecip[i][j][p]=0;
        sumprecipsq[i][j][p]=0;
      }
    }
  }


  /* Compute mean and std. deviations, ignorning missing values */

  for(p=0;p<MONTHS;p++){
    for(q=0;q<YEARS;q++){
      for(j=0;j<YPOINTS;j++){
        for(i=0;i<XPOINTS;i++){
          if(missingarray[i][j][p]==0){
            sumprecip[i][j][p]=sumprecip[i][j][p]+precip[i][j][p][q];
            sumprecipsq[i][j][p]=sumprecipsq[i][j][p]+precip[i][j][p][q]*precip[i][j][p][q];
	  }	  
	}
      }
    }
  }

  for(p=0;p<MONTHS;p++){
    printf("p=%d\n", p);
    for(j=0;j<YPOINTS;j++){
      for(i=0;i<XPOINTS;i++){
        if(missingarray[i][j][p]==0){
          meanprecip[i][j][p]=sumprecip[i][j][p]/years;
          stdevprecip[i][j][p]=sqrt((years*sumprecipsq[i][j][p]-(sumprecip[i][j][p]*sumprecip[i][j][p]))/(years*(years-1)));
	}
        if(missingarray[i][j][p]==1){
          meanprecip[i][j][p]=missingval;
          stdevprecip[i][j][p]=missingval;
	}
      }
    }
  }

  /* Write out mean and std. deviation to an asci file */

  for(p=0;p<MONTHS;p++){
    for(j=0;j<YPOINTS;j++){
      for(i=0;i<XPOINTS;i++){
	fprintf(outPtr1, "%.4f\n", meanprecip[i][j][p]);
        fprintf(outPtr2, "%.4f\n", stdevprecip[i][j][p]);
      }
    }
  }        

return 0;

}