bc9d664f by Roger Marquez

Merge branch 'feature/co2' into develop

2 parents 2aaf8ad9 53bc7534
......@@ -184,6 +184,10 @@ public class PotentialCropGrowthAdapter extends AnnotatedAdapter {
@Optional
@Input @Output public double BioYield;
@Description("the daily atmospheric level of CO2")
@Units("ppm")
@Input public double co2;
@Description("Amount N Added to Residue Pool After Harvesting")
@Units("kgN/ha")
@Output public double Addresidue_pooln;
......@@ -210,6 +214,9 @@ public class PotentialCropGrowthAdapter extends AnnotatedAdapter {
@Description("DLAI (Senescence)")
@Output public double dlai;
@Description("Co2Effect")
@Output public double co2effect;
private Map<HRU, PotentialCropGrowthSWAT> swatComponentMap = new ConcurrentHashMap<>();
@Override
......@@ -230,7 +237,9 @@ public class PotentialCropGrowthAdapter extends AnnotatedAdapter {
private void executeSWAT() {
PotentialCropGrowthSWAT component = swatComponentMap.get(hru);
if (component == null) {
//component = new PotentialCropGrowthSWAT();
component = new PotentialCropGrowthSWAT();
swatComponentMap.put(hru, component);
}
......@@ -288,6 +297,9 @@ public class PotentialCropGrowthAdapter extends AnnotatedAdapter {
component.FPHUact = FPHUact;
component.nfert = nfert;
// co2
component.co2 = co2;
component.execute();
plantStateReset = component.plantStateReset;
......@@ -316,6 +328,7 @@ public class PotentialCropGrowthAdapter extends AnnotatedAdapter {
topt = component.topt;
fr3N = component.fr3N;
dlai = component.dlai;
co2effect = component.co2eff;
}
private void executeUPGM() {
......
......@@ -21,6 +21,9 @@
package crop;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import java.util.logging.Logger;
public class PotentialCropGrowthSWAT {
......@@ -88,6 +91,20 @@ public class PotentialCropGrowthSWAT {
public double fracharvestn;
private double LAI_delta;
private static final boolean DO_CO2_BIOMASS = false;
// For co2 response.
public double co2;
public double co2eff;
// Assuming hardcoded values for initial testing will be ok. Eventually this could be user input.
// c4 crops
private static final List<Double> CO2_X_C4CROPS = Arrays.asList(0.0, 220.0, 280.0, 330.0, 400.0, 490.0, 570.0, 750.0, 990.0, 9999.0);
private static final List<Double> CO2_Y_C4CROPS = Arrays.asList(0.00, 0.85, 0.95, 1.00, 1.02, 1.04, 1.05, 1.06, 1.07, 1.08);
// c5 crops
private static final List<Double> CO2_X_C3CROPS = Arrays.asList(0.0, 220.0, 330.0, 440.0, 550.0, 660.0, 770.0, 880.0, 990.0, 9999.0);
private static final List<Double> CO2_Y_C3CROPS = Arrays.asList(0.00, 0.71, 1.00, 1.08, 1.17, 1.25, 1.32, 1.38, 1.43, 1.50);
public void execute() {
BioYield = 0;
......@@ -103,7 +120,12 @@ public class PotentialCropGrowthSWAT {
calc_phu();
calc_lai();
calc_biomass();
if (DO_CO2_BIOMASS) {
calc_co2_biomass();
} else {
co2eff = 1.0;
calc_biomass();
}
CanHeightAct = calc_canopy();
calc_root();
calc_nuptake();
......@@ -207,6 +229,43 @@ public class PotentialCropGrowthSWAT {
BioOpt_delta = (dormancy) ? 0 : BioOpt_delta;
}
private void calc_co2_biomass() {
co2eff = 1.0; // default = no impact.
int indexLess = 0;
int indexPlusOne = 0;
if (idc == 4) {
// c4 crops (?maybe?)
// either exact position of co2 value, or pos - 1 in array.
int indx = Collections.binarySearch(CO2_X_C4CROPS, co2);
if (indx < 0) {
indexLess = Math.abs(indx) - 2;
indexPlusOne = Math.abs(indx) - 1;
} else {
indexLess = indx;
indexPlusOne = Math.abs(indx) + 1;
}
co2eff = (co2 - CO2_X_C4CROPS.get(indexLess)) * (CO2_Y_C4CROPS.get(indexPlusOne) - CO2_Y_C4CROPS.get(indexLess))
/ (CO2_X_C4CROPS.get(indexPlusOne) - CO2_X_C4CROPS.get(indexLess)) + CO2_Y_C4CROPS.get(indexLess);
} else {
// either exact position of co2 value, or pos - 1 in array.
int indx = Collections.binarySearch(CO2_X_C3CROPS, co2);
if (indx < 0) {
indexLess = Math.abs(indx) - 2;
indexPlusOne = Math.abs(indx) - 1;
} else {
indexLess = indx;
indexPlusOne = Math.abs(indx) + 1;
}
co2eff = (co2 - CO2_X_C3CROPS.get(indexLess)) * (CO2_Y_C3CROPS.get(indexPlusOne) - CO2_Y_C3CROPS.get(indexLess))
/ (CO2_X_C3CROPS.get(indexPlusOne) - CO2_X_C3CROPS.get(indexLess)) + CO2_Y_C3CROPS.get(indexLess);
}
double Hphosyn = 0.5 * solRad * (1 - Math.exp(LExCoef * LAI)) * co2eff; // intercepted photosynthetically active radiation (MJ/m^2)
BioOpt_delta = rue * Hphosyn;
BioOpt_delta = (dormancy) ? 0 : BioOpt_delta;
}
// canopy cover
// canopy cover is expressed as a function of LAI
private double calc_canopy() {
......
/*
* This file is part of the AgroEcoSystem-Watershed (AgES-W) model component
* collection. AgES-W components are derived from multiple agroecosystem models
* including J2K and J2K-SN (FSU-Jena, DGHM, Germany), SWAT (USDA-ARS, USA),
* WEPP (USDA-ARS, USA), RZWQM2 (USDA-ARS, USA), and others.
*
* The AgES-W model is free software; you can redistribute the model and/or
* modify the components 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 potet;
import lib.Climate;
public class PenmanMonteithCO2TransRatio {
private static final double CP = 1.031E-3;
private static final double RSS = 150;
public String tempRes;
public double wind;
public double tmean;
public double rhum;
public double netRad;
public double actRsc0;
public double elevation;
public double area;
public double LAI;
public double actEffH;
public int idc;
public double co2;
public double potET;
public double actET;
public void execute() {
double tempFactor = 0;
if ("d".equals(tempRes)) {
tempFactor = 86400;
} else if ("h".equals(tempRes)) {
tempFactor = 3600;
}
double abs_temp = Climate.absTemp(tmean, "C");
double delta_s = Climate.slopeOfSaturationPressureCurve(tmean);
double pz = Climate.atmosphericPressure(elevation, abs_temp);
double est = Climate.saturationVapourPressure(tmean);
double ea = Climate.vapourPressure(rhum, est);
double latH = Climate.latentHeatOfVaporization(tmean);
double psy = Climate.psyConst(pz, latH);
double G = 0.1 * netRad;
double vT = Climate.virtualTemperature(abs_temp, pz, ea);
double pa = Climate.airDensityAtConstantPressure(vT, pz);
double Letp = calcETAllen(delta_s, netRad, G, pa, CP, est, ea, calcRa(actEffH, wind), calcRs(LAI, actRsc0, RSS),
psy, tempFactor, calcStomatalRatio(idc, co2, LAI));
potET = Letp / latH;
// convert mm to l
potET *= area;
// negative potET is not allowed
if (potET < 0) {
potET = 0;
}
actET = 0; // reset actET
}
private static double calcETAllen(double ds, double netRad, double G, double pa,
double CP, double est, double ea, double ra, double rs, double psy, double tempFactor, double stomatalRatio) {
return (ds * (netRad - G) + ((pa * CP * (est - ea) / ra) * tempFactor)) / (ds + psy * (1 + (rs / ra * stomatalRatio)));
}
private static double calcRa(double eff_height, double wind_speed) {
double ra;
if (wind_speed <= 0) {
wind_speed = 0.5;
}
if (eff_height < 10) {
ra = (9.5 / wind_speed) * Math.pow(Math.log(2 / (0.1 * eff_height)), 2);
} else {
ra = 20 / (0.1681 * wind_speed);
}
return ra;
}
private static double calcRs(double LAI, double rsc, double rss) {
double A = Math.pow(0.7, LAI);
return 1. / (((1 - A) / rsc) + ((A / rss)));
}
private static double calcStomatalRatio(int idc, double co2, double LAI) {
/*
Programmed using the approach outlined in the DSSAT-CSM 4.6 source file TRANS.FOR.
This approach was adapted for use with the PenmanMonteith PET calculation.
*/
// if there is no LAI, no plant, so no change
if (LAI < 0.01) {
return 1.0;
}
double rb = 10.0;
double rlf = 0.0;
double rlfc = 0.0;
/* DSSAT-CSM 4.6 Comments on following code:
!-----------------------------------------------------------------------
! RLF = Leaf stomatal resistance at 330.0 ppm CO2, s/m
! RLFC = Leaf stomatal resistance at other CO2 conc., s/m
! (Allen, 1986), Plant responses to rising CO2.
! CO2 = CO2 Conc of the increased ATM case
! CO2 = 330 Ambient CO2
!-----------------------------------------------------------------------
*/
// "IF (INDEX('MZ,ML,SG,SC,SW,BM,BH,BR,NP,SI',CROP) .GT. 0) THEN"
// Note: maize/corn/millet/sorghum seem to fall under IDC 4,
// I am not sure if the other crops listed as c-4 crops fall in this range too.
if (idc == 4) {
// c4 crops
// DSSAT-CSM 4.6-> "EQ 7 from Allen (1986) for corn."
rlf = (1.0 / (0.0328 - (5.49 * 10E-5 * 330.0) + (2.95E-8 * Math.pow(330.0, 2)))) + rb;
rlfc = (1.0 / (0.0328 - (5.49 * 10E-5 * co2) + (2.95E-8 * Math.pow(co2, 2)))) + rb;
} else {
// c3 crops
rlf = 9.72 + (0.0757 * 330.0) + rb;
rlfc = 9.72 + (0.0757 * co2) + rb;
}
// DSSAT-CSM assumes 2.88 for the reference crop LAI:
//" FAO-56 assumption of 2.88 LAI reference"
rlf = rlf / (0.5 * 2.88);
rlfc = rlfc / (0.5 * 2.88);
return rlfc / rlf;
}
}
/*
* This file is part of the AgroEcoSystem-Watershed (AgES-W) model component
* collection. AgES-W components are derived from multiple agroecosystem models
* including J2K and J2K-SN (FSU-Jena, DGHM, Germany), SWAT (USDA-ARS, USA),
* WEPP (USDA-ARS, USA), RZWQM2 (USDA-ARS, USA), and others.
*
* The AgES-W model is free software; you can redistribute the model and/or
* modify the components 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 potet;
import ages.types.HRU;
import annotations.Author;
import annotations.Bibliography;
import annotations.Documentation;
import annotations.Execute;
import annotations.Keywords;
import annotations.License;
import annotations.Role;
import static annotations.Role.PARAMETER;
import annotations.SourceInfo;
import annotations.Status;
import annotations.VersionInfo;
import crop.Crop;
import gov.usda.jcf.annotations.Description;
import gov.usda.jcf.annotations.Input;
import gov.usda.jcf.annotations.Output;
import gov.usda.jcf.annotations.Units;
import gov.usda.jcf.core.Context;
import gov.usda.jcf.core.adapters.AnnotatedAdapter;
@Description("PenmanMonteith PET calculation with CO2 impact through stomatal conductance ratio")
@Author(name = "Olaf David, Manfred Fink, James C. Ascough II", contact = "jim.ascough@ars.usda.gov")
@Keywords("Insert keywords")
@Bibliography("Insert bibliography here")
@VersionInfo("$ID:$")
// Nathan, do we need to change this?
@SourceInfo("http://javaforge.com/scm/file/3979484/default/src/potet/PenmanMonteith.java")
@License("http://www.gnu.org/licenses/gpl.html")
@Status(Status.TESTED)
@Documentation("src/potet/PenmanMonteith.xml")
public class PenmanMonteithCO2TransRatioAdapter extends AnnotatedAdapter {
@Description("daily or hourly time steps [d|h]")
@Units("d | h")
@Role(PARAMETER)
@Input public String tempRes;
@Description("wind")
@Input public double wind;
@Description("Mean Temperature")
@Units("C")
@Input public double tmean;
@Description("Relative Humidity")
@Input public double rhum;
@Description("Daily net radiation")
@Units("MJ/m2")
@Input public double netRad;
@Description("state variable rsc0")
@Input public double actRsc0;
@Description("HRU")
@Input public HRU hru;
@Description("LAI")
@Input public double LAI;
@Description("effective height")
@Input public double actEffH;
@Description("crop")
@Input public Crop crop;
@Description("the daily atmospheric level of CO2")
@Units("ppm")
@Input public double co2;
@Description("HRU potential Evapotranspiration")
@Units("mm")
@Output public double potET;
@Description("HRU actual Evapotranspiration")
@Units("mm")
@Output public double actET;
@Override
protected void run(Context context) {
PenmanMonteithCO2TransRatio component = new PenmanMonteithCO2TransRatio();
component.tempRes = tempRes;
component.wind = wind;
component.tmean = tmean;
component.rhum = rhum;
component.netRad = netRad;
component.actRsc0 = actRsc0;
component.elevation = hru.elevation;
component.area = hru.area;
component.LAI = LAI;
component.actEffH = actEffH;
// required for co2 impact
component.idc = crop.idc;
component.co2 = co2;
component.execute();
potET = component.potET;
actET = component.actET;
}
/**
* OMS execute statement for testing before JCF
*/
@Execute
public void execute() {
run(null);
}
}
Styling with Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!