| 1 | !------------------------ |
|---|
| 2 | ! I Initialization |
|---|
| 3 | ! I_a READ run.def |
|---|
| 4 | ! I_b READ of start_evol.nc and starfi_evol.nc |
|---|
| 5 | ! I_c Subslope parametrisation |
|---|
| 6 | ! I_d READ GCM data and convert to the physical grid |
|---|
| 7 | ! I_e Initialization of the PEM variable and soil |
|---|
| 8 | ! I_f Compute tendencies & Save initial situation |
|---|
| 9 | ! I_g Save initial PCM situation |
|---|
| 10 | ! I_h Read the PEMstart |
|---|
| 11 | ! I_i Compute orbit criterion |
|---|
| 12 | |
|---|
| 13 | ! II Run |
|---|
| 14 | ! II_a Update pressure, ice and tracers |
|---|
| 15 | ! II_b Evolution of the ice |
|---|
| 16 | ! II_c CO2 & H2O glaciers flows |
|---|
| 17 | ! II_d Update surface and soil temperatures |
|---|
| 18 | ! II_e Update the tendencies |
|---|
| 19 | ! II_f Checking the stopping criterion |
|---|
| 20 | |
|---|
| 21 | ! III Output |
|---|
| 22 | ! III_a Update surface value for the PCM start files |
|---|
| 23 | ! III_b Write restart_evol.nc and restartfi_evol.nc |
|---|
| 24 | ! III_c Write restartfi_PEM.nc |
|---|
| 25 | !------------------------ |
|---|
| 26 | |
|---|
| 27 | PROGRAM pem |
|---|
| 28 | |
|---|
| 29 | use phyetat0_mod, only: phyetat0 |
|---|
| 30 | use phyredem, only: physdem0, physdem1 |
|---|
| 31 | use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close |
|---|
| 32 | use turb_mod, only: q2, wstar |
|---|
| 33 | use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h |
|---|
| 34 | use logic_mod, only: iflag_phys |
|---|
| 35 | use mod_const_mpi, only: COMM_LMDZ |
|---|
| 36 | use comconst_mod, only: rad, g, cpp, pi |
|---|
| 37 | use infotrac |
|---|
| 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: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow |
|---|
| 42 | use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop |
|---|
| 43 | use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, & |
|---|
| 44 | m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial |
|---|
| 45 | use evol_co2_ice_s_mod, only: evol_co2_ice_s |
|---|
| 46 | use evol_h2o_ice_s_mod, only: evol_h2o_ice_s |
|---|
| 47 | use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, & |
|---|
| 48 | TI_PEM, inertiedat_PEM, & ! soil thermal inertia |
|---|
| 49 | tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths |
|---|
| 50 | fluxgeo, & ! Geothermal flux for the PEM and GCM |
|---|
| 51 | water_reservoir ! Water ressources |
|---|
| 52 | use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine |
|---|
| 53 | ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays |
|---|
| 54 | co2_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded |
|---|
| 55 | use time_evol_mod, only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini |
|---|
| 56 | use orbit_param_criterion_mod, only: orbit_param_criterion |
|---|
| 57 | use recomp_orb_param_mod, only: recomp_orb_param |
|---|
| 58 | use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, & |
|---|
| 59 | ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi |
|---|
| 60 | use soil_thermalproperties_mod, only: update_soil_thermalproperties |
|---|
| 61 | use time_phylmdz_mod, only: daysec, dtphys, day_end |
|---|
| 62 | |
|---|
| 63 | #ifndef CPP_STD |
|---|
| 64 | use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil |
|---|
| 65 | use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, & |
|---|
| 66 | albedodat, zmea, zstd, zsig, zgam, zthe, & |
|---|
| 67 | hmons, summit, base,albedo_h2o_frost, & |
|---|
| 68 | frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice, & |
|---|
| 69 | emisice, albedice |
|---|
| 70 | use dimradmars_mod, only: totcloudfrac, albedo |
|---|
| 71 | use dust_param_mod, only: tauscaling |
|---|
| 72 | use tracer_mod, only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses |
|---|
| 73 | use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master |
|---|
| 74 | use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit |
|---|
| 75 | use comcstfi_h, only: r, mugaz |
|---|
| 76 | use paleoclimate_mod, only: albedo_perenialco2 |
|---|
| 77 | #else |
|---|
| 78 | use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names |
|---|
| 79 | use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & |
|---|
| 80 | PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, & |
|---|
| 81 | ALBEDO_CO2_ICE_SPECTV, phys_state_var_init |
|---|
| 82 | use aerosol_mod, only: iniaerosol |
|---|
| 83 | use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit |
|---|
| 84 | use comcstfi_mod, only: r, mugaz |
|---|
| 85 | #endif |
|---|
| 86 | |
|---|
| 87 | #ifndef CPP_1D |
|---|
| 88 | use iniphysiq_mod, only: iniphysiq |
|---|
| 89 | use control_mod, only: iphysiq, day_step, nsplit_phys |
|---|
| 90 | #else |
|---|
| 91 | use time_phylmdz_mod, only: iphysiq, day_step |
|---|
| 92 | #endif |
|---|
| 93 | |
|---|
| 94 | #ifdef CPP_1D |
|---|
| 95 | use regular_lonlat_mod, only: init_regular_lonlat |
|---|
| 96 | use physics_distribution_mod, only: init_physics_distribution |
|---|
| 97 | use mod_grid_phy_lmdz, only: regular_lonlat |
|---|
| 98 | use init_phys_1d_mod, only: init_phys_1d |
|---|
| 99 | #endif |
|---|
| 100 | |
|---|
| 101 | IMPLICIT NONE |
|---|
| 102 | |
|---|
| 103 | include "dimensions.h" |
|---|
| 104 | include "paramet.h" |
|---|
| 105 | include "comgeom.h" |
|---|
| 106 | include "iniprint.h" |
|---|
| 107 | include "callkeys.h" |
|---|
| 108 | |
|---|
| 109 | integer ngridmx |
|---|
| 110 | parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) |
|---|
| 111 | |
|---|
| 112 | ! Same variable names as in the GCM |
|---|
| 113 | integer :: ngrid ! Number of physical grid points |
|---|
| 114 | integer :: nlayer ! Number of vertical layer |
|---|
| 115 | integer :: nq ! Number of tracer |
|---|
| 116 | integer :: day_ini ! First day of the simulation |
|---|
| 117 | real :: pday ! Physical day |
|---|
| 118 | real :: time_phys ! Same as GCM |
|---|
| 119 | real :: ptimestep ! Same as GCM |
|---|
| 120 | real :: ztime_fin ! Same as GCM |
|---|
| 121 | |
|---|
| 122 | ! Variables to read start.nc |
|---|
| 123 | character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM |
|---|
| 124 | |
|---|
| 125 | ! Dynamic variables |
|---|
| 126 | real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants |
|---|
| 127 | real :: teta(ip1jmp1,llm) ! temperature potentielle |
|---|
| 128 | real, dimension(:,:,:), allocatable :: q ! champs advectes |
|---|
| 129 | real :: ps(ip1jmp1) ! pression au sol |
|---|
| 130 | real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol |
|---|
| 131 | real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées |
|---|
| 132 | real :: masse(ip1jmp1,llm) ! masse d'air |
|---|
| 133 | real :: phis(ip1jmp1) ! geopotentiel au sol |
|---|
| 134 | real :: time_0 |
|---|
| 135 | |
|---|
| 136 | ! Variables to read starfi.nc |
|---|
| 137 | character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM |
|---|
| 138 | character*2 str2 |
|---|
| 139 | integer :: ncid, varid,status ! Variable for handling opening of files |
|---|
| 140 | integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files |
|---|
| 141 | integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files |
|---|
| 142 | integer :: apvarid, bpvarid ! Variable ID for Netcdf files |
|---|
| 143 | |
|---|
| 144 | ! Variables to read starfi.nc and write restartfi.nc |
|---|
| 145 | real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi |
|---|
| 146 | real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi |
|---|
| 147 | real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart |
|---|
| 148 | real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart |
|---|
| 149 | real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi |
|---|
| 150 | real :: Total_surface ! Total surface of the planet |
|---|
| 151 | |
|---|
| 152 | ! Variables for h2o_ice evolution |
|---|
| 153 | real :: ini_surf_h2o ! Initial surface of sublimating h2o ice |
|---|
| 154 | real :: ini_surf_co2 ! Initial surface of sublimating co2 ice |
|---|
| 155 | real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] |
|---|
| 156 | real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step |
|---|
| 157 | real :: global_ave_press_new ! constant: Global average pressure of current time step |
|---|
| 158 | real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] |
|---|
| 159 | real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] |
|---|
| 160 | real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] |
|---|
| 161 | real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step |
|---|
| 162 | logical :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? |
|---|
| 163 | logical :: STOPPING_1_water ! Logical : is there still water ice to sublimate? |
|---|
| 164 | logical :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? |
|---|
| 165 | logical :: STOPPING_pressure ! Logical : is the criterion (% of change in the surface pressure) reached? |
|---|
| 166 | integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param |
|---|
| 167 | real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] |
|---|
| 168 | real, allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] |
|---|
| 169 | real, allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM |
|---|
| 170 | 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] |
|---|
| 171 | real, allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] |
|---|
| 172 | integer :: timelen ! # time samples |
|---|
| 173 | real :: ave ! intermediate varibale to compute average |
|---|
| 174 | real, allocatable :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) |
|---|
| 175 | real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 |
|---|
| 176 | |
|---|
| 177 | ! Variables for slopes |
|---|
| 178 | real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] |
|---|
| 179 | real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] |
|---|
| 180 | real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] |
|---|
| 181 | real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] |
|---|
| 182 | real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] |
|---|
| 183 | real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice |
|---|
| 184 | real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state |
|---|
| 185 | real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state |
|---|
| 186 | real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year |
|---|
| 187 | 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 |
|---|
| 188 | real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial h2o ice |
|---|
| 189 | real, dimension(:,:), allocatable :: flag_co2flow(:,:) ! (ngrid,nslope): Flag where there is a CO2 glacier flow |
|---|
| 190 | real, dimension(:), allocatable :: flag_co2flow_mesh(:) ! (ngrid) : Flag where there is a CO2 glacier flow |
|---|
| 191 | real, dimension(:,:), allocatable :: flag_h2oflow(:,:) ! (ngrid,nslope): Flag where there is a H2O glacier flow |
|---|
| 192 | real, dimension(:), allocatable :: flag_h2oflow_mesh(:) ! (ngrid) : Flag where there is a H2O glacier flow |
|---|
| 193 | |
|---|
| 194 | ! Variables for surface and soil |
|---|
| 195 | real, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] |
|---|
| 196 | real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] |
|---|
| 197 | real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] |
|---|
| 198 | real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] |
|---|
| 199 | real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] |
|---|
| 200 | real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] |
|---|
| 201 | real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] |
|---|
| 202 | real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] |
|---|
| 203 | real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] |
|---|
| 204 | real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] |
|---|
| 205 | real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] |
|---|
| 206 | real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] |
|---|
| 207 | real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] |
|---|
| 208 | real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] |
|---|
| 209 | real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] |
|---|
| 210 | real, allocatable :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] |
|---|
| 211 | real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] |
|---|
| 212 | real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] |
|---|
| 213 | logical :: bool_sublim ! logical to check if there is sublimation or not |
|---|
| 214 | real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] |
|---|
| 215 | real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] |
|---|
| 216 | ! Some variables for the PEM run |
|---|
| 217 | real, parameter :: year_step = 1 ! timestep for the pem |
|---|
| 218 | integer :: year_iter ! number of iteration |
|---|
| 219 | integer :: year_iter_max ! maximum number of iterations before stopping |
|---|
| 220 | integer :: i_myear ! Global number of Martian years of the chained simulations |
|---|
| 221 | integer :: n_myear ! Maximum number of Martian years of the chained simulations |
|---|
| 222 | real :: timestep ! timestep [s] |
|---|
| 223 | real :: watercap_sum ! total mass of water cap [kg/m^2] |
|---|
| 224 | real :: water_sum ! total mass of water in the mesh [kg/m^2] |
|---|
| 225 | |
|---|
| 226 | #ifdef CPP_STD |
|---|
| 227 | real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice |
|---|
| 228 | real :: albedo_h2o_frost ! albedo of h2o frost |
|---|
| 229 | real, allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic |
|---|
| 230 | real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic |
|---|
| 231 | real, allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic |
|---|
| 232 | real, allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic |
|---|
| 233 | real, allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic |
|---|
| 234 | real, allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 235 | real, allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 236 | real, allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 237 | real, allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 238 | real, allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model |
|---|
| 239 | logical, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model |
|---|
| 240 | real, allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 241 | real, allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case |
|---|
| 242 | #endif |
|---|
| 243 | |
|---|
| 244 | #ifdef CPP_1D |
|---|
| 245 | integer :: nsplit_phys |
|---|
| 246 | integer, parameter :: jjm_value = jjm - 1 |
|---|
| 247 | integer :: ierr |
|---|
| 248 | #else |
|---|
| 249 | integer, parameter :: jjm_value = jjm |
|---|
| 250 | #endif |
|---|
| 251 | |
|---|
| 252 | ! Loop variables |
|---|
| 253 | integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap |
|---|
| 254 | |
|---|
| 255 | ! Parallel variables |
|---|
| 256 | #ifndef CPP_STD |
|---|
| 257 | is_sequential = .true. |
|---|
| 258 | is_parallel = .false. |
|---|
| 259 | is_mpi_root = .true. |
|---|
| 260 | is_omp_root = .true. |
|---|
| 261 | is_master = .true. |
|---|
| 262 | #endif |
|---|
| 263 | |
|---|
| 264 | day_ini = 0 ! test |
|---|
| 265 | time_phys = 0. ! test |
|---|
| 266 | |
|---|
| 267 | ! Some constants |
|---|
| 268 | ngrid = ngridmx |
|---|
| 269 | nlayer = llm |
|---|
| 270 | |
|---|
| 271 | A = (1/m_co2 - 1/m_noco2) |
|---|
| 272 | B = 1/m_noco2 |
|---|
| 273 | |
|---|
| 274 | year_day = 669 |
|---|
| 275 | daysec = 88775. |
|---|
| 276 | timestep = year_day*daysec/year_step |
|---|
| 277 | |
|---|
| 278 | !----------------------------- INITIALIZATION -------------------------- |
|---|
| 279 | !------------------------ |
|---|
| 280 | ! I Initialization |
|---|
| 281 | ! I_a READ run.def |
|---|
| 282 | !------------------------ |
|---|
| 283 | #ifndef CPP_1D |
|---|
| 284 | dtphys = 0 |
|---|
| 285 | call conf_gcm(99,.true.) |
|---|
| 286 | call infotrac_init |
|---|
| 287 | nq = nqtot |
|---|
| 288 | allocate(q(ip1jmp1,llm,nqtot)) |
|---|
| 289 | #else |
|---|
| 290 | ! load tracer names from file 'traceur.def' |
|---|
| 291 | open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr) |
|---|
| 292 | if (ierr /= 0) then |
|---|
| 293 | write(*,*) 'Cannot find required file "traceur.def"' |
|---|
| 294 | write(*,*) ' If you want to run with tracers, I need it' |
|---|
| 295 | write(*,*) ' ... might as well stop here ...' |
|---|
| 296 | stop |
|---|
| 297 | else |
|---|
| 298 | write(*,*) "pem1d: Reading file traceur.def" |
|---|
| 299 | ! read number of tracers: |
|---|
| 300 | read(90,*,iostat = ierr) nq |
|---|
| 301 | write(*,*) "nq",nq |
|---|
| 302 | nqtot = nq ! set value of nqtot (in infotrac module) as nq |
|---|
| 303 | if (ierr /= 0) then |
|---|
| 304 | write(*,*) "pem1d: error reading number of tracers" |
|---|
| 305 | write(*,*) " (first line of traceur.def) " |
|---|
| 306 | stop |
|---|
| 307 | endif |
|---|
| 308 | if (nq < 1) then |
|---|
| 309 | write(*,*) "pem1d: error number of tracers" |
|---|
| 310 | write(*,*) "is nq=",nq," but must be >=1!" |
|---|
| 311 | stop |
|---|
| 312 | endif |
|---|
| 313 | endif |
|---|
| 314 | nq = nqtot |
|---|
| 315 | allocate(q(ip1jmp1,llm,nqtot)) |
|---|
| 316 | allocate(ap(nlayer + 1)) |
|---|
| 317 | allocate(bp(nlayer + 1)) |
|---|
| 318 | call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp) |
|---|
| 319 | pi = 2.*asin(1.) |
|---|
| 320 | g = 3.72 |
|---|
| 321 | nsplit_phys = 1 |
|---|
| 322 | #endif |
|---|
| 323 | |
|---|
| 324 | call conf_pem(i_myear,n_myear) |
|---|
| 325 | |
|---|
| 326 | !------------------------ |
|---|
| 327 | ! I Initialization |
|---|
| 328 | ! I_b READ of start_evol.nc and starfi_evol.nc |
|---|
| 329 | !------------------------ |
|---|
| 330 | ! I_b.1 READ start_evol.nc |
|---|
| 331 | allocate(ps_start_GCM(ngrid)) |
|---|
| 332 | #ifndef CPP_1D |
|---|
| 333 | call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0) |
|---|
| 334 | |
|---|
| 335 | call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) |
|---|
| 336 | |
|---|
| 337 | call iniconst !new |
|---|
| 338 | call inigeom |
|---|
| 339 | |
|---|
| 340 | allocate(ap(nlayer + 1)) |
|---|
| 341 | allocate(bp(nlayer + 1)) |
|---|
| 342 | status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid) |
|---|
| 343 | status = nf90_inq_varid(ncid,"ap",apvarid) |
|---|
| 344 | status = nf90_get_var(ncid,apvarid,ap) |
|---|
| 345 | status = nf90_inq_varid(ncid,"bp",bpvarid) |
|---|
| 346 | status = nf90_get_var(ncid,bpvarid,bp) |
|---|
| 347 | status = nf90_close(ncid) |
|---|
| 348 | |
|---|
| 349 | call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, & |
|---|
| 350 | rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) |
|---|
| 351 | #else |
|---|
| 352 | ps_start_GCM(1) = ps(1) |
|---|
| 353 | #endif |
|---|
| 354 | |
|---|
| 355 | ! In the gcm, these values are given to the physic by the dynamic. |
|---|
| 356 | ! Here we simply read them in the startfi_evol.nc file |
|---|
| 357 | status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) |
|---|
| 358 | |
|---|
| 359 | allocate(longitude(ngrid)) |
|---|
| 360 | allocate(latitude(ngrid)) |
|---|
| 361 | allocate(cell_area(ngrid)) |
|---|
| 362 | |
|---|
| 363 | status = nf90_inq_varid(ncid,"longitude",lonvarid) |
|---|
| 364 | status = nf90_get_var(ncid,lonvarid,longitude) |
|---|
| 365 | |
|---|
| 366 | status = nf90_inq_varid(ncid,"latitude",latvarid) |
|---|
| 367 | status = nf90_get_var(ncid,latvarid,latitude) |
|---|
| 368 | |
|---|
| 369 | status = nf90_inq_varid(ncid,"area",areavarid) |
|---|
| 370 | status = nf90_get_var(ncid,areavarid,cell_area) |
|---|
| 371 | |
|---|
| 372 | status = nf90_inq_varid(ncid,"soildepth",sdvarid) |
|---|
| 373 | status = nf90_get_var(ncid,sdvarid,mlayer) |
|---|
| 374 | |
|---|
| 375 | status = nf90_close(ncid) |
|---|
| 376 | |
|---|
| 377 | ! I_b.2 READ startfi_evol.nc |
|---|
| 378 | ! First we read the initial state (starfi.nc) |
|---|
| 379 | #ifndef CPP_STD |
|---|
| 380 | call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, & |
|---|
| 381 | tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar, & |
|---|
| 382 | watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist) |
|---|
| 383 | |
|---|
| 384 | ! Remove unphysical values of surface tracer |
|---|
| 385 | do i = 1,ngrid |
|---|
| 386 | do nnq = 1,nqtot |
|---|
| 387 | do islope = 1,nslope |
|---|
| 388 | if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0. |
|---|
| 389 | enddo |
|---|
| 390 | enddo |
|---|
| 391 | enddo |
|---|
| 392 | |
|---|
| 393 | call surfini(ngrid,qsurf) |
|---|
| 394 | #else |
|---|
| 395 | call phys_state_var_init(nq) |
|---|
| 396 | if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) |
|---|
| 397 | call initracer(ngrid,nq) |
|---|
| 398 | call iniaerosol() |
|---|
| 399 | allocate(tsurf_read_generic(ngrid)) |
|---|
| 400 | allocate(qsurf_read_generic(ngrid,nq)) |
|---|
| 401 | allocate(tsoil_read_generic(ngrid,nsoilmx)) |
|---|
| 402 | allocate(emis_read_generic(ngrid)) |
|---|
| 403 | allocate(tsurf(ngrid,1)) |
|---|
| 404 | allocate(qsurf(ngrid,nq,1)) |
|---|
| 405 | allocate(tsoil(ngrid,nsoilmx,1)) |
|---|
| 406 | allocate(emis(ngrid,1)) |
|---|
| 407 | allocate(watercap(ngrid,1)) |
|---|
| 408 | allocate(watercaptag(ngrid)) |
|---|
| 409 | allocate(albedo_read_generic(ngrid,2)) |
|---|
| 410 | allocate(albedo(ngrid,2,1)) |
|---|
| 411 | allocate(inertiesoil(ngrid,nsoilmx,1)) |
|---|
| 412 | call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, & |
|---|
| 413 | tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & |
|---|
| 414 | qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & |
|---|
| 415 | tslab, tsea_ice,sea_ice) |
|---|
| 416 | call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) |
|---|
| 417 | |
|---|
| 418 | nslope = 1 |
|---|
| 419 | call ini_comslope_h(ngrid,1) |
|---|
| 420 | |
|---|
| 421 | qsurf(:,:,1) = qsurf_read_generic(:,:) |
|---|
| 422 | tsurf(:,1) = tsurf_read_generic(:) |
|---|
| 423 | tsoil(:,:,1) = tsoil_read_generic(:,:) |
|---|
| 424 | emis(:,1) = emis_read_generic(:) |
|---|
| 425 | watercap(:,1) = 0. |
|---|
| 426 | watercaptag(:) = .false. |
|---|
| 427 | albedo(:,1,1) = albedo_read_generic(:,1) |
|---|
| 428 | albedo(:,2,1) = albedo_read_generic(:,2) |
|---|
| 429 | inertiesoil(:,:,1) = inertiedat(:,:) |
|---|
| 430 | |
|---|
| 431 | if (nslope == 1) then |
|---|
| 432 | def_slope(1) = 0 |
|---|
| 433 | def_slope(2) = 0 |
|---|
| 434 | def_slope_mean = 0 |
|---|
| 435 | subslope_dist(:,1) = 1. |
|---|
| 436 | endif |
|---|
| 437 | |
|---|
| 438 | ! Remove unphysical values of surface tracer |
|---|
| 439 | do i = 1,ngrid |
|---|
| 440 | do nnq = 1,nqtot |
|---|
| 441 | qsurf(i,nnq,1) = qsurf_read_generic(i,nnq) |
|---|
| 442 | if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0. |
|---|
| 443 | enddo |
|---|
| 444 | enddo |
|---|
| 445 | #endif |
|---|
| 446 | |
|---|
| 447 | do nnq = 1,nqtot ! Why not using ini_tracer ? |
|---|
| 448 | if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq |
|---|
| 449 | if (noms(nnq) == "h2o_vap") then |
|---|
| 450 | igcm_h2o_vap = nnq |
|---|
| 451 | mmol(igcm_h2o_vap)=18. |
|---|
| 452 | endif |
|---|
| 453 | if (noms(nnq) == "co2") igcm_co2 = nnq |
|---|
| 454 | enddo |
|---|
| 455 | r = 8.314511*1000./mugaz |
|---|
| 456 | |
|---|
| 457 | !------------------------ |
|---|
| 458 | ! I Initialization |
|---|
| 459 | ! I_c Subslope parametrisation |
|---|
| 460 | !------------------------ |
|---|
| 461 | ! Define some slope statistics |
|---|
| 462 | iflat = 1 |
|---|
| 463 | do islope = 2,nslope |
|---|
| 464 | if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope |
|---|
| 465 | enddo |
|---|
| 466 | |
|---|
| 467 | write(*,*) 'Flat slope for islope = ',iflat |
|---|
| 468 | write(*,*) 'corresponding criterium = ',def_slope_mean(iflat) |
|---|
| 469 | |
|---|
| 470 | allocate(flag_co2flow(ngrid,nslope)) |
|---|
| 471 | allocate(flag_co2flow_mesh(ngrid)) |
|---|
| 472 | allocate(flag_h2oflow(ngrid,nslope)) |
|---|
| 473 | allocate(flag_h2oflow_mesh(ngrid)) |
|---|
| 474 | |
|---|
| 475 | flag_co2flow(:,:) = 0 |
|---|
| 476 | flag_co2flow_mesh(:) = 0 |
|---|
| 477 | flag_h2oflow(:,:) = 0 |
|---|
| 478 | flag_h2oflow_mesh(:) = 0 |
|---|
| 479 | |
|---|
| 480 | !------------------------ |
|---|
| 481 | ! I Initialization |
|---|
| 482 | ! I_d READ GCM data and convert to the physical grid |
|---|
| 483 | !------------------------ |
|---|
| 484 | ! 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 |
|---|
| 485 | call nb_time_step_GCM("data_GCM_Y1.nc",timelen) |
|---|
| 486 | |
|---|
| 487 | allocate(tsoil_ave(ngrid,nsoilmx,nslope)) |
|---|
| 488 | allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope)) |
|---|
| 489 | allocate(vmr_co2_gcm(ngrid,timelen)) |
|---|
| 490 | allocate(ps_timeseries(ngrid,timelen)) |
|---|
| 491 | allocate(min_co2_ice_1(ngrid,nslope)) |
|---|
| 492 | allocate(min_h2o_ice_1(ngrid,nslope)) |
|---|
| 493 | allocate(min_co2_ice_2(ngrid,nslope)) |
|---|
| 494 | allocate(min_h2o_ice_2(ngrid,nslope)) |
|---|
| 495 | allocate(tsurf_ave_yr1(ngrid,nslope)) |
|---|
| 496 | allocate(tsurf_ave(ngrid,nslope)) |
|---|
| 497 | allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen)) |
|---|
| 498 | allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen)) |
|---|
| 499 | allocate(q_co2_PEM_phys(ngrid,timelen)) |
|---|
| 500 | allocate(q_h2o_PEM_phys(ngrid,timelen)) |
|---|
| 501 | allocate(co2_ice_GCM(ngrid,nslope,timelen)) |
|---|
| 502 | allocate(watersurf_density_ave(ngrid,nslope)) |
|---|
| 503 | allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) |
|---|
| 504 | allocate(Tsurfave_before_saved(ngrid,nslope)) |
|---|
| 505 | allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) |
|---|
| 506 | allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) |
|---|
| 507 | allocate(delta_co2_adsorbded(ngrid)) |
|---|
| 508 | allocate(porefillingice_thickness_prev_iter(ngrid,nslope)) |
|---|
| 509 | allocate(delta_h2o_icetablesublim(ngrid)) |
|---|
| 510 | allocate(delta_h2o_adsorbded(ngrid)) |
|---|
| 511 | allocate(vmr_co2_pem_phys(ngrid,timelen)) |
|---|
| 512 | |
|---|
| 513 | write(*,*) "Downloading data Y1..." |
|---|
| 514 | 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, & |
|---|
| 515 | tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & |
|---|
| 516 | co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries) |
|---|
| 517 | write(*,*) "Downloading data Y1 done" |
|---|
| 518 | |
|---|
| 519 | ! 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 |
|---|
| 520 | write(*,*) "Downloading data Y2" |
|---|
| 521 | 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, & |
|---|
| 522 | tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & |
|---|
| 523 | co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries) |
|---|
| 524 | write(*,*) "Downloading data Y2 done" |
|---|
| 525 | |
|---|
| 526 | !------------------------ |
|---|
| 527 | ! I Initialization |
|---|
| 528 | ! I_e Initialization of the PEM variables and soil |
|---|
| 529 | !------------------------ |
|---|
| 530 | call end_comsoil_h_PEM |
|---|
| 531 | call ini_comsoil_h_PEM(ngrid,nslope) |
|---|
| 532 | call end_adsorption_h_PEM |
|---|
| 533 | call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) |
|---|
| 534 | call end_ice_table_porefilling |
|---|
| 535 | call ini_ice_table_porefilling(ngrid,nslope) |
|---|
| 536 | |
|---|
| 537 | if (soil_pem) then |
|---|
| 538 | call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) |
|---|
| 539 | do l=1,nsoilmx |
|---|
| 540 | tsoil_PEM(:,l,:)=tsoil_ave(:,l,:) |
|---|
| 541 | tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:) |
|---|
| 542 | watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:) |
|---|
| 543 | enddo |
|---|
| 544 | do l=nsoilmx+1,nsoilmx_PEM |
|---|
| 545 | tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:) |
|---|
| 546 | watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:) |
|---|
| 547 | enddo |
|---|
| 548 | watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen |
|---|
| 549 | endif !soil_pem |
|---|
| 550 | deallocate(tsoil_ave,tsoil_GCM_timeseries) |
|---|
| 551 | |
|---|
| 552 | !------------------------ |
|---|
| 553 | ! I Initialization |
|---|
| 554 | ! I_f Compute tendencies & Save initial situation |
|---|
| 555 | !------------------------ |
|---|
| 556 | allocate(tendencies_co2_ice(ngrid,nslope)) |
|---|
| 557 | allocate(tendencies_co2_ice_ini(ngrid,nslope)) |
|---|
| 558 | allocate(tendencies_h2o_ice(ngrid,nslope)) |
|---|
| 559 | |
|---|
| 560 | ! Compute the tendencies of the evolution of ice over the years |
|---|
| 561 | call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice) |
|---|
| 562 | tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:) |
|---|
| 563 | call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice) |
|---|
| 564 | |
|---|
| 565 | deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2) |
|---|
| 566 | |
|---|
| 567 | !------------------------ |
|---|
| 568 | ! I Initialization |
|---|
| 569 | ! I_g Save initial PCM situation |
|---|
| 570 | !------------------------ |
|---|
| 571 | allocate(initial_co2_ice_sublim(ngrid,nslope)) |
|---|
| 572 | allocate(initial_co2_ice(ngrid,nslope)) |
|---|
| 573 | allocate(initial_h2o_ice(ngrid,nslope)) |
|---|
| 574 | |
|---|
| 575 | ! We save the places where water ice is sublimating |
|---|
| 576 | ! We compute the surface of water ice sublimating |
|---|
| 577 | ini_surf_co2 = 0. |
|---|
| 578 | ini_surf_h2o = 0. |
|---|
| 579 | Total_surface = 0. |
|---|
| 580 | do i = 1,ngrid |
|---|
| 581 | Total_surface = Total_surface+cell_area(i) |
|---|
| 582 | do islope = 1,nslope |
|---|
| 583 | if (tendencies_co2_ice(i,islope) < 0) then |
|---|
| 584 | initial_co2_ice_sublim(i,islope) = 1. |
|---|
| 585 | ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope) |
|---|
| 586 | else |
|---|
| 587 | initial_co2_ice_sublim(i,islope) = 0. |
|---|
| 588 | endif |
|---|
| 589 | if (qsurf(i,igcm_co2,islope) > 0) then |
|---|
| 590 | initial_co2_ice(i,islope) = 1. |
|---|
| 591 | else |
|---|
| 592 | initial_co2_ice(i,islope) = 0. |
|---|
| 593 | endif |
|---|
| 594 | if (tendencies_h2o_ice(i,islope) < 0) then |
|---|
| 595 | initial_h2o_ice(i,islope) = 1. |
|---|
| 596 | ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope) |
|---|
| 597 | else |
|---|
| 598 | initial_h2o_ice(i,islope) = 0. |
|---|
| 599 | endif |
|---|
| 600 | enddo |
|---|
| 601 | enddo |
|---|
| 602 | |
|---|
| 603 | write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2 |
|---|
| 604 | write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o |
|---|
| 605 | write(*,*) "Total surface of the planet=", Total_surface |
|---|
| 606 | allocate(zplev_gcm(ngrid,nlayer + 1)) |
|---|
| 607 | |
|---|
| 608 | do l = 1,nlayer + 1 |
|---|
| 609 | do ig = 1,ngrid |
|---|
| 610 | zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) |
|---|
| 611 | enddo |
|---|
| 612 | enddo |
|---|
| 613 | |
|---|
| 614 | global_ave_press_old = 0. |
|---|
| 615 | do i = 1,ngrid |
|---|
| 616 | global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface |
|---|
| 617 | enddo |
|---|
| 618 | |
|---|
| 619 | global_ave_press_GCM = global_ave_press_old |
|---|
| 620 | global_ave_press_new = global_ave_press_old |
|---|
| 621 | write(*,*) "Initial global average pressure=", global_ave_press_GCM |
|---|
| 622 | |
|---|
| 623 | !------------------------ |
|---|
| 624 | ! I Initialization |
|---|
| 625 | ! I_h Read the PEMstart |
|---|
| 626 | !------------------------ |
|---|
| 627 | call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, & |
|---|
| 628 | porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, & |
|---|
| 629 | tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), & |
|---|
| 630 | qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, & |
|---|
| 631 | co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) |
|---|
| 632 | |
|---|
| 633 | delta_h2o_icetablesublim(:) = 0. |
|---|
| 634 | |
|---|
| 635 | do ig = 1,ngrid |
|---|
| 636 | do islope = 1,nslope |
|---|
| 637 | 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.) |
|---|
| 638 | enddo |
|---|
| 639 | enddo |
|---|
| 640 | |
|---|
| 641 | if (adsorption_pem) then |
|---|
| 642 | totmassco2_adsorbded = 0. |
|---|
| 643 | totmassh2o_adsorbded = 0. |
|---|
| 644 | do ig = 1,ngrid |
|---|
| 645 | do islope = 1, nslope |
|---|
| 646 | do l = 1,nsoilmx_PEM - 1 |
|---|
| 647 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & |
|---|
| 648 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
|---|
| 649 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & |
|---|
| 650 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
|---|
| 651 | enddo |
|---|
| 652 | enddo |
|---|
| 653 | enddo |
|---|
| 654 | |
|---|
| 655 | write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded |
|---|
| 656 | write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded |
|---|
| 657 | endif ! adsorption |
|---|
| 658 | deallocate(tsurf_ave_yr1) |
|---|
| 659 | |
|---|
| 660 | !------------------------ |
|---|
| 661 | ! I Initialization |
|---|
| 662 | ! I_i Compute orbit criterion |
|---|
| 663 | !------------------------ |
|---|
| 664 | #ifndef CPP_STD |
|---|
| 665 | call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
|---|
| 666 | #else |
|---|
| 667 | call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) |
|---|
| 668 | #endif |
|---|
| 669 | |
|---|
| 670 | if (evol_orbit_pem) then |
|---|
| 671 | call orbit_param_criterion(i_myear,year_iter_max) |
|---|
| 672 | else |
|---|
| 673 | year_iter_max = Max_iter_pem |
|---|
| 674 | endif |
|---|
| 675 | !-------------------------- END INITIALIZATION ------------------------- |
|---|
| 676 | |
|---|
| 677 | !-------------------------------- RUN ---------------------------------- |
|---|
| 678 | !------------------------ |
|---|
| 679 | ! II Run |
|---|
| 680 | ! II_a Update pressure, ice and tracers |
|---|
| 681 | !------------------------ |
|---|
| 682 | year_iter = 0 |
|---|
| 683 | |
|---|
| 684 | do while (year_iter < year_iter_max .and. i_myear < n_myear) |
|---|
| 685 | ! II.a.1. Compute updated global pressure |
|---|
| 686 | write(*,*) "Recomputing the new pressure..." |
|---|
| 687 | |
|---|
| 688 | do i = 1,ngrid |
|---|
| 689 | do islope = 1,nslope |
|---|
| 690 | 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 |
|---|
| 691 | enddo |
|---|
| 692 | enddo |
|---|
| 693 | write(*,*) 'Global average pressure old time step',global_ave_press_old |
|---|
| 694 | call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) |
|---|
| 695 | |
|---|
| 696 | if (adsorption_pem) then |
|---|
| 697 | do i = 1,ngrid |
|---|
| 698 | global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface |
|---|
| 699 | enddo |
|---|
| 700 | write(*,*) 'Global average pressure old time step',global_ave_press_old |
|---|
| 701 | write(*,*) 'Global average pressure new time step',global_ave_press_new |
|---|
| 702 | endif |
|---|
| 703 | |
|---|
| 704 | ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) |
|---|
| 705 | allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen)) |
|---|
| 706 | write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..." |
|---|
| 707 | do l = 1,nlayer + 1 |
|---|
| 708 | do ig = 1,ngrid |
|---|
| 709 | zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
|---|
| 710 | enddo |
|---|
| 711 | enddo |
|---|
| 712 | |
|---|
| 713 | ! II.a.3. Surface pressure timeseries |
|---|
| 714 | write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..." |
|---|
| 715 | do ig = 1,ngrid |
|---|
| 716 | ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old |
|---|
| 717 | enddo |
|---|
| 718 | |
|---|
| 719 | ! II.a.4. New pressure levels timeseries |
|---|
| 720 | allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen)) |
|---|
| 721 | write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..." |
|---|
| 722 | do l = 1,nlayer + 1 |
|---|
| 723 | do ig = 1,ngrid |
|---|
| 724 | zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
|---|
| 725 | enddo |
|---|
| 726 | enddo |
|---|
| 727 | |
|---|
| 728 | ! II.a.5. Tracers timeseries |
|---|
| 729 | write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..." |
|---|
| 730 | |
|---|
| 731 | l = 1 |
|---|
| 732 | do ig = 1,ngrid |
|---|
| 733 | do t = 1,timelen |
|---|
| 734 | 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))/ & |
|---|
| 735 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) |
|---|
| 736 | if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30 |
|---|
| 737 | if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1. |
|---|
| 738 | enddo |
|---|
| 739 | enddo |
|---|
| 740 | |
|---|
| 741 | do ig = 1,ngrid |
|---|
| 742 | do t = 1, timelen |
|---|
| 743 | 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))/ & |
|---|
| 744 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & |
|---|
| 745 | + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & |
|---|
| 746 | - (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ & |
|---|
| 747 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) |
|---|
| 748 | if (q_co2_PEM_phys(ig,t) < 0) then |
|---|
| 749 | q_co2_PEM_phys(ig,t) = 1.e-30 |
|---|
| 750 | elseif (q_co2_PEM_phys(ig,t) > 1) then |
|---|
| 751 | q_co2_PEM_phys(ig,t) = 1. |
|---|
| 752 | endif |
|---|
| 753 | mmean=1/(A*q_co2_PEM_phys(ig,t) + B) |
|---|
| 754 | vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 |
|---|
| 755 | enddo |
|---|
| 756 | enddo |
|---|
| 757 | |
|---|
| 758 | deallocate(zplev_new_timeseries,zplev_old_timeseries) |
|---|
| 759 | |
|---|
| 760 | !------------------------ |
|---|
| 761 | ! II Run |
|---|
| 762 | ! II_b Evolution of the ice |
|---|
| 763 | !------------------------ |
|---|
| 764 | write(*,*) "Evolution of h2o ice" |
|---|
| 765 | 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) |
|---|
| 766 | |
|---|
| 767 | write(*,*) "Evolution of co2 ice" |
|---|
| 768 | call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) |
|---|
| 769 | |
|---|
| 770 | do islope=1, nslope |
|---|
| 771 | write(str2(1:2),'(i2.2)') islope |
|---|
| 772 | call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) |
|---|
| 773 | call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) |
|---|
| 774 | call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) |
|---|
| 775 | call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) |
|---|
| 776 | enddo |
|---|
| 777 | |
|---|
| 778 | !------------------------ |
|---|
| 779 | ! II Run |
|---|
| 780 | ! II_c CO2 & H2O glaciers flows |
|---|
| 781 | !------------------------ |
|---|
| 782 | write(*,*) "CO2 glacier flows" |
|---|
| 783 | |
|---|
| 784 | if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, & |
|---|
| 785 | global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) |
|---|
| 786 | |
|---|
| 787 | write(*,*) "H2O glacier flows" |
|---|
| 788 | |
|---|
| 789 | 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) |
|---|
| 790 | |
|---|
| 791 | do islope = 1,nslope |
|---|
| 792 | write(str2(1:2),'(i2.2)') islope |
|---|
| 793 | call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) |
|---|
| 794 | call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) |
|---|
| 795 | call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) |
|---|
| 796 | enddo |
|---|
| 797 | |
|---|
| 798 | !------------------------ |
|---|
| 799 | ! II Run |
|---|
| 800 | ! II_d Update surface and soil temperatures |
|---|
| 801 | !------------------------ |
|---|
| 802 | ! II_d.1 Update Tsurf |
|---|
| 803 | write(*,*) "Updating the new Tsurf" |
|---|
| 804 | bool_sublim = .false. |
|---|
| 805 | Tsurfave_before_saved(:,:) = tsurf_ave(:,:) |
|---|
| 806 | do ig = 1,ngrid |
|---|
| 807 | do islope = 1,nslope |
|---|
| 808 | 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 |
|---|
| 809 | if (latitude_deg(ig) > 0) then |
|---|
| 810 | do ig_loop = ig,ngrid |
|---|
| 811 | do islope_loop = islope,iflat,-1 |
|---|
| 812 | if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then |
|---|
| 813 | tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) |
|---|
| 814 | bool_sublim = .true. |
|---|
| 815 | exit |
|---|
| 816 | endif |
|---|
| 817 | enddo |
|---|
| 818 | if (bool_sublim) exit |
|---|
| 819 | enddo |
|---|
| 820 | else |
|---|
| 821 | do ig_loop = ig,1,-1 |
|---|
| 822 | do islope_loop = islope,iflat |
|---|
| 823 | if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then |
|---|
| 824 | tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) |
|---|
| 825 | bool_sublim = .true. |
|---|
| 826 | exit |
|---|
| 827 | endif |
|---|
| 828 | enddo |
|---|
| 829 | if (bool_sublim) exit |
|---|
| 830 | enddo |
|---|
| 831 | endif |
|---|
| 832 | initial_co2_ice(ig,islope) = 0 |
|---|
| 833 | if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then |
|---|
| 834 | albedo(ig,1,islope) = albedo_h2o_frost |
|---|
| 835 | albedo(ig,2,islope) = albedo_h2o_frost |
|---|
| 836 | else |
|---|
| 837 | albedo(ig,1,islope) = albedodat(ig) |
|---|
| 838 | albedo(ig,2,islope) = albedodat(ig) |
|---|
| 839 | emis(ig,islope) = emissiv |
|---|
| 840 | endif |
|---|
| 841 | 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 |
|---|
| 842 | ave = 0. |
|---|
| 843 | do t = 1,timelen |
|---|
| 844 | if (co2_ice_GCM(ig,islope,t) > 1.e-3) then |
|---|
| 845 | ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) |
|---|
| 846 | else |
|---|
| 847 | ave = ave + tsurf_GCM_timeseries(ig,islope,t) |
|---|
| 848 | endif |
|---|
| 849 | enddo |
|---|
| 850 | tsurf_ave(ig,islope) = ave/timelen |
|---|
| 851 | ! set the surface albedo to be the ice albedo |
|---|
| 852 | if (latitude_deg(ig) > 0) then |
|---|
| 853 | icap = 1 |
|---|
| 854 | else |
|---|
| 855 | icap = 2 |
|---|
| 856 | endif |
|---|
| 857 | albedo(ig,1,islope) = albedice(icap) |
|---|
| 858 | albedo(ig,2,islope) = albedice(icap) |
|---|
| 859 | emis(ig,islope) = emisice(icap) |
|---|
| 860 | endif |
|---|
| 861 | enddo |
|---|
| 862 | enddo |
|---|
| 863 | |
|---|
| 864 | do t = 1,timelen |
|---|
| 865 | tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) |
|---|
| 866 | enddo |
|---|
| 867 | ! for the start |
|---|
| 868 | do ig = 1,ngrid |
|---|
| 869 | do islope = 1,nslope |
|---|
| 870 | tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope)) |
|---|
| 871 | enddo |
|---|
| 872 | enddo |
|---|
| 873 | |
|---|
| 874 | do islope = 1,nslope |
|---|
| 875 | write(str2(1:2),'(i2.2)') islope |
|---|
| 876 | call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) |
|---|
| 877 | enddo |
|---|
| 878 | |
|---|
| 879 | if (soil_pem) then |
|---|
| 880 | |
|---|
| 881 | ! II_d.2 Update soil temperature |
|---|
| 882 | allocate(TI_locslope(ngrid,nsoilmx_PEM)) |
|---|
| 883 | allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) |
|---|
| 884 | allocate(Tsurf_locslope(ngrid)) |
|---|
| 885 | write(*,*)"Updating soil temperature" |
|---|
| 886 | |
|---|
| 887 | ! Soil averaged |
|---|
| 888 | do islope = 1,nslope |
|---|
| 889 | TI_locslope(:,:) = TI_PEM(:,:,islope) |
|---|
| 890 | do t = 1,timelen |
|---|
| 891 | Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) |
|---|
| 892 | Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t) |
|---|
| 893 | call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) |
|---|
| 894 | call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) |
|---|
| 895 | tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) |
|---|
| 896 | do ig = 1,ngrid |
|---|
| 897 | do isoil = 1,nsoilmx_PEM |
|---|
| 898 | 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) |
|---|
| 899 | if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) |
|---|
| 900 | enddo |
|---|
| 901 | enddo |
|---|
| 902 | enddo |
|---|
| 903 | enddo |
|---|
| 904 | tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen |
|---|
| 905 | watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen |
|---|
| 906 | |
|---|
| 907 | write(*,*) "Update of soil temperature done" |
|---|
| 908 | |
|---|
| 909 | deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope) |
|---|
| 910 | write(*,*) "Compute ice table" |
|---|
| 911 | |
|---|
| 912 | ! II_d.3 Update the ice table |
|---|
| 913 | porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:) |
|---|
| 914 | call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) |
|---|
| 915 | |
|---|
| 916 | call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, & |
|---|
| 917 | delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere |
|---|
| 918 | |
|---|
| 919 | write(*,*) "Update soil propreties" |
|---|
| 920 | |
|---|
| 921 | ! II_d.4 Update the soil thermal properties |
|---|
| 922 | call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, & |
|---|
| 923 | porefillingice_depth,porefillingice_thickness,TI_PEM) |
|---|
| 924 | |
|---|
| 925 | ! II_d.5 Update the mass of the regolith adsorbded |
|---|
| 926 | if (adsorption_pem) then |
|---|
| 927 | call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, & |
|---|
| 928 | qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, & |
|---|
| 929 | q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) |
|---|
| 930 | |
|---|
| 931 | totmassco2_adsorbded = 0. |
|---|
| 932 | totmassh2o_adsorbded = 0. |
|---|
| 933 | do ig = 1,ngrid |
|---|
| 934 | do islope =1, nslope |
|---|
| 935 | do l = 1,nsoilmx_PEM - 1 |
|---|
| 936 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & |
|---|
| 937 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & |
|---|
| 938 | cell_area(ig) |
|---|
| 939 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & |
|---|
| 940 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & |
|---|
| 941 | cell_area(ig) |
|---|
| 942 | enddo |
|---|
| 943 | enddo |
|---|
| 944 | enddo |
|---|
| 945 | write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded |
|---|
| 946 | write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded |
|---|
| 947 | endif |
|---|
| 948 | endif !soil_pem |
|---|
| 949 | |
|---|
| 950 | !------------------------ |
|---|
| 951 | ! II Run |
|---|
| 952 | ! II_e Update the tendencies |
|---|
| 953 | !------------------------ |
|---|
| 954 | write(*,*) "Adaptation of the new co2 tendencies to the current pressure" |
|---|
| 955 | call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, & |
|---|
| 956 | global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) |
|---|
| 957 | |
|---|
| 958 | !------------------------ |
|---|
| 959 | ! II Run |
|---|
| 960 | ! II_f Checking the stopping criterion |
|---|
| 961 | !------------------------ |
|---|
| 962 | call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) |
|---|
| 963 | |
|---|
| 964 | call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & |
|---|
| 965 | initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope) |
|---|
| 966 | |
|---|
| 967 | year_iter = year_iter + dt_pem |
|---|
| 968 | i_myear = i_myear + dt_pem |
|---|
| 969 | |
|---|
| 970 | write(*,*) "Checking all the stopping criterion." |
|---|
| 971 | if (STOPPING_water) then |
|---|
| 972 | write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water |
|---|
| 973 | criterion_stop = 1 |
|---|
| 974 | endif |
|---|
| 975 | if (STOPPING_1_water) then |
|---|
| 976 | write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water |
|---|
| 977 | criterion_stop = 1 |
|---|
| 978 | endif |
|---|
| 979 | if (STOPPING_co2) then |
|---|
| 980 | write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 |
|---|
| 981 | criterion_stop = 2 |
|---|
| 982 | endif |
|---|
| 983 | if (STOPPING_pressure) then |
|---|
| 984 | write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure |
|---|
| 985 | criterion_stop = 3 |
|---|
| 986 | endif |
|---|
| 987 | if (year_iter >= year_iter_max) then |
|---|
| 988 | write(*,*) "STOPPING because maximum number of iterations reached" |
|---|
| 989 | criterion_stop = 4 |
|---|
| 990 | endif |
|---|
| 991 | if (i_myear >= n_myear) then |
|---|
| 992 | write(*,*) "STOPPING because maximum number of Martian years to be simulated reached" |
|---|
| 993 | criterion_stop = 5 |
|---|
| 994 | endif |
|---|
| 995 | |
|---|
| 996 | if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then |
|---|
| 997 | exit |
|---|
| 998 | else |
|---|
| 999 | write(*,*) "We continue!" |
|---|
| 1000 | write(*,*) "Number of iterations of the PEM: year_iter =", year_iter |
|---|
| 1001 | write(*,*) "Number of simulated Martian years: i_myear =", i_myear |
|---|
| 1002 | endif |
|---|
| 1003 | |
|---|
| 1004 | global_ave_press_old=global_ave_press_new |
|---|
| 1005 | |
|---|
| 1006 | enddo ! big time iteration loop of the pem |
|---|
| 1007 | !------------------------------ END RUN -------------------------------- |
|---|
| 1008 | |
|---|
| 1009 | !------------------------------- OUTPUT -------------------------------- |
|---|
| 1010 | !------------------------ |
|---|
| 1011 | ! III Output |
|---|
| 1012 | ! III_a Update surface value for the PCM start files |
|---|
| 1013 | !------------------------ |
|---|
| 1014 | ! III_a.1 Ice update (for startfi) |
|---|
| 1015 | |
|---|
| 1016 | ! H2O ice |
|---|
| 1017 | do ig = 1,ngrid |
|---|
| 1018 | if (watercaptag(ig)) then |
|---|
| 1019 | watercap_sum = 0. |
|---|
| 1020 | do islope = 1,nslope |
|---|
| 1021 | 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 |
|---|
| 1022 | 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 |
|---|
| 1023 | else |
|---|
| 1024 | ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap |
|---|
| 1025 | 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.)) |
|---|
| 1026 | qsurf(ig,igcm_h2o_ice,islope)=0. |
|---|
| 1027 | endif |
|---|
| 1028 | watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) |
|---|
| 1029 | watercap(ig,islope) = 0. |
|---|
| 1030 | enddo |
|---|
| 1031 | water_reservoir(ig) = water_reservoir(ig) + watercap_sum |
|---|
| 1032 | endif |
|---|
| 1033 | enddo |
|---|
| 1034 | |
|---|
| 1035 | do ig = 1,ngrid |
|---|
| 1036 | water_sum = 0. |
|---|
| 1037 | do islope = 1,nslope |
|---|
| 1038 | water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) |
|---|
| 1039 | enddo |
|---|
| 1040 | if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. |
|---|
| 1041 | 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 |
|---|
| 1042 | watercaptag(ig) = .true. |
|---|
| 1043 | 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 |
|---|
| 1044 | do islope = 1,nslope |
|---|
| 1045 | qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.) |
|---|
| 1046 | enddo |
|---|
| 1047 | endif |
|---|
| 1048 | else ! let's check that the infinite source of water disapear |
|---|
| 1049 | if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then |
|---|
| 1050 | watercaptag(ig) = .false. |
|---|
| 1051 | do islope = 1,nslope |
|---|
| 1052 | qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) |
|---|
| 1053 | enddo |
|---|
| 1054 | water_reservoir(ig) = 0. |
|---|
| 1055 | endif |
|---|
| 1056 | endif |
|---|
| 1057 | enddo |
|---|
| 1058 | |
|---|
| 1059 | ! CO2 ice |
|---|
| 1060 | do ig = 1,ngrid |
|---|
| 1061 | do islope = 1,nslope |
|---|
| 1062 | if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then |
|---|
| 1063 | perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) |
|---|
| 1064 | qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) |
|---|
| 1065 | albedo(ig,1,islope) = albedo_perenialco2 |
|---|
| 1066 | albedo(ig,2,islope) = albedo_perenialco2 |
|---|
| 1067 | endif |
|---|
| 1068 | enddo |
|---|
| 1069 | enddo |
|---|
| 1070 | |
|---|
| 1071 | ! III_a.2 Tsoil update (for startfi) |
|---|
| 1072 | if (soil_pem) then |
|---|
| 1073 | call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) |
|---|
| 1074 | tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) |
|---|
| 1075 | endif |
|---|
| 1076 | |
|---|
| 1077 | ! III_a.4 Pressure (for start) |
|---|
| 1078 | do i = 1,ip1jmp1 |
|---|
| 1079 | ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM |
|---|
| 1080 | enddo |
|---|
| 1081 | |
|---|
| 1082 | do i = 1,ngrid |
|---|
| 1083 | ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM |
|---|
| 1084 | enddo |
|---|
| 1085 | |
|---|
| 1086 | ! III_a.5 Tracer (for start) |
|---|
| 1087 | allocate(zplev_new(ngrid,nlayer + 1)) |
|---|
| 1088 | |
|---|
| 1089 | do l = 1,nlayer + 1 |
|---|
| 1090 | do ig = 1,ngrid |
|---|
| 1091 | zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) |
|---|
| 1092 | enddo |
|---|
| 1093 | enddo |
|---|
| 1094 | |
|---|
| 1095 | do nnq = 1,nqtot |
|---|
| 1096 | if (noms(nnq) /= "co2") then |
|---|
| 1097 | do l = 1,llm - 1 |
|---|
| 1098 | do ig = 1,ngrid |
|---|
| 1099 | 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)) |
|---|
| 1100 | enddo |
|---|
| 1101 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
|---|
| 1102 | enddo |
|---|
| 1103 | else |
|---|
| 1104 | do l = 1,llm - 1 |
|---|
| 1105 | do ig = 1,ngrid |
|---|
| 1106 | 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)) & |
|---|
| 1107 | + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ & |
|---|
| 1108 | (zplev_new(ig,l) - zplev_new(ig,l + 1)) |
|---|
| 1109 | enddo |
|---|
| 1110 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
|---|
| 1111 | enddo |
|---|
| 1112 | endif |
|---|
| 1113 | enddo |
|---|
| 1114 | |
|---|
| 1115 | ! Conserving the tracers mass for GCM start files |
|---|
| 1116 | do nnq = 1,nqtot |
|---|
| 1117 | do ig = 1,ngrid |
|---|
| 1118 | do l = 1,llm - 1 |
|---|
| 1119 | 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 |
|---|
| 1120 | extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) |
|---|
| 1121 | q(ig,l,nnq)=1. |
|---|
| 1122 | q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2)) |
|---|
| 1123 | write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" |
|---|
| 1124 | endif |
|---|
| 1125 | if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 |
|---|
| 1126 | enddo |
|---|
| 1127 | enddo |
|---|
| 1128 | enddo |
|---|
| 1129 | |
|---|
| 1130 | if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter) |
|---|
| 1131 | |
|---|
| 1132 | !------------------------ |
|---|
| 1133 | ! III Output |
|---|
| 1134 | ! III_b Write restart_evol.nc and restartfi_evol.nc |
|---|
| 1135 | !------------------------ |
|---|
| 1136 | ! III_b.1 Write restart_evol.nc |
|---|
| 1137 | ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys |
|---|
| 1138 | pday = day_ini |
|---|
| 1139 | ztime_fin = time_phys |
|---|
| 1140 | |
|---|
| 1141 | allocate(p(ip1jmp1,nlayer + 1)) |
|---|
| 1142 | #ifndef CPP_1D |
|---|
| 1143 | call pression (ip1jmp1,ap,bp,ps,p) |
|---|
| 1144 | call massdair(p,masse) |
|---|
| 1145 | call dynredem0("restart_evol.nc",day_ini,phis) |
|---|
| 1146 | call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps) |
|---|
| 1147 | write(*,*) "restart_evol.nc has been written" |
|---|
| 1148 | #else |
|---|
| 1149 | do nnq = 1, nqtot |
|---|
| 1150 | call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) |
|---|
| 1151 | enddo |
|---|
| 1152 | #endif |
|---|
| 1153 | |
|---|
| 1154 | ! III_b.2 Write restartfi_evol.nc |
|---|
| 1155 | #ifndef CPP_STD |
|---|
| 1156 | call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & |
|---|
| 1157 | nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & |
|---|
| 1158 | inertiedat,def_slope,subslope_dist) |
|---|
| 1159 | |
|---|
| 1160 | call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & |
|---|
| 1161 | ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & |
|---|
| 1162 | albedo,emis,q2,qsurf,tauscaling,totcloudfrac, & |
|---|
| 1163 | wstar,watercap,perenial_co2ice) |
|---|
| 1164 | #else |
|---|
| 1165 | call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & |
|---|
| 1166 | nlayer,nq,ptimestep,pday,time_phys,cell_area, & |
|---|
| 1167 | albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) |
|---|
| 1168 | |
|---|
| 1169 | call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & |
|---|
| 1170 | ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, & |
|---|
| 1171 | cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & |
|---|
| 1172 | tslab,tsea_ice,sea_ice) |
|---|
| 1173 | #endif |
|---|
| 1174 | write(*,*) "restartfi_evol.nc has been written" |
|---|
| 1175 | |
|---|
| 1176 | !------------------------ |
|---|
| 1177 | ! III Output |
|---|
| 1178 | ! III_c Write restartfi_PEM.nc |
|---|
| 1179 | !------------------------ |
|---|
| 1180 | call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & |
|---|
| 1181 | float(day_ini),0.,nslope,def_slope,subslope_dist) |
|---|
| 1182 | |
|---|
| 1183 | |
|---|
| 1184 | call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, & |
|---|
| 1185 | TI_PEM, porefillingice_depth,porefillingice_thickness, & |
|---|
| 1186 | co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) |
|---|
| 1187 | write(*,*) "restartfi_PEM.nc has been written" |
|---|
| 1188 | |
|---|
| 1189 | call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear) |
|---|
| 1190 | |
|---|
| 1191 | write(*,*) "The PEM has run for", year_iter, "Martian years." |
|---|
| 1192 | write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years." |
|---|
| 1193 | write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years." |
|---|
| 1194 | write(*,*) "LL & RV & JBC: so far, so good!" |
|---|
| 1195 | |
|---|
| 1196 | deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys) |
|---|
| 1197 | deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved) |
|---|
| 1198 | deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) |
|---|
| 1199 | deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) |
|---|
| 1200 | !----------------------------- END OUTPUT ------------------------------ |
|---|
| 1201 | |
|---|
| 1202 | END PROGRAM pem |
|---|
| 1203 | |
|---|
