write_debug.c

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

        free((char *)EVAP);
      }
      INIT_MOIST  = (double *)calloc(Nveg+1,sizeof(double));
      MOIST_ERROR = (double **)calloc(Nveg+1,sizeof(double*));
      INFLOW      = (double *)calloc(Nveg+1,sizeof(double));
      RUNOFF      = (double *)calloc(Nveg+1,sizeof(double));
      BASEFLOW    = (double *)calloc(Nveg+1,sizeof(double));
      EVAP        = (double *)calloc(Nveg+1,sizeof(double));
      for(i=0;i<=Nveg;i++)
        MOIST_ERROR[i] = (double *)calloc(options.Nlayer+3,sizeof(double));
    }
    if(rec==0 && dist==0) {
      for(band = 0; band < Nbands; band++) {
	INIT_MOIST[veg] = debug.store_moist[WET][band][options.Nlayer+2] * mu;
	INIT_MOIST[veg] += debug.store_moist[DRY][band][options.Nlayer+2] 
	  * (1. - mu);
      }
      for(i=0;i<options.Nlayer+3;i++) MOIST_ERROR[veg][i] = 0.;
    }
    if(dist==0) {
      INFLOW[veg] = 0.;
      RUNOFF[veg] = 0.;
      BASEFLOW[veg] = 0.;
      EVAP[veg] = 0.;
    }
    Evap       = (double *)calloc(options.Nlayer+3,sizeof(double));
    curr_moist = (double *)calloc(options.Nlayer+3,sizeof(double));
    for(band = 0; band < Nbands; band++) {
      if(soil_con->AreaFract[band]>0) {
	Evap[options.Nlayer+2] = 0.;
	curr_moist[options.Nlayer+2] = 0.;
	if(veg < Nveg) {
	  /** Vegetation Present **/
	  Evap[0]       = veg_var[band].canopyevap;
	  curr_moist[0] = veg_var[band].Wdew;
	  Evap[0]       += snow[band].canopy_vapor_flux * 1000.;
	  curr_moist[0] += (snow[band].snow_canopy) * 1000.;
	  Evap[options.Nlayer+2]       += Evap[0];
	  curr_moist[options.Nlayer+2] += curr_moist[0];
	}
	else {
	  /** No vegetation **/
	  Evap[0]       = 0.;
	  curr_moist[0] = 0.;
	}
	
	/** Snow **/
	Evap[1]                       = snow[band].vapor_flux * 1000.;
	Evap[options.Nlayer+2]       += Evap[1];
	curr_moist[1]                 = snow[band].swq * 1000.;
	curr_moist[options.Nlayer+2] += curr_moist[1];
	
	for(i = 0; i < options.Nlayer; i++) {
	  /** All Soil Layers **/
	  Evap[i+2]                     = cell[band].layer[i].evap;
	  Evap[options.Nlayer+2]       += Evap[i+2];
	  curr_moist[i+2]               = cell[band].layer[i].moist;
	  curr_moist[options.Nlayer+2] += curr_moist[i+2];
	}
	
	/** Compute Moisture Balance Error **/
	for(i = 0; i < options.Nlayer+3; i++) {
	  if(dist==0 && band==0) MOIST_ERROR[veg][i] = 0.;
	  MOIST_ERROR[veg][i] 
	    += (debug.inflow[dist][band][i] 
		- (debug.outflow[dist][band][i] + Evap[i]) 
		- (curr_moist[i] - debug.store_moist[dist][band][i])) 
	    * mu * soil_con->AreaFract[band];
	  if(fabs(MOIST_ERROR[veg][i]) > 1.e-4) {
	    fprintf(stderr,"WARNING: Debug Layer %i has a Moisture Balance Error of %f in rec %i, veg %i, band %i, precip dist %i\n",i,MOIST_ERROR[veg][i],rec,veg,band,dist);
	  }
	}

	/** Store Variables **/
	INFLOW[veg]   += atmos->prec[NR] * soil_con->Pfactor[band] 
	  / soil_con->AreaFract[band];
	BASEFLOW[veg] += cell[band].baseflow * mu / soil_con->AreaFract[band];
	RUNOFF[veg]   += cell[band].runoff * mu / soil_con->AreaFract[band];
	EVAP[veg]     += Evap[options.Nlayer+2] * mu / soil_con->AreaFract[band];
      }
    }

    if(rec==0 && veg==0 && dist==0) {
      fprintf(debug.fg_balance,"Date\tVeg Num\tPrecip\tRunoff");
      fprintf(debug.fg_balance,"\tBFlow\tEvap\tdStor\tError\n");
    }
    if((options.DIST_PRCP && dist==1) || !options.DIST_PRCP) {
      fprintf(debug.fg_balance,"%f\t%i\t%f\t%f\t%f",
	      (double)rec/24.0*(double)gp->dt,veg,INFLOW[veg],RUNOFF[veg],
	      BASEFLOW[veg]);
      fprintf(debug.fg_balance,"\t%f\t%f\t%f\n",
          EVAP[veg],curr_moist[options.Nlayer+2]-INIT_MOIST[veg],
          MOIST_ERROR[veg][options.Nlayer+2]);
    }
    free((char*)Evap);
    free((char*)curr_moist);
  }
 
  if(debug.PRT_TEMP) {
  
    /***** Temperature Profile Debugging Output *****/
 
    if(rec==0 && veg==0 && dist==0) {
      fprintf(debug.fg_temp,"%i\n",options.Nnode);
      fprintf(debug.fg_temp,"Date - hour(REC)\tveg\tband\tAir T\tFdpth\tTdpth");
      for(i=0;i<options.Nlayer;i++) 
	fprintf(debug.fg_temp,"\tLayer %i",i);
      sum=0.0;
      for(i=0;i<options.Nnode;i++) {
	fprintf(debug.fg_temp,"\tT%.0f",sum*100.0);
	sum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.0;
      }
      sum=0.0;
      for(i=0;i<options.Nnode;i++) {
	fprintf(debug.fg_temp,"\tM%.0f",sum*100.0);
	sum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.0;
      }
      sum=0.0;
      for(i=0;i<options.Nnode;i++) {
	fprintf(debug.fg_temp,"\tI%.0f",sum*100.0);
	sum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.0;
      }
      sum=0.0;
      for(i=0;i<options.Nnode;i++) {
	fprintf(debug.fg_temp,"\tK%.0f",sum*100.0);
	sum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.0;
      }
      sum=0.0;
      for(i=0;i<options.Nnode;i++) {
	fprintf(debug.fg_temp,"\tCs%.0f",sum*100.0);
	sum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.0;
      }
      fprintf(debug.fg_temp,"\n");
    }
    for(band = 0; band < Nbands; band++) {
      if(soil_con->AreaFract[band]>0) {
	fprintf(debug.fg_temp,"%02i/%02i/%04i - %02i\t%7f\t%i\t%i",
		dmy->day, dmy->month, dmy->year, dmy->hour,
		(float)rec/24.0*(float)gp->dt, veg, band);
	fprintf(debug.fg_temp,"\t%6.2f\t%6.2f\t%6.2f",
		atmos->air_temp[NR], energy->fdepth[0] * 100., 
		energy->tdepth[0] * 100.);
	/* print layer temperatures */
	for(i = 0; i < options.Nlayer; i++)
	  fprintf(debug.fg_temp,"\t%6.2f", cell[band].layer[i].T);
	/* print node temperatures */
	for(i = 0; i < options.Nnode; i++) 
	  fprintf(debug.fg_temp,"\t%6.2f", energy->T[i]);
	/* print node moisture contents */
	for(i = 0; i < options.Nnode; i++) 
	  fprintf(debug.fg_temp,"\t%6.4f", energy->moist[i]);
	/* print node ice contents */
	for(i = 0; i < options.Nnode; i++) 
	  fprintf(debug.fg_temp,"\t%6.4f", energy->ice[i]);
	/* print node thermal conductivities */
	for(i = 0; i < options.Nnode; i++) 
	  fprintf(debug.fg_temp,"\t%6.2f", energy->kappa_node[i]);
	/* print node volumetric heat capacities */
	for(i = 0;i < options.Nnode; i++) 
	  fprintf(debug.fg_temp,"\t%6.0f", energy->Cs_node[i]);
	fprintf(debug.fg_temp,"\n");
      }
    }
 
  }
 
  if(debug.PRT_MOIST) {

    /***** Moisture Profile Debugging Output *****/

    if(FIRST != -999) {
      fprintf(debug.fg_moist,"Date - hour(REC)        \tVeg Num\tDist Num");
      fprintf(debug.fg_moist,"\tT Air");
      fprintf(debug.fg_moist,"\tInflow\tRunoff");
      for(i=0;i<options.Nlayer;i++)
        fprintf(debug.fg_moist,"\t%i Moist\t%i Ice",i,i);
      fprintf(debug.fg_moist,"\n");
    }

    for(band = 0; band < Nbands; band++) {
      if(soil_con->AreaFract[band]>0) {

	fprintf(debug.fg_moist,"%02i/%02i/%04i - %02i\t%7f\t%i\t%i",
		dmy->day,dmy->month,dmy->year,dmy->hour,
		(float)rec/24.0*(float)gp->dt,veg,dist);
	fprintf(debug.fg_moist,"\t%6.2f\t%6.4f\t%6.4f",
		atmos->air_temp[NR], cell[band].inflow, cell[band].runoff);
	
	curr_moist = (double *)calloc(1,sizeof(double));
	for(i = 0; i < options.Nlayer; i++) {
	  curr_moist[0] = cell[band].layer[i].moist * (soil_con->depth[i]) 
	    / soil_con->depth[i];
	  fprintf(debug.fg_moist,"\t%6.4f",
		  curr_moist[0] / soil_con->depth[i] / 1000.);
	}
	fprintf(debug.fg_moist,"\n");
	free((char*)curr_moist);
      }
    }
  }
 
  if(debug.PRT_KAPPA) {
 
    /***** Soil Thermal Properties Profile Debugging Output *****/

    fprintf(debug.fg_kappa,"%02i/%02i/%04i - %02i\t%7f\t%i",
        dmy->day,dmy->month,dmy->year,dmy->hour,
        (float)rec/24.0*(float)gp->dt,veg);
 
    for(band = 0; band < Nbands; band++) {
      if(soil_con->AreaFract[band]>0) {

	for(i=0;i<options.Nlayer;i++) {
	  fprintf(debug.fg_kappa,"\t%6.2f\t%6.0f",
		  cell[band].layer[i].kappa,
		  cell[band].layer[i].Cs);
	}
	fprintf(debug.fg_kappa,"\n");
      }
    }
  }
 
  if(debug.PRT_GRID) {

    /***** Soil Themperature GMT Grided Profile Output *****/
 
    Zsum=0.;
    for(i=0;i<options.Nnode;i++) {
      fprintf(debug.fg_grid,"%7f\t%f\t%f\n",
          (float)rec/24.0*(float)gp->dt,
          Zsum,energy->T[i]);
      if(i<options.Nnode)
        Zsum+=(soil_con->dz_node[i]+soil_con->dz_node[i+1])/2.;
    }
 
  }

  FIRST = -999;
 
}
#endif