EtpPenmanMonteithFAO.java

/*
 * This file is part of Indicators.
 *
 * Indicators 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.
 *
 * Indicators 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 Indicators. If not, see <https://www.gnu.org/licenses/>.
 */
package fr.inrae.agroclim.indicators.model.data.climate;

import java.util.Arrays;
import java.util.HashSet;
import java.util.Set;

import fr.inrae.agroclim.indicators.model.data.Variable;
import lombok.EqualsAndHashCode;
import lombok.Setter;

/**
 * Calculator for ETP (aka ET0, ETR) according to Penman-Monteith, FAO (Allen et
 * al., 1998; Smith, 1992; Hargreaves, 1994). Rnl is neglected.
 *
 * http://www.fao.org/docrep/X0490E/x0490e06.htm#penman%20monteith%20equation
 *
 * Last changed : $Date$
 *
 * @author $Author$
 * @version $Revision$
 */
@EqualsAndHashCode(
        callSuper = false,
        of = {"windAt10m"}
)
public final class EtpPenmanMonteithFAO implements EtpCalculator {
    /**
     * UUID for Serializable.
     */
    private static final long serialVersionUID = 4872847709393688047L;

    /**
     * Correction for wind speed measured at 10m instead of 2m.
     */
    private static final double WIND_10M_TO_2M_COEFFICIENT = 0.7;

    /**
     * If wind speed is measured at 10m or 2m.
     *
     * Case for DRIAS, Meteo-France stations. Wind speed at INRAE stations is at
     * 2m.
     */
    @Setter
    private boolean windAt10m = true;

    /**
     * Constructor.
     */
    public EtpPenmanMonteithFAO() {
    }

    /**
     * Constructor.
     *
     * @param windSpeedAt10m
     *            True: wind speed is measured at 10 m. False: 2 m
     */
    public EtpPenmanMonteithFAO(final boolean windSpeedAt10m) {
        this.windAt10m = windSpeedAt10m;
    }

    @Override
    public EtpPenmanMonteithFAO clone() {
        EtpPenmanMonteithFAO clone = new EtpPenmanMonteithFAO();
        clone.windAt10m = windAt10m;
        return clone;
    }

    @Override
    public Double compute(final ClimaticDailyData data) throws IllegalArgumentException {
        if (data.getRadiation() == null) {
            throw new IllegalArgumentException("error.radiation.null");
        }
        if (data.getRh() == null) {
            throw new IllegalArgumentException("error.rh.null");
        }
        if (data.getWind() == null) {
            throw new IllegalArgumentException("error.wind.null");
        }
        final double rg = data.getRadiation();
        final double tmoy = data.getTmean();
        final double rh = data.getRh();
        final double wind;
        if (windAt10m) {
            wind = data.getWind() * WIND_10M_TO_2M_COEFFICIENT;
        } else {
            wind = data.getWind();
        }
        // albedo for green grass
        final double albedo = 0.23;

        // rayonnement net en Mj/m² alors qu'on le reçoit en W/m²
        // #7012 : radiation retrieved from climatic db in MJ/m2
        //         (converted by Season sql query if needed)

        final double rgNet = rg * (1 - albedo);

        // saturation vapour pressure (kPa)
        final double es = 0.6108 * Math.exp(17.27 * tmoy / (237.3 + tmoy));
        // actual vapour pressure deficit of the air (kPa)
        final double ea = es * rh / 100;
        // vapour pressure deficit
        final double dpv = es - ea;
        // slope of saturation vapour pressure curve (kPa/°C)
        final double delta = 4098 * es / ((237.3 + tmoy) * (237.3 + tmoy));
        // psychrometric constant (kPa/°C)
        final double psychro = 0.000665 * 101.3 * Math.pow(
                (293 - 0.0065 * 1) / 293, 5.26);
        final double numerator = 0.408 * delta * rgNet
                + psychro * 900 / (tmoy + 273) * wind * dpv;
        final double denominator = delta + psychro * (1 + 0.34 * wind);
        return numerator / denominator;
    }

    @Override
    public Set<Variable> getVariables() {
        Set<Variable> variables = new HashSet<>();
        variables.addAll(Arrays.asList(Variable.RADIATION, Variable.RH,
                Variable.TMEAN, Variable.WIND));
        return variables;
    }
}