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

    r3991 r4065  
    1414!-----------------------------------------------------------------------
    1515
    16 ! DECLARATION
    17 ! -----------
    18 implicit none
    19 
    20 ! MODULE VARIABLES
    21 ! ----------------
    22 integer                         :: iPCM_qh2o ! Index for H2O vapor tracer from PCM
    23 real, dimension(:), allocatable :: mmol      ! Molar masses of tracers [g/mol]
     16! DEPENDENCIES
     17! ------------
     18use numerics, only: dp, di, k4
     19
     20! DECLARATION
     21! -----------
     22implicit none
     23
     24! PARAMETERS
     25! ----------
     26character(11),                                parameter :: traceurdef_name = 'traceur.def'
     27integer(di),                                  protected :: nq        ! Number of tracers
     28character(30), dimension(:),     allocatable, protected :: qnames    ! Names of tracers
     29real(dp),      dimension(:,:,:), allocatable, protected :: q_PCM     ! Tracers in the "start.nc" at the beginning
     30integer(di),                                  protected :: iPCM_qh2o ! Index for H2O vapor tracer from PCM
     31real(dp),      dimension(:),     allocatable, protected :: mmol      ! Molar masses of tracers [g/mol]
    2432
    2533contains
    … …  
    2735
    2836!=======================================================================
    29 SUBROUTINE ini_tracers_id(nqtot,noms)
    30 !-----------------------------------------------------------------------
    31 ! NAME
    32 !     ini_tracers_id
    33 !
    34 ! DESCRIPTION
    35 !     Initialize tracer indices from PCM tracer names.
    36 !
    37 ! AUTHORS & DATE
    38 !     JB Clement, 12/2025
    39 !
    40 ! NOTES
    41 !
    42 !-----------------------------------------------------------------------
     37SUBROUTINE ini_tracers()
     38!-----------------------------------------------------------------------
     39! NAME
     40!     ini_tracers
     41!
     42! DESCRIPTION
     43!     Allocate tracer molar mass array.
     44!
     45! AUTHORS & DATE
     46!     JB Clement, 12/2025
     47!
     48! NOTES
     49!
     50!-----------------------------------------------------------------------
     51
     52! DEPENDENCIES
     53! ------------
     54use stoppage, only: stop_clean
     55use geometry, only: ngrid, nlayer
     56use display,  only: print_msg
     57use utility,  only: int2str
     58
     59! DECLARATION
     60! -----------
     61implicit none
     62
     63! LOCAL VARIABLES
     64! ---------------
     65logical(k4) :: here
     66integer(di) :: i, ierr, funit
     67
     68! CODE
     69! ----
     70inquire(file = traceurdef_name,exist = here)
     71if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//traceurdef_name//'"!',1)
     72call print_msg('> Reading "'//traceurdef_name//'"')
     73open(newunit = funit,file = traceurdef_name,status = 'old',form = 'formatted',action = 'read',iostat = ierr)
     74if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//traceurdef_name//'"!',ierr)
     75read(funit,*) nq
     76call print_msg('nq = '//int2str(nq))
     77
     78! Allocation
     79if (.not. allocated(mmol)) allocate(mmol(nq))
     80if (.not. allocated(qnames)) allocate(qnames(nq))
     81if (.not. allocated(q_PCM)) allocate(q_PCM(ngrid,nlayer,nq))
     82
     83! Initialize the names of tracers
     84do i = 1,nq
     85    read(funit,*,iostat = ierr) qnames(i)
     86    if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error reading the names of tracers in "'//traceurdef_name//'"!',1)
     87end do
     88close(funit)
     89
     90! Initialize the indices of tracers
     91iPCM_qh2o = -1
     92do i = 1,nq
     93    if (qnames(i) == "h2o_vap") then
     94        iPCM_qh2o = i
     95        mmol(i) = 18._dp
     96    end if
     97end do
     98
     99! Checking if everything has been found
     100if (iPCM_qh2o < 0) call stop_clean(__FILE__,__LINE__,'H2O vapour index not found!',1)
     101
     102END SUBROUTINE ini_tracers
     103!=======================================================================
     104
     105!=======================================================================
     106SUBROUTINE end_tracers()
     107!-----------------------------------------------------------------------
     108! NAME
     109!     end_tracers
     110!
     111! DESCRIPTION
     112!     Deallocate tracer molar mass array.
     113!
     114! AUTHORS & DATE
     115!     JB Clement, 12/2025
     116!
     117! NOTES
     118!
     119!-----------------------------------------------------------------------
     120
     121! DECLARATION
     122! -----------
     123implicit none
     124
     125! CODE
     126! ----
     127if (allocated(mmol)) deallocate(mmol)
     128if (allocated(qnames)) deallocate(qnames)
     129if (allocated(q_PCM)) deallocate(q_PCM)
     130
     131END SUBROUTINE end_tracers
     132!=======================================================================
     133
     134!=======================================================================
     135SUBROUTINE set_q_PCM(q_PCM_in,iq)
     136!-----------------------------------------------------------------------
     137! NAME
     138!     set_q_PCM
     139!
     140! DESCRIPTION
     141!     Setter for 'q_PCM'.
     142!
     143! AUTHORS & DATE
     144!     JB Clement, 12/2025
     145!
     146! NOTES
     147!
     148!-----------------------------------------------------------------------
     149
     150! DEPENDENCIES
     151! ------------
     152use geometry, only: dyngrd2vect, nlon, nlat, ngrid
    43153
    44154! DECLARATION
    … …  
    48158! ARGUMENTS
    49159! ---------
    50 integer,                        intent(in) :: nqtot ! Total number of tracers
    51 character(*), dimension(nqtot), intent(in) :: noms  ! Names of tracers
     160real(dp), dimension(:,:,:), intent(in) :: q_PCM_in
     161integer(di),                intent(in) :: iq
    52162
    53163! LOCAL VARIABLES
    54164! ---------------
    55 integer :: i
    56 
    57 ! CODE
    58 ! ----
    59 ! Allocation
    60 call ini_tracers(nqtot)
    61 
    62 ! Initialization
    63 iPCM_qh2o = -1
    64 
    65 ! Getting the index
    66 do i = 1,nqtot
    67     if (noms(i) == "h2o_vap") then
    68         iPCM_qh2o = i
    69         mmol(i) = 18.
    70     endif
    71 enddo
    72 
    73 ! Checking if everything has been found
    74 if (iPCM_qh2o < 0) error stop 'ini_frost_id: H2O vapour index not found!'
    75 
    76 END SUBROUTINE ini_tracers_id
    77 !=======================================================================
    78 
    79 !=======================================================================
    80 SUBROUTINE ini_tracers(nqtot)
    81 !-----------------------------------------------------------------------
    82 ! NAME
    83 !     ini_tracers
    84 !
    85 ! DESCRIPTION
    86 !     Allocate tracer molar mass array.
    87 !
    88 ! AUTHORS & DATE
    89 !     JB Clement, 12/2025
    90 !
    91 ! NOTES
    92 !
    93 !-----------------------------------------------------------------------
     165integer(di) :: l
     166
     167! CODE
     168! ----
     169do l = 1,size(q_PCM_in,3)
     170    call dyngrd2vect(nlon + 1,nlat,ngrid,q_PCM_in(:,:,l),q_PCM(:,l,iq))
     171end do
     172
     173END SUBROUTINE set_q_PCM
     174!=======================================================================
     175
     176!=======================================================================
     177SUBROUTINE adapt_tracers2pressure(ps_avg_glob_old,ps_avg_glob,ps_ts,q_h2o_ts,q_co2_ts)
     178!-----------------------------------------------------------------------
     179! NAME
     180!     adapt_tracers2pressure
     181!
     182! DESCRIPTION
     183!     Adapt the timeseries of tracers VMR to the new pressure.
     184!
     185! AUTHORS & DATE
     186!     JB Clement, 01/2026
     187!
     188! NOTES
     189!
     190!-----------------------------------------------------------------------
     191
     192! DEPENDENCIES
     193! ------------
     194use geometry,   only: ngrid, nday
     195use atmosphere, only: ap, bp
     196use display,    only: print_msg
    94197
    95198! DECLARATION
    … …  
    99202! ARGUMENTS
    100203! ---------
    101 integer, intent(in) :: nqtot ! Total number of tracers
    102 
    103 ! CODE
    104 ! ----
    105 if (.not. allocated(mmol)) allocate(mmol(nqtot))
    106 
    107 END SUBROUTINE ini_tracers
    108 !=======================================================================
    109 
    110 !=======================================================================
    111 SUBROUTINE end_tracers()
    112 !-----------------------------------------------------------------------
    113 ! NAME
    114 !     end_tracers
    115 !
    116 ! DESCRIPTION
    117 !     Deallocate tracer molar mass array.
    118 !
    119 ! AUTHORS & DATE
    120 !     JB Clement, 12/2025
    121 !
    122 ! NOTES
    123 !
    124 !-----------------------------------------------------------------------
    125 
    126 ! DECLARATION
    127 ! -----------
    128 implicit none
    129 
    130 ! CODE
    131 ! ----
    132 if (allocated(mmol)) deallocate(mmol)
    133 
    134 END SUBROUTINE end_tracers
     204real(dp),                        intent(in)    :: ps_avg_glob_old, ps_avg_glob
     205real(dp), dimension(ngrid,nday), intent(inout) :: ps_ts, q_h2o_ts, q_co2_ts
     206
     207! LOCAL VARIABLES
     208! ---------------
     209real(dp), dimension(ngrid,nday) :: plev1, plev2, plev1_old, plev2_old
     210
     211! CODE
     212! ----
     213call print_msg("> Adapting the timeseries of tracers VMR to the new pressure")
     214
     215! Build the first levels of pressure
     216plev1_old(:,:) = ap(1) + bp(1)*ps_ts(:,:)
     217plev2_old(:,:) = ap(2) + bp(2)*ps_ts(:,:)
     218
     219! Evolve timeseries of surface pressure
     220ps_ts(:,:) = ps_ts(:,:)*ps_avg_glob/ps_avg_glob_old
     221
     222! Build the first levels of pressure
     223plev1(:,:) = ap(1) + bp(1)*ps_ts(:,:)
     224plev2(:,:) = ap(2) + bp(2)*ps_ts(:,:)
     225
     226! Adapt the timeseries of VMR tracers
     227! H2O
     228q_h2o_ts(:,:) = q_h2o_ts(:,:)*(plev1_old(:,:) - plev2_old(:,:))/(plev1(:,:) - plev2(:,:))
     229where (q_h2o_ts < 0._dp) q_h2o_ts = 0._dp
     230where (q_h2o_ts > 1._dp) q_h2o_ts = 1._dp
     231! CO2
     232q_co2_ts(:,:) = q_co2_ts(:,:)*(plev1_old(:,:) - plev2_old(:,:))/(plev1(:,:) - plev2(:,:))
     233where (q_co2_ts < 0._dp) q_co2_ts = 0._dp
     234where (q_co2_ts > 1._dp) q_co2_ts = 1._dp
     235
     236END SUBROUTINE adapt_tracers2pressure
     237!=======================================================================
     238
     239!=======================================================================
     240SUBROUTINE build4PCM_tracers(ps4PCM,q4PCM)
     241!-----------------------------------------------------------------------
     242! NAME
     243!     build4PCM_tracers
     244!
     245! DESCRIPTION
     246!     Reconstructs the tracer VMRs for the PCM.
     247!
     248! AUTHORS & DATE
     249!     JB Clement, 12/2025
     250!     C. Metz, 01/2026
     251!
     252! NOTES
     253!
     254!-----------------------------------------------------------------------
     255
     256! DEPENDENCIES
     257! ------------
     258use atmosphere, only: ap, bp, ps_PCM
     259use geometry,   only: ngrid, nlayer
     260use display,    only: print_msg
     261
     262! DECLARATION
     263! -----------
     264implicit none
     265
     266! ARGUMENTS
     267! ---------
     268real(dp), dimension(:),     intent(in)  :: ps4PCM
     269real(dp), dimension(:,:,:), intent(out) :: q4PCM
     270
     271! LOCAL VARIABLES
     272! ---------------
     273integer(di)                           :: i, l, iq
     274real(dp)                              :: extra_mass
     275real(dp), dimension(:,:), allocatable :: dz_plev_PCM, dz_plev_new
     276
     277! CODE
     278! ----
     279call print_msg('> Building tracers for the PCM')
     280allocate(dz_plev_PCM(ngrid,nlayer),dz_plev_new(ngrid,nlayer))
     281
     282! Compute thickness of pressure levels (= atmospheric layers)
     283do l = 1,nlayer
     284    dz_plev_PCM(:,l) = ap(l) + bp(l)*ps_PCM(:) - ap(l + 1) - bp(l + 1)*ps_PCM(:)
     285    dz_plev_new(:,l) = ap(l) + bp(l)*ps_PCM(:) - ap(l + 1) - bp(l + 1)*ps4PCM(:)
     286end do
     287
     288! Reconstruct tracers
     289do iq = 1,nq
     290    do l = 1,nlayer
     291        do i = 1,ngrid
     292            ! Standard scaling
     293            q4PCM(i,l,iq) = q_PCM(i,l,iq)*dz_plev_PCM(i,l)/dz_plev_new(i,l)
     294
     295            ! CO2 logic: account for the dry air mass change
     296            if (trim(qnames(iq)) == "co2") q4PCM(i,l,iq) = q4PCM(i,l,iq) + (dz_plev_new(i,l) - dz_plev_PCM(i,l)) / dz_plev_new(i,l)
     297        end do
     298    end do
     299end do
     300
     301! Mass conservation and clipping
     302do iq = 1,nq
     303    do i = 1,ngrid
     304        do l = 1,nlayer
     305            ! Clipping negative values
     306            if (q4PCM(i,l,iq) < 0._dp) q4PCM(i,l,iq) = 0._dp
     307
     308            ! VMR limit check (ignore number density tracers)
     309            if (q4PCM(i,l,iq) > 1._dp .and.              &
     310                (qnames(iq) /= "dust_number") .and.      &
     311                (qnames(iq) /= "ccn_number") .and.       &
     312                (qnames(iq) /= "stormdust_number") .and. &
     313                (qnames(iq) /= "topdust_number")) then
     314
     315                if (l < nlayer) then
     316                    extra_mass = (q4PCM(i,l,iq) - 1.)*dz_plev_new(i,l) ! extra_mass units: [VMR*pressure]
     317                    q4PCM(i,l,iq) = 1._dp
     318
     319                    ! Transfer extra mass to the layer above
     320                    ! Note: We divide by dz_plev_new(i,l + 1) to convert back to VMR
     321                    q4PCM(i,l + 1,iq) = q4PCM(i,l + 1,iq) + extra_mass/dz_plev_new(i,l + 1)
     322                else
     323                    q4PCM(i,l,iq) = 1._dp
     324                end if
     325            end if
     326        end do
     327    end do
     328end do
     329
     330deallocate(dz_plev_PCM,dz_plev_new)
     331
     332END SUBROUTINE build4PCM_tracers
    135333!=======================================================================
    136334
Note: See TracChangeset for help on using the changeset viewer.