Commit a8733b7a authored by Beisman, James Joseph's avatar Beisman, James Joseph
Browse files

file renaming and organization

parent 79758d35
......@@ -9,7 +9,7 @@ It is not the same as the displa calculated in CanopyTemperature
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "clm_constants.h"
#include "landtype.h"
#include "frictionvelocity.h"
......
add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiation.cc SurfaceResistanceMod_functions.cc QSatMod_functions.cc CanopyTemperature.cc FrictionVelocityMod_functions.cc BareGroundFluxes.cc)
add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiation.cc SurfaceResistance.cc QSat.cc CanopyTemperature.cc FrictionVelocity.cc BareGroundFluxes.cc)
target_include_directories (physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
......@@ -6,7 +6,7 @@ functions derived from CanopyHydrologyMod.F90
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "clm_constants.h"
#include "landtype.h"
class CanopyHydrology {
......
......@@ -3,8 +3,8 @@
*/
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "QSatMod_functions.hh"
#include "clm_constants.h"
#include "qsat.h"
#include "landtype.h"
class CanopyTemperature {
......
......@@ -12,7 +12,7 @@ Vol. 11, 2628-2644.
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "clm_constants.h"
/* MoninObukIni()
DESCRIPTION: Initialization of the Monin-Obukhov length. The scheme is based on the work of
......@@ -71,7 +71,6 @@ void MoninObukIni(
double StabilityFunc1(double zeta)
{
double chik, chik2, retval;
chik2 = std::sqrt(1.0 - 16.0 * zeta);
chik = std::sqrt(chik2);
retval = 2.0 * std::log((1.0 + chik) * 0.5) + std::log((1.0 + chik2) * 0.5) - 2.0 * atan(chik) + ELM_PI * 0.5;
......@@ -82,7 +81,6 @@ double StabilityFunc1(double zeta)
double StabilityFunc2(double zeta)
{
double chik2, retval;
chik2 = std::sqrt(1.0 - 16.0 * zeta);
retval = 2.0 * std::log((1.0 + chik2) * 0.5);
return retval;
......@@ -108,7 +106,6 @@ void FrictionVelocityWind(
const double& um,
const double& obu,
const double& z0m,
double& ustar)
{
double zetam = 1.574; // transition point of flux-gradient relation (wind profile)
......@@ -146,7 +143,6 @@ void FrictionVelocityTemperature(
const double& displa,
const double& obu,
const double& z0h,
double& temp1)
{
double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
......@@ -189,7 +185,6 @@ void FrictionVelocityHumidity(
const double& z0h,
const double& z0q,
const double& temp1,
double& temp2)
{
double zetat = 0.465; // transition point of flux-gradient relation (temp. profile)
......
/*
*/
#include "clm_constants.h"
#include "qsat.h"
void QSat(
const double& T,
const double& p,
double& es,
double& esdT,
double& qs,
double& qsdT)
{
double td,vp,vp1,vp2,T_limit;
T_limit = T - tfrz;
if (T_limit > 100.0) { T_limit=100.0; }
if (T_limit < -75.0) { T_limit=-75.0; }
td = T_limit;
if (td >= 0.0) {
es = a0 + td *(a1 + td*(a2 + td*(a3 + td*(a4
+ td*(a5 + td*(a6 + td*(a7 + td*a8)))))));
esdT = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4
+ td*(b5 + td*(b6 + td*(b7 + td*b8)))))));
} else {
es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4
+ td*(c5 + td*(c6 + td*(c7 + td*c8)))))));
esdT = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4
+ td*(d5 + td*(d6 + td*(d7 + td*d8)))))));
}
es = es * 100.0; // pa
esdT = esdT * 100.0; // pa/K
vp = 1.0 / (p - 0.378 * es);
vp1 = 0.622 * vp;
vp2 = vp1 * vp;
qs = es * vp1; // kg/kg
qsdT = esdT * vp2 * p; // 1 / K
}
void QSat(const double& T, const double& p, double& es, double& esdT, double& qs, double& qsdT);
\ No newline at end of file
......@@ -10,7 +10,7 @@ SurfRadReflected
#include <cmath>
#include <assert.h>
#include "clm_constants.hh"
#include "clm_constants.h"
#include "landtype.h"
......
......@@ -2,7 +2,8 @@
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "clm_constants.h"
#include "landtype.h"
/* calc_soilevap_stress()
DESCRIPTION:
......@@ -12,8 +13,7 @@ This is the calc_beta_leepielke1992() function from CLM/ELM that gets called fro
code that gets called when use_vsfm == true is not currently included
INPUTS:
ltype [int] landunit type
lakpoi [bool] true => landunit is a lake point
Land [LandType] struct containing information about landtypet
frac_sno [double] fraction of ground covered by snow (0 to 1)
frac_h2osfc [double] fraction of ground covered by surface water (0 to 1)
watsat[nlevgrnd] [double] volumetric soil water at saturation (porosity)
......@@ -26,25 +26,24 @@ OUTPUTS:
soilbeta [double] factor that reduces ground evaporation
*/
void calc_soilevap_stress(
const int& ltype,
const bool& lakpoi,
const LandType& Land,
const double& frac_sno,
const double& frac_h2osfc,
const double* watsat,
const double* watfc,
const double* h2osoi_liq,
const double* h2osoi_ice,
const double* dz,
const double *watsat,
const double *watfc,
const double *h2osoi_liq,
const double *h2osoi_ice,
const double *dz,
double& soilbeta)
{
if (!lakpoi) {
if (!Land.lakpoi) {
// local variables
double fac, fac_fc, wx;
if (ltype != istwet && ltype != istice && ltype != istice_mec) {
if (ltype == istsoil || ltype == istcrop) {
if (Land.ltype != istwet && Land.ltype != istice && Land.ltype != istice_mec) {
if (Land.ltype == istsoil || Land.ltype == istcrop) {
wx = (h2osoi_liq[nlevsno] / denh2o + h2osoi_ice[nlevsno] / denice) / dz[nlevsno];
fac = std::min(1.0, wx / watsat[0]);
fac = std::max(fac, 0.01);
......@@ -58,11 +57,11 @@ void calc_soilevap_stress(
} else {
soilbeta = 1.0;
}
} else if (ltype == icol_road_perv) {
} else if (Land.ltype == icol_road_perv) {
soilbeta = 0.0;
} else if (ltype == icol_sunwall || ltype == icol_shadewall) {
} else if (Land.ltype == icol_sunwall || Land.ltype == icol_shadewall) {
soilbeta = 0.0;
} else if (ltype == icol_roof || ltype == icol_road_imperv) {
} else if (Land.ltype == icol_roof || Land.ltype == icol_road_imperv) {
soilbeta = 0.0;
}
} else {
......
......@@ -13,7 +13,7 @@ micro_sigma [double] microtopography pdf sigma (m)
#include <algorithm>
#include <cmath>
#include "clm_constants.hh"
#include "clm_constants.h"
void TopoSlopes(const double& topo_slope_in,
double& topo_slope_out)
......
/*
#ifndef QSAT_H_
#define QSAT_H_
void QSat(const double& T, const double& p, double& es, double& esdT, double& qs, double& qsdT);
*/
#include "clm_constants.hh"
void QSat(
const double& T,
const double& p,
double& es,
double& esdT,
double& qs,
double& qsdT)
{
// For water vapor (temperature range 0C-100C)
// For water vapor (temperature range 0C-100C)
const double a0 = 6.11213476;
const double a1 = 0.444007856;
const double a2 = 0.143064234e-01;
......@@ -54,31 +44,4 @@ void QSat(
const double d7 = 0.394116744e-13;
const double d8 = 0.498070196e-16;
double td,vp,vp1,vp2,T_limit;
T_limit = T - tfrz;
if (T_limit > 100.0) { T_limit=100.0; }
if (T_limit < -75.0) { T_limit=-75.0; }
td = T_limit;
if (td >= 0.0) {
es = a0 + td *(a1 + td*(a2 + td*(a3 + td*(a4
+ td*(a5 + td*(a6 + td*(a7 + td*a8)))))));
esdT = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4
+ td*(b5 + td*(b6 + td*(b7 + td*b8)))))));
} else {
es = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4
+ td*(c5 + td*(c6 + td*(c7 + td*c8)))))));
esdT = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4
+ td*(d5 + td*(d6 + td*(d7 + td*d8)))))));
}
es = es * 100.0; // pa
esdT = esdT * 100.0; // pa/K
vp = 1.0 / (p - 0.378 * es);
vp1 = 0.622 * vp;
vp2 = vp1 * vp;
qs = es * vp1; // kg/kg
qsdT = esdT * vp2 * p; // 1 / K
}
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment