estimate_daily_tstations.c

Go to the documentation of this file.

 
#include "prototypes_new.h"
 
void estimate_daily_tstations(int j)
 
{
 
extern int isom;
extern struct tdata tdata[10];
extern struct station tstation[1000];
extern int max_tstations;
int m,k,i,l,ii;
double lat1,lon1,fdist,fdata,fvalue[4],lat,lon,testdist,fmult,ftotal,isoh,isoh1,padj,stotal;
int num_missing;
double distlon,dist,df;
float conv=.0174;
int h;
 
if(tdata[j].data_time==0)
                     return;
 
for(m=0;m<max_tstations;m++) {
 
          /* dont estimate missing 24 hour tstations */
      
          for(h=0;h<4;h++)    
           tdata[j].stn[m].tlevel2[h].a=-99;
          
          for(h=0;h<6;h++) 
              if(tdata[j].stn[m].tlevel2[h].data == -99 ||
            (tdata[j].stn[m].tlevel2[h].qual!=0 &&
         tdata[j].stn[m].tlevel2[h].qual!=8 &&
                 tdata[j].stn[m].tlevel2[h].qual!=3 &&
         tdata[j].stn[m].tlevel2[h].qual!=2))
                                    break;
                    
      if(h!=6)
             continue;      
                   
      for(h=0;h<4;h++) {
      
           tdata[j].stn[m].tlevel2[h].a=
           (float)(tdata[j].stn[m].tlevel2[h].data-tdata[j].stn[m].tlevel2[5].data)/
           (float)(tdata[j].stn[m].tlevel2[4].data-tdata[j].stn[m].tlevel2[5].data);
                       
 
           }
           
           }
                         
for(m=0;m<max_tstations;m++) {
                  
          for(k=0;k<4;k++)         
               if(tdata[j].stn[m].tlevel2[k].data != -99 &&
              tdata[j].stn[m].tlevel2[k].qual != 6)
                                 break;
                     
          if(k != 4)
               continue;
           
      for(k=0;k<4;k++)         
                         tdata[j].stn[m].tlevel2[k].data = -99; 
                
      if(tdata[j].stn[m].tlevel2[4].data == -99 ||
             tdata[j].stn[m].tlevel2[5].data == -99)
                                          continue; 
                        
          if(tdata[j].stn[m].tlevel2[4].qual==1 ||            
             tdata[j].stn[m].tlevel2[4].qual==5 ||
         tdata[j].stn[m].tlevel2[5].qual==1 ||            
             tdata[j].stn[m].tlevel2[5].qual==5) 
                                           continue;  
     
          lat1=tstation[m].lat;
          lon1=tstation[m].lon;
      
          /*get isohyet */
 
      for(k=0;k<4;k++) {
 
             fdist=0.0;
             fdata=0.0;
 
             l=0;
 
             for(ii=0;ii<20;ii++){
 
                 i=tstation[m].index[ii]; 
 
                 /* dont estimate unless good or forced good */
 
                 if(tdata[j].stn[i].tlevel2[k].qual!=0 &&
                    tdata[j].stn[i].tlevel2[k].qual!=8 &&
                    tdata[j].stn[i].tlevel2[k].qual!=3 &&
            tdata[j].stn[i].tlevel2[k].qual!=2)
                                          continue;
         
                 /*dont use missing tstations */
                 
                 if(tdata[j].stn[i].tlevel2[k].data ==-99 ||
            tdata[j].stn[i].tlevel2[k].a < -98)
                                          continue;
 
                 lat=tstation[i].lat;
                 lon=tstation[i].lon;
         
         distlon=(lon1-lon)*cos((lat1+lat)/2*conv);
 
                 dist= pow((double)(lat1-lat),2) +
                       pow((double)(distlon),2);
               
             dist=pow(dist,.5)*60;
            
         df=50*abs(tstation[m].elev-tstation[i].elev)/5280; 
                           
                 dist=dist+df;
         
         if(dist==0.0)
                             dist=.000001;
         dist=1/dist;
 
         fdata=fdata+tdata[j].stn[i].tlevel2[k].a*dist;
         fdist=fdist+dist;
           
                 l++;
         
         if(l==5)
                break;
 
           }
     
 
             if(l < 5) {
 
             fdist=0.0;
             fdata=0.0;
 
             l=0;
 
             for(i=0;i<max_tstations;i++){
                  
                 if(i==m)
                      continue;
 
                 if(tdata[j].stn[i].tlevel2[k].qual!=0 &&
                    tdata[j].stn[i].tlevel2[k].qual!=8 &&
                    tdata[j].stn[i].tlevel2[k].qual!=3 &&
            tdata[j].stn[i].tlevel2[k].qual!=2)
                                          continue;
      
                 if(tdata[j].stn[i].tlevel2[k].data == -99 ||
            tdata[j].stn[i].tlevel2[k].a < -98)
                                          continue;
 
                 lat=tstation[i].lat;
                 lon=tstation[i].lon;
         
         distlon=(lon1-lon)*cos((lat1+lat)/2*conv);
 
                 dist= pow((double)(lat1-lat),2) +
                       pow((double)(distlon),2);
               
             dist=pow(dist,.5)*60;
            
         df=50*abs(tstation[m].elev-tstation[i].elev)/5280; 
                           
                 dist=dist+df;
         
                 if(dist==0.0)
                             dist=.000001;
 
         dist=1/dist;
 
         fdata=fdata+tdata[j].stn[i].tlevel2[k].a*dist;
         fdist=fdist+dist;
           
                 l++;
          
         }
           
         }
 
           if(l != 0){
         
           tdata[j].stn[m].tlevel2[k].a=fdata/fdist;
               tdata[j].stn[m].tlevel2[k].data=tdata[j].stn[m].tlevel2[k].a*
           (tdata[j].stn[m].tlevel2[4].data-tdata[j].stn[m].tlevel2[5].data) +
           tdata[j].stn[m].tlevel2[5].data;
               tdata[j].stn[m].tlevel2[k].qual=6;
           
           }
           
           
        
       }
 
            
}
}