diff --git a/esmvalcore/cmor/tables/custom/CMOR_alt.dat b/esmvalcore/cmor/tables/custom/CMOR_alt.dat new file mode 100644 index 0000000000..4d7475f9b8 --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_alt.dat @@ -0,0 +1,21 @@ +SOURCE: CMIP6 +!============ +variable_entry: alt +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: m +cell_methods: time: mean +cell_measures: area: areacella +long_name: Active Layer Thickness +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat index 490ddab6f5..8c42a9454e 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat @@ -65,6 +65,28 @@ must_have_bounds: no !---------------------------------- ! +!============ +axis_entry: sdepth +!============ +!---------------------------------- +! Axis attributes: +!---------------------------------- +standard_name: depth +units: m +axis: Z ! X, Y, Z, T (default: undeclared) +positive: down ! up or down (default: undeclared) +long_name: depth +!---------------------------------- +! Additional axis information: +!---------------------------------- +out_name: depth +valid_min: 0.0 +valid_max: 200.0 +stored_direction: increasing +type: double +must_have_bounds: yes +!---------------------------------- +! !============ axis_entry: time diff --git a/esmvalcore/cmor/tables/custom/CMOR_gtd.dat b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat new file mode 100644 index 0000000000..16248628da --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat @@ -0,0 +1,22 @@ +SOURCE: CMIP6 +!============ +variable_entry: gtd +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: K +cell_methods: time: mean +cell_measures: area: areacella +long_name: Permafrost Ground Temperature +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude sdepth time +out_name: gtd +type: real +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/custom/CMOR_pfr.dat b/esmvalcore/cmor/tables/custom/CMOR_pfr.dat new file mode 100644 index 0000000000..2900b0595a --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_pfr.dat @@ -0,0 +1,25 @@ +SOURCE: CMIP6 +!============ +variable_entry: pfr +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: % +cell_methods: time: mean +cell_measures: area: areacella +long_name: Permafrost Extent +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +valid_min: 0.0 +valid_max: 100.0 +ok_min_mean_abs: 0.0 +ok_max_mean_abs: 100.0 +!---------------------------------- +! diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py new file mode 100644 index 0000000000..182a598141 --- /dev/null +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -0,0 +1,143 @@ +"""Derivation of variable `pfr`.""" + +import logging + +import dask.array as da +import iris +import numpy as np +from iris import NameConstraint +from iris.time import PartialDateTime + +from ._baseclass import DerivedVariableBase + +logger = logging.getLogger(__name__) + +# Constants +THRESH_TEMPERATURE = 273.15 +FROZEN_MONTHS = 24 # valid range: 12-36 + + +class DerivedVariable(DerivedVariableBase): + """Derivation of variable `pfr` (permafrost extent).""" + + @staticmethod + def required(project): + """Declare the variables needed for derivation.""" + return [ + {"short_name": "tsl", "mip": "Lmon"}, + {"short_name": "sftlf", "mip": "fx"}, + {"short_name": "mrsos", "mip": "Lmon"}, + ] + + @staticmethod + def calculate(cubes): + """Compute permafrost extent. + + Permafrost is assumed if + - soil temperature in the deepest level is < 0°C + - for at least 24 consecutive months + - ice covered part of grid cell is excluded + Reference: Burke, E. J., Y. Zhang, and G. Krinner: + Evaluating permafrost physics in the Coupled Model + Intercomparison Project 6 (CMIP6) models and their + sensitivity to climate change, The Cryosphere, 14, + 3155-3174, doi: 10.5194/tc-14-3155-2020, 2020. + """ + # create a mask of land fraction (%) over ice-free grid cells + # use soil moisture as proxy for ice / ice-free grid cells + # 1) annual mean of fraction of grid cell covered with ice (%) + # assumption: top soil moisture = 0 --> ice covered + mrsos = cubes.extract_cube(NameConstraint(var_name="mrsos")) + iris.coord_categorisation.add_year(mrsos, "time") + mrsos_yr = mrsos.aggregated_by(["year"], iris.analysis.MEAN) + mrsos_yr.data = da.where(mrsos_yr.data < 0.001, 0.0, 1.0) + # 2) fraction of land cover of grid cell (%) (constant) + landfrac = cubes.extract_cube(NameConstraint(var_name="sftlf")) + # 3) create mask with fraction of ice-free land (%) + + # latitude/longitude coordinates of mrsos and sftlf sometimes + # differ by a very small amount for some models (probably because + # of rounding errors) preventing iris to do the math + # --> overwrite latitudes/longitudes in sftlf + + # fix longitudes if maximum differences are smaller than 1.0e-4 + x_coord1 = mrsos.coord(axis="X") + x_coord2 = landfrac.coord(axis="X") + delta_x_max = np.amax(x_coord1.core_points() - x_coord2.core_points()) + if delta_x_max != 0.0: + if abs(delta_x_max) < 1.0e-4: + x_coord2.points = x_coord1.points + x_coord2.bounds = x_coord1.bounds + else: + logger.error( + "Longitudes of mrsos and stflf fields differ (max = %f).", + delta_x_max, + ) + + # fix latitudes if maximum differences are smaller than 1.0e-4 + y_coord1 = mrsos.coord(axis="Y") + y_coord2 = landfrac.coord(axis="Y") + delta_y_max = np.amax(y_coord1.core_points() - y_coord2.core_points()) + if delta_y_max != 0.0: + if abs(delta_y_max) < 1.0e-4: + y_coord2.points = y_coord1.points + y_coord2.bounds = y_coord1.bounds + else: + logger.error( + "Latitudes of mrsos and stflf fields differ (max = %f).", + delta_y_max, + ) + + mask = iris.analysis.maths.multiply(mrsos_yr, landfrac) + + # extract deepest soil level + soiltemp = cubes.extract_cube(NameConstraint(var_name="tsl")) + z_coord = soiltemp.coord(axis="Z") + zmax = np.amax(z_coord.core_points()) + soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) + soiltemp.data = da.where( + soiltemp.core_data() < THRESH_TEMPERATURE, + 1, + 0, + ) + iris.coord_categorisation.add_year(soiltemp, "time") + + # prepare cube for permafrost extent with yearly time steps + pfr_yr = soiltemp.aggregated_by(["year"], iris.analysis.MEAN) + # get years to process + year_coord = pfr_yr.coord("year") + # calculate time period before and after current year to include + # in test for permafrost + test_period = (FROZEN_MONTHS - 12) / 2 + # loop over all years and test if frost is present throughout + # the whole test period, i.e. [year-test_period, year+test_period] + + for tidx, year in enumerate(year_coord.points): + # extract test period + pdt1 = PartialDateTime( + year=year - 1, + month=13 - test_period, + day=1, + ) + pdt2 = PartialDateTime(year=year + 1, month=test_period + 1, day=1) + daterange = iris.Constraint( + time=lambda cell, pdt1=pdt1, pdt2=pdt2: pdt1 + <= cell.point + < pdt2, + ) + soiltemp_window = soiltemp.extract(daterange) + # remove auxiliary coordinate 'year' to avoid lots of warnings + # from iris + soiltemp_window.remove_coord("year") + # calculate mean over test period + test_cube = soiltemp_window.collapsed("time", iris.analysis.MEAN) + # if all months in test period show soil tempeatures below zero + # then mark grid cell with "1" as permafrost and "0" otherwise + pfr_yr.data[tidx, :, :] = da.where(test_cube.data > 0.99, 1, 0) + + pfr_yr = pfr_yr * mask + pfr_yr.units = "%" + pfr_yr.rename("Permafrost extent") + pfr_yr.var_name = "pfr" + + return pfr_yr diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py new file mode 100644 index 0000000000..c7a9ddb13a --- /dev/null +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -0,0 +1,72 @@ +"""Test derivation of ``pfr``.""" + +import dask.array as da +import iris +import numpy as np +import pytest + +from esmvalcore.preprocessor._derive import pfr + +from .test_shared import get_cube + + +@pytest.fixture +def cubes(): + time_coord = iris.coords.DimCoord([0.0, 1.0, 2.0], standard_name="time") +# tsl_cube = iris.cube.Cube( +# [[[20, 10, 0], [10, 10, 10]], [[10, 10, 0], [10, 10, 20]], [[10, 10, 20], [10, 10, 10]]], +# units="K", +# standard_name="soil_temperature", +# var_name="tsl", +# dim_coords_and_dims=[(time_coord, 0, 1)], +# ) + tsl_cube = get_cube( + [[[[270.0]], [[260.0]], [[250.0]]]], + air_pressure_coord=False, + depth_coord=True, + standard_name="soil_temperature", + var_name="tsl", + units="K", + ) + sftlf_cube = get_cube( + [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], + air_pressure_coord=False, + depth_coord=False, + units="%", + standard_name="land_area_fraction", + var_name="sftlf", + dim_coords_and_dims=[(time_coord, 0)], + ) + mrsos_cube = iris.cube.Cube( + [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], + air_pressure_coord=False, + depth_coord=False, + units="kg m-2", + standard_name="mass_content_of_water_in_soil_layer", + var_name="mrsos", + dim_coords_and_dims=[(time_coord, 0)], + ) + return iris.cube.CubeList([tsl_cube, sftlf_cube, mrsos_cube]) + + +def test_pfr_calculation(cubes): + """Test function ``calculate``.""" + derived_var = pfr.DerivedVariable() + out_cube = derived_var.calculate(cubes) + assert out_cube.units == cf_units.Unit("%") + out_data = out_cube.data + expected = np.ma.ones_like(cubes[0].data) + expected[0][0][0] = 1.0 + np.testing.assert_array_equal(out_data.mask, expected.mask) + np.testing.assert_array_equal(out_data[0][0][0], expected[0][0][0]) + + +def test_pfr_required(): + """Test function ``required``.""" + derived_var = pfr.DerivedVariable() + output = derived_var.required(None) + assert output == [ + {"short_name": "tsl", "mip": "Lmon"}, + {"short_name": "sftlf", "mip": "fx"}, + {"short_name": "mrsos", "mip": "Lmon"}, + ] \ No newline at end of file