Ignore:
Timestamp:
Feb 12, 2026, 9:09:12 AM (3 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 edited

Legend:

Unmodified
Added
Removed

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

    r3991 r4065  
    1414!-----------------------------------------------------------------------
    1515
    16 ! DECLARATION
    17 ! -----------
    18 implicit none
     16! DEPENDENCIES
     17! ------------
     18use numerics, only: dp, di, k4, minieps
     19
     20! DECLARATION
     21! -----------
     22implicit none
     23
     24! PARAMATERS
     25! ----------
     26real(dp), dimension(:,:), allocatable, protected :: tsurf_PCM ! Surface temperature in the PCM at the beginning [K]
    1927
    2028contains
    21 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    22 
    23 !=======================================================================
    24 SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
    25 !-----------------------------------------------------------------------
    26 ! NAME
    27 !     update_tsurf_nearest_baresoil
     29!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     30
     31!=======================================================================
     32SUBROUTINE ini_surf_temp()
     33!-----------------------------------------------------------------------
     34! NAME
     35!     ini_surf_temp
     36!
     37! DESCRIPTION
     38!     Initialize the parameters of module 'surf_temp'.
     39!
     40! AUTHORS & DATE
     41!     JB Clement, 12/2025
     42!
     43! NOTES
     44!
     45!-----------------------------------------------------------------------
     46
     47! DEPENDENCIES
     48! ------------
     49use geometry, only: ngrid, nslope
     50
     51! DECLARATION
     52! -----------
     53implicit none
     54
     55! CODE
     56! ----
     57if (.not. allocated(tsurf_PCM)) allocate(tsurf_PCM(ngrid,nslope))
     58
     59END SUBROUTINE ini_surf_temp
     60!=======================================================================
     61
     62!=======================================================================
     63SUBROUTINE end_surf_temp()
     64!-----------------------------------------------------------------------
     65! NAME
     66!     end_surf_temp
     67!
     68! DESCRIPTION
     69!     Deallocate surf_temp arrays.
     70!
     71! AUTHORS & DATE
     72!     JB Clement, 12/2025
     73!
     74! NOTES
     75!
     76!-----------------------------------------------------------------------
     77
     78! DECLARATION
     79! -----------
     80implicit none
     81
     82! CODE
     83! ----
     84if (allocated(tsurf_PCM)) deallocate(tsurf_PCM)
     85
     86END SUBROUTINE end_surf_temp
     87!=======================================================================
     88
     89!=======================================================================
     90SUBROUTINE set_tsurf_PCM(tsurf_PCM_in)
     91!-----------------------------------------------------------------------
     92! NAME
     93!     set_tsurf_PCM
     94!
     95! DESCRIPTION
     96!     Setter for 'tsurf_PCM'.
     97!
     98! AUTHORS & DATE
     99!     JB Clement, 12/2025
     100!
     101! NOTES
     102!
     103!-----------------------------------------------------------------------
     104
     105! DECLARATION
     106! -----------
     107implicit none
     108
     109! ARGUMENTS
     110! ---------
     111real(dp), dimension(:,:), intent(in) :: tsurf_PCM_in
     112
     113! CODE
     114! ----
     115tsurf_PCM(:,:) = tsurf_PCM_in(:,:)
     116
     117END SUBROUTINE set_tsurf_PCM
     118!=======================================================================
     119
     120!=======================================================================
     121SUBROUTINE adapt_tsurf2disappearedice(surfice,is_ice_ini,ice_disappeared,tsurf_avg)
     122!-----------------------------------------------------------------------
     123! NAME
     124!     adapt_tsurf2disappearedice
    28125!
    29126! DESCRIPTION
    … …  
    40137! DEPENDENCIES
    41138! ------------
    42 use grid_conversion, only: vect2lonlat, lonlat2vect
     139use geometry, only: ngrid, nslope, nlon, nlat, latitudes, vect2lonlat, lonlat2vect
     140use display,  only: print_msg
     141use utility,  only: int2str
    43142
    44143! DECLARATION
    … …  
    48147! ARGUMENTS
    49148! ---------
    50 integer,                          intent(in)    :: nlon, nlat, nslope, ngrid ! Grid dimensions
    51 real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                   ! CO2 ice density
    52 real,    dimension(ngrid),        intent(in)    :: latitude                  ! Latitude
    53 logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_ini             ! Initial CO2 ice flag
    54 real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg                 ! Average surface temperature
    55 logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared        ! Ice disappeared flag
     149real(dp),    dimension(:,:), intent(in)    :: surfice
     150logical(k4), dimension(:,:), intent(in)    :: is_ice_ini
     151real(dp),    dimension(:,:), intent(inout) :: tsurf_avg
     152logical(k4), dimension(:,:), intent(inout) :: ice_disappeared
    56153
    57154! LOCAL VARIABLES
    58155! ---------------
    59 real, parameter                      :: eps = 1.e-10
    60 integer                              :: islope, i, j, k, radius, rmax, di, dj, ii, jj
    61 logical                              :: found
    62 real,    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
    63 real,    dimension(nlon,nlat)        :: latitude_ll
    64 real,    dimension(ngrid)            :: tmp
    65 integer, dimension(nslope - 1)       :: priority
     156integer(di)                              :: islope, i, j, k, radius, rmax, d_i, d_j, ii, jj
     157logical(k4)                              :: found
     158real(dp),    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, mask_co2ice_disappeared
     159real(dp),    dimension(nlon,nlat)        :: latitude_ll
     160real(dp),    dimension(ngrid)            :: tmp
     161integer(di), dimension(nslope - 1)       :: priority
    66162
    67163! CODE
    … …  
    70166if (ngrid == 1) return
    71167
    72 write(*,*) "> Updating surface temperature where ice disappeared"
     168call print_msg("> Adapting surface temperature where ice disappeared")
    73169! Convert from reduced grid to lon-lat grid
    74 call vect2lonlat(nlon,nlat,ngrid,latitude,latitude_ll)
     170call vect2lonlat(nlon,nlat,ngrid,latitudes,latitude_ll)
    75171do islope = 1,nslope
    76172    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
    77     call vect2lonlat(nlon,nlat,ngrid,co2_ice(:,islope),co2ice_ll(:,:,islope))
    78     call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
    79     call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
    80 enddo
     173    call vect2lonlat(nlon,nlat,ngrid,surfice(:,islope),co2ice_ll(:,:,islope))
     174    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,is_ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
     175    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,ice_disappeared(:,islope)),mask_co2ice_disappeared(:,:,islope))
     176end do
    81177
    82178! For each point where ice disappeared
    … …  
    85181    do i = 1,nlon
    86182        do islope = 1,nslope
    87             if (mask_co2ice_ini(i,j,islope) > 0.5 .and. co2ice_ll(i,j,islope) < eps .and. co2ice_disappeared_ll(i,j,islope) < 0.5) then
     183            if (mask_co2ice_ini(i,j,islope) > 0.5_dp .and. co2ice_ll(i,j,islope) < minieps .and. mask_co2ice_disappeared(i,j,islope) < 0.5_dp) then
    88184                found = .false.
    89                 co2ice_disappeared_ll(i,j,islope) = 1.
    90                 call get_slope_priority(latitude_ll(i,j),nslope,islope,priority)
     185                mask_co2ice_disappeared(i,j,islope) = 1._dp
     186                call get_slope_priority(latitude_ll(i,j),islope,priority)
    91187                do k = 1,nslope - 1
    92                     if (mask_co2ice_ini(i,j,priority(k)) < 0.5) then
     188                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5_dp) then
    93189                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
    94190                        found = .true.
    95191                        exit
    96                     endif
    97                 enddo
     192                    end if
     193                end do
    98194
    99195                radius = 1
    100196                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
    101                     do dj = -radius,radius
    102                         do di = -radius,radius
    103                             if (abs(di) + abs(dj) == radius) then
    104                                 ii = i + di
    105                                 jj = j + dj
     197                    do d_j = -radius,radius
     198                        do d_i = -radius,radius
     199                            if (abs(d_i) + abs(d_j) == radius) then
     200                                ii = i + d_i
     201                                jj = j + d_j
    106202                                ! Longitudinal periodicity
    107203                                if (ii < 1) then
    … …  
    109205                                else if (ii > nlon) then
    110206                                    ii = ii - nlon
    111                                 endif
     207                                end if
    112208                                ! Latitude boundaries
    113209                                if (jj >= 1 .and. jj <= nlat) then
    114                                     call get_slope_priority(latitude_ll(ii,jj),nslope,islope,priority)
     210                                    call get_slope_priority(latitude_ll(ii,jj),islope,priority)
    115211                                    do k = 1,nslope - 1
    116                                         if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5) then
     212                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5_dp) then
    117213                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
    118214                                            found = .true.
    119215                                            exit
    120                                         endif
    121                                     enddo
    122                                 endif
    123                             endif
     216                                        end if
     217                                    end do
     218                                end if
     219                            end if
    124220                            if (found) exit
    125                         enddo
     221                        end do
    126222                        if (found) exit
    127                     enddo
     223                    end do
    128224                    radius = radius + 1
    129                 enddo
    130                 if (.not. found) write(*,*) "WARNING: no bare soil found for ice disappeared on  point:",i,j,islope
    131             endif
    132         enddo
    133     enddo
    134 enddo
     225                end do
     226                if (.not. found) call print_msg("Warning: no bare soil found for ice disappeared at point ("//int2str(i)//','//int2str(j)//','//int2str(islope)//'!')
     227            end if
     228        end do
     229    end do
     230end do
    135231
    136232! Convert back from lon-lat grid to reduced grid
    137233do islope = 1,nslope
    138234    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
    139     call lonlat2vect(nlon,nlat,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
    140     where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
    141 enddo
    142 
    143 END SUBROUTINE update_tsurf_nearest_baresoil
    144 !=======================================================================
    145 
    146 !=======================================================================
    147 SUBROUTINE get_slope_priority(lat,nslope,islope,priority)
     235    call lonlat2vect(nlon,nlat,ngrid,mask_co2ice_disappeared(:,:,islope),tmp)
     236    where (tmp > 0.5_dp) ice_disappeared(:,islope) = .true.
     237end do
     238
     239END SUBROUTINE adapt_tsurf2disappearedice
     240!=======================================================================
     241
     242!=======================================================================
     243SUBROUTINE get_slope_priority(lat,islope,priority)
    148244!-----------------------------------------------------------------------
    149245! NAME
    … …  
    151247!
    152248! DESCRIPTION
    153 !     Determine slope priority based on latitude (equator-ward favored).
     249!     Determine slope priority based on latitudes (equator-ward favored).
    154250!
    155251! AUTHORS & DATE
    … …  
    159255!     Equator-ward slopes are most likely to hold no ice.
    160256!-----------------------------------------------------------------------
     257
     258! DEPENDENCIES
     259! ------------
     260use geometry, only: nslope
    161261
    162262! DECLARATION
    … …  
    166266! ARGUMENTS
    167267! ---------
    168 real,                           intent(in)  :: lat      ! Latitude [degrees]
    169 integer,                        intent(in)  :: nslope, islope
    170 integer, dimension(nslope - 1), intent(out) :: priority ! Priority ordering of slopes
     268real(dp),                  intent(in)  :: lat      ! Latitude [degrees]
     269integer(di),               intent(in)  :: islope
     270integer(di), dimension(:), intent(out) :: priority ! Priority ordering of slopes
    171271
    172272! LOCAL VARIABLES
    173273! ---------------
    174 integer :: i, k
     274integer(di) :: i, k
    175275
    176276! CODE
    … …  
    184284        priority(k) = i
    185285        k = k + 1
    186     enddo
     286    end do
    187287    ! Pole-ward slopes
    188288    do i = islope + 1,nslope
    189289        priority(k) = i
    190290        k = k + 1
    191     enddo
     291    end do
    192292else ! Southern hemisphere
    193293    ! Equator-ward slopes
    … …  
    195295        priority(k) = i
    196296        k = k + 1
    197     enddo
     297    end do
    198298    ! Pole-ward slopes
    199299    do i = islope - 1,1,-1
    200300        priority(k) = i
    201301        k = k + 1
    202     enddo
    203 endif
     302    end do
     303end if
    204304
    205305END SUBROUTINE get_slope_priority
    206306!=======================================================================
    207307
     308!=======================================================================
     309SUBROUTINE build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)
     310!-----------------------------------------------------------------------
     311! NAME
     312!     build4PCM_tsurf
     313!
     314! DESCRIPTION
     315!     Reconstructs surface temperature for the PCM.
     316!
     317! AUTHORS & DATE
     318!     JB Clement, 12/2025
     319!
     320! NOTES
     321!
     322!-----------------------------------------------------------------------
     323
     324! DEPENDENCIES
     325! ------------
     326use display, only: print_msg
     327
     328! DECLARATION
     329! -----------
     330implicit none
     331
     332! ARGUMENTS
     333! ---------
     334real(dp), dimension(:,:), intent(in)  :: tsurf_avg, tsurf_dev
     335real(dp), dimension(:,:), intent(out) :: tsurf4PCM
     336
     337! CODE
     338! ----
     339call print_msg('> Building surface temperature for the PCM')
     340tsurf4PCM(:,:) = tsurf_avg(:,:) + tsurf_dev(:,:)
     341
     342END SUBROUTINE build4PCM_tsurf
     343!=======================================================================
     344
    208345END MODULE surf_temp
Note: See TracChangeset for help on using the changeset viewer.