diff --git a/src/Zonstat.cc b/src/Zonstat.cc
index 53582e4c0fddc465463bb2d5ea840da17de6cef0..eb5fa624e2c46925cbaa2b2a057e75ea8831a835 100644
--- a/src/Zonstat.cc
+++ b/src/Zonstat.cc
@@ -50,19 +50,20 @@ static void
 add_operators(void)
 {
   // clang-format off
-  cdo_operator_add("zonmin",   FieldFunc_Min,   0, nullptr);
-  cdo_operator_add("zonmax",   FieldFunc_Max,   0, nullptr);
-  cdo_operator_add("zonrange", FieldFunc_Range, 0, nullptr);
-  cdo_operator_add("zonsum",   FieldFunc_Sum,   0, nullptr);
-  cdo_operator_add("zonmean",  FieldFunc_Mean,  0, nullptr);
-  cdo_operator_add("zonavg",   FieldFunc_Avg,   0, nullptr);
-  cdo_operator_add("zonvar",   FieldFunc_Var,   0, nullptr);
-  cdo_operator_add("zonvar1",  FieldFunc_Var1,  0, nullptr);
-  cdo_operator_add("zonstd",   FieldFunc_Std,   0, nullptr);
-  cdo_operator_add("zonstd1",  FieldFunc_Std1,  0, nullptr);
-  cdo_operator_add("zonskew",  FieldFunc_Skew,  0, nullptr);
-  cdo_operator_add("zonkurt",  FieldFunc_Kurt,  0, nullptr);
-  cdo_operator_add("zonpctl",  FieldFunc_Pctl,  0, nullptr);
+  cdo_operator_add("zonmin",    FieldFunc_Min,    0, nullptr);
+  cdo_operator_add("zonmax",    FieldFunc_Max,    0, nullptr);
+  cdo_operator_add("zonrange",  FieldFunc_Range,  0, nullptr);
+  cdo_operator_add("zonsum",    FieldFunc_Sum,    0, nullptr);
+  cdo_operator_add("zonmean",   FieldFunc_Mean,   0, nullptr);
+  cdo_operator_add("zonavg",    FieldFunc_Avg,    0, nullptr);
+  cdo_operator_add("zonvar",    FieldFunc_Var,    0, nullptr);
+  cdo_operator_add("zonvar1",   FieldFunc_Var1,   0, nullptr);
+  cdo_operator_add("zonstd",    FieldFunc_Std,    0, nullptr);
+  cdo_operator_add("zonstd1",   FieldFunc_Std1,   0, nullptr);
+  cdo_operator_add("zonskew",   FieldFunc_Skew,   0, nullptr);
+  cdo_operator_add("zonkurt",   FieldFunc_Kurt,   0, nullptr);
+  cdo_operator_add("zonmedian", FieldFunc_Median, 0, nullptr);
+  cdo_operator_add("zonpctl",   FieldFunc_Pctl,   0, nullptr);
   // clang-format on
 }
 
diff --git a/test/Zonstat.test.in b/test/Zonstat.test.in
index 48b065625a378da4f095ed75f17affd3bcb3687d..a19f51153984ea4ffaeed4d47ebd0346760e531c 100644
--- a/test/Zonstat.test.in
+++ b/test/Zonstat.test.in
@@ -1,13 +1,13 @@
 #! @BASH@
 #
-echo 1..10 # Number of tests to be executed.
+echo 1..13 # Number of tests to be executed.
 #
 . ./cdoTestFunctions.test
 #
 test -n "$CDO"      || CDO="@abs_top_builddir@/src/cdo $CDO_DEBUG"
 test -n "$DATAPATH" || DATAPATH="@abs_top_srcdir@/test/data/"
 #
-OPERATORS="zonmin zonmax zonrange zonsum zonavg zonmean zonstd zonstd1 zonvar zonvar1"
+OPERATORS="zonmin zonmax zonrange zonsum zonavg zonmean zonstd zonstd1 zonvar zonvar1 zonskew zonkurt zonmedian"
 #
 IFILE=$DATAPATH/t21_geosp_tsurf.grb
 #
diff --git a/test/data/Makefile.am b/test/data/Makefile.am
index 1a16c50b58f4deedbee96f151623ddd9ef6eba34..be7fdaf6ec76000ab87219db687f3452f859ff05 100644
--- a/test/data/Makefile.am
+++ b/test/data/Makefile.am
@@ -37,7 +37,7 @@ FLDSTAT2     = fldcor_ref fldcovar_ref
 FLDSTAT      = fldmin_ref fldmax_ref fldsum_ref fldavg_ref fldmean_ref fldstd_ref fldstd1_ref fldvar_ref fldvar1_ref fldrange_ref fldkurt_ref fldskew_ref fldmedian_ref
 FLDPSTAT     = fldpctl1_ref fldpctl20_ref fldpctl25_ref fldpctl33_ref fldpctl50_ref fldpctl66_ref fldpctl75_ref fldpctl80_ref fldpctl99_ref fldpctl100_ref
 MERSTAT      = mermin_ref mermax_ref mersum_ref meravg_ref mermean_ref merstd_ref merstd1_ref mervar_ref mervar1_ref merrange_ref
-ZONSTAT      = zonmin_ref zonmax_ref zonsum_ref zonavg_ref zonmean_ref zonstd_ref zonstd1_ref zonvar_ref zonvar1_ref zonrange_ref
+ZONSTAT      = zonmin_ref zonmax_ref zonsum_ref zonavg_ref zonmean_ref zonstd_ref zonstd1_ref zonvar_ref zonvar1_ref zonrange_ref zonskew_ref zonkurt_ref zonmedian_ref
 ENSSTAT      = ensmin_ref ensmax_ref enssum_ref ensavg_ref ensmean_ref ensstd_ref ensstd1_ref ensvar_ref ensvar1_ref ensrange_ref ensskew_ref enskurt_ref ensmedian_ref
 ENSMSTAT     = ensmmin_ref ensmmax_ref ensmsum_ref ensmavg_ref ensmmean_ref ensmstd_ref ensmstd1_ref ensmvar_ref ensmvar1_ref ensmrange_ref ensmskew_ref ensmkurt_ref ensmmedian_ref
 ENSPCTL      = enspctl1_ref enspctl20_ref enspctl25_ref enspctl33_ref enspctl50_ref enspctl66_ref enspctl75_ref enspctl80_ref enspctl99_ref enspctl100_ref
diff --git a/test/data/gen_refdata.sh b/test/data/gen_refdata.sh
index 3c52a0017fea9efbfe8deded4788a51c0eef7d5a..ecec384b3f6b97eeb202c9558d1481c8053a1e1e 100755
--- a/test/data/gen_refdata.sh
+++ b/test/data/gen_refdata.sh
@@ -6,6 +6,27 @@ FORMAT="-f srv -b F32"
 #
 ########################################################################
 #
+# Zonstat
+#
+STATS="min max range sum avg mean std std1 var var1 skew kurt median"
+STATS="skew kurt median"
+IFILE=t21_geosp_tsurf.grb
+for STAT in $STATS; do
+  $CDO $FORMAT zon$STAT $IFILE zon${STAT}_ref
+done
+exit
+########################################################################
+#
+# Merstat
+#
+STATS="min max range sum avg mean std std1 var var1 skew kurt median"
+IFILE=t21_geosp_tsurf.grb
+for STAT in $STATS; do
+  $CDO $FORMAT mer$STAT $IFILE mer${STAT}_ref
+done
+exit
+########################################################################
+#
 # Ensstat
 #
 IFILE=ts_mm_5years
@@ -509,26 +530,6 @@ for STAT in $STATS; do
 done
 ########################################################################
 #
-# Zonstat
-#
-STATS="min max sum avg mean std std1 var var1 range"
-IFILE=t21_geosp_tsurf.grb
-for STAT in $STATS; do
-  $CDO $FORMAT zon$STAT $IFILE zon${STAT}_ref
-done
-exit
-########################################################################
-#
-# Merstat
-#
-STATS="min max sum avg mean std std1 var var1"
-IFILE=t21_geosp_tsurf.grb
-for STAT in $STATS; do
-  $CDO $FORMAT mer$STAT $IFILE mer${STAT}_ref
-done
-exit
-########################################################################
-#
 # Arith
 #
 IFILE=arith1.srv
diff --git a/test/data/zonkurt_ref b/test/data/zonkurt_ref
new file mode 100644
index 0000000000000000000000000000000000000000..7056132ae5ffd74599021aee957728ef72a12b68
Binary files /dev/null and b/test/data/zonkurt_ref differ
diff --git a/test/data/zonmedian_ref b/test/data/zonmedian_ref
new file mode 100644
index 0000000000000000000000000000000000000000..747517e75aaef6e9865ee1e73a840d43ed186be3
Binary files /dev/null and b/test/data/zonmedian_ref differ
diff --git a/test/data/zonskew_ref b/test/data/zonskew_ref
new file mode 100644
index 0000000000000000000000000000000000000000..07b5d1501c88c445b2ff33c25b3847f151667b09
Binary files /dev/null and b/test/data/zonskew_ref differ