diff --git a/docs/INDEX.md b/docs/INDEX.md index 95b9bb7..c5a6376 100644 --- a/docs/INDEX.md +++ b/docs/INDEX.md @@ -6,6 +6,8 @@ For an overview of the needed atmospheric forcing variables for eCLM see [Overvi For the default creation of atmospheric forcing based on ERA5 data see [eCLM atmospheric forcing based on ERA5](era5forcing). +For the default creation of atmospheric forcing based on SEAS5 data see [eCLM atmospheric forcing based on SEAS5](seas5forcing). + For a Python script that validates eCLM atmospheric forcing files specified in `datm_in` namelist and stream files see [Checking Atmospheric Forcing](checkingforcing). diff --git a/docs/_toc.yml b/docs/_toc.yml index 4dfd047..ff10593 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -10,6 +10,8 @@ parts: title: Overview atmospheric forcing - file: users_guide/era5-forcing title: ERA5 atmospheric forcing + - file: users_guide/seas5-forcing + title: SEAS5 atmospheric forcing - file: users_guide/icon-forcing title: ICON output atmospheric forcing - file: users_guide/other-forcing diff --git a/docs/users_guide/era5-forcing.md b/docs/users_guide/era5-forcing.md index cea7bf8..abef0ab 100644 --- a/docs/users_guide/era5-forcing.md +++ b/docs/users_guide/era5-forcing.md @@ -45,6 +45,7 @@ Source the provided environment file source jsc.2024_Intel.sh ``` +(era5forcing-download)= ### Download of ERA5 data `download_ERA5_input.py` contains a prepared retrieval for the `cdsapi` python module. diff --git a/docs/users_guide/seas5-forcing.md b/docs/users_guide/seas5-forcing.md new file mode 100644 index 0000000..43af23d --- /dev/null +++ b/docs/users_guide/seas5-forcing.md @@ -0,0 +1,121 @@ +(seas5forcing)= +# eCLM atmospheric forcing based on SEAS5 + +Start by sourcing the provided environment file + +``` +source jsc.2024_Intel.sh +``` + +Provide a Python virtual environment including all necessary modules, +in particular `cdsapi` for downloading. It is important to load or +activate the pyvenv after sourcing the environment file, since the +environment file loads Python modules itself. + +First usage - creation of pyvenv: +``` + python -m venv python_cdsapi + source python_cdsapi/bin/activate + pip install cdsapi +``` + +Existing pyvenv - activating the pyvenv: +``` + source python_cdsapi/bin/activate +``` + +## Creation of forcing data from SEAS5 + +Creation of SEAS5 forcing files is adapted from the creation of ERA5 +forcing files. + +The folder `mkforcing/` contains the scripts that assist the SEAS5 +retrieval. + +``` + cd mkforcing +``` + +### Download of SEAS5 data + +`download_ERA5_input.py` contains a prepared retrieval for the cdsapi python module. + +More about cdsapi can be found in [Download of ERA5 +data](era5forcing-download). + +Usage: Three separate commands have to be executed, one for constant +variables (orography), one for daily variables (e.g. total +precipitation) and one for 6-hourly variables (e.g. temperature). + +For each type of output, a dedicated download directory is created. + +A main adaption from ERA5 download is the specification of a +customized CDSAPI request using `--request`. + +```bash + python download_ERA5_input.py --year --month --dirout cdsapidwn_SEAS5_const --request ../custom_request_SEAS5_const.py + python download_ERA5_input.py --year --month --dirout cdsapidwn_SEAS5_24h --request ../custom_request_SEAS5_24h.py + python download_ERA5_input.py --year --month --dirout cdsapidwn_SEAS5_06h --request ../custom_request_SEAS5_06h.py +``` + +**Note:** The wrapper script: `./download_ERA5_input_wrapper.sh` is +currently NOT SUPPORTED for SEAS5 download. + + +### Preparation of SEAS5 data: all variable to 06h + +Next, we want to have all variables available in 6-hourly interval + +- from constant: `z` has to be ported (same values as before) +- from daily: `strd`, `ssrd` (thermal, solar, each accumulated) and + `tp` (total precipitation), these values would have to be + distributed. For `tp` the value would be divided by four and + assigned to the four corresponding 6-hour intervals. For the + radiations, it would be zero in the first/last interval ("night") + and half the value in the middle two intervals, then converted + to flux (W/m²) by dividing by the time interval. + +The script outputs radiation directly as flux variables (`flds`, +`fsds`) in W/m². + +```bash + python seas5_daily_to_6hourly.py --const cdsapidwn_SEAS5_const/download_era5_2026_01.nc --daily cdsapidwn_SEAS5_24h/download_era5_2026_01.nc --hourly cdsapidwn_SEAS5_06h/download_era5_2026_01.nc --output cdsapidwn_SEAS5/download_era5_2026_01.nc --frequency 3 --include-hour-zero +``` + +### Preparation of SEAS5 data: correct input variables + +Steps for preparing SEAS5 data as eCLM input data + +1. Orography to elevation (adds a variable `elevation` to the netCDF + file) + +```bash + python orography_to_elevation.py cdsapidwn_SEAS5/download_era5_2026_01.nc +``` + +2. Mean sea level pressure to surface pressure + +```bash + python mslp_to_sp.py cdsapidwn_SEAS5/download_era5_2026_01.nc --elevation-var elevation +``` + +3. Humidity computed from dewpoint temperature and surface pressure + +``` + python dewpoint_to_specific_humidity.py cdsapidwn_SEAS5/download_era5_2026_01.nc +``` + +4. Temperature and Specific Humidity converted from 2m to 10m + +``` + python 2m_to_10m_conversion.py cdsapidwn_SEAS5/download_era5_2026_01.nc +``` + +### Preparation of SEAS5 data: Remapping, Data merging, CLM3.5 + + +Check inputs and replace according to your case. + +``` + sh prepare_SEAS5_input.sh lwgtdis=true lgriddes=true wgtcaf=../wgtdis_era5caf_to_DE-RuS-$(date +%y%m%d).nc griddesfile=../griddes_DE-RuS_$(date +%y%m%d).txt iyear=2026 imonth=01 author="Johannes KELLER" email="jo.keller@fz-juelich.de" tmpdir=tmpdir pathdata=../cdsapidwn_SEAS5 +``` diff --git a/mkforcing/custom_request_SEAS5_06h.py b/mkforcing/custom_request_SEAS5_06h.py new file mode 100644 index 0000000..8aa6d4b --- /dev/null +++ b/mkforcing/custom_request_SEAS5_06h.py @@ -0,0 +1,38 @@ +## Download for SEAS5 forecast variables +# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download + +# 6-HOURLY VARIABLES + +# Forecast for 7 months + +# import cdsapi + +dataset = "seasonal-original-single-levels" +request = { + "originating_centre": "ecmwf", + "system": "51", + "variable": [ + # # constant + # "orography", # used to convert mslp to sp + # # 24-hourly + # "surface_thermal_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + # "surface_solar_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + # "total_precipitation", + # 6-hourly + "mean_sea_level_pressure", # convert to surface pressure, use elevation (hypsometric formula), surface geopotential height (orography) + "10m_u_component_of_wind", + "10m_v_component_of_wind", + "2m_temperature", + "2m_dewpoint_temperature", + ], + "year": ["2025"], + "month": ["09"], + "day": ["01"], + "leadtime_hour": [str(h) for h in range(0, 5161, 6)], + "data_format": "netcdf", + "area": [50.870906+0.5, 6.4421445-0.5, 50.870906-0.5, 6.4421445+0.5] # Selhausen + # "area": [74, -42, 20, 69] # Europe +} + +# client = cdsapi.Client() +# client.retrieve(dataset, request).download() diff --git a/mkforcing/custom_request_SEAS5_24h.py b/mkforcing/custom_request_SEAS5_24h.py new file mode 100644 index 0000000..5a4b74e --- /dev/null +++ b/mkforcing/custom_request_SEAS5_24h.py @@ -0,0 +1,38 @@ +## Download for SEAS5 forecast variables +# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download + +# 24-HOURLY / DAILY VARIABLES + +# Forecast for 7 months + +# import cdsapi + +dataset = "seasonal-original-single-levels" +request = { + "originating_centre": "ecmwf", + "system": "51", + "variable": [ + # # constant + # "orography", # used to convert mslp to sp + # 24-hourly + "surface_thermal_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + "surface_solar_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + "total_precipitation", + # # 6-hourly + # "mean_sea_level_pressure", # convert to surface pressure, use elevation (hypsometric formula), surface geopotential height (orography) + # "10m_u_component_of_wind", + # "10m_v_component_of_wind", + # "2m_temperature", + # "2m_dewpoint_temperature", + ], + "year": ["2025"], + "month": ["09"], + "day": ["01"], + "leadtime_hour": [str(h) for h in range(0, 5161, 6)], + "data_format": "netcdf", + "area": [50.870906+0.5, 6.4421445-0.5, 50.870906-0.5, 6.4421445+0.5] # Selhausen + # "area": [74, -42, 20, 69] # Europe +} + +# client = cdsapi.Client() +# client.retrieve(dataset, request).download() diff --git a/mkforcing/custom_request_SEAS5_const.py b/mkforcing/custom_request_SEAS5_const.py new file mode 100644 index 0000000..c6203bf --- /dev/null +++ b/mkforcing/custom_request_SEAS5_const.py @@ -0,0 +1,38 @@ +## Download for SEAS5 forecast variables +# https://cds.climate.copernicus.eu/datasets/seasonal-original-single-levels?tab=download + +# CONSTANT VARIABLES + +# Forecast for 7 months + +# import cdsapi + +dataset = "seasonal-original-single-levels" +request = { + "originating_centre": "ecmwf", + "system": "51", + "variable": [ + # constant + "orography", # used to convert mslp to sp + # # 24-hourly + # "surface_thermal_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + # "surface_solar_radiation_downwards", # Unit conversion from accumulated value [J/m2] to mean rate [W/m2] + # "total_precipitation", + # # 6-hourly + # "mean_sea_level_pressure", # convert to surface pressure, use elevation (hypsometric formula), surface geopotential height (orography) + # "10m_u_component_of_wind", + # "10m_v_component_of_wind", + # "2m_temperature", + # "2m_dewpoint_temperature", + ], + "year": ["2025"], + "month": ["09"], + "day": ["01"], + "leadtime_hour": [str(h) for h in range(0, 5161, 6)], + "data_format": "netcdf", + "area": [50.870906+0.5, 6.4421445-0.5, 50.870906-0.5, 6.4421445+0.5] # Selhausen + # "area": [74, -42, 20, 69] # Europe +} + +# client = cdsapi.Client() +# client.retrieve(dataset, request).download() diff --git a/mkforcing/download_ERA5_input.py b/mkforcing/download_ERA5_input.py index 6bf1467..cdfbfbc 100755 --- a/mkforcing/download_ERA5_input.py +++ b/mkforcing/download_ERA5_input.py @@ -139,7 +139,11 @@ def generate_datarequest(year, monthstr, days, # Adapt year, month and day to input values request["year"] = [str(year)] request["month"] = [monthstr] - request["day"] = days + if dataset == "seasonal-original-single-levels": + # First day of month specified for SEAS5 + request["day"] = ["01"] + else: + request["day"] = days # Temporary filename w/o extension if not provided auto_detect_extension = target is None @@ -159,7 +163,7 @@ def generate_datarequest(year, monthstr, days, extension = detect_file_type(target) # Rename to final target with correct extension - final_target = f'{target}{extension}' + final_target = f'download_era5_{year}_{monthstr}{extension}' os.rename(target, final_target) target = final_target @@ -230,6 +234,8 @@ def generate_datarequest(year, monthstr, days, custom_request = None custom_dataset = args.dataset if args.request: + if not os.path.isfile(args.request): + raise FileNotFoundError(f"Custom request file not found: {args.request}") import importlib.util spec = importlib.util.spec_from_file_location("custom_request_module", args.request) custom_module = importlib.util.module_from_spec(spec) diff --git a/mkforcing/mslp_to_sp.py b/mkforcing/mslp_to_sp.py new file mode 100644 index 0000000..fcdfa89 --- /dev/null +++ b/mkforcing/mslp_to_sp.py @@ -0,0 +1,535 @@ +""" +Convert between Mean Sea Level Pressure (MSLP) and Surface Pressure +Using the Method from Stull's Practical Meteorology Chapter 9 + +This implementation follows the standard meteorological sea-level pressure +reduction method that accounts for the elevation and uses a fictitious +temperature for the imaginary air column between the surface and sea level. + +References (all open-access): +1. Stull, R., 2017: Practical Meteorology: An Algebra-based Survey of + Atmospheric Science. University of British Columbia, 940 pp. + ISBN 978-0-88865-283-6 + Available under Creative Commons License (CC BY-NC-SA 4.0) + URL: https://www.eoas.ubc.ca/books/Practical_Meteorology/ + See Chapter 9, Section "Sea-level Pressure Reduction" + +Mathematical Background: +The sea-level pressure reduction uses the hypsometric equation: + + P_MSL = P_surface * exp(z_stn / (a * T_v*)) + +where: +- P_MSL: Mean sea level pressure (Pa or hPa) +- P_surface: Surface pressure at the station (Pa or hPa) +- z_stn: Station elevation above sea level (m) +- a: Scale height parameter = R_d / g ≈ 29.3 m/K +- T_v*: Fictitious average virtual temperature (K) + +The fictitious temperature T_v* represents the temperature of an imaginary +air column between the surface and sea level. It is calculated as: + + T_v* = 0.5 * [T_v(t_0) + T_v(t_0 - 12h) + γ_sa * z_stn] + +where: +- T_v(t_0): Current virtual temperature at the surface +- T_v(t_0 - 12h): Virtual temperature 12 hours ago +- γ_sa: Standard atmosphere lapse rate = 0.0065 K/m +""" + +import argparse +import numpy as np +import netCDF4 + + +def mslp_to_surface_pressure( + mslp: float, + elevation: float, + temperature_current: float, + temperature_12h_ago: float = None, + pressure_units: str = "Pa", +) -> float: + """ + Convert mean sea level pressure to surface pressure at a given elevation. + + Uses the standard meteorological method with a fictitious temperature + for the imaginary air column between the surface and sea level. + + Parameters + ---------- + mslp : float + Mean sea level pressure (in units specified by pressure_units) + elevation : float + Station elevation above mean sea level (meters) + - Positive for locations above sea level + - Negative for locations below sea level + temperature_current : float + Current surface air temperature (Kelvin) + temperature_12h_ago : float, optional + Surface air temperature 12 hours ago (Kelvin) + If None, uses temperature_current (i.e., assumes steady conditions) + pressure_units : str, optional + Units for pressure ('Pa' or 'hPa'), default='Pa' + + Returns + ------- + float + Surface pressure at the given elevation (same units as input) + + Notes + ----- + 1. This follows the method described in Stull (2017), Chapter 9. + 2. The fictitious temperature accounts for: + - Current surface temperature + - Temperature 12 hours ago (for diurnal averaging) + - Lapse rate correction for the imaginary column below the surface + 3. For locations at sea level (elevation=0), returns mslp unchanged. + 4. The 12-hour averaging helps smooth out diurnal temperature variations. + + Examples + -------- + >>> # Sea level location + >>> surface_p = mslp_to_surface_pressure(101325, 0, 288.15, None, 'Pa') + >>> print(f"Surface pressure: {surface_p:.1f} Pa") + Surface pressure: 101325.0 Pa + + >>> # Mountain station (2000m) with temperature data + >>> surface_p = mslp_to_surface_pressure(1013.25, 2000, 278, 276, 'hPa') + >>> print(f"Surface pressure: {surface_p:.1f} hPa") + Surface pressure: 795.3 hPa + """ + + # Physical constants + R_d = 287.05 # Gas constant for dry air (J/(kg·K)) + g = 9.80665 # Gravitational acceleration (m/s²) + gamma_sa = 0.0065 # Standard atmosphere lapse rate (K/m) + + # Calculate scale height parameter + a = R_d / g # ≈ 29.3 m/K + + # If no 12-hour-ago temperature provided, use current temperature + if temperature_12h_ago is None: + temperature_12h_ago = temperature_current + + # Calculate fictitious average virtual temperature for the column + # This is the temperature of the imaginary air column between + # the surface and sea level + # T_v* = 0.5 * [T_v(now) + T_v(12h ago) + γ_sa * z_stn] + T_v_star = 0.5 * (temperature_current + temperature_12h_ago + gamma_sa * elevation) + + # Apply hypsometric equation (inverse direction) + # P_surface = P_MSL * exp(-z / (a * T_v*)) + exponent = -elevation / (a * T_v_star) + surface_pressure = mslp * np.exp(exponent) + + return surface_pressure + + +def surface_to_mslp( + surface_pressure: float, + elevation: float, + temperature_current: float, + temperature_12h_ago: float = None, + pressure_units: str = "Pa", +) -> float: + """ + Convert surface pressure to mean sea level pressure. + + This is the standard meteorological "sea-level pressure reduction" + operation used to create weather maps. + + Parameters + ---------- + surface_pressure : float + Pressure at the surface (in units specified by pressure_units) + elevation : float + Station elevation above mean sea level (meters) + temperature_current : float + Current surface air temperature (Kelvin) + temperature_12h_ago : float, optional + Surface air temperature 12 hours ago (Kelvin) + If None, uses temperature_current + pressure_units : str, optional + Units for pressure ('Pa' or 'hPa'), default='Pa' + + Returns + ------- + float + Mean sea level pressure (same units as input) + + Notes + ----- + This function uses the fictitious temperature method from Stull (2017): + + T_v* = 0.5 * [T_v(now) + T_v(12h ago) + γ_sa * z_stn] + + Then applies: + P_MSL = P_surface * exp(z / (a * T_v*)) + + Examples + -------- + >>> # Convert Denver surface pressure to MSLP + >>> mslp = surface_to_mslp(83500, 1600, 288.15, 285.15, 'Pa') + >>> print(f"MSLP: {mslp:.1f} Pa ({mslp/100:.1f} hPa)") + MSLP: 101362.4 Pa (1013.6 hPa) + + >>> # Mountain weather station + >>> mslp = surface_to_mslp(795.0, 2000, 278, 276, 'hPa') + >>> print(f"MSLP: {mslp:.2f} hPa") + MSLP: 1012.85 hPa + """ + + # Physical constants + R_d = 287.05 # Gas constant for dry air (J/(kg·K)) + g = 9.80665 # Gravitational acceleration (m/s²) + gamma_sa = 0.0065 # Standard atmosphere lapse rate (K/m) + + # Calculate scale height parameter + a = R_d / g # ≈ 29.3 m/K + + # If no 12-hour-ago temperature provided, use current temperature + if temperature_12h_ago is None: + temperature_12h_ago = temperature_current + + # Calculate fictitious average temperature for the imaginary column + T_v_star = 0.5 * (temperature_current + temperature_12h_ago + gamma_sa * elevation) + + # Apply hypsometric equation + # P_MSL = P_surface * exp(z / (a * T_v*)) + exponent = elevation / (a * T_v_star) + mslp = surface_pressure * np.exp(exponent) + + return mslp + + +def calculate_fictitious_temperature( + temperature_current: float, temperature_12h_ago: float, elevation: float +) -> float: + """ + Calculate the fictitious temperature for the imaginary air column. + + This is a helper function that computes T_v* according to the method + in Stull (2017), Chapter 9. + + Parameters + ---------- + temperature_current : float + Current surface temperature (Kelvin) + temperature_12h_ago : float + Temperature 12 hours ago (Kelvin) + elevation : float + Station elevation (meters) + + Returns + ------- + float + Fictitious temperature T_v* (Kelvin) + + Notes + ----- + The formula is: + T_v* = 0.5 * [T(now) + T(12h ago) + γ_sa * z] + + This represents: + 1. Average of current and 12-hour-ago temperatures (diurnal smoothing) + 2. Plus a lapse-rate correction for the imaginary column height + + Examples + -------- + >>> T_star = calculate_fictitious_temperature(288.15, 285.15, 1600) + >>> print(f"Fictitious temperature: {T_star:.2f} K ({T_star-273.15:.2f}°C)") + Fictitious temperature: 291.55 K (18.40°C) + """ + gamma_sa = 0.0065 # Standard atmosphere lapse rate (K/m) + + T_v_star = 0.5 * (temperature_current + temperature_12h_ago + gamma_sa * elevation) + + return T_v_star + + +def add_surface_pressure_to_netcdf(filename, elevation_var="z", temp_var="t2m", mslp_var="msl"): + """ + Read MSLP, temperature, and elevation from a netCDF file, + calculate surface pressure, and write it back to the file. + + Parameters: + ----------- + filename : str + Path to the netCDF file + elevation_var : str, optional + Name of the elevation variable in the file (default: 'z') + The geopotential variable will be converted to elevation (m) by dividing by g + temp_var : str, optional + Name of the temperature variable (default: 't2m') + mslp_var : str, optional + Name of the MSLP variable (default: 'msl') + + Returns: + -------- + None + Modifies the netCDF file in place by adding 'sp' variable + + Raises: + ------- + ValueError + If 'sp' variable already exists in the file + KeyError + If required variables are not found + """ + + # Open netCDF file in append mode + print(f"Opening {filename}...") + nc = netCDF4.Dataset(filename, "a") + + try: + # Check if sp already exists - if so, raise error and exit + if "sp" in nc.variables: + nc.close() + raise ValueError( + f"Variable 'sp' already exists in {filename}. " + "No changes made. Delete the variable first if you want to recalculate." + ) + + # Read the required variables + print(f"Reading {mslp_var}, {temp_var}, and {elevation_var}...") + msl = nc.variables[mslp_var][:] # Mean sea level pressure [Pa] + t2m = nc.variables[temp_var][:] # Temperature at 2m [K] + z_geopotential = nc.variables[elevation_var][:] # Geopotential [m^2/s^2] + + print(f"Data shapes - {mslp_var}: {msl.shape}, {temp_var}: {t2m.shape}, {elevation_var}: {z_geopotential.shape}") + + # Convert geopotential to elevation (meters) + g = 9.80665 # Standard gravity [m/s^2] + elevation = z_geopotential / g + + print(f"Elevation range: {np.min(elevation):.1f} to {np.max(elevation):.1f} m") + + # Determine time dimension for 12h offset + # Check for SEAS5 structure (number, forecast_reference_time, forecast_period, lat, lon) + # vs ERA5 structure (time, lat, lon) + dim_names = nc.variables[temp_var].dimensions + print(f"Dimension names: {dim_names}") + + # Find the time dimension - prefer 'forecast_period' for SEAS5, otherwise first dim + if "forecast_period" in dim_names: + time_dim = "forecast_period" + time_axis = dim_names.index("forecast_period") + elif "time" in dim_names: + time_dim = "time" + time_axis = dim_names.index("time") + else: + # Fall back to first dimension + time_dim = dim_names[0] + time_axis = 0 + + time_size = nc.dimensions[time_dim].size + print(f"Time dimension: {time_dim} (axis {time_axis}) with size {time_size}") + + # Calculate surface pressure + print("Calculating surface pressure...") + + # For SEAS5 data with forecast_period as time dimension, we need to + # process along the forecast_period axis while keeping other dimensions + # The mslp_to_surface_pressure function works element-wise via numpy + + # Create array for surface pressure + sp = np.zeros_like(msl) + + # For SEAS5 structure: (number, forecast_reference_time, forecast_period, lat, lon) + # time_axis = 2, so we need to slice along axis 2 + if time_axis == 2: + # SEAS5 structure + print("Detected SEAS5 data structure (ensemble data)...") + + # For the first timestep (index 0), we don't have 12h ago data + if time_size > 0: + print("Processing first timestep (no 12h-ago data available)...") + sp[:, :, 0, :, :] = mslp_to_surface_pressure( + msl[:, :, 0, :, :], + elevation[:, :, 0, :, :], + t2m[:, :, 0, :, :], + t2m[:, :, 0, :, :], # Use current temp as 12h-ago temp + "Pa" + ) + + # For 6-hourly data, 12h = 2 timesteps back + # For hourly data, 12h = 12 timesteps back + # Assume 6-hourly if time_size <= 8 per day, otherwise hourly + timesteps_12h = 2 if time_size <= 8 else 12 + + if time_size > timesteps_12h: + print(f"Processing remaining timesteps with 12h-ago data (offset={timesteps_12h})...") + for t in range(timesteps_12h, time_size): + sp[:, :, t, :, :] = mslp_to_surface_pressure( + msl[:, :, t, :, :], + elevation[:, :, t, :, :], + t2m[:, :, t, :, :], + t2m[:, :, t - timesteps_12h, :, :], + "Pa" + ) + + # For timesteps 1 to timesteps_12h-1, use current temp + if time_size > 1: + print(f"Processing timesteps 1-{min(timesteps_12h, time_size)-1} without 12h-ago data...") + for t in range(1, min(timesteps_12h, time_size)): + sp[:, :, t, :, :] = mslp_to_surface_pressure( + msl[:, :, t, :, :], + elevation[:, :, t, :, :], + t2m[:, :, t, :, :], + t2m[:, :, t, :, :], + "Pa" + ) + elif time_size > 1: + print("Processing all timesteps without 12h-ago data...") + for t in range(1, time_size): + sp[:, :, t, :, :] = mslp_to_surface_pressure( + msl[:, :, t, :, :], + elevation[:, :, t, :, :], + t2m[:, :, t, :, :], + t2m[:, :, t, :, :], + "Pa" + ) + else: + # ERA5 structure - time is first dimension + print("Detected ERA5 data structure...") + + # For the first timestep (index 0), we don't have 12h ago data + if time_size > 0: + print("Processing first timestep (no 12h-ago data available)...") + sp[0, ...] = mslp_to_surface_pressure( + msl[0, ...], + elevation[0, ...] if elevation.shape[0] == time_size else elevation, + t2m[0, ...], + t2m[0, ...], # Use current temp as 12h-ago temp + "Pa" + ) + + # For remaining timesteps, check if we can use 12h offset + # Assuming hourly data, 12h offset = 12 timesteps back + if time_size > 12: + print("Processing remaining timesteps with 12h-ago data...") + for t in range(12, time_size): + elev_t = elevation[t, ...] if elevation.shape[0] == time_size else elevation + sp[t, ...] = mslp_to_surface_pressure( + msl[t, ...], + elev_t, + t2m[t, ...], + t2m[t - 12, ...], + "Pa" + ) + + # For timesteps 1-11, use current temp + if time_size > 1: + print("Processing timesteps 1-11 without 12h-ago data...") + for t in range(1, min(12, time_size)): + elev_t = elevation[t, ...] if elevation.shape[0] == time_size else elevation + sp[t, ...] = mslp_to_surface_pressure( + msl[t, ...], + elev_t, + t2m[t, ...], + t2m[t, ...], + "Pa" + ) + elif time_size > 1: + print("Processing all timesteps without 12h-ago data...") + for t in range(1, time_size): + elev_t = elevation[t, ...] if elevation.shape[0] == time_size else elevation + sp[t, ...] = mslp_to_surface_pressure( + msl[t, ...], + elev_t, + t2m[t, ...], + t2m[t, ...], + "Pa" + ) + + # Create the new variable + print("Creating new variable 'sp'...") + + # Get dimension names from msl + dim_names = nc.variables[mslp_var].dimensions + + # Create the sp variable with same dimensions as msl + sp_var = nc.createVariable("sp", "f4", dim_names, zlib=True, complevel=4) + + # Add attributes + sp_var.units = "Pa" + sp_var.long_name = "Surface pressure" + sp_var.standard_name = "surface_air_pressure" + sp_var.description = ( + f"Calculated from {mslp_var}, {temp_var}, and {elevation_var} " + "using Stull (2017) Chapter 9 method" + ) + + # Write the data + sp_var[:] = sp + + print("Variable 'sp' created and written successfully!") + + # Print some statistics + print("\nStatistics:") + print(f" Surface pressure range: {np.min(sp):.1f} to {np.max(sp):.1f} Pa") + print(f" Surface pressure mean: {np.mean(sp):.1f} Pa ({np.mean(sp)/100:.2f} hPa)") + print(f" MSLP mean: {np.mean(msl):.1f} Pa ({np.mean(msl)/100:.2f} hPa)") + print(f" Mean difference: {np.mean(msl - sp):.1f} Pa ({np.mean(msl - sp)/100:.2f} hPa)") + + except KeyError as e: + print(f"Error: Required variable not found in netCDF file: {e}") + print(f"Available variables: {list(nc.variables.keys())}") + nc.close() + raise + + except ValueError as e: + # This catches the "sp already exists" error + print(f"Error: {e}") + raise + + except Exception as e: + print(f"Unexpected error occurred: {e}") + nc.close() + raise + + else: + # Only executes if no exception was raised + nc.close() + print(f"\nFile {filename} closed successfully.") + + +if __name__ == "__main__": + # Set up argument parser + parser = argparse.ArgumentParser( + description="Add surface pressure (sp) to a netCDF file based on MSLP, temperature, and elevation." + ) + parser.add_argument( + "filename", + type=str, + help="Path to the netCDF file containing MSLP, temperature, and elevation variables" + ) + parser.add_argument( + "--elevation-var", + type=str, + default="z", + help="Name of the elevation/geopotential variable (default: z)" + ) + parser.add_argument( + "--temp-var", + type=str, + default="t2m", + help="Name of the temperature variable (default: t2m)" + ) + parser.add_argument( + "--mslp-var", + type=str, + default="msl", + help="Name of the MSLP variable (default: msl)" + ) + + # Parse command-line arguments + args = parser.parse_args() + + # Process the file + add_surface_pressure_to_netcdf( + args.filename, + elevation_var=args.elevation_var, + temp_var=args.temp_var, + mslp_var=args.mslp_var + ) diff --git a/mkforcing/orography_to_elevation.py b/mkforcing/orography_to_elevation.py new file mode 100644 index 0000000..210435c --- /dev/null +++ b/mkforcing/orography_to_elevation.py @@ -0,0 +1,249 @@ +""" +Convert geopotential (orography) to elevation above mean sea level. + +Geopotential represents the gravitational potential energy per unit mass +at a given height. It is related to elevation through the gravitational +acceleration constant. + +Mathematical Background: +The relationship between geopotential and geometric height is: + + z = Φ / g + +where: +- z: Geometric height (elevation) above mean sea level (m) +- Φ: Geopotential (m²/s²) +- g: Standard gravitational acceleration = 9.80665 m/s² + +This conversion assumes a constant gravitational acceleration, which is +appropriate for typical atmospheric applications at Earth's surface. + +References: +1. WMO Guide to Meteorological Instruments and Methods of Observation + (WMO-No. 8), 2018 edition +2. ECMWF IFS Documentation - Part III: Dynamics and Numerical Procedures +""" + +import argparse +import numpy as np +import netCDF4 + + +def geopotential_to_elevation(geopotential, g=9.80665): + """ + Convert geopotential to elevation above mean sea level. + + Parameters + ---------- + geopotential : float or array + Geopotential [m²/s²] + g : float, optional + Standard gravitational acceleration [m/s²] (default: 9.80665) + + Returns + ------- + elevation : float or array + Elevation above mean sea level [m] + + Notes + ----- + The standard gravitational acceleration g = 9.80665 m/s² is defined + by ISO 80000-3:2006 and is used in meteorological applications. + + Examples + -------- + >>> # Sea level + >>> elevation = geopotential_to_elevation(0.0) + >>> print(f"Elevation: {elevation:.1f} m") + Elevation: 0.0 m + + >>> # Typical mountain height + >>> geopotential = 19613.3 # m²/s² + >>> elevation = geopotential_to_elevation(geopotential) + >>> print(f"Elevation: {elevation:.1f} m") + Elevation: 2000.0 m + """ + elevation = geopotential / g + return elevation + + +def elevation_to_geopotential(elevation, g=9.80665): + """ + Convert elevation to geopotential. + + Parameters + ---------- + elevation : float or array + Elevation above mean sea level [m] + g : float, optional + Standard gravitational acceleration [m/s²] (default: 9.80665) + + Returns + ------- + geopotential : float or array + Geopotential [m²/s²] + + Notes + ----- + This is the inverse operation of geopotential_to_elevation. + + Examples + -------- + >>> # Mount Everest + >>> geopotential = elevation_to_geopotential(8849) + >>> print(f"Geopotential: {geopotential:.1f} m²/s²") + Geopotential: 86763.3 m²/s² + + >>> # Below sea level (Dead Sea) + >>> geopotential = elevation_to_geopotential(-430) + >>> print(f"Geopotential: {geopotential:.1f} m²/s²") + Geopotential: -4216.6 m²/s² + """ + geopotential = elevation * g + return geopotential + + +def add_elevation_to_netcdf(filename, geopotential_var="z", elevation_var_name="elevation"): + """ + Read geopotential (orography) from a netCDF file, + calculate elevation, and write it back to the file. + + Parameters: + ----------- + filename : str + Path to the netCDF file + geopotential_var : str, optional + Name of the geopotential variable in the file (default: 'z') + elevation_var_name : str, optional + Name for the output elevation variable (default: 'elevation') + + Returns: + -------- + None + Modifies the netCDF file in place by adding elevation variable + + Raises: + ------- + ValueError + If elevation variable already exists in the file + KeyError + If required geopotential variable is not found + """ + + # Open netCDF file in append mode + print(f"Opening {filename}...") + nc = netCDF4.Dataset(filename, "a") + + try: + # Check if elevation variable already exists - if so, raise error and exit + if elevation_var_name in nc.variables: + nc.close() + raise ValueError( + f"Variable '{elevation_var_name}' already exists in {filename}. " + "No changes made. Delete the variable first if you want to recalculate." + ) + + # Read the geopotential variable + print(f"Reading geopotential variable '{geopotential_var}'...") + geopotential = nc.variables[geopotential_var][:] # Geopotential [m²/s²] + + print(f"Data shape - {geopotential_var}: {geopotential.shape}") + print(f"Geopotential range: {np.min(geopotential):.1f} to {np.max(geopotential):.1f} m²/s²") + + # Calculate elevation + print("Calculating elevation...") + elevation = geopotential_to_elevation(geopotential) + + print(f"Elevation range: {np.min(elevation):.1f} to {np.max(elevation):.1f} m") + + # Create the new variable + print(f"Creating new variable '{elevation_var_name}'...") + + # Get dimension names from geopotential variable + dim_names = nc.variables[geopotential_var].dimensions + + # Create the elevation variable with same dimensions as geopotential + elevation_var = nc.createVariable( + elevation_var_name, "f4", dim_names, zlib=True, complevel=4 + ) + + # Add attributes + elevation_var.units = "m" + elevation_var.long_name = "Elevation above mean sea level" + elevation_var.standard_name = "surface_altitude" + elevation_var.description = ( + f"Calculated from geopotential variable '{geopotential_var}' " + "using z = Φ / g with g = 9.80665 m/s²" + ) + + # Write the data + elevation_var[:] = elevation + + print(f"Variable '{elevation_var_name}' created and written successfully!") + + # Print some statistics + print("\nStatistics:") + print(f" Elevation range: {np.min(elevation):.1f} to {np.max(elevation):.1f} m") + print(f" Elevation mean: {np.mean(elevation):.1f} m") + print(f" Elevation std: {np.std(elevation):.1f} m") + + # Identify any interesting features + if np.min(elevation) < 0: + print(f" Below sea level: Yes (min: {np.min(elevation):.1f} m)") + if np.max(elevation) > 1000: + print(f" High elevation areas: Yes (max: {np.max(elevation):.1f} m)") + + except KeyError as e: + print(f"Error: Required variable not found in netCDF file: {e}") + print(f"Available variables: {list(nc.variables.keys())}") + nc.close() + raise + + except ValueError as e: + # This catches the "elevation already exists" error + print(f"Error: {e}") + raise + + except Exception as e: + print(f"Unexpected error occurred: {e}") + nc.close() + raise + + else: + # Only executes if no exception was raised + nc.close() + print(f"\nFile {filename} closed successfully.") + + +if __name__ == "__main__": + # Set up argument parser + parser = argparse.ArgumentParser( + description="Add elevation variable to a netCDF file based on geopotential (orography)." + ) + parser.add_argument( + "filename", + type=str, + help="Path to the netCDF file containing geopotential/orography variable" + ) + parser.add_argument( + "--geopotential-var", + type=str, + default="z", + help="Name of the geopotential variable (default: z)" + ) + parser.add_argument( + "--elevation-var", + type=str, + default="elevation", + help="Name for the output elevation variable (default: elevation)" + ) + + # Parse command-line arguments + args = parser.parse_args() + + # Process the file + add_elevation_to_netcdf( + args.filename, + geopotential_var=args.geopotential_var, + elevation_var_name=args.elevation_var + ) diff --git a/mkforcing/prepare_SEAS5_input.sh b/mkforcing/prepare_SEAS5_input.sh new file mode 100755 index 0000000..81a5a18 --- /dev/null +++ b/mkforcing/prepare_SEAS5_input.sh @@ -0,0 +1,162 @@ +#!/usr/bin/env bash +set -eo pipefail + +# default values of parameters +lrmp=true +lmerge=true +lwgtdis=false # Switch for creating wgtdis file in script +lgriddes=false # Switch for creating griddes file in script +ompthd=1 + +# TSMP2/eclm +pathdata=./ +domainfile=../domain.lnd.DE-RuS_DE-RuS.250926.nc +wgtcaf=/p/scratch/cslts/poll1/sim/euro-cordex/tsmp2_wfe_eur-11u/dta/rmp_gridwgts/wgtdis_era5caf_to_eur11u-189976.nc +# wgtmeteo=/p/scratch/cslts/poll1/sim/euro-cordex/tsmp2_wfe_eur-11u/dta/rmp_gridwgts/wgtdis_era5meteo_to_eur11u-189976.nc +griddesfile=/p/scratch/cslts/poll1/sim/euro-cordex/tsmp2_wfe_eur-11u/dta/rmp_gridwgts/griddes_eur-11u_189976.txt + +iyear=2017 +imonth=07 +tmpdir=tmpdir +wrkdir="" +author="Default AUTHOR" +email="d.fault@fz-juelich.de" +nens=51 # Number of ensemble members to process + +# Function to parse input +parse_arguments() { + for arg in "$@"; do + key="${arg%%=*}" + value="${arg#*=}" + + case "$key" in + lrmp) lrmp="$value" ;; + lmerge) lmerge="$value" ;; + lwgtdis) lwgtdis="$value" ;; + lgriddes) lgriddes="$value" ;; + ompthd) ompthd="$value" ;; + pathdata) pathdata="$value" ;; + wgtcaf) wgtcaf="$value" ;; + domainfile) domainfile="$value" ;; + # wgtmeteo) wgtmeteo="$value" ;; + griddesfile) griddesfile="$value" ;; + tmpdir) tmpdir="$value" ;; + wrkdir) wrkdir="$value" ;; + imonth) imonth="$value" ;; + iyear) iyear="$value" ;; + author) author="$value" ;; + email) email="$value" ;; + nens) nens="$value" ;; + *) echo "Warning: Unknown parameter: $key" ;; + esac + done +} + +# Call the function to parse the input arguments +# Users needs to make sure for consistent input +parse_arguments "$@" + +# +#cd $wrkdir +#mkdir -pv $tmpdir + +# +for year in ${iyear} +do +for month in ${imonth} +do + + # Go into working directory and create temporary directory + if [ -z ${wrkdir} ];then + wrkdir=${iyear}-${imonth} + fi + cd $wrkdir + mkdir -pv $tmpdir + + if $lrmp; then + + # Copy netCDF file to a template for remapping (preserves original with all ensemble members) + cp ${pathdata}/download_era5_${year}_${month}.nc ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + + # Extract the first ensemble member for the template + # The original file with all 51 ensemble members remains in pathdata for later processing + ncks --overwrite -d number,0 -O ${tmpdir}/download_era5_${year}_${month}_4rmp.nc ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + # Remove number and forecast_reference_time dimensions + ncwa --overwrite -a forecast_reference_time ${tmpdir}/download_era5_${year}_${month}_4rmp.nc ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + ncwa --overwrite -a number ${tmpdir}/download_era5_${year}_${month}_4rmp.nc ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + + # 1) Renaming variable 'valid_time' to 'time' in $file + # 2) Renaming dimension 'forecast_period' to 'valid_time' in $file + # Background: With this naming scheme the CDO remap command below will generate a time variable + # called "time" in "seconds since 1970-01-01" as for ERA5. + ncrename -v valid_time,time ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + ncrename -d forecast_period,valid_time ${tmpdir}/download_era5_${year}_${month}_4rmp.nc + + if $lwgtdis; then + cdo gendis,${domainfile} ${tmpdir}/download_era5_${year}_${month}_4rmp.nc ${wgtcaf} + fi + + if $lgriddes; then + cdo griddes ${domainfile} > ${griddesfile} + fi + + # TODO: Look into remapping! Now skipped here, but done in the loop below + # cdo -P ${ompthd} remap,${griddesfile},${wgtcaf} ${tmpdir}/download_era5_${year}_${month}_4rmp.nc ${tmpdir}/rmp_era5_${year}_${month}.nc + fi + + if $lmerge; then + + # Loop over all 51 ensemble members (indices 0-50) + for ens in $(seq 0 $((nens - 1))); do + # Format ensemble number as 5-digit with leading zeros (1-based: ens+1) + ens_num=$(printf "%05d" $((ens + 1))) + ens_dir=real_${ens_num} + mkdir -pv ${ens_dir} + + # Copy netCDF file for this ensemble member + cp ${pathdata}/download_era5_${year}_${month}.nc ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + + # Extract the specific ensemble member + ncks --overwrite -d number,${ens} -O ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + # Remove number and forecast_reference_time dimensions + ncwa --overwrite -a forecast_reference_time ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + ncwa --overwrite -a number ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + + # Rename dimensions/variables for CDO compatibility + ncrename -v valid_time,time ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + ncrename -d forecast_period,valid_time ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + + # Remap this ensemble member + cdo -P ${ompthd} remap,${griddesfile},${wgtcaf} ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/rmp_era5_${year}_${month}_ens${ens}.nc + + cdo -P ${ompthd} expr,'WIND=sqrt(u10^2+v10^2)' ${tmpdir}/rmp_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/${year}_${month}_temp.nc # Calculate WIND from u10 and v10 + cdo -f nc4c const,10,${tmpdir}/rmp_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/${year}_${month}_const.nc + ncpdq -U ${tmpdir}/rmp_era5_${year}_${month}_ens${ens}.nc ${tmpdir}/${year}_${month}_temp2.nc + + cdo merge ${tmpdir}/${year}_${month}_const.nc ${tmpdir}/${year}_${month}_temp2.nc \ + ${tmpdir}/${year}_${month}_temp.nc ${tmpdir}/${year}_${month}_temp4.nc + + ncks -Q -C -x -v hyai,hyam,hybi,hybm ${tmpdir}/${year}_${month}_temp4.nc ${tmpdir}/${year}_${month}_temp5.nc + + # Copy to ensemble-specific directory + cp ${tmpdir}/${year}_${month}_temp5.nc ${ens_dir}/${year}-${month}.nc + + # Rename variables + ncrename -v sp,PSRF -v fsds,FSDS -v flds,FLDS -v avg_tprate,PRECTmms -v const,ZBOT -v t10m,TBOT -v q10m,QBOT ${ens_dir}/${year}-${month}.nc + +# ncap2 -O -s 'where(FSDS<0.) FSDS=0' ${ens_dir}/${year}-${month}.nc + ncatted -O -a units,ZBOT,m,c,"m" ${ens_dir}/${year}-${month}.nc + + ncks -Q -O -h --glb author="${author}" ${ens_dir}/${year}-${month}.nc ${ens_dir}/${year}-${month}.nc + ncks -Q -O -h --glb contact="${email}" ${ens_dir}/${year}-${month}.nc ${ens_dir}/${year}-${month}.nc + + # Cleanup temporary files for this ensemble member + rm ${tmpdir}/download_era5_${year}_${month}_ens${ens}.nc + rm ${tmpdir}/rmp_era5_${year}_${month}_ens${ens}.nc + rm ${tmpdir}/${year}_${month}_temp*nc ${tmpdir}/${year}_${month}_const.nc + done + fi + +done +done + diff --git a/mkforcing/seas5_daily_to_6hourly.py b/mkforcing/seas5_daily_to_6hourly.py new file mode 100644 index 0000000..26ba778 --- /dev/null +++ b/mkforcing/seas5_daily_to_6hourly.py @@ -0,0 +1,716 @@ +#!/usr/bin/env python3 +""" +Convert SEAS5 constant, daily, and 6-hourly variables to a target time resolution. + +This script: +1. Adds constant `z` (orography) broadcast to all time steps +2. Expands 6-hourly variables (msl, u10, v10, t2m, d2m) to target frequency +3. Converts daily `tp` (total precipitation, m) to precipitation flux `avg_tprate` (kg m**-2 s**-1) +4. Converts daily `strd` (thermal radiation) to flux, distributed equally +5. Converts daily `ssrd` (solar radiation) to flux with bell-shaped diurnal cycle + (cosine distribution: zero at 6:00 and 18:00, peak at noon) + +Usage: + python seas5_daily_to_6hourly.py --const --daily \\ + --hourly <6h_file> [--output ] [--frequency ] + +If --output is not specified, the output file will be named based on frequency. +""" + +import argparse +import os +import sys + +import numpy as np +import xarray as xr + + +def forecast_period_to_hours(forecast_period): + """ + Convert forecast_period values to hours. + + xarray may store forecast_period as timedelta64[ns] (nanoseconds), + but we need hours for our calculations. + + Parameters + ---------- + forecast_period : array-like + Forecast period values (may be in nanoseconds or hours) + + Returns + ------- + np.ndarray + Forecast period values in hours + """ + values = np.asarray(forecast_period) + + # Check if values are in nanoseconds (very large numbers) + # 1 hour = 3600 seconds = 3.6e12 nanoseconds + if values.dtype.kind == "m": # timedelta type + # Convert timedelta64 to hours + return values.astype("timedelta64[h]").astype(float) + elif np.max(values) > 1e9: + # Likely nanoseconds, convert to hours + return values / (3600 * 1e9) + else: + # Assume already in hours + return values.astype(float) + + +def generate_target_forecast_periods(input_periods_hours, frequency_hours, include_hour_zero=False): + """ + Generate target forecast periods at the desired frequency. + + Parameters + ---------- + input_periods_hours : np.ndarray + Input forecast periods in hours (e.g., [6, 12, 18, 24, 30, 36, 42, 48]) + frequency_hours : int + Target frequency in hours (e.g., 1 for hourly) + include_hour_zero : bool + If True, include hour 0 as the first time step (for initial conditions) + + Returns + ------- + np.ndarray + Target forecast periods in hours + """ + start_hour = 0 if include_hour_zero else frequency_hours + end_hour = int(np.max(input_periods_hours)) + return np.arange(start_hour, end_hour + frequency_hours, frequency_hours, dtype=float) + + +def expand_6hourly_to_target(ds_6h, var_name, target_periods_hours, frequency_hours): + """ + Expand a 6-hourly variable to target frequency by repeating values. + + Each 6-hourly value is assigned to all sub-intervals within that 6-hour window. + For example, with hourly output, the value at hour 6 is assigned to hours 1-6. + + Parameters + ---------- + ds_6h : xarray.Dataset + Dataset containing 6-hourly variables + var_name : str + Variable name to expand + target_periods_hours : np.ndarray + Target forecast periods in hours + frequency_hours : int + Target frequency in hours + + Returns + ------- + np.ndarray + Variable data at target frequency + """ + var_data = ds_6h[var_name] + input_periods_hours = forecast_period_to_hours(ds_6h["forecast_period"].values) + + # Get dimensions + n_number = var_data.sizes.get("number", var_data.shape[0]) + n_ref_time = var_data.sizes.get("forecast_reference_time", 1) + n_lat = var_data.sizes.get("latitude", 1) + n_lon = var_data.sizes.get("longitude", 1) + n_target = len(target_periods_hours) + + # Create output array + output = np.zeros((n_number, n_ref_time, n_target, n_lat, n_lon), dtype=np.float32) + + # For each target period, find the corresponding 6-hourly period + for i, target_hour in enumerate(target_periods_hours): + # Find the 6-hourly period that contains this target hour + # Hour 0-6 -> 6h period at hour 6, Hour 7-12 -> 6h period at hour 12, etc. + containing_6h = int(np.ceil(target_hour / 6.0) * 6) + if containing_6h == 0: + containing_6h = 6 # Hour 0 uses the first available value (hour 6) + idx_6h = np.where(np.isclose(input_periods_hours, containing_6h))[0] + if len(idx_6h) > 0: + output[:, :, i, :, :] = var_data.isel(forecast_period=idx_6h[0]).values + + return output + + +def distribute_daily_to_target(ds_daily, var_name, target_periods_hours, frequency_hours): + """ + Distribute daily accumulated variable to target frequency intervals. + + Each daily value is divided equally among all intervals within that day. + + Parameters + ---------- + ds_daily : xarray.Dataset + Dataset containing daily variables + var_name : str + Variable name (e.g., 'tp') + target_periods_hours : np.ndarray + Target forecast periods in hours + frequency_hours : int + Target frequency in hours + + Returns + ------- + np.ndarray + Variable data at target frequency + """ + var_daily = ds_daily[var_name] + daily_periods_hours = forecast_period_to_hours(ds_daily["forecast_period"].values) + + # Get dimensions + n_number = var_daily.sizes.get("number", var_daily.shape[0]) + n_ref_time = var_daily.sizes.get("forecast_reference_time", 1) + n_lat = var_daily.sizes.get("latitude", 1) + n_lon = var_daily.sizes.get("longitude", 1) + n_target = len(target_periods_hours) + + # Number of intervals per day + intervals_per_day = int(24 / frequency_hours) + + # Create output array + output = np.zeros((n_number, n_ref_time, n_target, n_lat, n_lon), dtype=np.float32) + + # For each daily period, distribute to target intervals + for i, daily_period_hours in enumerate(daily_periods_hours): + # Find target intervals belonging to this day + # For first day, include hour 0 if present + day_start = daily_period_hours - 24 + if i > 0: + day_start += frequency_hours # Avoid overlap with previous day + day_end = daily_period_hours + + indices = [] + for j, target_hour in enumerate(target_periods_hours): + if day_start <= target_hour <= day_end: + indices.append(j) + + if len(indices) > 0: + # Divide daily value by number of intervals + daily_value = var_daily.isel(forecast_period=i).values + for idx in indices: + output[:, :, idx, :, :] = daily_value / len(indices) + + return output + + +def distribute_precipitation_to_flux(ds_daily, var_name, target_periods_hours, frequency_hours): + """ + Distribute daily precipitation to target frequency and convert to flux. + + The input data is accumulated since forecast start, so we first de-accumulate + by taking differences (day[i] - day[i-1]). The first day uses its value directly. + The de-accumulated daily values are then distributed evenly across all intervals. + Converts from meters [m] to precipitation rate [kg m**-2 s**-1] (= mm/s). + + Parameters + ---------- + ds_daily : xarray.Dataset + Dataset containing daily variables (accumulated since forecast start) + var_name : str + Variable name (typically 'tp') + target_periods_hours : np.ndarray + Target forecast periods in hours + frequency_hours : int + Target frequency in hours + + Returns + ------- + np.ndarray + Precipitation flux (kg m**-2 s**-1) at target frequency + """ + precip_daily = ds_daily[var_name] + daily_periods_hours = forecast_period_to_hours(ds_daily["forecast_period"].values) + + # Get dimensions + n_number = precip_daily.sizes.get("number", precip_daily.shape[0]) + n_ref_time = precip_daily.sizes.get("forecast_reference_time", 1) + n_lat = precip_daily.sizes.get("latitude", 1) + n_lon = precip_daily.sizes.get("longitude", 1) + n_target = len(target_periods_hours) + + # Time interval in seconds + dt_seconds = frequency_hours * 3600 + + # Create output array + flux = np.zeros((n_number, n_ref_time, n_target, n_lat, n_lon), dtype=np.float32) + + # For each daily period, distribute to target intervals + for i, daily_period_hours in enumerate(daily_periods_hours): + # For first day, include hour 0 if present + day_start = daily_period_hours - 24 + if i > 0: + day_start += frequency_hours # Avoid overlap with previous day + day_end = daily_period_hours + + indices = [] + for j, target_hour in enumerate(target_periods_hours): + if day_start <= target_hour <= day_end: + indices.append(j) + + if len(indices) > 0: + # Daily value in meters [m] + # De-accumulate: get actual daily value from accumulated values + # Input is accumulated since forecast start + accumulated_value = precip_daily.isel(forecast_period=i).values + if i == 0: + # First day: use accumulated value directly + daily_value = accumulated_value + else: + # Subsequent days: take difference with previous day + previous_accumulated = precip_daily.isel(forecast_period=i - 1).values + daily_value = accumulated_value - previous_accumulated + + # Convert: m -> mm (×1000), then divide by interval seconds to get rate + # flux [kg m**-2 s**-1] = (daily_value [m] × 1000) / (n_intervals × dt_seconds) + flux_value = (daily_value * 1000) / (len(indices) * dt_seconds) + for idx in indices: + flux[:, :, idx, :, :] = flux_value + + return flux + + +def distribute_thermal_radiation_to_flux(ds_daily, var_name, target_periods_hours, frequency_hours): + """ + Distribute daily thermal radiation to target frequency and convert to flux. + + The input data is accumulated since forecast start, so we first de-accumulate + by taking differences (day[i] - day[i-1]). The first day uses its value directly. + Thermal radiation is distributed evenly across all intervals (day and night). + + Parameters + ---------- + ds_daily : xarray.Dataset + Dataset containing daily variables (accumulated since forecast start) + var_name : str + Variable name (typically 'strd') + target_periods_hours : np.ndarray + Target forecast periods in hours + frequency_hours : int + Target frequency in hours + + Returns + ------- + np.ndarray + Radiation flux (W/m²) at target frequency + """ + rad_daily = ds_daily[var_name] + daily_periods_hours = forecast_period_to_hours(ds_daily["forecast_period"].values) + + # Get dimensions + n_number = rad_daily.sizes.get("number", rad_daily.shape[0]) + n_ref_time = rad_daily.sizes.get("forecast_reference_time", 1) + n_lat = rad_daily.sizes.get("latitude", 1) + n_lon = rad_daily.sizes.get("longitude", 1) + n_target = len(target_periods_hours) + + # Time interval in seconds + dt_seconds = frequency_hours * 3600 + + # Create output array + flux = np.zeros((n_number, n_ref_time, n_target, n_lat, n_lon), dtype=np.float32) + + # For each daily period, distribute to target intervals + for i, daily_period_hours in enumerate(daily_periods_hours): + # For first day, include hour 0 if present + day_start = daily_period_hours - 24 + if i > 0: + day_start += frequency_hours # Avoid overlap with previous day + day_end = daily_period_hours + + indices = [] + for j, target_hour in enumerate(target_periods_hours): + if day_start <= target_hour <= day_end: + indices.append(j) + + if len(indices) > 0: + # De-accumulate: get actual daily value from accumulated values + # Input is accumulated since forecast start (J/m²) + accumulated_value = rad_daily.isel(forecast_period=i).values + if i == 0: + # First day: use accumulated value directly + daily_value = accumulated_value + else: + # Subsequent days: take difference with previous day + previous_accumulated = rad_daily.isel(forecast_period=i - 1).values + daily_value = accumulated_value - previous_accumulated + + # Flux = (daily_value / n_intervals) / dt_seconds + flux_value = (daily_value / len(indices)) / dt_seconds + for idx in indices: + flux[:, :, idx, :, :] = flux_value + + return flux + + +def distribute_solar_radiation_to_flux(ds_daily, var_name, target_periods_hours, frequency_hours): + """ + Distribute daily solar radiation with bell-shaped diurnal cycle. + + The input data is accumulated since forecast start, so we first de-accumulate + by taking differences (day[i] - day[i-1]). The first day uses its value directly. + + Uses a cosine distribution: + - Zero radiation between 18:00 and 06:00 (night) + - Bell-shaped curve between 06:00 and 18:00, peak at noon + - weight(h) = max(0, cos((h - 12) * π / 12)) + + Parameters + ---------- + ds_daily : xarray.Dataset + Dataset containing daily variables (accumulated since forecast start) + var_name : str + Variable name (typically 'ssrd') + target_periods_hours : np.ndarray + Target forecast periods in hours + frequency_hours : int + Target frequency in hours + + Returns + ------- + np.ndarray + Radiation flux (W/m²) at target frequency + """ + rad_daily = ds_daily[var_name] + daily_periods_hours = forecast_period_to_hours(ds_daily["forecast_period"].values) + + # Get dimensions + n_number = rad_daily.sizes.get("number", rad_daily.shape[0]) + n_ref_time = rad_daily.sizes.get("forecast_reference_time", 1) + n_lat = rad_daily.sizes.get("latitude", 1) + n_lon = rad_daily.sizes.get("longitude", 1) + n_target = len(target_periods_hours) + + # Time interval in seconds + dt_seconds = frequency_hours * 3600 + + # Create output array + flux = np.zeros((n_number, n_ref_time, n_target, n_lat, n_lon), dtype=np.float32) + + # For each daily period, distribute with bell-shaped weights + for i, daily_period_hours in enumerate(daily_periods_hours): + # For first day, include hour 0 if present + day_start = daily_period_hours - 24 + if i > 0: + day_start += frequency_hours # Avoid overlap with previous day + day_end = daily_period_hours + + # Find indices and calculate weights for this day + indices = [] + weights = [] + for j, target_hour in enumerate(target_periods_hours): + if day_start <= target_hour <= day_end: + indices.append(j) + # Hour of day (0-24) - use center of interval + hour_of_day = ((target_hour - frequency_hours / 2) % 24) + # Bell-shaped weight: cosine centered at noon + # cos((h - 12) * π / 12) gives 0 at h=6 and h=18, 1 at h=12 + weight = max(0.0, np.cos((hour_of_day - 12) * np.pi / 12)) + weights.append(weight) + + if len(indices) > 0 and sum(weights) > 0: + # Normalize weights to sum to 1 + weights = np.array(weights) + weights = weights / weights.sum() + + # De-accumulate: get actual daily value from accumulated values + # Input is accumulated since forecast start (J/m²) + accumulated_value = rad_daily.isel(forecast_period=i).values + if i == 0: + # First day: use accumulated value directly + daily_value = accumulated_value + else: + # Subsequent days: take difference with previous day + previous_accumulated = rad_daily.isel(forecast_period=i - 1).values + daily_value = accumulated_value - previous_accumulated + + # Distribute according to weights and convert to flux + for idx, weight in zip(indices, weights): + # Energy for this interval = daily_value * weight + # Flux = energy / dt_seconds + flux[:, :, idx, :, :] = (daily_value * weight) / dt_seconds + + return flux + + +def main(): + parser = argparse.ArgumentParser( + description="Convert SEAS5 variables to target time resolution" + ) + parser.add_argument( + "--const", required=True, help="Path to netCDF file with constant variables (z)" + ) + parser.add_argument( + "--daily", + required=True, + help="Path to netCDF file with daily variables (strd, ssrd, tp)", + ) + parser.add_argument( + "--hourly", required=True, help="Path to netCDF file with 6-hourly variables" + ) + parser.add_argument( + "--output", + default=None, + help="Output file path (default: auto-generated based on frequency)", + ) + parser.add_argument( + "--frequency", + type=int, + default=6, + choices=[1, 2, 3, 6], + help="Target frequency in hours (default: 6). Must divide 24 evenly.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print what would be done without writing output", + ) + parser.add_argument( + "--include-hour-zero", + action="store_true", + help="Include hour 0 as initial time step (duplicates first available value)", + ) + + args = parser.parse_args() + + frequency_hours = args.frequency + + # Determine output file + if args.output: + output_file = args.output + else: + base = os.path.splitext(args.hourly)[0] + output_file = f"{base}_{frequency_hours}h.nc" + + # Check if output directory exists + output_dir = os.path.dirname(output_file) + if output_dir and not os.path.exists(output_dir): + print(f"Error: Output directory does not exist: {output_dir}") + print(f"Please create it with: mkdir -p {output_dir}") + sys.exit(1) + + print(f"Target frequency: {frequency_hours} hours") + print(f"Output file: {output_file}") + + # Load datasets + print(f"Loading constant file: {args.const}") + ds_const = xr.open_dataset(args.const).load() + + print(f"Loading daily file: {args.daily}") + ds_daily = xr.open_dataset(args.daily).load() + + print(f"Loading 6-hourly file: {args.hourly}") + ds_6h = xr.open_dataset(args.hourly).load() + + # Get input periods and generate target periods + input_periods_hours = forecast_period_to_hours(ds_6h["forecast_period"].values) + target_periods_hours = generate_target_forecast_periods( + input_periods_hours, frequency_hours, args.include_hour_zero + ) + n_target = len(target_periods_hours) + + daily_periods_hours = forecast_period_to_hours(ds_daily["forecast_period"].values) + + print(f"Input 6-hourly periods (hours): {input_periods_hours}") + print(f"Daily periods (hours): {daily_periods_hours}") + print(f"Target periods (hours): {target_periods_hours}") + + # 1. Expand constant z + print("Expanding constant z...") + z_const = ds_const["z"] + z_squeezed = z_const.squeeze("forecast_period", drop=True) + z_data = np.tile( + z_squeezed.values[:, :, np.newaxis, :, :], (1, 1, n_target, 1, 1) + ) + + # 2. Expand 6-hourly variables + vars_6h = ["msl", "u10", "v10", "t2m", "d2m"] + expanded_vars = {} + for var in vars_6h: + if var in ds_6h.data_vars: + print(f"Expanding {var} to {frequency_hours}-hourly...") + expanded_vars[var] = expand_6hourly_to_target( + ds_6h, var, target_periods_hours, frequency_hours + ) + + # 3. Distribute daily precipitation and convert to flux + print("Distributing daily precipitation and converting to flux...") + tprate_target = distribute_precipitation_to_flux( + ds_daily, "tp", target_periods_hours, frequency_hours + ) + + # 4. Distribute thermal radiation (evenly) + print("Distributing thermal radiation (flds)...") + flds_target = distribute_thermal_radiation_to_flux( + ds_daily, "strd", target_periods_hours, frequency_hours + ) + + # 5. Distribute solar radiation (bell-shaped) + print("Distributing solar radiation with bell-shaped diurnal cycle (fsds)...") + fsds_target = distribute_solar_radiation_to_flux( + ds_daily, "ssrd", target_periods_hours, frequency_hours + ) + + # Create output dataset + print("Creating output dataset...") + + dims = ["number", "forecast_reference_time", "forecast_period", "latitude", "longitude"] + + # Start with coordinates + ds_out = xr.Dataset( + coords={ + "number": ds_6h["number"], + "forecast_reference_time": ds_6h["forecast_reference_time"], + "forecast_period": target_periods_hours, + "latitude": ds_6h["latitude"], + "longitude": ds_6h["longitude"], + } + ) + + # Copy coordinate attributes + ds_out["forecast_period"].attrs = { + "long_name": "time since forecast_reference_time", + "standard_name": "forecast_period", + "units": "hours", + } + + # Compute valid_time = forecast_reference_time + forecast_period + # forecast_reference_time may be datetime64 or seconds since epoch + ref_time = ds_6h["forecast_reference_time"].values[0] + if np.issubdtype(type(ref_time), np.datetime64): + # Convert datetime64 to seconds since epoch + ref_time_seconds = (ref_time - np.datetime64('1970-01-01T00:00:00')) / np.timedelta64(1, 's') + else: + ref_time_seconds = float(ref_time) + valid_time_seconds = (ref_time_seconds + target_periods_hours * 3600).astype(np.int64) + # ds_out["valid_time"] = xr.DataArray( + # data=valid_time_seconds, + # dims=["forecast_period"], + # attrs={ + # "standard_name": "time", + # "long_name": "time", + # "units": "seconds since 1970-01-01", + # "calendar": "proleptic_gregorian", + # }, + # ) + + # Also add time in hours since 1900-01-01 (common format for climate data) + # Hours from 1900-01-01 to 1970-01-01: 613608 hours + hours_1900_to_1970 = 613608 + time_hours = (valid_time_seconds / 3600 + hours_1900_to_1970).astype(np.int32) + ds_out["valid_time"] = xr.DataArray( + data=time_hours, + dims=["forecast_period"], + attrs={ + "standard_name": "time", + "long_name": "time", + "units": "hours since 1900-01-01 00:00:00.0", + "calendar": "gregorian", + "axis": "T", + }, + ) + + # Add z + ds_out["z"] = xr.DataArray(data=z_data, dims=dims, attrs=z_const.attrs) + + # Add expanded 6-hourly variables + for var, data in expanded_vars.items(): + ds_out[var] = xr.DataArray(data=data, dims=dims, attrs=ds_6h[var].attrs) + + # Add precipitation flux + ds_out["avg_tprate"] = xr.DataArray( + data=tprate_target, + dims=dims, + attrs={ + "units": "kg m**-2 s**-1", + "long_name": "Time-mean total precipitation rate", + "standard_name": "precipitation_flux", + "description": f"Converted from daily tp [m], distributed equally to {frequency_hours}-hourly intervals", + }, + ) + + # Add radiation fluxes + ds_out["flds"] = xr.DataArray( + data=flds_target, + dims=dims, + attrs={ + "units": "W m-2", + "long_name": "Downward longwave radiation at surface", + "standard_name": "surface_downwelling_longwave_flux_in_air", + "description": f"Converted from daily strd, distributed equally to {frequency_hours}-hourly intervals", + }, + ) + + ds_out["fsds"] = xr.DataArray( + data=fsds_target, + dims=dims, + attrs={ + "units": "W m-2", + "long_name": "Downward shortwave radiation at surface", + "standard_name": "surface_downwelling_shortwave_flux_in_air", + "description": f"Converted from daily ssrd with bell-shaped diurnal cycle (cosine, peak at noon)", + }, + ) + + # Copy global attributes + ds_out.attrs = ds_6h.attrs.copy() + ds_out.attrs["frequency"] = f"{frequency_hours} hours" + + # Write output + if args.dry_run: + print(f"\n[DRY RUN] Would write output to: {output_file}") + print("[DRY RUN] Output dataset structure:") + print(ds_out) + else: + print(f"Writing output to: {output_file}") + + # Build encoding + valid_encodings = { + "compression", "dtype", "least_significant_digit", "zlib", + "_FillValue", "fletcher32", "complevel", "chunksizes", + "shuffle", "contiguous", "calendar", "units", + } + + def filter_encoding(enc): + return {k: v for k, v in enc.items() if k in valid_encodings} + + encoding = {} + + # Default encoding for all data variables + for var in ds_out.data_vars: + encoding[var] = {"zlib": True, "complevel": 4, "dtype": "float32"} + + # Preserve encoding from source files where applicable + for var in vars_6h: + if var in ds_6h.data_vars and ds_6h[var].encoding: + encoding[var].update(filter_encoding(ds_6h[var].encoding)) + + if ds_const["z"].encoding: + encoding["z"].update(filter_encoding(ds_const["z"].encoding)) + + if ds_daily["tp"].encoding: + encoding["avg_tprate"].update(filter_encoding(ds_daily["tp"].encoding)) + + # valid_time should be int64, time should be int32 + # encoding["valid_time"] = {"dtype": "int64"} + encoding["valid_time"] = {"dtype": "int32"} + + ds_out.to_netcdf(output_file, encoding=encoding) + + # Close datasets + ds_const.close() + ds_daily.close() + ds_6h.close() + + print("Done!") + + # Print summary + print(f"\nSummary (target frequency: {frequency_hours} hours):") + print(f" z: constant orography, broadcast to {n_target} time steps") + for var in vars_6h: + if var in expanded_vars: + print(f" {var}: expanded from 6-hourly by value repetition") + print(f" avg_tprate: precipitation rate [kg m**-2 s**-1] - converted from daily tp [m]") + print(f" flds: thermal radiation flux - distributed equally (day and night)") + print(f" fsds: solar radiation flux - bell-shaped (cosine, 6-18h, peak at noon)") + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..b7cfbfe --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,15 @@ +[project] +name = "eclm-atm-forcing-generator" +version = "0.1.0" +description = "Add your description here" +readme = "README.md" +requires-python = ">=3.11" +dependencies = [ + "cdsapi>=0.7.7", + "netcdf4>=1.7.4", + "numpy>=1.25.1", + "xarray>=2023.8.0", +] + +[tool.setuptools] +packages = [] diff --git a/run_atm_forcing_generator.sh b/run_atm_forcing_generator.sh new file mode 100755 index 0000000..2b7dfc5 --- /dev/null +++ b/run_atm_forcing_generator.sh @@ -0,0 +1,84 @@ +#!/usr/bin/env bash +set -euo pipefail +if [[ -z "$1" || -z "$2" || -z "$3" ]]; then + echo "Usage: $0 MODE YEAR MONTH" + exit 1 +fi + +MODE=$1 +YEAR=$2 +MONTH=$3 +DOMAINFILE="domain.lnd.DE-RuS_DE-RuS.250926.nc" +SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" + +source "${SCRIPT_DIR}/jsc.2024_Intel.sh" + +VENV_DIR="${SCRIPT_DIR}/pyvenv_eclm_atm_forcing_generator" +if [[ ! -d "$VENV_DIR" ]]; then + python -m venv "$VENV_DIR" +fi +source "${VENV_DIR}/bin/activate" +pip install "${SCRIPT_DIR}" + +mkdir -p "${SCRIPT_DIR}/${YEAR}-${MONTH}" +if [[ "$MODE" == "ERA5" ]]; then + echo "$MODE not available" + # mkdir -p data + # uv run mkforcing/download_ERA5_input.py \ + # --year $YEAR \ + # --month $MONTH \ + # --dirout data \ + # --request "${HOME}/mkforcing/custom_request_ERA5.py" \ + # --domainfile "${HOME}/domain.nc" + # unzip "data/download_era5_${YEAR}_${MONTH}.zip" -d data/ + # uv run mkforcing/dewpoint_to_specific_humidity.py data/data_stream-oper_stepType-instant.nc + # uv run mkforcing/2m_to_10m_conversion.py data/data_stream-oper_stepType-instant.nc + # mkforcing/prepare_ERA5_input.sh \ + # lrenametime=true lmeteo=false \ + # lunzip=false wgtcaf=../wgtdis_era5caf_to_domain.nc \ + # griddesfile=../domain_griddef.txt iyear=$YEAR \ + # imonth=$MONTH \ + # pathdata=../data \ + # lwgtdis=true lgriddes=true domainfile="${HOME}/domain.nc" + +else + python "${SCRIPT_DIR}/mkforcing/download_ERA5_input.py" \ + --year $YEAR --month ${MONTH} \ + --dirout "${SCRIPT_DIR}/cdsapidwn_SEAS5_const" \ + --request "${SCRIPT_DIR}/mkforcing/custom_request_SEAS5_const.py" # \ + # --domainfile $DOMAINFILE + python "${SCRIPT_DIR}/mkforcing/download_ERA5_input.py" \ + --year ${YEAR} --month ${MONTH} \ + --dirout "${SCRIPT_DIR}/cdsapidwn_SEAS5_24h" \ + --request "${SCRIPT_DIR}/mkforcing/custom_request_SEAS5_24h.py" # \ + # --domainfile $DOMAINFILE + python "${SCRIPT_DIR}/mkforcing/download_ERA5_input.py" \ + --year ${YEAR} --month ${MONTH} \ + --dirout "${SCRIPT_DIR}/cdsapidwn_SEAS5_06h" \ + --request "${SCRIPT_DIR}/mkforcing/custom_request_SEAS5_06h.py" # \ + # --domainfile $DOMAINFILE + mkdir -p "${SCRIPT_DIR}/cdsapidwn_SEAS5" + python "${SCRIPT_DIR}/mkforcing/seas5_daily_to_6hourly.py" \ + --const "${SCRIPT_DIR}/cdsapidwn_SEAS5_const/download_era5_${YEAR}_${MONTH}.nc" \ + --daily "${SCRIPT_DIR}/cdsapidwn_SEAS5_24h/download_era5_${YEAR}_${MONTH}.nc" \ + --hourly "${SCRIPT_DIR}/cdsapidwn_SEAS5_06h/download_era5_${YEAR}_${MONTH}.nc" \ + --output "${SCRIPT_DIR}/cdsapidwn_SEAS5/download_era5_${YEAR}_${MONTH}.nc" \ + --frequency 3 --include-hour-zero + python "${SCRIPT_DIR}/mkforcing/orography_to_elevation.py" \ + "${SCRIPT_DIR}/cdsapidwn_SEAS5/download_era5_${YEAR}_${MONTH}.nc" + python "${SCRIPT_DIR}/mkforcing/mslp_to_sp.py" \ + "${SCRIPT_DIR}/cdsapidwn_SEAS5/download_era5_${YEAR}_${MONTH}.nc" \ + --elevation-var elevation + python "${SCRIPT_DIR}/mkforcing/dewpoint_to_specific_humidity.py" \ + "${SCRIPT_DIR}/cdsapidwn_SEAS5/download_era5_${YEAR}_${MONTH}.nc" + python "${SCRIPT_DIR}/mkforcing/2m_to_10m_conversion.py" \ + "${SCRIPT_DIR}/cdsapidwn_SEAS5/download_era5_${YEAR}_${MONTH}.nc" + + "${SCRIPT_DIR}/mkforcing/prepare_SEAS5_input.sh" \ + lwgtdis=true lgriddes=true \ + domainfile="${SCRIPT_DIR}/${DOMAINFILE}" \ + griddesfile="${SCRIPT_DIR}/domain_griddef.txt" \ + wgtcaf="${SCRIPT_DIR}/wgtdis_era5caf_to_domain.nc" \ + iyear=${YEAR} imonth=${MONTH} \ + pathdata="${SCRIPT_DIR}/cdsapidwn_SEAS5" +fi