Ignore:
Timestamp:
Mar 26, 2026, 4:29:10 PM (2 weeks ago)
Author:
jbclement
Message:

PEM:

  • Separate variables ownership between the module "planet" for persistent climate state and the program "pem" for transient workflow logic. It provides a meaningful structure.
  • Add lifecycle helpers for clear allocation/deallocation logic.
  • Simplify string suffix for slopes variables.

JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4152 r4157  
    88!
    99! AUTHORS & DATE
    10 !     R. Vandemeulebrouck, 22/07/2022
     10!     R. Vandemeulebrouck, 22/07/2022 with r2778 & r2779
    1111!     L. Lange, 22/07/2022
    1212!     JB Clement, 2023-2025
    1313!
    1414! NOTES
    15 !
     15!     Ownership criterion for declarations:
     16!       - Declare in module "planet" variables that are persistent climate
     17!         state or persistent references used across multiple time steps.
     18!       - Declare in program "pem.F90" transient workflow/control varaibles,
     19!         one-shot initialization buffers and per-iteration working buffers.
    1620!-----------------------------------------------------------------------
    1721
    … …  
    2529use atmosphere,         only: ps_PCM, evolve_pressure, CO2cond_ps_PCM
    2630use backup,             only: save_clim_state, backup_rate
    27 use clim_state_init,    only: read_start, read_startfi, read_startpem
     31use clim_state_init,    only: read_start, read_startfi, read_startevo
    2832use config,             only: read_rundef, read_display_config
    2933use display,            only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN
    … …  
    3236use glaciers,           only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers
    3337use ice_table,          only: icetable_equilibrium, icetable_dynamic, evolve_ice_table
    34 use layered_deposits,   only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
     38use layered_deposits,   only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering
    3539use maths,              only: pi
    3640use numerics,           only: dp, qp, di, li, k4, minieps, minieps_qp
    3741use orbit,              only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit
    3842use output,             only: write_diagevo, dim_ngrid, dim_nsoil
     43use planet
    3944use physics,            only: g
    4045use slopes,             only: subslope_dist, def_slope_mean
    … …  
    6065! ---------------
    6166! Utility-related:
    62 integer(li)               :: cr     ! Number of clock ticks per second (count rate)
    63 integer(li)               :: c1, c2 ! Counts of processor clock
    64 character(:), allocatable :: num    ! To write slope variables
    65 integer(di)               :: i, islope
    66 ! Pressure-related:
    67 real(dp), dimension(:),   allocatable :: ps_avg          ! Average surface pressure [Pa]
    68 real(dp), dimension(:),   allocatable :: ps_dev          ! Deviation of surface pressure [Pa]
    69 real(dp), dimension(:,:), allocatable :: ps_ts           ! Surface pressure timeseries [Pa]
    70 real(dp)                              :: ps_avg_glob_ini ! Global average pressure at initialization [Pa]
    71 real(dp)                              :: ps_avg_glob_old ! Global average pressure of previous time step [Pa]
    72 real(dp)                              :: ps_avg_glob     ! Global average pressure of current time step [Pa]
     67integer(li)  :: cr     ! Number of clock ticks per second (count rate)
     68integer(li)  :: c1, c2 ! Counts of processor clock
     69character(8) :: num    ! Slope suffix to ouput variables
     70integer(di)  :: i, islope
    7371! Ice-related:
    74 real(dp),    dimension(:,:),   allocatable :: h2o_ice                    ! H2O ice [kg.m-2]
    75 real(dp),    dimension(:,:),   allocatable :: co2_ice                    ! CO2 ice [kg.m-2]
    76 real(dp)                                   :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2]
    77 real(dp)                                   :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2]
    78 logical(k4), dimension(:,:),   allocatable :: is_h2oice_ini              ! Initial location of H2O ice
    79 logical(k4), dimension(:,:),   allocatable :: is_co2ice_ini              ! Initial location of CO2 ice
    80 logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow             ! Flag for location of CO2 glacier flow
    81 logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow             ! Flag for location of H2O glacier flow
    82 real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice              ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
    83 real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice              ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
    84 real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost               ! Minimum of H2O frost over the last PCM year [kg.m-2]
    85 real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost               ! Minimum of CO2 frost over the last PCM year [kg.m-2]
    86 logical(k4), dimension(:,:),   allocatable :: is_co2ice_disappeared      ! Flag to check if CO2 ice disappeared at the previous timestep
     72logical(k4), dimension(:,:),   allocatable :: is_co2ice_flow   ! Flag for location of CO2 glacier flow
     73logical(k4), dimension(:,:),   allocatable :: is_h2oice_flow   ! Flag for location of H2O glacier flow
     74real(dp),    dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg.m-2]
     75real(dp),    dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2]
     76real(dp),    dimension(:,:,:), allocatable :: minPCM_h2ofrost  ! Minimum of H2O frost over the last PCM year [kg.m-2]
     77real(dp),    dimension(:,:,:), allocatable :: minPCM_co2frost  ! Minimum of CO2 frost over the last PCM year [kg.m-2]
    8778! Surface-related:
    88 real(dp), dimension(:,:), allocatable :: tsurf_avg           ! Average surface temperature [K]
    89 real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1       ! Average surface temperature of the second to last PCM run [K]
    90 real(dp), dimension(:,:), allocatable :: tsurf_dev           ! Deviation of surface temperature [K]
    91 real(dp), dimension(:,:), allocatable :: h2o_surfdensity_avg ! Average water surface density [kg/m^3]
    92 real(dp), dimension(:,:), allocatable :: zshift_surf         ! Elevation shift for the surface [m]
    93 real(dp), dimension(:,:), allocatable :: zlag                ! Newly built lag thickness [m]
    94 ! Soil-related:
    95 real(dp), dimension(:,:,:),   allocatable :: tsoil_avg                ! Average soil temperature [K]
    96 real(dp), dimension(:,:,:),   allocatable :: tsoil_dev                ! Deviation pf soil temperature [K]
    97 real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts                 ! Soil temperature timeseries [K]
    98 real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts_old             ! Soil temperature timeseries at the previous time step [K]
    99 real(dp), dimension(:,:,:),   allocatable :: h2o_soildensity_avg      ! Average of soil water soil density [kg/m^3]
    100 real(dp), dimension(:),       allocatable :: delta_co2_ads            ! Quantity of CO2 exchanged due to adsorption/desorption [kg/m^2]
    101 real(dp), dimension(:),       allocatable :: delta_h2o_ads            ! Quantity of H2O exchanged due to adsorption/desorption [kg/m^2]
    102 real(qp)                                  :: totmass_adsh2o           ! Total mass of H2O exchanged because of adsorption/desorption [kg]
    103 real(dp), dimension(:,:,:),   allocatable :: h2o_ads_reg              ! H2O adsorbed in the regolith [kg/m^2]
    104 real(dp), dimension(:,:,:),   allocatable :: co2_ads_reg              ! CO2 adsorbed in the regolith [kg/m^2]
    105 real(dp), dimension(:,:),     allocatable :: icetable_depth           ! Depth of the ice table [m]
    106 real(dp), dimension(:,:),     allocatable :: icetable_thickness       ! Thickness of the ice table [m]
    107 real(dp), dimension(:,:,:),   allocatable :: ice_porefilling          ! Amount of porefilling [m^3/m^3]
    108 real(dp), dimension(:,:),     allocatable :: icetable_depth_old       ! Old depth of the ice table [m]
    109 real(dp), dimension(:),       allocatable :: delta_icetable ! Total mass of the H2O that has sublimated / condenses from the ice table [kg]
    110 ! Tracer-related:
    111 real(dp), dimension(:,:),   allocatable :: q_co2_ts     ! CO2 mass mixing ratio in the first layer [kg/kg]
    112 real(dp), dimension(:,:),   allocatable :: q_co2_ts_ini ! Initial CO2 mass mixing ratio in the first layer [kg/kg]
    113 real(dp), dimension(:,:),   allocatable :: q_h2o_ts     ! H2O mass mixing ratio in the first layer [kg/kg]
     79real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K]
     80real(dp), dimension(:,:), allocatable :: zshift_surf   ! Elevation shift for the surface [m]
     81real(dp), dimension(:,:), allocatable :: zlag          ! Newly built lag thickness [m]
    11482! Tendency-related:
    115 real(dp), dimension(:,:), allocatable :: d_co2ice     ! Tendency of perennial CO2 ice [kg/m2/y]
    116 real(dp), dimension(:,:), allocatable :: d_co2ice_ini ! Tendency of perennial CO2 ice at the beginning [kg/m2/y]
    117 real(dp), dimension(:,:), allocatable :: d_h2oice     ! Tendency of perennial H2O ice [kg/m2/y]
    11883real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y]
    11984! Layering-related:
    120 type(layering), dimension(:,:), allocatable :: layerings_map    ! Layering for each grid point and slope
    12185type(ptrarray), dimension(:,:), allocatable :: current          ! Current active stratum in the layering
    12286real(dp),       dimension(:,:), allocatable :: h2oice_depth     ! Depth of subsurface ice layer
    … …  
    13094real(qp) :: totmass_co2ice_ini, totmass_atmco2_ini, totmass_adsco2_ini ! Initial total CO2 masses (surface ice|atmospheric|adsorbed)
    13195real(qp) :: totmass_ini                                                ! Initial total CO2 mass [kg]
     96real(qp) :: totmass_adsh2o                                             ! Current total adsorbed H2O mass [kg]
    13297real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm   ! Variables to balance H2O ice reservoirs
    13398
    … …  
    168133
    169134! Read the PCM data given by XIOS
    170 allocate(ps_avg(ngrid),ps_ts(ngrid,nday))
    171 allocate(tsurf_avg(ngrid,nslope),tsurf_avg_yr1(ngrid,nslope),h2o_surfdensity_avg(ngrid,nslope))
    172 allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope))
    173 allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday))
    174 allocate(minPCM_h2operice(ngrid,nslope,2),minPCM_co2perice(ngrid,nslope,2),minPCM_h2ofrost(ngrid,nslope,2),minPCM_co2frost(ngrid,nslope,2))
     135call allocate_xios_state()
     136allocate(tsurf_avg_yr1(ngrid,nslope))
     137allocate(minPCM_h2operice(ngrid,nslope,2))
     138allocate(minPCM_co2perice(ngrid,nslope,2))
     139allocate(minPCM_h2ofrost(ngrid,nslope,2))
     140allocate(minPCM_co2frost(ngrid,nslope,2))
    175141
    176142call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, &
    … …  
    181147
    182148! Compute the deviation from the average
    183 allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoil_PCM,nslope))
     149call allocate_deviation_state()
    184150ps_dev(:) = ps_PCM(:) - ps_avg(:)
    185151tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:)
    … …  
    193159
    194160! Read the "startevo.nc"
    195 allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,nslope))
    196 allocate(icetable_depth(ngrid,nslope),icetable_thickness(ngrid,nslope),ice_porefilling(ngrid,nsoil,nslope))
    197 allocate(h2o_ads_reg(ngrid,nsoil,nslope),co2_ads_reg(ngrid,nsoil,nslope),delta_h2o_ads(ngrid),delta_co2_ads(ngrid))
    198 allocate(layerings_map(ngrid,nslope))
    199 icetable_depth(:,:) = 0._dp
    200 icetable_thickness(:,:) = 0._dp
    201 ice_porefilling(:,:,:) = 0._dp
    202 delta_h2o_ads(:) = 0._dp
    203 delta_co2_ads(:) = 0._dp
    204 call read_startpem(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,h2o_ice,co2_ice, &
     161call allocate_startevo_state()
     162call read_startevo(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,h2o_ice,co2_ice, &
    205163                   tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map,       &
    206164                   h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
    … …  
    208166
    209167! Compute ice tendencies from yearly minima
    210 allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope))
     168call allocate_tendencies()
    211169call print_msg('> Computing surface ice tendencies',LVL_NFO)
    212170call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice)
    … …  
    218176! Save initial set-up useful for the next computations
    219177call print_msg('> Saving some initial climate state variables',LVL_NFO)
    220 allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday))
    221 allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
     178call allocate_initial_snapshots()
    222179ps_avg_glob_ini = ps_avg_glob
    223180d_co2ice_ini(:,:) = d_co2ice(:,:)
    … …  
    225182h2oice_sublim_coverage_ini = 0._dp
    226183co2ice_sublim_coverage_ini = 0._dp
    227 is_h2oice_ini(:,:) = .false.
    228 is_co2ice_ini(:,:) = .false.
    229184totmass_co2ice_ini = 0._qp
    230185totmass_atmco2_ini = 0._qp
    … …  
    268223    !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot
    269224end if
    270 if (nslope == 1) then ! No slope
    271     allocate(character(0) :: num)
    272 else ! Using slopes
    273     allocate(character(8) :: num)
    274 end if
    275 allocate(delta_icetable(ngrid),icetable_depth_old(ngrid,nslope),is_co2ice_disappeared(ngrid,nslope),tsoil_ts_old(ngrid,nsoil,nslope,nday))
    276 is_co2ice_disappeared(:,:) = .false.
    277 delta_icetable(:) = 0._dp
     225call allocate_loop_state()
    278226n_yr_run = 0
    279227idt = 0
    … …  
    353301    call write_diagevo('ps_avg_glob','Global average pressure','Pa',ps_avg_glob)
    354302    do islope = 1,nslope
    355         if (nslope /= 1) then
    356             num = '  '
    357             write(num,'(i2.2)') islope
    358             num = '_slope'//num
    359         end if
    360         call write_diagevo('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
    361         call write_diagevo('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
    362         call write_diagevo('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
    363         call write_diagevo('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
     303        if (nslope == 1) then
     304            num = ''
     305        else
     306            write(num,'("_slope",i2.2)') islope
     307        end if
     308        call write_diagevo('h2oice'//trim(num),'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/))
     309        call write_diagevo('co2ice'//trim(num),'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/))
     310        call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/))
     311        call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/))
    364312        if (co2ice_flow) then
    365             call write_diagevo('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
    366             deallocate(is_co2ice_flow)
     313            call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/))
    367314        end if
    368315        if (h2oice_flow) then
    369             call write_diagevo('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
    370             deallocate(is_h2oice_flow)
    371         end if
    372         call write_diagevo('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
     316            call write_diagevo('Flow_h2oice'//trim(num),'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/))
     317        end if
     318        call write_diagevo('tsurf'//trim(num),'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/))
    373319        if (do_soil) then
    374320            if (icetable_equilibrium) then
    375                 call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
    376                 call write_diagevo('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
     321                call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
     322                call write_diagevo('icetable_thick'//trim(num),'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/))
    377323            else if (icetable_dynamic) then
    378                 call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
    379                 call write_diagevo('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
     324                call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/))
     325                call write_diagevo('ice_porefilling'//trim(num),'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/))
    380326            end if
    381             call write_diagevo('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
    382             call write_diagevo('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
     327            call write_diagevo('tsoil_avg'//trim(num),'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/))
     328            call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/))
    383329            if (do_sorption) then
    384                 call write_diagevo('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
    385                 call write_diagevo('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
     330                call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
     331                call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/))
    386332            end if
    387333        end if
    388334    end do
     335    if (co2ice_flow .and. allocated(is_co2ice_flow)) deallocate(is_co2ice_flow)
     336    if (h2oice_flow .and. allocated(is_h2oice_flow)) deallocate(is_h2oice_flow)
    389337
    390338    ! Checking mass balance for CO2
    … …  
    458406    end if
    459407end do ! End of the evolution loop
    460 deallocate(is_co2ice_disappeared)
    461408
    462409! Finalization
    … …  
    474421
    475422! Deallocation
    476 deallocate(num)
    477423if (do_layering) then
    478424    deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current)
    … …  
    483429    end do
    484430end if
    485 deallocate(layerings_map)
    486 deallocate(h2o_ads_reg,co2_ads_reg)
    487 deallocate(h2o_ice,co2_ice,is_h2oice_ini,is_co2ice_ini)
    488 deallocate(ps_avg,ps_dev,ps_ts)
    489 deallocate(tsurf_avg,tsurf_dev,h2o_surfdensity_avg)
    490 deallocate(tsoil_avg,tsoil_dev,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg)
    491 deallocate(q_h2o_ts,q_co2_ts,q_co2_ts_ini)
    492 deallocate(d_h2oice,d_co2ice,d_co2ice_ini)
    493 deallocate(delta_h2o_ads,delta_co2_ads,delta_icetable,icetable_depth_old)
    494 deallocate(icetable_depth,icetable_thickness,ice_porefilling)
    495431call end_allocation()
    496432
Note: See TracChangeset for help on using the changeset viewer.