write_data.c

超强的大尺度水文模拟工具

      fwrite(tmp_fptr,1,sizeof(float),outfiles->snow);
    }
  }
  else if((options.FULL_ENERGY || options.FROZEN_SOIL)) {
    /***** Write ASCII full energy snow output file *****/
    fprintf(outfiles->snow  ,"%04i\t%02i\t%02i\t%02i\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\n",
	    dmy->year, dmy->month, dmy->day, dmy->hour,
	    out_data->swq[0], out_data->snow_depth[0], 
	    out_data->snow_canopy[0], /* out_data->coverage[0], */ 
	    out_data->advection[0], out_data->deltaCC[0], 
	    out_data->snow_flux[0], out_data->refreeze_energy[0]);
  }
  else {
    /***** Write ASCII water balance snow output file *****/
    fprintf(outfiles->snow  ,"%04i\t%02i\t%02i",
	    dmy->year, dmy->month, dmy->day);
    if(dt<24) 
      fprintf(outfiles->snow  ,"\t%02i",dmy->hour);
    fprintf(outfiles->snow  ,"\t%.4f\t%.4f\t%.4f\n",
	    out_data->swq[0], out_data->snow_depth[0], 
	    out_data->snow_canopy[0]);
  }

  /**************************************
    Ouput Snow Band Variables
  **************************************/
  if(options.PRT_SNOW_BAND && options.BINARY_OUTPUT) {
    /***** Write Binary Snow Band Output File *****/
    tmp_iptr[0] = dmy->year;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->snowband);
    tmp_iptr[0] = dmy->month;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->snowband);
    tmp_iptr[0] = dmy->day;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->snowband);
    if(dt<24) {
      tmp_iptr[0] = dmy->hour;
      fwrite(tmp_iptr,1,sizeof(int),outfiles->snowband);
    }
    for(band=0;band<options.SNOW_BAND;band++) {
      tmp_fptr[0] = (float)out_data->swq[band+1];
      fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
      tmp_fptr[0] = (float)out_data->snow_depth[band+1];
      fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
      tmp_fptr[0] = (float)out_data->snow_canopy[band+1];
      fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
      /*     tmp_fptr[0] = (float)out_data->coverage; */
      /*     fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband); */
      if(options.FULL_ENERGY || options.FROZEN_SOIL) {
	tmp_fptr[0] = (float)out_data->advection[band+1];
	fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
	tmp_fptr[0] = (float)out_data->deltaCC[band+1];
	fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
	tmp_fptr[0] = (float)out_data->snow_flux[band+1];
	fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
	tmp_fptr[0] = (float)out_data->refreeze_energy[band+1];
	fwrite(tmp_fptr,1,sizeof(float),outfiles->snowband);
      }
    }
  }
  else if((options.FULL_ENERGY || options.FROZEN_SOIL) 
	  && options.PRT_SNOW_BAND) {
    /***** Write ASCII full energy snow band output file *****/
    fprintf(outfiles->snowband  ,"%04i\t%02i\t%02i\t%02i",
	    dmy->year, dmy->month, dmy->day, dmy->hour);
    for(band=0;band<options.SNOW_BAND;band++) {
      fprintf(outfiles->snowband  ,"\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f",
	      out_data->swq[band+1], out_data->snow_depth[band+1], 
	      out_data->snow_canopy[band+1], /* out_data->coverage[band+1], */ 
	      out_data->advection[band+1], out_data->deltaCC[band+1], 
	      out_data->snow_flux[band+1], out_data->refreeze_energy[band+1]);
      fprintf(outfiles->snowband,"\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f", 
	      out_data->swband[band], out_data->lwband[band], 
	      out_data->albedoband[band], out_data->latentband[band], 
	      out_data->sensibleband[band], out_data->grndband[band]);
    }
    fprintf(outfiles->snowband,"\n");
  }
  else if(!options.FULL_ENERGY && options.PRT_SNOW_BAND) {
    /***** Write ASCII water balance snow band output file *****/
    fprintf(outfiles->snowband  ,"%04i\t%02i\t%02i",
	      dmy->year, dmy->month, dmy->day);
    if(dt<24) 
      fprintf(outfiles->snowband  ,"\t%02i", dmy->hour);
    for ( band = 0; band < options.SNOW_BAND; band++ ) {
      fprintf(outfiles->snowband  ,"\t%.4f\t%.4f\t%.4f",
	      out_data->swq[band+1], out_data->snow_depth[band+1], 
	      out_data->snow_canopy[band+1]); 
      fprintf(outfiles->snowband,"\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f", 
	      out_data->swband[band], out_data->lwband[band], 
	      out_data->albedoband[band], out_data->latentband[band], 
	      out_data->sensibleband[band], out_data->grndband[band]);
    }
    fprintf(outfiles->snowband,"\n");
  }

  /************************************
    Output Standard Energy and Moisture Flux Variables
  ************************************/
  if(options.BINARY_OUTPUT) {
    /***** Write Binary Fluxes File *****/
    tmp_iptr[0] = dmy->year;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->fluxes);
    tmp_iptr[0] = dmy->month;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->fluxes);
    tmp_iptr[0] = dmy->day;
    fwrite(tmp_iptr,1,sizeof(int),outfiles->fluxes);
    if(dt<24) {
      tmp_iptr[0] = dmy->hour;
      fwrite(tmp_iptr,1,sizeof(int),outfiles->fluxes);
    }
    tmp_fptr[0] = (float)out_data->prec;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->evap;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->runoff;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->baseflow;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->Wdew;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    for(j=0;j<options.Nlayer;j++) {    
      tmp_fptr[0] = (float)out_data->moist[j];
      fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    }
    if(options.FULL_ENERGY || options.FROZEN_SOIL) {
      tmp_fptr[0] = (float)out_data->rad_temp;
      fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    }
    tmp_fptr[0] = (float)out_data->net_short;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->r_net;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    if(options.FULL_ENERGY || options.FROZEN_SOIL) {
      tmp_fptr[0] = (float)out_data->latent;
      fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    }
    tmp_fptr[0] = (float)out_data->evap_canop;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->evap_veg;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->evap_bare;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->sub_canop;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->sub_snow;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    if(options.FULL_ENERGY || options.FROZEN_SOIL) {
      tmp_fptr[0] = (float)out_data->sensible;
      fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
      tmp_fptr[0] = (float)out_data->grnd_flux;
      fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    }
    tmp_fptr[0] = (float)out_data->aero_resist;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->surf_temp;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
    tmp_fptr[0] = (float)out_data->albedo;
    fwrite(tmp_fptr,1,sizeof(float),outfiles->fluxes);
  }
  else if(options.FULL_ENERGY || options.FROZEN_SOIL) {
    /***** Write ASCII energy balance fluxes file *****/
    fprintf(outfiles->fluxes,"%04i\t%02i\t%02i\t%02i\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f",
	    dmy->year, dmy->month, dmy->day, dmy->hour,
	    out_data->prec, out_data->evap, out_data->runoff,
	    out_data->baseflow, out_data->Wdew);
    for(j=0;j<options.Nlayer;j++) 
      fprintf(outfiles->fluxes,"\t%.4f", out_data->moist[j]);
    fprintf(outfiles->fluxes,"\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t\n",
	    out_data->rad_temp,
	    out_data->net_short, out_data->r_net, out_data->latent,
	    out_data->evap_canop, out_data->evap_veg,
	    out_data->evap_bare, out_data->sub_canop, 
	    out_data->sub_snow, out_data->sensible,
	    out_data->grnd_flux, out_data->aero_resist, out_data->surf_temp, 
	    out_data->albedo);
  }
  else {
    /***** Write ASCII Water Balance Fluxes File *****/
    fprintf(outfiles->fluxes,"%04i\t%02i\t%02i",
	    dmy->year, dmy->month, dmy->day);
    if(dt<24)
      fprintf(outfiles->fluxes,"\t%02i", dmy->hour);
    fprintf(outfiles->fluxes,"\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f",
	    out_data->prec, out_data->evap, out_data->runoff, 
	    out_data->baseflow, out_data->Wdew);
    for(j=0;j<options.Nlayer;j++) 
      fprintf(outfiles->fluxes,"\t%.4f", out_data->moist[j]);
    fprintf(outfiles->fluxes,"\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t\n",
	    out_data->net_short, out_data->r_net, out_data->evap_canop, 
	    out_data->evap_veg, out_data->evap_bare, 
	    out_data->sub_canop, out_data->sub_snow, 
	    out_data->aero_resist, out_data->surf_temp, out_data->albedo);
  }

  free((char *)tmp_iptr);
  free((char *)tmp_fptr);

#endif

}

void calc_water_balance_error(int    rec,
			      double inflow,
			      double outflow,
			      double storage) {
/***************************************************************
  calc_water_balance_error  Keith Cherkauer        April 1998

  This subroutine computes the overall model water balance, and 
  warns the model user if large errors are found.
***************************************************************/

  static double last_storage;
  static double cum_error;
  static double max_error;
  static int    error_cnt;
  static int    Nrecs;

  double error;

  if(rec<0) {
    last_storage = storage;
    cum_error    = 0.;
    max_error    = 0.;
    error_cnt    = 0;
    Nrecs        = -rec;
  }
  else {
    error = inflow - outflow - (storage - last_storage);
    cum_error += error;
    if(fabs(error)>fabs(max_error) && fabs(error)>1e-5) {
      max_error = error;
      fprintf(stderr,"Maximum Moist Error:\t%i\t%.5f\t%.5f\n",
	      rec,error,cum_error);
    }
    if(rec==Nrecs-1) {
      fprintf(stderr,"Total Cumulative Water Error for Grid Cell = %.4f\n",
	      cum_error);
    }
    last_storage = storage;
  }

}

void calc_energy_balance_error(int    rec,
			       double net_rad,
			       double latent,
			       double sensible,
			       double grnd_flux,
			       double snow_fluxes) {
/***************************************************************
  calc_energy_balance_error   Keith Cherkauer     April 1998

  This subroutine computes the overall model energy balance, and
  reports the maximum time step error above a thresehold to the
  user.  The total cumulative error for the grid cell is also 
  computed and reported at the end of the model run.
***************************************************************/

  static double cum_error;
  static double max_error;
  static int    Nrecs;

  double error;

  if(rec<0) {
    cum_error = 0;
    Nrecs     = -rec;
    max_error = 0;
  }
  else {
    error = net_rad - latent - sensible - grnd_flux + snow_fluxes;
    cum_error += error;
    if(fabs(error)>fabs(max_error) && fabs(error)>0.001) {
      max_error = error;
      if ( rec > 0 ) 
	fprintf(stderr,"Maximum Energy Error:\t%i\t%.4f\t%.4f\n",
		rec,error,cum_error/(double)rec);
      else 
	fprintf(stderr,"Maximum Energy Error:\t%i\t%.4f\t%.4f\n",
		rec,error,cum_error);
    }
    if(rec==Nrecs-1) {
      fprintf(stderr,"Total Cumulative Energy Error for Grid Cell = %.4f\n",
	      cum_error/(double)rec);
    }
  }
}