Changeset 4065 for trunk/LMDZ.COMMON/libf/evolution/tendencies.F90

Timestamp:

Feb 12, 2026, 9:09:12 AM (2 weeks ago)

Author:

jbclement

Message:

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

• Some PEM models were corrected or improved:
  • Frost/perennial ice semantics are clarified via renaming;
  • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
  • Geothermal flux for the PCM is computed based on the PEM state;
• New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
• New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
• Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
• Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
  • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
  • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
  • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
• Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
• Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
• Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
• Add module "slopes" to isolate slopes handling;
• Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
• Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
• Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
• The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
• Ice table logic is made restart-safe;
• 'Open'/'read' intrinsic logic is made safe and automatic;
• Improve discoverability and standardize the data handling (private vs protected vs public);
• Apply consistent documentation/header style already introduced;
• Update deftank scripts to reflect new names and launch wrappers;

This revision is a structural milestone aiming to be behavior-preserving where possible. It has been tested via compilation and short integration runs. However, due to extensive renames, moves, and API changes, full validation is still ongoing.
Note: the revision includes one (possibly two) easter egg hidden in the code for future archaeologists of the PEM. No physics were harmed.
JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/tendencies.F90 (modified) (8 diffs)

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

@@ -16 +16 @@
 !-----------------------------------------------------------------------
 
+! DEPENDENCIES
+! ------------
+use numerics, only: dp, di, minieps, tol
+
 ! DECLARATION
 ! -----------
@@ -24 +28 @@
 
 !=======================================================================
-SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice)
+SUBROUTINE compute_tend(min_ice,d_ice)
 !-----------------------------------------------------------------------
 ! NAME
@@ -47 +51 @@
 ! ARGUMENTS
 ! ---------
-integer,                         intent(in)  :: ngrid
-integer,                         intent(in)  :: nslope
-real, dimension(ngrid,nslope,2), intent(in)  :: min_ice ! Minima of ice at each point for the PCM years
-real, dimension(ngrid,nslope),   intent(out) :: d_ice   ! Evolution of perennial ice
+real(dp), dimension(:,:,:), intent(in)  :: min_ice
+real(dp), dimension(:,:),   intent(out) :: d_ice
 
 ! CODE
 ! ----
 ! We compute the difference
-d_ice = min_ice(:,:,2) - min_ice(:,:,1)
+d_ice(:,:) = min_ice(:,:,2) - min_ice(:,:,1)
 
 ! If the difference is too small, then there is no evolution
-where (abs(d_ice) < 1.e-10) d_ice = 0.
+where (abs(d_ice) < minieps) d_ice = 0._dp
 
 ! If the minimum over the last year is 0, then we have no perennial ice
-where (abs(min_ice(:,:,2)) < 1.e-10) d_ice = 0.
+where (abs(min_ice(:,:,2)) < minieps) d_ice = 0._dp
 
 END SUBROUTINE compute_tend
@@ -67 +69 @@
 
 !=======================================================================
-SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
-                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global)
-!-----------------------------------------------------------------------
-! NAME
-!     recomp_tend_co2
+SUBROUTINE evolve_tend_co2(d_co2ice_ini,co2ice,emissivity,q_co2_ts_ini,q_co2_ts,ps_PCM,ps_avg_glob_ini,ps_avg_global,d_co2ice)
+!-----------------------------------------------------------------------
+! NAME
+!     evolve_tend_co2
 !
 ! DESCRIPTION
@@ -86 +87 @@
 ! DEPENDENCIES
 ! ------------
-use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
+use geometry, only: ngrid, nslope, nday
+use planet,   only: alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, m_co2, A, B
+use orbit,    only: yr_len_sol, sol_len_s
+use display,  only: print_msg
+use utility,  only: real2str
 
 ! DECLARATION
@@ -94 +99 @@
 ! ARGUMENTS
 ! ---------
-integer,                        intent(in)    :: timelen, ngrid, nslope ! Time length, # of grid points and slopes
-real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PCM            ! CO2 VMR in PCM first layer
-real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PEM            ! CO2 VMR in PEM first layer
-real, dimension(ngrid,timelen), intent(in)    :: ps_PCM                 ! Surface pressure in PCM
-real,                           intent(in)    :: ps_avg_global_ini      ! Global average pressure (initial)
-real,                           intent(in)    :: ps_avg_global          ! Global average pressure (current)
-real, dimension(ngrid,nslope),  intent(in)    :: d_co2ice_ini           ! Initial CO2 ice evolution
-real, dimension(ngrid,nslope),  intent(in)    :: co2ice                 ! CO2 ice surface [kg/m^2]
-real, dimension(ngrid,nslope),  intent(in)    :: emissivity             ! Emissivity
-real, dimension(ngrid,nslope),  intent(inout) :: d_co2ice_phys          ! Updated CO2 ice evolution
+real(dp), dimension(:,:), intent(in)    :: q_co2_ts, q_co2_ts_ini, ps_PCM, d_co2ice_ini, co2ice, emissivity
+real(dp),                 intent(in)    :: ps_avg_global, ps_avg_glob_ini
+real(dp), dimension(:,:), intent(inout) :: d_co2ice
 
 ! LOCAL VARIABLES
 ! ---------------
-integer :: i, t, islope
-real    :: coef, avg
+integer(di)                     :: i, iday, islope
+real(dp)                        :: coef, avg
+real(dp), dimension(ngrid,nday) :: vmr_co2_ts, vmr_co2_ts_ini
 
 ! CODE
 ! ----
-write(*,*) "> Updating the CO2 ice tendency for the new pressure"
-
-! Evolution of the water ice for each physical point
+call print_msg("> Evolving the CO2 ice tendency according to the new pressure")
+call print_msg('Old CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
+
+! Compute volume mixing ratio
+vmr_co2_ts(:,:) = q_co2_ts(:,:)/(A*q_co2_ts(:,:) + B)/m_co2
+vmr_co2_ts_ini(:,:) = q_co2_ts_ini(:,:)/(A*q_co2_ts_ini(:,:) + B)/m_co2
+
+! Recompute the tendency
 do i = 1,ngrid
     do islope = 1,nslope
-        coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
-        avg = 0.
-        if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
-            do t = 1,timelen
-                avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 &
-                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4
-            enddo
-            if (avg < 1.e-4) avg = 0.
-            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen
-        endif
-    enddo
-enddo
-
-END SUBROUTINE recomp_tend_co2
-!=======================================================================
-
-!=======================================================================
-SUBROUTINE recomp_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice)
-!-----------------------------------------------------------------------
-! NAME
-!     recomp_tend_h2o
+        coef = yr_len_sol*sol_len_s*emissivity(i,islope)*sigmaB/Lco2
+        avg = 0._dp
+        if (co2ice(i,islope) > 1.e-4_dp .and. abs(d_co2ice(i,islope)) > 1.e-5_dp) then
+            do iday = 1,nday
+                avg = avg + (beta_clap_co2/(alpha_clap_co2 - log(q_co2_ts_ini(i,iday)*ps_PCM(i,iday)/100._dp)))**4 &
+                      - (beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_ts(i,iday)*ps_PCM(i,iday)*(ps_avg_global/ps_avg_glob_ini)/100._dp)))**4
+            end do
+            if (avg < 1.e-4_dp) avg = 0._dp
+            d_co2ice(i,islope) = d_co2ice_ini(i,islope) - coef*avg/nday
+        end if
+    end do
+end do
+call print_msg('New CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str( minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)))
+
+END SUBROUTINE evolve_tend_co2
+!=======================================================================
+
+!=======================================================================
+SUBROUTINE evolve_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_ts_old,tsoil_ts_new,d_h2oice)
+!-----------------------------------------------------------------------
+! NAME
+!     evolve_tend_h2o
 !
 ! DESCRIPTION
@@ -160 +165 @@
 ! ARGUMENTS
 ! ---------
-real,                 intent(in)    :: h2oice_depth_old ! Old H2O ice depth
-real,                 intent(in)    :: h2oice_depth_new ! New H2O ice depth
-real,                 intent(in)    :: tsurf            ! Surface temperature
-real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_old ! Old soil temperature time series
-real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_new ! New soil temperature time series
-real,                 intent(inout) :: d_h2oice         ! Evolution of perennial ice
+real(dp),                 intent(in)    :: h2oice_depth_old ! Old H2O ice depth
+real(dp),                 intent(in)    :: h2oice_depth_new ! New H2O ice depth
+real(dp),                 intent(in)    :: tsurf            ! Surface temperature
+real(dp), dimension(:,:), intent(in)    :: tsoil_ts_old     ! Old soil temperature time series
+real(dp), dimension(:,:), intent(in)    :: tsoil_ts_new     ! New soil temperature time series
+real(dp),                 intent(inout) :: d_h2oice         ! Evolution of perennial ice
 
 ! LOCAL VARIABLES
 ! ---------------
-real            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
-integer         :: t
-real, parameter :: coef_diff = 4.e-4 ! Diffusion coefficient
-real, parameter :: zcdv = 0.0325     ! Drag coefficient
+real(dp)            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
+integer(di)         :: iday
+real(dp), parameter :: coef_diff = 4.e-4 ! Diffusion coefficient
+real(dp), parameter :: zcdv = 0.0325     ! Drag coefficient
 
 ! CODE
@@ -182 +187 @@
 
 ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location
-psv_max_old = 0.
-psv_max_new = 0.
-do t = 1,size(tsoil_PEM_timeseries_old,2)
-    psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,h2oice_depth_old)))
-    psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,h2oice_depth_new)))
-enddo
+psv_max_old = 0._dp
+psv_max_new = 0._dp
+do iday = 1,size(tsoil_ts_old,2)
+    psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_ts_old(:,iday),tsurf,h2oice_depth_old)))
+    psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_ts_new(:,iday),tsurf,h2oice_depth_new)))
+end do
 
 ! Lower humidity due to growing lag layer (higher depth)
-if (abs(psv_max_old) < 1.e2*epsilon(1.)) then
-    hum_dec = 1.
+if (abs(psv_max_old) < tol) then
+    hum_dec = 1._dp
 else
     hum_dec = psv_max_new/psv_max_old ! Decrease because of lower water vapor pressure at the new depth
-endif
+end if
 
 ! Flux correction (decrease)
 d_h2oice = d_h2oice*R_dec*hum_dec
 
-END SUBROUTINE recomp_tend_h2o
+END SUBROUTINE evolve_tend_h2o
 !=======================================================================