PROGRAM pem !----------------------------------------------------------------------- ! NAME ! pem ! ! DESCRIPTION ! Main entry point for the Planetary Evolution Model (PEM). ! ! AUTHORS & DATE ! R. Vandemeulebrouck, 22/07/2022 ! L. Lange, 22/07/2022 ! JB Clement, 2023-2025 ! ! NOTES ! !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ ! Common modules use job_timelimit_mod, only: timelimit, antetime, timewall use parse_args_mod, only: parse_args ! PEM modules use allocation, only: ini_allocation, end_allocation use atmosphere, only: ps_PCM, evolve_pressure, CO2cond_ps_PCM, build4PCM_atmosphere use clim_state_init, only: read_start, read_startfi, read_startpem use clim_state_rec, only: write_restart, write_restartfi, write_restartpem use config, only: read_rundef use display, only: print_ini, print_end, print_msg use evolution, only: n_yr_run, n_yr_sim, ntot_yr_sim, nmax_yr_run, dt, idt, r_plnt2earth_yr, pem_ini_date use geometry, only: ngrid, nslope, nday, nsoil_PCM, nsoil, cell_area, total_surface, nlayer use glaciers, only: h2oice_flow, co2ice_flow, flow_co2glaciers, flow_h2oglaciers use ice_table, only: icetable_equilibrium, icetable_dynamic, icetable_depth, icetable_thickness, ice_porefilling, evolve_ice_table use layered_deposits, only: layering, do_layering, del_layering, evolve_layering, ptrarray, layering2surfice, surfice2layering use maths, only: pi use numerics, only: dp, qp, di, li, k4, minieps, minieps_qp use orbit, only: evol_orbit, read_orbitpm, compute_maxyr_orbit, update_orbit use output, only: write_diagevol, dim_ngrid, dim_nsoil use physics, only: g use slopes, only: subslope_dist, def_slope_mean use soil, only: do_soil, set_soil, TI, build4PCM_soil use soil_temp, only: tsoil_PCM, shift_tsoil2surf, evolve_soil_temp use soil_therm_inertia, only: update_soil_TI use sorption, only: do_sorption, compute_totmass_adsorbed, evolve_regolith_adsorption use stopping_crit, only: stopFlags, stopping_crit_pressure, stopping_crit_h2o, stopping_crit_h2oice, stopping_crit_co2ice use surface, only: emissivity_PCM, build4PCM_surf_rad_prop use surf_ice, only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs, build4PCM_perice use surf_temp, only: tsurf_PCM, adapt_tsurf2disappearedice, build4PCM_tsurf use tendencies, only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice use tracers, only: adapt_tracers2pressure, build4PCM_tracers, nq use utility, only: real2str use workflow_status, only: read_workflow_status, update_workflow_status use xios_data, only: load_xios_data ! DECLARATION ! ----------- implicit none ! LOCAL VARIABLES ! --------------- ! 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] real(dp), dimension(:), allocatable :: ps4PCM ! Surface pressure reconstruction to feed back into PCM [Pa] real(dp), dimension(:,:), allocatable :: teta4PCM ! Potential temperature reconstruction to feed back into PCM [K] real(dp), dimension(:,:), allocatable :: air_mass4PCM ! Air mass reconstruction to feed back into PCM [kg] real(dp) :: pa4PCM ! Transition pressure (for hybrid coord.) [Pa] real(dp) :: preff4PCM ! Reference pressure [Pa] ! 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] real(dp), dimension(:,:), allocatable :: h2o_ice4PCM ! H2O ice reconstruction to feed back into PCM [kg.m-2] real(dp), dimension(:,:), allocatable :: co2_ice4PCM ! CO2 ice reconstruction to feed back into PCM [kg.m-2] logical(k4), dimension(:), allocatable :: is_h2o_perice ! Location of H2O infinite reservoir, called 'watercaptag' in PCM logical(k4), dimension(:,:), allocatable :: is_co2ice_disappeared ! Flag to check if CO2 ice disappeared at the previous timestep ! 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] real(dp), dimension(:,:), allocatable :: albedo4PCM ! Albedo reconstruction to feed back into PCM real(dp), dimension(:,:), allocatable :: emissivity4PCM ! Emissivity reconstruction to feed back into PCM real(dp), dimension(:,:), allocatable :: tsurf4PCM ! Surface temperature reconstruction to feed back into PCM [K] ! 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 :: inertiesoil4PCM ! Soil thermal inertia reconstruction to feed back into PCM [SI] real(dp), dimension(:,:,:), allocatable :: tsoil4PCM ! Average soil temperature reconstruction to feed back into PCM [K] real(dp), dimension(:,:), allocatable :: flux_geo4PCM ! Geothermal flux reconstruction to feed back into PCM [W/m2] 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_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 :: q4PCM ! Tracers reconstruction to feed back into PCM [kg/kg] ! 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 real(dp), dimension(:,:), allocatable :: h2oice_depth_old ! Old depth of subsurface ice layer logical(k4), dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm ! Evolution-related: real(dp) :: nmax_yr_runorb ! Maximum number of years for the run due to orbital parameters type(stopFlags) :: stopcrit ! Stopping criteria ! Balance-related real(qp) :: totmass_co2ice, totmass_atmco2, totmass_adsco2 ! Current total CO2 masses (surface ice|atmospheric|adsorbed) 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) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O ice reservoirs ! CODE ! ---- ! Pre-processing step !~~~~~~~~~~~~~~~~~~~~ ! Elapsed time with system clock call system_clock(count_rate = cr) call system_clock(c1) ! Parse command-line options call parse_args() ! Initialization ! ~~~~~~~~~~~~~~ ! Header call print_ini() ! Allocate module arrays call ini_allocation() ! Read the duration information of the workflow call read_workflow_status() ! Read the PEM parameters call read_rundef() ! Read the orbital parameters call read_orbitpm(n_yr_sim) ! Read the "start.nc" call read_start() ! Read the "startfi.nc" call read_startfi() ! 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 load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, & q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost) ! Initiate soil settings and TI if (do_soil) call set_soil(TI) ! Compute the deviation from the average allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoil_PCM,nslope)) ps_dev(:) = ps_PCM(:) - ps_avg(:) tsurf_dev(:,:) = tsurf_PCM(:,:) - tsurf_avg(:,:) tsoil_dev(:,:,:) = tsoil_PCM(:,:,:) - tsoil_avg(:,1:nsoil_PCM,:) ! Compute global surface pressure ps_avg_glob = sum(cell_area*ps_avg)/total_surface ! Compute the accepted maximum number of years due to orbital parameters (if needed) call compute_maxyr_orbit(n_yr_sim,nmax_yr_runorb) ! Read the "startpem.nc" allocate(h2o_ice(ngrid,nslope),co2_ice(ngrid,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)) call read_startpem(tsurf_avg_yr1,tsurf_avg,ps_avg_glob_ini,ps_ts,q_co2_ts,q_h2o_ts,h2o_surfdensity_avg,d_h2oice,d_co2ice,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) deallocate(tsurf_avg_yr1) ! Compute ice tendencies from yearly minima allocate(d_h2oice(ngrid,nslope),d_co2ice(ngrid,nslope)) call print_msg('> Computing surface ice tendencies') call compute_tendice(minPCM_h2operice + minPCM_h2ofrost,h2o_ice(:,:) > 0._dp,d_h2oice) call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice))) call compute_tendice(minPCM_co2perice + minPCM_co2frost,co2_ice(:,:) > 0._dp,d_co2ice) call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice))) deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost) ! Save initial set-up useful for the next computations call print_msg('> Saving some initial climate state variables') allocate(d_co2ice_ini(ngrid,nslope),q_co2_ts_ini(ngrid,nday)) allocate(is_h2oice_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope)) ps_avg_glob_ini = ps_avg_glob d_co2ice_ini(:,:) = d_co2ice(:,:) q_co2_ts_ini(:,:) = q_co2_ts(:,:) 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 totmass_adsco2_ini = 0._qp totmass_adsh2o = 0._qp do i = 1,ngrid totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g do islope = 1,nslope totmass_co2ice_ini = totmass_co2ice_ini + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp) if (co2_ice(i,islope) > 0._dp) then is_co2ice_ini(i,islope) = .true. if (d_co2ice(i,islope) < 0._dp) co2ice_sublim_coverage_ini = co2ice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp) end if if (h2o_ice(i,islope) > 0._dp) then is_h2oice_ini(i,islope) = .true. if (d_h2oice(i,islope) < 0._dp) h2oice_sublim_coverage_ini = h2oice_sublim_coverage_ini + cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp) end if end do end do call print_msg("Initial global average pressure [Pa] = "//real2str(ps_avg_glob_ini)) call print_msg("Initial surface area of sublimating CO2 ice [m2] = "//real2str(co2ice_sublim_coverage_ini)) call print_msg("Initial surface area of sublimating H2O ice [m2] = "//real2str(h2oice_sublim_coverage_ini)) if (do_sorption) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2_ini,totmass_adsh2o) ! Main evolution loop ! ~~~~~~~~~~~~~~~~~~~ call print_msg('') call print_msg('********* Evolution *********') if (do_layering) then allocate(h2oice_depth(ngrid,nslope),h2oice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope)) new_str = .true. new_lag = .true. do islope = 1,nslope do i = 1,ngrid current(i,islope)%p => layerings_map(i,islope)%top end do end do ! Conversion to surface ice call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth) ! We put the sublimating tendency coming from subsurface ice into the overall tendency !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 n_yr_run = 0 idt = 0 do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim) call print_msg('**** Iteration of the PEM run [Planetary years]: '//real2str(n_yr_run + dt)) ! Evolve global surface pressure call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg) ! Adapt tracers according to global surface pressure call adapt_tracers2pressure(ps_avg_glob_old,ps_avg_glob,ps_ts,q_h2o_ts,q_co2_ts) ! Evolve surface ice allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) if (do_layering) then h2oice_depth_old(:,:) = h2oice_depth(:,:) do islope = 1,nslope do i = 1,ngrid call evolve_layering(layerings_map(i,islope),d_co2ice(i,islope),d_h2oice_new(i,islope),new_str(i,islope),zshift_surf(i,islope),new_lag(i,islope),zlag(i,islope),current(i,islope)%p) !call print_layering(layerings_map(i,islope)) end do end do ! Conversion to surface ice call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth) ! Balance H2O ice reservoirs allocate(d_h2oice_new(ngrid,nslope)) call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit) call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new) deallocate(d_h2oice_new) else zlag(:,:) = 0._dp call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit) call evolve_co2ice(co2_ice,d_co2ice,zshift_surf) end if ! Flow glaciers according to surface ice if (co2ice_flow) then allocate(is_co2ice_flow(ngrid,nslope)) call flow_co2glaciers(q_co2_ts,ps_ts,ps_avg_glob_old,ps_avg_glob,co2_ice,is_co2ice_flow) end if if (h2oice_flow) then allocate(is_h2oice_flow(ngrid,nslope)) call flow_h2oglaciers(tsurf_avg,h2o_ice,is_h2oice_flow) end if if (do_layering) call surfice2layering(h2o_ice,co2_ice,layerings_map) ! Adapt surface temperature if surface ice disappeared call adapt_tsurf2disappearedice(co2_ice,is_co2ice_ini,is_co2ice_disappeared,tsurf_avg) !call adapt_tsurf2disappearedice(h2o_ice,is_h2oice_ini,tsurf_avg) if (do_soil) then ! Shift soil temperature to surface call shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf_avg,tsoil_avg) ! Evolve soil temperature call evolve_soil_temp(tsoil_avg,tsurf_avg,tsoil_ts,tsoil_ts_old,h2o_soildensity_avg) ! Evolve ice table call evolve_ice_table(h2o_ice,h2o_surfdensity_avg,h2o_soildensity_avg,tsoil_avg,tsurf_avg,delta_icetable,q_h2o_ts, & ps_avg,icetable_depth,icetable_thickness,ice_porefilling,icetable_depth_old) ! Update soil thermal properties according to ice table and soil temperature call update_soil_TI(d_h2oice,h2o_ice,ps_avg_glob,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI) ! Evolve adsorbed species if (do_sorption) then call evolve_regolith_adsorption(d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_avg,TI,ps_ts,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads) call compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o) else totmass_adsh2o = 0._dp totmass_adsco2 = 0._dp end if end if ! do_soil deallocate(zshift_surf,zlag) ! Output the results call print_msg('> Standard outputs') call write_diagevol('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_diagevol('h2oice'//num,'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/)) call write_diagevol('co2ice'//num,'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/)) call write_diagevol('d_h2oice'//num,'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/)) call write_diagevol('d_co2ice'//num,'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/)) if (co2ice_flow) then call write_diagevol('Flow_co2ice'//num,'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/)) deallocate(is_co2ice_flow) end if if (h2oice_flow) then call write_diagevol('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_diagevol('tsurf'//num,'Surface temperature','K',tsurf_avg(:,islope),(/dim_ngrid/)) if (do_soil) then if (icetable_equilibrium) then call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) call write_diagevol('icetable_thick'//num,'Ice table thickness','m',icetable_thickness(:,islope),(/dim_ngrid/)) else if (icetable_dynamic) then call write_diagevol('icetable_depth'//num,'Ice table depth','m',icetable_depth(:,islope),(/dim_ngrid/)) call write_diagevol('ice_porefilling'//num,'Ice pore filling','-',ice_porefilling(:,:,islope),(/dim_ngrid,dim_nsoil/)) end if call write_diagevol('tsoil_avg'//num,'Soil temperature','K',tsoil_avg(:,:,islope),(/dim_ngrid,dim_nsoil/)) call write_diagevol('inertiesoil'//num,'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/)) if (do_sorption) then call write_diagevol('co2_ads_reg'//num,'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) call write_diagevol('h2o_ads_reg'//num,'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) end if end if end do ! Checking mass balance for CO2 if (abs(CO2cond_ps_PCM - 1._dp) < minieps) then totmass_co2ice = 0._dp totmass_atmco2 = 0._dp do i = 1,ngrid totmass_atmco2 = totmass_atmco2 + cell_area(i)*ps_avg(i)/g do islope = 1,nslope totmass_co2ice = totmass_co2ice + co2_ice(i,islope)*cell_area(i)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.) end do end do totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,minieps_qp) ! To avoid division by 0 call print_msg(" > Relative total CO2 mass balance = "//real2str(100._qp*(totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini)//' %') if (abs((totmass_atmco2 + totmass_co2ice + totmass_adsco2 - totmass_atmco2_ini - totmass_co2ice_ini - totmass_adsco2_ini)/totmass_ini) > 0.01_qp) then call print_msg(' /!\ Warning: mass balance is not conserved!') totmass_ini = max(totmass_atmco2_ini,minieps_qp) ! To avoid division by 0 call print_msg(' Atmospheric CO2 mass balance = '//real2str(100._qp*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini)//' %') totmass_ini = max(totmass_co2ice_ini,minieps_qp) ! To avoid division by 0 call print_msg(' CO2 ice mass balance = '//real2str(100._qp*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini)//' %') totmass_ini = max(totmass_adsco2_ini,minieps_qp) ! To avoid division by 0 call print_msg(' Adsorbed CO2 mass balance = '//real2str(100._qp*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini)//' %') end if end if ! Evolve the tendencies call evolve_tend_co2ice(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice) !~ if (do_layering) then !~ do i = 1,ngrid !~ do islope = 1,nslope !~ if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2oice(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope)) !~ end do !~ end do ! else ! do i = 1,ngrid ! do islope = 1,nslope ! call evolve_tend_h2oice(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope)) ! end do ! end do !~ end if ! Increment time n_yr_run = n_yr_run + dt n_yr_sim = n_yr_sim + dt idt = idt + 1 ! Check the stopping criteria call print_msg("> Checking the stopping criteria") call stopping_crit_pressure(ps_avg_glob_ini,ps_avg_glob,stopcrit) call stopping_crit_h2oice(h2oice_sublim_coverage_ini,h2o_ice,d_h2oice,stopcrit) call stopping_crit_co2ice(co2ice_sublim_coverage_ini,co2_ice,d_co2ice,stopcrit) call system_clock(c2) if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_run) stopcrit%nmax_yr_run_reached = .true. if (stopcrit%stop_code() == 0 .and. n_yr_run >= nmax_yr_runorb) stopcrit%nmax_yr_runorb_reached = .true. if (stopcrit%stop_code() == 0 .and. n_yr_sim >= ntot_yr_sim) stopcrit%nmax_yr_sim_reached = .true. if (stopcrit%stop_code() == 0 .and. timewall .and. real(c2 - c1,dp)/real(cr,dp) >= timelimit - antetime) stopcrit%time_limit_reached = .true. if (stopcrit%is_any_set()) then call print_msg(stopcrit%stop_message()) exit else call print_msg('**** This run has achieved '//real2str(n_yr_run)//' Planetary years.') call print_msg('**** The workflow has achieved '//real2str(n_yr_sim)//' Planetary years.') call print_msg('**** This run is continuing!') call print_msg('****') end if end do ! End of the evolution loop deallocate(is_co2ice_disappeared) ! Finalization ! ~~~~~~~~~~~~ call print_msg('') call print_msg('********* Finalization *********') ! Update orbital parameters if (evol_orbit) call update_orbit(n_yr_sim,n_yr_run) ! Build ice for the PCM allocate(h2o_ice4PCM(ngrid,nslope),co2_ice4PCM(ngrid,nslope),is_h2o_perice(ngrid)) call build4PCM_perice(h2o_ice,co2_ice,is_h2o_perice,h2o_ice4PCM,co2_ice4PCM) ! Build surface temperature for the PCM allocate(tsurf4PCM(ngrid,nslope)) call build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM) ! Build soil for the PCM if (do_soil) then allocate(tsoil4PCM(ngrid,nsoil_PCM,nslope),inertiesoil4PCM(ngrid,nsoil_PCM,nslope),flux_geo4PCM(ngrid,nslope)) call build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM) end if ! Build atmosphere for the PCM allocate(ps4PCM(ngrid),teta4PCM(ngrid,nlayer),air_mass4PCM(ngrid,nlayer)) call build4PCM_atmosphere(ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini,ps4PCM,pa4PCM,preff4PCM,teta4PCM,air_mass4PCM) ! Build tracers for the PCM allocate(q4PCM(ngrid,nlayer,nq)) call build4PCM_tracers(ps4PCM,q4PCM) ! Build surface radiative properties state for the PCM allocate(albedo4PCM(ngrid,nslope),emissivity4PCM(ngrid,nslope)) call build4PCM_surf_rad_prop(h2o_ice,co2_ice,albedo4PCM,emissivity4PCM) ! Write the "startpem.nc" call write_restartpem(h2o_ice,co2_ice,tsoil_avg,TI,icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map) ! Write the "startfi.nc" call write_restartfi(is_h2o_perice,h2o_ice4PCM,co2_ice4PCM,tsurf4PCM,tsoil4PCM,inertiesoil4PCM,albedo4PCM,emissivity4PCM,flux_geo4PCM) ! Write the "start.nc" call write_restart(ps4PCM,pa4PCM,preff4PCM,q4PCM,teta4PCM,air_mass4PCM) ! Update the duration information of the workflow call update_workflow_status(n_yr_run,stopcrit%stop_code(),n_yr_sim,ntot_yr_sim) ! Deallocation deallocate(emissivity4PCM,albedo4PCM) deallocate(q4PCM) deallocate(ps4PCM,teta4PCM,air_mass4PCM) if (do_soil) deallocate(tsoil4PCM,inertiesoil4PCM,flux_geo4PCM) deallocate(tsurf4PCM) deallocate(co2_ice4PCM,h2o_ice4PCM,is_h2o_perice) deallocate(num) if (do_layering) then deallocate(h2oice_depth,h2oice_depth_old,new_str,new_lag,current) do islope = 1,nslope do i = 1,ngrid call del_layering(layerings_map(i,islope)) end do 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) call end_allocation() ! Footer call system_clock(c2) call print_end(n_yr_run,real((c2 - c1),dp)/real(cr,dp),r_plnt2earth_yr,pem_ini_date,n_yr_sim) END PROGRAM pem