!------------------------ ! I Initialization ! I_a Read the "run.def" ! I_b Read the "start.nc" and "startfi.nc" ! I_c Subslope parametrisation ! I_d Read the PCM data and convert them to the physical grid ! I_e Initialization of the PEM variable and soil ! I_f Compute global surface pressure ! I_g Read the "startpem.nc" ! I_h Compute orbit criterion ! II Run ! II_a Update pressure, ice and tracers ! II_b Evolution of ice ! II_c Flow of glaciers ! II_d Update surface and soil temperatures ! II_e Outputs ! II_f Update the tendencies ! II_g Checking the stopping criterion ! III Output ! III_a Update surface values for the PCM start files ! III_b Write the "restart.nc" and "restartfi.nc" ! III_c Write the "restartpem.nc" !------------------------ PROGRAM pem !----------------------------------------------------------------------- ! NAME ! pem ! ! DESCRIPTION ! Main entry point for the Planetary Evolution Model (PEM). ! ! AUTHORS & DATE ! R. Vandemeulebrouck, 2022 ! L. Lange, 2022 ! JB Clement, 2023-2025 ! ! NOTES ! Drives initialization, run loop, and outputs for PEM/PCM coupling. !----------------------------------------------------------------------- ! DEPENDENCIES ! ------------ use phyetat0_mod, only: phyetat0 use phyredem, only: physdem0, physdem1 use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close use turb_mod, only: q2, wstar use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h use logic_mod, only: iflag_phys use mod_const_mpi, only: COMM_LMDZ use infotrac use geometry_mod, only: latitude_deg use config_pem, only: read_rundef use pemredem, only: pemdem0, pemdem1 use glaciers, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, computeTcondCO2 use stopping_crit, only: stopping_crit_h2o_ice, stopping_crit_co2, stopping_crit_h2o, stopFlags use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2 use surf_ice, only: evol_co2_ice, evol_h2o_ice, balance_h2oice_reservoirs, h2oice_cap_threshold, set_perice4PCM, h2o_ice, co2_ice, ini_surf_ice, end_surf_ice use soil, only: set_soil, do_soil, ini_soil, end_soil, nsoilmx_PEM, & TI_PEM, & ! Soil thermal inertia tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths fluxgeo ! Geothermal flux for the PEM and PCM use sorption, only: regolith_adsorption, do_sorption, & ! Bool to check if adsorption, main subroutine ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays co2_adsorbed_phys, h2o_adsorbed_phys ! Mass of co2 and h2O adsorbed use evolution, only: dt, evol_orbit, nyears_max, convert_years, year_bp_ini use orbit, only: orbit_param_criterion, set_new_orbitpm use ice_table, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, & ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi use soil_thermal_inertia, only: update_soil_TI use time_phylmdz_mod, only: daysec, dtphys use stop_pem, only: stop_clean use tendencies, only: compute_tend, recomp_tend_co2, recomp_tend_h2o use info_pem, only: update_info use pemetat0_mod, only: pemetat0 use xios_data, only: load_xios_data, get_timelen use soil_temp, only: compute_tsoil, shift_tsoil2surf use outputs, only: write_diagpem, write_diagsoilpem use co2condens_mod, only: CO2cond_ps use layered_deposits, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, & rho_co2ice, rho_h2oice, ptrarray, & stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str use subsurface_ice, only: dyn_ss_ice_m use parse_args_mod, only: parse_args use job_timelimit_mod, only: timelimit, antetime, timewall use paleoclimate_mod, only: h2o_ice_depth, zdqsdif_ssi_tot use surf_temp, only: update_tsurf_nearest_baresoil use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil use surfdat_h, only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h, & albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, & zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold, & watercap, watercaptag, perennial_co2ice, albedo_perennialco2 use dimradmars_mod, only: totcloudfrac, albedo use dust_param_mod, only: tauscaling use tracer_mod, only: noms ! Tracer names use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master use planete_h, only: year_day use surfini_mod, only: surfini use metamorphism, only: ini_frost_id, set_frost4PCM, iPCM_h2ofrost, iPCM_co2frost use tracers, only: mmol, iPCM_qh2o, ini_tracers_id, end_tracers use phys_constants, only: pi, g, r, mugaz, rad, cpp, rcp, read_constants #ifndef CPP_1D use iniphysiq_mod, only: iniphysiq use control_mod, only: iphysiq, day_step, nsplit_phys #else use time_phylmdz_mod, only: iphysiq, steps_per_sol use regular_lonlat_mod, only: init_regular_lonlat use physics_distribution_mod, only: init_physics_distribution use mod_grid_phy_lmdz, only: regular_lonlat use init_testphys1d_mod, only: init_testphys1d use comvert_mod, only: ap, bp use writerestart1D_mod, only: writerestart1D #endif ! DECLARATION ! ----------- implicit none include "dimensions.h" include "paramet.h" include "comgeom.h" include "iniprint.h" ! LOCAL VARIABLES ! --------------- ! Same variable names as in the PCM integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm integer, parameter :: nlayer = llm ! Number of vertical layer integer :: ngrid ! Number of physical grid points integer :: nq ! Number of tracer integer :: day_ini ! First day of the simulation real :: pday ! Physical day real :: time_phys ! Same as in PCM real :: ptimestep ! Same as in PCM real :: ztime_fin ! Same as in PCM ! Variables to read "start.nc" character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM ! Dynamic variables real, dimension(ip1jm,llm) :: vcov ! vents covariants real, dimension(ip1jmp1,llm) :: ucov ! vents covariants real, dimension(ip1jmp1,llm) :: teta ! Potential temperature real, dimension(:,:,:), allocatable :: q ! champs advectes real, dimension(:), allocatable :: pdyn ! pressure for the dynamic grid real, dimension(:), allocatable :: ps_start ! surface pressure in the start file real, dimension(:), allocatable :: ps_start0 ! surface pressure in the start file at the beginning real, dimension(:), allocatable :: ps_avg ! (ngrid) Average surface pressure real, dimension(:), allocatable :: ps_dev ! (ngrid x timelen) Surface pressure deviation real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure real, dimension(ip1jmp1,llm) :: masse ! Air mass real, dimension(ip1jmp1) :: phis ! geopotentiel au sol real :: time_0 ! Variables to read startfi.nc character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM character(2) :: str2 integer :: ncid, status ! Variable for handling opening of files integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files integer :: apvarid, bpvarid ! Variable ID for Netcdf files ! Variables to read starfi.nc and write restartfi.nc real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi real :: total_surface ! Total surface of the planet ! Variables for h2o ice evolution real, dimension(:,:), allocatable :: d_h2oice ! physical point x slope field: Tendency of evolution of perennial h2o ice real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice logical, dimension(:,:), allocatable :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating real :: ps_avg_global_ini ! constant: Global average pressure at initialization [Pa] real :: ps_avg_global_old ! constant: Global average pressure of previous time step real :: ps_avg_global_new ! constant: Global average pressure of current time step real, dimension(:,:), allocatable :: zplev_new ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] real, dimension(:,:), allocatable :: zplev_start0 ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2] real, dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] real, dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step type(stopFlags) :: stopCrit ! Stopping criteria real :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] integer :: timelen ! # time samples real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 real, dimension(:,:), allocatable :: d_h2oice_new ! Adjusted tendencies to keep balance between donor and recipient real, dimension(:,:), allocatable :: min_h2oice ! Yearly minimum of H2O ice for the last year of the PCM real :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O ! Variables for co2 ice evolution real, dimension(:,:), allocatable :: d_co2ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year real, dimension(:,:), allocatable :: d_co2ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM? real :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini ! Logical array to know if co2 ice is sublimating real, dimension(:,:), allocatable :: vmr_co2_PCM ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Grid points x Times co2 volume mixing ratio used in the PEM real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] real, dimension(:,:), allocatable :: min_co2ice ! Yearly minimum of CO2 ice for the last year of the PCM real(kind = 16) :: totmass_co2ice, totmass_atmco2 ! Current total CO2 masses real(kind = 16) :: totmass_co2ice_ini, totmass_atmco2_ini ! Initial total CO2 masses ! Variables for the evolution of layered layerings_map type(layering), dimension(:,:), allocatable :: layerings_map ! Layering for each grid point and slope type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering logical, dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm real, dimension(:,:), allocatable :: h2o_ice_depth_old ! Old depth of subsurface ice layer ! Variables for slopes integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow ! Variables for surface and soil real, dimension(:,:), allocatable :: tsurf_avg ! Grid points x Slope field: Average surface temperature [K] real, dimension(:,:), allocatable :: tsurf_dev ! Grid points x Slope field: Surface temperature deviation [K] real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Grid points x Slope field: Average surface temperature of first call of the PCM [K] real, dimension(:,:,:), allocatable :: tsoil_avg ! Grid points x Soil x Slope field: Average Soil Temperature [K] real, dimension(:,:), allocatable :: tsoil_avg_old ! Grid points x Soil field: Average Soil Temperature at the previous time step [K] real, dimension(:,:,:), allocatable :: tsoil_dev ! Grid points x Soil x Slope field: Soil temperature deviation [K] real, dimension(:,:,:,:), allocatable :: tsoil_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K] real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K] real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K] real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3] real, dimension(:,:), allocatable :: watersurf_density_avg ! Grid points x Slope: Average water surface density [kg/m^3] real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3] real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Grid points x Soil x Slopes: Average water soil density [kg/m^3] real, dimension(:), allocatable :: delta_co2_adsorbed ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] real, dimension(:), allocatable :: delta_h2o_adsorbed ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] real(kind = 16) :: totmass_adsco2, totmass_adsco2_ini ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] real :: totmass_adsh2o ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m] real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] real, dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table ! Some variables for the PEM run real, parameter :: year_step = 1 ! Timestep for the pem real :: i_myear_leg ! Number of iteration real :: nyears_maxorb ! Maximum number of iterations before stopping real :: i_myear ! Global number of Martian years of the chained simulations real :: nyears_tot ! Maximum number of Martian years of the chained simulations real :: timestep ! Timestep [s] integer(kind = 8) :: cr ! Number of clock ticks per second (count rate) integer(kind = 8) :: c1, c2 ! Counts of processor clock character(8) :: date character(10) :: time character(5) :: zone integer, dimension(8) :: values character(128) :: dir = ' ' character(32) :: logname = ' ' character(32) :: hostname = ' ' #ifdef CPP_1D integer :: nsplit_phys integer, parameter :: jjm_value = jjm - 1 integer :: day_step ! Dummy variables to use the subroutine 'init_testphys1d' logical :: therestart1D, therestartfi, ctrl_h2ovap, ctrl_h2oice integer :: ndt, day0 real :: ptif, pks, day, gru, grv, relaxtime_h2ovap, relaxtime_h2oice real, dimension(:), allocatable :: zqsat real, dimension(:,:,:), allocatable :: dq, dqdyn real, dimension(nlayer) :: play, w, qref_h2ovap, qref_h2oice real, dimension(nlayer + 1) :: plev #else integer, parameter :: jjm_value = jjm real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart real, dimension(:,:), allocatable :: p ! Grid points x Atmosphere: pressure to recompute and write in restart (ip1jmp1,llmp1) #endif integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap real :: totmass_ini logical :: num_str ! CODE ! ---- ! Elapsed time with system clock call system_clock(count_rate = cr) call system_clock(c1) ! Parse command-line options call parse_args() ! Header write(*,*) ' * . . + . * . + . . . ' write(*,*) ' + _______ ________ ____ ____ * + ' write(*,*) ' + . * |_ __ \|_ __ ||_ \ / _| . *' write(*,*) ' . . | |__) | | |_ \_| | \/ | * * . ' write(*,*) ' . | ___/ | _| _ | |\ /| | . . ' write(*,*) '. * * _| |_ _| |__/ | _| |_\/_| |_ * ' write(*,*) ' + |_____| |________||_____||_____| + . ' write(*,*) ' . * . * . + * . + .' ! Some user info call date_and_time(date,time,zone,values) call getcwd(dir) ! Current directory call getlog(logname) ! User name call hostnm(hostname) ! Machine/station name write(*,*) write(*,*) '********* PEM information *********' write(*,*) '+ User : '//trim(logname) write(*,*) '+ Machine : '//trim(hostname) write(*,*) '+ Directory: '//trim(dir) write(*,'(a,i2,a,i2,a,i4)') ' + Date : ',values(3),'/',values(2),'/',values(1) write(*,'(a,i2,a,i2,a,i2,a)') ' + Time : ',values(5),':',values(6),':',values(7) ! Parallel variables is_sequential = .true. is_parallel = .false. is_mpi_root = .true. is_omp_root = .true. is_master = .true. ! Some constants day_ini = 0 time_phys = 0. ngrid = ngridmx A = (1./m_co2 - 1./m_noco2) B = 1./m_noco2 year_day = 669 daysec = 88775. timestep = year_day*daysec*year_step !----------------------------- INITIALIZATION -------------------------- !------------------------ ! I Initialization ! I_a Read the "run.def" !------------------------ write(*,*) write(*,*) '********* PEM initialization *********' write(*,*) '> Reading "run.def" (PEM)' #ifndef CPP_1D dtphys = daysec/48. ! Dummy value (overwritten in phyetat0) call conf_gcm(99,.true.) call infotrac_init nq = nqtot allocate(q(ip1jmp1,llm,nqtot)) allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) #else allocate(q(1,llm,nqtot),pdyn(1)) allocate(longitude(1),latitude(1),cell_area(1)) therestart1D = .false. ! Default value inquire(file = 'start1D.txt',exist = therestart1D) if (.not. therestart1D) then write(*,*) 'There is no "start1D.txt" file!' error stop 'Initialization cannot be done for the 1D PEM.' endif therestartfi = .false. ! Default value inquire(file = 'startfi.nc',exist = therestartfi) if (.not. therestartfi) then write(*,*) 'There is no "startfi.nc" file!' error stop 'Initialization cannot be done for the 1D PEM.' endif write(*,*) '> Reading "start1D.txt" and "startfi.nc"' call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, & time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & play,plev,latitude,longitude,cell_area, & ctrl_h2ovap,relaxtime_h2ovap,qref_h2ovap,ctrl_h2oice,relaxtime_h2oice,qref_h2oice) nsplit_phys = 1 day_step = steps_per_sol #endif call read_rundef(i_myear,nyears_tot) !------------------------ ! I Initialization ! I_b Read of the "start.nc" and "starfi.nc" !------------------------ call read_constants(startfi_name) ! I_b.1 Read "start.nc" write(*,*) '> Reading "start.nc"' allocate(ps_start0(ngrid)) #ifndef CPP_1D allocate(pdyn(ip1jmp1)) call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0) call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0) call iniconst ! Initialization of dynamical constants (comconst_mod) call inigeom ! Initialization of geometry allocate(ap(nlayer + 1)) allocate(ap(nlayer + 1)) allocate(bp(nlayer + 1)) status = nf90_open(start_name,NF90_NOWRITE,ncid) status = nf90_inq_varid(ncid,"ap",apvarid) status = nf90_get_var(ncid,apvarid,ap) status = nf90_inq_varid(ncid,"bp",bpvarid) status = nf90_get_var(ncid,bpvarid,bp) status = nf90_close(ncid) ! Initialization of physics constants and variables call iniphysiq(iim,jjm,llm,(jjm - 1)*iim + 2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) #else ps_start0(1) = pdyn(1) #endif deallocate(pdyn) ! In the PCM, these values are given to the physic by the dynamic. ! Here we simply read them in the "startfi.nc" file status = nf90_open(startfi_name,NF90_NOWRITE,ncid) status = nf90_inq_varid(ncid,"longitude",lonvarid) status = nf90_get_var(ncid,lonvarid,longitude) status = nf90_inq_varid(ncid,"latitude",latvarid) status = nf90_get_var(ncid,latvarid,latitude) status = nf90_inq_varid(ncid,"area",areavarid) status = nf90_get_var(ncid,areavarid,cell_area) status = nf90_inq_varid(ncid,"soildepth",sdvarid) status = nf90_get_var(ncid,sdvarid,mlayer) status = nf90_close(ncid) ! I_b.2 Read the "startfi.nc" ! First we read the initial state (starfi.nc) write(*,*) '> Reading "startfi.nc"' call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) ! Remove unphysical values of surface tracer where (qsurf < 0.) qsurf = 0. call surfini(ngrid,nslope,qsurf) call ini_frost_id(nqtot,noms) call ini_tracers_id(nqtot,noms) !------------------------ ! I Initialization ! I_c Subslope parametrisation !------------------------ ! Define some slope statistics iflat = 1 do islope = 2,nslope if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope enddo write(*,*) 'Flat slope for islope = ',iflat write(*,*) 'Corresponding criterium = ',def_slope_mean(iflat) !------------------------ ! I Initialization ! I_d Read the PCM data and convert them to the physical grid !------------------------ ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value call get_timelen("Xoutdaily4pem_Y1.nc",timelen) allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen)) allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid)) allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope)) allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope)) allocate(min_co2ice(ngrid,nslope),min_h2oice(ngrid,nslope)) call load_xios_data(ngrid,nslope,nsoilmx,timelen,qsurf(:,iPCM_h2ofrost,:),qsurf(:,iPCM_co2frost,:),ps_avg,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_timeseries,watersurf_density_avg, & d_h2oice,d_co2ice,ps_timeseries,q_h2o_PEM_phys,q_co2_PEM_phys,watersoil_density_timeseries,min_h2oice,min_co2ice) vmr_co2_PCM = q_co2_PEM_phys/(A*q_co2_PEM_phys + B)/m_co2 d_co2ice_ini = d_co2ice ! Compute the deviation from the average allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope)) ps_dev = ps_start0 - ps_avg tsurf_dev = tsurf - tsurf_avg tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:) !------------------------ ! I Initialization ! I_e Initialization of the PEM variables and soil !------------------------ call end_soil call ini_soil(ngrid,nslope) call end_adsorption_h_PEM call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) call end_ice_table call ini_ice_table(ngrid,nslope,nsoilmx_PEM) allocate(tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen),watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) if (do_soil) then call set_soil(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries do l = nsoilmx + 1,nsoilmx_PEM tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:) watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:) tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:) enddo watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen endif !do_soil deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries) !------------------------ ! I Initialization ! I_f Compute global surface pressure !------------------------ total_surface = sum(cell_area) ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface ps_avg_global_new = ps_avg_global_ini write(*,*) "Total surface of the planet =", total_surface write(*,*) "Initial global average pressure =", ps_avg_global_ini !------------------------ ! I Initialization ! I_g Read the "startpem.nc" !------------------------ write(*,*) '> Reading "startpem.nc"' call ini_surf_ice(ngrid,nslope) allocate(layerings_map(ngrid,nslope)) allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid)) delta_h2o_icetablesublim = 0. call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, & tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice, & watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,layerings_map) deallocate(tsurf_avg_yr1) ! We save the places where h2o ice is sublimating ! We compute the surface of h2o ice sublimating allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope)) co2ice_sublim_surf_ini = 0. h2oice_ini_surf = 0. is_co2ice_sublim_ini = .false. is_h2oice_sublim_ini = .false. is_co2ice_ini = .false. co2ice_disappeared = .false. totmass_co2ice_ini = 0. totmass_atmco2_ini = 0. if (layering_algo) then do ig = 1,ngrid do islope = 1,nslope if (is_co2ice_str(layerings_map(ig,islope)%top)) then co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice else if (is_h2oice_str(layerings_map(ig,islope)%top)) then h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice else call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope)) endif enddo enddo ! We put the sublimating tendency coming from subsurface ice into the overall tendency where (h2o_ice_depth > 0. .and. zdqsdif_ssi_tot < -1.e-10) d_h2oice = zdqsdif_ssi_tot endif 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.) if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true. if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then is_co2ice_sublim_ini(i,islope) = .true. co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope) endif if (d_h2oice(i,islope) < 0.) then if (h2o_ice(i,islope) > 0.) then is_h2oice_sublim_ini(i,islope) = .true. h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) else if (h2o_ice_depth(i,islope) > 0.) then is_h2oice_sublim_ini(i,islope) = .true. endif endif enddo enddo write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf totmass_adsco2_ini = 0. totmass_adsh2o = 0. if (do_sorption) then do ig = 1,ngrid do islope = 1,nslope do l = 1,nsoilmx_PEM - 1 if (l == 1) then totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) else totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) endif enddo enddo enddo totmass_adsco2 = totmass_adsco2_ini write(*,*) "Tot mass of CO2 in the regolith =", totmass_adsco2 write(*,*) "Tot mass of H2O in the regolith =", totmass_adsh2o endif ! adsorption !------------------------ ! I Initialization ! I_h Compute orbit criterion !------------------------ if (evol_orbit) call orbit_param_criterion(i_myear,nyears_maxorb) !-------------------------- END INITIALIZATION ------------------------- !-------------------------------- RUN ---------------------------------- !------------------------ ! II Run ! II_a Update pressure, ice and tracers !------------------------ write(*,*) write(*,*) '********* PEM cycle *********' i_myear_leg = 0 if (layering_algo) then allocate(h2o_ice_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 ig = 1,ngrid current(ig,islope)%p => layerings_map(ig,islope)%top enddo enddo endif do while (i_myear_leg < min(nyears_maxorb,nyears_max) .and. i_myear < nyears_tot) ! II.a.1. Compute updated global pressure write(*,'(a,f10.2)') ' **** Iteration of the PEM leg (Martian years): ', i_myear_leg + dt write(*,*) "> Updating the surface pressure" ps_avg_global_old = ps_avg_global_new do i = 1,ngrid do islope = 1,nslope ps_avg_global_new = ps_avg_global_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface enddo enddo if (do_sorption) then do i = 1,ngrid ps_avg_global_new = ps_avg_global_new - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface enddo endif ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old write(*,*) 'Global average pressure old time step:',ps_avg_global_old write(*,*) 'Global average pressure new time step:',ps_avg_global_new ! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption) write(*,*) "> Updating the surface pressure timeseries for the new pressure" allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen)) do l = 1,nlayer + 1 do ig = 1,ngrid zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) enddo enddo ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old write(*,*) "> Updating the pressure levels timeseries for the new pressure" allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen)) do l = 1,nlayer + 1 do ig = 1,ngrid zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) enddo enddo ! II.a.3. Tracers timeseries write(*,*) "> Updating the tracer VMR timeseries for the new pressure" allocate(vmr_co2_PEM_phys(ngrid,timelen)) l = 1 do ig = 1,ngrid do t = 1,timelen ! H2O q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ & (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) if (q_h2o_PEM_phys(ig,t) < 0) then q_h2o_PEM_phys(ig,t) = 1.e-30 else if (q_h2o_PEM_phys(ig,t) > 1) then q_h2o_PEM_phys(ig,t) = 1. endif ! CO2 q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ & (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) & + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) & - (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ & (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) if (q_co2_PEM_phys(ig,t) < 0) then q_co2_PEM_phys(ig,t) = 1.e-30 else if (q_co2_PEM_phys(ig,t) > 1) then q_co2_PEM_phys(ig,t) = 1. endif mmean = 1./(A*q_co2_PEM_phys(ig,t) + B) vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 enddo enddo deallocate(zplev_timeseries_new,zplev_timeseries_old) !------------------------ ! II Run ! II_b Evolution of ice !------------------------ allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) if (layering_algo) then h2o_ice_depth_old = h2o_ice_depth allocate(d_h2oice_new(ngrid,nslope)) call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,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(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new) do islope = 1,nslope do ig = 1,ngrid call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice_new(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p) !call print_layering(layerings_map(ig,islope)) co2_ice(ig,islope) = 0. h2o_ice(ig,islope) = 0. h2o_ice_depth(ig,islope) = 0. if (is_co2ice_str(layerings_map(ig,islope)%top)) then co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice*rho_co2ice else if (is_h2oice_str(layerings_map(ig,islope)%top)) then h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice*rho_h2oice else call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope),co2_ice(ig,islope)) endif enddo enddo deallocate(d_h2oice_new) else zlag = 0. call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit) call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) endif !------------------------ ! II Run ! II_c Flow of glaciers !------------------------ allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope)) if (co2ice_flow .and. nslope > 1) then call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, & ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow) if (layering_algo) then do islope = 1,nslope do ig = 1,ngrid layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)/rho_co2ice enddo enddo endif endif if (h2oice_flow .and. nslope > 1) then call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow) if (layering_algo) then do islope = 1,nslope do ig = 1,ngrid layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)/rho_h2oice enddo enddo endif endif !------------------------ ! II Run ! II_d Update surface and soil temperatures !------------------------ ! II_d.1 Update Tsurf call update_tsurf_nearest_baresoil(ngrid,nslope,iim,jjm - 1,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared) if (do_soil) then ! II_d.2 Shifting soil temperature to surface call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) ! II_d.3 Update soil temperature write(*,*)"> Updating soil temperature profile" allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen)) tsoil_PEM_timeseries_old = tsoil_PEM_timeseries do islope = 1,nslope tsoil_avg_old = tsoil_PEM(:,:,islope) call compute_tsoil(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) call compute_tsoil(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) do t = 1,timelen do ig = 1,ngrid do isoil = 1,nsoilmx_PEM ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil) ! Update of watersoil density watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)**mmol(iPCM_qh2o)/(mugaz*r) if (isnan(tsoil_PEM(ig,isoil,islope))) call stop_clean("PEM - Update Tsoil","NaN detected in tsoil_PEM",1) enddo enddo enddo enddo watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen deallocate(tsoil_avg_old) ! II_d.4 Update the ice table allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope)) if (icetable_equilibrium) then write(*,*) "> Updating ice table (equilibrium method)" icetable_thickness_old = icetable_thickness call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere else if (icetable_dynamic) then write(*,*) "> Updating ice table (dynamic method)" ice_porefilling_old = ice_porefilling icetable_depth_old = icetable_depth allocate(porefill(nsoilmx_PEM)) do ig = 1,ngrid do islope = 1,nslope call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps_avg(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth) icetable_depth(ig,islope) = ssi_depth ice_porefilling(ig,:,islope) = porefill enddo enddo deallocate(porefill) call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere endif deallocate(icetable_thickness_old,ice_porefilling_old) ! II_d.5 Update the soil thermal properties call update_soil_TI(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,ps_avg_global_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) ! II_d.6 Update the mass of the regolith adsorbed totmass_adsco2 = 0. totmass_adsh2o = 0. if (do_sorption) then call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, & tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed) do ig = 1,ngrid do islope = 1,nslope do l = 1,nsoilmx_PEM if (l == 1) then totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) else totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) endif enddo enddo enddo write(*,*) "Total mass of CO2 in the regolith =", totmass_adsco2 write(*,*) "Total mass of H2O in the regolith =", totmass_adsh2o endif endif !do_soil deallocate(zshift_surf,zlag) !------------------------ ! II Run ! II_e Outputs !------------------------ call write_diagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/)) do islope = 1,nslope write(str2(1:2),'(i2.2)') islope call write_diagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope)) call write_diagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope)) call write_diagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope)) call write_diagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope)) call write_diagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope))) call write_diagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope))) call write_diagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope)) if (icetable_equilibrium) then call write_diagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) call write_diagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope)) else if (icetable_dynamic) then call write_diagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) endif if (do_soil) then call write_diagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope)) call write_diagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope)) if (icetable_dynamic) call write_diagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope)) if (do_sorption) then call write_diagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope)) call write_diagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope)) endif endif enddo deallocate(flag_co2flow,flag_h2oflow) ! Checking mass balance for CO2 if (abs(CO2cond_ps - 1.) < 1.e-10) then totmass_co2ice = 0. totmass_atmco2 = 0. do ig = 1,ngrid totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g do islope = 1,nslope totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)*cell_area(ig)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) enddo enddo totmass_ini = max(totmass_atmco2_ini + totmass_co2ice_ini + totmass_adsco2_ini,1.e-10) write(*,'(a,f8.3,a)') " > Relative total CO2 mass balance = ", 100.*(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) then write(*,*) ' /!\ Warning: mass balance is not conserved!' totmass_ini = max(totmass_atmco2_ini,1.e-10) write(*,*) ' Atmospheric CO2 mass balance = ', 100.*(totmass_atmco2 - totmass_atmco2_ini)/totmass_ini, ' %' totmass_ini = max(totmass_co2ice_ini,1.e-10) write(*,*) ' CO2 ice mass balance = ', 100.*(totmass_co2ice - totmass_co2ice_ini)/totmass_ini, ' %' totmass_ini = max(totmass_adsco2_ini,1.e-10) write(*,*) ' Adsorbed CO2 mass balance = ', 100.*(totmass_adsco2 - totmass_adsco2_ini)/totmass_ini, ' %' endif endif !------------------------ ! II Run ! II_f Update the tendencies !------------------------ call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new) write(*,*) "> Updating the H2O sub-surface ice tendency due to lag layer" if (layering_algo) then do ig = 1,ngrid do islope = 1,nslope if (is_h2oice_sublim_ini(ig,islope) .and. h2o_ice_depth(ig,islope) > 0.) call recomp_tend_h2o(h2o_ice_depth_old(ig,islope),h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) enddo enddo !~ else !~ do ig = 1,ngrid !~ do islope = 1,nslope !~ call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) !~ enddo !~ enddo endif if (do_soil) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old) deallocate(vmr_co2_PEM_phys) !------------------------ ! II Run ! II_g Checking the stopping criterion !------------------------ write(*,*) "> Checking the stopping criteria" call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid) call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopCrit,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope) i_myear_leg = i_myear_leg + dt i_myear = i_myear + dt if (stopCrit%stop_code() == 0 .and. i_myear_leg >= nyears_max) stopCrit%nyears_max_reached = .true. if (stopCrit%stop_code() == 0 .and. i_myear_leg >= nyears_maxorb) stopCrit%nyears_maxorb_reached = .true. if (stopCrit%stop_code() == 0 .and. i_myear >= nyears_tot) stopCrit%nyears_tot_reached = .true. call system_clock(c2) if (stopCrit%stop_code() == 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopCrit%time_limit_reached = .true. if (stopCrit%is_any_set()) then write(*,*) stopCrit%stop_message() exit else write(*,'(a,f10.2,a)') ' **** The chained simulation has run for ',i_myear,' Martian years.' write(*,*) '**** The PEM can continue!' write(*,*) '****' endif enddo ! big time iteration loop of the pem deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed) deallocate(watersoil_density_PEM_avg,watersurf_density_avg) deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries) deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim) deallocate(d_co2ice,d_co2ice_ini,d_h2oice) deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini) if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current) !------------------------------ END RUN -------------------------------- !------------------------------- OUTPUT -------------------------------- !------------------------ ! III Output ! III_a Update surface values for the PCM start files !------------------------ write(*,*) write(*,*) '********* PEM finalization *********' ! III_a.1 Ice update for start file call set_frost4PCM(qsurf) call set_perice4PCM(ngrid,nslope,qsurf,watercaptag,watercap,perennial_co2ice) ! III.a.3. Tsurf update for start file write(*,*) '> Reconstructing the surface temperature for the PCM' tsurf = tsurf_avg + tsurf_dev deallocate(tsurf_dev) ! III_a.4 Tsoil update for start file if (do_soil) then write(*,*) '> Reconstructing the soil temperature profile for the PCM' inertiesoil = TI_PEM(:,:nsoilmx,:) ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value do isoil = 1,nsoilmx tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:) enddo tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev flux_geo = fluxgeo endif deallocate(tsurf_avg,tsoil_dev) ! III_a.5 Pressure update for start file write(*,*) '> Reconstructing the pressure for the PCM' allocate(ps_start(ngrid)) ! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini deallocate(ps_avg,ps_dev) ! III_a.6 Tracers update for start file write(*,*) '> Reconstructing the tracer VMR for the PCM' allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1)) do l = 1,nlayer + 1 zplev_start0(:,l) = ap(l) + bp(l)*ps_start0 zplev_new(:,l) = ap(l) + bp(l)*ps_start enddo do nnq = 1,nqtot if (noms(nnq) /= "co2") then do l = 1,llm - 1 do ig = 1,ngrid q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) enddo q(:,llm,nnq) = q(:,llm - 1,nnq) enddo else do l = 1,llm - 1 do ig = 1,ngrid q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_start0(ig,l) - zplev_start0(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) enddo q(:,llm,nnq) = q(:,llm - 1,nnq) enddo endif enddo deallocate(zplev_start0) ! Conserving the tracers mass for start file do nnq = 1,nqtot do ig = 1,ngrid do l = 1,llm - 1 if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1)) q(ig,l,nnq) = 1. q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" endif if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 enddo enddo enddo deallocate(zplev_new) ! III_a.7 Albedo update for start file write(*,*) '> Reconstructing the albedo for the PCM' do ig = 1,ngrid if (latitude(ig) < 0.) then icap = 2 ! Southern hemisphere else icap = 1 ! Northern hemisphere endif do islope = 1,nslope ! Bare ground albedo(ig,:,islope) = albedodat(ig) emis(ig,islope) = emissiv ! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant ! H2O ice if (h2o_ice(ig,islope) > 0.) then albedo(ig,:,islope) = albedo_h2o_cap emis(ig,islope) = 1. endif ! CO2 ice if (co2_ice(ig,islope) > 0.) then albedo(ig,:,islope) = albedo_perennialco2(icap) emis(ig,islope) = emisice(icap) endif ! H2O frost if (qsurf(ig,iPCM_h2ofrost,islope) > 0.) then albedo(ig,:,islope) = albedo_h2o_frost emis(ig,islope) = 1. endif ! CO2 frost if (qsurf(ig,iPCM_co2frost,islope) > 0.) then albedo(ig,:,islope) = albedice(icap) emis(ig,islope) = emisice(icap) endif enddo enddo ! III_a.8 Orbital parameters update for start file write(*,*) '> Setting the new orbital parameters' if (evol_orbit) call set_new_orbitpm(i_myear,i_myear_leg) !------------------------ ! III Output ! III_b Write "restart.nc" and "restartfi.nc" !------------------------ ! III_b.1 Write "restart.nc" ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys pday = day_ini ztime_fin = time_phys #ifndef CPP_1D write(*,*) '> Writing "restart.nc"' ! Correction on teta due to surface pressure changes allocate(pdyn(ip1jmp1)) call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn) do i = 1,ip1jmp1 teta(i,:) = teta(i,:)*pdyn(i)**rcp enddo ! Correction on atmospheric pressure allocate(p(ip1jmp1,nlayer + 1)) call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn) call pression(ip1jmp1,ap,bp,pdyn,p) ! Correction on the mass of atmosphere call massdair(p,masse) call dynredem0("restart.nc",day_ini,phis) call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn) deallocate(ap,bp,p,pdyn) #else write(*,*) '> Writing "restart1D.txt"' call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q) #endif deallocate(ps_start0,ps_start) ! III_b.2 Write the "restartfi.nc" write(*,*) '> Writing "restartfi.nc"' call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & inertiedat,def_slope,subslope_dist) call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, & wstar,watercap,perennial_co2ice) !------------------------ ! III Output ! III_c Write the "restartpem.nc" !------------------------ write(*,*) '> Writing "restartpem.nc"' if (layering_algo) nb_str_max = get_nb_str_max(layerings_map,ngrid,nslope) ! Get the maximum number of "stratum" in the layerings_mapication (layerings) call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,co2_ice,layerings_map) call update_info(i_myear_leg,stopCrit%stop_code(),i_myear,nyears_tot) write(*,*) write(*,*) '****** PEM final information *******' write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years." write(*,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years." write(*,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years." write(*,*) "+ PEM: so far, so good!" write(*,*) '************************************' if (layering_algo) then do islope = 1,nslope do i = 1,ngrid call del_layering(layerings_map(i,islope)) enddo enddo endif deallocate(q,longitude,latitude,cell_area,tsoil_PEM) call end_surf_ice() deallocate(layerings_map) call end_tracers() !----------------------------- END OUTPUT ------------------------------ END PROGRAM pem