put_data.c

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

	  
	  /**********************************
            Record Energy Balance Variables
	  **********************************/

	  /** record surface radiative temperature **/
	  if(snow[veg][band].swq>0) 
	    rad_temp = snow[veg][band].surf_temp + KELVIN;
	  else
	    rad_temp = energy[veg][band].T[0] + KELVIN;
	  
	  /** record soil surface temperature **/
	  surf_temp = energy[veg][band].T[0];
	  // MODIFIED FOR USE WITH ROSEMOUNT SIMULATIONS
	  //surf_temp = (energy[veg][band].T[0] + energy[veg][band].T[1])/2.;
	  
	  /** record net shortwave radiation **/
	  out_data->net_short += energy[veg][band].shortwave
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->swband[band] = energy[veg][band].shortwave
	      * Cv;
	  
	  /** record net longwave radiation **/
	  out_data->net_long  += energy[veg][band].longwave
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->lwband[band] = energy[veg][band].longwave
	      * Cv;
	  
	  /** record incoming longwave radiation **/
	  out_data->in_long  += ((energy[veg][band].longwave + STEFAN_B 
				  * (rad_temp) * (rad_temp)
				  * (rad_temp) * (rad_temp))
				 * Cv * AreaFract[band] * TreeAdjustFactor[band]);
	  
	  /** record albedo **/
	  out_data->albedo    += energy[veg][band].albedo
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->albedoband[band] = energy[veg][band].albedo
	      * Cv;
	  
	  /** record latent heat flux **/
	  out_data->latent    -= energy[veg][band].latent
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->latentband[band] = energy[veg][band].latent
	      * Cv;

	  /** record sensible heat flux **/
	  out_data->sensible  -= energy[veg][band].sensible
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->sensibleband[band] = energy[veg][band].sensible
	      * Cv;

	  /** record ground heat flux (+ heat storage) **/
	  out_data->grnd_flux -= (energy[veg][band].grnd_flux
				  + energy[veg][band].deltaH)
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  if(options.PRT_SNOW_BAND) 
	    out_data->grndband[band] = -1.*(energy[veg][band].grnd_flux
					    + energy[veg][band].deltaH)
	      * Cv;

	  /** record heat storage **/
	  out_data->deltaH    -= energy[veg][band].deltaH
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record energy balance error **/
	  out_data->energy_error += energy[veg][band].error
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record radiative effective temperature [K], 
	      emissivities set = 1.0  **/
	  out_data->rad_temp += ((rad_temp) * (rad_temp) 
				 * (rad_temp) * (rad_temp)) 
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record mean surface temperature [C]  **/
	  out_data->surf_temp += surf_temp * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /*****************************
	    Record Snow Pack Variables 
	  *****************************/
	  
	  /** record snow water equivalence **/
	  out_data->swq[0]         
	    += snow[veg][band].swq * Cv * 1000. * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record snowpack depth **/
	  out_data->snow_depth[0]  
	    += snow[veg][band].depth * Cv * 100. * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record canopy intercepted snow **/
	  if ( veg < veg_con[0].vegetat_type_num )
#if LDAS_OUTPUT
	    out_data->swq[0] 
#else
	      out_data->snow_canopy[0] 
#endif
	      += (snow[veg][band].snow_canopy) 
	      * Cv * 1000. * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record snow cover fraction **/
	  out_data->coverage[0]    
	    += snow[veg][band].coverage * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record snowpack cold content **/
	  out_data->deltaCC[0]           
	    += energy[veg][band].deltaCC * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record snowpack advection **/
	  out_data->advection[0]         
	    += energy[veg][band].advection * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record snow energy flux **/
	  out_data->snow_flux[0]         
	    += energy[veg][band].snow_flux * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** record refreeze energy **/
	  out_data->refreeze_energy[0]   
	    += energy[veg][band].refreeze_energy 
	    * Cv * AreaFract[band] * TreeAdjustFactor[band];
	  
	  /** if snow elevation bands are to be printed separately **/
	  if(options.PRT_SNOW_BAND) {
	    
	    /** record band snow water equivalent **/
	    out_data->swq[band+1]         
	      += snow[veg][band].swq * Cv  * 1000.;
	    
	    /** record band snowpack depth **/
	    out_data->snow_depth[band+1]  
	      += snow[veg][band].depth * Cv * 100.;
	    
	    /** record band canopy intercepted snow **/
	    if ( veg < veg_con[0].vegetat_type_num )
#if LDAS_OUTPUT
	      out_data->swq[band+1]
#else
		out_data->snow_canopy[band+1]  
#endif
		+= (snow[veg][band].snow_canopy) * Cv * 1000.;
	    
	    /** record band snow coverage **/
	    out_data->coverage[band+1]    
	      += snow[veg][band].coverage * Cv;
	    
	    /** record band cold content **/
	    out_data->deltaCC[band+1]           
	      += energy[veg][band].deltaCC * Cv;
	    
	    /** record band advection **/
	    out_data->advection[band+1]         
	      += energy[veg][band].advection * Cv;
	    
	    /** record band snow flux **/
	    out_data->snow_flux[band+1]         
	      += energy[veg][band].snow_flux * Cv;
	    
	    /** record band refreeze energy **/
	    out_data->refreeze_energy[band+1]   
	      += energy[veg][band].refreeze_energy * Cv;
	    
	  }
	}
      }

#endif /* not OPTIMIZE */

    }
  }
  
#if !OPTIMIZE

  /** record radiative temperature **/
  out_data->rad_temp = pow(out_data->rad_temp,0.25);

  /** record net radiation **/
  out_data->r_net    = out_data->net_short + out_data->net_long;

  /********************
    Check Water Balance 
    ********************/
  inflow  = out_data->prec;
  outflow = out_data->evap+out_data->runoff+out_data->baseflow;
  storage = 0.;
  for(index=0;index<options.Nlayer;index++)
    if(options.MOISTFRACT)
      storage += (out_data->moist[index] + out_data->ice[index]) 
	* depth[index] * 1000;
    else
      storage += out_data->moist[index] + out_data->ice[index];
  storage += out_data->swq[0] + out_data->snow_canopy[0] + out_data->Wdew;
  calc_water_balance_error(rec,inflow,outflow,storage);
  if(options.FULL_ENERGY)
    calc_energy_balance_error(rec, out_data->net_short + out_data->net_long,
			      out_data->latent, out_data->sensible,
			      out_data->grnd_flux, out_data->advection[0] 
			      - out_data->deltaCC[0] - out_data->snow_flux[0]
			      + out_data->refreeze_energy[0]);

  /*************
    Write Data
  *************/

  if(rec >= skipyear)
    write_data(out_data, outfiles, dmy, dt);

#else

  if ( rec == 0 ) prtdt = 0;
  if ( prtdt == 0 ) {
    runoff = out_data->runoff;
    baseflow = out_data->baseflow;
    prtdt ++;
  }
  else {
    runoff += out_data->runoff;
    baseflow += out_data->baseflow;
    prtdt ++;
  }
  if ( prtdt == 24 / dt ) {
    out_data->runoff = runoff;
    out_data->baseflow = baseflow;
    write_data(out_data, outfiles, dmy, dt);
    prtdt = 0;
  } 

#endif

  free((char *)out_data); 

}