📄 calc_surf_energy_bal.c
字号:
#include <stdio.h>#include <stdlib.h>#include "vicNl.h"static char vcid[] = "$Id: calc_surf_energy_bal.c,v 4.2.2.2 2004/09/22 00:40:32 vicadmin Exp $";double calc_surf_energy_bal(int rec, int iveg, int nlayer, int Nveg, int dt, int Nnodes, int veg_class, int band, int hour, double ice0, double moist, double dp, double surf_atten, double T0, double shortwave, double longwave, double air_temp, double Le, double Ls, double mu, double displacement, double roughness, double ref_height, double snow_energy, double vapor_flux, double bare_albedo, double *aero_resist, double *wind, double *rainfall, double *ppt, float *root, atmos_data_struct *atmos, veg_var_struct *veg_var_wet, veg_var_struct *veg_var_dry, energy_bal_struct *energy, snow_data_struct *snow, layer_data_struct *layer_wet, layer_data_struct *layer_dry, soil_con_struct *soil_con, dmy_struct *dmy)/************************************************************** calc_surf_energy_bal.c Greg O'Donnell and Keith Cherkauer Sept 9 1997 This function calculates the surface temperature, in the case of no snow cover. Evaporation is computed using the previous ground heat flux, and then used to comput latent heat in the energy balance routine. Surface temperature is found using the Root Brent method (Numerical Recipies). modifications: 02-29-00 Included variables needed to compute energy flux through the snow pack. The ground surface energy balance will now be a mixture of snow covered and bare ground, controlled by the snow cover fraction set in solve_snow.c KAC 04-Jun-04 Placed "ERROR" at beginning of screen dump in error_print_surf_energy_bal. TJB 21-Sep-04 Added ErrorString to store error messages from root_brent. TJB***************************************************************/{ extern veg_lib_struct *veg_lib; extern option_struct options; char FIRST_SOLN[1]; double T2; double Ts_old; double T1_old; double Tair; double ra; double atmos_density; double albedo; double emissivity; double kappa1; double kappa2; double Cs1; double Cs2; double D1; double D2; double delta_t; double Vapor; double max_moist; double bubble; double expt; double T1; double snow_depth; double snow_density; double Tsnow_surf; double snow_cover_fraction; double snow_flux; int VEG; double Tsurf; double U; double error; double Wdew[2]; double *T_node; double Tnew_node[MAX_NODES]; double *dz_node; double *kappa_node; double *Cs_node; double *moist_node; double *bubble_node; double *expt_node; double *max_moist_node; double *ice_node; double *alpha; double *beta; double *gamma; layer_data_struct layer[MAX_LAYERS]; double T_lower, T_upper; char ErrorString[MAXSTRING]; /************************************************** Correct Aerodynamic Resistance for Stability **************************************************/ ra = aero_resist[0]; U = wind[0]; if (U > 0.0) ra /= StabilityCorrection(ref_height, displacement, T0, air_temp, U, roughness); else ra = HUGE_RESIST; /************************************************** Compute Evaporation and Transpiration **************************************************/ if(iveg!=Nveg) { if(veg_lib[veg_class].LAI[dmy->month-1] > 0.0) VEG = TRUE; else VEG = FALSE; } else VEG = FALSE; Vapor = vapor_flux / (double)dt / 3600.; /************************************************** Set All Variables For Use **************************************************/ T2 = energy->T[Nnodes-1]; Ts_old = energy->T[0]; T1_old = energy->T[1]; Tair = air_temp; atmos_density = atmos->density[hour]; albedo = bare_albedo; emissivity = 1.; kappa1 = energy->kappa[0]; kappa2 = energy->kappa[1]; Cs1 = energy->Cs[0]; Cs2 = energy->Cs[1]; D1 = soil_con->depth[0];/* D2 = soil_con->depth[1]; */ D2 = soil_con->depth[0]; delta_t = (double)dt*3600.; max_moist = soil_con->max_moist[0]/(soil_con->depth[0]*1000.); bubble = soil_con->bubble[0]; expt = soil_con->expt[0]; snow_depth = snow->depth; snow_density = snow->density; Tsnow_surf = snow->surf_temp; snow_cover_fraction = snow->coverage; Wdew[WET] = veg_var_wet->Wdew; if(options.DIST_PRCP) Wdew[DRY] = veg_var_dry->Wdew; FIRST_SOLN[0] = TRUE; /************************************************************* Prepare soil node variables for finite difference solution *************************************************************/ if(!options.QUICK_FLUX) { bubble_node = soil_con->bubble_node; expt_node = soil_con->expt_node; max_moist_node = soil_con->max_moist_node; alpha = soil_con->alpha; beta = soil_con->beta; gamma = soil_con->gamma; moist_node = energy->moist; kappa_node = energy->kappa_node; Cs_node = energy->Cs_node; T_node = energy->T; dz_node = soil_con->dz_node; ice_node = energy->ice; } else { bubble_node = NULL; expt_node = NULL; max_moist_node = NULL; alpha = NULL; beta = NULL; gamma = NULL; moist_node = NULL; kappa_node = NULL; Cs_node = NULL; T_node = NULL; dz_node = NULL; ice_node = NULL; } /************************************************** Find Surface Temperature Using Root Brent Method **************************************************/ if(options.FULL_ENERGY) { /** Added for temporary backwards compatability **/ if(snow->swq > 0) { T_lower = 0.; T_upper = energy->T[0]-SURF_DT; } else { T_lower = 0.5*(energy->T[0]+air_temp)-SURF_DT; T_upper = 0.5*(energy->T[0]+air_temp)+SURF_DT; }#if QUICK_FS Tsurf = root_brent(T_upper, T_lower, ErrorString, func_surf_energy_bal, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist, ice0, max_moist, bubble, expt, snow_depth, snow_density, Tsnow_surf, snow_cover_fraction, surf_atten, U, displacement, roughness, ref_height, (double)soil_con->elevation, soil_con->b_infilt, soil_con->max_infil, (double)dt, atmos->vpd[hour], snow_energy, mu, rainfall, Wdew, &energy->grnd_flux, &T1, &energy->latent, &energy->sensible, &energy->deltaH, &energy->snow_flux, &energy->error, soil_con->depth, soil_con->Wcr, soil_con->Wpwp, soil_con->resid_moist, T_node, Tnew_node, dz_node, kappa_node, Cs_node, moist_node, bubble_node, expt_node, max_moist_node, ice_node, alpha, beta, gamma, soil_con->ufwc_table_layer[0], soil_con->ufwc_table_node, root, layer_wet, layer_dry, veg_var_wet, veg_var_dry, VEG, veg_class, dmy->month, Nnodes, FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);#else Tsurf = root_brent(T_upper, T_lower, ErrorString, func_surf_energy_bal, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist, ice0, max_moist, bubble, expt, snow_depth, snow_density, Tsnow_surf, snow_cover_fraction, surf_atten, U, displacement, roughness, ref_height, (double)soil_con->elevation, soil_con->b_infilt, soil_con->max_infil, (double)dt, atmos->vpd[hour], snow_energy, mu, rainfall, Wdew, &energy->grnd_flux, &T1, &energy->latent, &energy->sensible, &energy->deltaH, &energy->snow_flux, &energy->error, soil_con->depth, soil_con->Wcr, soil_con->Wpwp, soil_con->resid_moist, T_node, Tnew_node, dz_node, kappa_node, Cs_node, moist_node, bubble_node, expt_node, max_moist_node, ice_node, alpha, beta, gamma, root, layer_wet, layer_dry, veg_var_wet, veg_var_dry, VEG, veg_class, dmy->month, Nnodes, FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);#endif if(Tsurf <= -9998) { #if QUICK_FS error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist, ice0, max_moist, bubble, expt, snow_depth, snow_density, Tsnow_surf, snow_cover_fraction, surf_atten, U, displacement, roughness, ref_height, (double)soil_con->elevation, soil_con->b_infilt, soil_con->max_infil, (double)dt, atmos->vpd[hour], snow_energy, mu, rainfall, Wdew, &energy->grnd_flux, &T1, &energy->latent, &energy->sensible, &energy->deltaH, &energy->snow_flux, &energy->error, soil_con->depth, soil_con->Wcr, soil_con->Wpwp, soil_con->resid_moist, T_node, Tnew_node, dz_node, kappa_node, Cs_node, moist_node, bubble_node, expt_node, max_moist_node, ice_node, alpha, beta, gamma, soil_con->ufwc_table_layer[0], soil_con->ufwc_table_node, root, layer_wet, layer_dry, veg_var_wet, veg_var_dry, VEG, veg_class, dmy->month, iveg, Nnodes, FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE, ErrorString);#else error = error_calc_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist, ice0, max_moist, bubble, expt, snow_depth, snow_density, Tsnow_surf, snow_cover_fraction, surf_atten, U, displacement, roughness, ref_height, (double)soil_con->elevation, soil_con->b_infilt, soil_con->max_infil, (double)dt, atmos->vpd[hour], snow_energy, mu, rainfall, Wdew, &energy->grnd_flux, &T1, &energy->latent, &energy->sensible, &energy->deltaH, &energy->snow_flux, &energy->error, soil_con->depth, soil_con->Wcr, soil_con->Wpwp, soil_con->resid_moist, T_node, Tnew_node, dz_node, kappa_node, Cs_node, moist_node, bubble_node, expt_node, max_moist_node, ice_node, alpha, beta, gamma, root, layer_wet, layer_dry, veg_var_wet, veg_var_dry, VEG, veg_class, dmy->month, iveg, Nnodes, FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE, ErrorString);#endif } } else { /** Frozen soil model run with no surface energy balance **/ Tsurf = Tair; } /************************************************** Recalculate Energy Balance Terms for Final Surface Temperature **************************************************/#if QUICK_FS error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, emissivity, kappa1, kappa2, Cs1, Cs2, D1, D2, dp, delta_t, Le, Ls, Vapor, moist, ice0, max_moist, bubble, expt, snow_depth, snow_density, Tsnow_surf, snow_cover_fraction, surf_atten, U, displacement, roughness, ref_height, (double)soil_con->elevation, soil_con->b_infilt, soil_con->max_infil, (double)dt, atmos->vpd[hour], snow_energy, mu, rainfall, Wdew, &energy->grnd_flux, &T1, &energy->latent, &energy->sensible, &energy->deltaH, &energy->snow_flux, &energy->error, soil_con->depth, soil_con->Wcr, soil_con->Wpwp, soil_con->resid_moist, T_node, Tnew_node, dz_node, kappa_node, Cs_node, moist_node, bubble_node, expt_node, max_moist_node, ice_node, alpha, beta, gamma, soil_con->ufwc_table_layer[0], soil_con->ufwc_table_node, root, layer_wet, layer_dry, veg_var_wet, veg_var_dry, VEG, veg_class, dmy->month, Nnodes, FIRST_SOLN, snow->snow, soil_con->FS_ACTIVE);#else error = solve_surf_energy_bal(Tsurf, T2, Ts_old, T1_old, Tair, ra, atmos_density, shortwave, longwave, albedo, ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -