Ignore:
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:
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/maths.F90

    r4064 r4065  
    1 module math_toolkit
    2 !-----------------------------------------------------------------------
    3 ! NAME
    4 !     math_toolkit
     1module maths
     2!-----------------------------------------------------------------------
     3! NAME
     4!     maths
    55!
    66! DESCRIPTION
    … …  
    1515!-----------------------------------------------------------------------
    1616
    17 ! DECLARATION
    18 ! -----------
    19 implicit none
     17! DEPENDENCIES
     18! ------------
     19use numerics, only: dp, di, k4, minieps
     20
     21! DECLARATION
     22! -----------
     23implicit none
     24
     25! PARAMETERS
     26! ----------
     27real(dp),    parameter :: pi = 4._dp*atan(1._dp) ! PI = 3.14159...
     28integer(di), parameter :: limiter_minmod  = 1
     29integer(di), parameter :: limiter_MC      = 2
     30integer(di), parameter :: limiter_vanLeer = 3
    2031
    2132contains
    … …  
    4758! ARGUMENTS
    4859! ---------
    49 integer,             intent(in)  :: nz       ! number of layer
    50 real, dimension(nz), intent(in)  :: z        ! depth layer
    51 real, dimension(nz), intent(in)  :: y        ! function which needs to be derived
    52 real,                intent(in)  :: y0, ybot ! boundary conditions
    53 real, dimension(nz), intent(out) :: dzY      ! derivative of y w.r.t depth
    54 
    55 ! LOCAL VARIABLES
    56 ! ---------------
    57 integer :: j
    58 real    :: hm, hp, c1, c2, c3
     60integer(di),             intent(in)  :: nz       ! number of layer
     61real(dp), dimension(nz), intent(in)  :: z        ! depth layer
     62real(dp), dimension(nz), intent(in)  :: y        ! function which needs to be derived
     63real(dp),                intent(in)  :: y0, ybot ! boundary conditions
     64real(dp), dimension(nz), intent(out) :: dzY      ! derivative of y w.r.t depth
     65
     66! LOCAL VARIABLES
     67! ---------------
     68integer(di) :: j
     69real(dp)    :: hm, hp, c1, c2, c3
    5970
    6071! CODE
    … …  
    6980    hm = z(j) - z(j - 1)
    7081    c1 = +hm/(hp*(z(j + 1) - z(j - 1)))
    71     c2 = 1/hm - 1/hp
     82    c2 = 1._dp/hm - 1._dp/hp
    7283    c3 = -hp/(hm*(z(j + 1) - z(j - 1)))
    7384    dzY(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
    74 enddo
    75 dzY(nz) = (ybot - y(nz - 1))/(2.*(z(nz) - z(nz - 1)))
     85end do
     86dzY(nz) = (ybot - y(nz - 1))/(2._dp*(z(nz) - z(nz - 1)))
    7687
    7788END SUBROUTINE deriv1
    … …  
    101112! ARGUMENTS
    102113! ---------
    103 integer,             intent(in)  :: nz
    104 real, dimension(nz), intent(in)  :: z, y
    105 real,                intent(in)  :: y0, yNp1
    106 real, dimension(nz), intent(out) :: yp2
    107 
    108 ! LOCAL VARIABLES
    109 ! ---------------
    110 integer :: j
    111 real    :: hm, hp, c1, c2, c3
    112 
    113 ! CODE
    114 ! ----
    115 c1 = +2./((z(2) - z(1))*z(2))
    116 c2 = -2./((z(2) - z(1))*z(1))
    117 c3 = +2./(z(1)*z(2))
     114integer(di),             intent(in)  :: nz
     115real(dp), dimension(nz), intent(in)  :: z, y
     116real(dp),                intent(in)  :: y0, yNp1
     117real(dp), dimension(nz), intent(out) :: yp2
     118
     119! LOCAL VARIABLES
     120! ---------------
     121integer(di) :: j
     122real(dp)    :: hm, hp, c1, c2, c3
     123
     124! CODE
     125! ----
     126c1 = +2._dp/((z(2) - z(1))*z(2))
     127c2 = -2._dp/((z(2) - z(1))*z(1))
     128c3 = +2._dp/(z(1)*z(2))
    118129yp2(1) = c1*y(2) + c2*y(1) + c3*y0
    119130do j = 2,nz - 1
    120131    hp = z(j + 1) - z(j)
    121132    hm = z(j) - z(j - 1)
    122     c1 = +2./(hp*(z(j + 1) - z(j - 1)))
    123     c2 = -2./(hp*hm)
    124     c3 = +2./(hm*(z(j + 1) - z(j-1)))
     133    c1 = +2._dp/(hp*(z(j + 1) - z(j - 1)))
     134    c2 = -2._dp/(hp*hm)
     135    c3 = +2._dp/(hm*(z(j + 1) - z(j - 1)))
    125136    yp2(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
    126 enddo
    127 yp2(nz) = (yNp1 - 2*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2
     137end do
     138yp2(nz) = (yNp1 - 2._dp*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2
    128139
    129140END SUBROUTINE deriv2_simple
    … …  
    154165! ARGUMENTS
    155166! ---------
    156 integer,             intent(in)  :: nz, j
    157 real, dimension(nz), intent(in)  :: z, y
    158 real,                intent(out) :: dy_zj
    159 
    160 ! LOCAL VARIABLES
    161 ! ---------------
    162 real :: h1, h2, c1, c2, c3
     167integer(di),             intent(in)  :: nz, j
     168real(dp), dimension(nz), intent(in)  :: z, y
     169real(dp),                intent(out) :: dy_zj
     170
     171! LOCAL VARIABLES
     172! ---------------
     173real(dp) :: h1, h2, c1, c2, c3
    163174
    164175! CODE
    165176! ----
    166177if (j < 1 .or. j > nz - 2) then
    167     dy_zj = -1.
     178    dy_zj = -1._dp
    168179else
    169180    h1 = z(j + 1) - z(j)
    170181    h2 = z(j + 2)- z(j + 1)
    171     c1 = -(2*h1 + h2)/(h1*(h1 + h2))
     182    c1 = -(2._dp*h1 + h2)/(h1*(h1 + h2))
    172183    c2 = (h1 + h2)/(h1*h2)
    173184    c3 = -h1/(h2*(h1 + h2))
    174185    dy_zj = c1*y(j) + c2*y(j + 1) + c3*y(j + 2)
    175 endif
     186end if
    176187
    177188END SUBROUTINE deriv1_onesided
    … …  
    201212! ARGUMENTS
    202213! ---------
    203 integer,             intent(in)  :: nz, i1, i2
    204 real, dimension(nz), intent(in)  :: y, z
    205 real,                intent(out) :: integral
    206 
    207 ! LOCAL VARIABLES
    208 ! ---------------
    209 integer             :: i
    210 real, dimension(nz) :: dz
     214integer(di),             intent(in)  :: nz, i1, i2
     215real(dp), dimension(nz), intent(in)  :: y, z
     216real(dp),                intent(out) :: integral
     217
     218! LOCAL VARIABLES
     219! ---------------
     220integer(di)             :: i
     221real(dp), dimension(nz) :: dz
    211222
    212223! CODE
    … …  
    215226do i = 2,nz - 1
    216227    dz(i) = (z(i + 1) - z(i - 1))/2.
    217 enddo
     228end do
    218229dz(nz) = z(nz) - z(nz - 1)
    219230integral = sum(y(i1:i2)*dz(i1:i2))
    … …  
    245256! ARGUMENTS
    246257! ---------
    247 real, intent(in)  :: y1, y2
    248 real, intent(in)  :: z1, z2
    249 real, intent(out) :: zr
     258real(dp), intent(in)  :: y1, y2
     259real(dp), intent(in)  :: z1, z2
     260real(dp), intent(out) :: zr
    250261
    251262! CODE
    … …  
    283294! ARGUMENTS
    284295! ---------
    285 integer,                intent(in)  :: n
    286 real, dimension(n),     intent(in)  :: b, d
    287 real, dimension(n - 1), intent(in)  :: a, c
    288 real, dimension(n),     intent(out) :: x
    289 integer,                intent(out) :: error
    290 
    291 ! LOCAL VARIABLES
    292 ! ---------------
    293 integer            :: i
    294 real               :: m
    295 real, dimension(n) :: cp, dp
     296integer(di),                intent(in)  :: n
     297real(dp), dimension(n),     intent(in)  :: b, d
     298real(dp), dimension(n - 1), intent(in)  :: a, c
     299real(dp), dimension(n),     intent(out) :: x
     300integer(di),                intent(out) :: error
     301
     302! LOCAL VARIABLES
     303! ---------------
     304integer(di)            :: i
     305real(dp)               :: m
     306real(dp), dimension(n) :: cp, dp
    296307
    297308! CODE
    … …  
    302313    error = 1
    303314    return
    304 endif
     315end if
    305316do i = 2,n - 1
    306317    if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then
    307318        error = 1
    308319        return
    309     endif
    310 enddo
     320    end if
     321end do
    311322if (abs(b(n)) < abs(a(n - 1))) then
    312323    error = 1
    313324    return
    314 endif
     325end if
    315326
    316327! Initialization
    … …  
    323334    cp(i) = c(i)/m
    324335    dp(i) = (d(i) - a(i - 1)*dp(i - 1))/m
    325 enddo
     336end do
    326337m = b(n) - a(n - 1)*cp(n - 1)
    327338dp(n) = (d(n) - a(n - 1)*dp(n - 1))/m
    … …  
    331342do i = n - 1,1,-1
    332343    x(i) = dp(i) - cp(i)*x(i + 1)
    333 enddo
     344end do
    334345
    335346END SUBROUTINE solve_tridiag
    … …  
    358369! DEPENDENCIES
    359370! ------------
    360 use stop_pem, only: stop_clean
    361 
    362 ! DECLARATION
    363 ! -----------
    364 implicit none
    365 
    366 ! ARGUMENTS
    367 ! ---------
    368 integer,                intent(in)  :: n
    369 real, dimension(n),     intent(in)  :: z, kappa
    370 real, dimension(n - 1), intent(in)  :: mz, mkappa
    371 real,                   intent(in)  :: T_left, q_right
    372 real, dimension(n),     intent(out) :: T
    373 
    374 ! LOCAL VARIABLES
    375 ! ---------------
    376 integer                :: i, error
    377 real, dimension(n)     :: b, d
    378 real, dimension(n - 1) :: a, c
     371use stoppage, only: stop_clean
     372
     373! DECLARATION
     374! -----------
     375implicit none
     376
     377! ARGUMENTS
     378! ---------
     379integer(di),                intent(in)  :: n
     380real(dp), dimension(n),     intent(in)  :: z, kappa
     381real(dp), dimension(n - 1), intent(in)  :: mz, mkappa
     382real(dp),                   intent(in)  :: T_left, q_right
     383real(dp), dimension(n),     intent(out) :: T
     384
     385! LOCAL VARIABLES
     386! ---------------
     387integer(di)                :: i, error
     388real(dp), dimension(n)     :: b, d
     389real(dp), dimension(n - 1) :: a, c
    379390
    380391! CODE
    381392! ----
    382393! Initialization
    383 a = 0.; b = 0.; c = 0.; d = 0.
     394a = 0._dp; b = 0._dp; c = 0._dp; d = 0._dp
    384395
    385396! Left boundary condition (Dirichlet: prescribed temperature)
    386 b(1) = 1.
     397b(1) = 1._dp
    387398d(1) = T_left
    388399
    … …  
    392403    c(i) = -mkappa(i)/((mz(i) - mz(i - 1))*(z(i + 1) - z(i)))
    393404    b(i) = -(a(i - 1) + c(i))
    394 enddo
     405end do
    395406
    396407! Right boundary condition (Neumann: prescribed temperature)
    … …  
    401412! Solve the tridiagonal linear system with the Thomas' algorithm
    402413call solve_tridiag(a,b,c,d,n,T,error)
    403 if (error /= 0) call stop_clean("solve_steady_heat","Unstable solving!",1)
     414if (error /= 0) call stop_clean(__FILE__,__LINE__,"unstable solving!",1)
    404415
    405416END SUBROUTINE solve_steady_heat
    406417!=======================================================================
    407418
    408 end module math_toolkit
     419end module maths
Note: See TracChangeset for help on using the changeset viewer.