!------------------------ ! I Initialization ! I_a READ run.def ! I_b READ of start_evol.nc and starfi_evol.nc ! I_c Subslope parametrisation ! I_d READ GCM data and convert to the physical grid ! I_e Initialization of the PEM variable and soil ! I_f Compute tendencies & Save initial situation ! I_g Save initial PCM situation ! I_h Read the PEMstart ! I_i Compute orbit criterion ! II Run ! II_a Update pressure, ice and tracers ! II_b Evolution of the ice ! II_c CO2 & H2O glaciers flows ! II_d Update surface and soil temperatures ! II_e Update the tendencies ! II_f Checking the stopping criterion ! III Output ! III_a Update surface value for the PCM start files ! III_b Write restart_evol.nc and restartfi_evol.nc ! III_c Write restartfi_PEM.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, major_slope, ini_comslope_h use logic_mod, only: iflag_phys use mod_const_mpi, only: COMM_LMDZ use comconst_mod, only: rad, g, cpp, pi use infotrac use geometry_mod, only: latitude_deg use conf_pem_mod, only: conf_pem use pemredem, only: pemdem0, pemdem1 use glaciers_mod, only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, & m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial use evol_co2_ice_s_mod, only: evol_co2_ice_s use evol_h2o_ice_s_mod, only: evol_h2o_ice_s use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, & TI_PEM, inertiedat_PEM, & ! soil thermal inertia tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths fluxgeo, & ! Geothermal flux for the PEM and GCM water_reservoir ! Water ressources 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: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, & ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi use soil_thermalproperties_mod, only: update_soil_thermalproperties use time_phylmdz_mod, only: daysec, dtphys, day_end #ifndef CPP_STD use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, & albedodat, zmea, zstd, zsig, zgam, zthe, & hmons, summit, base,albedo_h2o_frost, & frost_albedo_threshold, emissiv, watercaptag, perenial_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 use comcstfi_h, only: r, mugaz use paleoclimate_mod, only: albedo_perenialco2 #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: r, mugaz #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, day_step #endif #ifdef CPP_1D 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_phys_1d_mod, only: init_phys_1d #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 GCM integer :: ngrid ! Number of physical grid points integer :: nlayer ! Number of vertical layer integer :: nq ! Number of tracer integer :: day_ini ! First day of the simulation real :: pday ! Physical day real :: time_phys ! Same as GCM real :: ptimestep ! Same as GCM real :: ztime_fin ! Same as GCM ! Variables to read start.nc character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM ! Dynamic variables real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants real :: teta(ip1jmp1,llm) ! temperature potentielle real, dimension(:,:,:), allocatable :: q ! champs advectes real :: ps(ip1jmp1) ! pression au sol real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées real :: masse(ip1jmp1,llm) ! masse d'air real :: phis(ip1jmp1) ! geopotentiel au sol real :: time_0 ! Variables to read starfi.nc character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM character*2 str2 integer :: ncid, varid,status ! Variable for handling opening of files integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf 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 FILE_NAME and written in restartfi real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi real :: Total_surface ! Total surface of the planet ! Variables for h2o_ice evolution real :: ini_surf_h2o ! Initial surface of sublimating h2o ice real :: ini_surf_co2 ! Initial surface of sublimating co2 ice real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step real :: global_ave_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_gcm ! same but retrieved from the gcm [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 logical :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? logical :: STOPPING_1_water ! Logical : is there still water ice to sublimate? logical :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? logical :: STOPPING_pressure ! Logical : is the criterion (% of change in the surface pressure) reached? integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] real, allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] real, allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM real, allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg] real, allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] integer :: timelen ! # time samples real :: ave ! intermediate varibale to compute average real, 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 slopes real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year real, dimension(:,:), allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial 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, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] real, 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 real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] ! 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] real :: watercap_sum ! total mass of water cap [kg/m^2] real :: water_sum ! total mass of water in the mesh [kg/m^2] #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, allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic real, allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic real, allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic real, allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic real, allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case real, allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case real, allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case real, allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case real, allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model logical, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model real, allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case real, allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case #endif #ifdef CPP_1D integer :: nsplit_phys integer, parameter :: jjm_value = jjm - 1 integer :: ierr #else integer, parameter :: jjm_value = jjm #endif ! Loop variables integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap ! Parallel variables #ifndef CPP_STD is_sequential = .true. is_parallel = .false. is_mpi_root = .true. is_omp_root = .true. is_master = .true. #endif day_ini = 0 ! test time_phys = 0. ! test ! Some constants ngrid = ngridmx nlayer = llm 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 run.def !------------------------ #ifndef CPP_1D dtphys = 0 call conf_gcm(99,.true.) call infotrac_init nq = nqtot allocate(q(ip1jmp1,llm,nqtot)) #else ! load tracer names from file 'traceur.def' open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr) if (ierr /= 0) then write(*,*) 'Cannot find required file "traceur.def"' write(*,*) ' If you want to run with tracers, I need it' write(*,*) ' ... might as well stop here ...' stop else write(*,*) "pem1d: Reading file traceur.def" ! read number of tracers: read(90,*,iostat = ierr) nq write(*,*) "nq",nq nqtot = nq ! set value of nqtot (in infotrac module) as nq if (ierr /= 0) then write(*,*) "pem1d: error reading number of tracers" write(*,*) " (first line of traceur.def) " stop endif if (nq < 1) then write(*,*) "pem1d: error number of tracers" write(*,*) "is nq=",nq," but must be >=1!" stop endif endif nq = nqtot allocate(q(ip1jmp1,llm,nqtot)) allocate(ap(nlayer + 1)) allocate(bp(nlayer + 1)) call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp) pi = 2.*asin(1.) g = 3.72 nsplit_phys = 1 #endif call conf_pem(i_myear,n_myear) !------------------------ ! I Initialization ! I_b READ of start_evol.nc and starfi_evol.nc !------------------------ ! I_b.1 READ start_evol.nc allocate(ps_start_GCM(ngrid)) #ifndef CPP_1D call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0) call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) call iniconst !new call inigeom allocate(ap(nlayer + 1)) allocate(bp(nlayer + 1)) status = nf90_open(FILE_NAME_start,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_GCM(1) = ps(1) #endif ! In the gcm, these values are given to the physic by the dynamic. ! Here we simply read them in the startfi_evol.nc file status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) allocate(longitude(ngrid)) allocate(latitude(ngrid)) allocate(cell_area(ngrid)) 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 startfi_evol.nc ! First we read the initial state (starfi.nc) #ifndef CPP_STD call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, & tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar, & watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist) ! Remove unphysical values of surface tracer do i = 1,ngrid do nnq = 1,nqtot do islope = 1,nslope if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0. enddo enddo enddo call surfini(ngrid,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(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,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, & tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & qsurf_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 do i = 1,ngrid do nnq = 1,nqtot qsurf(i,nnq,1) = qsurf_read_generic(i,nnq) if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0. enddo enddo #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 GCM data and convert 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 GCM run, saving only the minimum value call nb_time_step_GCM("data_GCM_Y1.nc",timelen) allocate(tsoil_ave(ngrid,nsoilmx,nslope)) allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope)) allocate(vmr_co2_gcm(ngrid,timelen)) allocate(ps_timeseries(ngrid,timelen)) allocate(min_co2_ice_1(ngrid,nslope)) allocate(min_h2o_ice_1(ngrid,nslope)) allocate(min_co2_ice_2(ngrid,nslope)) allocate(min_h2o_ice_2(ngrid,nslope)) allocate(tsurf_ave_yr1(ngrid,nslope)) allocate(tsurf_ave(ngrid,nslope)) allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen)) allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(q_co2_PEM_phys(ngrid,timelen)) allocate(q_h2o_PEM_phys(ngrid,timelen)) allocate(co2_ice_GCM(ngrid,nslope,timelen)) allocate(watersurf_density_ave(ngrid,nslope)) allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) allocate(Tsurfave_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(porefillingice_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_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, & tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & co2_ice_GCM,watersurf_density_ave,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 GCM run, saving only the minimum value write(*,*) "Downloading data Y2" call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, & tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & co2_ice_GCM,watersurf_density_ave,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_porefilling call ini_ice_table_porefilling(ngrid,nslope) if (soil_pem) then call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) do l=1,nsoilmx tsoil_PEM(:,l,:)=tsoil_ave(:,l,:) tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:) watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:) enddo do l=nsoilmx+1,nsoilmx_PEM tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:) watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:) enddo watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen endif !soil_pem deallocate(tsoil_ave,tsoil_GCM_timeseries) !------------------------ ! I Initialization ! I_f Compute tendencies & Save initial situation !------------------------ allocate(tendencies_co2_ice(ngrid,nslope)) allocate(tendencies_co2_ice_ini(ngrid,nslope)) allocate(tendencies_h2o_ice(ngrid,nslope)) ! Compute the tendencies of the evolution of ice over the years call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice) tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:) call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice) deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2) !------------------------ ! I Initialization ! I_g Save initial PCM situation !------------------------ allocate(initial_co2_ice_sublim(ngrid,nslope)) allocate(initial_co2_ice(ngrid,nslope)) allocate(initial_h2o_ice(ngrid,nslope)) ! We save the places where water ice is sublimating ! We compute the surface of water ice sublimating ini_surf_co2 = 0. ini_surf_h2o = 0. Total_surface = 0. do i = 1,ngrid Total_surface = Total_surface+cell_area(i) do islope = 1,nslope if (tendencies_co2_ice(i,islope) < 0) then initial_co2_ice_sublim(i,islope) = 1. ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope) else initial_co2_ice_sublim(i,islope) = 0. endif if (qsurf(i,igcm_co2,islope) > 0) then initial_co2_ice(i,islope) = 1. else initial_co2_ice(i,islope) = 0. endif if (tendencies_h2o_ice(i,islope) < 0) then initial_h2o_ice(i,islope) = 1. ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope) else initial_h2o_ice(i,islope) = 0. endif enddo enddo write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2 write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o write(*,*) "Total surface of the planet=", Total_surface allocate(zplev_gcm(ngrid,nlayer + 1)) do l = 1,nlayer + 1 do ig = 1,ngrid zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) enddo enddo global_ave_press_old = 0. do i = 1,ngrid global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface enddo global_ave_press_GCM = global_ave_press_old global_ave_press_new = global_ave_press_old write(*,*) "Initial global average pressure=", global_ave_press_GCM !------------------------ ! I Initialization ! I_h Read the PEMstart !------------------------ call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, & porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, & tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), & qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, & co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) delta_h2o_icetablesublim(:) = 0. do ig = 1,ngrid do islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) enddo enddo if (adsorption_pem) then totmassco2_adsorbded = 0. totmassh2o_adsorbded = 0. do ig = 1,ngrid do islope = 1, nslope do l = 1,nsoilmx_PEM - 1 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - 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+1) - layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 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_ave_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 if (evol_orbit_pem) then call orbit_param_criterion(i_myear,year_iter_max) else year_iter_max = Max_iter_pem endif !-------------------------- END INITIALIZATION ------------------------- !-------------------------------- RUN ---------------------------------- !------------------------ ! II Run ! II_a Update pressure, ice and tracers !------------------------ year_iter = 0 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_ave_press_new = global_ave_press_new-g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface enddo enddo write(*,*) 'Global average pressure old time step',global_ave_press_old call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) if (adsorption_pem) then do i = 1,ngrid global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface enddo write(*,*) 'Global average pressure old time step',global_ave_press_old write(*,*) 'Global average pressure new time step',global_ave_press_new endif ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) 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_ave_press_new/global_ave_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) q_h2o_PEM_phys(ig,t) = 1.e-30 if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1. 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 elseif (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 the ice !------------------------ write(*,*) "Evolution of h2o ice" call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water) write(*,*) "Evolution of co2 ice" call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) do islope=1, nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) enddo !------------------------ ! II Run ! II_c CO2 & H2O glaciers flows !------------------------ write(*,*) "CO2 glacier flows" if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, & global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) write(*,*) "H2O glacier flows" if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh) do islope = 1,nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) enddo !------------------------ ! II Run ! II_d Update surface and soil temperatures !------------------------ ! II_d.1 Update Tsurf write(*,*) "Updating the new Tsurf" bool_sublim = .false. Tsurfave_before_saved(:,:) = tsurf_ave(:,:) do ig = 1,ngrid do islope = 1,nslope if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice if (latitude_deg(ig) > 0) then do ig_loop = ig,ngrid do islope_loop = islope,iflat,-1 if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then tsurf_ave(ig,islope) = tsurf_ave(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(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) bool_sublim = .true. exit endif enddo if (bool_sublim) exit enddo endif initial_co2_ice(ig,islope) = 0 if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,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 ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2 ave = 0. do t = 1,timelen if (co2_ice_GCM(ig,islope,t) > 1.e-3) then ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) else ave = ave + tsurf_GCM_timeseries(ig,islope,t) endif enddo tsurf_ave(ig,islope) = ave/timelen ! set the surface albedo to be the ice albedo if (latitude_deg(ig) > 0) then icap = 1 else icap = 2 endif albedo(ig,1,islope) = albedice(icap) albedo(ig,2,islope) = albedice(icap) emis(ig,islope) = emisice(icap) endif enddo enddo do t = 1,timelen tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) enddo ! for the start do ig = 1,ngrid do islope = 1,nslope tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope)) enddo enddo do islope = 1,nslope write(str2(1:2),'(i2.2)') islope call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) enddo if (soil_pem) then ! II_d.2 Update soil temperature allocate(TI_locslope(ngrid,nsoilmx_PEM)) allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) allocate(Tsurf_locslope(ngrid)) write(*,*)"Updating soil temperature" ! Soil averaged do islope = 1,nslope TI_locslope(:,:) = TI_PEM(:,:,islope) do t = 1,timelen Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t) call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) 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_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) enddo enddo enddo enddo tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen write(*,*) "Update of soil temperature done" deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope) write(*,*) "Compute ice table" ! II_d.3 Update the ice table porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:) call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, & delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere write(*,*) "Update soil propreties" ! II_d.4 Update the soil thermal properties call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, & porefillingice_depth,porefillingice_thickness,TI_PEM) ! II_d.5 Update the mass of the regolith adsorbded if (adsorption_pem) then call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, & qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),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 - 1 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - 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+1) - layer_PEM(l))* & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & cell_area(ig) 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 Update the tendencies !------------------------ write(*,*) "Adaptation of the new co2 tendencies to the current pressure" call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, & global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) !------------------------ ! II Run ! II_f Checking the stopping criterion !------------------------ call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope) year_iter = year_iter + dt_pem i_myear = i_myear + dt_pem write(*,*) "Checking all the stopping criterion." if (STOPPING_water) then write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water criterion_stop = 1 endif if (STOPPING_1_water) then write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water criterion_stop = 1 endif if (STOPPING_co2) then write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 criterion_stop = 2 endif if (STOPPING_pressure) then write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure criterion_stop = 3 endif if (year_iter >= year_iter_max) then write(*,*) "STOPPING because maximum number of iterations reached" criterion_stop = 4 endif if (i_myear >= n_myear) then write(*,*) "STOPPING because maximum number of Martian years to be simulated reached" criterion_stop = 5 endif if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then exit else write(*,*) "We continue!" write(*,*) "Number of iterations of the PEM: year_iter =", year_iter write(*,*) "Number of simulated Martian years: i_myear =", i_myear endif global_ave_press_old=global_ave_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) ! H2O ice do ig = 1,ngrid if (watercaptag(ig)) then watercap_sum = 0. do islope = 1,nslope if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost else ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) qsurf(ig,igcm_h2o_ice,islope)=0. endif watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) watercap(ig,islope) = 0. enddo water_reservoir(ig) = water_reservoir(ig) + watercap_sum endif enddo do ig = 1,ngrid water_sum = 0. do islope = 1,nslope water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) enddo if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1 watercaptag(ig) = .true. water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost do islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.) enddo endif else ! let's check that the infinite source of water disapear if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then watercaptag(ig) = .false. do islope = 1,nslope qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) enddo water_reservoir(ig) = 0. endif endif enddo ! CO2 ice do ig = 1,ngrid do islope = 1,nslope if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) albedo(ig,1,islope) = albedo_perenialco2 albedo(ig,2,islope) = albedo_perenialco2 endif enddo enddo ! III_a.2 Tsoil update (for startfi) if (soil_pem) then call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) endif ! III_a.4 Pressure (for start) do i = 1,ip1jmp1 ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM enddo do i = 1,ngrid ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM enddo ! III_a.5 Tracer (for start) allocate(zplev_new(ngrid,nlayer + 1)) do l = 1,nlayer + 1 do ig = 1,ngrid zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) enddo 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_gcm(ig,l) - zplev_gcm(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_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(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 GCM 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_evol.nc and restartfi_evol.nc !------------------------ ! III_b.1 Write restart_evol.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 call pression (ip1jmp1,ap,bp,ps,p) call massdair(p,masse) call dynredem0("restart_evol.nc",day_ini,phis) call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps) write(*,*) "restart_evol.nc has been written" #else do nnq = 1, nqtot call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) enddo #endif ! III_b.2 Write restartfi_evol.nc #ifndef CPP_STD call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & inertiedat,def_slope,subslope_dist) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & albedo,emis,q2,qsurf,tauscaling,totcloudfrac, & wstar,watercap,perenial_co2ice) #else call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & nlayer,nq,ptimestep,pday,time_phys,cell_area, & albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, & cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & tslab,tsea_ice,sea_ice) #endif write(*,*) "restartfi_evol.nc has been written" !------------------------ ! III Output ! III_c Write restartfi_PEM.nc !------------------------ call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & float(day_ini),0.,nslope,def_slope,subslope_dist) call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, & TI_PEM, porefillingice_depth,porefillingice_thickness, & co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) write(*,*) "restartfi_PEM.nc has been written" call info_run_PEM(year_iter,criterion_stop,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!" deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys) deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved) deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) !----------------------------- END OUTPUT ------------------------------ END PROGRAM pem