Changeset 4110 for trunk/LMDZ.COMMON/libf/evolution/pem.F90

Timestamp:

Mar 9, 2026, 10:29:53 AM (5 days ago)

Author:

jbclement

Message:

PEM:

• Introduction of a configurable display/logging system with options 'out2term', 'out2log', 'verbosity_lvl'. All messages now use verbosity levels ('LVL_NFO', 'LVL_WRN', 'LVL_ERR' and 'LVL_DBG').
• Code encapsulation improvements with systematic privacy/protection of module variables.
• Addition of workflow safety checks for required executables, dependencies (e.g. 'ncdump'), input files and callphys keys.
• Renaming of PEM starting and diagnostic files ("startevol.nc" into "startevo.nc", "diagevol.nc" into "diagevo.nc").

JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (14 diffs)

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -26 +26 @@
 use clim_state_init,    only: read_start, read_startfi, read_startpem
 use clim_state_rec,     only: write_restart, write_restartfi, write_restartpem
-use config,             only: read_rundef
-use display,            only: print_ini, print_end, print_msg
+use config,             only: read_rundef, read_display_config
+use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
 use evolution,          only: n_yr_run, n_yr_sim, ntot_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date
 use geometry,           only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface, nlayer
 use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
 use ice_table,          only: icetable_equilibrium, icetable_dynamic, icetable_depth, icetable_thickness, ice_porefilling, evolve_ice_table
-use layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering
+use layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
 use maths,              only: pi
 use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
-use orbit,              only: evol_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
-use output,             only: write_diagevol, dim_ngrid, dim_nsoil
+use orbit,              only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
+use output,             only: write_diagevo, dim_ngrid, dim_nsoil
 use physics,            only: g
 use slopes,             only: subslope_dist, def_slope_mean
@@ -158 +158 @@
 ! ~~~~~~~~~~~~~~
 ! Header
+call read_display_config()
 call print_ini()
 
@@ -203 +204 @@
 call compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb)
 
-! Read the "startevol.nc"
+! Read the "startevo.nc"
 allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,nslope))
 allocate(h2o_ads_reg(ngrid,nsoil,nslope),co2_ads_reg(ngrid,nsoil,nslope),delta_h2o_ads(ngrid),delta_co2_ads(ngrid))
@@ -214 +215 @@
 ! Compute ice tendencies from yearly minima
 allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
-call print_msg('> Computing surface ice tendencies')
+call print_msg('> Computing surface ice tendencies',LVL_NFO)
 call compute_tendice(minPCM_h2operice + minPCM_h2ofrost,h2o_ice(:,:) > 0._dp,d_h2oice)
-call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)))
+call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO)
 call compute_tendice(minPCM_co2perice + minPCM_co2frost,co2_ice(:,:) > 0._dp,d_co2ice)
-call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
+call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO)
 deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost)
 
 ! Save initial set-up useful for the next computations
-call print_msg('> Saving some initial climate state variables')
+call print_msg('> Saving some initial climate state variables',LVL_NFO)
 allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday))
 allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
@@ -250 +251 @@
     end do
 end do
-call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini))
-call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini))
-call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini))
+call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini),LVL_NFO)
+call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini),LVL_NFO)
+call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini),LVL_NFO)
 if (do_sorption) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2_ini,totmass_adsh2o)
 
 ! Main evolution loop
 ! ~~~~~~~~~~~~~~~~~~~
-call print_msg('')
-call print_msg('********* Evolution *********')
+call print_msg('',LVL_NFO)
+call print_msg('********* Evolution *********',LVL_NFO)
 if (do_layering) then
     allocate(h2oice_depth(ngrid,nslope),h2oice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
@@ -284 +285 @@
 idt = 0
 do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim)
-    call print_msg('**** Iteration of the PEM run [Planetary years]: '//real2str(n_yr_run + dt))
+    call print_msg('**** Iteration of the PEM run [Planetary years]: '//real2str(n_yr_run + dt),LVL_NFO)
     ! Evolve global surface pressure
     call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg)
@@ -298 +299 @@
             do i = 1,ngrid
                 call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p)
-                !call print_layering(layerings_map(i,islope))
+                call print_layering(layerings_map(i,islope))
             end do
         end do
@@ -355 +356 @@
 
     ! Output the results
-    call print_msg('> Standard outputs')
-    call write_diagevol('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
+    call print_msg('> Standard outputs',LVL_NFO)
+    call write_diagevo('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
     do islope = 1,nslope
         if (nslope /= 1) then
@@ -363 +364 @@
             num = '_slope'//num
         end if
-        call write_diagevol('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
-        call write_diagevol('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
-        call write_diagevol('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
-        call write_diagevol('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
+        call write_diagevo('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
+        call write_diagevo('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
         if (co2ice_flow) then
-            call write_diagevol('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
+            call write_diagevo('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
             deallocate(is_co2ice_flow)
         end if
         if (h2oice_flow) then
-            call write_diagevol('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
+            call write_diagevo('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
             deallocate(is_h2oice_flow)
         end if
-        call write_diagevol('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
+        call write_diagevo('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
         if (do_soil) then
             if (icetable_equilibrium) then
-                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
-                call write_diagevol('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
+                call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
+                call write_diagevo('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
             else if (icetable_dynamic) then
-                call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
-                call write_diagevol('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
+                call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
+                call write_diagevo('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
             end if
-            call write_diagevol('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
-            call write_diagevol('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
+            call write_diagevo('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
+            call write_diagevo('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
             if (do_sorption) then
-                call write_diagevol('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
-                call write_diagevol('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
+                call write_diagevo('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
+                call write_diagevo('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
             end if
         end if
@@ -404 +405 @@
         end do
         totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
-        call print_msg(" > Relative total CO2 mass balance = "//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %')
+        call print_msg(" > Relative total CO2 mass balance = "//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %',LVL_NFO)
         if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then
-            call print_msg('  /!\ Warning: mass balance is not conserved!')
+            call print_msg('Mass balance is not conserved!',LVL_WRN)
             totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
-            call print_msg('       Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %')
+            call print_msg('    Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %',LVL_WRN)
             totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0
-            call print_msg('       CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %')
+            call print_msg('    CO2 ice mass balance         = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %',LVL_WRN)
             totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
-            call print_msg('       Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %')
+            call print_msg('    Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %',LVL_WRN)
         end if
     end if
@@ -438 +439 @@
 
     ! Check the stopping criteria
-    call print_msg("> Checking the stopping criteria")
+    call print_msg("> Checking the stopping criteria",LVL_NFO)
     call stopping_crit_pressure(ps_avg_glob_ini,ps_avg_glob,stopcrit)
     call stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit)
@@ -448 +449 @@
     if (stopcrit%stop_code() == 0 .and. timewall .and. real(c2 - c1,dp)/real(cr,dp) >= timelimit - antetime) stopcrit%time_limit_reached = .true.
     if (stopcrit%is_any_set()) then
-        call print_msg(stopcrit%stop_message())
+        call print_msg(stopcrit%stop_message(),LVL_NFO)
         exit
     else
-        call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.')
-        call print_msg('**** The workflow has achieved '//real2str(n_yr_sim)//' Planetary years.')
-        call print_msg('**** This run is continuing!')
-        call print_msg('****')
+        call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.',LVL_NFO)
+        call print_msg('**** The workflow has achieved '//real2str(n_yr_sim)//' Planetary years.',LVL_NFO)
+        call print_msg('**** This run is continuing!',LVL_NFO)
+        call print_msg('****',LVL_NFO)
     end if
 end do ! End of the evolution loop
@@ -461 +462 @@
 ! Finalization
 ! ~~~~~~~~~~~~
-call print_msg('')
-call print_msg('********* Finalization *********')
+call print_msg('',LVL_NFO)
+call print_msg('********* Finalization *********',LVL_NFO)
 ! Update orbital parameters
-if (evol_orbit) call update_orbit(n_yr_sim,n_yr_run)
+if (evo_orbit) call update_orbit(n_yr_sim,n_yr_run)
 
 ! Build ice for the PCM
@@ -492 +493 @@
 call build4PCM_surf_rad_prop(h2o_ice,co2_ice,albedo4PCM,emissivity4PCM)
 
-! Write the "startevol.nc"
+! Write the "startevo.nc"
 call write_restartpem(h2o_ice,co2_ice,tsoil_avg,TI,icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map)