!------------------------ ! 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 tendencies ! I_g Save the initial situation ! I_h Read the "startpem.nc" ! I_i 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 value for the PCM start files ! III_b Write the "restart.nc" and "restartfi.nc" ! III_c Write the "restartpem.nc" !------------------------ PROGRAM pem 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 conf_pem_mod, only: conf_pem use pemredem, only: pemdem0, pemdem1 use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, & metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2 use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2 use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, 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 adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays co2_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded use time_evol_mod, only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini use orbit_param_criterion_mod, only: orbit_param_criterion use recomp_orb_param_mod, only: recomp_orb_param use ice_table_mod, 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_thermalproperties_mod, only: update_soil_thermalproperties use time_phylmdz_mod, only: daysec, dtphys use abort_pem_mod, only: abort_pem use soil_settings_PEM_mod, only: soil_settings_PEM use compute_tend_mod, only: compute_tend use info_PEM_mod, only: info_PEM use interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM use nb_time_step_PCM_mod, only: nb_time_step_PCM use pemetat0_mod, only: pemetat0 use read_data_PCM_mod, only: read_data_PCM use recomp_tend_co2_slope_mod, only: recomp_tend_co2_slope use compute_soiltemp_mod, only: compute_tsoil_pem use writediagpem_mod, only: writediagpem, writediagsoilpem use co2condens_mod, only: CO2cond_ps use layering_mod, only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo #ifndef CPP_STD use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & albedo_h2o_frost,frost_albedo_threshold, & emissiv, watercaptag, perennial_co2ice, emisice, albedice use dimradmars_mod, only: totcloudfrac, albedo use dust_param_mod, only: tauscaling use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit use comcstfi_h, only: pi, rad, g, mugaz, r use surfini_mod, only: surfini use comconst_mod, only: kappa, cpp #else use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, & ALBEDO_CO2_ICE_SPECTV, phys_state_var_init use aerosol_mod, only: iniaerosol use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit use comcstfi_mod, only: pi, rad, g, mugaz, r #endif #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 implicit none include "dimensions.h" include "paramet.h" include "comgeom.h" include "iniprint.h" include "callkeys.h" integer ngridmx parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) ! Same variable names as in the PCM 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 ! temperature potentielle real, dimension(:,:,:), allocatable :: q ! champs advectes real, dimension(ip1jmp1) :: ps ! pression au sol real, dimension(ip1jmp1) :: ps0 ! pression au sol initiale real, dimension(:), allocatable :: ps_start_PCM ! (ngrid) surface pressure real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure real, dimension(ip1jmp1,llm) :: masse ! masse d'air real, dimension(ip1jmp1) :: phis ! geopotentiel au sol real :: time_0 ! Variables to read starfi.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 :: h2o_ice ! h2o ice in the PEM real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2] real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice logical, dimension(:,:), allocatable :: ini_h2oice_sublim ! Logical array to know if h2o ice is sublimating real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa] real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step real :: global_avg_press_new ! constant: Global average pressure of current time step real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2] real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step integer :: stopPEM ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu real, save :: 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 ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] integer :: timelen ! # time samples real :: ave ! intermediate varibale to compute average real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 ! Variables for co2_ice evolution real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM? real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2] real :: co2ice_ini_surf ! Initial surface of sublimating co2 ice logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] ! Variables for stratification (layering) evolution type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm ! Variables for slopes real, dimension(:,:), allocatable :: tend_co2_ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year real, dimension(:,:), allocatable :: tend_co2_ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM real, dimension(:,:), allocatable :: tend_h2o_ice ! physical point x slope field: Tendency of evolution of perennial h2o ice real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow ! Variables for surface and soil real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] logical :: bool_sublim ! logical to check if there is sublimation or not logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep real, dimension(:,:), allocatable :: icetable_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table 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 ! Some variables for the PEM run real, parameter :: year_step = 1 ! Timestep for the pem integer :: year_iter ! Number of iteration integer :: year_iter_max ! Maximum number of iterations before stopping integer :: i_myear ! Global number of Martian years of the chained simulations integer :: n_myear ! Maximum number of Martian years of the chained simulations real :: timestep ! Timestep [s] character(20) :: job_id ! Job id provided as argument passed on the command line when the program was invoked logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job integer(kind=8) :: cr ! Number of clock ticks per second (count rate) integer(kind=8) :: c1, c2 ! Counts of processor clock character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command real :: timelimit ! Time limit for the PEM job in seconds real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default integer :: cstat, days, hours, minutes, seconds character(1) :: sep #ifdef CPP_STD real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice real :: albedo_h2o_frost ! albedo of h2o frost real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case #endif #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 integer :: ndt, day0 real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau real, dimension(:), allocatable :: zqsat real, dimension(:,:,:), allocatable :: dq, dqdyn real, dimension(nlayer) :: play, w 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 #endif ! Loop variables integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap ! Elapsed time with system clock call system_clock(count_rate = cr) call system_clock(c1) timewall = .true. timelimit = 86400 ! 86400 seconds = 24 h by default if (command_argument_count() > 0) then ! Read the job id passed as argument to the program call get_command_argument(1,job_id) ! Execute the system command call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat) if (cstat /= 0) then call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat) if (cstat > 0) then error stop 'pem: command execution failed!' else if (cstat < 0) then error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!' endif endif ! Read the output open(1,file = 'tmp_cmdout.txt',status = 'old') read(1,'(a)') chtimelimit close(1) chtimelimit = trim(chtimelimit) call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat) if (cstat > 0) then error stop 'pem: command execution failed!' else if (cstat < 0) then error stop 'pem: command execution not supported!' endif if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS" read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds timelimit = days*86400 + hours*3600 + minutes*60 + seconds else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS" read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds timelimit = hours*3600 + minutes*60 + seconds else ! 'chtimelimit' format is "MM:SS" read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds timelimit = minutes*60 + seconds endif else timewall = .false. endif ! Parallel variables #ifndef CPP_STD is_sequential = .true. is_parallel = .false. is_mpi_root = .true. is_omp_root = .true. is_master = .true. #endif ! Some constants day_ini = 0 ! test time_phys = 0. ! test 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" !------------------------ #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)) 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 call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, & time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) ps(2) = ps(1) nsplit_phys = 1 day_step = steps_per_sol #endif call conf_pem(i_myear,n_myear) !------------------------ ! I Initialization ! I_b Read of the "start.nc" and starfi_evol.nc !------------------------ ! I_b.1 Read "start.nc" allocate(ps_start_PCM(ngrid)) #ifndef CPP_1D call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0) call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM) call iniconst !new call inigeom 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) 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_start_PCM(1) = ps(1) #endif ! Save initial surface pressure ps0 = ps ! 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) #ifndef CPP_STD 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) #else call phys_state_var_init(nq) if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) call initracer(ngrid,nq) call iniaerosol() allocate(tsurf_read_generic(ngrid)) allocate(qsurf_read_generic(ngrid,nq)) allocate(tsoil_read_generic(ngrid,nsoilmx)) allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope)) allocate(emis_read_generic(ngrid)) allocate(tsurf(ngrid,1)) allocate(qsurf(ngrid,nq,1)) allocate(tsoil(ngrid,nsoilmx,1)) allocate(emis(ngrid,1)) allocate(watercap(ngrid,1)) allocate(watercaptag(ngrid)) allocate(albedo_read_generic(ngrid,2)) allocate(albedo(ngrid,2,1)) allocate(inertiesoil(ngrid,nsoilmx,1)) call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, & rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) nslope = 1 call ini_comslope_h(ngrid,1) qsurf(:,:,1) = qsurf_read_generic tsurf(:,1) = tsurf_read_generic tsoil(:,:,1) = tsoil_read_generic emis(:,1) = emis_read_generic watercap(:,1) = 0. watercaptag(:) = .false. albedo(:,1,1) = albedo_read_generic(:,1) albedo(:,2,1) = albedo_read_generic(:,2) inertiesoil(:,:,1) = inertiedat if (nslope == 1) then def_slope(1) = 0 def_slope(2) = 0 def_slope_mean = 0 subslope_dist(:,1) = 1. endif ! Remove unphysical values of surface tracer qsurf(:,:,1) = qsurf_read_generic where (qsurf < 0.) qsurf = 0. #endif do nnq = 1,nqtot ! Why not using ini_tracer ? if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq if (noms(nnq) == "h2o_vap") then igcm_h2o_vap = nnq mmol(igcm_h2o_vap) = 18. endif if (noms(nnq) == "co2") igcm_co2 = nnq enddo r = 8.314511*1000./mugaz !------------------------ ! 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) allocate(flag_co2flow(ngrid,nslope)) allocate(flag_co2flow_mesh(ngrid)) allocate(flag_h2oflow(ngrid,nslope)) allocate(flag_h2oflow_mesh(ngrid)) flag_co2flow = 0 flag_co2flow_mesh = 0 flag_h2oflow = 0 flag_h2oflow_mesh = 0 !------------------------ ! 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 nb_time_step_PCM("data_PCM_Y1.nc",timelen) allocate(tsoil_avg(ngrid,nsoilmx,nslope)) allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope)) allocate(vmr_co2_PCM(ngrid,timelen)) allocate(ps_timeseries(ngrid,timelen)) allocate(min_co2_ice(ngrid,nslope,2)) allocate(min_h2o_ice(ngrid,nslope,2)) allocate(tsurf_avg_yr1(ngrid,nslope)) allocate(tsurf_avg(ngrid,nslope)) allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen)) allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen)) allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) allocate(watersurf_density_avg(ngrid,nslope)) allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(Tsurfavg_before_saved(ngrid,nslope)) allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(delta_co2_adsorbded(ngrid)) allocate(co2ice_disappeared(ngrid,nslope)) allocate(icetable_thickness_prev_iter(ngrid,nslope)) allocate(delta_h2o_icetablesublim(ngrid)) allocate(delta_h2o_adsorbded(ngrid)) allocate(vmr_co2_PEM_phys(ngrid,timelen)) write(*,*) "Downloading data Y1..." call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), & tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & watersurf_density_avg,watersoil_density_timeseries) write(*,*) "Downloading data Y1 done!" ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value write(*,*) "Downloading data Y2..." call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), & tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & watersurf_density_avg,watersoil_density_timeseries) write(*,*) "Downloading data Y2 done!" !------------------------ ! I Initialization ! I_e Initialization of the PEM variables and soil !------------------------ call end_comsoil_h_PEM call ini_comsoil_h_PEM(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) if (soil_pem) then do t = 1,timelen tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - to recompute properly tsoil in the restart enddo call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries do l = nsoilmx + 1,nsoilmx_PEM tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:) watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:) enddo watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen endif !soil_pem deallocate(tsoil_avg) !------------------------ ! I Initialization ! I_f Compute tendencies !------------------------ allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope)) allocate(tend_co2_ice_ini(ngrid,nslope)) ! Compute the tendencies of the evolution of ice over the years call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice) call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice) tend_co2_ice_ini = tend_co2_ice deallocate(min_co2_ice,min_h2o_ice) !------------------------ ! I Initialization ! I_g Save the initial situation !------------------------ allocate(zplev_PCM(ngrid,nlayer + 1)) Total_surface = 0. do ig = 1,ngrid Total_surface = Total_surface + cell_area(ig) zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig) enddo global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface global_avg_press_PCM = global_avg_press_old global_avg_press_new = global_avg_press_old write(*,*) "Total surface of the planet =", Total_surface write(*,*) "Initial global average pressure =", global_avg_press_PCM !------------------------ ! I Initialization ! I_h Read the "startpem.nc" !------------------------ allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope)) allocate(stratif(ngrid,nslope)) if (layering_algo) then do islope = 1,nslope do i = 1,ngrid call ini_layering(stratif(i,islope)) enddo enddo endif 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,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice, & global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, & delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif) ! We save the places where h2o ice is sublimating ! We compute the surface of h2o ice sublimating allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope)) co2ice_ini_surf = 0. h2oice_ini_surf = 0. ini_co2ice_sublim = .false. ini_h2oice_sublim = .false. is_co2ice_ini = .false. co2ice_disappeared = .false. do i = 1,ngrid do islope = 1,nslope if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true. if (tend_co2_ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then ini_co2ice_sublim(i,islope) = .true. co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope) endif if (tend_h2o_ice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then ini_h2oice_sublim(i,islope) = .true. h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) endif enddo enddo write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf delta_h2o_icetablesublim = 0. if (adsorption_pem) then totmassco2_adsorbded = 0. totmassh2o_adsorbded = 0. do ig = 1,ngrid do islope = 1,nslope do l = 1,nsoilmx_PEM - 1 if (l==1) then totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) else totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_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) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_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(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded endif ! adsorption deallocate(tsurf_avg_yr1) !------------------------ ! I Initialization ! I_i Compute orbit criterion !------------------------ #ifndef CPP_STD call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) #else call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) #endif year_iter_max = Max_iter_pem if (evol_orbit_pem) call orbit_param_criterion(i_myear,year_iter_max) !-------------------------- END INITIALIZATION ------------------------- !-------------------------------- RUN ---------------------------------- !------------------------ ! II Run ! II_a Update pressure, ice and tracers !------------------------ year_iter = 0 stopPEM = 0 if (layering_algo) then allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope)) new_str = .true. new_lag1 = .true. new_lag2 = .true. do islope = 1,nslope do ig = 1,ngrid current1(ig,islope)%p => stratif(ig,islope)%top current2(ig,islope)%p => stratif(ig,islope)%top enddo enddo endif do while (year_iter < year_iter_max .and. i_myear < n_myear) ! II.a.1. Compute updated global pressure write(*,*) "Recomputing the new pressure..." do i = 1,ngrid do islope = 1,nslope global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface enddo enddo if (adsorption_pem) then do i = 1,ngrid global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface enddo endif write(*,*) 'Global average pressure old time step',global_avg_press_old write(*,*) 'Global average pressure new time step',global_avg_press_new ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption) allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen)) write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..." do l = 1,nlayer + 1 do ig = 1,ngrid zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) enddo enddo ! II.a.3. Surface pressure timeseries write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..." do ig = 1,ngrid ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old enddo ! II.a.4. New pressure levels timeseries allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen)) write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..." do l = 1,nlayer + 1 do ig = 1,ngrid zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) enddo enddo ! II.a.5. Tracers timeseries write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..." l = 1 do ig = 1,ngrid do t = 1,timelen q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(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 enddo enddo do ig = 1,ngrid do t = 1,timelen q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ & (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & - (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ & (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(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_new_timeseries,zplev_old_timeseries) !------------------------ ! II Run ! II_b Evolution of ice !------------------------ call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM) call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice) if (layering_algo) then do islope = 1,nslope do ig = 1,ngrid call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p) enddo enddo endif !------------------------ ! II Run ! II_c Flow of glaciers !------------------------ if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, & global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh) if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh) !------------------------ ! II Run ! II_d Update surface and soil temperatures !------------------------ ! II_d.1 Update Tsurf write(*,*) "Updating the new Tsurf" bool_sublim = .false. Tsurfavg_before_saved = tsurf_avg do ig = 1,ngrid do islope = 1,nslope if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then ! co2 ice disappeared, look for closest point without co2ice co2ice_disappeared(ig,islope) = .true. if (latitude_deg(ig) > 0) then do ig_loop = ig,ngrid do islope_loop = islope,iflat,-1 if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) bool_sublim = .true. exit endif enddo if (bool_sublim) exit enddo else do ig_loop = ig,1,-1 do islope_loop = islope,iflat if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) bool_sublim = .true. exit endif enddo if (bool_sublim) exit enddo endif if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then albedo(ig,1,islope) = albedo_h2o_frost albedo(ig,2,islope) = albedo_h2o_frost else albedo(ig,1,islope) = albedodat(ig) albedo(ig,2,islope) = albedodat(ig) emis(ig,islope) = emissiv endif else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2 ave = 0. do t = 1,timelen ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.)) enddo tsurf_avg(ig,islope) = ave/timelen endif enddo enddo do t = 1,timelen tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved enddo ! for the start do ig = 1,ngrid do islope = 1,nslope tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope)) enddo enddo if (soil_pem) then ! II_d.2 Update soil temperature write(*,*)"Updating soil temperature" allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) do islope = 1,nslope call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) do t = 1,timelen Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t) Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope) do ig = 1,ngrid do isoil = 1,nsoilmx_PEM watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r) if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1) enddo enddo enddo enddo watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen write(*,*) "Update of soil temperature done" deallocate(Tsoil_locslope) ! II_d.3 Update the ice table if (icetable_equilibrium) then write(*,*) "Compute ice table (equilibrium method)" icetable_thickness_prev_iter = 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_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere else if (icetable_dynamic) then write(*,*) "Compute ice table (dynamic method)" 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(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_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere endif ! II_d.4 Update the soil thermal properties call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbed if (adsorption_pem) then call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice, & h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) totmassco2_adsorbded = 0. totmassh2o_adsorbded = 0. do ig = 1,ngrid do islope = 1,nslope do l = 1,nsoilmx_PEM if (l == 1) then totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) else totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_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) totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_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(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded endif endif !soil_pem !------------------------ ! II Run ! II_e Outputs !------------------------ call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/)) do islope = 1,nslope write(str2(1:2),'(i2.2)') islope call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope)) call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope)) call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope)) call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope)) call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) if (icetable_equilibrium) then call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) else if (icetable_dynamic) then call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope)) endif if (soil_pem) then call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope)) call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope)) if (adsorption_pem) then call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope)) call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope)) endif endif enddo !------------------------ ! II Run ! II_f Update the tendencies !------------------------ call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, & global_avg_press_PCM,global_avg_press_new) !------------------------ ! II Run ! II_g Checking the stopping criterion !------------------------ write(*,*) "Checking the stopping criteria..." call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid) call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope) year_iter = year_iter + dt_pem i_myear = i_myear + dt_pem if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5 if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6 call system_clock(c2) if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7 if (stopPEM > 0) then select case (stopPEM) case(1) write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above." case(2) write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above." case(3) write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above." case(4) write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above." case(5) write(*,*) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM case(6) write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM case(7) write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM case(8) write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM case default write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM end select exit else write(*,*) "The PEM can continue!" write(*,*) "Number of iterations of the PEM: year_iter =", year_iter write(*,*) "Number of simulated Martian years: i_myear =", i_myear endif global_avg_press_old = global_avg_press_new enddo ! big time iteration loop of the pem !------------------------------ END RUN -------------------------------- !------------------------------- OUTPUT -------------------------------- !------------------------ ! III Output ! III_a Update surface value for the PCM start files !------------------------ ! III_a.1 Ice update (for startfi) watercap = 0. perennial_co2ice = co2_ice do ig = 1,ngrid ! H2O ice metamorphism if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold endif ! Is H2O ice still considered as an infinite reservoir for the PCM? if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then ! There is enough ice to be considered as an infinite reservoir watercaptag(ig) = .true. else ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost watercaptag(ig) = .false. qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:) h2o_ice(ig,:) = 0. endif ! CO2 ice metamorphism if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold qsurf(ig,igcm_co2,:) = metam_co2ice_threshold endif enddo ! III_a.2 Tsoil update (for startfi) if (soil_pem) then call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen) #ifndef CPP_STD flux_geo = fluxgeo #endif endif deallocate(tsoil_anom) ! III_a.4 Pressure (for start) ps = ps*global_avg_press_new/global_avg_press_PCM ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM ! III_a.5 Tracer (for start) allocate(zplev_new(ngrid,nlayer + 1)) do l = 1,nlayer + 1 zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM 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_PCM(ig,l) - zplev_PCM(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_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) enddo q(:,llm,nnq) = q(:,llm - 1,nnq) enddo endif enddo ! Conserving the tracers mass for PCM start files 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 if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter) !------------------------ ! 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 allocate(p(ip1jmp1,nlayer + 1)) #ifndef CPP_1D ! Correction on teta due to surface pressure changes do l = 1,nlayer do i = 1,ip1jmp1 teta(i,l)= teta(i,l)*(ps0(i)/ps(i))**kappa enddo enddo ! Correction on atmospheric pressure call pression(ip1jmp1,ap,bp,ps,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,ps) write(*,*) "restart.nc has been written" #else call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q) write(*,*) "restart1D.txt has been written" #endif ! III_b.2 Write the "restartfi.nc" #ifndef CPP_STD 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) #else call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & nlayer,nq,ptimestep,pday,time_phys,cell_area, & albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, & cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, & tsea_ice,sea_ice) #endif write(*,*) "restartfi.nc has been written" !------------------------ ! III Output ! III_c Write the "restartpem.nc" !------------------------ if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (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_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif) write(*,*) "restartpem.nc has been written" call info_PEM(year_iter,stopPEM,i_myear,n_myear) write(*,*) "The PEM has run for", year_iter, "Martian years." write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years." write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years." write(*,*) "LL & RV & JBC: so far, so good!" if (layering_algo) then do islope = 1,nslope do i = 1,ngrid call del_layering(stratif(i,islope)) enddo enddo deallocate(new_str,new_lag1,new_lag2,current1,current2) endif deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys) deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved) deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg) deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter) deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif) !----------------------------- END OUTPUT ------------------------------ END PROGRAM pem