water386_407a.cpp

处理激光雷达数据

/***************************************************************************/
/* PROG for calculating water vapor mixing ratio                           */
/*  1. N2 and Water Raman  signals                                         */
/*  2. Input *.407 and *.386 MCS data                                      */
/*                                                                         */
/* Edited on 2004-03-31                                                    */
/*                                                                         */
/***************************************************************************/

#include "READMCS_1.C"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <process.h>
#include <string.h>



void main()
{
  FILE *fpout,*fmol,*fact,*faer, *fr;
   int i,j,k,p;

//  char fname[20];

 float ratio1[1000],ratio2[1000];

 for(i=0;i<1000;i++)
  {
    ratio1[i]=0;ratio2[i]=0;
   }


/********************************************************************/
/*  Read MCS file  "*.407","*.386"
/*  Read channel geometrical correction MCS file, 407channel and
/*       386-channel.  calculate factor ,
/*   And make correction
/********************************************************************/

 float p407[1000],p386[1000],pnb[1000],pna[1000], temp[1000];
  char str1[20],str2[20],str3[20],filename[20];

 for(i=0;i<1000;i++)
  {
   p407[i]=0; p386[i]=0;pnb[i]=1;pna[i]=1; temp[i]=0; 
   }

 printf("Input date,such as: 040315\n");
 scanf("%s",str1);

 strcpy(str3,"");
 strcpy(str3,str1);
 strcat(str3,"sig407_386.dat");
 if( (fpout=fopen(str3,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str3);
     exit(1);
    }

 printf("\n\nInput MCS file (w?.407) group number,such as: 2 or 3 or 4....\n");
 scanf("%d",&k);


 for(i=1;i<k+1;i++)
  {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"w00");
     strcat(filename,str2);
     strcat(filename,".407");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"%s\n",filename);
     for(j=0;j<1000;j++)
      p407[j]=p407[j]+temp[j];
   }

  /* for(j=0;j<1000;j++)
     printf("%d  %.3f\n",j,p407[j]);   */

  printf("\n\n Input MCS file (n?.386) number,such as: 2 or 3 or 4....\n");
  scanf("%d",&k);


  for(i=1;i<k+1;i++)
  {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"n00");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"\n%s",filename);
     for(j=0;j<1000;j++)
        p386[j]=p386[j]+temp[j];
   }

 /*for(j=0;j<1000;j++)
   fprintf(fpout,"\n%d  %.3f",j,p386[j]);         */

 printf("\n\n Channel geometrical factor correction?\n");
 printf("     Yes, choise  '1' \n");
 printf("      No, choise  '0' \n");
 scanf("%d",&p);

 if(p==1)
  {
   for(i=0;i<1000;i++)
    {
      pnb[i]=0;
      pna[i]=0;
     }
   printf("\n Input channel geometric MCS file (a?.386) number\n");
   printf("       such as: 1 or 2 or 3... \n");
   scanf("%d",&k);

    fprintf(fpout,"\n");
   for(i=1;i<k+1;i++)
    {
     strcpy(filename,"");;
     itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"wa100");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"\n%s",filename);
     for(j=0;j<1000;j++)
        pna[j]=pna[j]+temp[j];
     }
   /*for(j=0;j<1000;j++)
       fprintf(fpout,"%d  %.3f\n",j,pna[j]);    */

   printf("\n Input channel geometric MCS file (b?.386) number\n");
   printf("       such as: 1 or 2 or 3... \n");
   scanf("%d",&k);

   fprintf(fpout,"\n");
   for(i=1;i<k+1;i++)
    {
     strcpy(filename,"");;
      itoa(i,str2,10);
     strcat(filename,str1);
     strcat(filename,"wb100");
     strcat(filename,str2);
     strcat(filename,".386");
     printf("%s   ",filename);

     readmcs_1(filename,temp);
     fprintf(fpout,"\n%s",filename);
     for(j=0;j<1000;j++)
        pnb[j]=pnb[j]+temp[j];
     }
  /*for(j=0;j<1000;j++)
      fprintf(fpout,"%d  %.3f\n",j,pnb[j]);    */
  }

 fprintf(fpout,"\n\n");
 for(j=0;j<1000;j++)
   fprintf(fpout,"%d  %.4e    %.4e    %.4e     %.4e\n",j,p407[j],p386[j],pna[j],pnb[j]);
fclose(fpout);

/*********************************************************
**
**   subtract background and data smoothing
**
**********************************************************/
    int start_bg, end_bg, gate, nn, newnum;
  float ave407=0, ave386=0, avena=0, avenb=0, spa_res, pp1, pp2, layer;
  float p1[1000],p2[1000],high[1000];



  printf("\n Input bakground range from channel ? to channel ?, such as: 700  900\n");
  scanf("%d%d",&start_bg,&end_bg);
  fflush(stdin);

  for(i=start_bg;i<end_bg;i++)
     {
      ave407=ave407+p407[i]/(end_bg-start_bg);
      ave386=ave386+p386[i]/(end_bg-start_bg);
      avena=avena+pna[i]/(end_bg-start_bg);
      avenb=avenb+pnb[i]/(end_bg-start_bg);
      }

  k=0;
  printf("\n Input start channel No(gate): ?  and spatial resolution(m),such as:30 \n");
  scanf("%d%f", &gate,&spa_res);      /* gate channel start channel */
  layer=spa_res/1000.;
  fflush(stdin);
  for(i=0;i<1000;i++)
     {
       if(i>=gate)
          {  //if(p407[i]<=0 || p386[i]<=0) break;
            high[k]=i*spa_res/1000.;
            p407[k]=p407[i]-ave407;
            p386[k]=p386[i]-ave386;
            pna[k]=pna[i]-avena;
            pnb[k]=pnb[i]-avenb;
 //           fprintf(fsig,"%6.3f    %6.3f    %6.3f\n",high[k],p407[k],p386[k]);
            k++;
         }
     }
  newnum=k;

/* to smooth lidar signal */
 for(i=0;i<newnum;i++)
    {
      p1[i]=p407[i]; p2[i]=p386[i];
      if(high[i]<2.5 ) nn=0.15/layer/2; /* lower areas */
      else nn=0.3/layer/2.0;            /* higher areas */

      if(i>=nn && i<(newnum-nn))
        { pp1=0;  pp2=0;
	       for(j=-nn;j<=nn;j++)
    	      { pp1=pp1+p407[i+j]/(2*nn+1);
              pp2=pp2+p386[i+j]/(2*nn+1);
             }
         p1[i]=pp1;  p2[i]=pp2;
       }
     if( p1[i]<=0 || p2[i]<0) break;
    }
 newnum=i;

 for(i=0;i<newnum;i++)
   {
    ratio1[i]=p1[i]/p2[i];
    p407[i]=p1[i];
    p386[i]=p2[i];
    }
/* to smooth channel sigal */

 for(i=0;i<1000;i++)
  {
   p1[i]=0;p2[i]=0;
   }

 for(i=0;i<newnum;i++)
    {
      p1[i]=pna[i]; p2[i]=pnb[i];
      if(high[i]<2.5 ) nn=0.15/layer/2; /* lower areas */
      else nn=0.3/layer/2.0;            /* higher areas */

      if(i>=nn && i<(newnum-nn))
        { pp1=0;  pp2=0;
	       for(j=-nn;j<=nn;j++)
    	      { pp1=pp1+pna[i+j]/(2*nn+1);
              pp2=pp2+pnb[i+j]/(2*nn+1);
             }
         p1[i]=pp1;  p2[i]=pp2;
       }
     if( p1[i]<=0 || p2[i]<0) break;
    }

 for(i=0;i<newnum;i++)
   if(high[i]<5)
     ratio2[i]=p1[i]/p2[i];

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

 float ave_cal=0, cal_high;

 printf("Input calibration range, such as: 2 or 3 or 4\n");
 scanf("%f",&cal_high);

 j=0;
 for(i=0;i<newnum;i++)
    if(high[i] > cal_high && high[i] < cal_high+1)
     {
       ave_cal=ave_cal+ratio2[i];
       j++  ;
      }
  ave_cal=ave_cal/j;

 strcpy(str3,"");
 strcat(str3,str1);
 strcat(str3,"fact_water.dat");

 if( (fact=fopen(str3,"w"))==NULL)
    {
     printf("\n cannot open %s\n",str3);
     exit(1);
    }
 for(i=0;i<newnum;i++)
   if(high[i]<5)
     {
      ratio2[i]=ratio2[i]/ave_cal;
      fprintf(fact,"%.3f   %.3f\n", high[i],ratio2[i]);
      }
  fclose(fact);


 for(i=0;i<newnum;i++)
    if(high[i]<cal_high)
      ratio1[i]=ratio1[i]/ratio2[i];

//*************************************************************************
//
//   read atmospheric model  "atmod355.dat".
//
//************************************************************************

   int nm;
 float mh[500], molex[500], mdensity, mp, mt, molHmax;

 if( (fmol=fopen("atmod355.dat","r"))==NULL ) // atmos.model with N(z),T(z)
    { printf("\n Cannot open atmod355.dat \n");
      exit(1);
    }

 fscanf(fmol,"%d\n", &nm); /* to read model N2 and T(K) data, ME90TN2.dat */
 for(i=0;i<nm;i++)
  {                                                             //读取氮气消光,密度
   fscanf(fmol,"%f %e %f %f %f\n",&mh[i],&molex[i],&mdensity,&mp,&mt);
  }     /* mh[i]/km,  molex[i]/ km-1  */
 molHmax=mh[i-1];
 fclose(fmol);

/* interpolate molecular data at Lidar height (km) */
 float am[1000];

 for(i=0;i<newnum;i++)
   { for(j=0;j<nm-1;j++)
       { if(high[i]>=mh[j]&& high[i]<mh[j+1])
           {am[i]=molex[j]+( molex[j+1]-molex[j])/(mh[j+1]-mh[j])*(high[i]-mh[j]);
            break;
           }
        }
      if (high[i]>=molHmax) break;
     }
  newnum=i;

 /* read 1998 mean aerosol data     */

                                                           // 首次猜测气溶胶文件名
  if( (faer=fopen("aer355_98.dat","r"))==NULL) /* Aerosol data file */
    { printf("\n Cannot open file aer355_98.dat \n");
      exit(1);                                                                   // !!!!
    }


  float ah[1000],aerex[1000],ap[1000];
   /* Aerosol extinction */

  i=0;
  while( (fscanf(faer,"%f %f\n",&ah[i],&aerex[i]))!=EOF)
      i++;     /* ah[i]/km,  aerex[i]/ km-1  */
 fclose(faer);

 /* interpolate aerosol data at Lidar height (km) */
 for(i=0;i<newnum;i++)
   { for(j=0;j<400;j++)
       { if(high[i]>=ah[j]&& high[i]<ah[j+1])
           {
            ap[i]=aerex[j]+( aerex[j+1]-aerex[j])/(ah[j+1]-ah[j])*(high[i]-ah[j]);
            break;
           }
        }
   }

  /******************************************************************************/
/* to calculate molcular optical depth (0-z) */

   float taum[1000],scale[1000], tmp1, tmp2, taua[1000];
  for(i=0;i<1000;i++)
  {
	  taum[i]=0;
    scale[i]=0;
     taua[i]=0;
  }

  for(i=1;i<newnum;i++)
  {
    taum[i]=taum[i-1]+(am[i]+am[i-1])/2.0*(high[i]-high[i-1]);
    tmp1=taum[i]*(1-pow(354.7/386,4.0))-taum[i]*(1-pow(354.7/407,4.0));

    taua[i]=taua[i-1]+(ap[i]+ap[i-1])/2.0*(high[i]-high[i-1]);
    tmp2=taua[i]*(1-354.7/386)-taua[i]*(1-354.7/407);

    scale[i]=ratio1[i]*exp(-tmp1-tmp2);
  }

/*******************************************/
/* to calculte  water vapor mixing ratio   */
/*******************************************/
  float cons;

  printf(" \n If calibrated by radiosonde, choise '1'.\n");
  printf(" else  with a constant, choise '2'.\n");
  scanf("%d",&k);

  if(k==1)
     printf("ok\n");
    else
     {
      printf("Input calibration constant:?\n");
      scanf("%f",&cons);
      }

  for(i=1;i<newnum;i++)
     {
      scale[i]=scale[i]*cons;
     }

  /*  smooth and output results (0.45-1.0 km)  */

 strcpy(filename,"");
 strcat(filename,str1);
 strcat(filename,"W_Mix_R.dat");
 if( (fr=fopen(filename,"w"))==NULL)
    {
     printf("\n cannot open %s\n", filename);
     exit(1);
    }

 fprintf(fr,"Height/km   Water_vapor_mixing_ratio/g/kg\n");
 for(i=1; i<newnum; i++)
   { tmp1=scale[i];
     if(high[i]>=1.5)
      { if(high[i]>=1.5 && high[i]<=3.0) nn=0.150/layer/2;
        else if(high[i]>3.0 && high[i]<=6.0) nn=0.3/layer/2;
        else nn=0.3/layer/2;
        if(i>=nn&& i<newnum-nn)
      	 { tmp1=0;
	         for(j=-nn;j<=nn;j++)
      	      tmp1=tmp1+scale[i+j]/(2*nn+1);
           }
      }
    scale[i]=tmp1;

    if(high[i]<=13)
      fprintf(fr,"%4.3f   %6.4f\n",high[i],scale[i]);
  }

 fclose(fr);
 printf("\n ! water vapor mixing ratio: %s\n", filename);


 }