Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4152) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4157) @@ -8,10 +8,14 @@ ! ! AUTHORS & DATE -! R. Vandemeulebrouck, 22/07/2022 +! R. Vandemeulebrouck, 22/07/2022 with r2778 & r2779 ! L. Lange, 22/07/2022 ! JB Clement, 2023-2025 ! ! NOTES -! +! Ownership criterion for declarations: +! - Declare in module "planet" variables that are persistent climate +! state or persistent references used across multiple time steps. +! - Declare in program "pem.F90" transient workflow/control varaibles, +! one-shot initialization buffers and per-iteration working buffers. !----------------------------------------------------------------------- @@ -25,5 +29,5 @@ use atmosphere, only: ps_PCM, evolve_pressure, CO2cond_ps_PCM use backup, only: save_clim_state, backup_rate -use clim_state_init, only: read_start, read_startfi, read_startpem +use clim_state_init, only: read_start, read_startfi, read_startevo use config, only: read_rundef, read_display_config use display, only: print_ini, print_end, print_msg, LVL_NFO, LVL_WRN @@ -32,9 +36,10 @@ use glaciers, only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers use ice_table, only: icetable_equilibrium, icetable_dynamic, evolve_ice_table -use layered_deposits, only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering +use layered_deposits, only: do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering, print_layering use maths, only: pi use numerics, only: dp, qp, di, li, k4, minieps, minieps_qp use orbit, only: evo_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit use output, only: write_diagevo, dim_ngrid, dim_nsoil +use planet use physics, only: g use slopes, only: subslope_dist, def_slope_mean @@ -60,63 +65,22 @@ ! --------------- ! Utility-related: -integer(li) :: cr ! Number of clock ticks per second (count rate) -integer(li) :: c1, c2 ! Counts of processor clock -character(:), allocatable :: num ! To write slope variables -integer(di) :: i, islope -! Pressure-related: -real(dp), dimension(:), allocatable :: ps_avg ! Average surface pressure [Pa] -real(dp), dimension(:), allocatable :: ps_dev ! Deviation of surface pressure [Pa] -real(dp), dimension(:,:), allocatable :: ps_ts ! Surface pressure timeseries [Pa] -real(dp) :: ps_avg_glob_ini ! Global average pressure at initialization [Pa] -real(dp) :: ps_avg_glob_old ! Global average pressure of previous time step [Pa] -real(dp) :: ps_avg_glob ! Global average pressure of current time step [Pa] +integer(li) :: cr ! Number of clock ticks per second (count rate) +integer(li) :: c1, c2 ! Counts of processor clock +character(8) :: num ! Slope suffix to ouput variables +integer(di) :: i, islope ! Ice-related: -real(dp), dimension(:,:), allocatable :: h2o_ice ! H2O ice [kg.m-2] -real(dp), dimension(:,:), allocatable :: co2_ice ! CO2 ice [kg.m-2] -real(dp) :: h2oice_sublim_coverage_ini ! Initial surface area of sublimating H2O ice [m2] -real(dp) :: co2ice_sublim_coverage_ini ! Initial surface area of sublimating CO2 ice [m2] -logical(k4), dimension(:,:), allocatable :: is_h2oice_ini ! Initial location of H2O ice -logical(k4), dimension(:,:), allocatable :: is_co2ice_ini ! Initial location of CO2 ice -logical(k4), dimension(:,:), allocatable :: is_co2ice_flow ! Flag for location of CO2 glacier flow -logical(k4), dimension(:,:), allocatable :: is_h2oice_flow ! Flag for location of H2O glacier flow -real(dp), dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_h2ofrost ! Minimum of H2O frost over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_co2frost ! Minimum of CO2 frost over the last PCM year [kg.m-2] -logical(k4), dimension(:,:), allocatable :: is_co2ice_disappeared ! Flag to check if CO2 ice disappeared at the previous timestep +logical(k4), dimension(:,:), allocatable :: is_co2ice_flow ! Flag for location of CO2 glacier flow +logical(k4), dimension(:,:), allocatable :: is_h2oice_flow ! Flag for location of H2O glacier flow +real(dp), dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg.m-2] +real(dp), dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2] +real(dp), dimension(:,:,:), allocatable :: minPCM_h2ofrost ! Minimum of H2O frost over the last PCM year [kg.m-2] +real(dp), dimension(:,:,:), allocatable :: minPCM_co2frost ! Minimum of CO2 frost over the last PCM year [kg.m-2] ! Surface-related: -real(dp), dimension(:,:), allocatable :: tsurf_avg ! Average surface temperature [K] -real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K] -real(dp), dimension(:,:), allocatable :: tsurf_dev ! Deviation of surface temperature [K] -real(dp), dimension(:,:), allocatable :: h2o_surfdensity_avg ! Average water surface density [kg/m^3] -real(dp), dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] -real(dp), dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] -! Soil-related: -real(dp), dimension(:,:,:), allocatable :: tsoil_avg ! Average soil temperature [K] -real(dp), dimension(:,:,:), allocatable :: tsoil_dev ! Deviation pf soil temperature [K] -real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts ! Soil temperature timeseries [K] -real(dp), dimension(:,:,:,:), allocatable :: tsoil_ts_old ! Soil temperature timeseries at the previous time step [K] -real(dp), dimension(:,:,:), allocatable :: h2o_soildensity_avg ! Average of soil water soil density [kg/m^3] -real(dp), dimension(:), allocatable :: delta_co2_ads ! Quantity of CO2 exchanged due to adsorption/desorption [kg/m^2] -real(dp), dimension(:), allocatable :: delta_h2o_ads ! Quantity of H2O exchanged due to adsorption/desorption [kg/m^2] -real(qp) :: totmass_adsh2o ! Total mass of H2O exchanged because of adsorption/desorption [kg] -real(dp), dimension(:,:,:), allocatable :: h2o_ads_reg ! H2O adsorbed in the regolith [kg/m^2] -real(dp), dimension(:,:,:), allocatable :: co2_ads_reg ! CO2 adsorbed in the regolith [kg/m^2] -real(dp), dimension(:,:), allocatable :: icetable_depth ! Depth of the ice table [m] -real(dp), dimension(:,:), allocatable :: icetable_thickness ! Thickness of the ice table [m] -real(dp), dimension(:,:,:), allocatable :: ice_porefilling ! Amount of porefilling [m^3/m^3] -real(dp), dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table [m] -real(dp), dimension(:), allocatable :: delta_icetable ! Total mass of the H2O that has sublimated / condenses from the ice table [kg] -! Tracer-related: -real(dp), dimension(:,:), allocatable :: q_co2_ts ! CO2 mass mixing ratio in the first layer [kg/kg] -real(dp), dimension(:,:), allocatable :: q_co2_ts_ini ! Initial CO2 mass mixing ratio in the first layer [kg/kg] -real(dp), dimension(:,:), allocatable :: q_h2o_ts ! H2O mass mixing ratio in the first layer [kg/kg] +real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K] +real(dp), dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] +real(dp), dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] ! Tendency-related: -real(dp), dimension(:,:), allocatable :: d_co2ice ! Tendency of perennial CO2 ice [kg/m2/y] -real(dp), dimension(:,:), allocatable :: d_co2ice_ini ! Tendency of perennial CO2 ice at the beginning [kg/m2/y] -real(dp), dimension(:,:), allocatable :: d_h2oice ! Tendency of perennial H2O ice [kg/m2/y] real(dp), dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendency of perennial H2O ice to keep balance between donor and recipient [kg/m2/y] ! Layering-related: -type(layering), dimension(:,:), allocatable :: layerings_map ! Layering for each grid point and slope type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering real(dp), dimension(:,:), allocatable :: h2oice_depth ! Depth of subsurface ice layer @@ -130,4 +94,5 @@ real(qp) :: totmass_co2ice_ini, totmass_atmco2_ini, totmass_adsco2_ini ! Initial total CO2 masses (surface ice|atmospheric|adsorbed) real(qp) :: totmass_ini ! Initial total CO2 mass [kg] +real(qp) :: totmass_adsh2o ! Current total adsorbed H2O mass [kg] real(qp) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O ice reservoirs @@ -168,9 +133,10 @@ ! Read the PCM data given by XIOS -allocate(ps_avg(ngrid),ps_ts(ngrid,nday)) -allocate(tsurf_avg(ngrid,nslope),tsurf_avg_yr1(ngrid,nslope),h2o_surfdensity_avg(ngrid,nslope)) -allocate(tsoil_avg(ngrid,nsoil,nslope),tsoil_ts(ngrid,nsoil,nslope,nday),h2o_soildensity_avg(ngrid,nsoil,nslope)) -allocate(q_h2o_ts(ngrid,nday),q_co2_ts(ngrid,nday)) -allocate(minPCM_h2operice(ngrid,nslope,2),minPCM_co2perice(ngrid,nslope,2),minPCM_h2ofrost(ngrid,nslope,2),minPCM_co2frost(ngrid,nslope,2)) +call allocate_xios_state() +allocate(tsurf_avg_yr1(ngrid,nslope)) +allocate(minPCM_h2operice(ngrid,nslope,2)) +allocate(minPCM_co2perice(ngrid,nslope,2)) +allocate(minPCM_h2ofrost(ngrid,nslope,2)) +allocate(minPCM_co2frost(ngrid,nslope,2)) call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, & @@ -181,5 +147,5 @@ ! Compute the deviation from the average -allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoil_PCM,nslope)) +call allocate_deviation_state() ps_dev(:) = ps_PCM(:) - ps_avg(:) tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:) @@ -193,14 +159,6 @@ ! Read the "startevo.nc" -allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,nslope)) -allocate(icetable_depth(ngrid,nslope),icetable_thickness(ngrid,nslope),ice_porefilling(ngrid,nsoil,nslope)) -allocate(h2o_ads_reg(ngrid,nsoil,nslope),co2_ads_reg(ngrid,nsoil,nslope),delta_h2o_ads(ngrid),delta_co2_ads(ngrid)) -allocate(layerings_map(ngrid,nslope)) -icetable_depth(:,:) = 0._dp -icetable_thickness(:,:) = 0._dp -ice_porefilling(:,:,:) = 0._dp -delta_h2o_ads(:) = 0._dp -delta_co2_ads(:) = 0._dp -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, & +call allocate_startevo_state() +call 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, & tsoil_avg,h2o_soildensity_avg,icetable_depth,icetable_thickness,ice_porefilling,layerings_map, & h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) @@ -208,5 +166,5 @@ ! Compute ice tendencies from yearly minima -allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope)) +call allocate_tendencies() call print_msg('> Computing surface ice tendencies',LVL_NFO) call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice) @@ -218,6 +176,5 @@ ! Save initial set-up useful for the next computations call print_msg('> Saving some initial climate state variables',LVL_NFO) -allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday)) -allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope)) +call allocate_initial_snapshots() ps_avg_glob_ini = ps_avg_glob d_co2ice_ini(:,:) = d_co2ice(:,:) @@ -225,6 +182,4 @@ h2oice_sublim_coverage_ini = 0._dp co2ice_sublim_coverage_ini = 0._dp -is_h2oice_ini(:,:) = .false. -is_co2ice_ini(:,:) = .false. totmass_co2ice_ini = 0._qp totmass_atmco2_ini = 0._qp @@ -268,12 +223,5 @@ !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot end if -if (nslope == 1) then ! No slope - allocate(character(0) :: num) -else ! Using slopes - allocate(character(8) :: num) -end if -allocate(delta_icetable(ngrid),icetable_depth_old(ngrid,nslope),is_co2ice_disappeared(ngrid,nslope),tsoil_ts_old(ngrid,nsoil,nslope,nday)) -is_co2ice_disappeared(:,:) = .false. -delta_icetable(:) = 0._dp +call allocate_loop_state() n_yr_run = 0 idt = 0 @@ -353,38 +301,38 @@ call write_diagevo('ps_avg_glob','Global average pressure','Pa',ps_avg_glob) do islope = 1,nslope - if (nslope /= 1) then - num = ' ' - write(num,'(i2.2)') islope - num = '_slope'//num - end if - call write_diagevo('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/)) - call write_diagevo('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/)) - call write_diagevo('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/)) - call write_diagevo('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/)) + if (nslope == 1) then + num = '' + else + write(num,'("_slope",i2.2)') islope + end if + call write_diagevo('h2oice'//trim(num),'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/)) + call write_diagevo('co2ice'//trim(num),'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/)) + call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/)) + call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/)) if (co2ice_flow) then - call write_diagevo('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/)) - deallocate(is_co2ice_flow) + call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/)) end if if (h2oice_flow) then - call write_diagevo('Flow_h2oice'//num,'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/)) - deallocate(is_h2oice_flow) - end if - call write_diagevo('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/)) + call write_diagevo('Flow_h2oice'//trim(num),'H2O ice flow location','T/F',merge(1._dp,0._dp,is_h2oice_flow(:,islope)),(/dim_ngrid/)) + end if + call write_diagevo('tsurf'//trim(num),'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/)) if (do_soil) then if (icetable_equilibrium) then - call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) - call write_diagevo('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/)) + call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) + call write_diagevo('icetable_thick'//trim(num),'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/)) else if (icetable_dynamic) then - call write_diagevo('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) - call write_diagevo('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('icetable_depth'//trim(num),'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) + call write_diagevo('ice_porefilling'//trim(num),'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/)) end if - call write_diagevo('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/)) - call write_diagevo('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('tsoil_avg'//trim(num),'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/)) if (do_sorption) then - call write_diagevo('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) - call write_diagevo('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) end if end if end do + if (co2ice_flow .and. allocated(is_co2ice_flow)) deallocate(is_co2ice_flow) + if (h2oice_flow .and. allocated(is_h2oice_flow)) deallocate(is_h2oice_flow) ! Checking mass balance for CO2 @@ -458,5 +406,4 @@ end if end do ! End of the evolution loop -deallocate(is_co2ice_disappeared) ! Finalization @@ -474,5 +421,4 @@ ! Deallocation -deallocate(num) if (do_layering) then deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current) @@ -483,14 +429,4 @@ end do end if -deallocate(layerings_map) -deallocate(h2o_ads_reg,co2_ads_reg) -deallocate(h2o_ice,co2_ice,is_h2oice_ini,is_co2ice_ini) -deallocate(ps_avg,ps_dev,ps_ts) -deallocate(tsurf_avg,tsurf_dev,h2o_surfdensity_avg) -deallocate(tsoil_avg,tsoil_dev,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg) -deallocate(q_h2o_ts,q_co2_ts,q_co2_ts_ini) -deallocate(d_h2oice,d_co2ice,d_co2ice_ini) -deallocate(delta_h2o_ads,delta_co2_ads,delta_icetable,icetable_depth_old) -deallocate(icetable_depth,icetable_thickness,ice_porefilling) call end_allocation()