source: trunk/LMDZ.COMMON/libf/evolution/numerics.F90 @ 4065

Last change on this file since 4065 was 4065, checked in by jbclement, 2 weeks ago

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 size: 6.1 KB
Line 
1MODULE numerics
2!-----------------------------------------------------------------------
3! NAME
4!     numerics
5!
6! DESCRIPTION
7!     Management of numerical precision and special values (NaN and Inf)
8!     throughout the code.
9!
10! AUTHORS & DATE
11!     JB Clement, 01/2026
12!
13! NOTES
14!     The intrinsic function storage_size() returns the storage size of
15!     argument in bits.
16!-----------------------------------------------------------------------
17
18! DEPENDENCIES
19! ------------
20use, intrinsic :: iso_fortran_env, only : int8, int16, int32, int64, real32, real64, real128
21!~ use, intrinsic :: ieee_arithmetic
22
23! DECLARATION
24! -----------
25implicit none
26
27! PARAMETERS
28! ----------
29! Integers
30integer, parameter :: ti = int8                   ! Tiny integer      = 8 bits
31integer, parameter :: si = int16                  ! Short integer     = 16 bits
32integer, parameter :: di = int32                  ! Standard integer  = 32 bits (default)
33integer, parameter :: li = int64                  ! Long integer      = 64 bits
34integer, parameter :: lli = selected_int_kind(30) ! Long long integer = 128 bits
35integer, parameter :: wi = di                     ! Working integer
36
37! Reals
38integer, parameter :: sp = real32  ! Simple precision    = 32  bits = 1-bit sign + 8-bits  exponent + 23-bits  significand
39integer, parameter :: dp = real64  ! Double precision    = 64  bits = 1-bit sign + 11-bits exponent + 52-bits  significand
40integer, parameter :: qp = real128 ! Quadruple precision = 128 bits = 1-bit sign + 15-bits exponent + 112-bits significand
41integer, parameter :: wp = dp      ! Working precision
42
43! Logicals
44! The precision is given by the compiler to be the same as the standard integer precision (32 bits) for memory reasons.
45! It can be interesting to have smaller logical to reduce memory load.
46! But be careful: definition is compiler-dependent and processor-dependent! So the following kinds are not guaranteed for logicals!
47integer, parameter :: k1 = 1 ! kind = 1 should be 8  bits
48integer, parameter :: k2 = 2 ! kind = 2 should be 16 bits
49integer, parameter :: k4 = 4 ! kind = 4 should be 32 bits (default)
50
51! Smallest strictly positive representable real number
52real(sp), parameter :: smallest_sp = tiny(1._sp)
53real(dp), parameter :: smallest_dp = tiny(1._dp)
54real(qp), parameter :: smallest_qp = tiny(1._qp)
55real(wp), parameter :: smallest_nb = tiny(1._wp)
56
57! Largest representable number
58integer(ti),  parameter :: largest_ti = huge(1_ti)
59integer(si),  parameter :: largest_si = huge(1_si)
60integer(di),  parameter :: largest_di = huge(1_di)
61integer(li),  parameter :: largest_li = huge(1_li)
62integer(lli), parameter :: largest_lli = huge(1_lli)
63integer(wi),  parameter :: largest_wi = huge(1_wi)
64real(sp),     parameter :: largest_sp = huge(1._sp)
65real(dp),     parameter :: largest_dp = huge(1._dp)
66real(qp),     parameter :: largest_qp = huge(1._qp)
67real(wp),     parameter :: largest_nb = huge(1._wp)
68
69! Machine epsilon: smallest positive real number such that 1 + epsilon > 1
70real(sp), parameter :: minieps_sp = epsilon(1._sp)
71real(dp), parameter :: minieps_dp = epsilon(1._dp)
72real(qp), parameter :: minieps_qp = epsilon(1._qp)
73real(wp), parameter :: minieps = epsilon(1._wp)
74real(wp), parameter :: tol = 100._wp*minieps
75
76! Minimum number of significant decimal digits
77integer, parameter :: prec_sp = precision(1._sp)
78integer, parameter :: prec_dp = precision(1._dp)
79integer, parameter :: prec_qp = precision(1._qp)
80integer, parameter :: prec = precision(1._wp)
81
82!~ contains
83!~ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84
85!~ !=======================================================================
86!~ FUNCTION is_finite(var) RESULT(flag)
87!~ !-----------------------------------------------------------------------
88!~ ! NAME
89!~ !     is_finite
90!~ !
91!~ ! DESCRIPTION
92!~ !     Returns whether array values are finite.
93!~ !
94!~ ! AUTHORS & DATE
95!~ !     JB Clement, 01/2026
96!~ !
97!~ ! NOTES
98!~ !     Finite values are neither NaN nor Infinite.
99!~ !-----------------------------------------------------------------------
100
101!~ ! DECLARATION
102!~ ! -----------
103!~ implicit none
104
105!~ ! ARGUMENTS
106!~ ! ---------
107!~ real(dp), dimension(..), intent(in)  :: var
108
109!~ ! LOCAL VARIABLES
110!~ ! ---------------
111!~ logical(k4) :: flag
112
113!~ ! CODE
114!~ ! ----
115!~ flag = .false.
116!~ if (all(ieee_is_finite(var))) flag = .true.
117
118!~ END FUNCTION is_finite
119!~ !=======================================================================
120
121!~ !=======================================================================
122!~ FUNCTION is_nan(var) RESULT(flag)
123!~ !-----------------------------------------------------------------------
124!~ ! NAME
125!~ !     is_nan
126!~ !
127!~ ! DESCRIPTION
128!~ !     Returns whether array values has a Not-a-Number (NaN).
129!~ !
130!~ ! AUTHORS & DATE
131!~ !     JB Clement, 01/2026
132!~ !
133!~ ! NOTES
134!~ !
135!~ !-----------------------------------------------------------------------
136
137!~ ! DECLARATION
138!~ ! -----------
139!~ implicit none
140
141!~ ! ARGUMENTS
142!~ ! ---------
143!~ real(dp), dimension(..), intent(in)  :: var
144
145!~ ! LOCAL VARIABLES
146!~ ! ---------------
147!~ logical(k4) :: flag
148
149!~ ! CODE
150!~ ! ----
151!~ flag = .false.
152!~ if (any(ieee_is_nan(var))) flag = .true.
153
154!~ END FUNCTION is_nan
155!~ !=======================================================================
156
157!~ !=======================================================================
158!~ FUNCTION is_normal(var) RESULT(flag)
159!~ !-----------------------------------------------------------------------
160!~ ! NAME
161!~ !     is_normal
162!~ !
163!~ ! DESCRIPTION
164!~ !     Returns whether array values is normal.
165!~ !
166!~ ! AUTHORS & DATE
167!~ !     JB Clement, 01/2026
168!~ !
169!~ ! NOTES
170!~ !     A normal value is finite, non-zero and not subnormal/denormal (very
171!~ !     small finite number with reduced precision).
172!~ !-----------------------------------------------------------------------
173
174!~ ! DECLARATION
175!~ ! -----------
176!~ implicit none
177
178!~ ! ARGUMENTS
179!~ ! ---------
180!~ real(dp), dimension(..), intent(in)  :: var
181
182!~ ! LOCAL VARIABLES
183!~ ! ---------------
184!~ logical(k4) :: flag
185
186!~ ! CODE
187!~ ! ----
188!~ flag = .false.
189!~ if (all(ieee_is_normal(var))) flag = .true.
190
191!~ END FUNCTION is_normal
192!~ !=======================================================================
193
194END MODULE numerics
195
Note: See TracBrowser for help on using the repository browser.