[2779] | 1 | !------------------------ |
---|
[3028] | 2 | ! I Initialization |
---|
[3161] | 3 | ! I_a Read the "run.def" |
---|
[3317] | 4 | ! I_b Read the "start.nc" and "startfi.nc" |
---|
[2835] | 5 | ! I_c Subslope parametrisation |
---|
[3161] | 6 | ! I_d Read the PCM data and convert them to the physical grid |
---|
[3028] | 7 | ! I_e Initialization of the PEM variable and soil |
---|
[3143] | 8 | ! I_f Compute tendencies |
---|
[3571] | 9 | ! I_g Compute global surface pressure |
---|
[3384] | 10 | ! I_h Read the "startpem.nc" |
---|
[2835] | 11 | ! I_i Compute orbit criterion |
---|
[2779] | 12 | |
---|
| 13 | ! II Run |
---|
[3028] | 14 | ! II_a Update pressure, ice and tracers |
---|
[3149] | 15 | ! II_b Evolution of ice |
---|
| 16 | ! II_c Flow of glaciers |
---|
[2835] | 17 | ! II_d Update surface and soil temperatures |
---|
[3088] | 18 | ! II_e Outputs |
---|
| 19 | ! II_f Update the tendencies |
---|
| 20 | ! II_g Checking the stopping criterion |
---|
[2779] | 21 | |
---|
| 22 | ! III Output |
---|
[2835] | 23 | ! III_a Update surface value for the PCM start files |
---|
[3317] | 24 | ! III_b Write the "restart.nc" and "restartfi.nc" |
---|
[3161] | 25 | ! III_c Write the "restartpem.nc" |
---|
[2779] | 26 | !------------------------ |
---|
| 27 | |
---|
| 28 | PROGRAM pem |
---|
| 29 | |
---|
[3161] | 30 | use phyetat0_mod, only: phyetat0 |
---|
| 31 | use phyredem, only: physdem0, physdem1 |
---|
| 32 | use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close |
---|
| 33 | use turb_mod, only: q2, wstar |
---|
[3206] | 34 | use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h |
---|
[3161] | 35 | use logic_mod, only: iflag_phys |
---|
| 36 | use mod_const_mpi, only: COMM_LMDZ |
---|
[3028] | 37 | use infotrac |
---|
[3161] | 38 | use geometry_mod, only: latitude_deg |
---|
| 39 | use conf_pem_mod, only: conf_pem |
---|
| 40 | use pemredem, only: pemdem0, pemdem1 |
---|
| 41 | use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, & |
---|
[3571] | 42 | metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2 |
---|
[3161] | 43 | use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2 |
---|
| 44 | use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2 |
---|
| 45 | use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice |
---|
[3571] | 46 | use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, & |
---|
| 47 | TI_PEM, & ! Soil thermal inertia |
---|
[3330] | 48 | tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths |
---|
| 49 | fluxgeo ! Geothermal flux for the PEM and PCM |
---|
[3161] | 50 | use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine |
---|
| 51 | ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays |
---|
[3571] | 52 | co2_adsorbed_phys, h2o_adsorbed_phys ! Mass of co2 and h2O adsorbed |
---|
[3498] | 53 | use time_evol_mod, only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini |
---|
[3161] | 54 | use orbit_param_criterion_mod, only: orbit_param_criterion |
---|
| 55 | use recomp_orb_param_mod, only: recomp_orb_param |
---|
[3493] | 56 | use ice_table_mod, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, & |
---|
| 57 | ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi |
---|
[3161] | 58 | use soil_thermalproperties_mod, only: update_soil_thermalproperties |
---|
[3206] | 59 | use time_phylmdz_mod, only: daysec, dtphys |
---|
[3161] | 60 | use abort_pem_mod, only: abort_pem |
---|
| 61 | use soil_settings_PEM_mod, only: soil_settings_PEM |
---|
| 62 | use compute_tend_mod, only: compute_tend |
---|
| 63 | use info_PEM_mod, only: info_PEM |
---|
[3571] | 64 | use get_timelen_PCM_mod, only: get_timelen_PCM |
---|
[3161] | 65 | use pemetat0_mod, only: pemetat0 |
---|
| 66 | use read_data_PCM_mod, only: read_data_PCM |
---|
[3610] | 67 | use recomp_tend_mod, only: recomp_tend_co2, recomp_tend_h2o |
---|
[3553] | 68 | use compute_soiltemp_mod, only: compute_tsoil_pem, shift_tsoil2surf |
---|
[3181] | 69 | use writediagpem_mod, only: writediagpem, writediagsoilpem |
---|
[3207] | 70 | use co2condens_mod, only: CO2cond_ps |
---|
[3498] | 71 | use layering_mod, only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo |
---|
[3527] | 72 | use dyn_ss_ice_m_mod, only: dyn_ss_ice_m |
---|
[3578] | 73 | use version_info_mod, only: print_version_info |
---|
[3638] | 74 | use paleoclimate_mod, only: h2o_ice_depth, zdqsdif_ssi_tot |
---|
[2985] | 75 | |
---|
[2842] | 76 | #ifndef CPP_STD |
---|
[3206] | 77 | use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil |
---|
[3584] | 78 | use surfdat_h, only: tsurf, qsurf, emis, emissiv, emisice, ini_surfdat_h, & |
---|
| 79 | albedodat, albedice, albedo_h2o_frost, albedo_h2o_cap, & |
---|
| 80 | zmea, zstd, zsig, zgam, zthe, frost_albedo_threshold, & |
---|
| 81 | watercap, watercaptag, perennial_co2ice, albedo_perennialco2 |
---|
[3028] | 82 | use dimradmars_mod, only: totcloudfrac, albedo |
---|
| 83 | use dust_param_mod, only: tauscaling |
---|
[3149] | 84 | use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses |
---|
[3028] | 85 | use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master |
---|
[3096] | 86 | use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit |
---|
[3143] | 87 | use surfini_mod, only: surfini |
---|
[3609] | 88 | use comconst_mod, only: pi, rad, g, cpp, kappa |
---|
[3618] | 89 | use comcstfi_h, only: mugaz |
---|
[2842] | 90 | #else |
---|
[3028] | 91 | use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names |
---|
| 92 | use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & |
---|
| 93 | PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, & |
---|
| 94 | ALBEDO_CO2_ICE_SPECTV, phys_state_var_init |
---|
[3039] | 95 | use aerosol_mod, only: iniaerosol |
---|
| 96 | use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit |
---|
[3442] | 97 | use comcstfi_mod, only: pi, rad, g, mugaz, r |
---|
[2842] | 98 | #endif |
---|
[2985] | 99 | |
---|
[3028] | 100 | #ifndef CPP_1D |
---|
[3618] | 101 | use comconst_mod, only: r |
---|
[3076] | 102 | use iniphysiq_mod, only: iniphysiq |
---|
| 103 | use control_mod, only: iphysiq, day_step, nsplit_phys |
---|
[3019] | 104 | #else |
---|
[3618] | 105 | use comcstfi_h, only: r |
---|
[3386] | 106 | use time_phylmdz_mod, only: iphysiq, steps_per_sol |
---|
[3028] | 107 | use regular_lonlat_mod, only: init_regular_lonlat |
---|
| 108 | use physics_distribution_mod, only: init_physics_distribution |
---|
| 109 | use mod_grid_phy_lmdz, only: regular_lonlat |
---|
[3065] | 110 | use init_testphys1d_mod, only: init_testphys1d |
---|
| 111 | use comvert_mod, only: ap, bp |
---|
[3076] | 112 | use writerestart1D_mod, only: writerestart1D |
---|
[2980] | 113 | #endif |
---|
[2835] | 114 | |
---|
[3076] | 115 | implicit none |
---|
[2980] | 116 | |
---|
[3028] | 117 | include "dimensions.h" |
---|
| 118 | include "paramet.h" |
---|
| 119 | include "comgeom.h" |
---|
| 120 | include "iniprint.h" |
---|
[3039] | 121 | include "callkeys.h" |
---|
[2779] | 122 | |
---|
[3028] | 123 | integer ngridmx |
---|
| 124 | parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) |
---|
[2794] | 125 | |
---|
[3096] | 126 | ! Same variable names as in the PCM |
---|
[3065] | 127 | integer, parameter :: nlayer = llm ! Number of vertical layer |
---|
| 128 | integer :: ngrid ! Number of physical grid points |
---|
| 129 | integer :: nq ! Number of tracer |
---|
| 130 | integer :: day_ini ! First day of the simulation |
---|
| 131 | real :: pday ! Physical day |
---|
[3149] | 132 | real :: time_phys ! Same as in PCM |
---|
| 133 | real :: ptimestep ! Same as in PCM |
---|
[3571] | 134 | real :: ztime_fin ! Same as in PCM |
---|
[2794] | 135 | |
---|
[3571] | 136 | ! Variables to read "start.nc" |
---|
[3317] | 137 | character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM |
---|
[2779] | 138 | |
---|
[3028] | 139 | ! Dynamic variables |
---|
[3065] | 140 | real, dimension(ip1jm,llm) :: vcov ! vents covariants |
---|
| 141 | real, dimension(ip1jmp1,llm) :: ucov ! vents covariants |
---|
[3571] | 142 | real, dimension(ip1jmp1,llm) :: teta ! Potential temperature |
---|
[3065] | 143 | real, dimension(:,:,:), allocatable :: q ! champs advectes |
---|
[3603] | 144 | real, dimension(:), allocatable :: pdyn ! pressure for the dynamic grid |
---|
[3571] | 145 | real, dimension(:), allocatable :: ps_start ! surface pressure in the start file |
---|
| 146 | real, dimension(:), allocatable :: ps_start0 ! surface pressure in the start file at the beginning |
---|
[3673] | 147 | real, dimension(:), allocatable :: ps_avg ! (ngrid) Average surface pressure |
---|
[3571] | 148 | real, dimension(:), allocatable :: ps_dev ! (ngrid x timelen) Surface pressure deviation |
---|
| 149 | real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure |
---|
| 150 | real, dimension(ip1jmp1,llm) :: masse ! Air mass |
---|
[3065] | 151 | real, dimension(ip1jmp1) :: phis ! geopotentiel au sol |
---|
[3028] | 152 | real :: time_0 |
---|
[2779] | 153 | |
---|
[3028] | 154 | ! Variables to read starfi.nc |
---|
[3317] | 155 | character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM |
---|
[3143] | 156 | character(2) :: str2 |
---|
[3206] | 157 | integer :: ncid, status ! Variable for handling opening of files |
---|
| 158 | integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files |
---|
| 159 | integer :: apvarid, bpvarid ! Variable ID for Netcdf files |
---|
[2794] | 160 | |
---|
[3028] | 161 | ! Variables to read starfi.nc and write restartfi.nc |
---|
[3143] | 162 | real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi |
---|
| 163 | real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi |
---|
| 164 | real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi |
---|
[3571] | 165 | real :: total_surface ! Total surface of the planet |
---|
[2897] | 166 | |
---|
[3028] | 167 | ! Variables for h2o_ice evolution |
---|
[3327] | 168 | real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM |
---|
[3571] | 169 | real, dimension(:,:), allocatable :: d_h2oice ! physical point x slope field: Tendency of evolution of perennial h2o ice |
---|
[3327] | 170 | real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2] |
---|
| 171 | real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice |
---|
[3571] | 172 | logical, dimension(:,:), allocatable :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating |
---|
| 173 | real :: ps_avg_global_ini ! constant: Global average pressure at initialization [Pa] |
---|
| 174 | real :: ps_avg_global_old ! constant: Global average pressure of previous time step |
---|
| 175 | real :: ps_avg_global_new ! constant: Global average pressure of current time step |
---|
| 176 | real, dimension(:,:), allocatable :: zplev_new ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] |
---|
| 177 | real, dimension(:,:), allocatable :: zplev_start0 ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2] |
---|
| 178 | real, dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] |
---|
| 179 | real, dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step |
---|
[3327] | 180 | 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 |
---|
[3571] | 181 | real :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] |
---|
| 182 | real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] |
---|
[3327] | 183 | integer :: timelen ! # time samples |
---|
| 184 | real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 |
---|
[2779] | 185 | |
---|
[3130] | 186 | ! Variables for co2_ice evolution |
---|
[3571] | 187 | real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM |
---|
| 188 | real, dimension(:,:), allocatable :: d_co2ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year |
---|
| 189 | real, dimension(:,:), allocatable :: d_co2ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM |
---|
| 190 | logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM? |
---|
| 191 | real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2] |
---|
| 192 | real :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice |
---|
| 193 | logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini ! Logical array to know if co2 ice is sublimating |
---|
| 194 | real, dimension(:,:), allocatable :: vmr_co2_PCM ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] |
---|
| 195 | real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Grid points x Times co2 volume mixing ratio used in the PEM |
---|
| 196 | real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] |
---|
[3130] | 197 | |
---|
[3297] | 198 | ! Variables for stratification (layering) evolution |
---|
| 199 | type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope |
---|
| 200 | type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering |
---|
| 201 | logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm |
---|
| 202 | |
---|
[3028] | 203 | ! Variables for slopes |
---|
[3591] | 204 | integer(kind = 1), dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow |
---|
| 205 | integer(kind = 1), dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow |
---|
[2779] | 206 | |
---|
[3028] | 207 | ! Variables for surface and soil |
---|
[3673] | 208 | real, dimension(:,:), allocatable :: tsurf_avg ! Grid points x Slope field: Average surface temperature [K] |
---|
[3598] | 209 | real, dimension(:,:), allocatable :: tsurf_dev ! Grid points x Slope field: Surface temperature deviation [K] |
---|
[3673] | 210 | real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Grid points x Slope field: Average surface temperature of first call of the PCM [K] |
---|
| 211 | real, dimension(:,:,:), allocatable :: tsoil_avg ! Grid points x Soil x Slope field: Average Soil Temperature [K] |
---|
| 212 | real, dimension(:,:), allocatable :: tsoil_avg_old ! Grid points x Soil field: Average Soil Temperature at the previous time step [K] |
---|
[3571] | 213 | real, dimension(:,:,:), allocatable :: tsoil_dev ! Grid points x Soil x Slope field: Soil temperature deviation [K] |
---|
| 214 | real, dimension(:,:,:,:), allocatable :: tsoil_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K] |
---|
| 215 | real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K] |
---|
[3610] | 216 | real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K] |
---|
[3571] | 217 | real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3] |
---|
[3673] | 218 | real, dimension(:,:), allocatable :: watersurf_density_avg ! Grid points x Slope: Average water surface density [kg/m^3] |
---|
[3571] | 219 | real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3] |
---|
[3673] | 220 | real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Grid points x Soil x Slopes: Average water soil density [kg/m^3] |
---|
[3554] | 221 | real, dimension(:), allocatable :: delta_co2_adsorbed ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] |
---|
| 222 | real, dimension(:), allocatable :: delta_h2o_adsorbed ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] |
---|
| 223 | real :: totmassco2_adsorbed ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] |
---|
| 224 | real :: totmassh2o_adsorbed ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] |
---|
[3493] | 225 | logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep |
---|
[3525] | 226 | real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] |
---|
| 227 | real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m] |
---|
[3493] | 228 | real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] |
---|
| 229 | real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table |
---|
| 230 | real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table |
---|
[3553] | 231 | real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] |
---|
| 232 | real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] |
---|
[3610] | 233 | real, dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table |
---|
[3065] | 234 | |
---|
[3028] | 235 | ! Some variables for the PEM run |
---|
[3616] | 236 | real, parameter :: year_step = 1 ! Timestep for the pem |
---|
| 237 | real :: i_myear_leg ! Number of iteration |
---|
| 238 | real :: n_myear_leg ! Maximum number of iterations before stopping |
---|
| 239 | real :: i_myear ! Global number of Martian years of the chained simulations |
---|
| 240 | real :: n_myear ! Maximum number of Martian years of the chained simulations |
---|
| 241 | real :: timestep ! Timestep [s] |
---|
| 242 | character(100) :: arg ! To read command-line arguments program was invoked |
---|
| 243 | logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job |
---|
| 244 | integer(kind = 8) :: cr ! Number of clock ticks per second (count rate) |
---|
| 245 | integer(kind = 8) :: c1, c2 ! Counts of processor clock |
---|
| 246 | character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command |
---|
| 247 | real :: timelimit ! Time limit for the PEM job in seconds |
---|
[3680] | 248 | real, parameter :: antetime = 3600 ! Anticipation time to prevent reaching the job time limit: 3600 s by default (it should cover the computing time of the reshaping tool) |
---|
[3616] | 249 | integer :: cstat, days, hours, minutes, seconds |
---|
| 250 | character(1) :: sep |
---|
| 251 | character(8) :: date |
---|
| 252 | character(10) :: time |
---|
| 253 | character(5) :: zone |
---|
| 254 | integer, dimension(8) :: values |
---|
| 255 | character(128) :: dir = ' ' |
---|
| 256 | character(32) :: logname = ' ' |
---|
| 257 | character(32) :: hostname = ' ' |
---|
[2779] | 258 | |
---|
[2842] | 259 | #ifdef CPP_STD |
---|
[3571] | 260 | real :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice |
---|
| 261 | real :: albedo_h2o_frost ! Albedo of h2o frost |
---|
[3143] | 262 | real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
| 263 | real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
| 264 | real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
| 265 | real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
| 266 | real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
[3065] | 267 | real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case |
---|
| 268 | real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case |
---|
| 269 | real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case |
---|
| 270 | real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case |
---|
| 271 | real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model |
---|
[3068] | 272 | logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model |
---|
[3065] | 273 | real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case |
---|
| 274 | real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case |
---|
[2842] | 275 | #endif |
---|
| 276 | |
---|
[2980] | 277 | #ifdef CPP_1D |
---|
[3143] | 278 | integer :: nsplit_phys |
---|
| 279 | integer, parameter :: jjm_value = jjm - 1 |
---|
[3386] | 280 | integer :: day_step |
---|
[3065] | 281 | |
---|
| 282 | ! Dummy variables to use the subroutine 'init_testphys1d' |
---|
[3742] | 283 | logical :: therestart1D, therestartfi, prescribed_h2ovap |
---|
[3068] | 284 | integer :: ndt, day0 |
---|
[3742] | 285 | real :: ptif, pks, day, gru, grv, h2ovap_relax_time |
---|
[3068] | 286 | real, dimension(:), allocatable :: zqsat |
---|
| 287 | real, dimension(:,:,:), allocatable :: dq, dqdyn |
---|
[3742] | 288 | real, dimension(nlayer) :: play, w, q_prescribed_h2o_vap |
---|
[3068] | 289 | real, dimension(nlayer + 1) :: plev |
---|
[2980] | 290 | #else |
---|
[3591] | 291 | integer, parameter :: jjm_value = jjm |
---|
| 292 | real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart |
---|
| 293 | real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart |
---|
[3628] | 294 | real, dimension(:,:), allocatable :: p ! Grid points x Atmosphere: pressure to recompute and write in restart (ip1jmp1,llmp1) |
---|
[2980] | 295 | #endif |
---|
| 296 | |
---|
[3028] | 297 | ! Loop variables |
---|
[3584] | 298 | integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap |
---|
[3599] | 299 | logical :: num_str |
---|
[2779] | 300 | |
---|
[3649] | 301 | write(*,*) ' * . . + . * . + . . . ' |
---|
[3603] | 302 | write(*,*) ' + _______ ________ ____ ____ * + ' |
---|
| 303 | write(*,*) ' + . * |_ __ \|_ __ ||_ \ / _| . *' |
---|
| 304 | write(*,*) ' . . | |__) | | |_ \_| | \/ | * * . ' |
---|
| 305 | write(*,*) ' . | ___/ | _| _ | |\ /| | . . ' |
---|
| 306 | write(*,*) '. * * _| |_ _| |__/ | _| |_\/_| |_ * ' |
---|
[3649] | 307 | write(*,*) ' + |_____| |________||_____||_____| + . ' |
---|
| 308 | write(*,*) ' . * . * . + * . + .' |
---|
[3603] | 309 | |
---|
[3363] | 310 | ! Elapsed time with system clock |
---|
| 311 | call system_clock(count_rate = cr) |
---|
| 312 | call system_clock(c1) |
---|
[3394] | 313 | timewall = .true. |
---|
[3363] | 314 | timelimit = 86400 ! 86400 seconds = 24 h by default |
---|
[3574] | 315 | if (command_argument_count() > 0) then ! Get the number of command-line arguments |
---|
| 316 | call get_command_argument(1,arg) ! Read the argument given to the program |
---|
[3599] | 317 | num_str = .true. |
---|
| 318 | do i = 1,len_trim(arg) |
---|
| 319 | if (arg(i:i) < '0' .or. arg(i:i) > '9') then |
---|
| 320 | num_str = .false. |
---|
| 321 | exit |
---|
| 322 | endif |
---|
| 323 | enddo |
---|
[3738] | 324 | |
---|
[3599] | 325 | if (num_str) then ! This is a numeric sting so we considerer this is the job id |
---|
[3574] | 326 | ! Execute the system command |
---|
| 327 | call execute_command_line('squeue -j '//trim(adjustl(arg))//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat) |
---|
| 328 | if (cstat /= 0) then |
---|
| 329 | call execute_command_line('qstat -f '//trim(adjustl(arg))//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat) |
---|
| 330 | if (cstat > 0) then |
---|
| 331 | error stop 'pem: command execution failed!' |
---|
| 332 | else if (cstat < 0) then |
---|
| 333 | error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!' |
---|
| 334 | endif |
---|
| 335 | endif |
---|
| 336 | ! Read the output |
---|
| 337 | open(1,file = 'tmp_cmdout.txt',status = 'old') |
---|
| 338 | read(1,'(a)') chtimelimit |
---|
| 339 | close(1) |
---|
| 340 | chtimelimit = trim(adjustl(chtimelimit)) |
---|
| 341 | call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat) |
---|
| 342 | if (cstat > 0) then |
---|
| 343 | error stop 'pem: command execution failed!' |
---|
| 344 | else if (cstat < 0) then |
---|
| 345 | error stop 'pem: command execution not supported!' |
---|
| 346 | endif |
---|
| 347 | if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS" |
---|
| 348 | read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds |
---|
| 349 | timelimit = days*86400 + hours*3600 + minutes*60 + seconds |
---|
| 350 | else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS" |
---|
| 351 | read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds |
---|
| 352 | timelimit = hours*3600 + minutes*60 + seconds |
---|
| 353 | else ! 'chtimelimit' format is "MM:SS" |
---|
| 354 | read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds |
---|
| 355 | timelimit = minutes*60 + seconds |
---|
| 356 | endif |
---|
[3599] | 357 | else ! Arg is not a numeric string |
---|
| 358 | select case (trim(adjustl(arg))) |
---|
| 359 | case('version') |
---|
| 360 | call print_version_info() |
---|
| 361 | stop |
---|
| 362 | case default |
---|
| 363 | error stop "The argument given to the program is unknown!" |
---|
| 364 | end select |
---|
| 365 | endif |
---|
| 366 | else |
---|
| 367 | timewall = .false. |
---|
[3363] | 368 | endif |
---|
| 369 | |
---|
[3616] | 370 | ! Some user info |
---|
| 371 | call date_and_time(date,time,zone,values) |
---|
| 372 | call getcwd(dir) ! Current directory |
---|
| 373 | call getlog(logname) ! User name |
---|
| 374 | call hostnm(hostname) ! Machine/station name |
---|
| 375 | write(*,*) |
---|
[3673] | 376 | write(*,*) '********* PEM information *********' |
---|
[3616] | 377 | write(*,*) '+ User : '//trim(logname) |
---|
| 378 | write(*,*) '+ Machine : '//trim(hostname) |
---|
| 379 | write(*,*) '+ Directory: '//trim(dir) |
---|
| 380 | write(*,'(a,i2,a,i2,a,i4)') ' + Date : ',values(3),'/',values(2),'/',values(1) |
---|
| 381 | write(*,'(a,i2,a,i2,a,i2,a)') ' + Time : ',values(5),':',values(6),':',values(7) |
---|
| 382 | |
---|
[3028] | 383 | ! Parallel variables |
---|
[2842] | 384 | #ifndef CPP_STD |
---|
[3028] | 385 | is_sequential = .true. |
---|
| 386 | is_parallel = .false. |
---|
| 387 | is_mpi_root = .true. |
---|
| 388 | is_omp_root = .true. |
---|
| 389 | is_master = .true. |
---|
[2842] | 390 | #endif |
---|
[2779] | 391 | |
---|
[3065] | 392 | ! Some constants |
---|
[3591] | 393 | day_ini = 0 |
---|
| 394 | time_phys = 0. |
---|
[3028] | 395 | ngrid = ngridmx |
---|
[3591] | 396 | A = (1./m_co2 - 1./m_noco2) |
---|
| 397 | B = 1./m_noco2 |
---|
[3028] | 398 | year_day = 669 |
---|
| 399 | daysec = 88775. |
---|
[3426] | 400 | timestep = year_day*daysec*year_step |
---|
[2794] | 401 | |
---|
[3028] | 402 | !----------------------------- INITIALIZATION -------------------------- |
---|
[2779] | 403 | !------------------------ |
---|
[3028] | 404 | ! I Initialization |
---|
[3161] | 405 | ! I_a Read the "run.def" |
---|
[2779] | 406 | !------------------------ |
---|
[3616] | 407 | write(*,*) |
---|
[3673] | 408 | write(*,*) '********* PEM initialization *********' |
---|
[3603] | 409 | write(*,*) '> Reading "run.def" (PEM)' |
---|
[2980] | 410 | #ifndef CPP_1D |
---|
[3383] | 411 | dtphys = daysec/48. ! Dummy value (overwritten in phyetat0) |
---|
[3028] | 412 | call conf_gcm(99,.true.) |
---|
| 413 | call infotrac_init |
---|
| 414 | nq = nqtot |
---|
| 415 | allocate(q(ip1jmp1,llm,nqtot)) |
---|
[3065] | 416 | allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) |
---|
[2980] | 417 | #else |
---|
[3603] | 418 | allocate(q(1,llm,nqtot),pdyn(1)) |
---|
[3065] | 419 | allocate(longitude(1),latitude(1),cell_area(1)) |
---|
[3129] | 420 | |
---|
[3143] | 421 | therestart1D = .false. ! Default value |
---|
[3363] | 422 | inquire(file = 'start1D.txt',exist = therestart1D) |
---|
[3129] | 423 | if (.not. therestart1D) then |
---|
[3363] | 424 | write(*,*) 'There is no "start1D.txt" file!' |
---|
[3129] | 425 | error stop 'Initialization cannot be done for the 1D PEM.' |
---|
| 426 | endif |
---|
[3143] | 427 | therestartfi = .false. ! Default value |
---|
[3317] | 428 | inquire(file = 'startfi.nc',exist = therestartfi) |
---|
[3129] | 429 | if (.not. therestartfi) then |
---|
[3317] | 430 | write(*,*) 'There is no "startfi.nc" file!' |
---|
[3129] | 431 | error stop 'Initialization cannot be done for the 1D PEM.' |
---|
| 432 | endif |
---|
| 433 | |
---|
[3603] | 434 | write(*,*) '> Reading "start1D.txt" and "startfi.nc"' |
---|
[3609] | 435 | call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, & |
---|
| 436 | time_0,pdyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & |
---|
[3742] | 437 | play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap) |
---|
[3028] | 438 | nsplit_phys = 1 |
---|
[3399] | 439 | day_step = steps_per_sol |
---|
[2980] | 440 | #endif |
---|
[2779] | 441 | |
---|
[3039] | 442 | call conf_pem(i_myear,n_myear) |
---|
[2779] | 443 | |
---|
[2835] | 444 | !------------------------ |
---|
[3028] | 445 | ! I Initialization |
---|
[3571] | 446 | ! I_b Read of the "start.nc" and "starfi.nc" |
---|
[3028] | 447 | !------------------------ |
---|
[3317] | 448 | ! I_b.1 Read "start.nc" |
---|
[3603] | 449 | write(*,*) '> Reading "start.nc"' |
---|
[3591] | 450 | allocate(ps_start0(ngrid)) |
---|
[2980] | 451 | #ifndef CPP_1D |
---|
[3603] | 452 | allocate(pdyn(ip1jmp1)) |
---|
| 453 | call dynetat0(start_name,vcov,ucov,teta,q,masse,pdyn,phis,time_0) |
---|
[2779] | 454 | |
---|
[3603] | 455 | call gr_dyn_fi(1,iip1,jjp1,ngridmx,pdyn,ps_start0) |
---|
[2897] | 456 | |
---|
[3589] | 457 | call iniconst ! Initialization of dynamical constants (comconst_mod) |
---|
| 458 | call inigeom ! Initialization of geometry |
---|
[2980] | 459 | |
---|
[3028] | 460 | allocate(ap(nlayer + 1)) |
---|
| 461 | allocate(bp(nlayer + 1)) |
---|
[3143] | 462 | status = nf90_open(start_name,NF90_NOWRITE,ncid) |
---|
[3028] | 463 | status = nf90_inq_varid(ncid,"ap",apvarid) |
---|
| 464 | status = nf90_get_var(ncid,apvarid,ap) |
---|
| 465 | status = nf90_inq_varid(ncid,"bp",bpvarid) |
---|
| 466 | status = nf90_get_var(ncid,bpvarid,bp) |
---|
| 467 | status = nf90_close(ncid) |
---|
[2779] | 468 | |
---|
[3591] | 469 | ! Initialization of physics constants and variables (comcstfi_h) |
---|
[3571] | 470 | 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) |
---|
[2980] | 471 | #else |
---|
[3628] | 472 | ps_start0(1) = pdyn(1) |
---|
[2980] | 473 | #endif |
---|
[3603] | 474 | deallocate(pdyn) |
---|
[2980] | 475 | |
---|
[3096] | 476 | ! In the PCM, these values are given to the physic by the dynamic. |
---|
[3317] | 477 | ! Here we simply read them in the "startfi.nc" file |
---|
[3571] | 478 | status = nf90_open(startfi_name,NF90_NOWRITE,ncid) |
---|
[3028] | 479 | status = nf90_inq_varid(ncid,"longitude",lonvarid) |
---|
| 480 | status = nf90_get_var(ncid,lonvarid,longitude) |
---|
| 481 | status = nf90_inq_varid(ncid,"latitude",latvarid) |
---|
| 482 | status = nf90_get_var(ncid,latvarid,latitude) |
---|
| 483 | status = nf90_inq_varid(ncid,"area",areavarid) |
---|
| 484 | status = nf90_get_var(ncid,areavarid,cell_area) |
---|
| 485 | status = nf90_inq_varid(ncid,"soildepth",sdvarid) |
---|
| 486 | status = nf90_get_var(ncid,sdvarid,mlayer) |
---|
| 487 | status = nf90_close(ncid) |
---|
[2963] | 488 | |
---|
[3317] | 489 | ! I_b.2 Read the "startfi.nc" |
---|
[2779] | 490 | ! First we read the initial state (starfi.nc) |
---|
[2842] | 491 | #ifndef CPP_STD |
---|
[3603] | 492 | write(*,*) '> Reading "startfi.nc"' |
---|
[3143] | 493 | call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & |
---|
[3149] | 494 | tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & |
---|
[3130] | 495 | watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) |
---|
[2779] | 496 | |
---|
[3070] | 497 | ! Remove unphysical values of surface tracer |
---|
| 498 | where (qsurf < 0.) qsurf = 0. |
---|
[2885] | 499 | |
---|
[3143] | 500 | call surfini(ngrid,nslope,qsurf) |
---|
[2842] | 501 | #else |
---|
[3028] | 502 | call phys_state_var_init(nq) |
---|
| 503 | if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) |
---|
| 504 | call initracer(ngrid,nq) |
---|
| 505 | call iniaerosol() |
---|
| 506 | allocate(tsurf_read_generic(ngrid)) |
---|
| 507 | allocate(qsurf_read_generic(ngrid,nq)) |
---|
| 508 | allocate(tsoil_read_generic(ngrid,nsoilmx)) |
---|
[3114] | 509 | allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope)) |
---|
[3028] | 510 | allocate(emis_read_generic(ngrid)) |
---|
[3591] | 511 | allocate(albedo_read_generic(ngrid,2)) |
---|
| 512 | allocate(qsurf(ngrid,nq,1)) |
---|
[3028] | 513 | allocate(tsurf(ngrid,1)) |
---|
| 514 | allocate(tsoil(ngrid,nsoilmx,1)) |
---|
| 515 | allocate(emis(ngrid,1)) |
---|
| 516 | allocate(watercap(ngrid,1)) |
---|
| 517 | allocate(watercaptag(ngrid)) |
---|
| 518 | allocate(albedo(ngrid,2,1)) |
---|
| 519 | allocate(inertiesoil(ngrid,nsoilmx,1)) |
---|
[3143] | 520 | call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & |
---|
[3149] | 521 | tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & |
---|
| 522 | qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, & |
---|
[3114] | 523 | rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) |
---|
[3065] | 524 | call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) |
---|
[2842] | 525 | |
---|
[3028] | 526 | nslope = 1 |
---|
| 527 | call ini_comslope_h(ngrid,1) |
---|
[2842] | 528 | |
---|
[3149] | 529 | qsurf(:,:,1) = qsurf_read_generic |
---|
| 530 | tsurf(:,1) = tsurf_read_generic |
---|
| 531 | tsoil(:,:,1) = tsoil_read_generic |
---|
| 532 | emis(:,1) = emis_read_generic |
---|
[3028] | 533 | watercap(:,1) = 0. |
---|
| 534 | watercaptag(:) = .false. |
---|
| 535 | albedo(:,1,1) = albedo_read_generic(:,1) |
---|
| 536 | albedo(:,2,1) = albedo_read_generic(:,2) |
---|
[3149] | 537 | inertiesoil(:,:,1) = inertiedat |
---|
[2842] | 538 | |
---|
[3028] | 539 | if (nslope == 1) then |
---|
| 540 | def_slope(1) = 0 |
---|
| 541 | def_slope(2) = 0 |
---|
| 542 | def_slope_mean = 0 |
---|
| 543 | subslope_dist(:,1) = 1. |
---|
| 544 | endif |
---|
[2842] | 545 | |
---|
[3070] | 546 | ! Remove unphysical values of surface tracer |
---|
[3149] | 547 | qsurf(:,:,1) = qsurf_read_generic |
---|
[3070] | 548 | where (qsurf < 0.) qsurf = 0. |
---|
[3591] | 549 | |
---|
| 550 | deallocate(tsurf_read_generic,qsurf_read_generic,qsoil_read_generic,emis_read_generic) |
---|
[2842] | 551 | #endif |
---|
| 552 | |
---|
[3603] | 553 | do nnq = 1,nqtot ! Why not using ini_tracer? |
---|
[3028] | 554 | if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq |
---|
| 555 | if (noms(nnq) == "h2o_vap") then |
---|
| 556 | igcm_h2o_vap = nnq |
---|
[3143] | 557 | mmol(igcm_h2o_vap) = 18. |
---|
[3028] | 558 | endif |
---|
| 559 | if (noms(nnq) == "co2") igcm_co2 = nnq |
---|
[3065] | 560 | enddo |
---|
[3028] | 561 | |
---|
[2835] | 562 | !------------------------ |
---|
[3028] | 563 | ! I Initialization |
---|
[2835] | 564 | ! I_c Subslope parametrisation |
---|
| 565 | !------------------------ |
---|
[3028] | 566 | ! Define some slope statistics |
---|
| 567 | iflat = 1 |
---|
| 568 | do islope = 2,nslope |
---|
| 569 | if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope |
---|
| 570 | enddo |
---|
| 571 | write(*,*) 'Flat slope for islope = ',iflat |
---|
[3603] | 572 | write(*,*) 'Corresponding criterium = ',def_slope_mean(iflat) |
---|
[2794] | 573 | |
---|
| 574 | !------------------------ |
---|
[3028] | 575 | ! I Initialization |
---|
[3161] | 576 | ! I_d Read the PCM data and convert them to the physical grid |
---|
[3028] | 577 | !------------------------ |
---|
[3096] | 578 | ! 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 |
---|
[3571] | 579 | call get_timelen_PCM("data_PCM_Y1.nc",timelen) |
---|
[2794] | 580 | |
---|
[3591] | 581 | allocate(vmr_co2_PCM(ngrid,timelen),q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen)) |
---|
| 582 | allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid)) |
---|
[3553] | 583 | allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2)) |
---|
[3591] | 584 | allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope)) |
---|
[3599] | 585 | allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_timeseries(ngrid,nsoilmx,nslope,timelen)) |
---|
[3591] | 586 | allocate(watersurf_density_avg(ngrid,nslope),watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) |
---|
[2794] | 587 | |
---|
[3571] | 588 | call read_data_PCM("data_PCM_Y1.nc","data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, & |
---|
| 589 | tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries) |
---|
[2985] | 590 | |
---|
[3571] | 591 | ! Compute the deviation from the average |
---|
[3591] | 592 | allocate(ps_dev(ngrid),tsurf_dev(ngrid,nslope),tsoil_dev(ngrid,nsoilmx,nslope)) |
---|
| 593 | ps_dev = ps_start0 - ps_avg |
---|
[3571] | 594 | tsurf_dev = tsurf - tsurf_avg |
---|
| 595 | tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:) |
---|
[2794] | 596 | |
---|
[2835] | 597 | !------------------------ |
---|
[3028] | 598 | ! I Initialization |
---|
| 599 | ! I_e Initialization of the PEM variables and soil |
---|
[2835] | 600 | !------------------------ |
---|
[3028] | 601 | call end_comsoil_h_PEM |
---|
| 602 | call ini_comsoil_h_PEM(ngrid,nslope) |
---|
| 603 | call end_adsorption_h_PEM |
---|
| 604 | call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) |
---|
[3493] | 605 | call end_ice_table |
---|
| 606 | call ini_ice_table(ngrid,nslope,nsoilmx_PEM) |
---|
[2794] | 607 | |
---|
[3599] | 608 | allocate(tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen),watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) |
---|
[3028] | 609 | if (soil_pem) then |
---|
| 610 | call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) |
---|
[3367] | 611 | tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg |
---|
[3149] | 612 | watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries |
---|
[3571] | 613 | tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries |
---|
[3070] | 614 | do l = nsoilmx + 1,nsoilmx_PEM |
---|
[3367] | 615 | tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:) |
---|
[3070] | 616 | watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:) |
---|
[3571] | 617 | tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:) |
---|
[3028] | 618 | enddo |
---|
[3367] | 619 | watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen |
---|
[3028] | 620 | endif !soil_pem |
---|
[3571] | 621 | deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries) |
---|
[2794] | 622 | |
---|
[2779] | 623 | !------------------------ |
---|
[3028] | 624 | ! I Initialization |
---|
[3143] | 625 | ! I_f Compute tendencies |
---|
[3028] | 626 | !------------------------ |
---|
[3591] | 627 | allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope),d_co2ice_ini(ngrid,nslope)) |
---|
[3498] | 628 | call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice) |
---|
| 629 | call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice) |
---|
| 630 | d_co2ice_ini = d_co2ice |
---|
[3365] | 631 | deallocate(min_co2_ice,min_h2o_ice) |
---|
[2895] | 632 | |
---|
[2835] | 633 | !------------------------ |
---|
[3028] | 634 | ! I Initialization |
---|
[3571] | 635 | ! I_g Compute global surface pressure |
---|
[3028] | 636 | !------------------------ |
---|
[3571] | 637 | total_surface = sum(cell_area) |
---|
| 638 | ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface |
---|
| 639 | ps_avg_global_new = ps_avg_global_ini |
---|
[3673] | 640 | write(*,*) "Total surface of the planet =", total_surface |
---|
| 641 | write(*,*) "Initial global average pressure =", ps_avg_global_ini |
---|
[3384] | 642 | |
---|
| 643 | !------------------------ |
---|
| 644 | ! I Initialization |
---|
| 645 | ! I_h Read the "startpem.nc" |
---|
| 646 | !------------------------ |
---|
[3603] | 647 | write(*,*) '> Reading "startpem.nc"' |
---|
[3591] | 648 | allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope)) |
---|
| 649 | allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid)) |
---|
[3319] | 650 | if (layering_algo) then |
---|
| 651 | do islope = 1,nslope |
---|
| 652 | do i = 1,ngrid |
---|
| 653 | call ini_layering(stratif(i,islope)) |
---|
| 654 | enddo |
---|
[3297] | 655 | enddo |
---|
[3319] | 656 | endif |
---|
[3297] | 657 | |
---|
[3591] | 658 | delta_h2o_icetablesublim = 0. |
---|
[3602] | 659 | call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, & |
---|
| 660 | tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice, & |
---|
| 661 | watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif) |
---|
[3571] | 662 | deallocate(tsurf_avg_yr1) |
---|
[2779] | 663 | |
---|
[3384] | 664 | ! We save the places where h2o ice is sublimating |
---|
| 665 | ! We compute the surface of h2o ice sublimating |
---|
[3591] | 666 | allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope),co2ice_disappeared(ngrid,nslope)) |
---|
[3571] | 667 | co2ice_sublim_surf_ini = 0. |
---|
[3384] | 668 | h2oice_ini_surf = 0. |
---|
[3571] | 669 | is_co2ice_sublim_ini = .false. |
---|
| 670 | is_h2oice_sublim_ini = .false. |
---|
[3384] | 671 | is_co2ice_ini = .false. |
---|
| 672 | co2ice_disappeared = .false. |
---|
| 673 | do i = 1,ngrid |
---|
| 674 | do islope = 1,nslope |
---|
| 675 | if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true. |
---|
[3498] | 676 | if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then |
---|
[3571] | 677 | is_co2ice_sublim_ini(i,islope) = .true. |
---|
| 678 | co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope) |
---|
[3384] | 679 | endif |
---|
[3498] | 680 | if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then |
---|
[3571] | 681 | is_h2oice_sublim_ini(i,islope) = .true. |
---|
[3384] | 682 | h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) |
---|
| 683 | endif |
---|
| 684 | enddo |
---|
| 685 | enddo |
---|
[3603] | 686 | write(*,*) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini |
---|
| 687 | write(*,*) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf |
---|
[3384] | 688 | |
---|
[3028] | 689 | if (adsorption_pem) then |
---|
[3553] | 690 | totmassco2_adsorbed = 0. |
---|
| 691 | totmassh2o_adsorbed = 0. |
---|
[3028] | 692 | do ig = 1,ngrid |
---|
[3070] | 693 | do islope = 1,nslope |
---|
[3028] | 694 | do l = 1,nsoilmx_PEM - 1 |
---|
[3571] | 695 | if (l == 1) then |
---|
[3553] | 696 | totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
[3028] | 697 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
[3553] | 698 | totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
[3028] | 699 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
[3264] | 700 | else |
---|
[3553] | 701 | totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & |
---|
[3264] | 702 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
[3553] | 703 | totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & |
---|
[3264] | 704 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
| 705 | endif |
---|
[3028] | 706 | enddo |
---|
[2961] | 707 | enddo |
---|
[3028] | 708 | enddo |
---|
[3553] | 709 | write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed |
---|
| 710 | write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed |
---|
[3028] | 711 | endif ! adsorption |
---|
[2794] | 712 | |
---|
[2835] | 713 | !------------------------ |
---|
[3028] | 714 | ! I Initialization |
---|
[2835] | 715 | ! I_i Compute orbit criterion |
---|
[3028] | 716 | !------------------------ |
---|
[2842] | 717 | #ifndef CPP_STD |
---|
[3050] | 718 | call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
---|
[2842] | 719 | #else |
---|
[3050] | 720 | call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) |
---|
[2842] | 721 | #endif |
---|
[2794] | 722 | |
---|
[3498] | 723 | n_myear_leg = Max_iter_pem |
---|
| 724 | if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg) |
---|
[3403] | 725 | |
---|
[3028] | 726 | !-------------------------- END INITIALIZATION ------------------------- |
---|
[2794] | 727 | |
---|
[3028] | 728 | !-------------------------------- RUN ---------------------------------- |
---|
[2794] | 729 | !------------------------ |
---|
| 730 | ! II Run |
---|
[3065] | 731 | ! II_a Update pressure, ice and tracers |
---|
[2794] | 732 | !------------------------ |
---|
[3616] | 733 | write(*,*) |
---|
[3673] | 734 | write(*,*) '********* PEM cycle *********' |
---|
[3498] | 735 | i_myear_leg = 0 |
---|
[3149] | 736 | stopPEM = 0 |
---|
[3319] | 737 | if (layering_algo) then |
---|
| 738 | allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope)) |
---|
| 739 | new_str = .true. |
---|
| 740 | new_lag1 = .true. |
---|
| 741 | new_lag2 = .true. |
---|
| 742 | do islope = 1,nslope |
---|
| 743 | do ig = 1,ngrid |
---|
| 744 | current1(ig,islope)%p => stratif(ig,islope)%top |
---|
| 745 | current2(ig,islope)%p => stratif(ig,islope)%top |
---|
| 746 | enddo |
---|
[3297] | 747 | enddo |
---|
[3319] | 748 | endif |
---|
[2794] | 749 | |
---|
[3498] | 750 | do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear) |
---|
[2835] | 751 | ! II.a.1. Compute updated global pressure |
---|
[3673] | 752 | write(*,'(a,f10.2)') ' **** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1 |
---|
[3603] | 753 | write(*,*) "> Updating the surface pressure" |
---|
[3571] | 754 | ps_avg_global_old = ps_avg_global_new |
---|
[3028] | 755 | do i = 1,ngrid |
---|
| 756 | do islope = 1,nslope |
---|
[3611] | 757 | ps_avg_global_new = ps_avg_global_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface |
---|
[3028] | 758 | enddo |
---|
| 759 | enddo |
---|
| 760 | if (adsorption_pem) then |
---|
| 761 | do i = 1,ngrid |
---|
[3611] | 762 | ps_avg_global_new = ps_avg_global_new - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface |
---|
[3050] | 763 | enddo |
---|
[3028] | 764 | endif |
---|
[3680] | 765 | ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old |
---|
[3603] | 766 | write(*,*) 'Global average pressure old time step:',ps_avg_global_old |
---|
| 767 | write(*,*) 'Global average pressure new time step:',ps_avg_global_new |
---|
[2835] | 768 | |
---|
[3571] | 769 | ! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption) |
---|
[3603] | 770 | write(*,*) "> Updating the surface pressure timeseries for the new pressure" |
---|
[3571] | 771 | allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen)) |
---|
[3028] | 772 | do l = 1,nlayer + 1 |
---|
| 773 | do ig = 1,ngrid |
---|
[3571] | 774 | zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
---|
[3028] | 775 | enddo |
---|
| 776 | enddo |
---|
[3673] | 777 | ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old |
---|
[3603] | 778 | write(*,*) "> Updating the pressure levels timeseries for the new pressure" |
---|
[3571] | 779 | allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen)) |
---|
[3028] | 780 | do l = 1,nlayer + 1 |
---|
| 781 | do ig = 1,ngrid |
---|
[3571] | 782 | zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
---|
[3028] | 783 | enddo |
---|
| 784 | enddo |
---|
[2779] | 785 | |
---|
[3571] | 786 | ! II.a.3. Tracers timeseries |
---|
[3603] | 787 | write(*,*) "> Updating the tracer VMR timeseries for the new pressure" |
---|
[3591] | 788 | allocate(vmr_co2_PEM_phys(ngrid,timelen)) |
---|
[3028] | 789 | l = 1 |
---|
| 790 | do ig = 1,ngrid |
---|
| 791 | do t = 1,timelen |
---|
[3571] | 792 | ! H2O |
---|
| 793 | q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ & |
---|
| 794 | (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) |
---|
[3143] | 795 | if (q_h2o_PEM_phys(ig,t) < 0) then |
---|
| 796 | q_h2o_PEM_phys(ig,t) = 1.e-30 |
---|
| 797 | else if (q_h2o_PEM_phys(ig,t) > 1) then |
---|
| 798 | q_h2o_PEM_phys(ig,t) = 1. |
---|
| 799 | endif |
---|
[3571] | 800 | ! CO2 |
---|
| 801 | q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ & |
---|
| 802 | (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) & |
---|
| 803 | + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) & |
---|
| 804 | - (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ & |
---|
| 805 | (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) |
---|
[3028] | 806 | if (q_co2_PEM_phys(ig,t) < 0) then |
---|
| 807 | q_co2_PEM_phys(ig,t) = 1.e-30 |
---|
[3143] | 808 | else if (q_co2_PEM_phys(ig,t) > 1) then |
---|
[3028] | 809 | q_co2_PEM_phys(ig,t) = 1. |
---|
| 810 | endif |
---|
[3571] | 811 | mmean = 1./(A*q_co2_PEM_phys(ig,t) + B) |
---|
[3149] | 812 | vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 |
---|
[3028] | 813 | enddo |
---|
| 814 | enddo |
---|
[3571] | 815 | deallocate(zplev_timeseries_new,zplev_timeseries_old) |
---|
[2794] | 816 | |
---|
[3028] | 817 | !------------------------ |
---|
[2835] | 818 | ! II Run |
---|
[3149] | 819 | ! II_b Evolution of ice |
---|
[3028] | 820 | !------------------------ |
---|
[3591] | 821 | allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) |
---|
[3553] | 822 | call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) |
---|
| 823 | call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) |
---|
[3319] | 824 | if (layering_algo) then |
---|
| 825 | do islope = 1,nslope |
---|
| 826 | do ig = 1,ngrid |
---|
[3553] | 827 | call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p) |
---|
[3319] | 828 | enddo |
---|
[3297] | 829 | enddo |
---|
[3553] | 830 | else |
---|
| 831 | zlag = 0. |
---|
[3319] | 832 | endif |
---|
[2794] | 833 | |
---|
| 834 | !------------------------ |
---|
| 835 | ! II Run |
---|
[3149] | 836 | ! II_c Flow of glaciers |
---|
[2794] | 837 | !------------------------ |
---|
[3591] | 838 | allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope)) |
---|
[3571] | 839 | if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, & |
---|
[3591] | 840 | ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow) |
---|
| 841 | if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow) |
---|
[3065] | 842 | |
---|
[2794] | 843 | !------------------------ |
---|
| 844 | ! II Run |
---|
[2835] | 845 | ! II_d Update surface and soil temperatures |
---|
[2794] | 846 | !------------------------ |
---|
[2835] | 847 | ! II_d.1 Update Tsurf |
---|
[3603] | 848 | write(*,*) "> Updating surface temperature" |
---|
[3028] | 849 | do ig = 1,ngrid |
---|
| 850 | do islope = 1,nslope |
---|
[3571] | 851 | ! CO2 ice disappeared so we look for the closest point without CO2 ice |
---|
| 852 | if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then |
---|
[3331] | 853 | co2ice_disappeared(ig,islope) = .true. |
---|
[3616] | 854 | if (latitude_deg(ig) > 0.) then ! North hemisphere |
---|
| 855 | outer1: do ig_loop = ig,ngrid ! Go towards equator |
---|
| 856 | do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be equator-ward) |
---|
[3330] | 857 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
[3367] | 858 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
[3571] | 859 | exit outer1 |
---|
[3028] | 860 | endif |
---|
| 861 | enddo |
---|
[3616] | 862 | do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be pole-ward) |
---|
| 863 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
| 864 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
| 865 | exit outer1 |
---|
| 866 | endif |
---|
| 867 | enddo |
---|
[3571] | 868 | enddo outer1 |
---|
[3616] | 869 | else ! South hemisphere |
---|
| 870 | outer2: do ig_loop = ig,1,-1 ! Go towards equator |
---|
| 871 | do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be equator-ward) |
---|
[3330] | 872 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
[3367] | 873 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
[3571] | 874 | exit outer2 |
---|
[3028] | 875 | endif |
---|
| 876 | enddo |
---|
[3616] | 877 | do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be pole-ward) |
---|
| 878 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
| 879 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
| 880 | exit outer2 |
---|
| 881 | endif |
---|
| 882 | enddo |
---|
[3571] | 883 | enddo outer2 |
---|
[2835] | 884 | endif |
---|
[3616] | 885 | else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as Tcond CO2 |
---|
[3571] | 886 | call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg) |
---|
[2835] | 887 | endif |
---|
| 888 | enddo |
---|
[3028] | 889 | enddo |
---|
[2794] | 890 | |
---|
[3028] | 891 | if (soil_pem) then |
---|
[3553] | 892 | ! II_d.2 Shifting soil temperature to surface |
---|
[3673] | 893 | !call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) |
---|
[2794] | 894 | |
---|
[3553] | 895 | ! II_d.3 Update soil temperature |
---|
[3603] | 896 | write(*,*)"> Updating soil temperature profile" |
---|
[3610] | 897 | allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen)) |
---|
| 898 | tsoil_PEM_timeseries_old = tsoil_PEM_timeseries |
---|
[3426] | 899 | do islope = 1,nslope |
---|
[3571] | 900 | tsoil_avg_old = tsoil_PEM(:,:,islope) |
---|
[3426] | 901 | call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) |
---|
| 902 | call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) |
---|
[2794] | 903 | |
---|
[3028] | 904 | do t = 1,timelen |
---|
| 905 | do ig = 1,ngrid |
---|
| 906 | do isoil = 1,nsoilmx_PEM |
---|
[3609] | 907 | ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries |
---|
| 908 | tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil) |
---|
| 909 | ! Update of watersoil density |
---|
[3571] | 910 | watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)*mmol(igcm_h2o_vap)/(mugaz*r) |
---|
[3426] | 911 | if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1) |
---|
[3028] | 912 | enddo |
---|
| 913 | enddo |
---|
| 914 | enddo |
---|
| 915 | enddo |
---|
[3367] | 916 | watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen |
---|
[3571] | 917 | deallocate(tsoil_avg_old) |
---|
[2888] | 918 | |
---|
[3553] | 919 | ! II_d.4 Update the ice table |
---|
[3610] | 920 | allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope)) |
---|
[3170] | 921 | if (icetable_equilibrium) then |
---|
[3603] | 922 | write(*,*) "> Updating ice table (equilibrium method)" |
---|
[3525] | 923 | icetable_thickness_old = icetable_thickness |
---|
[3493] | 924 | call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) |
---|
[3525] | 925 | call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere |
---|
[3490] | 926 | else if (icetable_dynamic) then |
---|
[3603] | 927 | write(*,*) "> Updating ice table (dynamic method)" |
---|
[3525] | 928 | ice_porefilling_old = ice_porefilling |
---|
[3610] | 929 | icetable_depth_old = icetable_depth |
---|
[3493] | 930 | allocate(porefill(nsoilmx_PEM)) |
---|
| 931 | do ig = 1,ngrid |
---|
| 932 | do islope = 1,nslope |
---|
[3571] | 933 | call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps_avg(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth) |
---|
[3493] | 934 | icetable_depth(ig,islope) = ssi_depth |
---|
| 935 | ice_porefilling(ig,:,islope) = porefill |
---|
| 936 | enddo |
---|
| 937 | enddo |
---|
| 938 | deallocate(porefill) |
---|
[3673] | 939 | call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere |
---|
[3486] | 940 | endif |
---|
[3591] | 941 | deallocate(icetable_thickness_old,ice_porefilling_old) |
---|
[3486] | 942 | |
---|
[3553] | 943 | ! II_d.5 Update the soil thermal properties |
---|
[3571] | 944 | call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,ps_avg_global_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
[2794] | 945 | |
---|
[3553] | 946 | ! II_d.6 Update the mass of the regolith adsorbed |
---|
[3028] | 947 | if (adsorption_pem) then |
---|
[3571] | 948 | call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, & |
---|
| 949 | tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & |
---|
[3553] | 950 | h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed) |
---|
[2794] | 951 | |
---|
[3553] | 952 | totmassco2_adsorbed = 0. |
---|
| 953 | totmassh2o_adsorbed = 0. |
---|
[3028] | 954 | do ig = 1,ngrid |
---|
[3493] | 955 | do islope = 1,nslope |
---|
[3532] | 956 | do l = 1,nsoilmx_PEM |
---|
[3493] | 957 | if (l == 1) then |
---|
[3553] | 958 | totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
[3264] | 959 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
[3553] | 960 | totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
[3264] | 961 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
| 962 | else |
---|
[3553] | 963 | totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & |
---|
[3264] | 964 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
[3553] | 965 | totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & |
---|
[3264] | 966 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
| 967 | endif |
---|
[3028] | 968 | enddo |
---|
| 969 | enddo |
---|
| 970 | enddo |
---|
[3603] | 971 | write(*,*) "Total mass of CO2 in the regolith =", totmassco2_adsorbed |
---|
| 972 | write(*,*) "Total mass of H2O in the regolith =", totmassh2o_adsorbed |
---|
[3028] | 973 | endif |
---|
| 974 | endif !soil_pem |
---|
[3666] | 975 | deallocate(zshift_surf,zlag) |
---|
[3028] | 976 | |
---|
[2794] | 977 | !------------------------ |
---|
| 978 | ! II Run |
---|
[3088] | 979 | ! II_e Outputs |
---|
[2794] | 980 | !------------------------ |
---|
[3571] | 981 | call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/)) |
---|
[3088] | 982 | do islope = 1,nslope |
---|
| 983 | write(str2(1:2),'(i2.2)') islope |
---|
[3181] | 984 | call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope)) |
---|
| 985 | call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope)) |
---|
[3498] | 986 | call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope)) |
---|
| 987 | call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope)) |
---|
[3591] | 988 | call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope))) |
---|
| 989 | call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope))) |
---|
[3571] | 990 | call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope)) |
---|
[3339] | 991 | if (icetable_equilibrium) then |
---|
[3493] | 992 | call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) |
---|
[3537] | 993 | call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope)) |
---|
[3490] | 994 | else if (icetable_dynamic) then |
---|
[3493] | 995 | call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) |
---|
[3486] | 996 | endif |
---|
| 997 | |
---|
[3339] | 998 | if (soil_pem) then |
---|
[3571] | 999 | call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope)) |
---|
| 1000 | call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope)) |
---|
[3543] | 1001 | if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope)) |
---|
[3171] | 1002 | if (adsorption_pem) then |
---|
[3553] | 1003 | call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope)) |
---|
| 1004 | call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope)) |
---|
[3532] | 1005 | endif |
---|
[3171] | 1006 | endif |
---|
[3088] | 1007 | enddo |
---|
[3591] | 1008 | deallocate(flag_co2flow,flag_h2oflow) |
---|
[3088] | 1009 | |
---|
| 1010 | !------------------------ |
---|
| 1011 | ! II Run |
---|
| 1012 | ! II_f Update the tendencies |
---|
| 1013 | !------------------------ |
---|
[3571] | 1014 | call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new) |
---|
[3707] | 1015 | write(*,*) "> Updating the H2O sub-surface ice tendency due to lag layer" |
---|
| 1016 | do ig = 1,ngrid |
---|
| 1017 | do islope = 1,nslope |
---|
| 1018 | call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) |
---|
| 1019 | enddo |
---|
| 1020 | enddo |
---|
[3666] | 1021 | if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old) |
---|
| 1022 | deallocate(vmr_co2_PEM_phys) |
---|
[3711] | 1023 | write(*,*) "> Updating the H2O sub-surface ice depth" |
---|
| 1024 | do ig = 1,ngrid |
---|
| 1025 | do islope = 1,nslope |
---|
| 1026 | if (icetable_depth(ig,islope) > 0) then |
---|
| 1027 | icetable_depth(ig,islope)=icetable_depth(ig,islope) + d_h2oice(ig,islope)*0.01 !!! 0.01 is for 1 percent dust, needs to be updated !!! |
---|
| 1028 | endif |
---|
| 1029 | enddo |
---|
| 1030 | enddo |
---|
[2794] | 1031 | |
---|
[3711] | 1032 | |
---|
[2835] | 1033 | !------------------------ |
---|
| 1034 | ! II Run |
---|
[3088] | 1035 | ! II_g Checking the stopping criterion |
---|
[2835] | 1036 | !------------------------ |
---|
[3603] | 1037 | write(*,*) "> Checking the stopping criteria" |
---|
[3711] | 1038 | h2o_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) |
---|
| 1039 | co2_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) |
---|
[3571] | 1040 | call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid) |
---|
| 1041 | call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope) |
---|
[3498] | 1042 | i_myear_leg = i_myear_leg + dt |
---|
| 1043 | i_myear = i_myear + dt |
---|
| 1044 | if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5 |
---|
[3430] | 1045 | if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6 |
---|
[3389] | 1046 | call system_clock(c2) |
---|
[3430] | 1047 | if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7 |
---|
[3432] | 1048 | if (stopPEM > 0) then |
---|
[3149] | 1049 | select case (stopPEM) |
---|
| 1050 | case(1) |
---|
[3673] | 1051 | write(*,'(a,i0,a)') " **** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above." |
---|
[3149] | 1052 | case(2) |
---|
[3673] | 1053 | write(*,'(a,i0,a)') " **** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above." |
---|
[3149] | 1054 | case(3) |
---|
[3673] | 1055 | write(*,'(a,i0,a)') " **** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above." |
---|
[3149] | 1056 | case(4) |
---|
[3673] | 1057 | write(*,'(a,i0,a)') " **** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above." |
---|
[3149] | 1058 | case(5) |
---|
[3673] | 1059 | write(*,'(a,i0)') " **** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM |
---|
[3149] | 1060 | case(6) |
---|
[3673] | 1061 | write(*,'(a,i0)') " **** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM |
---|
[3363] | 1062 | case(7) |
---|
[3673] | 1063 | write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM |
---|
[3430] | 1064 | case(8) |
---|
[3673] | 1065 | write(*,'(a,i0)') " **** STOPPING because the layering algorithm met an hasty end: ", stopPEM |
---|
[3149] | 1066 | case default |
---|
[3673] | 1067 | write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM |
---|
[3149] | 1068 | end select |
---|
[2779] | 1069 | exit |
---|
[3028] | 1070 | else |
---|
[3673] | 1071 | write(*,'(a,f10.2,a)') ' **** The chained simulation has run for ',i_myear,' Martian years.' |
---|
| 1072 | write(*,*) '**** The PEM can continue!' |
---|
| 1073 | write(*,*) '****' |
---|
[3028] | 1074 | endif |
---|
[3149] | 1075 | enddo ! big time iteration loop of the pem |
---|
[3591] | 1076 | deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed) |
---|
[3598] | 1077 | deallocate(watersoil_density_PEM_avg,watersurf_density_avg) |
---|
[3591] | 1078 | deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries) |
---|
| 1079 | deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim) |
---|
| 1080 | deallocate(d_co2ice,d_co2ice_ini,d_h2oice) |
---|
| 1081 | deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini) |
---|
[3666] | 1082 | if (layering_algo) deallocate(new_str,new_lag1,new_lag2,current1,current2) |
---|
[3028] | 1083 | !------------------------------ END RUN -------------------------------- |
---|
[2779] | 1084 | |
---|
[3028] | 1085 | !------------------------------- OUTPUT -------------------------------- |
---|
[2794] | 1086 | !------------------------ |
---|
| 1087 | ! III Output |
---|
[2835] | 1088 | ! III_a Update surface value for the PCM start files |
---|
[2794] | 1089 | !------------------------ |
---|
[3571] | 1090 | ! III_a.1 Ice update for start file |
---|
[3616] | 1091 | write(*,*) |
---|
[3673] | 1092 | write(*,*) '********* PEM finalization *********' |
---|
[3603] | 1093 | write(*,*) '> Reconstructing perennial ice and frost for the PCM' |
---|
[3149] | 1094 | watercap = 0. |
---|
[3705] | 1095 | perennial_co2ice =sum(stratif(ig,islope)%top%co2ice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) |
---|
[3028] | 1096 | do ig = 1,ngrid |
---|
[3159] | 1097 | ! H2O ice metamorphism |
---|
[3705] | 1098 | !if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then |
---|
| 1099 | ! h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold |
---|
| 1100 | ! qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold |
---|
| 1101 | !endif |
---|
[3159] | 1102 | |
---|
| 1103 | ! Is H2O ice still considered as an infinite reservoir for the PCM? |
---|
[3705] | 1104 | !if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then |
---|
| 1105 | if (sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then |
---|
| 1106 | ! There is enough ice to be considered as an infinite reservoir |
---|
[3149] | 1107 | watercaptag(ig) = .true. |
---|
| 1108 | else |
---|
[3159] | 1109 | ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost |
---|
[3149] | 1110 | watercaptag(ig) = .false. |
---|
[3705] | 1111 | |
---|
| 1112 | ! qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:) |
---|
| 1113 | |
---|
| 1114 | qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) |
---|
[3149] | 1115 | h2o_ice(ig,:) = 0. |
---|
[3028] | 1116 | endif |
---|
[3159] | 1117 | |
---|
| 1118 | ! CO2 ice metamorphism |
---|
[3705] | 1119 | !if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then |
---|
| 1120 | ! perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold |
---|
| 1121 | ! qsurf(ig,igcm_co2,:) = metam_co2ice_threshold |
---|
| 1122 | !endif |
---|
[3028] | 1123 | enddo |
---|
[2888] | 1124 | |
---|
[3571] | 1125 | ! III.a.2. Tsurf update for start file |
---|
[3603] | 1126 | write(*,*) '> Reconstructing the surface temperature for the PCM' |
---|
[3571] | 1127 | tsurf = tsurf_avg + tsurf_dev |
---|
[3591] | 1128 | deallocate(tsurf_dev) |
---|
[3571] | 1129 | |
---|
| 1130 | ! III_a.3 Tsoil update for start file |
---|
[3028] | 1131 | if (soil_pem) then |
---|
[3603] | 1132 | write(*,*) '> Reconstructing the soil temperature profile for the PCM' |
---|
[3532] | 1133 | inertiesoil = TI_PEM(:,:nsoilmx,:) |
---|
[3571] | 1134 | ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value |
---|
| 1135 | do isoil = 1,nsoilmx |
---|
| 1136 | tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:) |
---|
| 1137 | enddo |
---|
| 1138 | tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev |
---|
[3172] | 1139 | #ifndef CPP_STD |
---|
[3181] | 1140 | flux_geo = fluxgeo |
---|
[3172] | 1141 | #endif |
---|
[3028] | 1142 | endif |
---|
[3591] | 1143 | deallocate(tsurf_avg,tsoil_dev) |
---|
[2779] | 1144 | |
---|
[3571] | 1145 | ! III_a.4 Pressure update for start file |
---|
[3603] | 1146 | write(*,*) '> Reconstructing the pressure for the PCM' |
---|
[3591] | 1147 | allocate(ps_start(ngrid)) |
---|
[3680] | 1148 | ! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop |
---|
| 1149 | ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini |
---|
[3591] | 1150 | deallocate(ps_avg,ps_dev) |
---|
[2794] | 1151 | |
---|
[3571] | 1152 | ! III_a.5 Tracers update for start file |
---|
[3603] | 1153 | write(*,*) '> Reconstructing the tracer VMR for the PCM' |
---|
[3571] | 1154 | allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1)) |
---|
[3028] | 1155 | do l = 1,nlayer + 1 |
---|
[3571] | 1156 | zplev_start0(:,l) = ap(l) + bp(l)*ps_start0 |
---|
| 1157 | zplev_new(:,l) = ap(l) + bp(l)*ps_start |
---|
[3028] | 1158 | enddo |
---|
[2835] | 1159 | |
---|
[3028] | 1160 | do nnq = 1,nqtot |
---|
| 1161 | if (noms(nnq) /= "co2") then |
---|
| 1162 | do l = 1,llm - 1 |
---|
| 1163 | do ig = 1,ngrid |
---|
[3571] | 1164 | q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
[3028] | 1165 | enddo |
---|
| 1166 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
---|
| 1167 | enddo |
---|
| 1168 | else |
---|
| 1169 | do l = 1,llm - 1 |
---|
| 1170 | do ig = 1,ngrid |
---|
[3571] | 1171 | q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & |
---|
| 1172 | + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_start0(ig,l) - zplev_start0(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
[3028] | 1173 | enddo |
---|
| 1174 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
---|
| 1175 | enddo |
---|
| 1176 | endif |
---|
| 1177 | enddo |
---|
[3591] | 1178 | deallocate(zplev_start0) |
---|
[2835] | 1179 | |
---|
[3571] | 1180 | ! Conserving the tracers mass for start file |
---|
[3028] | 1181 | do nnq = 1,nqtot |
---|
| 1182 | do ig = 1,ngrid |
---|
| 1183 | do l = 1,llm - 1 |
---|
| 1184 | 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 |
---|
[3065] | 1185 | extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
| 1186 | q(ig,l,nnq) = 1. |
---|
| 1187 | q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) |
---|
[3028] | 1188 | write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" |
---|
[2835] | 1189 | endif |
---|
[3028] | 1190 | if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 |
---|
| 1191 | enddo |
---|
| 1192 | enddo |
---|
| 1193 | enddo |
---|
[3591] | 1194 | deallocate(zplev_new) |
---|
[2779] | 1195 | |
---|
[3584] | 1196 | ! III_a.6 Albedo update for start file |
---|
[3603] | 1197 | write(*,*) '> Reconstructing the albedo for the PCM' |
---|
[3584] | 1198 | do ig = 1,ngrid |
---|
| 1199 | if (latitude(ig) < 0.) then |
---|
| 1200 | icap = 2 ! Southern hemisphere |
---|
| 1201 | else |
---|
| 1202 | icap = 1 ! Northern hemisphere |
---|
| 1203 | endif |
---|
[3603] | 1204 | do islope = 1,nslope |
---|
[3584] | 1205 | ! Bare ground |
---|
| 1206 | albedo(ig,:,islope) = albedodat(ig) |
---|
| 1207 | emis(ig,islope) = emissiv |
---|
| 1208 | |
---|
| 1209 | ! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant |
---|
| 1210 | ! H2O ice |
---|
| 1211 | if (h2o_ice(ig,islope) > 0.) then |
---|
| 1212 | albedo(ig,:,islope) = albedo_h2o_cap |
---|
| 1213 | emis(ig,islope) = 1. |
---|
| 1214 | endif |
---|
| 1215 | ! CO2 ice |
---|
| 1216 | if (co2_ice(ig,islope) > 0.) then |
---|
| 1217 | albedo(ig,:,islope) = albedo_perennialco2(icap) |
---|
| 1218 | emis(ig,islope) = emisice(icap) |
---|
| 1219 | endif |
---|
| 1220 | ! H2O frost |
---|
| 1221 | if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then |
---|
| 1222 | albedo(ig,:,islope) = albedo_h2o_frost |
---|
| 1223 | emis(ig,islope) = 1. |
---|
| 1224 | endif |
---|
| 1225 | ! CO2 frost |
---|
| 1226 | if (qsurf(ig,igcm_co2,islope) > 0.) then |
---|
| 1227 | albedo(ig,:,islope) = albedice(icap) |
---|
| 1228 | emis(ig,islope) = emisice(icap) |
---|
| 1229 | endif |
---|
| 1230 | enddo |
---|
| 1231 | enddo |
---|
| 1232 | |
---|
| 1233 | ! III_a.7 Orbital parameters update for start file |
---|
[3603] | 1234 | write(*,*) '> Setting the new orbital parameters' |
---|
[3498] | 1235 | if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg) |
---|
[2779] | 1236 | |
---|
| 1237 | !------------------------ |
---|
[3028] | 1238 | ! III Output |
---|
[3317] | 1239 | ! III_b Write "restart.nc" and "restartfi.nc" |
---|
[3028] | 1240 | !------------------------ |
---|
[3317] | 1241 | ! III_b.1 Write "restart.nc" |
---|
[3042] | 1242 | ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys |
---|
[3028] | 1243 | pday = day_ini |
---|
[3042] | 1244 | ztime_fin = time_phys |
---|
[2980] | 1245 | #ifndef CPP_1D |
---|
[3603] | 1246 | write(*,*) '> Writing "restart.nc"' |
---|
[3442] | 1247 | ! Correction on teta due to surface pressure changes |
---|
[3603] | 1248 | allocate(pdyn(ip1jmp1)) |
---|
| 1249 | call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn) |
---|
| 1250 | do i = 1,ip1jmp1 |
---|
| 1251 | teta(i,:) = teta(i,:)*pdyn(i)**kappa |
---|
[3442] | 1252 | enddo |
---|
| 1253 | ! Correction on atmospheric pressure |
---|
[3603] | 1254 | allocate(p(ip1jmp1,nlayer + 1)) |
---|
| 1255 | call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn) |
---|
| 1256 | call pression(ip1jmp1,ap,bp,pdyn,p) |
---|
[3442] | 1257 | ! Correction on the mass of atmosphere |
---|
[3028] | 1258 | call massdair(p,masse) |
---|
[3317] | 1259 | call dynredem0("restart.nc",day_ini,phis) |
---|
[3628] | 1260 | call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn) |
---|
[3603] | 1261 | deallocate(ap,bp,p,pdyn) |
---|
[2980] | 1262 | #else |
---|
[3603] | 1263 | write(*,*) '> Writing "restart1D.txt"' |
---|
[3571] | 1264 | call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q) |
---|
[2980] | 1265 | #endif |
---|
[3591] | 1266 | deallocate(ps_start0,ps_start) |
---|
[2980] | 1267 | |
---|
[3317] | 1268 | ! III_b.2 Write the "restartfi.nc" |
---|
[3603] | 1269 | write(*,*) '> Writing "restartfi.nc"' |
---|
[2842] | 1270 | #ifndef CPP_STD |
---|
[3317] | 1271 | call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & |
---|
[3327] | 1272 | nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & |
---|
[3028] | 1273 | inertiedat,def_slope,subslope_dist) |
---|
[3327] | 1274 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & |
---|
[3114] | 1275 | ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & |
---|
| 1276 | albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, & |
---|
[3130] | 1277 | wstar,watercap,perennial_co2ice) |
---|
[2842] | 1278 | #else |
---|
[3591] | 1279 | if (allocated(noms)) deallocate(noms) |
---|
| 1280 | deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil) |
---|
[3317] | 1281 | call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & |
---|
[3327] | 1282 | nlayer,nq,ptimestep,pday,time_phys,cell_area, & |
---|
[3028] | 1283 | albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) |
---|
[3327] | 1284 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & |
---|
[3114] | 1285 | ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, & |
---|
| 1286 | cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, & |
---|
| 1287 | tsea_ice,sea_ice) |
---|
[2842] | 1288 | #endif |
---|
| 1289 | |
---|
[2794] | 1290 | !------------------------ |
---|
| 1291 | ! III Output |
---|
[3161] | 1292 | ! III_c Write the "restartpem.nc" |
---|
[2794] | 1293 | !------------------------ |
---|
[3603] | 1294 | write(*,*) '> Writing "restartpem.nc"' |
---|
[3319] | 1295 | if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings) |
---|
[3206] | 1296 | call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) |
---|
[3537] | 1297 | call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & |
---|
[3553] | 1298 | co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif) |
---|
[2779] | 1299 | |
---|
[3498] | 1300 | call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear) |
---|
[3149] | 1301 | |
---|
[3616] | 1302 | write(*,*) |
---|
[3738] | 1303 | write(*,*) '****** PEM final information *******' |
---|
[3616] | 1304 | write(*,'(a,f10.2,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years." |
---|
[3638] | 1305 | write(*,'(a,f10.2,a,f10.2,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years." |
---|
[3616] | 1306 | write(*,'(a,f10.2,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years." |
---|
| 1307 | write(*,*) "+ PEM: so far, so good!" |
---|
[3673] | 1308 | write(*,*) '************************************' |
---|
[2794] | 1309 | |
---|
[3319] | 1310 | if (layering_algo) then |
---|
| 1311 | do islope = 1,nslope |
---|
| 1312 | do i = 1,ngrid |
---|
| 1313 | call del_layering(stratif(i,islope)) |
---|
| 1314 | enddo |
---|
[3297] | 1315 | enddo |
---|
[3319] | 1316 | endif |
---|
[3591] | 1317 | deallocate(q,longitude,latitude,cell_area,tsoil_PEM) |
---|
| 1318 | deallocate(co2_ice,h2o_ice,stratif) |
---|
[3028] | 1319 | !----------------------------- END OUTPUT ------------------------------ |
---|
[2897] | 1320 | |
---|
[2779] | 1321 | END PROGRAM pem |
---|