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/geometry.F90

    r4064 r4065  
    1 MODULE grid_conversion
    2 !-----------------------------------------------------------------------
    3 ! NAME
    4 !     grid_conversion
    5 !
    6 ! DESCRIPTION
    7 !     Provides tools to convert data between lon x lat grid format
    8 !     and vector format.
    9 !
    10 ! AUTHORS & DATE
    11 !     JB Clement, 12/2025
    12 !
    13 ! NOTES
    14 !     Handles pole duplication differences between the grid format
    15 !     and vector format.
    16 !-----------------------------------------------------------------------
    17 
    18 ! DECLARATION
    19 ! -----------
    20 implicit none
     1MODULE geometry
     2!-----------------------------------------------------------------------
     3! NAME
     4!     geometry
     5!
     6! DESCRIPTION
     7!     Grid metrics/parameters and basic tools to convert between
     8!     different grids.
     9!
     10! AUTHORS & DATE
     11!     JB Clement, 12/2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DEPENDENCIES
     18! ------------
     19use numerics, only: dp, di, k4
     20
     21! DECLARATION
     22! -----------
     23implicit none
     24
     25! PARAMETERS
     26! ----------
     27integer(di),                         protected :: nlon               ! Number of longitudes
     28integer(di),                         protected :: nlat               ! Number of latitudes
     29integer(di),                         protected :: ngrid              ! Number of grid points
     30integer(di),                         protected :: nlayer             ! Number of atmospheric layers
     31integer(di),                         protected :: nsoil_PCM          ! Number of soil layers in the PCM
     32integer(di),                         parameter :: nsoil = 68_di      ! Number of soil layers in the PEM
     33integer(di),                         protected :: nslope             ! Number of slopes
     34integer(di),                         protected :: nday               ! Number of sols in a year
     35real(dp), dimension(:), allocatable, protected :: longitudes         ! Longitudes
     36real(dp), dimension(:), allocatable, protected :: latitudes          ! Latitudes
     37real(dp), dimension(:), allocatable, protected :: cell_area          ! Cell area
     38real(dp),                            protected :: total_surface      ! Total surface of the grid/planet
     39logical(k4),                             protected :: dim_init = .false. ! Flag to true if dimensions are well initialized
    2140
    2241contains
    23 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     42!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     43
     44!=======================================================================
     45SUBROUTINE ini_geometry()
     46!-----------------------------------------------------------------------
     47! NAME
     48!     ini_geometry
     49!
     50! DESCRIPTION
     51!     Initialize the parameters of module 'geometry'.
     52!
     53! AUTHORS & DATE
     54!     JB Clement, 12/2025
     55!
     56! NOTES
     57!
     58!-----------------------------------------------------------------------
     59
     60! DEPENDENCIES
     61! ------------
     62use stoppage,   only: stop_clean
     63use io_netcdf,  only: startfi_name, xios_day_name1, open_nc, close_nc, get_dim_nc, get_var_nc
     64use display,    only: print_msg
     65use utility,    only: int2str
     66
     67! DECLARATION
     68! -----------
     69implicit none
     70
     71! LOCAL VARIABLES
     72! ---------------
     73logical(k4) :: here, found
     74
     75! CODE
     76! ----
     77! Reading from "startfi.nc" file
     78inquire(file = startfi_name,exist = here)
     79if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//startfi_name//'"!',1)
     80call open_nc(startfi_name,'read')
     81call get_dim_nc('physical_points',ngrid)
     82call get_dim_nc('subsurface_layers',nsoil_PCM)
     83call get_dim_nc('nlayer',nlayer)
     84call get_dim_nc('nslope',nslope,found)
     85if (.not. found) nslope = 1
     86if (.not. allocated(longitudes)) allocate(longitudes(ngrid))
     87if (.not. allocated(latitudes)) allocate(latitudes(ngrid))
     88if (.not. allocated(cell_area)) allocate(cell_area(ngrid))
     89call get_var_nc('longitude',longitudes)
     90call get_var_nc('latitude',latitudes)
     91call get_var_nc('area',cell_area)
     92call close_nc(startfi_name)
     93
     94! Reading from XIOS daily output file
     95inquire(file = xios_day_name1,exist = here)
     96if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//xios_day_name1//'"!',1)
     97call open_nc(xios_day_name1,'read')
     98call get_dim_nc('lon',nlon)
     99call get_dim_nc('lat',nlat)
     100call get_dim_nc('time_counter',nday)
     101call close_nc(xios_day_name1)
     102
     103! Total surface
     104total_surface = sum(cell_area)
     105
     106! State that dimentions are well initialized
     107dim_init = .true.
     108call print_msg('> Reading dimensions')
     109call print_msg('nlon      = '//int2str(nlon))
     110call print_msg('nlat      = '//int2str(nlat))
     111call print_msg('ngrid     = '//int2str(ngrid))
     112call print_msg('nslope    = '//int2str(nslope))
     113call print_msg('nlayer    = '//int2str(nlayer))
     114call print_msg('nsoil_PCM = '//int2str(nsoil_PCM))
     115call print_msg('nsoil     = '//int2str(nsoil))
     116call print_msg('nday      = '//int2str(nday))
     117
     118! Check consistency of grids
     119if (ngrid /= 2 + nlon*(nlat - 2)) call stop_clean(__FILE__,__LINE__,'geometry inconsistency between nlon/nlat and ngrid!',1)
     120if (nsoil_PCM > nsoil) call stop_clean(__FILE__,__LINE__,'geometry inconsistency between nsoil and nsoil_PCM! nsoil must be greater than nsoil_PCM.',1)
     121
     122END SUBROUTINE ini_geometry
     123!=======================================================================
     124
     125!=======================================================================
     126SUBROUTINE end_geometry()
     127!-----------------------------------------------------------------------
     128! NAME
     129!     end_geometry
     130!
     131! DESCRIPTION
     132!     Deallocate geometry arrays.
     133!
     134! AUTHORS & DATE
     135!     JB Clement, 12/2025
     136!
     137! NOTES
     138!
     139!-----------------------------------------------------------------------
     140
     141! DECLARATION
     142! -----------
     143implicit none
     144
     145! CODE
     146! ----
     147if (allocated(longitudes)) deallocate(longitudes)
     148if (allocated(latitudes)) deallocate(latitudes)
     149if (allocated(cell_area)) deallocate(cell_area)
     150
     151END SUBROUTINE end_geometry
     152!=======================================================================
    24153
    25154!=======================================================================
    … …  
    39168!     The longitudes -180 and +180 are not duplicated like in the PCM
    40169!     dynamics.
    41 !
    42 !-----------------------------------------------------------------------
    43 
    44 ! DECLARATION
    45 ! -----------
    46 implicit none
    47 
    48 ! Arguments
     170!-----------------------------------------------------------------------
     171
     172! DEPENDENCIES
     173! ------------
     174use stoppage, only: stop_clean
     175
     176! DECLARATION
     177! -----------
     178implicit none
     179
     180! ARGUMENTS
    49181!----------
    50 integer,                    intent(in)  :: nlon, nlat, ngrid
    51 real, dimension(nlon,nlat), intent(in)  :: v_ll
    52 real, dimension(ngrid),     intent(out) :: v_vect
    53 
    54 ! Local variables
     182integer(di),                    intent(in)  :: nlon, nlat, ngrid
     183real(dp), dimension(nlon,nlat), intent(in)  :: v_ll
     184real(dp), dimension(ngrid),     intent(out) :: v_vect
     185
     186! LOCAL VARIABLES
    55187!----------------
    56 integer :: i, j
    57 
    58 ! Code
     188integer(di) :: i, j
     189
     190! CODE
    59191!-----
    60 ! 1D case
    61 #ifdef CPP_1D
     192if (ngrid == 1) then ! 1D case
    62193    v_vect(1) = v_ll(1,1)
    63194    return
    64 #else
     195else ! 3D
    65196    ! Check
    66     if (ngrid /= nlon*(nlat - 2) + 2) error stop 'lonlat2vect: problem of dimensions!'
     197    if (ngrid /= nlon*(nlat - 2) + 2) call stop_clean(__FILE__,__LINE__,'problem of dimensions!',1)
    67198
    68199    ! Initialization
    … …  
    77208        do i = 1,nlon
    78209            v_vect(1 + i + (j - 2)*nlon) = v_ll(i,j)
    79         enddo
    80     enddo
    81 #endif
     210        end do
     211    end do
     212end if
    82213
    83214END SUBROUTINE lonlat2vect
    … …  
    102233!-----------------------------------------------------------------------
    103234
    104 ! DECLARATION
    105 ! -----------
    106 implicit none
    107 
    108 ! Arguments
     235! DEPENDENCIES
     236! ------------
     237use stoppage, only: stop_clean
     238
     239! DECLARATION
     240! -----------
     241implicit none
     242
     243! ARGUMENTS
    109244!----------
    110 integer,                    intent(in)  :: nlon, nlat, ngrid
    111 real, dimension(ngrid),     intent(in)  :: v_vect
    112 real, dimension(nlon,nlat), intent(out) :: v_ll
    113 
    114 ! Local variables
     245integer(di),                    intent(in)  :: nlon, nlat, ngrid
     246real(dp), dimension(ngrid),     intent(in)  :: v_vect
     247real(dp), dimension(nlon,nlat), intent(out) :: v_ll
     248
     249! LOCAL VARIABLES
    115250!----------------
    116 integer :: i, j
    117 
    118 ! Code
     251integer(di) :: i, j
     252
     253! CODE
    119254!-----
    120 ! 1D case
    121 #ifdef CPP_1D
     255if (ngrid == 1) then ! 1D case
    122256    v_ll(1,1) = v_vect(1)
    123257    return
    124 #else
     258else ! 3D
    125259    ! Check
    126     if (ngrid /= nlon*(nlat - 2) + 2) error stop 'vect2lonlat: problem of dimensions!'
     260    if (ngrid /= nlon*(nlat - 2) + 2) call stop_clean(__FILE__,__LINE__,'problem of dimensions!',1)
    127261
    128262    ! Initialization
    … …  
    137271        do i = 1,nlon
    138272            v_ll(i,j) = v_vect(1 + i + (j - 2)*nlon)
    139         enddo
    140     enddo
    141 #endif
     273        end do
     274    end do
     275end if
    142276
    143277END SUBROUTINE vect2lonlat
    144278!=======================================================================
    145279
    146 END MODULE grid_conversion
     280!=======================================================================
     281SUBROUTINE dyngrd2vect(ni,nj,ngrid,v_dyn,v_vect)
     282!-----------------------------------------------------------------------
     283! NAME
     284!     dyngrd2vect
     285!
     286! DESCRIPTION
     287!     Convert data from dynamical lon x lat grid (where values at the
     288!     poles are duplicated and longitude -180°/+180° is duplicated) into
     289!     a vector.
     290!
     291! AUTHORS & DATE
     292!     JB Clement, 12/2025
     293!
     294! NOTES
     295!     The longitudes -180 and +180 are duplicated (PCM dynamics).
     296!-----------------------------------------------------------------------
     297
     298! DEPENDENCIES
     299! ------------
     300use stoppage, only: stop_clean
     301
     302! DECLARATION
     303! -----------
     304implicit none
     305
     306! ARGUMENTS
     307!----------
     308integer(di),                intent(in)  :: ni, nj, ngrid
     309real(dp), dimension(ni,nj), intent(in)  :: v_dyn
     310real(dp), dimension(ngrid), intent(out) :: v_vect
     311
     312! LOCAL VARIABLES
     313!----------------
     314integer(di) :: i, j
     315
     316! CODE
     317!-----
     318if (ngrid == 1) then ! 1D case
     319    v_vect(1) = v_dyn(1,1)
     320    return
     321else ! 3D
     322    ! Check
     323    if (ngrid /= (ni - 1)*(nj - 2) + 2) call stop_clean(__FILE__,__LINE__,'problem of dimensions!',1)
     324
     325    ! Initialization
     326    v_vect = 0.
     327
     328    ! Treatment of the poles
     329    v_vect(1) = v_dyn(1,1)
     330    v_vect(ngrid) = v_dyn(1,nj)
     331
     332    ! Treatment of regular points
     333    do j = 2,nj - 1
     334        do i = 1,ni - 1
     335            v_vect(1 + i + (j - 2)*(ni - 1)) = v_dyn(i,j)
     336        end do
     337    end do
     338end if
     339
     340END SUBROUTINE dyngrd2vect
     341!=======================================================================
     342
     343!=======================================================================
     344SUBROUTINE vect2dyngrd(ni,nj,ngrid,v_vect,v_dyn)
     345!-----------------------------------------------------------------------
     346! NAME
     347!     vect2dyngrd
     348!
     349! DESCRIPTION
     350!     Convert data from a vector into lon x lat grid (where values at the
     351!     poles are duplicated and longitude -180°/+180° is duplicated).
     352!
     353! AUTHORS & DATE
     354!     JB Clement, 12/2025
     355!
     356! NOTES
     357!     The longitudes -180 and +180 are not duplicated like in the PCM
     358!     dynamics.
     359!-----------------------------------------------------------------------
     360
     361! DEPENDENCIES
     362! ------------
     363use stoppage, only: stop_clean
     364
     365! DECLARATION
     366! -----------
     367implicit none
     368
     369! ARGUMENTS
     370!----------
     371integer(di),                intent(in)  :: ni, nj, ngrid
     372real(dp), dimension(ngrid), intent(in)  :: v_vect
     373real(dp), dimension(ni,nj), intent(out) :: v_dyn
     374
     375! LOCAL VARIABLES
     376!----------------
     377integer(di) :: i, j
     378
     379! CODE
     380!-----
     381if (ngrid == 1) then ! 1D case
     382    v_dyn(1,1) = v_vect(1)
     383    return
     384else ! 3D
     385    ! Check
     386    if (ngrid /= (ni - 1)*(nj - 2) + 2) call stop_clean(__FILE__,__LINE__,'problem of dimensions!',1)
     387
     388    ! Initialization
     389    v_dyn = 0.
     390
     391    ! Treatment of the poles
     392    v_dyn(:,1) = v_vect(1)
     393    v_dyn(:,nj) = v_vect(ngrid)
     394
     395    ! Treatment of regular points
     396    do j = 2,nj - 1
     397        do i = 1,ni - 1
     398            v_dyn(i,j) = v_vect(1 + i + (j - 2)*(ni - 1))
     399        end do
     400        v_dyn(ni,j) =  v_dyn(1,j)
     401    end do
     402end if
     403
     404END SUBROUTINE vect2dyngrd
     405!=======================================================================
     406
     407END MODULE geometry
Note: See TracChangeset for help on using the changeset viewer.