From 4d2c4edd39071e8a9dc8e0a22fe83edd4163640c Mon Sep 17 00:00:00 2001
From: Bruno Ringeval <bruno.ringeval@inrae.fr>
Date: Tue, 18 Feb 2025 16:39:20 +0100
Subject: [PATCH] progress Recombi output

---
 GPCROP/README.txt                             |  3 ++
 GPCROP/Recombine_gridpack.py                  | 38 ++++++++++---------
 GPCROP/comp_simul_vs_obs_on_site.py           | 18 ---------
 ....py => plot_multi_simul_on_site_wodata.py} |  0
 GPCROP/sbatch_Recombi.sh                      | 13 ++++---
 GPCROP/sbatch_STATS.sh                        | 27 -------------
 6 files changed, 31 insertions(+), 68 deletions(-)
 rename GPCROP/{plot_multi_simul_on_site.py => plot_multi_simul_on_site_wodata.py} (100%)
 delete mode 100755 GPCROP/sbatch_STATS.sh

diff --git a/GPCROP/README.txt b/GPCROP/README.txt
index e754e35..be19c3b 100644
--- a/GPCROP/README.txt
+++ b/GPCROP/README.txt
@@ -85,3 +85,6 @@ scp -r output_main/main_soilmodelv1.1_begin1900_end2018_writebegin2018_writeend2
 #
 #3 simulations: 10111, 11001, 11111
 #./launch_GPCROP_Global.sh
+#
+#4) mkdir good directories for output then sbatch_Recombi.sh
+#scp sbatch_Recombi.sh Recombine_gridpack.py $dist:$scratch_dist/GPCROP/
diff --git a/GPCROP/Recombine_gridpack.py b/GPCROP/Recombine_gridpack.py
index 43b6e88..04a9880 100755
--- a/GPCROP/Recombine_gridpack.py
+++ b/GPCROP/Recombine_gridpack.py
@@ -27,10 +27,11 @@ yri_simul          = int(sys.argv[1]) #only one year in the current version
 ini_nbgrid         = int(sys.argv[2])
 max_nbgrid         = int(sys.argv[3])
 sizepack_nbgrid    = int(sys.argv[4])
+selectvar          = int(sys.argv[5])
 
-config=1
+config=3
 ini_nmonte=0
-max_nmonte=30
+max_nmonte=72
 
 #create lat and lon
 lat_grid_1D = np.zeros((max_nbgrid),np.float32)#np.zeros((max_nbgrid-1),np.float32)
@@ -41,13 +42,13 @@ lon_grid_1D = np.zeros((max_nbgrid),np.float32)#np.zeros((max_nbgrid-1),np.float
 nmonte_len=max_nmonte-ini_nmonte
 time_len=750
 
-name_output="/11111_2009_config1_nmonte0_30/"
-name_simul="11111"
+#name_output="/11111_2009_config1_nmonte0_72/"
+#name_simul="11111"
 
-#name_output="/11001_2009_config1_nmonte0_30/"
-#name_simul="11001"
+name_output="/11001_2009_config1_nmonte0_72/"
+name_simul="11001"
 
-#name_output="/10000_2009_config1_nmonte0_30/"
+#name_output="/10000_2009_config1_nmonte0_72/"
 #name_simul="10000"
 
 ###########################
@@ -64,7 +65,7 @@ while (deb_grid < max_nbgrid):
     print(deb_grid, fin_grid,flush=True)
         
     #open the netcdf file and fill the local variables with variable within netcdf file
-    fich= Dataset(path_output+name_output+"/SIM+GPASOIL_"+name_simul+"_"+str(yri_simul)+"_config"+str(config)+"_nmonte"+str(ini_nmonte)+"_"+str(max_nmonte)+"_grid"+(len(str(max_nbgrid))-len(str(deb_grid)))*'0'+str(deb_grid)+"_"+(len(str(max_nbgrid))-len(str(fin_grid)))*'0'+str(fin_grid)+"_spinup20d.nc",'r',format='NETCDF4')
+    fich= Dataset(path_output+name_output+"/SIM+GPASOIL_"+name_simul+"_"+str(yri_simul)+"_config"+str(config)+"_nmonte"+str(ini_nmonte)+"_"+str(max_nmonte)+"_grid"+(len(str(max_nbgrid))-len(str(deb_grid)))*'0'+str(deb_grid)+"_"+(len(str(max_nbgrid))-len(str(fin_grid)))*'0'+str(fin_grid)+".nc",'r',format='NETCDF4')
 
     if deb_grid==ini_nbgrid:
         list_var = list(fich.variables)
@@ -72,14 +73,14 @@ while (deb_grid < max_nbgrid):
         list_var.remove("lon")
 
         #subsample of var
-        list_var_copy=list_var.copy()
-        for var in list_var_copy:
-            #if 'hour' in var :
-            #    if 'Pilab' not in var:
-            #        list_var.remove(var)
-            #if var[:5] not in ["grain","nleaf","succe"]:
-            if var not in ["grain_data","nleaf_data","nleafwosene_data","success_optim_data",]:
-                list_var.remove(var)
+        if selectvar:
+            list_var_copy=list_var.copy()
+            for var in list_var_copy:
+                #if 'hour' in var :
+                #    if 'Pilab' not in var:
+                #        list_var.remove(var)
+                if var not in ["grain_data","nleaf_data","nleafwosene_data","success_optim_data","RSR","shootwosene_data","root_data","concstar_shootwograin_data","Psupply_day","Pdemand_day"]: #shootwosene_data to compute RSRwosene, Psupply_day and Pdemand_day only useful for no-intera
+                    list_var.remove(var)
         #end subsample of var
         
         for var in list_var:
@@ -114,7 +115,10 @@ while (deb_grid < max_nbgrid):
 ####################
 #create netcdf file#
 ####################
-new_file= Dataset(path_output+name_output+"/SIM+GPASOIL_"+name_simul+"_"+str(yri_simul)+"_config"+str(config)+"_nmonte"+str(ini_nmonte)+"_"+str(max_nmonte)+"_grid_allgrid_spinup20d_selectvar_grainandother.nc",mode='w')
+if selectvar:
+    new_file= Dataset(path_output+name_output+"/SIM+GPASOIL_"+name_simul+"_"+str(yri_simul)+"_config"+str(config)+"_nmonte"+str(ini_nmonte)+"_"+str(max_nmonte)+"_grid_allgrid_selectvar.nc",mode='w')
+else:
+    new_file= Dataset(path_output+name_output+"/SIM+GPASOIL_"+name_simul+"_"+str(yri_simul)+"_config"+str(config)+"_nmonte"+str(ini_nmonte)+"_"+str(max_nmonte)+"_grid_allgrid.nc",mode='w')    
 #define axis of output netcdf file
 grid_dim  = new_file.createDimension('grid', len(lat_grid_1D))
 time_dim  = new_file.createDimension('time', 0)   # unlimited axis
diff --git a/GPCROP/comp_simul_vs_obs_on_site.py b/GPCROP/comp_simul_vs_obs_on_site.py
index 8a732c3..185ad5f 100644
--- a/GPCROP/comp_simul_vs_obs_on_site.py
+++ b/GPCROP/comp_simul_vs_obs_on_site.py
@@ -72,8 +72,6 @@ for var in dic_corres_namedata_namevar.keys(): #["LAI","Nbre"]:
     for treat in ["P0","P1.5","P3"]:
         prefix_var = dic_corres_namedata_namevar[var]+"_obs_"+dic_corres_treatdata_treatvar[treat]
         print(prefix_var)
-        #exec(""+prefix_var+"_avg = np.zeros(len(day_JAS))")
-        #exec(""+prefix_var+"_std = np.zeros(len(day_JAS))")
         exec(""+prefix_var+"     = np.zeros((4,len(day_JAS)),np.float32)") #keep the 4 repeat
         index_d = 0
         for d in day_JAS:
@@ -81,29 +79,13 @@ for var in dic_corres_namedata_namevar.keys(): #["LAI","Nbre"]:
             #TEST: df[df.JAS.isin(["22","152"])][df['Trait']==treat].reset_index()['CPepi'].to_numpy() #select more than 1 value for JAS
             select = ma.masked_where(select!=select,select) #mask
             if len(select.compressed())!=0:
-                #exec(""+prefix_var+"_avg[index_d] = ma.average(select)")
-                #exec(""+prefix_var+"_std[index_d] = ma.std(select)")
                 exec(""+prefix_var+"[:,index_d]   = select")
             else:
-                #exec(""+prefix_var+"_avg[index_d] = -9999.")
-                #exec(""+prefix_var+"_std[index_d] = -9999.")
                 exec(""+prefix_var+"[:,index_d]   = -9999.")
             index_d = index_d + 1
-        #exec(""+prefix_var+"_avg = ma.masked_where("+prefix_var+"_avg[:]==-9999.,"+prefix_var+"_avg[:])")
-        #exec(""+prefix_var+"_std = ma.masked_where("+prefix_var+"_std[:]==-9999.,"+prefix_var+"_std[:])")
         exec(""+prefix_var+"     = ma.masked_where("+prefix_var+"[:]==-9999.,"+prefix_var+"[:])")
 
 
-###create stemleaf variable for obs. not used any more
-#for treat in [dic_corres_treatdata_treatvar[x] for x in dic_corres_treatdata_treatvar.keys()]:
-#    stat="avg"
-#    exec("stemleaf_obs_"+treat+"_"+stat+"  = stem_obs_"+treat+"_"+stat+"  + leaf_obs_"+treat+"_"+stat+"")
-#    exec("stemleafP_obs_"+treat+"_"+stat+" = stemP_obs_"+treat+"_"+stat+" + leafP_obs_"+treat+"_"+stat+"")
-#    #
-#    exec("stemleaf_obs_"+treat+"_std   = ma.sqrt(stem_obs_"+treat+"_std**2  + leaf_obs_"+treat+"_std**2)")
-#    exec("stemleafP_obs_"+treat+"_std  = ma.sqrt(stemP_obs_"+treat+"_std**2  + leafP_obs_"+treat+"_std**2)")
-
-
 ####################
 #read GPCROP output#
 ####################
diff --git a/GPCROP/plot_multi_simul_on_site.py b/GPCROP/plot_multi_simul_on_site_wodata.py
similarity index 100%
rename from GPCROP/plot_multi_simul_on_site.py
rename to GPCROP/plot_multi_simul_on_site_wodata.py
diff --git a/GPCROP/sbatch_Recombi.sh b/GPCROP/sbatch_Recombi.sh
index 0d37cd5..cf3ffd1 100755
--- a/GPCROP/sbatch_Recombi.sh
+++ b/GPCROP/sbatch_Recombi.sh
@@ -7,24 +7,25 @@
 ##SBATCH --account=ispa
 
 #to be launched with
-#sbatch --exclude=muse241,muse098,muse146,muse189,muse108,muse135,muse239,muse099,muse113,muse114,muse115,muse116 --nodes=1 --cpus-per-task=1 -J Recombi -o Recombi.out sbatch_Recombi.sh
+#sbatch --nodes=1 --cpus-per-task=1 -J Recombi -o Recombi_select.out sbatch_Recombi.sh
+#or
+#sbatch --nodes=1 --cpus-per-task=1 -J Recombi -o Recombi.out sbatch_Recombi.sh
 
 path_rel="./"
 source $path_rel/path.txt
 
 yri=2009
 ini_nbgrid=0
-max_nbgrid=14765
+max_nbgrid=14714
 sizepack_nbgrid=200
-
+select_var=0
 
 server=1
 
-
 if [[ $server -eq 0 ]];then
-    singularity exec $path_rel$path_general/BRpython_gp_pulp /opt/miniconda3/bin/python $path_rel/$path_script/Recombine_gridpack.py $yri $ini_nbgrid $max_nbgrid $sizepack_nbgrid
+    singularity exec $path_rel$path_general/BRpython_gp_pulp /opt/miniconda3/bin/python $path_rel/$path_script/Recombine_gridpack.py $yri $ini_nbgrid $max_nbgrid $sizepack_nbgrid $select_var
 else
     ls /dev/loop*
     module load singularity
-    singularity exec $path_rel$path_general/BRpython_gp_pulp /opt/miniconda3/bin/python $path_rel/$path_script/Recombine_gridpack.py $yri $ini_nbgrid $max_nbgrid $sizepack_nbgrid > outRecombi.txt
+    singularity exec $path_rel$path_general/BRpython_gp_pulp /opt/miniconda3/bin/python $path_rel/$path_script/Recombine_gridpack.py $yri $ini_nbgrid $max_nbgrid $sizepack_nbgrid $select_var > outRecombi_"$select_var".txt
 fi
diff --git a/GPCROP/sbatch_STATS.sh b/GPCROP/sbatch_STATS.sh
deleted file mode 100755
index a6f7eeb..0000000
--- a/GPCROP/sbatch_STATS.sh
+++ /dev/null
@@ -1,27 +0,0 @@
-#!/bin/bash
-#SBATCH -N 1 # number of nodes
-#SBATCH -n 1 # number of cores
-#SBATCH --mem=30G
-#SBATCH -t 0-05:00 # time (D-HH:MM)
-##SBATCH --partition=defq
-##SBATCH --account=ispa
-
-#to be launched with
-#sbatch --nodes=1 --cpus-per-task=1 -J STAT -o STAT.out sbatch_STATS.sh
-
-path_rel="./"
-source $path_rel/path.txt
-
-yri=2000
-nmonte_init=799
-nmonte_end=800
-server=0
-
-
-if [[ $server -eq 0 ]];then
-    singularity exec $path_rel$path_general/BRpython_gp /opt/miniconda3/bin/python $path_rel/$path_script/post_treatGPCROP.py $yri $nmonte_init $nmonte_end
-else
-    ls /dev/loop*
-    module load singularity
-    singularity exec $path_rel$path_general/BRpython_gp /opt/miniconda3/bin/python $path_rel/$path_script/post_treatGPCROP.py $yri $nmonte_init $nmonte_end 2>&1
-fi
-- 
GitLab