write_debug.c

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <vicNl.h>

static char vcid[] = "$Id: write_debug.c,v 4.1 2000/05/16 21:07:16 vicadmin Exp $";

#if LINK_DEBUG
void write_debug(atmos_data_struct    *atmos,
                 soil_con_struct      *soil_con,
                 cell_data_struct     *cell,
                 energy_bal_struct    *energy,
                 snow_data_struct     *snow,
                 veg_var_struct       *veg_var,
                 dmy_struct           *dmy,
                 global_param_struct  *gp,
                 double                out_short,
                 double                mu,
                 int                   Nveg,
                 int                   veg,
                 int                   rec,
                 int                   gridcell,
                 int                   dist,
                 char                  NEWCELL) {
/**********************************************************************
  write_debug		Keith Cherkauer		October 8, 1997

  This subroutine controls the output of selected debug data files.
  Which debugging files are output is controlled by flags set in the
  model control file.

  This subroutine is not essential to the operation of the model, and 
  can be excluded from the compiled code as long as references to it 
  from full_energy.c are removed.  Due to the number of static variables
  in this routine, it may contribut significantly to the size of the
  compiled code.

  Modifications:
  07-15-98 modified to work with elevation bands                 KAC

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

  extern option_struct options;
  extern debug_struct debug;

  static short int   FIRST;
  static double    **MOIST_ERROR;
  static double     *INIT_MOIST;
  static double     *ENERGY_ERROR;
  static double     *ENERGY_ERROR_CALC;
  static double     *INFLOW;
  static double     *RUNOFF;
  static double     *BASEFLOW;
  static double     *EVAP;
  static double     *INSHORT;
  static double     *OUTSHORT;
  static double     *INLONG;
  static double     *OUTLONG;
  static double     *SENSIBLE;
  static double     *LATENT;
  static double     *GRND_FLUX;
  static double     *ADVECTION;
  static double     *DELTA_CC;
  static double     *SNOW_FLUX;
  static double     *REFREEZEENERGY;
  static double     *DELTA_H;

  int     i;
  int     Ntemp;
  int     band;
  int     Nbands;
  double  sum;
  double *Evap;
  double *curr_moist;
  double  curr_ice;
  double  Zsum;
  double  grnd_flux;
  double  advection;
  double  deltaH;
  double  deltaCC;
  double  Qnet;

  Nbands = options.SNOW_BAND;

  if(debug.PRT_FLUX && options.FULL_ENERGY) {

    /***** Record Hourly Energy Balance Terms *****/

    if(NEWCELL && dist==0) {
      if(gridcell>0) {
        free((char *)ENERGY_ERROR);
        free((char *)ENERGY_ERROR_CALC);
        free((char *)INSHORT);
        free((char *)OUTSHORT);
        free((char *)INLONG);
        free((char *)OUTLONG);
        free((char *)SENSIBLE);
        free((char *)LATENT);
        free((char *)GRND_FLUX);
        free((char *)ADVECTION);
        free((char *)DELTA_H);
        free((char *)DELTA_CC);
        free((char *)SNOW_FLUX);
        free((char *)REFREEZEENERGY);
      }
      ENERGY_ERROR      = (double *)calloc(Nveg+1,sizeof(double));
      ENERGY_ERROR_CALC = (double *)calloc(Nveg+1,sizeof(double));
      INSHORT           = (double *)calloc(Nveg+1,sizeof(double));
      OUTSHORT          = (double *)calloc(Nveg+1,sizeof(double));
      INLONG            = (double *)calloc(Nveg+1,sizeof(double));
      OUTLONG           = (double *)calloc(Nveg+1,sizeof(double));
      SENSIBLE          = (double *)calloc(Nveg+1,sizeof(double));
      LATENT            = (double *)calloc(Nveg+1,sizeof(double));
      GRND_FLUX         = (double *)calloc(Nveg+1,sizeof(double));
      ADVECTION         = (double *)calloc(Nveg+1,sizeof(double));
      DELTA_H           = (double *)calloc(Nveg+1,sizeof(double));
      DELTA_CC          = (double *)calloc(Nveg+1,sizeof(double));
      SNOW_FLUX         = (double *)calloc(Nveg+1,sizeof(double));
      REFREEZEENERGY    = (double *)calloc(Nveg+1,sizeof(double));
    }
    if(rec==0 && dist==0) {
      ENERGY_ERROR[veg]      = 0.;
      ENERGY_ERROR_CALC[veg] = 0.;
    }
    if(dist==0) {
      INSHORT[veg]        = 0.;
      OUTSHORT[veg]       = 0.;
      INLONG[veg]         = 0.;
      OUTLONG[veg]        = 0.;
      SENSIBLE[veg]       = 0.;
      LATENT[veg]         = 0.;
      GRND_FLUX[veg]      = 0.;
      ADVECTION[veg]      = 0.;
      DELTA_H[veg]        = 0.;
      DELTA_CC[veg]       = 0.;
      SNOW_FLUX[veg]      = 0.;
      REFREEZEENERGY[veg] = 0.;
    } 

    for(band = 0; band < Nbands; band++) {
      deltaCC   = energy[band].deltaCC;
      Qnet      = snow[band].Qnet / soil_con->AreaFract[band]; 
      advection = energy[band].advection;
      deltaH    = energy[band].deltaH;
      
      ENERGY_ERROR[veg] += energy[band].error * mu;
      ENERGY_ERROR_CALC[veg] 
	+= ((1.-energy[band].albedo)*energy[band].shortwave
	    + energy[band].longwave + energy[band].grnd_flux
	    + energy[band].latent + energy[band].sensible 
	    + energy[band].deltaH - energy[band].deltaCC 
	    - energy[band].snow_flux + energy[band].refreeze_energy 
	    + energy[band].advection) * mu / soil_con->AreaFract[band];
      
      INSHORT[veg]        += (1.-energy[band].albedo)*atmos->shortwave[NR] 
	* mu / soil_con->AreaFract[band];
      INLONG[veg]         += atmos->longwave[NR] * mu 
	/ soil_con->AreaFract[band];
      SENSIBLE[veg]       += energy[band].sensible * mu 
	/ soil_con->AreaFract[band];
      LATENT[veg]         += energy[band].latent * mu 
	/ soil_con->AreaFract[band];
      GRND_FLUX[veg]      += energy[band].grnd_flux * mu 
	/ soil_con->AreaFract[band];
      ADVECTION[veg]      += energy[band].advection * mu 
	/ soil_con->AreaFract[band];
      DELTA_H[veg]        += energy[band].deltaH * mu 
	/ soil_con->AreaFract[band];
      DELTA_CC[veg]       += energy[band].deltaCC * mu 
	/ soil_con->AreaFract[band];
      SNOW_FLUX[veg]      += energy[band].snow_flux * mu 
	/ soil_con->AreaFract[band];
      REFREEZEENERGY[veg] += energy[band].refreeze_energy * mu 
	/ soil_con->AreaFract[band];
    }

    if(rec==0 && veg==0 && dist==0) {
      fprintf(debug.fg_energy,"DATE\tNET SHT\tNET LNG\t");
      fprintf(debug.fg_energy,"GRND F\tLATENT\tSENSBL\tADVEC\tdel H\t");
      fprintf(debug.fg_energy,"del CC\tSNWFLX\tMELT\t");
      fprintf(debug.fg_energy,"ERROR\tERR CAL\tGRND T\tT_1\tWIND\n");
    }
    if((options.DIST_PRCP && dist==1) || !options.DIST_PRCP) {
      fprintf(debug.fg_energy,"%7.4f\t%7.4f\t%7.4f",
	      (float)rec/24.0*(float)gp->dt, INSHORT[veg],
	      INLONG[veg]);
      fprintf(debug.fg_energy,"\t%7.4f\t%7.4f\t%7.4f",
	      -GRND_FLUX[veg], LATENT[veg], SENSIBLE[veg]);
      fprintf(debug.fg_energy,"\t%7.4f\t%7.4f\t%7.4f\t%7.4f\t%7.4f",
	      ADVECTION[veg], DELTA_H[veg], DELTA_CC[veg], SNOW_FLUX[veg], 
	      REFREEZEENERGY[veg]);
      fprintf(debug.fg_energy,"\t%7.4f\t%7.4f",
	      ENERGY_ERROR[veg],ENERGY_ERROR_CALC[veg]);
      fprintf(debug.fg_energy,"\t%7.4f\t%7.4f\t%7.4f\n",
	      energy[0].T[0], energy[0].T[1],
	      atmos->wind[NR]);
    }
  }
 
  if(debug.PRT_SNOW) {

    /***** Record Hourly Snow Terms *****/
    
    for(band = 0; band < Nbands; band++) {
      
      if(snow[band].snow) grnd_flux = energy->grnd_flux;
      else grnd_flux = 0.;
 
      if(rec==0 && veg==0 && dist==0 && band==0) {
	/** Print File Header **/
	fprintf(debug.fg_snow,"Date\tBand\tSWE TOT\tSWE SRF\tSWE PCK\t");
	fprintf(debug.fg_snow,"GRND T\tlyr1 T\tSURF T\tPACK T\tMELT\t");
	fprintf(debug.fg_snow,"VPR FLX\tAIR T\tSNOW\tRAIN\tGRNDFLX\t");
	fprintf(debug.fg_snow,"DEPTH\tKAPPA\tCANOPY\tCNPYFLUX\n");
      }
      fprintf(debug.fg_snow,"%7.4f\t%7.4f\t%7.4f\t%7.4f\t%7.3f",
	      (float)rec/24.0*(float)gp->dt,snow[band].swq*1000.,
	      snow[band].surf_water*1000.,
	      snow[band].pack_water*1000.,energy->T[0]);
      fprintf(debug.fg_snow,"\t%7.3f\t%7.3f\t%7.3f\t%7.3f",
        energy->T[1],snow[band].surf_temp,snow[band].pack_temp,
        snow[band].vapor_flux*1000.);
      fprintf(debug.fg_snow,"\t%7.3f\t%7.4f",
	      atmos->air_temp[NR], grnd_flux);
      fprintf(debug.fg_snow,"\t%7.4f\t%7.4f\t%7.4f\t%7.4f\n",
	      snow[band].depth,snow[band].density,
	      snow[band].snow_canopy*1000.,
	      snow[band].canopy_vapor_flux*1000.);
    }
  }
 
  if(debug.PRT_BALANCE) {
 
    /***** Compute Water Balance Error *****/

    if(NEWCELL && dist==0) {
      if(gridcell>0) {
        for(i=0;i<=Nveg;i++) free((char *)MOIST_ERROR[i]);
        free((char *)INIT_MOIST);
        free((char *)MOIST_ERROR);
        free((char *)INFLOW);
        free((char *)RUNOFF);
        free((char *)BASEFLOW);