// based on ascii2cdf_multi
// specific to reading in files of 
// OPRA data 
// specific for files created during MAP
// SEY 28 August 2000
//
// SEY 12 Dec 2000 added V field with thresholding on NCP
//
// SEY 18 Jan 2001 added first guess calc of w using methods of Joss and MH
//
//
// puts in the tr tube correction
// there are 100 gates in the MAP data
//
// 18 Dec 1998 according to Grant, the first 3 gates are always bad.
// first gate is in the radar
// second gate is 0-150 m up
//
// testing
// cd /home/disk/hail/yuter/mtn/rib
// /home/disk/rain2/dt/yuter/src/aiw-src/opra2cdf 1461 980516.1340.1440.trim
// /home/disk/rain2/dt/yuter/src/aiw-src/opra2cdf 20 980516.1340.1440.trim.20
//
// SEY 1 Dec 1998 revision to put parameters on command line
// /home/disk/rain2/dt/yuter/src/aiw-src/opra2cdf 100 0.5 0.0 47.65 -122.3 /home/disk/corona/kpearl/opradata/opr81121.13z 0 0 7
//
///home/disk/rain2/dt/yuter/src/aiw-src/opra2cdf 100 0.5 0.0 47.65 -122.3 /home/disk/corona/kpearl/opradata/opr81125.12z -25 0 7
// /home/echo/yuter/src/cdf-src/mapopra2cdf 

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sunmath.h>
#include <stdlib.h>                       /* standard library functions      */
//#include "volume.h"  /* general header for volume   */
//#include "mudras.h" 
#include "defs.h" 
#include "zebcdflib.h" 

#define GATERESKM 0.15 // 150 meters
#define MAXTIMES 60*60*2 // 60 min * 60 sec * 2
#define MISSING -99.0
#define MAXLINE 1050

float computeMeanDropVelocity (float Z, float temp, float pres, 
			  float No, float m);
float estimateTemp(float surftemp, float gatenum, float lowestlevelkm);
float estimatePressure(float surpres, float gatenum, float lowestlevelkm);
float estimateMHVtrain (float Z, float height);
float estimateMHVtsnow (float Z, float height);
float computeMassflux (float w, float height);
int main(int argc, char **argv)

{
  if (argc !=7)  {
     fprintf(stdout, 
  "usage: %s Zfilename Vfilename Nfilename firstgate topheighttoprocesskm deriveflag (0=no 1=yes)\n",
           argv[0]);
     exit(1);
  }
  int q;
  static char date[10];
  static char time[10];
  static char holdstring[8];
  char syear[3];
  char smonth[3];
  char sday[3];
  char shour[3];
  char sminute[3];
  char ssecond[3];
  char fieldname[20];
  static char rest_of_line[MAXLINE+1]; // this is like a loose cannon variable
  int year, month, day, hour, minute, second; // start time
  int lyear, lmonth, lday, lhour, lminute, lsecond; // time for a given line
  int numgrids = 6;
  int lineindex= 0;
  float trcorr[]={49.98, 49.98, 69.50, 61.40, 57.50, 55.30, 53.80, 52.90,52.20,
		   51.70, 51.30,51.05, 50.80, 50.66, 50.50, 50.38, 50.26,
		   50.20, 50.12, 50.03, 50.00};
  int numzap=0;
     float pressure, temperature, Z, Vr;
     int l, i;

    // the correction at each gate is trcorr[i]-trcorr[0]
  /* char *fn   = argv[]; */
  //strcpy(date,argv[1]);
  //strcpy(time,argv[2]);
  
  // fill in the stuff that is fixed for MAP
  int numgates = 100; // not fixed for OPRA MAP is 84 in some circumstances
  float timeressec = 1; // fixed for OPRA MAP
  float radarheight =0.307; // fixed for OPRA MAP height of radar AGL
  
  float origin_lat=46.17;
  float origin_lon=8.75;
  
  char **gridfnlist = new char*[numgrids];
  gridfnlist[0]=argv[1];
  gridfnlist[1]=argv[2];
  char ncpfn[MAX_LINE];
  strcpy(ncpfn,argv[3]);
  float Zoffset = 0; // no offset
  int firstgate = atoi(argv[4]);
  float topheight = atof(argv[5]);
  int deriveflag = atoi(argv[6]);

  if (deriveflag==0) numgrids=2; // just do maxdz and radial velocity

  float lowestlevelkm = 0.307; // 0.382; // height of radar AGL
  // first gate is centered at this height 382-75=307
  // since the first gate is in the radar hack adjust lowest level
  // adjust zeb display accordingly
  lowestlevelkm=lowestlevelkm-GATERESKM;

  // figure out top gate to process
  int numgatesproc = irint((topheight-lowestlevelkm)/GATERESKM);
  fprintf(stderr,"Number of gates writing to cdf file is %d\n",
	  numgatesproc);
  
  printf("number of grids %d\n",numgrids);
  // create the arrays of pointers 
  char **gridnames = new char*[numgrids]; 
  float **gridlist = new float*[numgrids];
  // read in the fieldnames and their assoc files
  gridnames[0]="maxdz";
  gridnames[1]="radialvelocity";
  // computing two versions of w for debugging
  if (deriveflag){
  gridnames[2]="w";
  gridnames[3]="fallspeed";
  gridnames[4]="w2";
  gridnames[5]="massflux";
  }
  
  fprintf(stderr,"Zfilename is %s\n",gridfnlist[0]);
  fprintf(stderr,"Vfilename is %s\n",gridfnlist[1]);
  fprintf(stderr,"Nfilename is %s\n",ncpfn);
   
  // this will be the prefix for the cdf file
  char *p;
  p="opra99";
  strcpy(fieldname,p);
 
  float clat = origin_lat; 
  float clon = origin_lon;

  FILE *fp;
  FILE *ncpfp;
  float* grid;

  // first read in the ncp data, that is used for qc on the V field
  float* ncpgrid = new float[numgatesproc*MAXTIMES];
     float ncp;
     ncpfp = fopen(ncpfn, "r");
     int j=0;
     while (fscanf(ncpfp, "%f", &ncp)!=EOF) {
	if (j==(numgatesproc+2)) {
	 j=-1; // rest x index and incre lineindex
	 //fprintf(stderr,"about to set lineindex from %d\n",lineindex);
       lineindex++;
       fgets(rest_of_line, MAXLINE, ncpfp);
       }

	if (j>=2) {
	 int dataarrayindex=numgatesproc*lineindex+j-2;
	 ncpgrid[dataarrayindex] = ncp;
	}
	j++;
     }
     
     // debug
     for (int j=0;j<numgatesproc*2;j++){
        fprintf(stderr,"%.1f ",ncpgrid[j]);
      if (j==(numgatesproc-1)) fprintf(stderr,"\n");
     }
       
     fprintf(stderr,"\n ");
  // for each input grid read in the data 
 float* offsettimes = new float[MAXTIMES];
 float* gridin1 = new float[numgatesproc*MAXTIMES];
 float* gridin2 = new float[numgatesproc*MAXTIMES];
 gridlist[0]=gridin1;
 gridlist[1]=gridin2;

 for (q=0;q<2;q++){
     /* FILE *fp = fopen(fn, "r");*/
   //fprintf(stderr,"q=%d\n",q);
     fp = fopen(gridfnlist[q], "r");
     fprintf(stderr,"reading in %s\n",gridfnlist[q]);
     if (!fp)
       {
	   fprintf(stderr, "Can't open file %s.\n",gridfnlist[q]);
	   exit(1);
       }
     float v;
     lineindex=0;
     
     // start reading in the file
     while (fscanf(fp, "%s", &holdstring)!=EOF) {
       v=atof(holdstring);
       if (j==(numgatesproc+2)) { // end of line to read
	 j=-1; // rest x index and incre lineindex
	 //fprintf(stderr,"about to set lineindex from %d\n",lineindex);
       lineindex++;
       fgets(rest_of_line, MAXLINE, fp);
       }
	if (lineindex > MAXTIMES) {
	 fprintf(stderr,
		 "number of lines in file %d exceeds MAXTIMES can handle.. Aborting\n",lineindex);
	 exit(1);
       }
	if (j==0) { // date
	  strcpy(date,holdstring);
	  if (q==0){ // save the date
	  sprintf(syear,&date[0]);
	  syear[2] = '\0';
	  lyear   = atoi(syear);

	  sprintf(smonth,&date[2]);
	  smonth[2] = '\0';
	  lmonth  = atoi(smonth);

	  sprintf(sday,&date[4]);
	  sday[2] = '\0';
	  lday    = atoi(sday);
	  
	  if (lineindex==0) { // this first time is the basetime
	    year=lyear;
	    month=lmonth;
	    day=lday;
	  }}}
	if (j==1) { // time
	  strcpy(time,holdstring);
	  if (q==0){ // save the time
	  sprintf(shour,&time[0]);
	  shour[2] = '\0';
	  lhour   = atoi(shour);

	  sprintf(sminute,&time[2]);
	  sminute[2] = '\0';
	  lminute = atoi(sminute);

	  sprintf(ssecond,&time[4]);
	  ssecond[2] = '\0';
	  lsecond = atoi(ssecond);
	
	  if (lineindex==0){ // this first time is the basetime
	    hour=lhour;
	    minute=lminute;
	    second=lsecond;

       fprintf(stderr, "mapopra2cdf: Date = %d %d %d and time = %d %d %d\n",
	       year,month,day,hour, minute,second);
	  }}}

	  
	if ((j==1) && (lineindex < 10)
	    && (q==0)) {// just print out the time for the
	   // first few so that it does not slow down the processing. 
	   fprintf(stderr,"time is %s\n",time);
	   fprintf(stderr,"offsettime is %d\n",
		   computeTimeOffsetsecs(year,month,day,hour,minute,second,
				lyear,lmonth,lday,lhour,lminute,lsecond));
	  }
	  if ((j==1) && (q==0)) { // compute the offset seconds
	    offsettimes[lineindex]=(float) computeTimeOffsetsecs(
				year,month,day,hour,minute,second,
				lyear,lmonth,lday,lhour,lminute,lsecond);
	  }
     
	// deal with the data
	int dataarrayindex=numgatesproc*lineindex+j-2;
	if (j>=2) { // is a gate
	  if (v == MISSING) v=v-Zoffset;
	  if ((j-2)<firstgate) gridlist[q][dataarrayindex] = MISSING;
	  else if (((j-2)<20) && (q==0)){ // apply tr correction to Z
	    gridlist[q][dataarrayindex] = v+trcorr[j-2]-trcorr[0]+Zoffset;
	    /* fprintf(stderr,"correction applied %.2f yielding %.2f\n",
	       trcorr[j-2]-trcorr[0],grid[j-2]); */
	  }
	  else gridlist[q][dataarrayindex] = v+Zoffset;
	
	if ((q==1) && (ncpgrid[dataarrayindex]<0.15)){ // Vr data thresholding
	  numzap++;
	  //fprintf(stderr,"zapping V value of %.2f\n",grid[dataarrayindex]);
	  gridlist[q][dataarrayindex] = MISSING;
	}}
	j++;
     }

     fprintf(stderr,"lineindex processed = %d, q= %d\n",lineindex+1,q);
     fclose(fp);
 } // end for q loop  // end loop for eAch file

 
 // for derived fields
 // output grids
 float* gridout1 = new float[numgatesproc*MAXTIMES];
 float* gridout2 = new float[numgatesproc*MAXTIMES];
 float* gridout3 = new float[numgatesproc*MAXTIMES];
 float* gridout4 = new float[numgatesproc*MAXTIMES];

if (deriveflag) {

fprintf(stderr,"now derive the fall and air velocities for %d gates\n",numgatesproc);

 // debug
 fprintf(stderr,"DEBUG meandropvelocity = %.4f\n",
	 computeMeanDropVelocity(26.08,0,1013,2179,2.38));

 for (l=0;l<(lineindex+1);l++){
   for (i=0;i<numgatesproc;i++){
     temperature=estimateTemp(21,i,lowestlevelkm);
     pressure=estimatePressure(1013,i,lowestlevelkm);
     float height=lowestlevelkm+(i*GATERESKM);
     if (l<1){
     fprintf(stderr,"gate=%d,height=%.2f,temp=%.2f,pres=%.2f\n",i,height,temperature,pressure);
     }
     Z= gridlist[0][l*numgatesproc+i];
     Vr= gridlist[1][l*numgatesproc+i];
     //Z=26.08;
     //Vr=0;
     if ((Vr==MISSING) || (Z==MISSING)){ 
       gridout1[l*numgatesproc+i]=MISSING;
       gridout2[l*numgatesproc+i]=MISSING;
       gridout3[l*numgatesproc+i]=MISSING;
       gridout4[l*numgatesproc+i]=MISSING;
     }
     else {
       gridout1[l*numgatesproc+i]=Vr+computeMeanDropVelocity(Z,temperature,pressure,8000,0);
     
       //gridout2[l*numgatesproc+i]=Vr+computeMeanDropVelocity(Z,temperature,pressure,2179,2.38);

       // zero Doppler velocity, so just the fall speed
       gridout2[l*numgatesproc+i]=0+computeMeanDropVelocity(Z,temperature,pressure,8000,0);

       // testing MH methods
       if (temperature > 0) // in rain
	 gridout3[l*numgatesproc+i]=Vr+estimateMHVtrain(Z,height);
       // in snow
       else gridout3[l*numgatesproc+i]=Vr+estimateMHVtsnow(Z,height);

       // mass flux
       gridout4[l*numgatesproc+i]=
	 computeMassflux(gridout1[l*numgatesproc+i],height);
     }

   }}

 gridlist[2]=gridout1;
 gridlist[3]=gridout2;
 gridlist[4]=gridout3;
 gridlist[5]=gridout4;
 }

 fprintf(stderr,"last time processed %s %s\n",date,time);
 fprintf(stderr,"numgatesproc = %d, lines processed = %d\n",
	 numgatesproc,lineindex+1);
 // since lineindex counts from 0 total number of lines is
 // lineindex +1

writeCDFxtime_multi(numgatesproc,lineindex,GATERESKM,timeressec,year,month,day,
		   hour, minute, second,
		   numgrids,offsettimes,gridnames, gridlist,
		   1.0,
		    MISSING, // missing val
		   clat, clon,
		   lowestlevelkm, // altitude of the grid 
		    filename(".", fieldname,
                    year, month, day, hour, minute));


  fprintf(stderr,"Number of V values below NCP threshold =%d\n",numzap); 

  delete[] gridout1;
  gridout1=NULL;
  delete[] gridout2;
  gridout2=NULL;
  delete[] gridout3;
  gridout3=NULL;
  delete[] gridout4;
  gridout4=NULL;
  delete[] gridin1;
  gridin1=NULL;
  delete[] gridin2;
  gridin2=NULL;
  delete[] gridnames;
  gridnames=NULL;
  delete[] gridfnlist;
  gridfnlist=NULL;
  delete[] gridlist;
  gridlist=NULL;
  delete[] ncpgrid;
  ncpgrid=NULL;
  return(0);
}

// OPRA QC
// SEY 1 Dec 1998
// if adjacent gates are different by more than x dB then set
// the gate at higher altitude to 0

// compute mean drop velocity in m/s using Joss' algorithm
// Z in dBZ, temp in C, pressure in hPa
float computeMeanDropVelocity (float Z, float temp, float pres, 
			  float No, float m) {

float a, k, b, Vt;

 if (Z==MISSING) Vt=MISSING;
 else{
   a=pow(((temp+273)/pres),0.45)*pow((No),-0.12);
   k=0.4099*a+0.00558/a;
   b=-0.15+15.48*a-1.3229*(log10(m+1)) + 0.03975*m;
   Vt=(b+k*Z); // fall speed is defined as positive
 }

 return(Vt); 
}

// use standard lapse rate to calc temp in deg C from surface temp
float estimateTemp(float surftemp, float gatenum, float lowestlevelkm){

float height, esttemp;

height=lowestlevelkm+(gatenum*GATERESKM);
// standard lapse rate is 7 deg /km
// MAP lapse rate of 6 deg/km, modified 3 May 2001 SEY

esttemp=surftemp-height*6.0;

return(esttemp);
}

float estimatePressure(float surpres, float gatenum, float lowestlevelkm){

float height, estpress;
height=lowestlevelkm+(gatenum*GATERESKM);
// WH page 56
// p2=p1exp((z1-z2)/H)
// where H is the scale height Tv*29.3
// assuming Tv=288, then H=8.4384 km
// for isothermal, dry atmosphere
estpress=surpres*exp(-height/8.4384);

return(estpress);
  }

// Marks and Houze 1987 from Joss and Waldvogel 1970
float estimateMHVtrain (float Z, float height) {

float rhocorr, Vt;

 if (Z==MISSING) Vt=MISSING;
 else{
   rhocorr=pow((1/(exp(height/10.0))),0.45);
   Vt=rhocorr*2.6*pow(Z,0.107);
 }

 return(Vt); 
}

float estimateMHVtsnow (float Z, float height) {

float rhocorr, Vt;

 if (Z==MISSING) Vt=MISSING;
 else{
   rhocorr=pow((1/(exp(height/10.0))),0.45);
   Vt=rhocorr*0.817*pow(Z,0.063);
 }

 return(Vt); 
}

// compute mass flux in kg s-1 m-2 i.e. dont normalize by beamwidth
float computeMassflux (float w, float height){

float density;
float massflux;

density=1.275*exp(-height/10.0);

// fprintf(stderr,"height = %.2f density =%.2f\n",height, density);

massflux=w*density;

return(massflux);

}