diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 540d4617563552d17aedcf3df14aaf5f830a68df..9932cd8157e2674af333efd3a89ccf764843a33f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -13,7 +13,6 @@ cache:
paths:
- .cache/pip
-
.config_freva: &config_freva
image: ghcr.io/freva-clint/freva-dev:latest
tags:
@@ -42,34 +41,42 @@ func-test:
<<: *config_freva
stage: test
script:
- # âš ï¸ To make the following jobs work, please have a look at the `climate-change-profile` plugin at the link below:
- # https://gitlab.dkrz.de/ch1187/plugins4freva/climate-change-profile/-/blob/main/.gitlab-ci.yml?ref_type=heads#L53
- export SCRATCH=/tmp
- mkdir -p /home/freva
- cd /home/freva
- tar -xJvf $CI_PROJECT_DIR/data/btu-data.tar.xz
+ - tar -xJvf $CI_PROJECT_DIR/data/cmip-data.tar.xz
- cd $CI_PROJECT_DIR
- make prepare-freva # This command prepares the configuration files for the plugin
- ./data-crawler solr add -ds nukleus-btu
- export EVALUATION_SYSTEM_PLUGINS=$CI_PROJECT_DIR,pet-wrapper-api
- freva-plugin --list
- # TODO: After generalization to all avaiable dataset, we need to provide a func-test here
- - freva-plugin pet product=ceu-3 institute=clmcom-btu model=ec-earth-consortium-ec-earth3-veg-clmcom-btu-icon-2-6-5-rc-nukleus-x2yn2-v1 experiment=historical ensemble=r1i1p1f1 start_date=1967-01-01 end_date=1969-12-31 region=2,6,46,54 output_dir=$CI_PROJECT_DIR/output debug=True -vvv
+ - freva-plugin pet project=nukleus product=ceu-3 institute=clmcom-btu model=ec-earth-consortium-ec-earth3-veg-clmcom-btu-icon-2-6-5-rc-nukleus-x2yn2-v1 experiment=historical ensemble=r1i1p1f1 start_date=1967-01-01 end_date=1969-12-31 region=2,6,46,54 output_dir=$CI_PROJECT_DIR/output debug=True -vvv
- mkdir -p build/func_test
- FREVA_JOB_ID=$(ls plugin_env/freva_output/freva/freva/output/pet/ | head -n 1)
- cd plugin_env/freva_output/freva/freva/output/pet/$FREVA_JOB_ID
- |
- desired_out_files=("APOTEVAP_S_1967010100-1968010100.ncz" "APOTEVAP_S_1968010100-1969010100.ncz" "APOTEVAP_S_1969010100-1970010100.ncz")
+ desired_out_files=("evspsblpot_1967010100-1968010100.ncz" "evspsblpot_1968010100-1969010100.ncz" "evspsblpot_1969010100-1970010100.ncz")
for file in ${desired_out_files[@]}; do
- if [ ! -f APOTEVAP_S/$file ]; then
+ if [ ! -f evspsblpot/$file ]; then
+ exit 1
+ fi
+ done
+ - cd $CI_PROJECT_DIR
+ - freva-plugin pet project=cmip6 product=cmip institute=awi model=awi-esm-1-recom experiment=esm-hist ensemble=r1i1p1f1 start_date=1900-01-01 end_date=1902-12-31 region=5,15,47,55 output_dir=$CI_PROJECT_DIR/output debug=True -vvv
+ - FREVA_JOB_ID=$(ls plugin_env/freva_output/freva/freva/output/pet/ | head -n 2 | tail -1)
+ - cd plugin_env/freva_output/freva/freva/output/pet/$FREVA_JOB_ID
+ - |
+ desired_out_files=("evspsblpot_1900010100-1901010100.ncz" "evspsblpot_1901010100-1902010100.ncz" "evspsblpot_1902010100-1903010100.ncz")
+ for file in ${desired_out_files[@]}; do
+ if [ ! -f evspsblpot/$file ]; then
exit 1
fi
done
- echo '<html><body>' >> index.html
- echo '<h1>PET plugin Results:</h1>' >> index.html
- - echo '<a href="'$CI_PROJECT_URL'/-/jobs/'$CI_JOB_ID'/artifacts/browse/build/func_test/APOTEVAP_S/" target="_blank">Link to data</a>' >> index.html
+ - echo '<a href="'$CI_PROJECT_URL'/-/jobs/'$CI_JOB_ID'/artifacts/browse/build/func_test/evspsblpot/" target="_blank">Link to data</a>' >> index.html
- echo '</body></html>' >> index.html
- - ls -la
- cp -r * $CI_PROJECT_DIR/build/func_test
artifacts:
when: always
@@ -91,8 +98,8 @@ prod-test:
- module load cdo/2.2.2-gcc-11.2.0 nco/5.0.6-gcc-11.2.0
- export EVALUATION_SYSTEM_PLUGINS=$CI_PROJECT_DIR,pet-wrapper-api
- freva-plugin --list
- # TODO: After generalization to all avaiable dataset, we need to provide a prod-test here
- - freva-plugin pet product=ceu-3 institute=clmcom-btu model=ec-earth-consortium-ec-earth3-veg-clmcom-btu-icon-2-6-5-rc-nukleus-x2yn2-v1 experiment=historical ensemble=r1i1p1f1 start_date=1970-01-01 end_date=1971-12-31 region=2,6,46,54 caption="PET production test from CI piepline https://gitlab.dkrz.de/ch1187/plugins4freva/pet/-/jobs/$CI_JOB_ID"-vvv
+ - freva-plugin pet project=nukleus product=ceu-3 institute=clmcom-btu model=ec-earth-consortium-ec-earth3-veg-clmcom-btu-icon-2-6-5-rc-nukleus-x2yn2-v1 experiment=historical ensemble=r1i1p1f1 start_date=1970-01-01 end_date=1971-12-31 region=2,6,46,54 caption="PET production test from CI piepline https://gitlab.dkrz.de/ch1187/plugins4freva/pet/-/jobs/$CI_JOB_ID" -vvv
+ - freva-plugin pet project=cmip6 product=cmip institute=awi model=awi-esm-1-recom experiment=esm-hist ensemble=r1i1p1f1 start_date=1900-01-01 end_date=1902-12-31 region=5,15,47,55 caption="PET production test from CI piepline https://gitlab.dkrz.de/ch1187/plugins4freva/pet/-/jobs/$CI_JOB_ID" -vvv
- mkdir -p build/prod_test
- REGIKLIM_JOB_ID=$(ls -t /work/$ACCOUNT/regiklim-work/$GITLAB_USER_LOGIN/regiklim-ces/output/pet/ | head -n 1)
- cd /work/$ACCOUNT/regiklim-work/$GITLAB_USER_LOGIN/regiklim-ces/output/pet/$REGIKLIM_JOB_ID
diff --git a/CHANGELOG.md b/CHANGELOG.md
index dcc1b18a5214fdd4a7a78cbeb426433e7e3d0e12..94cd01d2265b0ec3487c3da0c96da6e1e909bd1a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,3 +22,16 @@ SPDX-License-Identifier: CC0-1.0
* Outsourced definition of calculation functions into seperate file: `src/pet/pet_functions.sh`
* Outsourced basic configurations into seperate file: `src/pet/utils.sh`
* Outsourced additional expression files into own directory: `src/pet/helper`
+
+## v2025.01.28: Intelligent variable handling
+
+### Added
+ * Added support for projects besides NUKLEUS via project facet
+ * Added calculation of daily vapour pressure (pvap) if data not available: `src/pet/helper/calc_pvap`
+ * Added project-based selection of orographic data if no orographic data is given for selected collection of facet options
+ * Added feature checking availability of required variables based on selection of facet options and providing suggestions of suitable facet options if varaibles are missing
+ * Added prod-test and func-test for CMIP6 data to test calculation of pvap
+ * Added test data for CMIP6
+### Changed
+ * Set of varaibles are now located within a dictionary taking into account that different projects have different capitalization of varaible names
+
diff --git a/data/btu-data.tar.xz b/data/btu-data.tar.xz
index 3045c169865b61cb346d98e1d78b4f2e28008a26..ef2ea8675bf3135219c4f673775c53e7ee532eb7 100644
Binary files a/data/btu-data.tar.xz and b/data/btu-data.tar.xz differ
diff --git a/data/cmip-data.tar.xz b/data/cmip-data.tar.xz
new file mode 100644
index 0000000000000000000000000000000000000000..9c2b9aaad1002e2a043051b7c5ecc479ea5aaa1c
Binary files /dev/null and b/data/cmip-data.tar.xz differ
diff --git a/data/cmip-data.tar.xz.license b/data/cmip-data.tar.xz.license
new file mode 100644
index 0000000000000000000000000000000000000000..07c923760f9ef68e5243b6df768b8ce71313d240
--- /dev/null
+++ b/data/cmip-data.tar.xz.license
@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2025 Brandenburgische Technische Universität Cottbus-Senftenberg
+
+SPDX-License-Identifier: CC0-1.0
diff --git a/pet-wrapper-api.py b/pet-wrapper-api.py
index e8e2f721dfb2a9d716492ac8a94b20612585be28..2ef69fef5e01d5138e8a50c7ac365b638529bc0d 100644
--- a/pet-wrapper-api.py
+++ b/pet-wrapper-api.py
@@ -11,9 +11,12 @@ PET plugin API wrapper.
# libraries as needed
######################################################
+import sys
import os
import json
+import re
from pathlib import Path
+from datetime import datetime
import subprocess
# from evaluation_system.api.exceptions import PluginError
@@ -40,9 +43,9 @@ class PET(plugin.PluginAbstract):
# Define the details for the plugin
######################################################
- __version__ = (2025, 1, 10)
+ __version__ = (2025, 1, 28)
__short_description__: str = (
- "Calculation of daily Potential Evatransporation (only NUKLEUS dataset)"
+ "Calculation of daily Potential Evatransporation"
)
__category__: Optional[str] = "postproc"
__tags__: Optional[List[str]] = ["Potential Evatranspiration"]
@@ -52,23 +55,30 @@ class PET(plugin.PluginAbstract):
}
__long_description__: Optional[
str
- ] = "This plugin calculates daily Potential Evapotranspiration (PET) based on the methodology outlined in FAO Irrigation and Drainage Paper No. 56 (Crop Evapotranspiration) by Allen, Pereira et al., 1998.\n\nFor more details, refer to the official publication: https://www.fao.org/4/x0490e/x0490e00.htm \n\n ATTENTION: At the moment this plugin only works with NUKLEUS project dataset. But in the near future it would be generalized to current available projects in data-browser."
+ ] = "This plugin calculates daily Potential Evapotranspiration (PET) based on the methodology outlined in FAO Irrigation and Drainage Paper No. 56 (Crop Evapotranspiration) by Allen, Pereira et al., 1998.\n\nFor more details, refer to the official publication: https://www.fao.org/4/x0490e/x0490e00.htm"
######################################################
# Define the parameters for the plugin * Mandatory
######################################################
__parameters__ = ParameterDictionary(
+ SolrField(
+ name="project",
+ facet="project",
+ mandatory=True,
+ help="Project facet for input data search.",
+ predefined_facets={
+ "variable": ["hurs", "sfcWind"],
+ "time_frequency": ["1day"],
+ "project": ["cmip6", "cordex", "nukleus", "reanalysis"],
+ }),
SolrField(
name="product",
facet="product",
mandatory=True,
- default="ceu-3",
help="Product facet for input data search.",
predefined_facets={
- "project": ["nukleus"],
- "product": ["ceu-3", "eur-11"],
- "variable": ["rsds"],
+ "variable": ["hurs", "sfcWind"],
"time_frequency": ["1day"],
}),
SolrField(
@@ -77,8 +87,7 @@ class PET(plugin.PluginAbstract):
mandatory=True,
help="Institute facet for input data search.",
predefined_facets={
- "project": ["nukleus"],
- "variable": ["rsds"],
+ "variable": ["hurs", "sfcWind"],
"time_frequency": ["1day"],
}),
SolrField(
@@ -87,8 +96,7 @@ class PET(plugin.PluginAbstract):
mandatory=True,
help="Model facet for input data search.",
predefined_facets={
- "project": ["nukleus"],
- "variable": ["rsds"],
+ "variable": ["hurs", "sfcWind"],
"time_frequency": ["1day"],
}),
SolrField(
@@ -97,8 +105,7 @@ class PET(plugin.PluginAbstract):
mandatory=True,
help="Experiment facet for input data search.",
predefined_facets={
- "project": ["nukleus"],
- "variable": ["rsds"],
+ "variable": ["hurs", "sfcWind"],
"time_frequency": ["1day"],
}),
SolrField(
@@ -107,8 +114,7 @@ class PET(plugin.PluginAbstract):
mandatory=True,
help="ensemble facet for input data search.",
predefined_facets={
- "project": ["nukleus"],
- "variable": ["rsds"],
+ "variable": ["hurs", "sfcWind"],
"time_frequency": ["1day"],
}),
Date(
@@ -140,6 +146,59 @@ class PET(plugin.PluginAbstract):
),
)
+ # dictionary containing variable names for all projects (necessary because of different capitalization)
+ project_vars = ["hurs", "tasmax", "tasmin", "sfcWind", "rsds"]
+
+ def check_variables(self, config_dict):
+ '''
+ Checks if given set of facet values provides data for all required variables and suggestion more suitable options if this is not the case.
+
+ Parameters:
+ config_dict (dict): Dictionary containing facet values to search database
+ '''
+
+ facets = ["project", "product", "institute", "model", "experiment", "ensemble"]
+
+ facet_suggestions = {}
+ missing_vars = list(filter(lambda x: len(config_dict["files"][x]) == 0, config_dict["files"]))
+
+ search_args = {key: config_dict[key] for key in facets}
+
+ if (set(self.project_vars).issubset(missing_vars)):
+ sys.exit(logger.error("There is no data for all required variables. Please check if the chosen data covers the defined time period."))
+ elif (len(missing_vars) > 0):
+
+ # remove facet stepwise and filter for possible alternative facet values
+ removed_facets = []
+ for facet in reversed(facets):
+ del search_args[facet]
+ removed_facets.append(facet)
+ freva_results = dict(freva.facet_search(**search_args, facet="variable"))["variable"]
+
+ # if alternative facets found covering all variables stop the loop and provide suggestion
+ if (set(self.project_vars).issubset(freva_results)):
+
+ # check all possible entries of latest removed facet, all entries covering all variables will be saved as suggestions
+ facet_suggestions = []
+ if (facet == "project"):
+ facet_entries = ["cmip6", "cordex", "nukleus", "reanalysis"] # only search in relevant project datasets
+ else:
+ facet_entries = list(freva.facet_search(**search_args,facet=facet)[facet])
+
+ for sugg in facet_entries:
+ test_search_args = search_args.copy()
+ test_search_args[facet] = sugg
+ facet_variables = dict(freva.facet_search(**test_search_args, facet="variable"))["variable"]
+
+ if (not set(self.project_vars).issubset(facet_variables)):
+ facet_suggestions.append(sugg)
+
+ # if at least one enrty covers all required variables, stop the search, otherwise continue
+ if (len(facet_suggestions) > 0):
+ error_msg = "There is data for some of the required variables missing (" + ", ".join(missing_vars) + ") ! Consider choose a different " + facet + ". \n Suggestions: " + ", ".join(facet_suggestions)
+ sys.exit(logger.error(error_msg))
+ break
+
def find_orog(self, config_dict):
'''
Retrieves a list of orographic files based on given facets and returns first entry as suggestion
@@ -150,8 +209,7 @@ class PET(plugin.PluginAbstract):
orog_files[0] (string): First entry of list of found files
'''
facets = ["product", "institute", "model", "experiment", "ensemble"]
- #search_args = {"project": config_dict["project"], "variable": "orog"} for later if project facte is given
- search_args = {"project": "nukleus", "variable": "orog"}
+ search_args = {"project": config_dict["project"], "variable": "orog"}
orog_files = list(freva.databrowser(**search_args))
for facet in facets:
@@ -166,6 +224,37 @@ class PET(plugin.PluginAbstract):
return orog_files[0]
+ def extract_timeperiod(self, config_dict):
+ '''
+ Retrieves start and end date based on given datasets and exchanges provided start and/or end date in configuration if data for users selected time period is not available
+
+ Parameters:
+ config_dict (dict): Dictionary containing facet values as well as user given time period
+ '''
+
+ # all variables should cover the same time period, we only need one to extract time period -> tasmax
+ config_dict["files"]["tasmax"].sort(key=lambda f: int(''.join(filter(str.isdigit, f)))) #sorting list of files based on given start year in file name
+
+ # cmorized daily data includes time period of data file as 'YYYYMMDD-YYYYMMDD' within the file name
+ date1 = re.search(r'\d{4}\d{2}\d{2}-\d{4}\d{2}\d{2}', config_dict["files"]["tasmax"][0]).group(0).split("-",1)[0]
+ date2 = re.search(r'\d{4}\d{2}\d{2}-\d{4}\d{2}\d{2}', config_dict["files"]["tasmax"][-1]).group(0).split("-",1)[1]
+
+ start_date = datetime.strptime(date1, '%Y%m%d').strftime('%Y-%m-%d')
+ end_date = datetime.strptime(date2, '%Y%m%d').strftime('%Y-%m-%d')
+
+ time_adapted = False
+
+ if (start_date > config_dict["start_date"]):
+ config_dict["start_date"] = start_date
+ time_adapted = True
+
+ if (end_date < config_dict["end_date"]):
+ config_dict["end_date"] = end_date
+ time_adapted = True
+
+ if (time_adapted):
+ print('INFO Adapted time period based on data availability to ',config_dict["start_date"],'-', config_dict["end_date"])
+
def find_files(self, config_dict, variables: List[str], time_frequency: str, time: Optional[str] = None):
'''
Selects data files based on the given parameters using the freva databrowser
@@ -180,7 +269,7 @@ class PET(plugin.PluginAbstract):
'''
search_keys = (
- #"project",
+ "project",
"product",
"institute",
"model",
@@ -198,7 +287,7 @@ class PET(plugin.PluginAbstract):
search_args["variable"] = variable
files = list(freva.databrowser(**search_args))
file_dict[variable] = files
-
+
return file_dict
def prepare_html_output(self, path_output):
@@ -334,13 +423,15 @@ class PET(plugin.PluginAbstract):
with open(html_output_path, "w") as f:
f.write(html_content)
-
def run_tool(self, config_dict=None):
- config_dict["files"] = self.find_files(config_dict, ["tasmax", "tasmin", "sfcWind", "pvap", "rsds"], "day", config_dict["start_date"] + " to " + config_dict["end_date"])
+ config_dict["files"] = self.find_files(config_dict, self.project_vars, "day", config_dict["start_date"] + " to " + config_dict["end_date"])
config_dict["alti_file"] = self.find_orog(config_dict)
if "alti_file" not in config_dict or not config_dict["alti_file"]:
logger.error("Can't find altitude file")
+
+ self.check_variables(config_dict)
+ self.extract_timeperiod(config_dict)
out_dir = self._special_variables.substitute({"d": "$USER_OUTPUT_DIR"})
out_dir = Path(out_dir["d"]) / str(self.rowid)
diff --git a/src/pet/calculate_pet.sh b/src/pet/calculate_pet.sh
index 86098d8fa6ff43ab8842b0a34e4374681252f3c5..57542e94beaaf15633635e8d3ccba9de91bea72a 100755
--- a/src/pet/calculate_pet.sh
+++ b/src/pet/calculate_pet.sh
@@ -52,9 +52,9 @@ DEBUG=$(jq '.debug' $1 | sed 's/\"//g')
MAXPP=30
WIND_H_NAME=Wind_height # variable name for height of wind speed (needed for log. wind profile)
-PET_NAME=APOTEVAP_S
-PET_LNAME="surface potential evapotranspiration amount"
-PET_SNAME="water_potential_evapotranspiration_amount"
+PET_NAME=evspsblpot
+PET_LNAME="potential daily evapotranspiration"
+PET_SNAME="potential_daily_evapotranspiration"
PET_UNIT="kg m-2"
##################################################################################################
@@ -113,12 +113,16 @@ log_debug 'Checking availability of all necessary variables...'
VARIABLES=$(jq '.files | keys | .[]' $1)
for VAR in ${VARIABLES[@]}; do
- if [ $(jq ".files.$VAR | length " $1) > 0 ]; then
+ if [ $(jq ".files.$VAR | length " $1) -gt 0 ]; then
log_debug "${VAR//\"}...OK"
else
- log_error 'ERROR: No data for '${VAR}' available!'
- log_error 'Please select a dataset containg all necessary varaiables ('${VARIABLES}')!'
- exit 1
+ if [ ${VAR//\"} == "pvap" ]; then
+ log_debug 'No data for '${VAR}' available. Data for '${VAR}' will be calculated.'
+ else
+ log_error 'ERROR: No data for '${VAR}' available!'
+ log_error 'Please select a dataset containg all necessary varaiables!'
+ exit 1
+ fi
fi
done
@@ -127,6 +131,12 @@ echo ''
log_debug 'Preparing files for checks...'
for VAR in ${VARIABLES[@]}; do
+ # skip preparation if no pvap data available
+ if [ ${VAR//\"} == "pvap" ] && [ $(jq ".files.${VAR} | length" $1) == 0 ]; then
+ log_debug "Skipping preparation of missing pvap data"
+ continue
+ fi
+
# requirement: allfiles for a specific variable are located within one directory
VAR_FILES=$(jq ".files.${VAR//\"}[]" $1 | sed 's/\"//g' | grep '.nc')
TMP_FILE=${TMP_DIR}'/'${VAR//\"}'_'${START_YEAR}'-'${END_YEAR}'.nc'
@@ -136,7 +146,7 @@ for VAR in ${VARIABLES[@]}; do
log_debug 'Selecting time period...'
TIME_FILE=${TMP_DIR}'/'${VAR//\"}'_timeperiod.nc'
if [ ! ${START_DATE} == null ] && [ ! ${END_DATE} == null ]; then
- cdo seldate,${START_DATE}'T00:00:00',${END_DATE}'T00:00:00' -mergetime ${VAR_FILES} ${TIME_FILE}
+ cdo -s seldate,${START_DATE}'T00:00:00',${END_DATE}'T00:00:00' -mergetime ${VAR_FILES} ${TIME_FILE}
fi
if [ ! -f ${TIME_FILE} ]; then
@@ -151,7 +161,7 @@ for VAR in ${VARIABLES[@]}; do
REGION_FILE=${TMP_DIR}'/'${VAR//\"}'_region.nc'
log_debug 'Selecting region...'
- cdo -sellonlatbox,${REGION} ${TMP_FILE} ${REGION_FILE}
+ cdo -s -sellonlatbox,${REGION} ${TMP_FILE} ${REGION_FILE}
if [ ! -f ${REGION_FILE} ]; then
log_error "Error selecting region and creating: ${REGION_FILE}"
@@ -166,6 +176,13 @@ done
echo ''
log_debug 'Checking variable names, heights and units...'
for VAR in ${VARIABLES[@]}; do
+
+ # skipping if pvap data is missing
+ if [ ${VAR//\"} == "pvap" ] && [ $(jq ".files.${VAR} | length" $1) == 0 ]; then
+ log_debug "Skipping pvap"
+ continue
+ fi
+
VAR_FILE_NAME=${TMP_DIR}'/'${VAR//\"}'_'${START_YEAR}'-'${END_YEAR}'.nc'
X_SIZE=$(get_x_size ${VAR_FILE_NAME})
@@ -187,6 +204,22 @@ for VAR in ${VARIABLES[@]}; do
convert_unit ${VAR_FILE_NAME} ${VAR//\"}
done
+### CALCULATION OF PVAP ###
+if [ $(jq ".files.pvap | length" $1) == 0 ]; then
+ log_debug "Calculating pvap"
+ TASMAX_FILE=${TMP_DIR}'/tasmax_'${START_YEAR}'-'${END_YEAR}'.nc'
+ TASMIN_FILE=${TMP_DIR}'/tasmin_'${START_YEAR}'-'${END_YEAR}'.nc'
+ HURS_FILE=${TMP_DIR}'/hurs_'${START_YEAR}'-'${END_YEAR}'.nc'
+ CALC_TMP_FILE=${TMP_DIR}'/tmp_calc_data.nc'
+
+ ncks -A ${TASMAX_FILE} ${CALC_TMP_FILE}
+ ncks -A ${TASMIN_FILE} ${CALC_TMP_FILE}
+ ncks -A ${HURS_FILE} ${CALC_TMP_FILE}
+
+ PVAP_MYEXP_FILE=$(dirname "$0")/helper/calc_pvap
+ cdo -s -exprf,${PVAP_MYEXP_FILE} ${CALC_TMP_FILE} ${TMP_DIR}'/pvap_'${START_YEAR}'-'${END_YEAR}'.nc'
+fi
+
### ALTITUDE ###
ALTITUDE_NAME=$(cdo codetab ${ALTI_FILE} | awk '{print $2;}')
ALTITUDE_UNIT=$(cdo partab ${ALTI_FILE} | grep "units" | grep -o -P '(?<=").*(?=")')
@@ -196,7 +229,7 @@ log_debug "${ALTITUDE_NAME} [${ALTITUDE_UNIT}]"
case "${ALTITUDE_UNIT}" in
"m")
- cdo -setmissval,nan -select,name=${ALTITUDE_NAME} ${ALTI_FILE} ${TMP_ALT_FILE}
+ cdo -s -setmissval,nan -select,name=${ALTITUDE_NAME} ${ALTI_FILE} ${TMP_ALT_FILE}
;;
*)
log_error "ERROR: Altitude unit must be 'm'"
@@ -205,7 +238,7 @@ case "${ALTITUDE_UNIT}" in
esac
if [ ! ${REGION} == null ]; then
- cdo sellonlatbox,${REGION} ${TMP_ALT_FILE} 'regional_alt_data.nc'
+ cdo -s sellonlatbox,${REGION} ${TMP_ALT_FILE} 'regional_alt_data.nc'
mv 'regional_alt_data.nc' ${TMP_ALT_FILE}
fi
@@ -229,17 +262,20 @@ for ((YEAR=$START_YEAR;YEAR<END_YEAR+1;YEAR+=1)); do
TMP_FILE=${TMP_DIR_YEAR}'/tmp.nc'
# TMAX
- cdo -selyear,${YEAR} ${TMP_DIR}'/tasmax_'${START_YEAR}'-'${END_YEAR}'.nc' ${INPUT_FILE}
+ cdo -s -selyear,${YEAR} ${TMP_DIR}'/tasmax_'${START_YEAR}'-'${END_YEAR}'.nc' ${INPUT_FILE}
# TMIN
- cdo -selyear,${YEAR} ${TMP_DIR}'/tasmin_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
+ cdo -s -selyear,${YEAR} ${TMP_DIR}'/tasmin_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
ncks -A ${TMP_FILE} ${INPUT_FILE}
# WIND
- cdo -selyear,${YEAR} ${TMP_DIR}'/sfcWind_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
+ cdo -s -selyear,${YEAR} ${TMP_DIR}'/sfcWind_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
ncks -A ${TMP_FILE} ${INPUT_FILE}
# add Wind_height as new variable (required for log. windprofile ...m to 2 m)
TMP_FILE2=${TMP_DIR_YEAR}'/tmp2.nc'
- cdo -setmissval,nan -setattribute,${WIND_H_NAME}@units="m" -chname,sfcWind,${WIND_H_NAME} -setcindexbox,${WIND_HEIGHT},1,${X_SIZE},1,${Y_SIZE} ${TMP_FILE} ${TMP_FILE2} # windspeed -> wind_height=10m
+ WIND_HEIGHT=$(get_height ${TMP_FILE})
+ X_SIZE=$(get_x_size ${TMP_FILE})
+ Y_SIZE=$(get_y_size ${VAR_FILE_NAME})
+ cdo -s -setmissval,nan -setattribute,${WIND_H_NAME}@units="m" -chname,sfcWind,${WIND_H_NAME} -setcindexbox,${WIND_HEIGHT},1,${X_SIZE},1,${Y_SIZE} ${TMP_FILE} ${TMP_FILE2} # windspeed -> wind_height=10m
#cdo(2) setcindexbox (Warning): First latitude index out of range, set to 1!
#cdo(2) setcindexbox (Warning): Last longitude index out of range, set to 1!
@@ -248,11 +284,11 @@ for ((YEAR=$START_YEAR;YEAR<END_YEAR+1;YEAR+=1)); do
ncks -A ${TMP_FILE2} ${INPUT_FILE}
# VAP
- cdo -selyear,${YEAR} ${TMP_DIR}'/pvap_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
+ cdo -s -selyear,${YEAR} ${TMP_DIR}'/pvap_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
ncks -A ${TMP_FILE} ${INPUT_FILE}
# SOL
- cdo -selyear,${YEAR} ${TMP_DIR}'/rsds_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
+ cdo -s -selyear,${YEAR} ${TMP_DIR}'/rsds_'${START_YEAR}'-'${END_YEAR}'.nc' ${TMP_FILE}
ncks -A ${TMP_FILE} ${INPUT_FILE}
# ALTITUDE
@@ -263,13 +299,13 @@ for ((YEAR=$START_YEAR;YEAR<END_YEAR+1;YEAR+=1)); do
rm ${TMP_FILE2}
MYEXP_FILE=$(dirname "$0")/helper/myexpr
- cdo exprf,${MYEXP_FILE} ${INPUT_FILE} ${TMP_DIR_YEAR}'/res_'${YEAR}'.nc'
+ cdo -s exprf,${MYEXP_FILE} ${INPUT_FILE} ${TMP_DIR_YEAR}'/res_'${YEAR}'.nc'
if [ ! -f ${TMP_DIR_YEAR}'/res_'${YEAR}'.nc' ]; then
log_error 'Error create: '${TMP_DIR_YEAR}'/res_'${YEAR}'.nc'
exit 1
fi
- cdo -setmissval,-1.0E20 -chname,ET0,${PET_NAME} -selname,ET0 ${TMP_DIR_YEAR}'/res_'${YEAR}'.nc' ${TMP_DIR}'/'${PET_NAME}'_'${YEAR}'_tmp.nc'
+ cdo -s -setmissval,-1.0E20 -chname,ET0,${PET_NAME} -selname,ET0 ${TMP_DIR_YEAR}'/res_'${YEAR}'.nc' ${TMP_DIR}'/'${PET_NAME}'_'${YEAR}'_tmp.nc'
ncatted -O -a units,${PET_NAME},c,c,"""${PET_UNIT}""" ${TMP_DIR}'/'${PET_NAME}'_'${YEAR}'_tmp.nc'
ncatted -O -a long_name,${PET_NAME},c,c,"""${PET_LNAME}""" ${TMP_DIR}'/'${PET_NAME}'_'${YEAR}'_tmp.nc'
ncatted -O -a standard_name,${PET_NAME},c,c,"""${PET_SNAME}""" ${TMP_DIR}'/'${PET_NAME}'_'${YEAR}'_tmp.nc'
diff --git a/src/pet/helper/calc_pvap b/src/pet/helper/calc_pvap
new file mode 100644
index 0000000000000000000000000000000000000000..87ae5827dd6c786dc5a420a115596c10487566c0
--- /dev/null
+++ b/src/pet/helper/calc_pvap
@@ -0,0 +1,10 @@
+######################################################
+#
+# Calculation of daily vapor pressure
+# according to: FAO Irrigation and Drainage Paper No. 56 (Crop Evapotranspiration), Allan, Pereira 1998
+#
+######################################################
+
+### CALCULATION ###
+_pvap_sat=(0.6108*exp(17.27*tasmax/(tasmax+237.3))+0.6108*exp(17.27*tasmin/(tasmin+237.3)))/2; # (Equation 12)
+pvap=((hurs/100)*_pvap_sat); # water vapour pressure in Pa (Equation 10)
diff --git a/src/pet/helper/calc_pvap.license b/src/pet/helper/calc_pvap.license
new file mode 100644
index 0000000000000000000000000000000000000000..1e66d6f06c2798c809741da503ad3fd87439d165
--- /dev/null
+++ b/src/pet/helper/calc_pvap.license
@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2025 Brandenburgische Technische Universität Cottbus-Senftenberg
+
+SPDX-License-Identifier: BSD-3-Clause
diff --git a/src/pet/helper/myexpr b/src/pet/helper/myexpr
index 3b167e106276d54d2b35cbae008f90f89514c318..76c2ecdaa4814d142e792177f6a849fb001a8aaa 100644
--- a/src/pet/helper/myexpr
+++ b/src/pet/helper/myexpr
@@ -18,14 +18,14 @@
_var_tmax=tasmax;
_var_tmin=tasmin;
e_a=pvap;
-_var_wind=sfcWind;
R_s=rsds;
_var_altitude=orog;
_var_wind_height=Wind_height;
+_var_wind=sfcWind;
### NATURAL CONSTANTS ###
-c_p=1.013*10^(-3); # specific heat at constant pressure [MJ kg-1 °C-1]
+c_p=1.013*10^(-3); # specific heat at constant pressure [MJ kg-1 �C-1]
epsilon=0.622; # ratio molecular weight of water vapour/dry air
lambda=2.45; # latent heat of vaporization [MJ kg-1]
G_sc=0.0820; # solar constant [MJ m-2 d-1]
@@ -36,18 +36,18 @@ pi=3.14159265358979323846; # Pi
### CALCULATIONS ###
-# MEAN TEMPERATURE: daily mean temperature [°C]
+# MEAN TEMPERATURE: daily mean temperature [�C]
T_mean=(_var_tmax+_var_tmin)/2; #(Eq.9)
# VAPOR PRESSURE DEFICIT [kPa]
pressure=101.3*(((293-0.0065*_var_altitude)/293)^5.26); # atmospheric pressure at altitude z [kPa] (Eq.7) 2dim
-gamma=(c_p*pressure/(epsilon*lambda)); # psychrometric constant [kPa °C-1] (Eq.8) 3dim? daily values equal
+gamma=(c_p*pressure/(epsilon*lambda)); # psychrometric constant [kPa �C-1] (Eq.8) 3dim? daily values equal
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
##where are time und time_bnds in output file from? input file
e_sat=(0.6108*exp(17.27*_var_tmax/(_var_tmax+237.3)) + 0.6108*exp(17.27*_var_tmin/(_var_tmin+237.3)))/2; # mean saturation vapour pressure [kPa] (Eq.11,12)
e_deficit=e_sat-e_a; # vapour pressure deficit
-# SLOPE OF SATURATION VAPOUR PRESSURE CURVE [kPa °C-1]
+# SLOPE OF SATURATION VAPOUR PRESSURE CURVE [kPa �C-1]
delta=(4098*(0.6108*exp(17.27*T_mean/(T_mean+237.3))))/((T_mean+237.3)^2); #(Eq.13)
# WIND SPEED [m s-1]
@@ -61,7 +61,7 @@ R_ns=(1-a)*R_s; #(Eq.38)
# extraterrestrial radiation
lat=clat(_var_tmax); # get latitudes as variable
##cdi warning (cdfCheckVarname): Changed double entry of variable name 'lat' to 'lat_2'!
-lat_rad=pi/180*(lat); # convert latitude from [°] to [rad] (Eq.22)
+lat_rad=pi/180*(lat); # convert latitude from [�] to [rad] (Eq.22)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
#day_of_year=ctimestep(); # timestep number = day of year (for yearly files, not necassary to set with cdo: days since 1.1.) not in file
@@ -72,10 +72,10 @@ N1=floor(275 * month / 9);
N2=floor((month + 9) / 12);
N3=(1 + floor((year - 4 * floor(year / 4) + 2) / 3));
day_of_year=N1-(N2 * N3) + day - 30;
-### fehler -1 Tag für 1900,2100 (kein schaltjahr)
+### fehler -1 Tag f�r 1900,2100 (kein schaltjahr)
dist_earth_sun=1+0.033*cos((2*pi/365)*day_of_year); # inverse relative distance Earth-Sun (Eq.23) not in file
-solar_dec=0.409*sin((2*pi/365)*day_of_year-1.39); # solar declination [rad], must be -23.45° <= x <= 23.45° (Eq.24) not in file
+solar_dec=0.409*sin((2*pi/365)*day_of_year-1.39); # solar declination [rad], must be -23.45� <= x <= 23.45� (Eq.24) not in file
test_var=(-tan(lat_rad)*tan(solar_dec)); #(Eq.25)
##cdi warning (vlistDefVarTiles): Unexpected time type -1, set to TIME_VARYING!
# calculate sunset_angle only for valid values of cosine function
@@ -105,7 +105,7 @@ test_var3=(test_var2>1)+(test_var2<=1)*test_var2;
R_nl=sigma*((T_max_K^4+T_min_K^4)/2)*(0.34-0.14*sqrt(e_a))*(1.35*test_var3-0.35); #(Eq.39)
# net radiation [MJ m-2 d-1]
-R_n=R_ns-R_nl; #(Eq.40)
+R_n=R_ns-R_nl; #(Eq.4����0)
G=0; # estimated soil heat flux for daily periods [MJ m-2] (Eq.42)
# daily potential evapotranspiration ET0 [mm d-1]
diff --git a/src/pet/out.json b/src/pet/out.json
new file mode 100644
index 0000000000000000000000000000000000000000..3f82f2faeed2bbc3fd695a71cd4478176d754e09
--- /dev/null
+++ b/src/pet/out.json
@@ -0,0 +1,35 @@
+{
+ "project": "cmip6",
+ "product": "c4mip",
+ "institute": "mohc",
+ "model": "ukesm1-0-ll",
+ "experiment": "hist-bgc",
+ "ensemble": "r4i1p1f2",
+ "start_date": "1900-01-01",
+ "end_date": "1910-12-30",
+ "region": "5,15,47,55",
+ "output_dir": "/work/ch1187/regiklim-work/k202206/regiklim-ces/output/pet/6488",
+ "debug": false,
+ "files": {
+ "hurs": [
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/hurs/gn/v20190829/hurs_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_19500101-20141230.nc",
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/hurs/gn/v20190829/hurs_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_18500101-19491230.nc"
+ ],
+ "tasmax": [
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/tasmax/gn/v20190829/tasmax_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_18500101-19491230.nc",
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/tasmax/gn/v20190829/tasmax_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_19500101-20141230.nc"
+ ],
+ "tasmin": [
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/tasmin/gn/v20190829/tasmin_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_19500101-20141230.nc",
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/tasmin/gn/v20190829/tasmin_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_18500101-19491230.nc"
+ ],
+ "sfcWind": [
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/sfcWind/gn/v20190829/sfcWind_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_19500101-20141230.nc",
+ "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MOHC/UKESM1-0-LL/hist-bgc/r4i1p1f2/day/sfcWind/gn/v20190829/sfcWind_day_UKESM1-0-LL_hist-bgc_r4i1p1f2_gn_18500101-19491230.nc"
+ ],
+ "pvap": [],
+ "rsds": []
+ },
+ "alti_file": "/work/ik1017/CMIP6/data/CMIP6/C4MIP/MIROC/MIROC-ES2L/hist-bgc/r1i1p1f2/fx/orog/gn/v20191129/orog_fx_MIROC-ES2L_hist-bgc_r1i1p1f2_gn.nc",
+ "cache_dir": "/scratch/k/k202206/pet/6488"
+}
diff --git a/src/pet/pet_functions.sh b/src/pet/pet_functions.sh
index 330266330aff53a61979afa6fc1df915ec44eba0..e236e84c128ec3d122d0a923f8ac5ff9d042bb20 100755
--- a/src/pet/pet_functions.sh
+++ b/src/pet/pet_functions.sh
@@ -96,7 +96,7 @@ convert_unit() {
case ${UNIT} in
"K" | "Kelvin" | "KELVIN")
- cdo -s -setmissval,nan -setattribute,$3@units="Celsius" -addc,-237.15 $1 $1'_new.nc'
+ cdo -s -setmissval,nan -setattribute,$2@units="Celsius" -addc,-273.15 $1 $1'_new.nc'
mv $1'_new.nc' $1
;;
"°C" | "Celsius" | "CELSIUS")
@@ -104,7 +104,7 @@ convert_unit() {
mv $1'_new.nc' $1
;;
"km h-1")
- cdo -s -setmissval,nan -setattribute,$3@units="m s-1" -divc,3.6 $1 $1'_new.nc'
+ cdo -s -setmissval,nan -setattribute,$2@units="m s-1" -divc,3.6 $1 $1'_new.nc'
mv $1'_new.nc' $1
;;
"m s-1")
@@ -112,7 +112,7 @@ convert_unit() {
mv $1'_new.nc' $1
;;
"Pa")
- cdo -s -setmissval,nan -setattribute,$3@units="kPa" -divc,1000 $1 $1'_new.nc'
+ cdo -s -setmissval,nan -setattribute,$2@units="kPa" -divc,1000 $1 $1'_new.nc'
mv $1'_new.nc' $1
;;
"kPa")
@@ -120,13 +120,17 @@ convert_unit() {
mv $1'_new.nc' $1
;;
"W m-2")
- cdo -s -setmissval,nan -setattribute,$3@units="MJ m-2 d-1" -mulc,0.0864 $1 $1'_new.nc' # 60*60*24 / 1000000
+ cdo -s -setmissval,nan -setattribute,$2@units="MJ m-2 d-1" -mulc,0.0864 $1 $1'_new.nc' # 60*60*24 / 1000000
mv $1'_new.nc' $1
;;
"MJ m-2 d-1")
cdo -s -setmissval,nan $1 $1'_new.nc'
mv $1'_new.nc' $1
;;
+ "%")
+ cdo -s -setmissval,nan $1 $1'_new.nc'
+ mv $1'_new.nc' $1
+ ;;
*)
log_error "Error: The data doesnt match the given units (temperature: [°C]/[K], wind: [km h-1]/[m s-1], water vapour pressure: [Pa]/[kPa], radiation: [W m-2]/[MJ m-2 d-1])"
exit 1 ### exit /b
diff --git a/src/tests/conftest.py b/src/tests/conftest.py
index a900250fe30de6be25b505f90bcd432dc37006cb..c95e20c902f544a7a188140ec2e7db3f69160d3a 100644
--- a/src/tests/conftest.py
+++ b/src/tests/conftest.py
@@ -37,6 +37,12 @@ def update_drs_config():
}
config["nukleus-btu.defaults"] = {"project": "nukleus"}
+ config["cmip6"] = {
+ "root_path": "/home/freva/work/ch1187/freva-regiklim/data/model/global/cmip6/CMIP/output",
+ "drs_format": "cmip6",
+ }
+ config["cmip6.defaults"] = {"project": "cmip6"}
+
with open(drs_config_file, "w") as file:
toml.dump(config, file)