Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3362)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3363)
@@ -348,2 +348,7 @@
 The launching script can now operate a relaunch from any already computed PCM/PEM run. The starting point is asked to the user interactively on the terminal.
 Some features might not work very well yet since not every stituation has not been tested.
+
+== 10/06/2024 == JBC
+- The PEM can now stop itself cleanly before the SLURM time limit for the job is reached and it continues the simulation with a new cycle.
+- Update of "jobPEM.slurm" in the deftank to guarantee enough memory space to run the job.
+- Few minor cleanings.
Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/jobPEM.slurm
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/deftank/jobPEM.slurm	(revision 3362)
+++ /trunk/LMDZ.COMMON/libf/evolution/deftank/jobPEM.slurm	(revision 3363)
@@ -5,5 +5,5 @@
 #SBATCH --constraint=GENOA
 ### Number of Nodes to use
-#SBATCH --nodes=1
+#SBATCH --nodes=4 # to run with enough memory
 #SBATCH --ntasks-per-node=1
 #SBATCH --cpus-per-task=1
@@ -61,3 +61,9 @@
 fi
 
-./launchPEM.sh new
+# Launch the next cycle
+#if [ "$(awk 'END{print $NF}' info_PEM.txt)" -eq 7 ]; then
+#    read i_myear n_myear convert_years iPCM iPEM nPCM nPCM_ini < info_PEM.txt
+#    ./launchPEM.sh cont # Continue the PEM run if it stopped because of job time limit
+#else
+    ./launchPEM.sh new
+#fi
Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh	(revision 3362)
+++ /trunk/LMDZ.COMMON/libf/evolution/deftank/launchPEM.sh	(revision 3363)
@@ -64,4 +64,5 @@
             errlaunch
         fi
+        echo "The relaunch is initialized with a specific previous successful run."
         while true; do
             echo "Do you want to relaunch from a 'PCM' or 'PEM' run?"
@@ -106,4 +107,14 @@
             relaunchPEM
         fi
+
+    # Continuing the PEM run
+    elif [ $1 = "cont" ]; then
+        exec >> log_launchPEM.txt 2>&1
+        echo
+        echo "This is a continuation of the previous PEM run."
+        date
+        submitPEM
+
+    # Default case: error
     else
         echo "Error: given argument '$1' for the launching script is unknown!"
Index: unk/LMDZ.COMMON/libf/evolution/deftank/output_layering.py
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/deftank/output_layering.py	(revision 3362)
+++ 	(revision )
@@ -1,104 +1,0 @@
-####################################################################################
-### Python script to output the stratification data from the "startpem.nc" files ###
-####################################################################################
-
-import netCDF4 as nc
-import numpy as np
-import matplotlib.pyplot as plt
-import sys
-import os.path
-
-##############################
-### Parameters to fill in
-filename = 'startpem9.nc' # File name
-igrid = 1  # Grid point
-islope = 1 # Slope number
-istr = 4   # Stratum number
-##############################
-
-### Open the NetCDF file
-if os.path.isfile(filename):
-    nc_file = nc.Dataset(filename,'r')
-else:
-    sys.exit('The file \"' + filename + '\" does not exist!')
-
-### Get the dimensions
-Time = len(nc_file.dimensions['Time'])
-ngrid = len(nc_file.dimensions['physical_points'])
-nslope = len(nc_file.dimensions['nslope'])
-nb_str_max = len(nc_file.dimensions['nb_str_max'])
-if igrid > ngrid or igrid < 1:
-    sys.exit('Asked grid point is not possible!')
-if islope > nslope or islope < 1:
-    sys.exit('Asked slope number is not possible!')
-if istr > nb_str_max or istr < 1:
-    sys.exit('Asked stratum number is not possible!')
-
-### Get the stratification properties
-stratif_thickness = []
-stratif_top_elevation = []
-stratif_co2ice_volfrac = []
-stratif_h2oice_volfrac = []
-stratif_dust_volfrac = []
-stratif_air_volfrac = []
-for i in range(1,nslope + 1):
-    stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:])
-    stratif_top_elevation.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_top_elevation'][:])
-    stratif_co2ice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_co2ice_volfrac'][:])
-    stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:])
-    stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:])
-    stratif_air_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_air_volfrac'][:])
-
-### Display the data
-igrid = igrid - 1
-islope = islope - 1
-labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air'
-contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1])
-height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1)
-height[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:]
-height[1:] = stratif_top_elevation[islope][0,:,igrid]
-for i in range(len(stratif_top_elevation[islope][0,:,igrid])):
-    contents[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid]
-    contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]
-    contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]
-    contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid]
-contents[:,0] = contents[:,1]
-
-# Simple subplots for a layering
-fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,layout = 'constrained',sharey = True)
-fig.suptitle('Simple content profiles for the layering')
-ax1.step(contents[0,:],height,where = 'post')
-ax2.step(contents[1,:],height,where = 'post')
-ax3.step(contents[2,:],height,where = 'post')
-ax4.step(contents[3,:],height,where = 'post')
-ax1.set_ylabel('Elevation [m]')
-ax1.set_xlabel('Volume fraction [m3/m3]')
-ax2.set_xlabel('Volume fraction [m3/m3]')
-ax3.set_xlabel('Volume fraction [m3/m3]')
-ax4.set_xlabel('Volume fraction [m3/m3]')
-ax1.set_title(labels[0])
-ax2.set_title(labels[1])
-ax3.set_title(labels[2])
-ax4.set_title(labels[3])
-plt.savefig('layering_simpleprofiles.png')
-
-# Stackplot for a layering
-fig, ax = plt.subplots()
-ax.fill_betweenx(height,0,contents[0,:],label = labels[0],color = 'r',step = 'pre')
-ax.fill_betweenx(height,contents[0,:],sum(contents[0:2,:]),label = labels[1],color = 'b',step = 'pre')
-ax.fill_betweenx(height,sum(contents[0:2,:]),sum(contents[0:3,:]),label = labels[2],color = 'y',step = 'pre')
-ax.fill_betweenx(height,sum(contents[0:3,:]),sum(contents),label = labels[3],color = 'g',step = 'pre')
-plt.vlines(x = 0.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
-plt.vlines(x = 1.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
-for i in range(len(height)):
-    plt.hlines(y = height[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--')
-ax.set_title('Stack content profiles for the layering')
-plt.xlabel('Volume fraction [m3/m3]')
-plt.ylabel('Elevation [m]')
-ax.legend(loc = 'center left',bbox_to_anchor = (1,0.5))
-plt.savefig('layering_stackprofiles.png')
-
-plt.show()
-
-### Close the NetCDF file
-nc_file.close()
Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py	(revision 3363)
+++ /trunk/LMDZ.COMMON/libf/evolution/deftank/visu_layering.py	(revision 3363)
@@ -0,0 +1,104 @@
+####################################################################################
+### Python script to output the stratification data from the "startpem.nc" files ###
+####################################################################################
+
+import netCDF4 as nc
+import numpy as np
+import matplotlib.pyplot as plt
+import sys
+import os.path
+
+##############################
+### Parameters to fill in
+filename = 'startpem9.nc' # File name
+igrid = 1  # Grid point
+islope = 1 # Slope number
+istr = 4   # Stratum number
+##############################
+
+### Open the NetCDF file
+if os.path.isfile(filename):
+    nc_file = nc.Dataset(filename,'r')
+else:
+    sys.exit('The file \"' + filename + '\" does not exist!')
+
+### Get the dimensions
+Time = len(nc_file.dimensions['Time'])
+ngrid = len(nc_file.dimensions['physical_points'])
+nslope = len(nc_file.dimensions['nslope'])
+nb_str_max = len(nc_file.dimensions['nb_str_max'])
+if igrid > ngrid or igrid < 1:
+    sys.exit('Asked grid point is not possible!')
+if islope > nslope or islope < 1:
+    sys.exit('Asked slope number is not possible!')
+if istr > nb_str_max or istr < 1:
+    sys.exit('Asked stratum number is not possible!')
+
+### Get the stratification properties
+stratif_thickness = []
+stratif_top_elevation = []
+stratif_co2ice_volfrac = []
+stratif_h2oice_volfrac = []
+stratif_dust_volfrac = []
+stratif_air_volfrac = []
+for i in range(1,nslope + 1):
+    stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:])
+    stratif_top_elevation.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_top_elevation'][:])
+    stratif_co2ice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_co2ice_volfrac'][:])
+    stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:])
+    stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:])
+    stratif_air_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_air_volfrac'][:])
+
+### Display the data
+igrid = igrid - 1
+islope = islope - 1
+labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air'
+contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1])
+height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1)
+height[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:]
+height[1:] = stratif_top_elevation[islope][0,:,igrid]
+for i in range(len(stratif_top_elevation[islope][0,:,igrid])):
+    contents[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid]
+    contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]
+    contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]
+    contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid]
+contents[:,0] = contents[:,1]
+
+# Simple subplots for a layering
+fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,layout = 'constrained',sharey = True)
+fig.suptitle('Simple content profiles for the layering')
+ax1.step(contents[0,:],height,where = 'post')
+ax2.step(contents[1,:],height,where = 'post')
+ax3.step(contents[2,:],height,where = 'post')
+ax4.step(contents[3,:],height,where = 'post')
+ax1.set_ylabel('Elevation [m]')
+ax1.set_xlabel('Volume fraction [m3/m3]')
+ax2.set_xlabel('Volume fraction [m3/m3]')
+ax3.set_xlabel('Volume fraction [m3/m3]')
+ax4.set_xlabel('Volume fraction [m3/m3]')
+ax1.set_title(labels[0])
+ax2.set_title(labels[1])
+ax3.set_title(labels[2])
+ax4.set_title(labels[3])
+plt.savefig('layering_simpleprofiles.png')
+
+# Stackplot for a layering
+fig, ax = plt.subplots()
+ax.fill_betweenx(height,0,contents[0,:],label = labels[0],color = 'r',step = 'pre')
+ax.fill_betweenx(height,contents[0,:],sum(contents[0:2,:]),label = labels[1],color = 'b',step = 'pre')
+ax.fill_betweenx(height,sum(contents[0:2,:]),sum(contents[0:3,:]),label = labels[2],color = 'y',step = 'pre')
+ax.fill_betweenx(height,sum(contents[0:3,:]),sum(contents),label = labels[3],color = 'g',step = 'pre')
+plt.vlines(x = 0.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
+plt.vlines(x = 1.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-')
+for i in range(len(height)):
+    plt.hlines(y = height[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--')
+ax.set_title('Stack content profiles for the layering')
+plt.xlabel('Volume fraction [m3/m3]')
+plt.ylabel('Elevation [m]')
+ax.legend(loc = 'center left',bbox_to_anchor = (1,0.5))
+plt.savefig('layering_stackprofiles.png')
+
+plt.show()
+
+### Close the NetCDF file
+nc_file.close()
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3362)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3363)
@@ -228,10 +228,18 @@
 
 ! Some variables for the PEM run
-real, parameter :: year_step = 1 ! timestep for the pem
-integer         :: year_iter     ! number of iteration
-integer         :: year_iter_max ! maximum number of iterations before stopping
-integer         :: i_myear       ! Global number of Martian years of the chained simulations
-integer         :: n_myear       ! Maximum number of Martian years of the chained simulations
-real            :: timestep      ! timestep [s]
+real, parameter :: year_step = 1   ! Timestep for the pem
+integer         :: year_iter       ! Number of iteration
+integer         :: year_iter_max   ! Maximum number of iterations before stopping
+integer         :: i_myear         ! Global number of Martian years of the chained simulations
+integer         :: n_myear         ! Maximum number of Martian years of the chained simulations
+real            :: timestep        ! Timestep [s]
+character(20)   :: job_id          ! Job id provided as argument passed on the command line when the program was invoked
+integer(kind=8) :: cr              ! Number of clock ticks per second (count rate)
+integer(kind=8) :: c1, c2          ! Counts of processor clock
+character(100)  :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
+real            :: timelimit       ! Time limit for the PEM job in seconds
+real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time time limit: 1200 s = 20 min by default
+integer         :: cstat, days, hours, minutes, seconds
+character(1)    :: sep
 
 #ifdef CPP_STD
@@ -274,4 +282,39 @@
 integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
 
+! Elapsed time with system clock
+call system_clock(count_rate = cr)
+call system_clock(c1)
+timelimit = 86400 ! 86400 seconds = 24 h by default
+if (command_argument_count() > 0) then
+    ! Read the job id passed as argument to the program
+    call get_command_argument(1,job_id)
+    ! Execute the system command
+    call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
+    if (cstat > 0) then
+        error stop 'pem: command execution failed!'
+    else if (cstat < 0) then
+        error stop 'pem: command execution not supported!'
+    endif
+    ! Read the output
+    open(1,file = 'tmp_cmdout.txt',status = 'old')
+    read(1,'(a)') chtimelimit
+    close(1)
+    chtimelimit = trim(chtimelimit)
+    call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
+    if (cstat > 0) then
+        error stop 'pem: command execution failed!'
+    else if (cstat < 0) then
+        error stop 'pem: command execution not supported!'
+    endif
+    ! Assume 'chtimelimit' format is "D-HH:MM:SS" or "HH:MM:SS"
+    if (index(chtimelimit,'-') > 0) then
+        read(chtimelimit, '(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
+        timelimit = days*86400 + hours*3600 + minutes*60 + seconds
+    else
+        read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
+        timelimit = hours*3600 + minutes*60 + seconds
+    endif
+endif
+
 ! Parallel variables
 #ifndef CPP_STD
@@ -310,7 +353,7 @@
 
     therestart1D = .false. ! Default value
-    inquire(file = 'start1D_evol.txt',exist = therestart1D)
+    inquire(file = 'start1D.txt',exist = therestart1D)
     if (.not. therestart1D) then
-        write(*,*) 'There is no "start1D_evol.txt" file!'
+        write(*,*) 'There is no "start1D.txt" file!'
         error stop 'Initialization cannot be done for the 1D PEM.'
     endif
@@ -322,6 +365,6 @@
     endif
 
-    call init_testphys1d('start1D_evol.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, &
-                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,     &
+    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,        &
+                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
                          play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
     ps(2) = ps(1)
@@ -505,5 +548,5 @@
 write(*,*) "Downloading data Y1..."
 call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
-                   tsurf_avg_yr1,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                      &
+                   tsurf_avg_yr1,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                     &
                    co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
 write(*,*) "Downloading data Y1 done!"
@@ -958,4 +1001,6 @@
 !    II_g Checking the stopping criterion
 !------------------------
+    call system_clock(c2)
+    if (timelimit - real((c2 - c1)/cr) <= antetime) stopPEM = 7
     call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
     call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
@@ -981,4 +1026,6 @@
             case(6)
                 write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
+            case(7)
+                write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
             case default
                 write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
@@ -1104,6 +1151,6 @@
     write(*,*) "restart.nc has been written"
 #else
-    call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
-    write(*,*) "restart1D_evol.txt has been written"
+    call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
+    write(*,*) "restart1D.txt has been written"
 #endif
 
Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3362)
+++ /trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3363)
@@ -81,5 +81,5 @@
 call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
 
-! Dowload the data from the file
+! Download the data from the file
 call get_var3("ps",ps_PCM)
 write(*,*) "Data for surface pressure downloaded!"
@@ -172,4 +172,8 @@
 #endif
 endif
+
+! Close the NetCDF file
+write(*,*) "Closing "//filename//"..."
+call error_msg(nf90_close(fID),"close",filename)
 
 ! Compute the minimum over the year for each point
@@ -336,4 +340,5 @@
     case('inq');   msg="Field <"//trim(nam)//"> is missing"
     case('get');   msg="Reading failed for <"//trim(nam)//">"
+    case('put');   msg="Writing failed for <"//trim(nam)//">"
     case('open');  msg="File opening failed for <"//trim(nam)//">"
     case('close'); msg="File closing failed for <"//trim(nam)//">"
