From e0899ef348bea2f3c39cc079ca23ff05c108f485 Mon Sep 17 00:00:00 2001 From: ftomei Date: Sat, 9 Dec 2023 16:12:09 +0100 Subject: [PATCH 1/4] update snow --- crop/root.cpp | 11 ++-- crop/root.h | 2 +- snow/snow.cpp | 129 ++++++++++++++++++++++++---------------------- snow/snow.h | 4 +- snow/snowMaps.cpp | 8 +-- snow/snowMaps.h | 4 +- 6 files changed, 81 insertions(+), 77 deletions(-) diff --git a/crop/root.cpp b/crop/root.cpp index 93a0d6707..02ddc3f8d 100644 --- a/crop/root.cpp +++ b/crop/root.cpp @@ -493,7 +493,7 @@ namespace root } - bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, int nrLayers, + bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, unsigned int nrLayers, const std::vector &layerDepth, const std::vector &layerThickness) { // check soil @@ -506,7 +506,7 @@ namespace root // Initialize myCrop->roots.rootDensity.clear(); - myCrop->roots.rootDensity.resize(nrLayers); + myCrop->roots.rootDensity.resize(nrLayers); for (unsigned int i = 0; i < nrLayers; i++) { myCrop->roots.rootDensity[i] = 0.0; @@ -554,14 +554,13 @@ namespace root numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers); } - double maxLayerDepth = layerDepth[nrLayers-1] + layerThickness[nrLayers-1] * 0.5; int atom = 0; - double currentDepth = atom / 100.; // [m] + double currentDepth = double(atom) * 0.01; // [m] double rootDensitySum = 0.; while (currentDepth <= maxLayerDepth && atom < nrAtoms) { - for (int l = 0; l < nrLayers; l++) + for (unsigned int l = 0; l < nrLayers; l++) { double upperDepth = layerDepth[l] - layerThickness[l] * 0.5; double lowerDepth = layerDepth[l] + layerThickness[l] * 0.5; @@ -574,7 +573,7 @@ namespace root } atom++; - currentDepth = atom / 100.; // [m] + currentDepth = double(atom) * 0.01; // [m] } if (rootDensitySum <= EPSILON) diff --git a/crop/root.h b/crop/root.h index 019d9b295..19caf0ee9 100644 --- a/crop/root.h +++ b/crop/root.h @@ -54,7 +54,7 @@ double getRootLengthDD(const Crit3DRoot &myRoot, double currentDD, double emergenceDD); bool computeRootDensity(Crit3DCrop* myCrop, const std::vector &soilLayers); - bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, int nrLayers, + bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, unsigned int nrLayers, const std::vector &layerDepth, const std::vector &layerThickness); int getNrAtoms(const std::vector &soilLayers, double &minThickness, std::vector &atoms); diff --git a/snow/snow.cpp b/snow/snow.cpp index 85e050705..5938279ca 100644 --- a/snow/snow.cpp +++ b/snow/snow.cpp @@ -43,7 +43,7 @@ Crit3DSnowParameters::Crit3DSnowParameters() void Crit3DSnowParameters::initialize() { // default values - snowSkinThickness = 0.02; /*!< [m] */ // LC: VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice) + skinThickness = 0.02; /*!< [m] */ // LC: VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice) soilAlbedo = 0.2; /*!< [-] bare soil - 20% */ snowVegetationHeight = 1; snowWaterHoldingCapacity = 0.05; @@ -151,8 +151,8 @@ void Crit3DSnow::computeSnowBrooksModel() double prevIceContent, prevLWaterContent; /*!< [mm] */ double currentRatio; - double AirActualVapDensity; /*!< [kg m-3] */ - double WaterActualVapDensity; /*!< [kg m-3] */ + double airActualVapDensity; /*!< [kg m-3] */ + double waterActualVapDensity; /*!< [kg m-3] */ double longWaveAtmEmissivity; /*!< [-] */ double albedo; /*!< [-] */ @@ -215,10 +215,10 @@ void Crit3DSnow::computeSnowBrooksModel() vegetationShadowing = std::max(std::min(heightVegetation / maxVegetationHeight, 1.), 0.); solarRadTot = _globalRadiation - _beamRadiation * vegetationShadowing; - double prevSurfacetemp = _snowSurfaceTemp; double previousSWE = _snowWaterEquivalent; double prevInternalEnergy = _internalEnergy; double prevSurfaceEnergy = _surfaceEnergy; + double prevSurfaceTemp = _snowSurfaceTemp; /*-------------------------------------------------------------------- // controlli di coerenza per eventuali modifiche manuali su mappa SWE @@ -238,8 +238,8 @@ void Crit3DSnow::computeSnowBrooksModel() // stato fittizio:: neve recente prossima alla fusione, con una settimana di età _ageOfSnow = 7; - prevSurfacetemp = std::min(prevSurfacetemp, -0.1); - prevSurfaceEnergy = std::min(prevSurfaceEnergy, -0.1); + prevSurfaceTemp = std::min(prevSurfaceTemp, 0.); + prevSurfaceEnergy = std::min(prevSurfaceEnergy, 0.); } /*! check on sum */ @@ -268,21 +268,22 @@ void Crit3DSnow::computeSnowBrooksModel() // pag. 52 (3.20) // source: Jensen et al. (1990) and Tetens (1930) // saturated vapor density - AirActualVapDensity = double(exp((16.78 * dewPoint - 116.9) / (dewPoint + 237.3)) + airActualVapDensity = double(exp((16.78 * dewPoint - 116.9) / (dewPoint + 237.3)) / ((ZEROCELSIUS + dewPoint) * THERMO_WATER_VAPOR) ); - // over water - WaterActualVapDensity = double( exp((16.78 * prevSurfacetemp - 116.9) / (prevSurfacetemp + 237.3)) - / ((ZEROCELSIUS + prevSurfacetemp) * THERMO_WATER_VAPOR) ); + // over snow + waterActualVapDensity = double( exp((16.78 * prevSurfaceTemp - 116.9) / (prevSurfaceTemp + 237.3)) + / ((ZEROCELSIUS + prevSurfaceTemp) * THERMO_WATER_VAPOR) ); // over ice // LC: non trovo riferimenti a questa formula /* - if (prevInternalEnergy <= 0) + if (previousSWE > SNOW_MINIMUM_HEIGHT) { WaterActualVapDensity *= exp(MH2O * LATENT_HEAT_FUSION * prevSurfacetemp * 1000 - / (R_GAS * pow(prevSurfacetemp + ZEROCELSIUS, 2))); - }*/ + / (R_GAS * pow(prevSurfacetemp + ZEROCELSIUS, 2)) ); + } + */ /*! * \brief Atmospheric Emissivity Calculations for Longwave Radiation @@ -304,12 +305,13 @@ void Crit3DSnow::computeSnowBrooksModel() /*! \brief Incoming Energy Fluxes */ // pag. 52 (3.22) considerando i 2 contributi invece che solo uno - QPrecipW = (HEAT_CAPACITY_WATER / 1000.) * (_precRain / 1000.) * (std::max(0., _airT) - prevSurfacetemp); - QPrecipS = (HEAT_CAPACITY_SNOW / 1000.) * (_precSnow / 1000.) * (std::min(0., _airT) - prevSurfacetemp); + QPrecipW = (HEAT_CAPACITY_WATER / 1000.) * (_precRain / 1000.) * (std::max(0., _airT) - prevSurfaceTemp); + QPrecipS = (HEAT_CAPACITY_SNOW / 1000.) * (_precSnow / 1000.) * (std::min(0., _airT) - prevSurfaceTemp); QPrecip = QPrecipW + QPrecipS; // temperatura dell'acqua: almeno 1 grado - QWaterHeat = (HEAT_CAPACITY_WATER / 1000.) * (_surfaceWaterContent / 1000.) * (std::max(1., (prevSurfacetemp + _airT) / 2.) - prevSurfacetemp); + QWaterHeat = (HEAT_CAPACITY_WATER / 1000.) * (_surfaceWaterContent / 1000.) + * (std::max(1., (prevSurfaceTemp + _airT) / 2.) - prevSurfaceTemp); // energia acqua libera QWaterKinetic = 0; @@ -346,15 +348,15 @@ void Crit3DSnow::computeSnowBrooksModel() // pag. 50 (3.15) QLongWave = double(STEFAN_BOLTZMANN * 3.6 * (longWaveAtmEmissivity * pow(_airT + ZEROCELSIUS, 4.0) - - surfaceEmissivity * pow (prevSurfacetemp + ZEROCELSIUS, 4.0))); + - surfaceEmissivity * pow (prevSurfaceTemp + ZEROCELSIUS, 4.0))); // sensible heat pag. 50 (3.17) - QTempGradient = 3600. * (HEAT_CAPACITY_AIR / 1000.) * (_airT - prevSurfacetemp) / aerodynamicResistance; + QTempGradient = 3600. * (HEAT_CAPACITY_AIR / 1000.) * (_airT - prevSurfaceTemp) / aerodynamicResistance; // latent heat pag. 51 (3.19) // FT tolta WATER_DENSITY dall'eq. (non corrispondevano le unità di misura) QVaporGradient = 3600. * (LATENT_HEAT_VAPORIZATION + LATENT_HEAT_FUSION) - * (AirActualVapDensity - WaterActualVapDensity) / aerodynamicResistance; + * (airActualVapDensity - waterActualVapDensity) / aerodynamicResistance; // TODO serve formula diversa quando non c'è neve if (previousSWE < SNOW_MINIMUM_HEIGHT) @@ -362,64 +364,68 @@ void Crit3DSnow::computeSnowBrooksModel() QVaporGradient *= 0.4; } - /*! \brief Energy Balance */ + /*! Energy Balance */ QTotal = QSolar + QPrecip + QLongWave + QTempGradient + QVaporGradient + QWaterHeat + QWaterKinetic; _sensibleHeat = QTempGradient; _latentHeat = QVaporGradient; - /*! \brief Condensation/Evaporation */ - sublimation = 0; + /*! Condensation (positive) or evaporation (negative) */ + sublimation = 0; // [mm] + _evaporation = 0; // [mm] if (previousSWE > SNOW_MINIMUM_HEIGHT) { - // pag. 51 (3.21) - sublimation = QVaporGradient / (WATER_DENSITY * (LATENT_HEAT_FUSION + LATENT_HEAT_VAPORIZATION)); // [m] - sublimation *= 1000.; // [mm] + // pag. 51 (3.21) [mm] + sublimation = QVaporGradient / (LATENT_HEAT_FUSION + LATENT_HEAT_VAPORIZATION); if (sublimation < 0) { - // controllo aggiunto: può evaporare al massimo la neve presente + /*! Evaporation [mm] + * controllo aggiunto (può evaporare solo la neve presente) + */ sublimation = -std::min(fabs(sublimation), previousSWE + _precSnow); + _evaporation = -sublimation; } } - /*! sign of evaporation is negative */ - if (sublimation < 0) - _evaporation = -sublimation; - else - _evaporation = 0; - /*! refreeze or melt */ double freeze_melt = 0; // [mm] freeze (positive) or melt (negative) - // pag.53 (3.25) /*! net amount of liquid water that freezes (heat added to the snow pack) - * and ice that melts (heat removed from the snow pack) */ - double w = (prevInternalEnergy + QTotal) / (LATENT_HEAT_FUSION * WATER_DENSITY); // [m] + * and ice that melts (heat removed from the snow pack) + * pag.53 (3.25) + */ + double w = (prevInternalEnergy + QTotal) / (LATENT_HEAT_FUSION * WATER_DENSITY); // [m] if (w < 0) { - // controllo aggiunto - if (prevSurfacetemp <= 0) + /*! freeze + * add check on surface temperatures + */ + if (prevSurfaceTemp <= 0) { - /*! refreeze */ freeze_melt = std::min(prevLWaterContent + _precRain, -w * 1000.); // [mm] } } - else + else if (w > 0) { /*! melt */ freeze_melt = -std::min(prevIceContent + _precSnow + sublimation, w * 1000.); // [mm] } - // pag.53 (3.23) modificata (errore nel denominatore) + /*! Snowmelt or refreeze [mm] - source/sink for Criteria3D */ + _snowMelt = -freeze_melt; + + /*! latent heat exchange in the snow pack [kJ m-2] + * pag.53 (3.23) modificata (errore nel denominatore) + */ double Qr = (freeze_melt / 1000.) * LATENT_HEAT_FUSION * WATER_DENSITY; - /*! Internal energy */ + /*! Internal energy [kJ m-2] */ _internalEnergy = prevInternalEnergy + QTotal + Qr; /*! Snow Pack Mass */ /*! Ice content */ - if (_internalEnergy >= EPSILON) + if (_internalEnergy > EPSILON) { _iceContent = 0; } @@ -433,7 +439,7 @@ void Crit3DSnow::computeSnowBrooksModel() / (1 - snowParameters.snowWaterHoldingCapacity); // [%] /*! Liquid water content */ - if (_internalEnergy >= EPSILON) + if (_internalEnergy > EPSILON) { _liquidWaterContent = 0; } @@ -450,37 +456,36 @@ void Crit3DSnow::computeSnowBrooksModel() if (_snowWaterEquivalent < EPSILON) _ageOfSnow = NODATA; else if (_ageOfSnow == NODATA || _precSnow > EPSILON) - _ageOfSnow = ONE_HOUR; + _ageOfSnow = 0; else _ageOfSnow += ONE_HOUR; - /*! Snowmelt (or refreeze) - source/sink for Criteria3D */ - _snowMelt = previousSWE + _precSnow + sublimation - _snowWaterEquivalent; - - /*! Snow surface energy */ + /*! Snow surface energy [kJ m-2] */ + double swe_m = _snowWaterEquivalent / 1000.; // [m] if (_snowWaterEquivalent > EPSILON) { - if (fabs(_internalEnergy) < EPSILON) - _surfaceEnergy = 0.; - else + // snow + double ratio = (std::min(swe_m, snowParameters.skinThickness) / SNOW_DAMPING_DEPTH) + + (std::max(snowParameters.skinThickness - swe_m, 0.) / SOIL_DAMPING_DEPTH); + _surfaceEnergy = prevSurfaceEnergy + (QTotal + Qr) * ratio; + // check + if (_snowWaterEquivalent >= SNOW_MINIMUM_HEIGHT && _internalEnergy <= EPSILON && _surfaceEnergy > 0) { - _surfaceEnergy = std::min(0., prevSurfaceEnergy + (QTotal + Qr) - * (std::min(_snowWaterEquivalent / 1000., snowParameters.snowSkinThickness) / SNOW_DAMPING_DEPTH - + std::max(snowParameters.snowSkinThickness - (_snowWaterEquivalent / 1000.), 0.) / SOIL_DAMPING_DEPTH)); + _surfaceEnergy *= 0.5; } } else { // soil - _surfaceEnergy = prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.snowSkinThickness / SOIL_DAMPING_DEPTH); + _surfaceEnergy = prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SOIL_DAMPING_DEPTH); } - _snowSurfaceTemp = _surfaceEnergy / ((WATER_DENSITY * SNOW_SPECIFIC_HEAT * std::min(snowParameters.snowSkinThickness, _snowWaterEquivalent / 1000.)) - + (DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * std::max(0., snowParameters.snowSkinThickness - _snowWaterEquivalent / 1000.))); + _snowSurfaceTemp = _surfaceEnergy / ((WATER_DENSITY * SNOW_SPECIFIC_HEAT * std::min(snowParameters.skinThickness, swe_m)) + + (DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * std::max(0., snowParameters.skinThickness - swe_m))); - // controllo aggiunto per problemi numerici - _snowSurfaceTemp = std::min(_snowSurfaceTemp, prevSurfacetemp + 10); - _snowSurfaceTemp = std::max(_snowSurfaceTemp, prevSurfacetemp - 10); + // controlli aggiunti per problemi numerici + _snowSurfaceTemp = std::min(_snowSurfaceTemp, prevSurfaceTemp + 10); + _snowSurfaceTemp = std::max(_snowSurfaceTemp, prevSurfaceTemp - 10); } @@ -622,8 +627,8 @@ double computeInternalEnergy(double initSoilPackTemp,int bulkDensity, double ini } -double computeSurfaceEnergy(double initSnowSurfaceTemp, double snowSkinThickness) +double computeSurfaceEnergy(double initSnowSurfaceTemp, double skinThickness) { - return initSnowSurfaceTemp * WATER_DENSITY * SNOW_SPECIFIC_HEAT * snowSkinThickness; + return initSnowSurfaceTemp * WATER_DENSITY * SNOW_SPECIFIC_HEAT * skinThickness; } diff --git a/snow/snow.h b/snow/snow.h index 85346dc58..c6a929506 100644 --- a/snow/snow.h +++ b/snow/snow.h @@ -29,7 +29,7 @@ class Crit3DSnowParameters { public: - double snowSkinThickness; /*!< [m] */ + double skinThickness; /*!< [m] */ double soilAlbedo; /*!< [-] bare soil */ double snowVegetationHeight; /*!< [m] height of vegetation */ double snowWaterHoldingCapacity; /*!< [-] percentuale di acqua libera che il manto nevoso può trattenere */ @@ -113,7 +113,7 @@ double aerodynamicResistanceCampbell77(bool isSnow , double zRefWind, double windSpeed, double vegetativeHeight); double computeInternalEnergy(double initSoilPackTemp,int bulkDensity, double initSWE); - double computeSurfaceEnergy(double initSnowSurfaceTemp, double snowSkinThickness); + double computeSurfaceEnergy(double initSnowSurfaceTemp, double skinThickness); #endif // SNOW_H diff --git a/snow/snowMaps.cpp b/snow/snowMaps.cpp index 19405cbb2..cf9c17c31 100644 --- a/snow/snowMaps.cpp +++ b/snow/snowMaps.cpp @@ -77,7 +77,7 @@ void Crit3DSnowMaps::clear() } -void Crit3DSnowMaps::initialize(const gis::Crit3DRasterGrid &dtm, double snowSkinThickness) +void Crit3DSnowMaps::initialize(const gis::Crit3DRasterGrid &dtm, double skinThickness) { _snowFallMap->initializeGrid(dtm); _snowMeltMap->initializeGrid(dtm); @@ -99,7 +99,7 @@ void Crit3DSnowMaps::initialize(const gis::Crit3DRasterGrid &dtm, double snowSki // initialize with zero values _snowWaterEquivalentMap->setConstantValueWithBase(0, dtm); - resetSnowModel(snowSkinThickness); + resetSnowModel(skinThickness); isInitialized = true; } @@ -151,7 +151,7 @@ void Crit3DSnowMaps::setPoint(Crit3DSnow &snowPoint, int row, int col) } -void Crit3DSnowMaps::resetSnowModel(double snowSkinThickness) +void Crit3DSnowMaps::resetSnowModel(double skinThickness) { float initSWE; /*!< [mm] */ int surfaceBulkDensity; /*!< [kg/m^3] */ @@ -172,7 +172,7 @@ void Crit3DSnowMaps::resetSnowModel(double snowSkinThickness) _snowSurfaceTempMap->value[row][col] = float(_initSnowSurfaceTemp); - _surfaceEnergyMap->value[row][col] = float(computeSurfaceEnergy(_initSnowSurfaceTemp, snowSkinThickness)); + _surfaceEnergyMap->value[row][col] = float(computeSurfaceEnergy(_initSnowSurfaceTemp, skinThickness)); _internalEnergyMap->value[row][col] = float(computeInternalEnergy(_initSoilPackTemp, surfaceBulkDensity, initSWE/1000.)); diff --git a/snow/snowMaps.h b/snow/snowMaps.h index a032cdcd5..5cf20baa3 100644 --- a/snow/snowMaps.h +++ b/snow/snowMaps.h @@ -17,8 +17,8 @@ ~Crit3DSnowMaps(); void clear(); - void initialize(const gis::Crit3DRasterGrid &dtm, double snowSkinThickness); - void resetSnowModel(double snowSkinThickness); + void initialize(const gis::Crit3DRasterGrid &dtm, double skinThickness); + void resetSnowModel(double skinThickness); void updateMap(Crit3DSnow &snowPoint, int row, int col); void setPoint(Crit3DSnow &snowPoint, int row, int col); From dcf61f67ee2f934aaa1ee7d10d0215aa3385f3ce Mon Sep 17 00:00:00 2001 From: Fausto Tomei Date: Sat, 9 Dec 2023 18:53:57 +0100 Subject: [PATCH 2/4] fix date 3D --- project/project.cpp | 7 ++++--- snow/snow.cpp | 8 ++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/project/project.cpp b/project/project.cpp index deed22571..958122551 100644 --- a/project/project.cpp +++ b/project/project.cpp @@ -851,11 +851,12 @@ Crit3DTime Project::getCrit3DCurrentTime() QDateTime Project::getCurrentTime() { - QDateTime myTime; - myTime.setDate(this->currentDate); - return myTime.addSecs(this->currentHour * HOUR_SECONDS); + QDateTime myDateTime; + myDateTime.setDate(this->currentDate); + return myDateTime.addSecs(this->currentHour * HOUR_SECONDS); } + void Project::getMeteoPointsRange(float& minimum, float& maximum, bool useNotActivePoints) { minimum = NODATA; diff --git a/snow/snow.cpp b/snow/snow.cpp index 5938279ca..d6c1a4353 100644 --- a/snow/snow.cpp +++ b/snow/snow.cpp @@ -43,8 +43,8 @@ Crit3DSnowParameters::Crit3DSnowParameters() void Crit3DSnowParameters::initialize() { // default values - skinThickness = 0.02; /*!< [m] */ // LC: VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice) - soilAlbedo = 0.2; /*!< [-] bare soil - 20% */ + skinThickness = 0.02; /*!< [m] */ // LC: VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice) + soilAlbedo = 0.2; /*!< [-] bare soil */ snowVegetationHeight = 1; snowWaterHoldingCapacity = 0.05; snowMaxWaterContent = 0.1; @@ -205,8 +205,8 @@ void Crit3DSnow::computeSnowBrooksModel() else cloudCover = 0.1; - // ombreggiamento per vegetazione (4m sopra manto nevoso: ombreggiamento completo) - // TODO improve - use LAI when available + // vegetation shadowing (4 m sopra manto nevoso: ombreggiamento completo) + // TODO: use LAI when available double maxSnowDensity = 10; // 1 mm snow = 1 cm water double maxVegetationHeight = 4; // [m] double vegetationShadowing; // [-] From 01ac71434933d159261e02e6455b72019dda75e5 Mon Sep 17 00:00:00 2001 From: Fausto Tomei Date: Sun, 10 Dec 2023 20:03:08 +0100 Subject: [PATCH 3/4] update snow --- mathFunctions/commonConstants.h | 2 +- snow/snow.cpp | 86 +++++++++++++++------------------ snow/snow.h | 40 +++++++-------- 3 files changed, 59 insertions(+), 69 deletions(-) diff --git a/mathFunctions/commonConstants.h b/mathFunctions/commonConstants.h index a2d4d260c..c2dbb3db3 100644 --- a/mathFunctions/commonConstants.h +++ b/mathFunctions/commonConstants.h @@ -175,7 +175,7 @@ // [W m-2 K-4] Stefan-Boltzmann constant #define STEFAN_BOLTZMANN 5.670373E-8 - // [] Von Kármán constant + // [-] Von Kármán constant #define VON_KARMAN_CONST 0.41 // [J kg-1 K-1] specific heat at constant pressure #define CP 1013. diff --git a/snow/snow.cpp b/snow/snow.cpp index d6c1a4353..5f6e5c497 100644 --- a/snow/snow.cpp +++ b/snow/snow.cpp @@ -31,9 +31,6 @@ #include "snow.h" -// Warning: commonConstant unit is [J m-3 K-1] -// check division by 1000 to transform in [kJ] - Crit3DSnowParameters::Crit3DSnowParameters() { initialize(); @@ -87,7 +84,7 @@ void Crit3DSnow::initialize() _liquidWaterContent = NODATA; _internalEnergy = NODATA; _surfaceEnergy = NODATA; - _snowSurfaceTemp = NODATA; + _surfaceTemp = NODATA; _ageOfSnow = NODATA; } @@ -114,7 +111,7 @@ bool Crit3DSnow::checkValidPoint() || int(_globalRadiation) == int(NODATA) || int(_beamRadiation) == int(NODATA) || int(_snowWaterEquivalent) == int(NODATA) - || int(_snowSurfaceTemp) == int(NODATA) ) + || int(_surfaceTemp) == int(NODATA) ) { return false; } @@ -184,7 +181,7 @@ void Crit3DSnow::computeSnowBrooksModel() _liquidWaterContent = NODATA; _snowWaterEquivalent = NODATA; _surfaceEnergy = NODATA; - _snowSurfaceTemp = NODATA; + _surfaceTemp = NODATA; _ageOfSnow = NODATA; _precSnow = NODATA; @@ -218,7 +215,7 @@ void Crit3DSnow::computeSnowBrooksModel() double previousSWE = _snowWaterEquivalent; double prevInternalEnergy = _internalEnergy; double prevSurfaceEnergy = _surfaceEnergy; - double prevSurfaceTemp = _snowSurfaceTemp; + double prevSurfaceTemp = _surfaceTemp; /*-------------------------------------------------------------------- // controlli di coerenza per eventuali modifiche manuali su mappa SWE @@ -278,10 +275,10 @@ void Crit3DSnow::computeSnowBrooksModel() // over ice // LC: non trovo riferimenti a questa formula /* - if (previousSWE > SNOW_MINIMUM_HEIGHT) + if (previousSWE > EPSILON) { - WaterActualVapDensity *= exp(MH2O * LATENT_HEAT_FUSION * prevSurfacetemp * 1000 - / (R_GAS * pow(prevSurfacetemp + ZEROCELSIUS, 2)) ); + waterActualVapDensity *= exp(MH2O * LATENT_HEAT_FUSION * prevSurfaceTemp + / (R_GAS * 0.001 * pow(prevSurfaceTemp + ZEROCELSIUS, 2)) ); } */ @@ -359,7 +356,7 @@ void Crit3DSnow::computeSnowBrooksModel() * (airActualVapDensity - waterActualVapDensity) / aerodynamicResistance; // TODO serve formula diversa quando non c'è neve - if (previousSWE < SNOW_MINIMUM_HEIGHT) + if (previousSWE < EPSILON) { QVaporGradient *= 0.4; } @@ -372,7 +369,7 @@ void Crit3DSnow::computeSnowBrooksModel() /*! Condensation (positive) or evaporation (negative) */ sublimation = 0; // [mm] _evaporation = 0; // [mm] - if (previousSWE > SNOW_MINIMUM_HEIGHT) + if (previousSWE > EPSILON) { // pag. 51 (3.21) [mm] sublimation = QVaporGradient / (LATENT_HEAT_FUSION + LATENT_HEAT_VAPORIZATION); @@ -436,7 +433,7 @@ void Crit3DSnow::computeSnowBrooksModel() } double waterHoldingCapacity = snowParameters.snowWaterHoldingCapacity - / (1 - snowParameters.snowWaterHoldingCapacity); // [%] + / (1 - snowParameters.snowWaterHoldingCapacity); // [%] /*! Liquid water content */ if (_internalEnergy > EPSILON) @@ -452,40 +449,33 @@ void Crit3DSnow::computeSnowBrooksModel() /*! Snow water equivalent */ _snowWaterEquivalent = _iceContent + _liquidWaterContent; - /*! _ageOfSnow */ - if (_snowWaterEquivalent < EPSILON) - _ageOfSnow = NODATA; - else if (_ageOfSnow == NODATA || _precSnow > EPSILON) - _ageOfSnow = 0; - else - _ageOfSnow += ONE_HOUR; + /*! surface energy [kJ m-2] and surface temeprature [°C] */ + // all snow + double surfaceEnergySnow = std::min(0., prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SNOW_DAMPING_DEPTH) ); + double surfaceTempSnow = surfaceEnergySnow / (WATER_DENSITY * SNOW_SPECIFIC_HEAT * snowParameters.skinThickness); + + // all soil + double surfaceEnergySoil = prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SOIL_DAMPING_DEPTH); + double surfaceTempSoil = surfaceEnergySoil / (DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * snowParameters.skinThickness); + + double snowDepthRatio = 5.; + double snowFraction = std::min(_snowWaterEquivalent * snowDepthRatio / 1000., snowParameters.skinThickness) / snowParameters.skinThickness; - /*! Snow surface energy [kJ m-2] */ - double swe_m = _snowWaterEquivalent / 1000.; // [m] - if (_snowWaterEquivalent > EPSILON) + _surfaceEnergy = (surfaceEnergySnow * snowFraction) + surfaceEnergySoil * (1 - snowFraction); + _surfaceTemp = (surfaceTempSnow * snowFraction) + surfaceTempSoil * (1 - snowFraction); + + /*! _ageOfSnow [days] */ + if (_snowWaterEquivalent < EPSILON) { - // snow - double ratio = (std::min(swe_m, snowParameters.skinThickness) / SNOW_DAMPING_DEPTH) - + (std::max(snowParameters.skinThickness - swe_m, 0.) / SOIL_DAMPING_DEPTH); - _surfaceEnergy = prevSurfaceEnergy + (QTotal + Qr) * ratio; - // check - if (_snowWaterEquivalent >= SNOW_MINIMUM_HEIGHT && _internalEnergy <= EPSILON && _surfaceEnergy > 0) - { - _surfaceEnergy *= 0.5; - } + _ageOfSnow = NODATA; } else { - // soil - _surfaceEnergy = prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SOIL_DAMPING_DEPTH); + if (_ageOfSnow == NODATA || _precSnow > EPSILON) + _ageOfSnow = 0; + else + _ageOfSnow += ONE_HOUR; } - - _snowSurfaceTemp = _surfaceEnergy / ((WATER_DENSITY * SNOW_SPECIFIC_HEAT * std::min(snowParameters.skinThickness, swe_m)) - + (DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * std::max(0., snowParameters.skinThickness - swe_m))); - - // controlli aggiunti per problemi numerici - _snowSurfaceTemp = std::min(_snowSurfaceTemp, prevSurfaceTemp + 10); - _snowSurfaceTemp = std::max(_snowSurfaceTemp, prevSurfaceTemp - 10); } @@ -560,12 +550,12 @@ void Crit3DSnow::setSurfaceEnergy(float value) double Crit3DSnow::getSnowSurfaceTemp() { - return _snowSurfaceTemp; + return _surfaceTemp; } void Crit3DSnow::setSnowSurfaceTemp(float value) { - _snowSurfaceTemp = double(value); + _surfaceTemp = double(value); } double Crit3DSnow::getAgeOfSnow() @@ -581,10 +571,10 @@ void Crit3DSnow::setAgeOfSnow(float value) /*! * \brief Computes aerodynamic Resistance (resistance to heat transfer) - * Brooks pag.51 (3.18) - * \param zRefWind [m] heights of windspeed measurements - * \param windSpeed [m s-1] wind speed - * \param vegetativeHeight [m] height of the vegetation + * \param zRefWind [m] heights of windspeed measurements + * \param windSpeed [m s-1] wind speed + * \param vegetativeHeight [m] height of the vegetation + * \ref Brooks pag.51 (3.18) * \return aerodynamic Resistance [s m-1] */ double aerodynamicResistanceCampbell77(bool isSnow , double zRefWind, double windSpeed, double vegetativeHeight) @@ -593,7 +583,7 @@ double aerodynamicResistanceCampbell77(bool isSnow , double zRefWind, double win double momentumRoughness; /*! [m] momentum roughness parameter (for snow = 0.001m, for vegetative cover zm = 0.13 times the height of the vegetation) */ double zRefTemp = 2; /*! [m] heights of temperature measurements */ - /*! check on wind speed [m/s] */ + /*! check on wind speed [m s-1] */ windSpeed = std::max(windSpeed, 0.05); windSpeed = std::min(windSpeed, 10.); diff --git a/snow/snow.h b/snow/snow.h index c6a929506..8ec4183eb 100644 --- a/snow/snow.h +++ b/snow/snow.h @@ -15,27 +15,27 @@ /*! * heat of fusion for ice at 0 °C */ - #define LATENT_HEAT_FUSION 335. /*!< [kJ kg-1] */ - #define LATENT_HEAT_VAPORIZATION 2500. /*!< [kJ kg-1] */ + #define LATENT_HEAT_FUSION 335. /*!< [kJ kg-1] */ + #define LATENT_HEAT_VAPORIZATION 2500. /*!< [kJ kg-1] */ #define SNOW_SPECIFIC_HEAT 2.1 /*!< [KJ kg-1 °C-1] */ #define SOIL_SPECIFIC_HEAT 1.4 /*!< [KJ kg-1 °C-1] wet soil */ #define DEFAULT_BULK_DENSITY 1350 /*!< [kg m-3] */ #define SOIL_DAMPING_DEPTH 0.3 /*!< [m] */ #define SNOW_DAMPING_DEPTH 0.05 /*!< [m] */ - #define SNOW_MINIMUM_HEIGHT 2. /*!< [mm] */ + #define SNOW_MINIMUM_HEIGHT 1. /*!< [mm] */ class Crit3DSnowParameters { public: - double skinThickness; /*!< [m] */ - double soilAlbedo; /*!< [-] bare soil */ - double snowVegetationHeight; /*!< [m] height of vegetation */ - double snowWaterHoldingCapacity; /*!< [-] percentuale di acqua libera che il manto nevoso può trattenere */ - double snowMaxWaterContent; /*!< [m] acqua libera (torrenti, laghetti) */ - double tempMaxWithSnow; /*!< [°C] */ - double tempMinWithRain; /*!< [°C] */ + double skinThickness; /*!< [m] */ + double soilAlbedo; /*!< [-] bare soil */ + double snowVegetationHeight; /*!< [m] height of vegetation */ + double snowWaterHoldingCapacity; /*!< [-] percentuale di acqua libera che il manto nevoso può trattenere */ + double snowMaxWaterContent; /*!< [m] acqua libera (torrenti, laghetti) */ + double tempMaxWithSnow; /*!< [°C] */ + double tempMinWithRain; /*!< [°C] */ Crit3DSnowParameters(); @@ -82,15 +82,15 @@ private: // input - double _clearSkyTransmissivity; /*!< [-] */ - double _transmissivity; /*!< [-] */ - double _globalRadiation; /*!< [W m-2] */ - double _beamRadiation; /*!< [W m-2] */ - double _prec; /*!< [mm] */ - double _airT; /*!< [°C] */ - double _airRH; /*!< [%] */ - double _windInt; /*!< [m/s] */ - double _surfaceWaterContent; /*!< [mm] */ + double _clearSkyTransmissivity; /*!< [-] */ + double _transmissivity; /*!< [-] */ + double _globalRadiation; /*!< [W m-2] */ + double _beamRadiation; /*!< [W m-2] */ + double _prec; /*!< [mm] */ + double _airT; /*!< [°C] */ + double _airRH; /*!< [%] */ + double _windInt; /*!< [m/s] */ + double _surfaceWaterContent; /*!< [mm] */ // output double _evaporation; /*!< [mm] */ @@ -106,7 +106,7 @@ double _liquidWaterContent; /*!< [mm] */ double _internalEnergy; /*!< [kJ m-2] */ double _surfaceEnergy; /*!< [kJ m-2] */ - double _snowSurfaceTemp; /*!< [°C] */ + double _surfaceTemp; /*!< [°C] */ double _ageOfSnow; /*!< [days] */ }; From aefa16d967e09a4258ad64ece18e94eab5443e95 Mon Sep 17 00:00:00 2001 From: ftomei Date: Mon, 11 Dec 2023 16:02:34 +0100 Subject: [PATCH 4/4] update netcdf --- netcdfHandler/netcdfHandler.cpp | 17 +++++++++++++++++ netcdfHandler/netcdfHandler.h | 2 ++ 2 files changed, 19 insertions(+) diff --git a/netcdfHandler/netcdfHandler.cpp b/netcdfHandler/netcdfHandler.cpp index e0e409669..9a6331bbc 100644 --- a/netcdfHandler/netcdfHandler.cpp +++ b/netcdfHandler/netcdfHandler.cpp @@ -149,6 +149,8 @@ void NetCDFHandler::clear() variables.clear(); metadata.str(""); timeUnit = ""; + + missingValue = NODATA; } @@ -525,11 +527,25 @@ bool NetCDFHandler::readProperties(string fileName) else if (ncTypeId == NC_INT) { nc_get_att(ncId, v, attrName, &valueInt); + + // no data + if (lowerCase(string(attrName)) == "missing_value" || lowerCase(string(attrName)) == "nodata") + { + missingValue = double(valueInt); + } + metadata << attrName << " = " << valueInt << endl; } else if (ncTypeId == NC_DOUBLE) { nc_get_att(ncId, v, attrName, &value); + + // no data + if (lowerCase(string(attrName)) == "missing_value" || lowerCase(string(attrName)) == "nodata") + { + missingValue = value; + } + metadata << attrName << " = " << value << endl; } } @@ -1131,6 +1147,7 @@ bool NetCDFHandler::extractVariableMap(int idVar, const Crit3DTime& myTime, std: } default: { + delete[] values; errorStr = "Wrong variable type."; return false; } diff --git a/netcdfHandler/netcdfHandler.h b/netcdfHandler/netcdfHandler.h index 294f6a2f2..1d8a284d9 100644 --- a/netcdfHandler/netcdfHandler.h +++ b/netcdfHandler/netcdfHandler.h @@ -54,6 +54,8 @@ bool isHourly; bool isDaily; + double missingValue; + gis::Crit3DLatLonHeader latLonHeader; NetCDFHandler();