get_frzlvl_data.c

Go to the documentation of this file.

#include "prototypes_new.h"
 
int maxi,maxj;
float latitude[50][50],longitude[50][50],frzlvl[50][50],top[50][50];
void get_elev(float,float,float *, float *);
 
 
void get_frzlvl_data(int ii)
   
{
 
extern struct pdata pdata[10];
extern struct station station[3000];
extern int max_stations;
extern char freezing_level_directory[100];
FILE *fp;
char ibuf[150],*sname[3],*p;
float tval[50][50],val[50];
int j,k,m,iflag,ier;
char hb5[10],parm[10];
int tip,elev;
float lat,lon;
float conv;
char fbuf[50],savebuf[50];
time_t tget;
struct tm *gmtim;
int i;
float temp,felev,diff,tempo;
float lapse=3.3;
 
for(k=0;k<max_stations;k++) {
 
                  for(m=0;m<5;m++) 
                               pdata[ii].stn[k].frzlvl[m]=-99;
 
        }
 
conv=360/(2*3.14159);
 
sname[0]="latitude";
sname[1]="longitude";
sname[2]="topo";
 
for(m=0;m<3;m++) {
 
    strcpy(fbuf,freezing_level_directory);
    strcat(fbuf,"/");
    strcat(fbuf,sname[m]);
    fp=fopen(fbuf,"r");
 
    if(fp==NULL) {
 
              printf("could not open %s\n",fbuf);
              return;
 
        }
 
    j=-1;
    iflag=0;
    maxi=0;
    maxj=0;
 
    for(;;) {
 
            p=fgets(ibuf,100,fp);
 
            if(p==NULL)
                      break;
 
            if(strstr(ibuf,"ROW ")!=NULL) {
 
                   iflag=1;
                   p=p+8;
                   k=0;
                   j++;
 
         }
 
            if(iflag==1) {
 
 
                   ier=sscanf(p,"%f %f %f %f %f %f %f %f",
                              &val[0],&val[1],&val[2],&val[3],
                              &val[4],&val[5],&val[6],&val[7]);
 
                   for(i=0;i<ier;i++)  {
                             
                              tval[j][i+k]=val[i];
 
                              if(m==0)
                                 tval[j][i+k]=tval[j][i+k]*conv/100;
 
                              if(m==1)
                                 tval[j][i+k]=-tval[j][i+k]*conv/10;
 
                              if(m==2)
                                 tval[j][i+k]=tval[j][i+k];
 
                }
 
                   k=k+ier;
 
                   if(k > maxi)
                         maxi=k;
 
         }
 
      }
 
  maxj=j;
 
  if(m==0)
        memcpy(latitude,tval,2500*sizeof(float));
 
  if(m==1)
        memcpy(longitude,tval,2500*sizeof(float));
 
  if(m==2)
        memcpy(top,tval,2500*sizeof(float));
 
  fclose(fp);
 
  }
 
for(m=0;m<4;m++) {
 
/* get correct date */
 
                 tget=pdata[ii].data_time;
 
                 if(m < 2) 
                      tget=tget-86400;
          
                 gmtim=gmtime(&tget);
 
                 sprintf(fbuf,"%s/sfctemp.%02d-%02d-%02d.",
                         freezing_level_directory,
                         gmtim->tm_mon+1,gmtim->tm_mday,gmtim->tm_year);
 
                 if(m==0)
                       strcat(fbuf,"15");
           
                 else if(m==1)
                       strcat(fbuf,"21");
 
                 else if(m==2)
                       strcat(fbuf,"03");
           
                 else if(m==3)
                       strcat(fbuf,"09");
 
                 fp=fopen(fbuf,"r");
 
                 if(fp==NULL) {
 
                      printf("%s not found\n",fbuf);
 
                      continue;
 
                    }
 
 
                 j=-1;
                 iflag=0;
                
                 for(;;) {
 
                   p=fgets(ibuf,100,fp);
 
                   if(p==NULL)
                      break;
 
                   if(strstr(ibuf,"ROW ")!=NULL) {
 
                      iflag=1;
                      p=p+8;
                      k=0;
                      j++;
 
         }
 
                  if(iflag==1) {
 
 
                     ier=sscanf(p,"%f %f %f %f %f %f %f %f",
                                &val[0],&val[1],&val[2],&val[3],
                                &val[4],&val[5],&val[6],&val[7]);
 
                     for(i=0;i<ier;i++)  {
 
                              tval[j][i+k]=val[i];
 
                }
 
                   k=k+ier;
 
         }
 
      }
           
        memcpy(frzlvl,tval,2500*sizeof(float));
 
        for(k=0;k<max_stations;k++) {
 
            if(station[k].tip==1)
                   continue;
 
            get_elev(station[k].lat,station[k].lon,&felev,&temp);
 
            diff=felev - station[k].elev;
 
            tempo=temp+diff*lapse/1000;
            tempo= tempo* 9./5. + 32.;
            pdata[ii].stn[k].frzlvl[m]=(int)tempo;
 
      }
 
    fclose(fp);   
 
    }
 
}
 
void get_elev(float lat, float lon,float *felev, float *temp)
 
{
 
int i,j;
double pgon[4][2],point[2];
double distance,value;
int ii,jj;
double dist, dist1,dist2;
int flag,elev;
double value2,felev2;
 
point[0]=lat;
point[1]=lon;
 
for(j=0;j<maxj;j++) {
for(i=0;i<maxi;i++) {
 
  pgon[0][0]=latitude[j][i];
  pgon[0][1]=longitude[j][i];
  pgon[1][0]=latitude[j][i+1];
  pgon[1][1]=longitude[j][i+1];
  pgon[2][0]=latitude[j+1][i+1];
  pgon[2][1]=longitude[j+1][i+1];
  pgon[3][0]=latitude[j+1][i];
  pgon[3][1]=longitude[j+1][i];
 
  flag=macmartintest(pgon,4,point);
 
  if(flag==1) 
           goto cont;
 
 
}
 
}
 
printf("not found\n");
return;
 
cont:
 
distance=0.0;
value=0.0;
value2=0.0;
 
for(jj=j;jj<=j+1;jj++) {
      
for(ii=i;ii<=i+1;ii++) {
 
           dist1=point[0]-latitude[jj][ii];
           dist2=point[1]-longitude[jj][ii];
 
           dist=pow(dist1,2)+pow(dist2,2);
 
           if(dist==0.0)
                      dist=1.0;
 
           dist=1/dist;
 
           value=value+frzlvl[jj][ii]*dist;
           value2=value2+top[jj][ii]*dist;
 
           distance=distance+dist;
 
             }
             }
    
*temp=value/distance;      
*felev=value2/distance * 3.048;
 
return;
 
}