Index: /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/README.txt
===================================================================
--- /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/README.txt	(revision 6267)
+++ /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/README.txt	(revision 6267)
@@ -0,0 +1,37 @@
+Reconstruction of the historical PFT maps is done with the following steps:
+
+1. Aggregation of the ESA CCI land cover data for 2010 into generic PFT map (0.25 degree resolution, 15 PFT classes), using:
+   - ESA CCI land cover map for 2010 (v2.0.7b, 300 meters resolution, 38 land cover classes)
+   - aggregation tool "lc-user-tools" (v3.13)
+   - additional Koeppen-Geiger climate zone map
+   - default cross-walking table (v2.4), matching ESA CCI classes with generic PFT classes
+   - additional cross-walking table (v2.4), refining the conversion of several ESA CCI classes depending on the climate zone
+
+2. Aggregation of the generic PFT classes for 2010 into 15 ORCHIDEE PFTs plus 4 additional classes (urban, inland water, coastal water, permanent snow/ice), using:
+   - python script "CCI_to_PFT19.py"
+   - additional partitioning C3/C4 vegetation map (Still et al., 2009)
+   - Earth surface topography data at 5 minute resolution (ETOPO5.nc)
+
+3. PFT map reconstruction for by merging ORCHIDEE PFT map for 2010 with LUH2v2h or scenarios data:
+   One want to use LUHv2 for the evolution of crops and pasture lands. On the other side, one want to benefit of the ESA map for defining PFT mosaic within each PFT. Based on reference year (here we use 2010), we substract the fractions of the LUH2 that we identify as managed to the fractions of crops and grasslands in the ESA map. The remaining land fractions in Ethe SA map are assumed to correspond to the mosaic of vegetation which forms the natural lands of each tile. From this. This PFT fractions related to the natural land is then used to split, each year, the total lands of LUH2 corresponding to the sum of Primary forested land, Primary non-forested land, Secondary forested land and Secondary non-forested land. The remaining vegetation classes from LUH2 are more easily allocated into C3 grass, C4 grass, C3 crop and C4 crop. C3/C4 grass split uses the information contained into the ESA map. 
+
+  For historical period (850-2014), use of the following script and data
+   - ferret script "process_landusemap_19PFT_quarterdeg_final_temp.jnl"
+   - LUH data : /home/basc/dnarayan/LandUsemap_ORCHIDEE/LUH2v2h/states.nc : they are the same than those stored on input4mip
+   - ESA: /home/satellites7/vbastri/LC/orchidee/ESA2.0.7b/19PFT_quarterdeg_final/PFTmap_2010.nc
+   
+  For the scenarios (2015-2100)
+   - process_landusemap_19PFT_quarterdeg_final_temp_future.jnl : generation of a reference map for year 2010 with the PFT fractions corresponding to the natural vegetation (fractions_from_ESA2.0.7b_19PFT_quarterdeg_PFTmap_2010.nc). This is the same that the one used internally in the script process_landusemap_19PFT_quarterdeg_final_temp.jnl for the historical period.
+   - process_landusemap_19PFT_quarterdeg_final_temp_future_2dstep.jnl : generation of the land-use map for a specific scenario given as argument to the script (for instance, AIM-ssp370-2-1-f or GCAM-ssp434-2-1-f) 
+   - ScenarioMIP data : states variables for the different scenarios (SSP-RCP) stored under /home/surface7/vuichard/LandUseMap_CMIP6/(they are the same than those stored on input4mip)
+   - ESA: /home/satellites7/vbastri/LC/orchidee/ESA2.0.7b/19PFT_quarterdeg_final/PFTmap_2010.nc
+
+
+4. Post-processing of the historical maps (adding inland water as bare soil; extension of land points into the ocean by 8 pixels; formatting netcdf-files according to ORCHIDEE requirements), using:
+   - python script "post_PFT.py"
+
+ESA CCI land cover data can be found at https://www.esa-landcover-cci.org/
+All other used data and tools can be found in PFTMAPS/SOURCE directory
+To run python scripts Python 2.7 or higher with numpy and netCDF4 libraries is required
+For additional information please visit https://orchidas.lsce.ipsl.fr/dev/LCCCI.php
+
Index: /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future.jnl
===================================================================
--- /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future.jnl	(revision 6267)
+++ /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future.jnl	(revision 6267)
@@ -0,0 +1,148 @@
+! This is a map with 15 tiles
+! The first 13 tiles correspond to those considered in ORCHIDEE
+! The 14th is inland water
+! The 15th us urban land
+use "/home/satellites7/vbastri/LC/orchidee/ESA2.0.7b/19PFT_quarterdeg_final/PFTmap_2010.nc"
+
+! Renormalized maxvegetfrac over land area excluding inland water
+let mymaxvegetfrac=maxvegetfrac
+
+! This is the reference states file from LUH2v2h
+! from 850 to 2015
+use "/home/basc/dnarayan/LandUsemap_ORCHIDEE/LUH2v2h/states.nc"
+
+
+let primf_cor=primf[l=1161,d=2]
+let primn_cor=primn[l=1161,d=2]
+let secdf_cor=secdf[l=1161,d=2]
+let secdn_cor=secdn[l=1161,d=2]
+let urban_cor=urban[l=1161,d=2]
+let c3ann_cor=c3ann[l=1161,d=2]
+let c4ann_cor=c4ann[l=1161,d=2]
+let c3per_cor=c3per[l=1161,d=2]
+let c4per_cor=c4per[l=1161,d=2]
+let c3nfx_cor=c3nfx[l=1161,d=2]
+let pastr_cor=pastr[l=1161,d=2]
+let range_cor=range[l=1161,d=2]
+
+
+define axis/t="1-JAN-2010:00:00":"1-JAN-2010:00:00":1/calendar=noleap myaxis
+
+! We compute the total fraction of anthropogenic tiles in LUH2 that 
+! we can atttribute to specific PFTs into ORCHIDEE
+! crop is either C3 or C4 crops, pasture is grasslands
+! on the opposite, there is no correspondance for urban or rangeland, 
+! so, these tiles are not accounted for, here.
+let anthropogenic=pastr_cor+c3ann_cor+c4ann_cor+c3per_cor+c4per_cor+c3nfx_cor
+
+! pft10z11z14z15 will correspond to the fraction of 'natural grasslands'
+let pft10z11z14z15=MAX(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn]-(anthropogenic[l=1]-mymaxvegetfrac[k=12:13@sum,d=1,gt=myaxis@asn]),0.)
+! fraction of PFT10 and PFT11 (fraction of natural C3 and C4 grasslands) is obtained by splitting pft10z11z12z14z15 
+! using the relative contribution of C3grass and C4grass in the ESA-CCI map
+let pft10=if((mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn]) gt 0.) then mymaxvegetfrac[k=10,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])*pft10z11z14z15 else 0.
+let pft11=if((mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn]) gt 0.) then mymaxvegetfrac[k=11,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])*pft10z11z14z15 else 0.
+let pft14=if((mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn]) gt 0.) then mymaxvegetfrac[k=14,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])*pft10z11z14z15 else 0.
+let pft15=if((mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn]) gt 0.) then mymaxvegetfrac[k=15,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])*pft10z11z14z15 else 0.
+
+
+! We renormalize fractions for PFT1 to PFT9 and for variables pft10 and pft11
+! to get the sum equal to 1
+let pftnat=mymaxvegetfrac[k=1:9@sum,d=1,gt=myaxis@asn]+pft10+pft11+pft14+pft15
+let frac_pft1=mymaxvegetfrac[k=1,d=1,gt=myaxis@asn]/pftnat
+let frac_pft2=mymaxvegetfrac[k=2,d=1,gt=myaxis@asn]/pftnat
+let frac_pft3=mymaxvegetfrac[k=3,d=1,gt=myaxis@asn]/pftnat
+let frac_pft4=mymaxvegetfrac[k=4,d=1,gt=myaxis@asn]/pftnat
+let frac_pft5=mymaxvegetfrac[k=5,d=1,gt=myaxis@asn]/pftnat
+let frac_pft6=mymaxvegetfrac[k=6,d=1,gt=myaxis@asn]/pftnat
+let frac_pft7=mymaxvegetfrac[k=7,d=1,gt=myaxis@asn]/pftnat
+let frac_pft8=mymaxvegetfrac[k=8,d=1,gt=myaxis@asn]/pftnat
+let frac_pft9=mymaxvegetfrac[k=9,d=1,gt=myaxis@asn]/pftnat
+let frac_pft10=pft10/pftnat
+let frac_pft11=pft11/pftnat
+let frac_pft14=pft14/pftnat
+let frac_pft15=pft15/pftnat
+
+
+! We want a vector from k=1:13, not one scalar for each pft 
+let test_frac_pft1=k[k=1:15]/k[k=1:15]*frac_pft1
+let test_frac_pft2=k[k=1:15]/k[k=1:15]*frac_pft2
+let test_frac_pft3=k[k=1:15]/k[k=1:15]*frac_pft3
+let test_frac_pft4=k[k=1:15]/k[k=1:15]*frac_pft4
+let test_frac_pft5=k[k=1:15]/k[k=1:15]*frac_pft5
+let test_frac_pft6=k[k=1:15]/k[k=1:15]*frac_pft6
+let test_frac_pft7=k[k=1:15]/k[k=1:15]*frac_pft7
+let test_frac_pft8=k[k=1:15]/k[k=1:15]*frac_pft8
+let test_frac_pft9=k[k=1:15]/k[k=1:15]*frac_pft9
+let test_frac_pft10=k[k=1:15]/k[k=1:15]*frac_pft10
+let test_frac_pft11=k[k=1:15]/k[k=1:15]*frac_pft11
+let test_frac_pft14=k[k=1:15]/k[k=1:15]*frac_pft14
+let test_frac_pft15=k[k=1:15]/k[k=1:15]*frac_pft15
+
+let test2_frac_pft1=if(k[k=1:15] eq 1) then test_frac_pft1 else 0*test_frac_pft1
+let test2_frac_pft2=if(k[k=1:15] eq 2) then test_frac_pft2 else 0*test_frac_pft2
+let test2_frac_pft3=if(k[k=1:15] eq 3) then test_frac_pft3 else 0*test_frac_pft3
+let test2_frac_pft4=if(k[k=1:15] eq 4) then test_frac_pft4 else 0*test_frac_pft4
+let test2_frac_pft5=if(k[k=1:15] eq 5) then test_frac_pft5 else 0*test_frac_pft5
+let test2_frac_pft6=if(k[k=1:15] eq 6) then test_frac_pft6 else 0*test_frac_pft6
+let test2_frac_pft7=if(k[k=1:15] eq 7) then test_frac_pft7 else 0*test_frac_pft7
+let test2_frac_pft8=if(k[k=1:15] eq 8) then test_frac_pft8 else 0*test_frac_pft8
+let test2_frac_pft9=if(k[k=1:15] eq 9) then test_frac_pft9 else 0*test_frac_pft9
+let test2_frac_pft10=if(k[k=1:15] eq 10) then test_frac_pft10 else 0*test_frac_pft10
+let test2_frac_pft11=if(k[k=1:15] eq 11) then test_frac_pft11 else 0*test_frac_pft11
+let test2_frac_pft14=if(k[k=1:15] eq 14) then test_frac_pft14 else 0*test_frac_pft14
+let test2_frac_pft15=if(k[k=1:15] eq 15) then test_frac_pft15 else 0*test_frac_pft15
+
+! mask_cont always equals 1
+let mask_cont=mymaxvegetfrac[k=1:19@sum,d=1,gt=myaxis@asn]/mymaxvegetfrac[k=1:19@sum,d=1,gt=myaxis@asn]
+
+
+! This is the vector for k=1:15 which give a description of the PFT distribution of natural lands in each pixel
+! fracnat[k=12] always equals 0, the same for fracnat[k=13]
+let fracnat_init=test2_frac_pft1+test2_frac_pft2+test2_frac_pft3+test2_frac_pft4+test2_frac_pft5+test2_frac_pft6+test2_frac_pft7+test2_frac_pft8+test2_frac_pft9+test2_frac_pft10+test2_frac_pft11+test2_frac_pft14+test2_frac_pft15
+let fracnat=fill_xy(fracnat_init,mask_cont,100)*k[k=1:15]/k[k=1:15]
+
+
+
+! These are just masks
+! 1 for k=13 (or 12), else 0
+let fraccropc4=if(k[k=1:15] eq 13) then mask_cont else 0*mask_cont
+let fraccropc3=if(k[k=1:15] eq 12) then mask_cont else 0*mask_cont
+let fracbare=if(k[k=1:15] eq 1) then mask_cont else 0*mask_cont
+
+! To compute the relative proportion of C3 and C4 grasses in the (C3+C4) grass mix
+! A priori, it's 0 or 1
+let splitgrassPFT10=mymaxvegetfrac[k=10,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])
+let splitgrassPFT11=mymaxvegetfrac[k=11,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])
+let splitgrassPFT14=mymaxvegetfrac[k=14,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])
+let splitgrassPFT15=mymaxvegetfrac[k=15,d=1,gt=myaxis@asn]/(mymaxvegetfrac[k=10:11@sum,d=1,gt=myaxis@asn]+mymaxvegetfrac[k=14:15@sum,d=1,gt=myaxis@asn])
+
+let splitgrassPFT10_fill=fill_xy(splitgrassPFT10,mask_cont,100)
+let splitgrassPFT11_fill=fill_xy(splitgrassPFT11,mask_cont,100)
+let splitgrassPFT14_fill=fill_xy(splitgrassPFT14,mask_cont,100)
+let splitgrassPFT15_fill=fill_xy(splitgrassPFT15,mask_cont,100)
+let splitgrassTOT_fill=splitgrassPFT10_fill+splitgrassPFT11_fill+splitgrassPFT14_fill+splitgrassPFT15_fill
+let fracgrassPFT10=if(k[k=1:15] eq 10) then splitgrassPFT10_fill/splitgrassTOT_fill*mask_cont else 0*mask_cont
+let fracgrassPFT11=if(k[k=1:15] eq 11) then splitgrassPFT11_fill/splitgrassTOT_fill*mask_cont else 0*mask_cont
+let fracgrassPFT14=if(k[k=1:15] eq 14) then splitgrassPFT14_fill/splitgrassTOT_fill*mask_cont else 0*mask_cont
+let fracgrassPFT15=if(k[k=1:15] eq 15) then splitgrassPFT15_fill/splitgrassTOT_fill*mask_cont else 0*mask_cont
+
+
+
+let frac_iceandinlandwater_temp=if(maxvegetfrac[d=1,k=1,gt=myaxis@asn] NE `maxvegetfrac[d=1,k=1,gt=myaxis@asn], RETURN=bad`) then maxvegetfrac[d=1,k=17:18@sum,gt=myaxis@asn]/maxvegetfrac[d=1,k=17:19@sum,gt=myaxis@asn] else 0
+let frac_iceandinlandwater=if(frac_iceandinlandwater_temp NE `frac_iceandinlandwater_temp, RETURN=bad`) then frac_iceandinlandwater_temp else 1
+
+let frac_iceandinlandwater_fill=fill_xy(frac_iceandinlandwater,mask_cont,100)
+let frac_iceandinlandwaterPFT1=if(k[k=1:15] eq 1) then frac_iceandinlandwater_fill*mask_cont else 0*mask_cont
+
+
+
+
+
+
+
+set memory/size=1200
+
+save/file="/home/surface7/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final_future/fractions_from_ESA2.0.7b_19PFT_quarterdeg_PFTmap_2010.nc" fracnat[l=1],fraccropc3[l=1],fraccropc4[l=1],fracgrassPFT10[l=1],fracgrassPFT11[l=1],fracgrassPFT14[l=1],fracgrassPFT15[l=1],frac_iceandinlandwaterPFT1[l=1],fracbare[l=1]
+
+							 
+
Index: /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future_2dstep.jnl
===================================================================
--- /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future_2dstep.jnl	(revision 6267)
+++ /trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/process_landusemap_19PFT_quarterdeg_final_temp_future_2dstep.jnl	(revision 6267)
@@ -0,0 +1,27 @@
+use "/home/surface7/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final_future/fractions_from_ESA2.0.7b_19PFT_quarterdeg_PFTmap_2010.nc"
+
+let file={SPAWN: "find /home/surface7/vuichard/LandUseMap_CMIP6/UofMD-landState-$1/land/yr/multiple-states/gn/ -name multiple-states_input4MIPs_landState_ScenarioMIP_UofMD-$1_gn_2015-2100.nc"}
+
+list file
+
+use "`file`"
+
+
+
+! l=1161 corresponds to year 2010 - The same used for ESA-CCI file
+let tot_luh=primf+primn+secdf+secdn+urban+c3ann+c4ann+c3per+c4per+c3nfx+pastr+range
+
+
+let yearlymaxvegetfrac=(primf[d=2]+primn[d=2]+secdf[d=2]+secdn[d=2]+range[d=2])*fracnat[l=1,d=1]+urban[d=2]*fracbare[l=1,d=1]+(c3ann[d=2]+c3per[d=2]+c3nfx[d=2])*fraccropc3[l=1,d=1]+(c4ann[d=2]+c4per[d=2])*fraccropc4[l=1,d=1]+pastr[d=2]*(fracgrassPFT10[l=1,d=1]+fracgrassPFT11[l=1,d=1]+fracgrassPFT14[l=1,d=1]+fracgrassPFT15[l=1,d=1])+(1.-tot_luh[d=2])*frac_iceandinlandwaterPFT1[l=1,d=1]
+
+
+let yearlymaxvegetfrac_normalize=yearlymaxvegetfrac/yearlymaxvegetfrac[k=@sum]
+
+set memory/size=1200
+
+
+
+repeat/l=1:86 save/file="/home/surface7/vuichard/LandUseMap_CMIP6/UofMD-landState-$1/outputs/$1_`l+2014`.nc" yearlymaxvegetfrac_normalize[l=`l`]
+
+
+
