Java tutorial
/* * JGrass - Free Open Source Java GIS http://www.jgrass.org * (C) HydroloGIS - www.hydrologis.com * * This library is free software; you can redistribute it and/or modify it under * the terms of the GNU Library General Public License as published by the Free * Software Foundation; either version 2 of the License, or (at your option) any * later version. * * This library 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 Library General Public License for more * details. * * You should have received a copy of the GNU Library General Public License * along with this library; if not, write to the Free Foundation, Inc., 59 * Temple Place, Suite 330, Boston, MA 02111-1307 USA */ package org.jgrasstools.hortonmachine.modules.hydrogeomorphology.energybalance; import static java.lang.Math.PI; import static java.lang.Math.abs; import static java.lang.Math.acos; import static java.lang.Math.asin; import static java.lang.Math.cos; import static java.lang.Math.exp; import static java.lang.Math.log; import static java.lang.Math.pow; import static java.lang.Math.sin; import static java.lang.Math.tan; import static org.jgrasstools.gears.libs.modules.JGTConstants.C_ice; import static org.jgrasstools.gears.libs.modules.JGTConstants.C_liq; import static org.jgrasstools.gears.libs.modules.JGTConstants.Isc; import static org.jgrasstools.gears.libs.modules.JGTConstants.Lf; import static org.jgrasstools.gears.libs.modules.JGTConstants.Lv; import static org.jgrasstools.gears.libs.modules.JGTConstants.Tf; import static org.jgrasstools.gears.libs.modules.JGTConstants.isNovalue; import static org.jgrasstools.gears.libs.modules.JGTConstants.ka; import static org.jgrasstools.gears.libs.modules.JGTConstants.omega; import static org.jgrasstools.gears.libs.modules.JGTConstants.rho_i; import static org.jgrasstools.gears.libs.modules.JGTConstants.rho_w; import static org.jgrasstools.gears.libs.modules.JGTConstants.sigma; import static org.jgrasstools.gears.libs.modules.JGTConstants.tk; import java.io.File; import java.io.FileInputStream; import java.io.FileOutputStream; import java.io.IOException; import java.io.ObjectInputStream; import java.io.ObjectOutputStream; import java.util.ArrayList; import java.util.Arrays; import java.util.HashMap; import java.util.List; import java.util.Set; import oms3.annotations.Author; import oms3.annotations.Description; import oms3.annotations.Execute; import oms3.annotations.Finalize; 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.geotools.data.simple.SimpleFeatureCollection; import org.geotools.feature.FeatureIterator; import org.jgrasstools.gears.io.eicalculator.EIAreas; import org.jgrasstools.gears.io.eicalculator.EIEnergy; 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; import org.opengis.feature.simple.SimpleFeature; import org.opengis.feature.simple.SimpleFeatureType; import com.vividsolutions.jts.geom.Geometry; import static org.jgrasstools.hortonmachine.i18n.HortonMessages.*; @Description(OMSENERGYBALANCE_DESCRIPTION) @Author(name = OMSENERGYBALANCE_AUTHORNAMES, contact = OMSENERGYBALANCE_AUTHORCONTACTS) @Keywords(OMSENERGYBALANCE_KEYWORDS) @Label(OMSENERGYBALANCE_LABEL) @Name(OMSENERGYBALANCE_NAME) @Status(OMSENERGYBALANCE_STATUS) @License(OMSENERGYBALANCE_LICENSE) public class OmsEnergyBalance extends JGTModel { private static final int GLACIER_SWE = 2000; @Description(OMSENERGYBALANCE_inBasins_DESCRIPTION) @In public SimpleFeatureCollection inBasins; @Description(OMSENERGYBALANCE_fBasinid_DESCRIPTION) @In public String fBasinid = null; @Description(OMSENERGYBALANCE_fBasinlandcover_DESCRIPTION) @In public String fBasinlandcover = null; @Description(OMSENERGYBALANCE_pTrain_DESCRIPTION) @Unit("degrees") @In public double pTrain = 2.0; @Description(OMSENERGYBALANCE_pTsnow_DESCRIPTION) @Unit("degrees") @In public double pTsnow = 0.0; @Description(OMSENERGYBALANCE_pInternaltimestep_DESCRIPTION) @In public double pInternaltimestep = 1; @Description(OMSENERGYBALANCE_pRhosnow_DESCRIPTION) @Unit("kg/m3") @In public double pRhosnow = 400.0; @Description(OMSENERGYBALANCE_pInitswe_DESCRIPTION) @In public double pInitswe = -9999.0; @Description(OMSENERGYBALANCE_pCanopyh_DESCRIPTION) @In public double pCanopyh = 0.0; @Description(OMSENERGYBALANCE_pGlacierid_DESCRIPTION) @In public double pGlacierid = -9999.0; @Description(OMSENERGYBALANCE_pSnowrefv_DESCRIPTION) @In public double pSnowrefv = 0.85; @Description(OMSENERGYBALANCE_pSnowrefir_DESCRIPTION) @In public double pSnowrefir = 0.65; @Description(OMSENERGYBALANCE_tTimestep_DESCRIPTION) @In public int tTimestep; @Description(OMSENERGYBALANCE_tCurrent_DESCRIPTION) @In public String tCurrent = null; @Description(OMSENERGYBALANCE_inRain_DESCRIPTION) @In public HashMap<Integer, double[]> inRain; @Description(OMSENERGYBALANCE_inTemp_DESCRIPTION) @In public HashMap<Integer, double[]> inTemp; @Description(OMSENERGYBALANCE_inWind_DESCRIPTION) @In public HashMap<Integer, double[]> inWind; @Description(OMSENERGYBALANCE_inPressure_DESCRIPTION) @In public HashMap<Integer, double[]> inPressure; @Description(OMSENERGYBALANCE_inRh_DESCRIPTION) @In public HashMap<Integer, double[]> inRh; @Description(OMSENERGYBALANCE_inDtday_DESCRIPTION) @In public HashMap<Integer, double[]> inDtday; @Description(OMSENERGYBALANCE_inDtmonth_DESCRIPTION) @In public HashMap<Integer, double[]> inDtmonth; @Description(OMSENERGYBALANCE_inEnergy_DESCRIPTION) @In public List<EIEnergy> inEnergy = null; @Description(OMSENERGYBALANCE_inAreas_DESCRIPTION) @In public List<EIAreas> inAreas = null; @Description(OMSENERGYBALANCE_pInitsafepoint_DESCRIPTION) @UI(JGTConstants.FILEIN_UI_HINT) @In public String pInitsafepoint; @Description(OMSENERGYBALANCE_pEndsafepoint_DESCRIPTION) @UI(JGTConstants.FILEOUT_UI_HINT) @In public String pEndsafepoint; @Description(OMSENERGYBALANCE_outPnet_DESCRIPTION) @Out public HashMap<Integer, double[]> outPnet; @Description(OMSENERGYBALANCE_outPrain_DESCRIPTION) @Out public HashMap<Integer, double[]> outPrain; @Description(OMSENERGYBALANCE_outPsnow_DESCRIPTION) @Out public HashMap<Integer, double[]> outPsnow; @Description(OMSENERGYBALANCE_outSwe_DESCRIPTION) @Out public HashMap<Integer, double[]> outSwe; @Description(OMSENERGYBALANCE_outNetradiation_DESCRIPTION) @Out public HashMap<Integer, double[]> outNetradiation; @Description(OMSENERGYBALANCE_outNetshortradiation_DESCRIPTION) @Out public HashMap<Integer, double[]> outNetshortradiation; /* * Model's variables definition */ private double[] averageTemperature; // public double[] fullAdigeData; private double defaultTollU0 = 1000.0; private double defaultTollU = 500.0; private double defaultTollW0 = 10.0; private double defaultTollW = 5.0; private double defaultTollTs = 1.0; private double defaultUWiter = 10.0; private double defaultTsiter = 5.0; private int num_ES = -9999; private int num_EI = -9999; /** * latitude latitude in rad. */ private double latitude; /** * longitude longitude in rad. */ private double longitude; /** * standard_time difference from the UTM [hour]. */ private double standard_time; private double E0; private double alpha; private int basinNum = -1; private SafePoint safePoint; private ArrayList<SimpleFeature> basinsFeatures; private double[] Abasin; private HashMap<Integer, Integer> basinid2BasinindexMap; private HashMap<Integer, Integer> basinindex2BasinidMap; private double[][][] EI; private double[][][] A; /* * Full adige vector data contains for every basin in the following order: * 1. net precipitation * 2. net radiation * 3. net short radiation * 4. temperature * 5. relative humidity * 6. wind speed * 7. air pressure */ // private double[] fullAdigeData; private int usoFieldIndex = -1; private List<Integer> usoList; private DateTimeFormatter formatter = JGTConstants.utcDateFormatterYYYYMMDDHHMM; @Execute public void process() throws Exception { outPnet = new HashMap<Integer, double[]>(); outPrain = new HashMap<Integer, double[]>(); outPsnow = new HashMap<Integer, double[]>(); outSwe = new HashMap<Integer, double[]>(); outNetradiation = new HashMap<Integer, double[]>(); outNetshortradiation = new HashMap<Integer, double[]>(); if (safePoint == null) safePoint = new SafePoint(); // retrieve number of bands num_EI = 0; for (EIEnergy energy : inEnergy) { int j = energy.energeticBandId; if (j > num_EI) { num_EI = j; } } num_EI++; num_ES = 0; for (EIAreas area : inAreas) { int j = area.altimetricBandId; if (j > num_ES) { num_ES = j; } } num_ES++; if (basinid2BasinindexMap == null) { // get basin features from feature link basinsFeatures = new ArrayList<SimpleFeature>(); FeatureIterator<SimpleFeature> featureIterator = inBasins.features(); basinNum = inBasins.size(); SimpleFeatureType featureType = inBasins.getSchema(); int basinIdFieldIndex = featureType.indexOf(fBasinid); if (basinIdFieldIndex == -1) { throw new IllegalArgumentException( "The field of the basin id couldn't be found in the supplied basin data."); } if (fBasinlandcover != null) { usoFieldIndex = featureType.indexOf(fBasinlandcover); if (usoFieldIndex == -1) { throw new IllegalArgumentException( "The field of the soil type (usofield) couldn't be found in the supplied basin data."); } } basinid2BasinindexMap = new HashMap<Integer, Integer>(); basinindex2BasinidMap = new HashMap<Integer, Integer>(); pm.beginTask("Read basins data.", inBasins.size()); int index = 0; Abasin = new double[basinNum]; while (featureIterator.hasNext()) { pm.worked(1); SimpleFeature feature = featureIterator.next(); basinsFeatures.add(feature); basinid2BasinindexMap.put(((Number) feature.getAttribute(basinIdFieldIndex)).intValue(), index); basinindex2BasinidMap.put(index, ((Number) feature.getAttribute(basinIdFieldIndex)).intValue()); Geometry basinGeometry = (Geometry) feature.getDefaultGeometry(); Abasin[index] = basinGeometry.getArea() / 1000000.0; // area in km2 as the input // area for energetic and // altimetric bands index++; // read land cover if requested if (usoFieldIndex != -1) { if (usoList == null) { usoList = new ArrayList<Integer>(); } int uso = ((Number) feature.getAttribute(usoFieldIndex)).intValue(); usoList.add(uso); } } featureIterator.close(); pm.done(); } // get rain from scalar link double[] rain = new double[basinNum]; Set<Integer> basinIdSet = inRain.keySet(); pm.beginTask("Read rain data.", basinIdSet.size()); for (Integer basinId : basinIdSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { continue; } double[] value = inRain.get(basinId); if (!isNovalue(value[0])) { rain[index] = value[0]; } else { rain[index] = 0.0; } } pm.done(); // get energy values from scalar link ([12][num_EI][basinNum]) 12 == // 0,1,2,3,4,5,5,4,3,2,1,0 ones at the beginning of the simulation if (EI == null) { EI = new double[12][num_EI][basinNum]; pm.beginTask("Read energy index data.", inEnergy.size()); for (EIEnergy energy : inEnergy) { pm.worked(1); int tempId = energy.basinId; Integer index = basinid2BasinindexMap.get(tempId); if (index == null) { // basinid2BasinindexMap.remove(tempId); continue; } int j = energy.energeticBandId; int k = energy.virtualMonth; int kInverse = 11 - k; EI[k][j][index] = energy.energyValue; EI[kInverse][j][index] = energy.energyValue; } pm.done(); } // get area bande fascie from scalar link ([num_ES][num_EI][basinNum]) ones at the // beginning of the simulation if (A == null) { A = new double[num_ES][num_EI][basinNum]; pm.beginTask("Read area per heigth and band data.", inAreas.size()); HashMap<Integer, HashMap<Integer, HashMap<Integer, Double>>> idbasinMap = new HashMap<Integer, HashMap<Integer, HashMap<Integer, Double>>>(); for (EIAreas area : inAreas) { int idBasin = area.basinId; HashMap<Integer, HashMap<Integer, Double>> idfasceMap = idbasinMap.get(idBasin); if (idfasceMap == null) { idfasceMap = new HashMap<Integer, HashMap<Integer, Double>>(); idbasinMap.put(idBasin, idfasceMap); } int idAltimetricBand = area.altimetricBandId; HashMap<Integer, Double> idbandeMap = idfasceMap.get(idAltimetricBand); if (idbandeMap == null) { idbandeMap = new HashMap<Integer, Double>(); idfasceMap.put(idAltimetricBand, idbandeMap); } int idEnergeticBand = area.energyBandId; double areaValue = area.areaValue; idbandeMap.put(idEnergeticBand, areaValue); pm.worked(1); } pm.done(); for (int i = 0; i < basinNum; i = i + 1) { Integer index = basinindex2BasinidMap.get(i); if (index == null) { basinid2BasinindexMap.remove(i); continue; } HashMap<Integer, HashMap<Integer, Double>> fasceMap = idbasinMap.get(index); for (int j = 0; j < num_ES; j++) { HashMap<Integer, Double> bandeMap = fasceMap.get(j); for (int k = 0; k < num_EI; k++) { A[j][k][i] = bandeMap.get(k); } } } } // get T (temperatures per basin per band) from scalar input link at each time step double[][] T = null; if (inTemp != null) { T = new double[basinNum][num_ES]; pm.beginTask("Read temperature data.", inTemp.size()); Set<Integer> basinIdsSet = inTemp.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inTemp.get(basinId); T[index] = values; } pm.done(); } // get V (wind speed per basin per band) from scalar link at each time step double[][] V = null; if (inWind != null) { V = new double[basinNum][num_ES]; pm.beginTask("Read wind speed data.", inWind.size()); Set<Integer> basinIdsSet = inWind.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inWind.get(basinId); V[index] = values; } } // get P (pressure per basin per band) from scalar link at each time step double[][] P = null; if (inPressure != null) { P = new double[basinNum][num_ES]; pm.beginTask("Read pressure data.", inPressure.size()); Set<Integer> basinIdsSet = inPressure.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inPressure.get(basinId); P[index] = values; } pm.done(); } // get RH (relative humidity per basin per band) from scalar link at each time step double[][] RH = null; if (inRh != null) { RH = new double[basinNum][num_ES]; pm.beginTask("Read humidity data.", inRh.size()); Set<Integer> basinIdsSet = inRh.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inRh.get(basinId); RH[index] = values; } pm.done(); } // get dtday (daily temperature range per basin per band) from scalar link at each time // step double[][] DTd = null; if (inDtday != null) { DTd = new double[basinNum][num_ES]; pm.beginTask("Read dtday data.", inDtday.size()); Set<Integer> basinIdsSet = inDtday.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inDtday.get(basinId); DTd[index] = values; } pm.done(); } // get dtmonth (monthly temperature range per basin per band) from scalar link at each // time step double[][] DTm = null; if (inDtmonth != null) { DTm = new double[basinNum][num_ES]; pm.beginTask("Read dtday data.", inDtmonth.size()); Set<Integer> basinIdsSet = inDtmonth.keySet(); for (Integer basinId : basinIdsSet) { pm.worked(1); Integer index = basinid2BasinindexMap.get(basinId); if (index == null) { // data for a basin that is not considered, ignore it continue; } double[] values = inDtmonth.get(basinId); DTm[index] = values; } pm.done(); } /* * set the current time: day, month and hour */ DateTime currentDatetime = formatter.parseDateTime(tCurrent); int currentMonth = currentDatetime.getMonthOfYear(); int currentDay = currentDatetime.getDayOfMonth(); int currentMinute = currentDatetime.getMinuteOfDay(); double hour = currentMinute / 60.0; System.out.println("ora: " + hour); if (averageTemperature == null) { averageTemperature = new double[2 * basinNum]; } else { Arrays.fill(averageTemperature, 0.0); } /* * these have to be taken from initial values */ if (safePoint.SWE == null) { if (pInitsafepoint != null && new File(pInitsafepoint).exists()) { safePoint = getSafePointData(); } else { safePoint.SWE = new double[num_ES][num_EI][basinNum]; if (pInitswe == -9999.0) { pInitswe = 0.0; } for (int i = 0; i < basinNum; i++) { double sweTmp = pInitswe; if (usoList != null) { int usoTmp = usoList.get(i); if (usoTmp == pGlacierid) { sweTmp = GLACIER_SWE; } } for (int k = 0; k < num_ES; k++) { for (int j = 0; j < num_EI; j++) { safePoint.SWE[j][k][i] = sweTmp; } } } safePoint.U = new double[num_ES][num_EI][basinNum]; safePoint.SnAge = new double[num_ES][num_EI][basinNum]; safePoint.Ts = new double[num_ES][num_EI][basinNum]; } } // this has to be taken from a file, scalarreader // TODO add the input canopyLink for the canopy height for each altimetric band /* * if there is no canopy input matrix for the model create an empty canopy matrix for each elevation band and for each basin */ double[][] canopy = new double[num_ES][basinNum]; for (int i = 0; i < canopy.length; i++) { for (int j = 0; j < canopy[0].length; j++) { canopy[i][j] = pCanopyh; } } checkParametersAndRunEnergyBalance(rain, T, V, P, RH, currentMonth, currentDay, hour, Abasin, A, EI, DTd, DTm, canopy); } private SafePoint getSafePointData() { FileInputStream fis = null; ObjectInputStream in = null; try { fis = new FileInputStream(pInitsafepoint); in = new ObjectInputStream(fis); SafePoint readSafePoint = (SafePoint) in.readObject(); in.close(); return readSafePoint; } catch (IOException ex) { ex.printStackTrace(); } catch (ClassNotFoundException ex) { ex.printStackTrace(); } return null; } @Finalize public void writeSafePoint() { if (pEndsafepoint != null && new File(pEndsafepoint).getParentFile() != null) { FileOutputStream fos = null; ObjectOutputStream out = null; try { fos = new FileOutputStream(pEndsafepoint); out = new ObjectOutputStream(fos); out.writeObject(safePoint); out.close(); } catch (IOException ex) { ex.printStackTrace(); } } } /** * Method to check the input parameters. * * <p> * Due to the huge amount of parameters, this method is used to do * necessary checks and set default values. This is made so the initialize * method doesn't get flooded. * </p> * * @param dtData * @param rain * @param T matrix of temperatures of every altimetric band for every basin.[basin][altim. band] * @param V * @param P * @param RH * @param month * @param day * @param hour * @param Abasin area of the different basins (as defined through the features/geometries) * @param A area for altim and energetic bands. Coming from eicalculator. * @param EI energy index matrix, coming from eicalculator. * @param DTd daily temperature range. * @param DTm monthly temperature range. * @param canopy */ private void checkParametersAndRunEnergyBalance(double[] rain, double[][] T, double[][] V, double[][] P, double[][] RH, double month, double day, double hour, double[] Abasin, double[][][] A, double[][][] EI, double[][] DTd, double[][] DTm, double[][] canopy) { double Dt = ((double) tTimestep / (double) pInternaltimestep) * 60.0; /* * some hardcoded variables */ boolean hasNoStations = false; double zmes_T = 2.0; // quota misura temperatura,pressione e umidita' double zmes_U = 2.0; // quota misura velocita' vento [m] double z0T = 0.005; // [m] roughness length della temperatura double z0U = 0.05; // [m] roughness length del vento double K = ka * ka / ((log(zmes_U / z0U)) * (log(zmes_T / z0T))); double eps = 0.98; // emissivita' neve double Lc = 0.05; // ritenzione capillare double Ksat = 3.0; // 5.55; // conducibilita' idraulica della neve a saturazione double Ks = 5.55E-5; // conducibilita' termica superficiale della neve double aep = 50.0; // albedo extinction parameter (kg/m2==mm) double rho_g = 1600; // densita' del suolo [kg/m3] double De = 0.4; // suolo termicamente attivo double C_g = 890.0; // capacita' termica del suolo [J/(kg K)] double albedo_land = 0.2; double Ts_min = -20.0; double Ts_max = 20.0; // TODO check parameters and add to the model parameter latitude = 46.6 * Math.PI / 180.0; // [rad] longitude = 10.88 * Math.PI / 180.0; // [rad] standard_time = -1.0; // differenza rispetto a UMT [h] /* * start calculations */ sun(hour, day); for (int i = 0; i < basinNum; i++) { calculateEnergyBalance(i, month, hasNoStations, V[i], canopy, T[i], P[i], RH[i], rain, pTrain, pTsnow, Dt, A, Abasin, EI, DTd[i], DTm[i], K, eps, Lc, pRhosnow, Ksat, rho_g, De, C_g, aep, albedo_land, Ks, Ts_min, Ts_max); } } /** * @param i index for the basins list. * @param month * @param hasNoStations * @param windSpeed vettore della velocita' del vento sulle fasce altimetriche per il bacino considerato. * @param canopy * @param T vettore della temperatura sulle fasce altimetriche per il bacino considerato. * @param P vettore della pressione sulle fasce altimetriche per il bacino considerato. * @param RH vettore dell'umidita' relativa sulle fasce altimetriche per il bacino considerato. * @param rain vettore della pioggia * @param Train * @param Tsnow * @param Dt * @param A * @param Abasin * @param EI * @param DTd * @param DTm * @param K * @param eps snow emissivity. * @param Lc capillar ritention. * @param rho_sn snow density [kg/m3] * @param Ksat conducibilita' idraulica della neve a saturazione [kg/(m2 s)]. * @param rho_g densita' del suolo [kg/m3]. * @param De suolo termicamente attivo [m]. * @param C_g capacita' termica del suolo [J/(kg K)]. * @param aep albedo extinction parameter (kg/m2==mm). * @param albedo_land * @param Ks conducibilita' superficiale della neve [m/s]. * @param Ts_min * @param Ts_max * @param SWE snow water equivalent della [banda altimetrica][banda energetica][bacino]. * @param U * @param SnAge snow age [banda altimetrica][banda energetica][bacino] * @param Ts surface temperature */ private void calculateEnergyBalance(int i, double month, boolean hasNoStations, double[] windSpeed, double[][] canopy, double[] T, double[] P, double[] RH, double[] rain, double Train, double Tsnow, double Dt, double[][][] A, double[] Abasin, double[][][] EI, double[] DTd, double[] DTm, double K, double eps, double Lc, double rho_sn, double Ksat, double rho_g, double De, double C_g, double aep, double albedo_land, double Ks, double Ts_min, double Ts_max) { double rho, cp, ea, Psnow, T_snow, Prain, T_rain, Qp, Pnet; double[] tausn = new double[1]; double[] Rsw = new double[1]; double[] Rlwin = new double[1]; double[] netRadiation = new double[1]; double[] netShortRadiation = new double[1]; double[] Wice = new double[1]; double[] Tin = new double[1]; double[] Fliq = new double[1]; double Tsur, Se, H0, L0, R0, M0, U0, W0, H1, L1, R1, M1, U1, W1, U2, W2; int tol, cont; int[] conv = new int[1]; double[][][] SWE = safePoint.SWE; double[][][] SnAge = safePoint.SnAge; double[][][] Ts = safePoint.Ts; double[][][] U = safePoint.U; /* * set first value to the id of the basin */ Integer basinId = basinindex2BasinidMap.get(i); averageTemperature[2 * i] = basinId; // valori medi per bacino // creo un vettore contenente un valore di Snow Water Equivalent per // ogni bacino double tmpPnet = 0.0; double tmpPrain = 0.0; double tmpPsnow = 0.0; double tmpSwe = 0.0; double tmpNetradiation = 0.0; double tmpNetShortRadiation = 0.0; for (int j = 0; j < num_ES; j++) { // per tutte le BANDE // ALTIMETRICHE // riduzione della velocita' del vento dovuta alla canopy windSpeed[j] *= (1.0 - 0.8 * canopy[j][i]); // printf("\n j=%4ld T=%10.5f",j,V[j]); // densita' dell'aria (kg/m3) rho = 1.2922 * (tk / (T[j] + tk)) * (P[j] / 1013.25); // calore specifico a pressione costante (J/(kg K)) cp = 1005.00 + (T[j] + 23.15) * (T[j] + 23.15) / 3364.0; // pressione di vapore dell'aria (mbar) // BEFORE ea = 0.01 * RH[j] * pvap(T[j]); ea = 0.01 * RH[j] * pVap(T[j], P[j]); // precipitazione liquida e solida // se la temperatura della fascia altimetrica maggiore di Train // (parameters) // tutta la precipitazione pioggia if (T[j] > Train) { Prain = rain[i]; Psnow = 0.0; // se la temperatura della fascia altimetrica compresa tra // Train e Tsnow // ci sar una parte di pioggia e una di neve // si trascura l'effetto del vento } else if (T[j] <= Train && T[j] >= Tsnow) { Prain = rain[i] * (T[j] - Tsnow) / (Train - Tsnow); Psnow = rain[i] - Prain; // se la temperatura della fascia altimetrica minore di Tsnow // tutto neve } else { Prain = 0.0; Psnow = rain[i]; } // temperatura della neve T_snow = T[j]; // se la temperatura della neve e' maggiore dello zero assegno a // T_snow lo zero if (T_snow > Tf) T_snow = Tf; T_rain = T[j]; // se la temperatura dell'acqua minore dello zero assegno a T_rain // lo zero if (T_rain < Tf) T_rain = Tf; // calore trasportato dalla precipitazione Qp = (Psnow * C_ice * T_snow + Prain * (Lf + C_liq * T_rain)) / Dt; // Qp[W/m^2] // P[kg/m^2] // c[J/(kg*K)] // Lf[J/kg] for (int k = 0; k < num_EI; k++) { // per tutte le BANDE // ENERGETICHE // se il SWE della banda altimetrica, banda energetica del // bacino i o la prec // nevosa sono maggiori di zero if (SWE[j][k][i] > 0 || Psnow > 0) { // calcolo contenuto di ghiaccio e temperatura della neve da // U calculateTemp(Wice, Tin, Fliq, 1.0E3 * U[j][k][i], SWE[j][k][i], rho_g, De, C_g); // radiazione e albedo tausn[0] = SnAge[j][k][i]; // et della neve // adimensionale // BEFORE radiation(parameters, EI[month][k][i], // alpha, E0, // Ts[j][k][i], Wice[0], T[j], ea, Psnow, // DTd[j], DTm[j], tausn, Rsw, Rlwin); calculateRadiation(EI[(int) month][k][i], Ts[j][k][i], Wice[0], T[j], ea, P[j], Psnow, DTd[j], DTm[j], tausn, Rsw, Rlwin, Dt, aep, albedo_land, netRadiation, netShortRadiation, pSnowrefv, pSnowrefir); // double Dt, double aep, double albedo_land SnAge[j][k][i] = tausn[0]; // riduzione della canopy sui flussi radiativi Rsw[0] *= (1.0 - canopy[j][i]); Rlwin[0] *= (1.0 - canopy[j][i]); // PREDICTOR // temperatura della superficie Tsur = calculateSurfaceTemperature(conv, Ts[j][k][i], Tin[0], T[j], P[j], windSpeed[j], ea, Rsw[0], Rlwin[0], Qp, rho, cp, canopy[j][i], K, rho_sn, Ks, eps, Ts_min, Ts_max); if (Double.isNaN(Tsur) || conv[0] != 1) Tsur = T[j]; // controllo in caso di non // convergenza if (Tsur > Tf) Tsur = Tf; // vedi appunti OK // calore sensibile H0 = K * windSpeed[j] * rho * cp * (T[j] - Tsur); // calore latente L0 = Lv * K * windSpeed[j] * rho * 0.622 * (ea - pVap(Tsur, P[j])) / P[j]; // radiazione onde lunghe uscente R0 = (1.0 - canopy[j][i]) * eps * sigma * pow(Tsur + tk, 4.0); // scioglimento Se = (Fliq[0] / (1.0 - Fliq[0]) - Lc) / (rho_w / rho_sn - rho_w / rho_i - Lc); if (Se < 0) Se = 0.0; M0 = Ksat * pow(Se, 3.0); // flusso di acqua uscente if (Double.isNaN(Se)) M0 = 0.0; if (M0 * Dt > SWE[j][k][i] - Wice[0]) M0 = (SWE[j][k][i] - Wice[0]) / Dt; // aggiornamento W1 = SWE[j][k][i] + rain[i] + (L0 / Lv - M0) * Dt; U1 = U[j][k][i] + 1.0E-3 * (Rsw[0] + Rlwin[0] + Qp - Lf * M0 - R0 + H0 + L0) * Dt; U0 = U1; // aggiorno i parametri di U e W con i valori // calcolati W0 = W1; // di primo tentativo // contatori e controlli cont = 0; tol = 0; // CORRECTOR do { // calcolo Wice, Fliq e Tin da U1 calculateTemp(Wice, Tin, Fliq, 1.0E3 * U1, W1, rho_g, De, C_g); // temperatura della superficie Tsur = calculateSurfaceTemperature(conv, Tsur, Tin[0], T[j], P[j], windSpeed[j], ea, Rsw[0], Rlwin[0], Qp, rho, cp, canopy[j][i], K, rho_sn, Ks, eps, Ts_min, Ts_max); if (Double.isNaN(Tsur) || conv[0] != 1) Tsur = T[j]; // controllo in caso di non // convergenza if (Tsur > Tf) Tsur = Tf; // calore sensibile H1 = K * windSpeed[j] * rho * cp * (T[j] - Tsur); // calore latente L1 = Lv * K * windSpeed[j] * rho * 0.622 * (ea - pVap(Tsur, P[j])) / P[j]; // radiazione onde lunghe uscente R1 = (1.0 - canopy[j][i]) * eps * sigma * pow(Tsur + tk, 4.0); // scioglimento if (Fliq[0] == 1) { M1 = W1 / Dt; U2 = 0.0; W2 = 0.0; tol = 3; } else { Se = (Fliq[0] / (1.0 - Fliq[0]) - Lc) / (rho_w / rho_sn - rho_w / rho_i - Lc); if (Se < 0) Se = 0.0; M1 = Ksat * pow(Se, 3.0); if (Double.isNaN(Se)) M1 = 0.0; if (M1 * Dt > W1 - Wice[0]) M1 = (W1 - Wice[0]) / Dt; // aggiornamento W2 = SWE[j][k][i] + rain[i] + (0.5 * (L0 / Lv - M0) + 0.5 * (L1 / Lv - M1)) * Dt; U2 = U[j][k][i] + 1.0E-3 * (Rsw[0] + Rlwin[0] + Qp + 0.5 * (-Lf * M0 - R0 + H0 + L0) + 0.5 * (-Lf * M1 - R1 + H1 + L1)) * Dt; cont += 1; // controllo convergenza if (cont == 1 && abs(U2 - U1) < defaultTollU0 && abs(W2 - W1) < defaultTollW0) tol = 1; if (cont > 1 && abs(U2 - U1) < defaultTollU && abs(W2 - W1) < defaultTollW) tol = 2; if (conv[0] != 1) tol = 0; // se Tsur non converge, non va bene // la soluzione // aggiornamento U1 = U2; W1 = W2; } } while (tol == 0 && cont <= defaultUWiter); // se non c'e' convergenza cerco una soluzione esplicita // (approssimata) if (tol == 0) { // solo frazione liquida: ho trovato i // valori e aggiorno i dati U[j][k][i] = U0; SWE[j][k][i] = W0; Ts[j][k][i] = Tsur; Pnet = M0 * Dt; } else if (tol == 3) { // non ho convergenza e setto a zero // tutto U[j][k][i] = 0.0; SWE[j][k][i] = 0.0; Ts[j][k][i] = 0.0; Pnet = M1 * Dt; } else { // frazione liquida e solida: ho trovato i // valori e aggiorno i dati U[j][k][i] = U1; SWE[j][k][i] = W1; Ts[j][k][i] = Tsur; Pnet = (0.5 * M0 + 0.5 * M1) * Dt; } // se tutto lo SWE e' liquido o se c' troppo poco SWE // metto tutto lo SWE nella Pn if (1.0E3 * U[j][k][i] >= Lf * SWE[j][k][i] || SWE[j][k][i] <= 0) { if (SWE[j][k][i] < 0) SWE[j][k][i] = 0.0; Pnet += SWE[j][k][i]; if (Pnet < 0) Pnet = 0.0; SWE[j][k][i] = 0.0; U[j][k][i] = 0.0; } } else { Pnet = Prain; Tsur = T[j]; Ts[j][k][i] = Tsur; tausn[0] = 0.0; calculateRadiation(EI[(int) month][k][i], Ts[j][k][i], 0.0, T[j], ea, P[j], Psnow, DTd[j], DTm[j], tausn, Rsw, Rlwin, Dt, aep, albedo_land, netRadiation, netShortRadiation, pSnowrefv, pSnowrefir); // for( m = 2; m <= 40; m++ ) { // logg[m] = -1.0; // } } safePoint.SWE[j][k][i] = SWE[j][k][i]; safePoint.U[j][k][i] = U[j][k][i]; safePoint.Ts[j][k][i] = Ts[j][k][i]; safePoint.SnAge[j][k][i] = SnAge[j][k][i]; // calcolo i valori medi per bacino di SWE, Pnet, Prain, Psnow tmpSwe = tmpSwe + SWE[j][k][i] * (A[j][k][i] / Abasin[i]); tmpPnet = tmpPnet + Pnet * (A[j][k][i] / Abasin[i]); tmpPrain = tmpPrain + Prain * (A[j][k][i] / Abasin[i]); tmpPsnow = tmpPsnow + Psnow * (A[j][k][i] / Abasin[i]); tmpNetradiation = tmpNetradiation + netRadiation[0] * (A[j][k][i] / Abasin[i]); tmpNetShortRadiation = tmpNetShortRadiation + netShortRadiation[0] * (A[j][k][i] / Abasin[i]); // System.out.println("swe = " + fullAdigeData[8 * i + 8]); } averageTemperature[2 * i + 1] += T[j]; } outSwe.put(basinId, new double[] { tmpSwe }); outPnet.put(basinId, new double[] { tmpPnet }); outPrain.put(basinId, new double[] { tmpPrain }); outPsnow.put(basinId, new double[] { tmpPsnow }); outNetradiation.put(basinId, new double[] { tmpNetradiation }); outNetshortradiation.put(basinId, new double[] { tmpNetShortRadiation }); // System.out.println("rad media= " + fullAdigeData[8 * i + 2]); // System.out.println("short media= " + fullAdigeData[8 * i + 3]); averageTemperature[2 * i + 1] /= num_ES; } private void calculateTemp(double[] Wice, double[] Tin, double[] Fliq, double U, double SWE, double rho_g, double De, double C_g) { if (U <= 0) { Wice[0] = SWE; Tin[0] = U / (SWE * C_ice + rho_g * De * C_g); Fliq[0] = 0.0; } else if (U > 0 && U <= Lf * SWE) { Wice[0] = SWE - U / Lf; Tin[0] = 0.0; Fliq[0] = (SWE - Wice[0]) / SWE; } else { Wice[0] = 0.0; Tin[0] = (U - Lf * SWE) / (rho_g * De * C_g + SWE * C_liq); Fliq[0] = 1.0; } } private void calculateRadiation(double EI, double Ts, double Wice, double Ta, double ea, double P, double Psnow, double DTd, double DTm, double[] tausn, double[] Rsw, double[] Rlwin, double Dt, double aep, double albedo_land, double[] netRadiation, double[] netShortRadiation, double avo, double airo) { double coszen, r1, r2, r3, fzen, fage, avd, avis, aird, anir, albedo, rr, AtmTrans; double eps_clsky, CF, eps; double bb = 2.0, cv = 0.2, cr = 0.5, a = 0.8, c = 2.4, b; double diff2glob; // COSINE OF ZENITHAL ANGLE coszen = EI * sin(alpha); // ALBEDO // effect snow surface temperature r1 = exp(5000.0 * (1.0 / (Tf + tk) - 1.0 / (Ts + tk))); // effect melt and refreezing r2 = pow(r1, 10); if (r2 > 1.0) r2 = 1.0; // effect of dirt r3 = 0.03; // non-dimensional snow age: 10 mm of snow precipitation restore snow // age Dt(s) tausn[0] = (tausn[0] + (r1 + r2 + r3) * Dt * 1.0E-6) * (1.0 - Psnow / 10.0); if (tausn[0] < 0.0) tausn[0] = 0.0; // dipendence from solar angle if (coszen < 0.5) { fzen = 1.0 / bb * ((bb + 1.0) / (1.0 + 2.0 * bb * coszen) - 1.0); } else { fzen = 0.0; } // dipendence from snow age fage = tausn[0] / (1.0 + tausn[0]); // diffuse visible albedo avd = (1.0 - cv * fage) * avo; // global visible albedo avis = avd + 0.4 * fzen * (1.0 - avd); // diffuse near infared albedo aird = (1.0 - cr * fage) * airo; // global near infared albedo anir = aird + 0.4 * fzen * (1.0 - aird); // albedo is taken as average albedo = (avis + anir) / 2.0; // Linear transition from snow albedo to bare ground albedo if (Psnow == 0) { albedo = albedo_land; } else if (Wice < aep) { rr = (1.0 - Wice / aep) * exp(-Wice * 0.5 / aep); albedo = rr * albedo_land + (1.0 - rr) * albedo; } // NET SHORTWAVE RADIATION if (DTm < 0) DTm = 0.0; if (DTd < 0) DTd = 0.0; // ADDED a = 0.48 + 0.29 * (1013.25 / P) * sin(alpha); if (a > 1) a = 1.0; b = 0.036 * exp(-0.154 * DTm); AtmTrans = a * (1.0 - exp(-b * pow(DTd, c))); // ADDED // ratio diffuse to global radiation (Erbs et al., 1982) if (AtmTrans <= 0.22) { diff2glob = 1.0 - 0.09 * AtmTrans; } else if (AtmTrans > 0.22 && AtmTrans <= 0.80) { diff2glob = 0.9511 - 0.1604 * AtmTrans + 4.388 * pow(AtmTrans, 2.0) - 16.638 * pow(AtmTrans, 3.0) + 12.336 * pow(AtmTrans, 4.0); } else { diff2glob = 0.165; } // Rsw=direct+diffuse Rsw[0] = Isc * E0 * AtmTrans * (1.0 - albedo) * ((1.0 - diff2glob) * coszen + diff2glob * sin(alpha)); // INCOMING LONGWAVE RADIATION eps_clsky = 1.08 * (1.0 - exp(-pow(ea, (Ta + tk) / 2016.0))); CF = 1 - AtmTrans / a; eps = (1.0 - CF) * eps_clsky + CF; Rlwin[0] = eps * sigma * pow(Ta + tk, 4.0); // net radiation netRadiation[0] = Rsw[0] + Rlwin[0] - albedo * Rsw[0] - eps * sigma * pow(Ts + tk, 4.0); // System.out.println("Albedo " + albedo); // System.out.println("Radiazione netta: " + netRadiation[0]); netShortRadiation[0] = (1 - albedo) * Rsw[0]; // System.out.println("Net Short: " + netShortRadiation[0]); } /** * Calcola la temperatura della superficie della neve (C). * * <p> * Risolve il bilancio di energia alla superficie linearizzando e * iterando, secondo il metodo di Tarboton. * </p> * * @param conv * @param Ts temperatura della superficie nevosa di primo tentativo (o dell'istante precedente). * @param Tin temperatura della neve di primo tentativo (o dell'istante precedente). * @param Ta temperatura dell'aria. * @param P pressione (mbar). * @param V velocita' del vento (m/s). * @param ea pressione di vapore in aria (mbar). * @param Rsw * @param Rlwin * @param Qp * @param rho * @param cp * @param Fcanopy * @param K * @param rho_sn * @param Ks * @param eps * @param Ts_min * @param Ts_max * @return */ private double calculateSurfaceTemperature(int[] conv, double Ts, double Tin, double Ta, double P, double V, double ea, double Rsw, double Rlwin, double Qp, double rho, double cp, double Fcanopy, double K, double rho_sn, double Ks, double eps, double Ts_min, double Ts_max) { double a, b, Ts0; short cont; // coefficienti non dipendenti dalla temperatura della superficie a = Rsw + Rlwin + Qp + K * V * Ta * rho * cp + rho_sn * C_ice * Ks * Tin + 0.622 * K * V * Lv * rho * ea / P; b = rho_sn * C_ice * Ks + K * V * rho * cp; cont = 0; do { Ts0 = Ts; Ts = (a - 0.622 * K * V * Lv * rho * (pVap(Ts0, P) - Ts * dpVap(Ts0, P)) / P + 3.0 * Fcanopy * sigma * eps * pow(Ts0 + tk, 4.0)) / (b + 0.622 * dpVap(Ts0, P) * K * V * Lv * rho / P + 4.0 * Fcanopy * eps * sigma * pow(Ts0 + tk, 3.0)); cont += 1; } while (abs(Ts - Ts0) > defaultTollTs && cont <= defaultTsiter); // controlli if (abs(Ts - Ts0) > defaultTollTs) { conv[0] = 0; // non converge } else { conv[0] = 1; // converge } if (Ts < Ts_min || Ts > Ts_max) conv[0] = -1; // fuori dai limiti di ammissibilita' Ts0 = Ts; return Ts0; } /** * calcola la pressione di vapore [mbar] a saturazione in dipendenza * dalla temperatura [gradi Celsius]. * * @param T * @param P * @return la pressione di vapore [mbar] a saturazione */ private double pVap(double T, double P) { double A = 6.1121 * (1.0007 + 3.46E-6 * P); double b = 17.502; double c = 240.97; double e = A * exp(b * T / (c + T)); return e; } /** * Calcola la derivata della pressione di vapore [mbar] a saturazione * rispetto dalla temperatura [gradi Celsius]. * * @param T * @param P * @return derivata della pressione di vapore. */ private double dpVap(double T, double P) { double A = 6.1121 * (1.0007 + 3.46E-6 * P); double b = 17.502; double c = 240.97; double De = (A * exp(b * T / (c + T))) * (b / (c + T) - b * T / pow(c + T, 2.0)); return De; } /** * @param hour * @param day * @param E0 earth-sun distance correction. * @param alpha solar height (complementar to zenith angle), [rad]. */ private void sun(double hour, double day) { // standard latitude according to standard time double lst = standard_time * PI / 12.0; // correction sideral time double G = 2.0 * PI * (day - 1.0) / 365.0; double Et = 0.000075 + 0.001868 * cos(G) - 0.032077 * sin(G) - 0.014615 * cos(2 * G) - 0.04089 * sin(2 * G); // local time double lh = hour + (longitude - lst) / omega + Et / omega; // earth-sun distance correction E0 = 1.00011 + 0.034221 * cos(G) + 0.00128 * sin(G) + 0.000719 * cos(2 * G) + 0.000077 * sin(2 * G); // solar declination double D = 0.006918 - 0.399912 * cos(G) + 0.070257 * sin(G) - 0.006758 * cos(2 * G) + 0.000907 * sin(2 * G) - 0.002697 * cos(3 * G) + 0.00148 * sin(3 * G); // Sunrise and sunset with respect to 12pm [hour] double Thr = (acos(-tan(D) * tan(latitude))) / omega; if (lh >= 12 - Thr && lh <= 12 + Thr) { // alpha: solar height (complementar to zenith angle), [rad] alpha = asin(sin(latitude) * sin(D) + cos(latitude) * cos(D) * cos(omega * (12.0 - lh))); } else { alpha = 0.0; } } }