Java tutorial
/* * This file is part of JGrasstools (http://www.jgrasstools.org) * (C) HydroloGIS - www.hydrologis.com * * JGrasstools is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package org.jgrasstools.hortonmachine.modules.hydrogeomorphology.etp; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_AUTHORCONTACTS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_AUTHORNAMES; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_KEYWORDS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_LABEL; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_LICENSE; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_NAME; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_STATUS; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_UI; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inNetradiation_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inPressure_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inRh_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inShortradiation_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inSwe_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inTemp_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inVegetation_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_inWind_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_outEtp_DESCRIPTION; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSPENMANETP_tCurrent_DESCRIPTION; import java.util.HashMap; import java.util.Map.Entry; import java.util.Set; import oms3.annotations.Author; import oms3.annotations.Description; import oms3.annotations.Execute; import oms3.annotations.In; import oms3.annotations.Keywords; import oms3.annotations.Label; import oms3.annotations.License; import oms3.annotations.Name; import oms3.annotations.Out; import oms3.annotations.Status; import oms3.annotations.UI; import oms3.annotations.Unit; import org.jgrasstools.gears.io.adige.VegetationLibraryRecord; import org.jgrasstools.gears.libs.exceptions.ModelsIllegalargumentException; import org.jgrasstools.gears.libs.modules.JGTConstants; import org.jgrasstools.gears.libs.modules.JGTModel; import org.joda.time.DateTime; import org.joda.time.format.DateTimeFormatter; @Description(OMSPENMANETP_DESCRIPTION) @Author(name = OMSPENMANETP_AUTHORNAMES, contact = OMSPENMANETP_AUTHORCONTACTS) @Keywords(OMSPENMANETP_KEYWORDS) @Label(OMSPENMANETP_LABEL) @Name(OMSPENMANETP_NAME) @Status(OMSPENMANETP_STATUS) @License(OMSPENMANETP_LICENSE) @UI(OMSPENMANETP_UI) public class OmsPenmanEtp extends JGTModel { // @Description("Baricenter elevation of the HillSlope for every basin on which to calculate.") // @Unit("m") // @In // public HashMap<Integer, double[]> inElevations; @Description(OMSPENMANETP_inVegetation_DESCRIPTION) @In public HashMap<Integer, VegetationLibraryRecord> inVegetation; @Description(OMSPENMANETP_inNetradiation_DESCRIPTION) @Unit("W/m2") @In public HashMap<Integer, double[]> inNetradiation; @Description(OMSPENMANETP_inShortradiation_DESCRIPTION) @In public HashMap<Integer, double[]> inShortradiation; @Description(OMSPENMANETP_inTemp_DESCRIPTION) @Unit("C") @In public HashMap<Integer, double[]> inTemp; @Description(OMSPENMANETP_inRh_DESCRIPTION) @In public HashMap<Integer, double[]> inRh; @Description(OMSPENMANETP_inWind_DESCRIPTION) @In public HashMap<Integer, double[]> inWind; @Description(OMSPENMANETP_inPressure_DESCRIPTION) @In public HashMap<Integer, double[]> inPressure; @Description(OMSPENMANETP_inSwe_DESCRIPTION) @In public HashMap<Integer, double[]> inSwe; @Description(OMSPENMANETP_tCurrent_DESCRIPTION) @In public String tCurrent; @Description(OMSPENMANETP_outEtp_DESCRIPTION) @Unit("mm/day") @Out public HashMap<Integer, double[]> outEtp; private DateTimeFormatter formatter = JGTConstants.utcDateFormatterYYYYMMDDHHMM; private static final double Z0_SNOW = 0.010; private static final double CLOSURE = 4000.0; private static final double RSMAX = 5000.0; // Pa private static final double VPDMINFACTOR = 0.1; private static final double A_SVP = 0.61078; private static final double B_SVP = 17.269; private static final double C_SVP = 237.3; private static final double CP_PM = 1013.0; /* specific heat of moist air J/kg/C (Handbook of Hydrology) */ // private static final double PS_PM = 101300.0; /* sea level air pressure in Pa */ // private static final double LAPSE_PM = -0.006; /* environmental lapse rate in C/m */ // private static final int SECPHOUR = 3600; /* seconds per hour */ private static final int SEC_PER_DAY = 86400; /* seconds per day */ private static final double HUGE_RESIST = 1.e20; /* largest allowable double number */ private static final double VON_K = 0.41; /* Von Karman constant for evapotranspiration */ private static final double ZREF = 2.0; // reference height for wind speed @Execute public void penman() { checkNull(inPressure, inTemp, inRh, inWind, inSwe, inVegetation, inShortradiation, inNetradiation); outEtp = new HashMap<Integer, double[]>(); DateTime currentTimestamp = formatter.parseDateTime(tCurrent); int monthOfYear = currentTimestamp.getMonthOfYear(); Set<Entry<Integer, double[]>> elevSet = inTemp.entrySet(); for (Entry<Integer, double[]> entry : elevSet) { Integer basinId = entry.getKey(); // double elevation = inElevations.get(basinId)[0]; double tair = entry.getValue()[0]; double pressure = inPressure.get(basinId)[0]; double relativeHumidity = inRh.get(basinId)[0]; double wind = inWind.get(basinId)[0]; double snowWaterEquivalent = inSwe.get(basinId)[0]; double shortRadiation = inShortradiation.get(basinId)[0]; double netRadiation = inNetradiation.get(basinId)[0]; VegetationLibraryRecord vegetation = inVegetation.get(basinId); double displacement = vegetation.getDisplacement(monthOfYear); double roughness = vegetation.getRoughness(monthOfYear); double rs = vegetation.getMinStomatalResistance(); double RGL = vegetation.getRgl(); double lai = vegetation.getLai(monthOfYear); double rarc = vegetation.getArchitecturalResistance(); /* * set the pressure in Pascal instead of in hPa as the input link gives */ pressure = pressure / 100; double vpd = svp(tair) - (relativeHumidity * 100 / svp(tair)); // vpd in KPa double ra = calcAerodynamic(displacement, roughness, ZREF, wind, snowWaterEquivalent); // CONSIDER THE SOIL RESISTANCE NULL // // calculate gsm_inv: soil moisture stress factor // double criticalSoilMoisture = 0.33; // fraction of soil moisture content at the // critical // // point // double wiltingPointSoilMoisture = 0.133; // double waterContentCriticalPoint = criticalSoilMoisture * maxMoisture; // double waterContentWiltingPoint = wiltingPointSoilMoisture * maxMoisture; // double gsm_inv; // soil moisture stress factor // if (soilMoisture >= waterContentCriticalPoint) { // gsm_inv = 1.0; // } else if (soilMoisture >= waterContentWiltingPoint) { // gsm_inv = (soilMoisture - waterContentWiltingPoint) / (waterContentCriticalPoint - // waterContentWiltingPoint); // } else { // gsm_inv = 0.0; // } /* calculate the slope of the saturated vapor pressure curve in Pa/K */ double slope = svp_slope(tair) * 1000.0; /* factor for canopy resistance based on photosynthesis */ double dayFactor; /* calculate resistance factors (Wigmosta et al., 1994) */ double f = 0.0; if (rs > 0.) { if (RGL < 0) { throw new ModelsIllegalargumentException("Invalid value of RGL for the current class.", this); } else if (RGL == 0) { f = shortRadiation; } else { f = shortRadiation / RGL; } dayFactor = (1. + f) / (f + rs / RSMAX); } else dayFactor = 1.; /* factor for canopy resistance based on temperature */ double tFactor = .08 * tair - 0.0016 * tair * tair; tFactor = (tFactor <= 0.0) ? 1e-10 : tFactor; /* factor for canopy resistance based on vpd */ double vpdFactor = 1 - vpd / CLOSURE; vpdFactor = (vpdFactor < VPDMINFACTOR) ? VPDMINFACTOR : vpdFactor; /* calculate canopy resistance in s/m */ // double rc = rs / (lai * gsm_inv * tFactor * vpdFactor) * dayFactor; double rc = rs / (lai * tFactor * vpdFactor) * dayFactor; rc = (rc > RSMAX) ? RSMAX : rc; // double h; /* scale height in the atmosphere (m) */ // /* calculate scale height based on average temperature in the column */ // h = 287 / 9.81 * ((tair + 273.15) + 0.5 * (double) elevation * LAPSE_PM); // // /* use hypsometric equation to calculate p_z, assume that virtual temperature is // equal // air_temp */ // pz = PS_PM * Math.exp(-(double) elevation / h); // instead of calculating the pressure in h.adige it is possible to read the // interpolated // pressure from input link // pz is the surface air pressure /* calculate latent heat of vaporization. Eq. 4.2.1 in Handbook of Hydrology, assume Ts is Tair */ double lv = 2501000 - 2361 * tair; /* calculate the psychrometric constant gamma (Pa/C). Eq. 4.2.28. Handbook of Hydrology */ double gamma = 1628.6 * pressure / lv; /* calculate factor to be applied to rc/ra */ /* calculate the air density (in kg/m3), using eq. 4.2.4 Handbook of Hydrology */ double r_air = 0.003486 * pressure / (275 + tair); /* calculate the Penman-Monteith evaporation in mm/day (by not dividing by * the density of water (~1000 kg/m3)), the result ends up being in mm instead of m */ double evap = ((slope * netRadiation + r_air * CP_PM * vpd / ra) / (lv * (slope + gamma * (1 + (rc + rarc) / ra))) * SEC_PER_DAY) / 24.0; if (vpd >= 0.0 && evap < 0.0) evap = 0.0; outEtp.put(basinId, new double[] { evap }); } } /** * This routine computes the gradient of d(svp)/dT using Handbook of Hydrology eqn 4.2.3 * @param tair * @return saturated vapor pressure slope */ private double svp_slope(double temp) { double satVaporPressureSlope = (B_SVP * C_SVP) / ((C_SVP + temp) * (C_SVP + temp)) * svp(temp); return satVaporPressureSlope; } /** * This routine computes the saturated vapor pressure using Handbook of Hydrology eqn 4.2.2 (Pressure in kPa) * * @param temp temperature. * @return saturated vapor pressure. */ private double svp(double temp) { double SVP; SVP = A_SVP * Math.exp((B_SVP * temp) / (C_SVP + temp)); if (temp < 0) SVP *= 1.0 + .00972 * temp + .000042 * temp * temp; return (SVP); } /** * Calculates the aerodynamic resistance for the vegetation layer. * * <p>Calculates the aerodynamic resistance for the vegetation layer, and * the wind 2m above the layer boundary.</p> * <p>The values are normalized based on a reference height wind * speed, Uref, of 1 m/s. To get wind speeds and aerodynamic resistances for * other values of Uref, you need to multiply the here calculated wind * speeds by Uref and divide the here calculated aerodynamic resistances * by Uref</p> * * @param displacement * @param roughness * @param Zref reference height for windspeed. * @param windSpeed * @return the aerodynamic resistance for the vegetation layer. */ private double calcAerodynamic(double displacement, double roughness, double Zref, double windSpeed, double snowWaterEquivalent) { double ra = 0.0; double d_Lower; double K2; double Z0_Lower; double tmp_wind; // only a value of these quantities are input of the method: // - wind speed // - relative humidity // - Zref // - ra tmp_wind = windSpeed; K2 = VON_K * VON_K; if (displacement > Zref) Zref = displacement + Zref + roughness; /* No OverStory, thus maximum one soil layer */ // for bare soil // if (iveg == Nveg) { // Z0_Lower = Z0_SOIL; //Z0_SOIL is the soil roughness // d_Lower = 0; // } else { //with vegetation Z0_Lower = roughness; d_Lower = displacement; // } if cycle is deleted thinking that bare soil is a vegetation type with roughness and // displacement /* With snow on the surface */ if (snowWaterEquivalent > 0) { windSpeed = Math.log((2. + Z0_SNOW) / Z0_SNOW) / Math.log(Zref / Z0_SNOW); ra = Math.log((2. + Z0_SNOW) / Z0_SNOW) * Math.log(Zref / Z0_SNOW) / K2; } else { /* No snow on the surface*/ windSpeed = Math.log((2. + Z0_Lower) / Z0_Lower) / Math.log((Zref - d_Lower) / Z0_Lower); ra = Math.log((2. + (1.0 / 0.63 - 1.0) * d_Lower) / Z0_Lower) * Math.log((2. + (1.0 / 0.63 - 1.0) * d_Lower) / (0.1 * Z0_Lower)) / K2; } if (tmp_wind > 0.) { windSpeed *= tmp_wind; ra /= tmp_wind; } else { windSpeed *= tmp_wind; ra = HUGE_RESIST; pm.message("Aerodinamic resistance is set to the maximum value!"); } return ra; } }