From f5d472003bd7672c12abb73fe4db6a8a1b2d7353 Mon Sep 17 00:00:00 2001
From: Bruno Ringeval <bruno.ringeval@inrae.fr>
Date: Wed, 12 Mar 2025 17:38:41 +0100
Subject: [PATCH] correct issue in senescence

---
 GPCROP/GPCROP.py                                  | 6 +++---
 GPCROP/README.txt                                 | 2 +-
 SIM_vsMueller_2023/modules/compute_consistency.py | 7 -------
 SIM_vsMueller_2023/modules/compute_senescence.py  | 4 +++-
 4 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/GPCROP/GPCROP.py b/GPCROP/GPCROP.py
index 054ba73..1551c7a 100755
--- a/GPCROP/GPCROP.py
+++ b/GPCROP/GPCROP.py
@@ -748,7 +748,7 @@ for yr in np.arange(yri_simul,yri_simul+1): #pour l instant qu'une seule annee.
         #################
         #end uncertainty#
         #################
-
+        
         #write text file to keep track of the param
         if flag_grid and prescribe_gridcond:
             #no uncertainty about soil P in that case
@@ -1040,7 +1040,7 @@ for yr in np.arange(yri_simul,yri_simul+1): #pour l instant qu'une seule annee.
                         success_optim=0
                         problem.writeMPS(path_rel+"/MPS_optim_nosuccess.mps")
                         problem.writeLP(path_rel+"/MPS_optim_nosuccess.lp")
-                        print("no optim for:", grid_select[grid], flush=True)
+                        print("no optim for:", day, grid_select[grid], flush=True)
                         for param in list_param:
                             exec("print(param,"+param+"[grid],flush=True)")
                         for param in list_paramGPCROP:
@@ -1165,7 +1165,7 @@ for yr in np.arange(yri_simul,yri_simul+1): #pour l instant qu'une seule annee.
             path_rel=path_scriptSIM
             script=path_rel+"/path.txt";exec(open(script).read())
             script=path_rel+"/modules/compute_senescence.py";exec(open(script).read()) #
-
+            
             #remobilization: add corresponding P to reserve P pool. without remobilization, P lost due to senescence is not added to the reserve pool
             if yes_remob:
                 reserveP_data[day,:] = reserveP_data[day,:] + 1e+4 * 1e-3 * (shoot_data_saved[:] - shoot_data[day,:])  * concstar_shootwograin_data[day,:] #all change in biom is here attributed to senescence of shootwograin. conc of shootwograin is constar
diff --git a/GPCROP/README.txt b/GPCROP/README.txt
index 45792f0..b1b0da1 100644
--- a/GPCROP/README.txt
+++ b/GPCROP/README.txt
@@ -39,7 +39,7 @@ mkdir -p GPASOIL/output_main/
 cd GPCROP/
 scp -r GPCROP.py RSR_conc.py param_RSR_conc.py Psupplydemand_precomputation_optim.py Psupplydemand_core_optim.py Psupplydemand_postcomputation_optim.py  $dist:$scratch_dist/GPCROP/.
 scp -r launch_GPCROP_Tartas.sh sbatch_GPCROP_init.sh path.txt  $dist:$scratch_dist/GPCROP/.
-scp -r GPCROPinputfile_for_GPCROPuncertainty_Tartasprescribed_LAILIN.txt GPCROPinputfile_for_GPCROPuncertainty_global.txt $dist:$scratch_dist/GPCROP/.
+scp -r GPCROPinputfile_for_GPCROPuncertainty_Tartasprescribed_LAILIN_nmonteline.txt GPCROPinputfile_for_GPCROPuncertainty_global.txt $dist:$scratch_dist/GPCROP/.
 
 #
 cd GPASOIL/
diff --git a/SIM_vsMueller_2023/modules/compute_consistency.py b/SIM_vsMueller_2023/modules/compute_consistency.py
index 8330895..778fa25 100644
--- a/SIM_vsMueller_2023/modules/compute_consistency.py
+++ b/SIM_vsMueller_2023/modules/compute_consistency.py
@@ -17,10 +17,3 @@ elif LAIrelation=="LIN":
     nleaf_data[day,:] = ma.where(success_optim_data[day,:] == 0, LAI_data[day,:]/SD, nleaf_data[day,:])
 
 nleafwosene_data[day,:] = ma.where(success_optim_data[day,:] == 0, ma.amin( ma.array([maxnleaf*ma.ones(nleafwosene_data[day-1,:].shape, np.float32), (nleafwosene_data[day-1,:] + (nleaf_data[day,:] - nleaf_data[day-1,:]))]) ,axis=0), nleafwosene_data[day,:])
-
-#if day==121:
-#    stop
-#if day>=100:
-#    print(day," after consis:", shoot_data[day,0],nleaf_data[day,0],LAI_data[day,0])
-#if day>=130:
-#    stop
diff --git a/SIM_vsMueller_2023/modules/compute_senescence.py b/SIM_vsMueller_2023/modules/compute_senescence.py
index 2826f3d..fd2a426 100644
--- a/SIM_vsMueller_2023/modules/compute_senescence.py
+++ b/SIM_vsMueller_2023/modules/compute_senescence.py
@@ -40,6 +40,7 @@ if LAIrelation=="EXP":
     exec("LAI_"+model+"[day,:]  = ma.where(nleaf_data_keep[:]!=nleaf_data[day,:], aLAI * (ma.exp( bLAI * nleaf_"+model+"[day,:]) - 1), LAI_"+model+"[day,:])")
 elif LAIrelation=="LIN":
     exec("LAI_"+model+"[day,:]  = ma.where(nleaf_data_keep[:]!=nleaf_data[day,:], SD * nleaf_"+model+"[day,:], LAI_"+model+"[day,:])")
+    
 
 #reduce shoot according to senescence
 shoot_data_saved = shoot_data[day,:].copy()
@@ -52,6 +53,7 @@ else: #solu 2: decrease based on SLA
 exec("shootwograin_"+model+"[day,:] = shoot_"+model+"[day,:] - grain_"+model+"[day,:]")
 
 #senescence_data=1 if senescence starts!
-senescence_data[:] = ma.where((senescence_data[:]==0).astype(np.int32) + (nleafwosene_data[day,:]>nleaf_data[day,:]).astype(np.int32) == 2, 1, senescence_data[:]) 
+#senescence_data[:] = ma.where((senescence_data[:]==0).astype(np.int32) + (nleafwosene_data[day,:]>nleaf_data[day,:]).astype(np.int32) == 2, 1, senescence_data[:]) 
+senescence_data[:] = ma.where((senescence_data[:]==0).astype(np.int32) + (nleafwosene_data[day,:]>nleaf_data[day,:]).astype(np.int32) + (ma.fabs(nleafwosene_data[day,:]-nleaf_data[day,:])>1e-4).astype(np.int32) == 3, 1, senescence_data[:]) # put a condition with fabs to prevent error due to diff in precision between nleafwosene_data and nleaf_data
 senescenceday_data[day,:] = senescence_data[:] #daily variable
 
-- 
GitLab