Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/cwt.py
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/cwt.py	(revision 7881)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/cwt.py	(revision 7881)
@@ -0,0 +1,274 @@
+# Cross-walking to join 31 Koeppen-Geiger climate zones into 7
+cwt_kg31to7 = [
+    ["KG31/KG7", "tropical", "arid_warm", "arid_cool", "temp_warm", "temp_cool", "boreal_warm", "boreal_cool"],
+    ["Af",       1,          0,           0,           0,           0,           0,             0],
+    ["Am",       1,          0,           0,           0,           0,           0,             0],
+    ["As",       1,          0,           0,           0,           0,           0,             0],
+    ["Aw",       1,          0,           0,           0,           0,           0,             0],
+    ["BSh",      0,          1,           0,           0,           0,           0,             0],
+    ["BSk",      0,          0,           1,           0,           0,           0,             0],
+    ["BWh",      0,          1,           0,           0,           0,           0,             0],
+    ["BWk",      0,          0,           1,           0,           0,           0,             0],
+    ["Cfa",      0,          0,           0,           1,           0,           0,             0],
+    ["Cfb",      0,          0,           0,           0,           1,           0,             0],
+    ["Cfc",      0,          0,           0,           0,           1,           0,             0],
+    ["Csa",      0,          0,           0,           1,           0,           0,             0],
+    ["Csb",      0,          0,           0,           0,           1,           0,             0],
+    ["Csc",      0,          0,           0,           0,           1,           0,             0],
+    ["Cwa",      0,          0,           0,           1,           0,           0,             0],
+    ["Cwb",      0,          0,           0,           0,           1,           0,             0],
+    ["Cwc",      0,          0,           0,           0,           1,           0,             0],
+    ["Dfa",      0,          0,           0,           0,           1,           0,             0],
+    ["Dfb",      0,          0,           0,           0,           0,           1,             0],
+    ["Dfc",      0,          0,           0,           0,           0,           0,             1],
+    ["Dfd",      0,          0,           0,           0,           0,           0,             1],
+    ["Dsa",      0,          0,           0,           0,           0,           1,             0],
+    ["Dsb",      0,          0,           0,           0,           0,           1,             0],
+    ["Dsc",      0,          0,           0,           0,           0,           0,             1],
+    ["Dsd",      0,          0,           0,           0,           0,           0,             1],
+    ["Dwa",      0,          0,           0,           0,           0,           1,             0],
+    ["Dwb",      0,          0,           0,           0,           0,           1,             0],
+    ["Dwc",      0,          0,           0,           0,           0,           0,             1],
+    ["Dwd",      0,          0,           0,           0,           0,           0,             1],
+    ["EF",       0,          0,           0,           0,           0,           0,             1],
+    ["ET",       0,          0,           0,           0,           0,           0,             1],
+]
+
+# Cross-walking to join 31 Koeppen-Geiger climate zones into 3
+cwt_kg31to3 = [
+    ["KG31/KG3", "tropical", "temperate", "boreal"],
+    ["Af",       1,          0,           0],
+    ["Am",       1,          0,           0],
+    ["As",       1,          0,           0],
+    ["Aw",       1,          0,           0],
+    ["BSh",      1,          0,           0],
+    ["BSk",      0,          1,           0],
+    ["BWh",      1,          0,           0],
+    ["BWk",      0,          1,           0],
+    ["Cfa",      0,          1,           0],
+    ["Cfb",      0,          1,           0],
+    ["Cfc",      0,          1,           0],
+    ["Csa",      0,          1,           0],
+    ["Csb",      0,          1,           0],
+    ["Csc",      0,          1,           0],
+    ["Cwa",      0,          1,           0],
+    ["Cwb",      0,          1,           0],
+    ["Cwc",      0,          1,           0],
+    ["Dfa",      0,          1,           0],
+    ["Dfb",      0,          0,           1],
+    ["Dfc",      0,          0,           1],
+    ["Dfd",      0,          0,           1],
+    ["Dsa",      0,          0,           1],
+    ["Dsb",      0,          0,           1],
+    ["Dsc",      0,          0,           1],
+    ["Dsd",      0,          0,           1],
+    ["Dwa",      0,          0,           1],
+    ["Dwb",      0,          0,           1],
+    ["Dwc",      0,          0,           1],
+    ["Dwd",      0,          0,           1],
+    ["EF",       0,          0,           1],
+    ["ET",       0,          0,           1],
+]
+
+# Cross-walking to transform HRLC classes to generic PFTs for Amazon/Africa
+cwt_hrlc2gen_tropics = [
+    ["HRLC/GEN",   "TrBrEv", "TrBrDe", "TrNeEv", "TrNeDe", "ShBrEv", "ShBrDe", "ShNeEv", "ShNeDe", "NatGr", "Crops", "BS",  "Water", "SnowIce", "Urban", "No Data"],
+    ["NoData",      0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       1],
+    ["TrEvBr",      1,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrEvNe",      0,        0,        1,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrDeBr",      0,        1,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrDeNe",      0,        0,        0,        1,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["ShEv",        1,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["ShDe",        0,        1,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["Grass",       0,        0,        0,        0,        0,        0,        0,        0,        1,       0,       0,     0,       0,        0,       0],
+    ["Crop",        0,        0,        0,        0,        0,        0,        0,        0,        0,       1,       0,     0,       0,        0,       0],
+    ["Wood/flood",  0.8,      0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0.2,     0,        0,       0],
+    ["Grass/flood", 0,        0,        0,        0,        0,        0,        0,        0,        0.8,     0,       0,     0.2,     0,        0,       0],
+    ["LichMoss",    0,        0,        0,        0,        0,        0,        0,        0,        0.85,    0,       0.15,  0,       0,        0,       0],
+    ["Bare",        0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       1,     0,       0,        0,       0],
+    ["BuiltUp",     0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        1,       0],
+    ["WaterSeason", 0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     1,       0,        0,       0],
+    ["WaterPerm",   0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     1,       0,        0,       0],
+    ["SnowIce",     0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       1,        0,       0],
+]
+
+# Cross-walking to transform HRLC classes to generic PFTs for Siberia
+cwt_hrlc2gen_siberia = [
+    ["HRLC/GEN",   "TrBrEv", "TrBrDe", "TrNeEv", "TrNeDe", "ShBrEv", "ShBrDe", "ShNeEv", "ShNeDe", "NatGr", "Crops", "BS",  "Water", "SnowIce", "Urban", "No Data"],
+    ["NoData",      0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       1],
+    ["TrEvBr",      1,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrEvNe",      0,        0,        1,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrDeBr",      0,        1,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["TrDeNe",      0,        0,        0,        1,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["ShEv",        0,        0,        1,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["ShDe",        0,        0.25,     0,        0.75,     0,        0,        0,        0,        0,       0,       0,     0,       0,        0,       0],
+    ["Grass",       0,        0,        0,        0,        0,        0,        0,        0,        1,       0,       0,     0,       0,        0,       0],
+    ["Crop",        0,        0,        0,        0,        0,        0,        0,        0,        0,       1,       0,     0,       0,        0,       0],
+    ["Wood/flood",  0,        0,        0.8,      0,        0,        0,        0,        0,        0,       0,       0,     0.2,     0,        0,       0],
+    ["Grass/flood", 0,        0,        0,        0,        0,        0,        0,        0,        0.8,     0,       0,     0.2,     0,        0,       0],
+    ["LichMoss",    0,        0,        0,        0,        0,        0,        0,        0,        0.85,    0,       0.15,  0,       0,        0,       0],
+    ["Bare",        0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       1,     0,       0,        0,       0],
+    ["BuiltUp",     0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       0,        1,       0],
+    ["WaterSeason", 0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     1,       0,        0,       0],
+    ["WaterPerm",   0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     1,       0,        0,       0],
+    ["SnowIce",     0,        0,        0,        0,        0,        0,        0,        0,        0,       0,       0,     0,       1,        0,       0],
+]
+
+# Cross-walking table to transform generic PFTs to ORCHIDEE PFTs
+# Change from 20/03/2022 : Tr(Sh)NeEv for tropical and arid_warm from PFT2 to PFT4
+cwt_gen2orc = [
+    ["GENERIC", "KG/Still/LUH",  "PFT1", "PFT2", "PFT3", "PFT4", "PFT5", "PFT6", "PFT7", "PFT8", "PFT9", "PFT10", "PFT11", "PFT12", "PFT13", "PFT14", "PFT15"],
+    ["TrBrEv",  "tropical",       0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "arid_warm",      0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "arid_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "temp_warm",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "temp_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "boreal_warm",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "temp_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "boreal_warm",    0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      1,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "tropical",       0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "arid_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "arid_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "temp_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "temp_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "boreal_warm",    0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "temp_cool",      0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "boreal_warm",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "tropical",       0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "arid_warm",      0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "arid_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "temp_warm",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "temp_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "boreal_warm",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "temp_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "boreal_warm",    0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      1,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "tropical",       0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "arid_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "arid_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "temp_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "temp_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "boreal_warm",    0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "temp_cool",      0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "boreal_warm",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0],
+    ["NatGr",   "tropical/C3",    0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       1,       0],
+    ["NatGr",   "arid_warm/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       1,       0],
+    ["NatGr",   "arid_cool/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0],
+    ["NatGr",   "temp_warm/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0],
+    ["NatGr",   "temp_cool/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0],
+    ["NatGr",   "boreal_warm/C3", 0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       1],
+    ["NatGr",   "boreal_cool/C3", 0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       1],
+    ["NatGr",   "C4",             0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       1,       0,       0,       0,       0],
+    ["Crops",   "luhC3",          0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       1,       0,       0,       0],
+    ["Crops",   "luhC4",          0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       1,       0,       0],
+    ["BS",      "",               1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["Water",   "inland",         1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["SnowIce", "",               1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0],
+    ["Urban",   "tropical/C3",    0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0.2,     0],
+    ["Urban",   "arid_warm/C3",   0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0.2,     0],
+    ["Urban",   "arid_cool/C3",   0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0.2,     0,       0,       0,       0,       0],
+    ["Urban",   "temp_warm/C3",   0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0.2,     0,       0,       0,       0,       0],
+    ["Urban",   "temp_cool/C3",   0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0.2,     0,       0,       0,       0,       0],
+    ["Urban",   "boreal_warm/C3", 0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0.2],
+    ["Urban",   "boreal_cool/C3", 0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0.2],
+    ["Urban",   "C4",             0.8,    0,      0,      0,      0,      0,      0,      0,      0,      0,       0.2,     0,       0,       0,       0],
+]
+
+# Cross-walking table to transform generic PFTs to ORCHIDEE PFTs with additional urban type (PFT16)
+cwt_gen2orc16 = [
+    ["GENERIC", "KG/Still/LUH",  "PFT1", "PFT2", "PFT3", "PFT4", "PFT5", "PFT6", "PFT7", "PFT8", "PFT9", "PFT10", "PFT11", "PFT12", "PFT13", "PFT14", "PFT15", "PFT16"],
+    ["TrBrEv",  "tropical",       0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "arid_warm",      0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "arid_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "temp_warm",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "temp_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "boreal_warm",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "temp_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "boreal_warm",    0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrBrDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      1,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "tropical",       0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "arid_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "arid_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "temp_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "temp_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "boreal_warm",    0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "temp_cool",      0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "boreal_warm",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["TrNeDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "tropical",       0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "arid_warm",      0,      1,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "arid_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "temp_warm",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "temp_cool",      0,      0,      0,      0,      1,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "boreal_warm",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "temp_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "boreal_warm",    0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShBrDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      1,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "tropical",       0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "arid_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "arid_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "temp_warm",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "temp_cool",      0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "boreal_warm",    0,      0,      0,      1,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeEv",  "boreal_cool",    0,      0,      0,      0,      0,      0,      1,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "tropical",       0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "arid_warm",      0,      0,      1,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "arid_cool",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "temp_warm",      0,      0,      0,      0,      0,      1,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "temp_cool",      0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "boreal_warm",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["ShNeDe",  "boreal_cool",    0,      0,      0,      0,      0,      0,      0,      0,      1,      0,       0,       0,       0,       0,       0,       0],
+    ["NatGr",   "tropical/C3",    0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       1,       0,       0],
+    ["NatGr",   "arid_warm/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       1,       0,       0],
+    ["NatGr",   "arid_cool/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0,       0],
+    ["NatGr",   "temp_warm/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0,       0],
+    ["NatGr",   "temp_cool/C3",   0,      0,      0,      0,      0,      0,      0,      0,      0,      1,       0,       0,       0,       0,       0,       0],
+    ["NatGr",   "boreal_warm/C3", 0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       1,       0],
+    ["NatGr",   "boreal_cool/C3", 0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       1,       0],
+    ["NatGr",   "C4",             0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       1,       0,       0,       0,       0,       0],
+    ["Crops",   "luhC3",          0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       1,       0,       0,       0,       0],
+    ["Crops",   "luhC4",          0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       1,       0,       0,       0],
+    ["BS",      "",               1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["Water",   "inland",         1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["SnowIce", "",               1,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       0],
+    ["Urban",   "",               0,      0,      0,      0,      0,      0,      0,      0,      0,      0,       0,       0,       0,       0,       0,       1],
+]
+
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/gen2orc.py
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/gen2orc.py	(revision 7881)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/gen2orc.py	(revision 7881)
@@ -0,0 +1,177 @@
+import argparse, datetime, numpy as np, netCDF4, sys
+from grid import Grid
+from cwt import cwt_kg31to7, cwt_gen2orc, cwt_gen2orc16
+from basics import printStat, writePFTMAP
+CWT1 = cwt_kg31to7
+CWT2 = cwt_gen2orc # cwt_gen2orc for standard 15 PFTs, cwt_gen2orc16 for additional urbain PFT
+debug = False
+
+parser = argparse.ArgumentParser()
+parser.add_argument("year", action="store", type=int, help="year to process)")
+parser.add_argument("nlat", action="store", type=int, help="number of latitude steps (143 for LMDZ resolution)")
+parser.add_argument("-part", action="store", type=int, default=0, help="part (1 or 2 for nlat=14400 to split in two parts", required = False)
+args = parser.parse_args()
+year = args.year
+nlat = args.nlat
+part = args.part
+nlon = 144 if nlat == 143 else nlat * 2
+if part == 0: lat1=0; lat2=nlat
+elif part == 1: lat1=0; lat2=int(nlat/2)
+elif part == 2: lat1=int(nlat/2); lat2=nlat
+step = 180. / nlat
+luh_year0 = 1992
+esa = "v2.0.7b" if year <= 2015 else "v2.1.1"
+print("GEN2ORC for %s, %sx%s, %s:%s" % (year, nlat, nlon, lat1, lat2))
+
+atts = dict(description = "ORCHIDEE PFT map derived from ESACCI Land Cover map %(esa)s for %(year)s, aggregated using user-tool v4.3, cross-walking table v2.4, C3/C4 for grasses from C.Still (May 2018), C3/C4 for crops from LUH2 (2022)" % dict(esa = esa, year = year),
+            PFT1  = "bare soil",
+            PFT2  = "tropical broadleaf evergreen",
+            PFT3  = "tropical broadleaf raingreen",
+            PFT4  = "temperate needleleaf evergreen",
+            PFT5  = "temperate broadleaf evergreen",
+            PFT6  = "temperate broadleaf summergreen",
+            PFT7  = "boreal needleleaf evergreen",
+            PFT8  = "boreal broadleaf summergreen",
+            PFT9  = "boreal needleleaf summergreen",
+            PFT10 = "temperate C3 grass",
+            PFT11 = "C4 grass",
+            PFT12 = "C3 agriculture",
+            PFT13 = "C4 agriculture",
+            PFT14 = "tropical C3 grass",
+            PFT15 = "boreal C3 grass",
+            # PFT16 = "urban", # add this attribute for 16 PFTs (with additional urbain PFT)
+            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (peylin@lsce.ipsl.fr)",
+            web = "https://orchidas.lsce.ipsl.fr/dev/lccci/",
+            date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"))
+
+path_gen = "/home/surface5/vbastri/BDD/ESACCI-LC/GENERIC_PFT/%(step)s/ESACCI-LC-L4-LCCS-Map-300m-P1Y-aggregated-%(step).6fDeg-%(year)s-%(esa)s.nc" % dict(step = step, year = year, esa = esa)
+path_luh = "/home/surface5/vbastri/BDD/LUH2/TRENDY2022/LUH2v11_1992-2022_%(nlat)sx%(nlon)s.nc" % dict(nlat = nlat, nlon = nlon)
+path_kg  = "/home/surface5/vbastri/BDD/KG2017/KG_%(nlat)sx%(nlon)s.nc" % dict(nlat = nlat, nlon = nlon)
+path_c4  = "/home/surface5/vbastri/BDD/C4FRAC/C4frac_%(nlat)sx%(nlon)s.nc" % dict(nlat = nlat, nlon = nlon)
+path_wb  = "/home/surface5/vbastri/BDD/ESACCI-LC/WATER/WB_%(nlat)sx%(nlon)s.nc" % dict(nlat = nlat, nlon = nlon)
+path_out = "PFTmap_%(year)s_%(nlat)sx%(nlon)s%(end)s.nc" % dict(year = year, nlat = nlat, nlon = nlon, end = "_%(lat1)s-%(lat2)s" % dict(lat1 = lat1, lat2 = lat2) if part!=0 else "")
+GEN_LIST = ["Tree_Broadleaf_Evergreen", "Tree_Broadleaf_Deciduous", "Tree_Needleleaf_Evergreen", "Tree_Needleleaf_Deciduous", "Shrub_Broadleaf_Evergreen", "Shrub_Broadleaf_Deciduous", "Shrub_Needleleaf_Evergreen", "Shrub_Needleleaf_Deciduous", "Natural_Grass", "Crops", "Bare_soil", "Water", "Snow_Ice", "Urban", "NoData"]
+GEN_KEYS = ["TrBrEv", "TrBrDe", "TrNeEv", "TrNeDe", "ShBrEv", "ShBrDe", "ShNeEv", "ShNeDe", "NatGr", "Crops", "BS",  "Water", "SnowIce", "Urban", "NoData"]
+KG_KEYS = CWT1[0][1:]
+ORC_KEYS = CWT2[0][2:]
+
+# Function to safely divide two streams: if denominator equals zero, division result is set to zero
+# During division zero values in denominator are changed to 1 so that python does not warn on errors
+def divide(num, den): return np.where(den == 0, 0, num / np.where(den == 0, 1, den))
+
+# Function to make plots for data streams
+def plot(**kwargs):
+    for name, data in kwargs.items():
+        for idx, item in enumerate(data if len(data.shape)>2 else [data]):
+            grid.plotmap(np.ma.masked_where(item==0, item), dpi=150, vmin=0, vmax=1,
+                         filename="%s%s.png" % (name.lower(), idx if len(data.shape)>2 else ""), title="%s%s" % (name.upper(), idx if len(data.shape)>2 else ""),
+                         info="min=%.3g\nmax=%.3g\nmean=%.3g\nnzpts=%s" % (item.min(), item.max(), grid.meanmask(np.ma.masked_where(item==0, item)), (item>0).sum()))
+
+# Read and normalize generic PFTmap
+print("\nREAD GENERIC PFTMAP")
+nc = netCDF4.Dataset(path_gen)
+lat = nc.variables["lat"][lat1:lat2]
+lon = nc.variables["lon"][:]
+gen = np.zeros((len(GEN_LIST), len(lat), len(lon)))
+for idx, item in enumerate(GEN_LIST): gen[idx] = nc.variables[item][lat1:lat2].astype("f8").filled(0)
+nc.close()
+printStat(gen=gen)
+gen = divide(gen, gen.sum(axis=0))
+printStat(genn=gen)
+if debug: grid = Grid(lat=lat, lon=lon)
+if debug: plot(gen=gen, gensum=gen.sum(axis=0))
+
+# Read and normalize Koeppen-Geiger map with 31 classes
+print("\nREAD KG")
+nc = netCDF4.Dataset(path_kg)
+if not np.allclose(lat, nc.variables["lat"][lat1:lat2]) or not np.allclose(lon, nc.variables["lon"][:]): print("Different grid"); sys.exit()
+kg31 = nc.variables["kg"][:,lat1:lat2].astype("f8").filled(0)
+nc.close()
+printStat(kg31=kg31)
+kg31 = divide(kg31, kg31.sum(axis=0))
+printStat(kg31n=kg31)
+
+# Aggregate KG into 7 groups
+print("Check CWT1 kg31to7 = %s" % np.all([np.array(item[1:]).sum() == 1 for item in CWT1[1:]]))
+kg = np.zeros((len(KG_KEYS), len(lat), len(lon)))
+for idx, rule in enumerate(CWT1[1:]): kg = kg + kg31[idx] * np.array(rule[1:])[:,None,None]
+printStat(kg=kg)
+if debug: plot(kg=kg)
+# Fill non-defined values in KG with the dominant class across each latitude (or last class if no data in a latitude)
+domclass = np.where(kg.sum(axis=(0,2)) == 0, len(kg)-1, kg.sum(axis=2).argmax(axis=0))
+for idx in range(len(lat)): kg[domclass[idx], idx, kg[:,idx,:].sum(axis=0) == 0] = 1
+printStat(kgf=kg)
+if debug: plot(kgf=kg)
+
+print("\nREAD C4")
+nc = netCDF4.Dataset(path_c4)
+if not np.allclose(lat, nc.variables["lat"][lat1:lat2]) or not np.allclose(lon, nc.variables["lon"][:]): print("Different grid"); sys.exit()
+c4 = nc.variables["C4frac"][lat1:lat2].astype("f8").filled(0)
+nc.close()
+printStat(c4=c4)
+if debug: plot(c4=c4)
+
+print("\nREAD WB")
+nc = netCDF4.Dataset(path_wb)
+if not np.allclose(lat, nc.variables["lat"][lat1:lat2]) or not np.allclose(lon, nc.variables["lon"][:]): print("Different grid"); sys.exit()
+wb = nc.variables["wb"][:,lat1:lat2].astype("f8").filled(0)
+nc.close()
+printStat(wb=wb)
+wb = divide(wb, wb.sum(axis=0))
+printStat(wbn=wb)
+ocean = np.where(wb[0] > 0.999, 1, 0) # ocean 0/1-mask
+coast = np.where((wb[0] > 0) & (wb[0] <= 0.999), 1, 0) # coast 0/1-mask
+inland = np.where(wb[0] == 0, wb[2], 0) # inland water fraction excluding coasts
+print("ocean_npts=%s, coast_npts=%s, inland_npts=%s" % (ocean.sum(), coast.sum(), (inland>0).sum()))
+if debug: plot(wb=wb, ocean=ocean, coast=coast, inland=inland)
+
+# read C4 map from LUH2 for crops
+print("\nREAD LUH")
+nc = netCDF4.Dataset(path_luh)
+if not np.allclose(lat, nc.variables["lat"][lat1:lat2]) or not np.allclose(lon, nc.variables["lon"][:]): print("Different grid"); sys.exit()
+c3ann = nc.variables["c3ann"][year-luh_year0, lat1:lat2].astype("f8").filled(0)
+c4ann = nc.variables["c4ann"][year-luh_year0, lat1:lat2].astype("f8").filled(0)
+c3per = nc.variables["c3per"][year-luh_year0, lat1:lat2].astype("f8").filled(0)
+c4per = nc.variables["c4per"][year-luh_year0, lat1:lat2].astype("f8").filled(0)
+c3nfx = nc.variables["c3nfx"][year-luh_year0, lat1:lat2].astype("f8").filled(0)
+total = c3ann + c4ann + c3per + c4per + c3nfx
+luhc4 = divide(c4ann + c4per, total)
+nc.close()
+printStat(luhc4=luhc4)
+if debug: plot(luhc4=luhc4)
+
+# aggregate generic PFTs to ORCHIDEE PFTs
+print("\nAGGREGATE GENERIC PFTs to ORCHIDEE PFTs")
+print("Check CWT2 gen2orc = %s" % np.all([np.array(item[2:]).sum() == 1 for item in CWT2[1:]]))
+orc = np.zeros((len(ORC_KEYS), len(lat), len(lon)))
+for rule in CWT2[1:]:
+    idx = GEN_KEYS.index(rule[0])
+    cond = rule[1]
+    fracs = np.array(rule[2:])
+    if cond in KG_KEYS: coef = kg[KG_KEYS.index(cond)]
+    elif cond in ["C3", "C4"]: coef = c4 if cond == "C4" else 1-c4
+    elif cond in ["luhC3", "luhC4"]: coef = luhc4 if cond == "luhC4" else 1-luhc4
+    elif "/" in cond: coef = kg[KG_KEYS.index(cond.split("/")[0])] * (c4 if cond.split("/")[1] == "C4" else 1-c4)
+    elif cond == "inland": coef = inland
+    elif cond == "": coef = 1
+    else: print("Unknown condition: %s" % cond); sys.exit()
+    orc = orc + gen[idx] * coef * fracs[:,None,None]
+printStat(orc=orc)
+
+# set small fractions to zero (< 0.001)
+orc = np.where(orc < 0.001, 0, orc)
+# set 1 to bare soil on the coast if no other PFTs
+orc[0] = np.where((coast == 1) & (orc[1:].sum(axis=0) < 0.001), 1, orc[0])
+# normalize
+orc = divide(orc, orc.sum(axis=0))
+# apply land-sea mask
+orc = np.ma.array(orc)
+orc.mask = (ocean == 1)
+# filter out Antarctica (lat < S60)
+orc[:, lat < -60, :] = np.ma.masked
+
+printStat(orcn=orc)
+if debug: plot(orc=orc, orcsum=orc.sum(axis=0))
+
+writePFTMAP(orc, lat=lat, lon=lon, year=year, filename=path_out, atts=atts)
+
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/mergeLUH.py
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/mergeLUH.py	(revision 7881)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/mergeLUH.py	(revision 7881)
@@ -0,0 +1,237 @@
+import datetime, numpy as np, netCDF4, sys
+from grid import Grid
+from cwt import cwt_kg31to3
+from basics import printStat, writePFTMAP
+debug = False
+
+path_orc = "/home/surface5/vbastri/BDD/ESACCI-LC/ORCHIDEE_NEW/0.25/PFTmap_2010.nc"
+path_luh = "/home/surface5/vbastri/BDD/LUH2/GCB2022/states.nc"
+path_kg  = "/home/surface5/vbastri/BDD/KG2017/KG_720x1440.nc"
+path_c4  = "/home/surface5/vbastri/BDD/C4FRAC/C4frac_720x1440.nc"
+path_wb  = "/home/surface5/vbastri/BDD/ESACCI-LC/WATER/WB_720x1440.nc"
+path_out = "netcdf/PFTmap_%(year)s.nc"
+luh_year0 = 850
+merge_year = 2010
+build_years = np.arange(850, 2022+1)
+
+atts = dict(description = "ORCHIDEE PFT map derived from ESACCI Land Cover map v2.0.7b for 2010, aggregated using user-tool v4.3, cross-walking table v2.4, C3/C4 for grasses from C.Still (May 2018), C3/C4 for crops from LUH2-GCB2022, followed by historical reconstrution using G. Hurtt data LUH2-GCB2022 with merging rules v2.0",
+            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (peylin@lsce.ipsl.fr)",
+            web = "https://orchidas.lsce.ipsl.fr/dev/lccci/",
+            date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"))
+
+# ORCHIDEE PFT indices
+BS    = 0
+TREBF = 1
+TRDBF = 2
+TEENF = 3
+TEEBF = 4
+TEDBF = 5
+BOENF = 6
+BODBF = 7
+BODNF = 8
+TEGR  = 9
+C4GR  = 10
+C3CR  = 11
+C4CR  = 12
+TRGR  = 13
+BOGR  = 14
+
+# LUH class indices
+PRIMF = 0
+PRIMN = 1
+SECDF = 2
+SECDN = 3
+URBAN = 4
+C3ANN = 5
+C4ANN = 6
+C3PER = 7
+C4PER = 8
+C3NFX = 9
+PASTR = 10
+RANGE = 11
+ICWTR = 12
+LUH_LIST = ["primf", "primn", "secdf", "secdn", "urban", "c3ann", "c4ann", "c3per", "c4per", "c3nfx", "pastr", "range"]
+
+# Function to safely divide two streams: if denominator equals zero, division result is set to zero
+# During division denominator zeros are changed to units, so that python does not warn on errors
+def divide(num, den): return np.where(den == 0, 0, num / np.where(den == 0, 1, den))
+
+# Function to make plots for data streams
+def plot(**kwargs):
+    for name, data in kwargs.items():
+        for idx, item in enumerate(data if len(data.shape)>2 else [data]):
+            grid.plotmap(np.ma.masked_where(item==0, item), dpi=150, vmin=0, vmax=1,
+                         filename="%s%s.png" % (name.lower(), idx if len(data.shape)>2 else ""), title="%s%s" % (name.upper(), idx if len(data.shape)>2 else ""),
+                         info="min=%.3g\nmax=%.3g\nmean=%.3g\nnzpts=%s" % (item.min(), item.max(), grid.meanmask(np.ma.masked_where(item==0, item)), (item>0).sum()))
+
+# Function to extend land by 1 pixel outside land-sea mask
+def extend(veg):
+    newveg = np.ma.masked_all((len(veg), len(lat), len(lon)))
+    loc_area = np.where(veg[0].mask == False, grid.area, 0)
+    wveg = veg * loc_area
+    for ilat in range(len(lat)):
+        for ilon in range(len(lon)):
+            if veg[0,ilat,ilon] is not np.ma.masked:
+                newveg[:,ilat,ilon] = veg[:,ilat,ilon]
+            else:
+                if wveg[:, max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].count() == 0: continue
+                newveg[:,ilat,ilon] = wveg[:, max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].sum(axis=-1).sum(axis=-1) / loc_area[max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].sum()
+    return newveg
+
+# Read present-day ORCHIDEE map to be merged with LUH
+print("\nREAD ORCHIDEE")
+nc = netCDF4.Dataset(path_orc)
+lat = nc.variables["lat"][:]
+lon = nc.variables["lon"][:]
+orc = nc.variables["maxvegetfrac"][0].astype("f8").filled(0.)
+nc.close()
+printStat(orc=orc)
+orc = divide(orc, orc.sum(axis=0))
+printStat(orcn=orc)
+grid = Grid(lat=lat, lon=lon)
+if debug: plot(orc=orc, orcsum=orc.sum(axis=0))
+
+# Read Koeppen-Geiger map with 31 classes
+print("\nREAD KG")
+nc = netCDF4.Dataset(path_kg)
+if not np.all(lat == nc.variables["lat"][:]) or not np.all(lon == nc.variables["lon"][:]): print("Different grid"); sys.exit()
+kg31 = nc.variables["kg"][:].astype("f8").filled(0)
+nc.close()
+printStat(kg31=kg31)
+kg31 = divide(kg31, kg31.sum(axis=0))
+printStat(kg31n=kg31)
+
+# Aggregate KG into 3 groups
+kg = np.zeros((len(cwt_kg31to3[0])-1, len(lat), len(lon)))
+for idx, rule in enumerate(cwt_kg31to3[1:]): kg = kg + kg31[idx] * np.array(rule[1:])[:,None,None]
+printStat(kg3=kg)
+if debug: plot(kg=kg)
+# Fill KG where undefined with the dominant class across each latitude (or last class if no data in a latitude)
+domclass = np.where(kg.sum(axis=(0,2)) == 0, len(kg)-1, kg.sum(axis=2).argmax(axis=0))
+for idx in range(len(lat)): kg[domclass[idx], idx, kg[:,idx,:].sum(axis=0) == 0] = 1
+printStat(kg3f=kg)
+if debug: plot(kgf=kg)
+
+print("\nREAD C4")
+nc = netCDF4.Dataset(path_c4)
+if not np.all(lat == nc.variables["lat"][:]) or not np.all(lon == nc.variables["lon"][:]): print("Different grid"); sys.exit()
+c4 = nc.variables["C4frac"][:].astype("f8").filled(0)
+nc.close()
+printStat(c4=c4)
+if debug: plot(c4=c4)
+
+print("\nREAD WB")
+nc = netCDF4.Dataset(path_wb)
+if not np.all(lat == nc.variables["lat"][:]) or not np.all(lon == nc.variables["lon"][:]): print("Different grid"); sys.exit()
+wb = nc.variables["wb"][:].astype("f8").filled(0)
+nc.close()
+printStat(wb=wb)
+lsm = np.where(wb[1] > 0.001, 1, 0)
+print("land pts=%s" % lsm.sum())
+if debug: plot(wb=wb, lsm=lsm)
+
+# Read LUH for the merging year
+print("\nREAD LUH")
+nc = netCDF4.Dataset(path_luh)
+if not np.all(lat == nc.variables["lat"][:]) or not np.all(lon == nc.variables["lon"][:]): print("Different grid"); sys.exit()
+luh = np.zeros((len(LUH_LIST)+1, len(lat), len(lon))) # array for all LUH maps plus ICWTR
+for idx, item in enumerate(LUH_LIST): luh[idx] = nc.variables[item][merge_year-luh_year0].astype("f8") # do not fill with 0 to calculate icwtr over land (ice and water)
+luh[ICWTR] = np.clip(1 - luh.sum(axis=0), 0, 1)
+luh = np.ma.masked_invalid(luh).filled(0)
+luh = divide(luh, luh.sum(axis=0))
+nc.close()
+printStat(luh=luh)
+if debug: plot(luh=luh, luhsum=luh.sum(axis=0))
+
+onlyorc = np.where((orc.sum(axis=0) > 0) & (luh.sum(axis=0) == 0), 1, 0)
+onlyluh = np.where((orc.sum(axis=0) == 0) & (luh.sum(axis=0) > 0), 1, 0)
+print("ORC>0 LUH=0", onlyorc.sum())
+print("ORC=0 LUH>0", onlyluh.sum())
+if debug: plot(onlyorc=onlyorc, onlyluh=onlyluh)
+
+print("\nBUILD FRACTIONS")
+luh_C3  = luh[C3ANN] + luh[C3PER] + luh[C3NFX]
+luh_C4  = luh[C4ANN] + luh[C4PER]
+luh_ant = luh_C3 + luh_C4 + luh[PASTR]
+
+orc_gr = orc[TEGR] + orc[TRGR] + orc[BOGR] + orc[C4GR]
+orc_cr = orc[C3CR] + orc[C4CR]
+gr_ant = np.clip(luh_ant - orc_cr, 0, np.inf)
+gr_nat = np.clip(orc_gr  - gr_ant, 0, np.inf)
+gr_tr = gr_nat * divide(orc[TRGR], orc_gr)
+gr_te = gr_nat * divide(orc[TEGR], orc_gr)
+gr_bo = gr_nat * divide(orc[BOGR], orc_gr)
+gr_C4 = gr_nat * divide(orc[C4GR], orc_gr)
+total_nat = orc[BS] + orc[TREBF] + orc[TRDBF] + orc[TEENF] + orc[TEEBF] + orc[TEDBF] + orc[BOENF] + orc[BODBF] + orc[BODNF] + gr_tr + gr_te + gr_bo + gr_C4
+
+frac = divide(orc, total_nat)
+frac[TRGR] = divide(gr_tr, total_nat)
+frac[TEGR] = divide(gr_te, total_nat)
+frac[BOGR] = divide(gr_bo, total_nat)
+frac[C4GR] = divide(gr_C4, total_nat)
+frac[C3CR] = 0
+frac[C4CR] = 0
+printStat(frac=frac)
+print("ORC>0 FRAC=0", ((orc.sum(axis=0) > 0) & (frac.sum(axis=0) == 0)).sum()) # these are pixels with crops in ORCHIDEE (161)
+print("LUH>0 FRAC=0", ((luh.sum(axis=0) > 0) & (frac.sum(axis=0) == 0)).sum()) # these are mostly ICWTR pixels in LUH, few as PRIMN (361)
+print("ORC>0 LUH>0 FRAC=0", ((orc.sum(axis=0) > 0) & (luh.sum(axis=0) > 0) & (frac.sum(axis=0) == 0)).sum()) # 40
+if debug: plot(frac=frac, fracsum=frac.sum(axis=0))
+
+nc = netCDF4.Dataset(path_luh)
+for year in build_years:
+    print("\nREAD LUH FOR YEAR %s" % year)
+    yluh = np.zeros((len(LUH_LIST)+1, len(lat), len(lon)))
+    for idx, item in enumerate(LUH_LIST): yluh[idx] = nc.variables[item][year-luh_year0].astype("f8")
+    yluh[ICWTR] = np.clip(1 - yluh.sum(axis=0), 0, 1)
+    yluh = np.ma.masked_invalid(yluh).filled(0)
+    yluh = divide(yluh, yluh.sum(axis=0))
+    printStat(yluh=yluh)
+
+    yluh_nat = yluh[PRIMF] + yluh[PRIMN] + yluh[SECDF] + yluh[SECDN] + yluh[RANGE] + yluh[URBAN] + yluh[ICWTR]
+    yluh_C3  = yluh[C3ANN] + yluh[C3PER] + yluh[C3NFX]
+    yluh_C4  = yluh[C4ANN] + yluh[C4PER]
+    
+    print("\nBUILD ORC FOR YEAR %s" % year)    
+    yorc = frac * yluh_nat
+    yorc[TRGR] = yorc[TRGR] + kg[0] * (1 - c4) * yluh[PASTR]
+    yorc[TEGR] = yorc[TEGR] + kg[1] * (1 - c4) * yluh[PASTR]
+    yorc[BOGR] = yorc[BOGR] + kg[2] * (1 - c4) * yluh[PASTR]
+    yorc[C4GR] = yorc[C4GR] + c4 * yluh[PASTR]
+    yorc[C3CR] = yluh_C3
+    yorc[C4CR] = yluh_C4
+    printStat(yorc=yorc)
+    print("0<YORC<0.999", ((yorc.sum(axis=0)>0) & (yorc.sum(axis=0)<0.999)).sum())
+    print("ORC>0 LUH>0 YORC<0.999", ((orc.sum(axis=0)>0) & (luh.sum(axis=0)>0) & (yorc.sum(axis=0)<0.999)).sum())
+    print("ORC>0 LUH=0 YORC=0", ((orc.sum(axis=0)>0) & (luh.sum(axis=0)==0) & (yorc.sum(axis=0)==0)).sum())
+    print("ORC>0 YORC<0.999", ((orc.sum(axis=0)>0) & (yorc.sum(axis=0)<0.999)).sum())
+    if debug: plot(yorc=yorc)
+    
+    # POST PROCESS
+    # - drop pixels where ORC has no info
+    # - use data from ORC where LUH has no info or problematic balance
+    # - drop fractions below 0.001
+    # - normalize to 1
+    # - if OCEAN>0.999 apply mask
+    # - filter out Antarctica (lat < S60)
+    yorc = np.where((orc.sum(axis=0)==0), 0, yorc)
+    yorc = np.where((orc.sum(axis=0)>0) & (yorc.sum(axis=0)<0.999), orc, yorc)
+    yorc = np.where(yorc < 0.001, 0, yorc)
+    yorc = divide(yorc, yorc.sum(axis=0))
+    yorc = np.ma.where(wb[0] < 0.999, yorc, np.ma.masked)
+    yorc[:, lat < -60, :] = np.ma.masked
+    printStat(orcn=yorc)
+    if debug: plot(yorcn=yorc, yorcnsum=yorc.sum(axis=0))
+
+    # Check that resulting LSM is identical to ORC original
+    m1 = np.where(orc.sum(axis=0)>0, 1, 0)
+    m2 = np.where(yorc.sum(axis=0)>0, 1, 0)
+    print(np.all(m1==m2))
+    
+    # Extend by 1 pixel
+    yorc = extend(yorc)
+    printStat(orcext=yorc)
+    
+    writePFTMAP(yorc, lat=lat, lon=lon, year=year, filename=path_out % dict(year=year), atts=atts)
+
+nc.close()
+
