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

Timestamp:

Apr 10, 2026, 7:17:55 PM (6 hours ago)

Author:

jbclement

Message:

PEM:

• Rework layering-related logic, especially clarify interactions between surface and subsurface water tendencies and disable CO2 lag resistance (inconsistent without updating pressure and mass balance + PCM).
• Prevent simultaneous activation of layering and ice flows.
• Add warning when flux_geo /= 0 while soil is disabled.
• Add new utility function "is_lvl_enabled" for displaying.
• Replace deprecated 'minieps' with 'eps'/'tol'.

JBC

File:

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

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

@@ -31 +31 @@
 use clim_state_init,    only: read_start, read_startfi, read_startevo
 use config,             only: read_rundef, read_display_config
-use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
+use display,            only: print_ini, print_end, print_msg, is_lvl_enabled, LVL_NFO, LVL_WRN, LVL_DBG
 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, nsoil_PCM, nsoil, cell_area, total_surface
@@ -38 +38 @@
 use layered_deposits,   only: 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 numerics,           only: dp, qp, di, li, k4, eps, eps_qp
 use orbit,              only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
 use output,             only: write_diagevo, dim_ngrid, dim_nsoil
@@ -233 +233 @@
     if (do_layering) then
         h2oice_depth_old(:,:) = h2oice_depth(:,:)
-        ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency
-        where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:)
         do islope = 1,nslope
             do i = 1,ngrid
-                call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(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 evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice(i,islope),h2oice_depth(i,islope),flux_ssice_avg(i,islope), &
+                                     new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p,i,islope)
+                if (is_lvl_enabled(LVL_DBG)) call print_layering(layerings_map(i,islope))
             end do
         end do
@@ -255 +254 @@
         call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit)
         call evolve_co2ice(co2_ice,d_co2ice,zshift_surf)
-    end if
-
-    ! Flow glaciers according to surface ice
-    if (co2ice_flow) then
-        allocate(is_co2ice_flow(ngrid,nslope))
-        call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
-    end if
-    if (h2oice_flow) then
-        allocate(is_h2oice_flow(ngrid,nslope))
-        call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
-    end if
-    if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map)
+
+        ! Flow glaciers according to surface ice
+        if (co2ice_flow) then
+            allocate(is_co2ice_flow(ngrid,nslope))
+            call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow)
+        end if
+        if (h2oice_flow) then
+            allocate(is_h2oice_flow(ngrid,nslope))
+            call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow)
+        end if
+    end if ! do_layering
 
     ! Adapt surface temperature if surface ice disappeared
@@ -337 +335 @@
 
     ! Checking mass balance for CO2
-    if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then
+    if (abs(CO2cond_ps_PCM - 1._dp) < eps) then
         totmass_co2ice = 0._dp
         totmass_atmco2 = 0._dp
@@ -346 +344 @@
             end do
         end do
-        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0
+        totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,eps_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)//' %',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('Mass balance is not conserved!',LVL_WRN)
-            totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0
+            totmass_ini = max(totmass_atmco2_ini,eps_qp) ! To avoid division by 0
             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
+            totmass_ini = max(totmass_co2ice_ini,eps_qp) ! To avoid division by 0
             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
+            totmass_ini = max(totmass_adsco2_ini,eps_qp) ! To avoid division by 0
             call print_msg('    Adsorbed CO2 mass balance    = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %',LVL_WRN)
         end if
@@ -412 +410 @@
 call update_workflow_status(n_yr_run,stopcrit%stop_code(),n_yr_sim,ntot_yr_sim)
 
-! Deallocation
-if (do_layering) then
-    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
-    do islope = 1,nslope
-        do i = 1,ngrid
-            call del_layering(layerings_map(i,islope))
-        end do
-    end do
-end if
+! Clean-up
+if (do_layering) deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
 call end_allocation()