- Timestamp:
- Aug 14, 2023, 3:22:41 PM (16 months ago)
- Location:
- trunk
- Files:
-
- 2 edited
-
LMDZ.COMMON/libf/evolution/pem.F90 (modified) (2 diffs)
-
LMDZ.MARS/changelog.txt (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3022 r3028 1 2 !------------------------ 3 4 ! I Initialisation 5 ! I_a READ run.def 1 !------------------------ 2 ! I Initialization 3 ! I_a READ run.def 6 4 ! I_b READ of start_evol.nc and starfi_evol.nc 7 5 ! I_c Subslope parametrisation 8 ! I_d READ GCM data and convert to the physical grid 9 ! I_e Initiali sation of the PEM variable and soil6 ! I_d READ GCM data and convert to the physical grid 7 ! I_e Initialization of the PEM variable and soil 10 8 ! I_f Compute tendencies & Save initial situation 11 9 ! I_g Save initial PCM situation … … 14 12 15 13 ! II Run 16 ! II_a update pressure, ice and tracers14 ! II_a Update pressure, ice and tracers 17 15 ! II_b Evolution of the ice 18 ! II_c CO2 glaciers flows16 ! II_c CO2 & H2O glaciers flows 19 17 ! II_d Update surface and soil temperatures 20 ! II_e Update the tendencies 18 ! II_e Update the tendencies 21 19 ! II_f Checking the stopping criterion 22 20 23 21 ! III Output 24 22 ! III_a Update surface value for the PCM start files 25 ! III_b Write start and starfi.nc 26 ! III_c Write start_pem 27 23 ! III_b Write restart_evol.nc and restartfi_evol.nc 24 ! III_c Write restartfi_PEM.nc 28 25 !------------------------ 29 26 30 27 PROGRAM pem 31 28 32 ! 1: Modules needed for reading and writting startfi: 33 use phyetat0_mod, only: phyetat0 34 use phyredem, only: physdem0, physdem1 35 use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, & 36 nf90_get_var, nf90_inq_varid, nf90_inq_dimid, & 37 nf90_inquire_dimension,nf90_close 38 ! netcdf is needed to read info like lat and lon in the physiq file 39 use turb_mod, only: q2, wstar 40 ! 1a: Modules specific from the Marsian physiq 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 temps_mod_evol, only: dt_pem, evol_orbit_pem, Max_iter_pem 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 60 use soil_thermalproperties_mod, only: update_soil_thermalproperties 61 41 62 #ifndef CPP_STD 42 use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil43 use surfdat_h, only: tsurf, emis,&44 qsurf,watercap, ini_surfdat_h,&45 albedodat, zmea, zstd, zsig, zgam, zthe,&46 hmons, summit, base,albedo_h2o_frost, &47 frost_albedo_threshold,emissiv,watercaptag,perenial_co2ice48 use dimradmars_mod, only: totcloudfrac, albedo49 use dust_param_mod, only: tauscaling50 use tracer_mod, only: noms,igcm_h2o_ice,igcm_co2,mmol,igcm_h2o_vap ! tracer names and molar masses51 #else 52 ! 1b: Modules specific from the Generic physiq 53 use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa 54 use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &55 emissiv56 use tracer_h, only: noms,igcm_h2o_ice,igcm_co2 ! tracer names57 use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &58 PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, &59 ALBEDO_BAREGROUND,ALBEDO_CO2_ICE_SPECTV, phys_state_var_init60 use aerosol_mod, only : iniaerosol63 use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil 64 use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, & 65 albedodat, zmea, zstd, zsig, zgam, zthe, & 66 hmons, summit, base,albedo_h2o_frost, & 67 frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice 68 use dimradmars_mod, only: totcloudfrac, albedo 69 use dust_param_mod, only: tauscaling 70 use tracer_mod, only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses 71 use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master 72 use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit 73 use comcstfi_h, only: r, mugaz 74 #else 75 use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names 76 use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & 77 PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, & 78 ALBEDO_CO2_ICE_SPECTV, phys_state_var_init 79 use aerosol_mod, only: iniaerosol 80 use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit 81 use comcstfi_mod, only: r, mugaz 61 82 #endif 62 ! 1c: Modules specific from the 1d marsian physiq 83 63 84 #ifndef CPP_1D 64 USE iniphysiq_mod, ONLY: iniphysiq65 USE control_mod, ONLY: iphysiq, day_step,nsplit_phys85 use iniphysiq_mod, only: iniphysiq 86 use control_mod, only: iphysiq, day_step, nsplit_phys 66 87 #else 67 use time_phylmdz_mod, only: daysec, iphysiq,day_step88 use time_phylmdz_mod, only: daysec, iphysiq, day_step, dtphys 68 89 #endif 69 90 91 #ifdef CPP_1D 92 use regular_lonlat_mod, only: init_regular_lonlat 93 use physics_distribution_mod, only: init_physics_distribution 94 use mod_grid_phy_lmdz, only: regular_lonlat 95 use init_phys_1d_mod, only: init_phys_1d 96 #endif 97 98 IMPLICIT NONE 99 100 include "dimensions.h" 101 include "paramet.h" 102 include "comgeom.h" 103 include "iniprint.h" 104 105 integer ngridmx 106 parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) 107 108 ! Same variable names as in the GCM 109 integer :: ngrid ! Number of physical grid points 110 integer :: nlayer ! Number of vertical layer 111 integer :: nq ! Number of tracer 112 integer :: day_ini ! First day of the simulation 113 real :: pday ! Physical day 114 real :: time_phys ! Same as GCM 115 real :: ptimestep ! Same as GCM 116 real :: ztime_fin ! Same as GCM 117 118 ! Variables to read start.nc 119 character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM 120 121 ! Dynamic variables 122 real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants 123 real :: teta(ip1jmp1,llm) ! temperature potentielle 124 real, dimension(:,:,:), allocatable :: q ! champs advectes 125 real :: ps(ip1jmp1) ! pression au sol 126 real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol 127 real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées 128 real :: masse(ip1jmp1,llm) ! masse d'air 129 real :: phis(ip1jmp1) ! geopotentiel au sol 130 real :: time_0 131 132 ! Variables to read starfi.nc 133 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM 134 character*2 str2 135 integer :: ncid, varid,status ! Variable for handling opening of files 136 integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files 137 integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files 138 integer :: apvarid, bpvarid ! Variable ID for Netcdf files 139 140 ! Variables to read starfi.nc and write restartfi.nc 141 real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi 142 real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi 143 real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart 144 real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart 145 real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi 146 real :: Total_surface ! Total surface of the planet 147 148 ! Variables for h2o_ice evolution 149 real :: ini_surf_h2o ! Initial surface of sublimating h2o ice 150 real :: ini_surf_co2 ! Initial surface of sublimating co2 ice 151 real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] 152 real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step 153 real :: global_ave_press_new ! constant: Global average pressure of current time step 154 real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] 155 real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] 156 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] 157 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step 158 logical :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? 159 logical :: STOPPING_1_water ! Logical : is there still water ice to sublimate? 160 logical :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? 161 logical :: STOPPING_pressure ! Logical : is the criterion (% of change in the surface pressure) reached? 162 integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param 163 real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 164 real, allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] 165 real, allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM 166 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] 167 real, allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] 168 integer :: timelen ! # time samples 169 real :: ave ! intermediate varibale to compute average 170 real, allocatable :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 171 real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 172 173 ! Variables for slopes 174 real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] 175 real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] 176 real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] 177 real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] 178 real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] 179 real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice 180 real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state 181 real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state 182 real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year 183 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 184 real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial h2o ice 185 real, dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope): Flag where there is a CO2 glacier flow 186 real, dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) : Flag where there is a CO2 glacier flow 187 real, dimension(:,:), allocatable :: flag_h2oflow(:,:) !(ngrid,nslope): Flag where there is a H2O glacier flow 188 real, dimension(:), allocatable :: flag_h2oflow_mesh(:) !(ngrid) : Flag where there is a H2O glacier flow 189 190 ! Variables for surface and soil 191 real, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] 192 real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] 193 real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] 194 real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 195 real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 196 real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] 197 real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] 198 real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] 199 real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 200 real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 201 real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 202 real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 203 real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] 204 real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] 205 real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 206 real, allocatable :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 207 real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 208 real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 209 logical :: bool_sublim ! logical to check if there is sublimation or not 210 211 ! Some variables for the PEM run 212 real, parameter :: year_step = 1 ! timestep for the pem 213 integer :: year_iter ! number of iteration 214 integer :: year_iter_max ! maximum number of iterations before stopping 215 real :: timestep ! timestep [s] 216 real :: watercap_sum ! total mass of water cap [kg/m^2] 217 real :: water_sum ! total mass of water in the mesh [kg/m^2] 218 219 #ifdef CPP_STD 220 real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice 221 real :: albedo_h2o_frost ! albedo of h2o frost 222 real, allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 223 real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 224 real, allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 225 real, allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 226 real, allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 227 real, allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case 228 real, allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case 229 real, allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case 230 real, allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case 231 real, allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model 232 logical, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model 233 real, allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case 234 real, allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case 235 #endif 236 237 #ifdef CPP_1D 238 integer :: nsplit_phys 239 integer, parameter :: jjm_value = jjm - 1 240 integer :: ierr 241 #else 242 integer, parameter :: jjm_value = jjm 243 #endif 244 245 ! Loop variables 246 integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil 247 248 ! Parallel variables 70 249 #ifndef CPP_STD 71 use mod_phys_lmdz_para, only: is_parallel, is_sequential, & 72 is_mpi_root, is_omp_root, & 73 is_master 74 use planete_h, only: aphelie, periheli, year_day, peri_day, & 75 obliquit 250 is_sequential = .true. 251 is_parallel = .false. 252 is_mpi_root = .true. 253 is_omp_root = .true. 254 is_master = .true. 255 #endif 256 257 day_ini = 0 ! test 258 time_phys = 0. ! test 259 260 ! Some constants 261 ngrid = ngridmx 262 nlayer = llm 263 264 A = (1/m_co2 - 1/m_noco2) 265 B = 1/m_noco2 266 267 year_day = 669 268 daysec = 88775. 269 timestep = year_day*daysec/year_step 270 271 !----------------------------- INITIALIZATION -------------------------- 272 !------------------------ 273 ! I Initialization 274 ! I_a READ run.def 275 !------------------------ 276 #ifndef CPP_1D 277 dtphys = 0 278 call conf_gcm(99,.true.) 279 call infotrac_init 280 nq = nqtot 281 allocate(q(ip1jmp1,llm,nqtot)) 76 282 #else 77 use planete_mod, only: apoastr, periastr, year_day, peri_day, & 78 obliquit 79 #endif 80 81 USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & 82 subslope_dist,iflat, & 83 major_slope,ini_comslope_h 84 85 #ifndef CPP_STD 86 USE comcstfi_h, only: r, mugaz 87 #else 88 USE comcstfi_mod, only: r, mugaz 89 #endif 90 91 USE logic_mod, ONLY: iflag_phys 92 USE mod_const_mpi, ONLY: COMM_LMDZ 93 use time_phylmdz_mod, only: daysec,dtphys 94 USE comconst_mod, ONLY: rad,g,cpp,pi 95 USE infotrac 96 USE geometry_mod, only: latitude_deg 97 98 use conf_pem_mod, only: conf_pem 99 use pemredem, only: pemdem0,pemdem1 100 use glaciers_mod, only: co2glaciers_evol,h2oglaciers_evol,co2glaciersflow,h2oglaciersflow 101 use criterion_pem_stop_mod, only: criterion_waterice_stop,criterion_co2_stop 102 use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, & 103 m_co2,m_noco2,threshold_water_frost2perenial,threshold_co2_frost2perenial 104 use evol_co2_ice_s_mod, only: evol_co2_ice_s 105 use evol_h2o_ice_s_mod, only: evol_h2o_ice_s 106 use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, & 107 TI_PEM,inertiedat_PEM, & ! soil thermal inertia 108 tsoil_PEM, mlayer_PEM,layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths 109 fluxgeo, & ! geothermal flux for the PEM and GCM 110 water_reservoir ! Water ressources 111 use adsorption_mod, only : regolith_adsorption,adsorption_pem, & ! bool to check if adsorption, main subroutine 112 ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! allocate arrays 113 co2_adsorbded_phys, h2o_adsorbded_phys ! mass of co2 and h2O adsorbded 114 USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem 115 use orbit_param_criterion_mod, only : orbit_param_criterion 116 use recomp_orb_param_mod, only: recomp_orb_param 117 use ice_table_mod, only: porefillingice_depth,porefillingice_thickness,& 118 end_ice_table_porefilling,ini_ice_table_porefilling, & 119 computeice_table_equilibrium 120 use soil_thermalproperties_mod, only: update_soil_thermalproperties 121 122 #ifdef CPP_1D 123 use regular_lonlat_mod, only: init_regular_lonlat 124 use physics_distribution_mod, only: init_physics_distribution 125 use mod_grid_phy_lmdz, only : regular_lonlat 126 use init_phys_1d_mod, only : init_phys_1d 127 #endif 128 129 130 IMPLICIT NONE 131 132 include "dimensions.h" 133 include "paramet.h" 134 135 INTEGER ngridmx 136 PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm ) 137 138 include "comgeom.h" 139 include "iniprint.h" 140 ! Same variable's name as in the GCM 141 142 INTEGER :: ngrid !Number of physical grid points 143 INTEGER :: nlayer !Number of vertical layer 144 INTEGER :: nq !Number of tracer 145 INTEGER :: day_ini !First day of the simulation 146 REAL :: pday !Physical day 147 REAL :: time_phys !Same as GCM 148 REAL :: ptimestep !Same as GCM 149 REAL :: ztime_fin !Same as GCM 150 151 ! Variable for reading start.nc 152 character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM 153 ! variables dynamiques 154 REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants 155 REAL teta(ip1jmp1,llm) ! temperature potentielle 156 REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes 157 REAL ps(ip1jmp1) ! pression au sol 158 159 REAL, dimension(:),allocatable :: ps_start_GCM !(ngrid) pression au sol 160 REAL, dimension(:,:),allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées 161 162 REAL masse(ip1jmp1,llm) ! masse d'air 163 REAL phis(ip1jmp1) ! geopotentiel au sol 164 REAL time_0 165 166 ! Variable for reading starfi.nc 167 168 character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM 169 character*2 str2 170 integer :: ncid, varid,status !Variable for handling opening of files 171 integer :: phydimid, subdimid, nlayerdimid, nqdimid !Variable ID for Netcdf files 172 integer :: lonvarid, latvarid, areavarid,sdvarid !Variable ID for Netcdf files 173 integer :: apvarid,bpvarid !Variable ID for Netcdf files 174 175 ! Variable for reading starfi.nc and writting restartfi.nc 176 177 REAL, dimension(:),allocatable :: longitude !Longitude read in FILE_NAME and written in restartfi 178 REAL, dimension(:),allocatable :: latitude !Latitude read in FILE_NAME and written in restartfi 179 REAL, dimension(:),allocatable :: ap !Coefficient ap read in FILE_NAME_start and written in restart 180 REAL, dimension(:),allocatable :: bp !Coefficient bp read in FILE_NAME_start and written in restart 181 REAL, dimension(:),allocatable :: cell_area !Cell_area read in FILE_NAME and written in restartfi 182 REAL :: Total_surface !Total surface of the planet 183 184 ! Variable for h2o_ice evolution 185 REAL :: ini_surf_h2o ! Initial surface of sublimating h2o ice 186 REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice 187 188 REAL :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa] 189 REAL :: global_ave_press_old ! constant: Global average pressure of initial/previous time step 190 REAL :: global_ave_press_new ! constant: Global average pressure of current time step 191 192 REAL , dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2] 193 REAL , dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2] 194 REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] 195 REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step 196 197 LOGICAL :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? 198 LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate? 199 LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? 200 LOGICAL :: STOPPING_pressure! Logical : is the criterion (% of change in the surface pressure) reached? 201 INTEGER :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param 202 203 real,save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 204 real ,allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] 205 real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM 206 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] 207 REAL ,allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg] 208 integer :: timelen ! # time samples 209 REAL :: ave ! intermediate varibale to compute average 210 211 REAL, ALLOCATABLE :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 212 REAL :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 213 214 215 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE 216 REAL, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2] 217 REAL, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2] 218 REAL, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2] 219 REAL, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2] 220 REAL, dimension(:,:,:),allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2] 221 REAL, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice 222 REAL, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state 223 REAL, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state 224 REAL, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year 225 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 226 REAL, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial h2o ice 227 REAL, dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope): Flag where there is a CO2 glacier flow 228 REAL, dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) : Flag where there is a CO2 glacier flow 229 REAL, dimension(:,:), allocatable :: flag_h2oflow(:,:) !(ngrid,nslope): Flag where there is a H2O glacier flow 230 REAL, dimension(:), allocatable :: flag_h2oflow_mesh(:) !(ngrid) : Flag where there is a H2O glacier flow 231 232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL 233 234 REAL, ALLOCATABLE :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K] 235 REAL, ALLOCATABLE :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K] 236 REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K] 237 238 REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 239 REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K] 240 REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K] 241 242 REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI] 243 REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K] 244 REAL,ALLOCATABLE :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] 245 REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] 246 REAL,ALLOCATABLE :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3] 247 REAL,ALLOCATABLE :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] 248 REAL,ALLOCATABLE :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3] 249 REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K] 250 REAL, ALLOCATABLE :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 251 REAL, ALLOCATABLE :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 252 REAL :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 253 REAL :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 254 LOGICAL :: bool_sublim ! logical to check if there is sublimation or not 255 256 257 !! Some parameters for the PEM run 258 REAL, PARAMETER :: year_step = 1 ! timestep for the pem 259 INTEGER :: year_iter ! number of iteration 260 INTEGER :: year_iter_max ! maximum number of iterations before stopping 261 REAL :: timestep ! timestep [s] 262 REAL :: watercap_sum ! total mass of water cap [kg/m^2] 263 REAL :: water_sum ! total mass of water in the mesh [kg/m^2] 264 265 #ifdef CPP_STD 266 REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice 267 REAL :: albedo_h2o_frost ! albedo of h2o frost 268 REAL,allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 269 REAL,allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 270 REAL,allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 271 REAL,allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 272 REAL,allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic 273 REAL,allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case 274 REAL,allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case 275 REAL,allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case 276 REAL,allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case 277 REAL,allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model 278 LOGICAL, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model 279 REAL,allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case 280 REAL,allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case 281 #endif 282 283 #ifdef CPP_1D 284 INTEGER :: nsplit_phys 285 integer,parameter:: jjm_value=jjm-1 286 integer :: ierr 287 #else 288 integer,parameter:: jjm_value=jjm 289 #endif 290 291 !!!!!!!!!!!!!!!!!!!!!!!!!!!! 292 293 ! Loop variable 294 INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil 295 #ifndef CPP_STD 296 ! Parallel variables 297 is_sequential=.true. 298 is_parallel=.false. 299 is_mpi_root=.true. 300 is_omp_root=.true. 301 is_master=.true. 302 #endif 303 304 day_ini=0 !test 305 time_phys=0. !test 306 307 ! Some constants 308 309 ngrid=ngridmx 310 nlayer=llm 311 312 A =(1/m_co2 - 1/m_noco2) 313 B=1/m_noco2 314 315 year_day=669 316 daysec=88775. 317 timestep=year_day*daysec/year_step 318 319 !------------------------ 320 321 ! I Initialisation 322 ! I_a READ run.def 323 324 !------------------------ 325 326 !----------------------------READ run.def --------------------- 327 328 #ifndef CPP_1D 329 dtphys=0 330 CALL conf_gcm( 99, .TRUE. ) 331 call infotrac_init 332 nq=nqtot 333 allocate(q(ip1jmp1,llm,nqtot)) 334 #else 335 ! load tracer names from file 'traceur.def' 336 open(90,file='traceur.def',status='old',form='formatted',& 337 iostat=ierr) 338 if (ierr.ne.0) then 339 write(*,*) 'Cannot find required file "traceur.def"' 340 write(*,*) ' If you want to run with tracers, I need it' 341 write(*,*) ' ... might as well stop here ...' 342 stop 343 else 344 write(*,*) "pem1d: Reading file traceur.def" 345 ! read number of tracers: 346 read(90,*,iostat=ierr) nq 347 print *, "nq",nq 348 nqtot=nq ! set value of nqtot (in infotrac module) as nq 349 if (ierr.ne.0) then 283 ! load tracer names from file 'traceur.def' 284 open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr) 285 if (ierr /= 0) then 286 write(*,*) 'Cannot find required file "traceur.def"' 287 write(*,*) ' If you want to run with tracers, I need it' 288 write(*,*) ' ... might as well stop here ...' 289 stop 290 else 291 write(*,*) "pem1d: Reading file traceur.def" 292 ! read number of tracers: 293 read(90,*,iostat = ierr) nq 294 write(*,*) "nq",nq 295 nqtot = nq ! set value of nqtot (in infotrac module) as nq 296 if (ierr /= 0) then 350 297 write(*,*) "pem1d: error reading number of tracers" 351 298 write(*,*) " (first line of traceur.def) " 352 299 stop 353 354 if (nq<1) then300 endif 301 if (nq < 1) then 355 302 write(*,*) "pem1d: error number of tracers" 356 303 write(*,*) "is nq=",nq," but must be >=1!" 357 304 stop 358 endif359 305 endif 360 nq=nqtot361 allocate(q(ip1jmp1,llm,nqtot))362 allocate(ap(nlayer+1))363 allocate(bp(nlayer+1))364 call init_phys_1d(llm,nqtot,vcov,ucov, &365 teta,q,ps, time_0,ap,bp)366 pi=2.E+0*asin(1.E+0)367 g=3.72368 nsplit_phys=1306 endif 307 nq = nqtot 308 allocate(q(ip1jmp1,llm,nqtot)) 309 allocate(ap(nlayer + 1)) 310 allocate(bp(nlayer + 1)) 311 call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp) 312 pi = 2.*asin(1.) 313 g = 3.72 314 nsplit_phys = 1 369 315 #endif 370 CALL conf_pem 371 372 !------------------------ 373 374 ! I Initialisation 375 ! I_a READ run.def 316 317 call conf_pem 318 319 !------------------------ 320 ! I Initialization 376 321 ! I_b READ of start_evol.nc and starfi_evol.nc 377 378 !------------------------ 379 380 !----------------------------Initialisation : READ some constant of startfi_evol.nc --------------------- 381 382 !----------------------------READ start.nc --------------------- 383 384 322 !------------------------ 323 ! I_b.1 READ start_evol.nc 324 allocate(ps_start_GCM(ngrid)) 385 325 #ifndef CPP_1D 386 CALL dynetat0(FILE_NAME_start,vcov,ucov, & 387 teta,q,masse,ps,phis, time_0) 388 #endif 389 390 allocate(ps_start_GCM(ngrid)) 391 #ifndef CPP_1D 392 call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) 393 394 CALL iniconst !new 395 CALL inigeom 396 397 allocate(ap(nlayer+1)) 398 allocate(bp(nlayer+1)) 399 status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid) 400 status = nf90_inq_varid(ncid, "ap", apvarid) 401 status = nf90_get_var(ncid, apvarid, ap) 402 status = nf90_inq_varid(ncid, "bp", bpvarid) 403 status = nf90_get_var(ncid, bpvarid, bp) 404 status =nf90_close(ncid) 405 326 call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0) 327 328 call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM) 329 330 call iniconst !new 331 call inigeom 332 333 allocate(ap(nlayer + 1)) 334 allocate(bp(nlayer + 1)) 335 status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid) 336 status = nf90_inq_varid(ncid,"ap",apvarid) 337 status = nf90_get_var(ncid,apvarid,ap) 338 status = nf90_inq_varid(ncid,"bp",bpvarid) 339 status = nf90_get_var(ncid,bpvarid,bp) 340 status = nf90_close(ncid) 341 342 call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, & 343 rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) 406 344 #else 407 408 ps_start_GCM(1)=ps(1) 409 410 #endif 411 412 413 #ifndef CPP_1D 414 CALL iniphysiq(iim,jjm,llm, & 415 (jjm-1)*iim+2,comm_lmdz, & 416 daysec,day_ini,dtphys/nsplit_phys, & 417 rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, & 418 iflag_phys) 345 ps_start_GCM(1) = ps(1) 419 346 #endif 420 347 421 348 ! In the gcm, these values are given to the physic by the dynamic. 422 349 ! Here we simply read them in the startfi_evol.nc file 423 status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) 424 425 allocate(longitude(ngrid)) 426 allocate(latitude(ngrid)) 427 allocate(cell_area(ngrid)) 428 429 status = nf90_inq_varid(ncid, "longitude", lonvarid) 430 status = nf90_get_var(ncid, lonvarid, longitude) 431 432 status = nf90_inq_varid(ncid, "latitude", latvarid) 433 status = nf90_get_var(ncid, latvarid, latitude) 434 435 status = nf90_inq_varid(ncid, "area", areavarid) 436 status = nf90_get_var(ncid, areavarid, cell_area) 437 438 status = nf90_inq_varid(ncid, "soildepth", sdvarid) 439 status = nf90_get_var(ncid, sdvarid, mlayer) 440 441 status =nf90_close(ncid) 442 443 !----------------------------READ startfi.nc --------------------- 444 350 status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid) 351 352 allocate(longitude(ngrid)) 353 allocate(latitude(ngrid)) 354 allocate(cell_area(ngrid)) 355 356 status = nf90_inq_varid(ncid,"longitude",lonvarid) 357 status = nf90_get_var(ncid,lonvarid,longitude) 358 359 status = nf90_inq_varid(ncid,"latitude",latvarid) 360 status = nf90_get_var(ncid,latvarid,latitude) 361 362 status = nf90_inq_varid(ncid,"area",areavarid) 363 status = nf90_get_var(ncid,areavarid,cell_area) 364 365 status = nf90_inq_varid(ncid,"soildepth",sdvarid) 366 status = nf90_get_var(ncid,sdvarid,mlayer) 367 368 status = nf90_close(ncid) 369 370 ! I_b.2 READ startfi_evol.nc 445 371 ! First we read the initial state (starfi.nc) 446 447 372 #ifndef CPP_STD 448 CALL phyetat0 (FILE_NAME,0,0, & 449 nsoilmx,ngrid,nlayer,nq, & 450 day_ini,time_phys, & 451 tsurf,tsoil,albedo,emis, & 452 q2,qsurf,tauscaling,totcloudfrac,wstar, & 453 watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist) 373 call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, & 374 tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar, & 375 watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist) 454 376 455 377 ! Remove unphysical values of surface tracer 456 DO i=1,ngrid 457 DO nnq=1,nqtot 458 DO islope=1,nslope 459 if(qsurf(i,nnq,islope).LT.0) then 460 qsurf(i,nnq,islope)=0. 461 endif 462 enddo 463 enddo 464 enddo 465 466 call surfini(ngrid,qsurf) 467 378 do i = 1,ngrid 379 do nnq = 1,nqtot 380 do islope = 1,nslope 381 if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0. 382 enddo 383 enddo 384 enddo 385 386 call surfini(ngrid,qsurf) 468 387 #else 469 call phys_state_var_init(nq) 470 IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) 471 call initracer(ngrid,nq) 472 call iniaerosol() 473 allocate(tsurf_read_generic(ngrid)) 474 allocate(qsurf_read_generic(ngrid,nq)) 475 allocate(tsoil_read_generic(ngrid,nsoilmx)) 476 allocate(emis_read_generic(ngrid)) 477 allocate(tsurf(ngrid,1)) 478 allocate(qsurf(ngrid,nq,1)) 479 allocate(tsoil(ngrid,nsoilmx,1)) 480 allocate(emis(ngrid,1)) 481 allocate(watercap(ngrid,1)) 482 allocate(watercaptag(ngrid)) 483 allocate(albedo_read_generic(ngrid,2)) 484 allocate(albedo(ngrid,2,1)) 485 allocate(inertiesoil(ngrid,nsoilmx,1)) 486 call phyetat0(.true., & 487 ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, & 488 day_ini,time_phys,tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,qsurf_read_generic, & 489 cloudfrac,totcloudfrac,hice, & 490 rnat,pctsrf_sic,tslab, tsea_ice,sea_ice) 491 call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 492 493 nslope=1 494 call ini_comslope_h(ngrid,1) 495 496 qsurf(:,:,1)=qsurf_read_generic(:,:) 497 tsurf(:,1)=tsurf_read_generic(:) 498 tsoil(:,:,1)=tsoil_read_generic(:,:) 499 emis(:,1)=emis_read_generic(:) 500 watercap(:,1)=0. 501 watercaptag(:)=.false. 502 albedo(:,1,1)=albedo_read_generic(:,1) 503 albedo(:,2,1)=albedo_read_generic(:,2) 504 inertiesoil(:,:,1)=inertiedat(:,:) 505 506 if (nslope.eq.1) then 507 def_slope(1) = 0 508 def_slope(2) = 0 509 def_slope_mean=0 510 subslope_dist(:,1) = 1. 511 endif 388 call phys_state_var_init(nq) 389 if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) 390 call initracer(ngrid,nq) 391 call iniaerosol() 392 allocate(tsurf_read_generic(ngrid)) 393 allocate(qsurf_read_generic(ngrid,nq)) 394 allocate(tsoil_read_generic(ngrid,nsoilmx)) 395 allocate(emis_read_generic(ngrid)) 396 allocate(tsurf(ngrid,1)) 397 allocate(qsurf(ngrid,nq,1)) 398 allocate(tsoil(ngrid,nsoilmx,1)) 399 allocate(emis(ngrid,1)) 400 allocate(watercap(ngrid,1)) 401 allocate(watercaptag(ngrid)) 402 allocate(albedo_read_generic(ngrid,2)) 403 allocate(albedo(ngrid,2,1)) 404 allocate(inertiesoil(ngrid,nsoilmx,1)) 405 call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, & 406 tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & 407 qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & 408 tslab, tsea_ice,sea_ice) 409 call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 410 411 nslope = 1 412 call ini_comslope_h(ngrid,1) 413 414 qsurf(:,:,1) = qsurf_read_generic(:,:) 415 tsurf(:,1) = tsurf_read_generic(:) 416 tsoil(:,:,1) = tsoil_read_generic(:,:) 417 emis(:,1) = emis_read_generic(:) 418 watercap(:,1) = 0. 419 watercaptag(:) = .false. 420 albedo(:,1,1) = albedo_read_generic(:,1) 421 albedo(:,2,1) = albedo_read_generic(:,2) 422 inertiesoil(:,:,1) = inertiedat(:,:) 423 424 if (nslope == 1) then 425 def_slope(1) = 0 426 def_slope(2) = 0 427 def_slope_mean = 0 428 subslope_dist(:,1) = 1. 429 endif 512 430 513 431 ! Remove unphysical values of surface tracer 514 DO i=1,ngrid 515 DO nnq=1,nqtot 516 qsurf(i,nnq,1)=qsurf_read_generic(i,nnq) 517 if(qsurf(i,nnq,1).LT.0) then 518 qsurf(i,nnq,1)=0. 519 endif 520 enddo 521 enddo 432 do i = 1,ngrid 433 do nnq = 1,nqtot 434 qsurf(i,nnq,1) = qsurf_read_generic(i,nnq) 435 if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0. 436 enddo 437 enddo 522 438 #endif 523 439 524 DO nnq=1,nqtot ! Why not using ini_tracer ? 525 if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq 526 if(noms(nnq).eq."h2o_vap") then 527 igcm_h2o_vap = nnq 528 mmol(igcm_h2o_vap)=18. 529 endif 530 if(noms(nnq).eq."co2") igcm_co2 = nnq 531 ENDDO 532 r= 8.314511E+0 *1000.E+0/mugaz 533 !------------------------ 534 535 ! I Initialisation 536 ! I_a READ run.def 537 ! I_b READ of start_evol.nc and starfi_evol.nc 440 do nnq = 1,nqtot ! Why not using ini_tracer ? 441 if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq 442 if (noms(nnq) == "h2o_vap") then 443 igcm_h2o_vap = nnq 444 mmol(igcm_h2o_vap)=18. 445 endif 446 if (noms(nnq) == "co2") igcm_co2 = nnq 447 enddo 448 r = 8.314511E+0*1000.E+0/mugaz 449 450 !------------------------ 451 ! I Initialization 538 452 ! I_c Subslope parametrisation 539 540 !------------------------ 541 542 !----------------------------Subslope parametrisation definition --------------------- 543 544 ! Define some slope statistics 545 iflat=1 546 DO islope=2,nslope 547 IF(abs(def_slope_mean(islope)).lt. & 548 abs(def_slope_mean(iflat))) THEN 549 iflat = islope 550 ENDIF 551 ENDDO 552 553 PRINT*,'Flat slope for islope = ',iflat 554 PRINT*,'corresponding criterium = ',def_slope_mean(iflat) 555 556 allocate(flag_co2flow(ngrid,nslope)) 557 allocate(flag_co2flow_mesh(ngrid)) 558 allocate(flag_h2oflow(ngrid,nslope)) 559 allocate(flag_h2oflow_mesh(ngrid)) 560 561 flag_co2flow(:,:) = 0 562 flag_co2flow_mesh(:) = 0 563 flag_h2oflow(:,:) = 0 564 flag_h2oflow_mesh(:) = 0 565 !---------------------------- READ GCM data --------------------- 566 567 ! I Initialisation 568 ! I_a READ run.def 569 ! I_b READ of start_evol.nc and starfi_evol.nc 570 ! I_c Subslope parametrisation 571 ! I_d READ GCM data 572 573 !------------------------ 574 453 !------------------------ 454 ! Define some slope statistics 455 iflat = 1 456 do islope = 2,nslope 457 if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope 458 enddo 459 460 write(*,*) 'Flat slope for islope = ',iflat 461 write(*,*) 'corresponding criterium = ',def_slope_mean(iflat) 462 463 allocate(flag_co2flow(ngrid,nslope)) 464 allocate(flag_co2flow_mesh(ngrid)) 465 allocate(flag_h2oflow(ngrid,nslope)) 466 allocate(flag_h2oflow_mesh(ngrid)) 467 468 flag_co2flow(:,:) = 0 469 flag_co2flow_mesh(:) = 0 470 flag_h2oflow(:,:) = 0 471 flag_h2oflow_mesh(:) = 0 472 473 !------------------------ 474 ! I Initialization 475 ! I_d READ GCM data and convert to the physical grid 476 !------------------------ 575 477 ! 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 576 577 call nb_time_step_GCM("data_GCM_Y1.nc",timelen) 578 579 allocate(tsoil_ave(ngrid,nsoilmx,nslope)) 580 allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope)) 581 582 allocate(vmr_co2_gcm(ngrid,timelen)) 583 allocate(ps_timeseries(ngrid,timelen)) 584 allocate(min_co2_ice_1(ngrid,nslope)) 585 allocate(min_h2o_ice_1(ngrid,nslope)) 586 allocate(min_co2_ice_2(ngrid,nslope)) 587 allocate(min_h2o_ice_2(ngrid,nslope)) 588 allocate(tsurf_ave_yr1(ngrid,nslope)) 589 allocate(tsurf_ave(ngrid,nslope)) 590 591 allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen)) 592 allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen)) 593 allocate(q_co2_PEM_phys(ngrid,timelen)) 594 allocate(q_h2o_PEM_phys(ngrid,timelen)) 595 allocate(co2_ice_GCM(ngrid,nslope,timelen)) 596 allocate(watersurf_density_ave(ngrid,nslope)) 597 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) 598 599 allocate(Tsurfave_before_saved(ngrid,nslope)) 600 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 601 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 602 allocate(delta_co2_adsorbded(ngrid)) 603 allocate(delta_h2o_adsorbded(ngrid)) 604 allocate(vmr_co2_pem_phys(ngrid,timelen)) 605 606 print *, "Downloading data Y1..." 607 608 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,& 609 tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & 610 watersurf_density_ave,watersoil_density_timeseries) 478 call nb_time_step_GCM("data_GCM_Y1.nc",timelen) 479 480 allocate(tsoil_ave(ngrid,nsoilmx,nslope)) 481 allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope)) 482 allocate(vmr_co2_gcm(ngrid,timelen)) 483 allocate(ps_timeseries(ngrid,timelen)) 484 allocate(min_co2_ice_1(ngrid,nslope)) 485 allocate(min_h2o_ice_1(ngrid,nslope)) 486 allocate(min_co2_ice_2(ngrid,nslope)) 487 allocate(min_h2o_ice_2(ngrid,nslope)) 488 allocate(tsurf_ave_yr1(ngrid,nslope)) 489 allocate(tsurf_ave(ngrid,nslope)) 490 allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen)) 491 allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen)) 492 allocate(q_co2_PEM_phys(ngrid,timelen)) 493 allocate(q_h2o_PEM_phys(ngrid,timelen)) 494 allocate(co2_ice_GCM(ngrid,nslope,timelen)) 495 allocate(watersurf_density_ave(ngrid,nslope)) 496 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) 497 allocate(Tsurfave_before_saved(ngrid,nslope)) 498 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 499 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 500 allocate(delta_co2_adsorbded(ngrid)) 501 allocate(delta_h2o_adsorbded(ngrid)) 502 allocate(vmr_co2_pem_phys(ngrid,timelen)) 503 504 write(*,*) "Downloading data Y1..." 505 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, & 506 tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 507 co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries) 508 write(*,*) "Downloading data Y1 done" 611 509 612 510 ! 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 613 614 print *, "Downloading data Y1 done" 615 print *, "Downloading data Y2" 616 617 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, & 618 tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, & 619 watersurf_density_ave,watersoil_density_timeseries) 620 621 print *, "Downloading data Y2 done" 622 623 !------------------------ 624 625 ! I Initialisation 626 ! I_a READ run.def 627 ! I_b READ of start_evol.nc and starfi_evol.nc 628 ! I_c Subslope parametrisation 629 ! I_d READ GCM data and convert to the physical grid 630 ! I_e Initialisation of the PEM variable and soil 631 632 !------------------------ 633 634 !---------------------------- Initialisation of the PEM soil and values --------------------- 635 636 call end_comsoil_h_PEM 637 call ini_comsoil_h_PEM(ngrid,nslope) 638 call end_adsorption_h_PEM 639 call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) 640 call end_ice_table_porefilling 641 call ini_ice_table_porefilling(ngrid,nslope) 642 643 if(soil_pem) then 644 call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) 645 DO l=1,nsoilmx 511 write(*,*) "Downloading data Y2" 512 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, & 513 tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 514 co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries) 515 write(*,*) "Downloading data Y2 done" 516 517 !------------------------ 518 ! I Initialization 519 ! I_e Initialization of the PEM variables and soil 520 !------------------------ 521 call end_comsoil_h_PEM 522 call ini_comsoil_h_PEM(ngrid,nslope) 523 call end_adsorption_h_PEM 524 call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) 525 call end_ice_table_porefilling 526 call ini_ice_table_porefilling(ngrid,nslope) 527 528 if (soil_pem) then 529 call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) 530 do l=1,nsoilmx 646 531 tsoil_PEM(:,l,:)=tsoil_ave(:,l,:) 647 532 tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:) 648 533 watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:) 649 ENDDO650 DOl=nsoilmx+1,nsoilmx_PEM534 enddo 535 do l=nsoilmx+1,nsoilmx_PEM 651 536 tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:) 652 537 watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:) 653 ENDDO 654 watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen 655 endif !soil_pem 656 deallocate(tsoil_ave) 657 deallocate(tsoil_GCM_timeseries) 658 659 !------------------------ 660 661 ! I Initialisation 662 ! I_a READ run.def 663 ! I_b READ of start_evol.nc and starfi_evol.nc 664 ! I_c Subslope parametrisation 665 ! I_d READ GCM data and convert to the physical grid 666 ! I_e Initialisation of the PEM variable and soil 538 enddo 539 watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen 540 endif !soil_pem 541 deallocate(tsoil_ave,tsoil_GCM_timeseries) 542 543 !------------------------ 544 ! I Initialization 667 545 ! I_f Compute tendencies & Save initial situation 668 669 !----- Compute tendencies from the PCM run 670 671 allocate(tendencies_co2_ice(ngrid,nslope)) 672 allocate(tendencies_co2_ice_ini(ngrid,nslope)) 673 allocate(tendencies_h2o_ice(ngrid,nslope)) 674 675 ! Compute the tendencies of the evolution of ice over the years 676 677 call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,& 678 min_co2_ice_2,tendencies_co2_ice) 679 680 tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:) 681 682 call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,& 683 min_h2o_ice_2,tendencies_h2o_ice) 684 685 deallocate(min_co2_ice_1) 686 deallocate(min_co2_ice_2) 687 deallocate(min_h2o_ice_1) 688 deallocate(min_h2o_ice_2) 689 690 !------------------------ 691 692 ! I Initialisation 693 ! I_a READ run.def 694 ! I_b READ of start_evol.nc and starfi_evol.nc 695 ! I_c Subslope parametrisation 696 ! I_d READ GCM data and convert to the physical grid 697 ! I_e Initialisation of the PEM variable and soil 698 ! I_f Compute tendencies & Save initial situation 546 !------------------------ 547 allocate(tendencies_co2_ice(ngrid,nslope)) 548 allocate(tendencies_co2_ice_ini(ngrid,nslope)) 549 allocate(tendencies_h2o_ice(ngrid,nslope)) 550 551 ! Compute the tendencies of the evolution of ice over the years 552 call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice) 553 tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:) 554 call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice) 555 556 deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2) 557 558 !------------------------ 559 ! I Initialization 699 560 ! I_g Save initial PCM situation 700 701 !---------------------------- Save initial PCM situation --------------------- 702 703 allocate(initial_co2_ice_sublim(ngrid,nslope)) 704 allocate(initial_co2_ice(ngrid,nslope)) 705 allocate(initial_h2o_ice(ngrid,nslope)) 561 !------------------------ 562 allocate(initial_co2_ice_sublim(ngrid,nslope)) 563 allocate(initial_co2_ice(ngrid,nslope)) 564 allocate(initial_h2o_ice(ngrid,nslope)) 706 565 707 566 ! We save the places where water ice is sublimating 708 567 ! We compute the surface of water ice sublimating 709 ini_surf_co2=0. 710 ini_surf_h2o=0. 711 Total_surface=0. 712 do i=1,ngrid 713 Total_surface=Total_surface+cell_area(i) 714 do islope=1,nslope 715 if (tendencies_co2_ice(i,islope).LT.0) then 716 initial_co2_ice_sublim(i,islope)=1. 717 ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope) 718 else 719 initial_co2_ice_sublim(i,islope)=0. 720 endif 721 if (qsurf(i,igcm_co2,islope).GT.0) then 722 initial_co2_ice(i,islope)=1. 723 else 724 initial_co2_ice(i,islope)=0. 725 endif 726 if (tendencies_h2o_ice(i,islope).LT.0) then 727 initial_h2o_ice(i,islope)=1. 728 ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope) 729 else 730 initial_h2o_ice(i,islope)=0. 731 endif 732 enddo 733 enddo 734 735 print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2 736 print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o 737 print *, "Total surface of the planet=", Total_surface 738 739 allocate(zplev_gcm(ngrid,nlayer+1)) 740 741 DO l=1,nlayer+1 742 DO ig=1,ngrid 743 zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) 744 ENDDO 745 ENDDO 746 747 global_ave_press_old=0. 748 do i=1,ngrid 749 global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(i)/Total_surface 750 enddo 751 752 global_ave_press_GCM=global_ave_press_old 753 global_ave_press_new=global_ave_press_old 754 print *, "Initial global average pressure=", global_ave_press_GCM 755 756 !------------------------ 757 758 ! I Initialisation 759 ! I_a READ run.def 760 ! I_b READ of start_evol.nc and starfi_evol.nc 761 ! I_c Subslope parametrisation 762 ! I_d READ GCM data and convert to the physical grid 763 ! I_e Initialisation of the PEM variable and soil 764 ! I_f Compute tendencies & Save initial situation 765 ! I_g Save initial PCM situation 568 ini_surf_co2 = 0. 569 ini_surf_h2o = 0. 570 Total_surface = 0. 571 do i = 1,ngrid 572 Total_surface = Total_surface+cell_area(i) 573 do islope = 1,nslope 574 if (tendencies_co2_ice(i,islope) < 0) then 575 initial_co2_ice_sublim(i,islope) = 1. 576 ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope) 577 else 578 initial_co2_ice_sublim(i,islope) = 0. 579 endif 580 if (qsurf(i,igcm_co2,islope) > 0) then 581 initial_co2_ice(i,islope) = 1. 582 else 583 initial_co2_ice(i,islope) = 0. 584 endif 585 if (tendencies_h2o_ice(i,islope) < 0) then 586 initial_h2o_ice(i,islope) = 1. 587 ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope) 588 else 589 initial_h2o_ice(i,islope) = 0. 590 endif 591 enddo 592 enddo 593 594 write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2 595 write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o 596 write(*,*) "Total surface of the planet=", Total_surface 597 allocate(zplev_gcm(ngrid,nlayer + 1)) 598 599 do l = 1,nlayer + 1 600 do ig = 1,ngrid 601 zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) 602 enddo 603 enddo 604 605 global_ave_press_old = 0. 606 do i = 1,ngrid 607 global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface 608 enddo 609 610 global_ave_press_GCM = global_ave_press_old 611 global_ave_press_new = global_ave_press_old 612 write(*,*) "Initial global average pressure=", global_ave_press_GCM 613 614 !------------------------ 615 ! I Initialization 766 616 ! I_h Read the PEMstart 767 768 !---------------------------- Read the PEMstart --------------------- 769 770 call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth,porefillingice_thickness, & 771 tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,& 772 tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,& 773 watersurf_density_ave,watersoil_density_PEM_ave, & 774 co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) 617 !------------------------ 618 call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, & 619 porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, & 620 tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), & 621 qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, & 622 co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir) 623 624 do ig = 1,ngrid 625 do islope = 1,nslope 626 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.) 627 enddo 628 enddo 629 630 if (adsorption_pem) then 631 totmassco2_adsorbded = 0. 632 totmassh2o_adsorbded = 0. 633 do ig = 1,ngrid 634 do islope = 1, nslope 635 do l = 1,nsoilmx_PEM - 1 636 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 637 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 638 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 639 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 640 enddo 641 enddo 642 enddo 643 644 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded 645 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded 646 endif ! adsorption 647 deallocate(tsurf_ave_yr1) 648 649 !------------------------ 650 ! I Initialization 651 ! I_i Compute orbit criterion 652 !------------------------ 653 #ifndef CPP_STD 654 call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) 655 #else 656 call iniorbit(apoastr, periastr, year_day, peri_day,obliquit) 657 #endif 658 659 if (evol_orbit_pem) then 660 call orbit_param_criterion(year_iter_max) 661 else 662 year_iter_max = Max_iter_pem 663 endif 664 !-------------------------- END INITIALIZATION ------------------------- 665 666 !-------------------------------- RUN ---------------------------------- 667 !------------------------ 668 ! II Run 669 ! II_a Update pressure, ice and tracers 670 !------------------------ 671 year_iter = 0 672 do while (year_iter < year_iter_max) 673 674 ! II.a.1. Compute updated global pressure 675 write(*,*) "Recomputing the new pressure..." 676 677 do i = 1,ngrid 678 do islope = 1,nslope 679 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 680 enddo 681 enddo 682 write(*,*) 'Global average pressure old time step',global_ave_press_old 683 call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) 684 685 if (adsorption_pem) then 686 do i = 1,ngrid 687 global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface 688 enddo 689 write(*,*) 'Global average pressure old time step',global_ave_press_old 690 write(*,*) 'Global average pressure new time step',global_ave_press_new 691 endif 692 693 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) 694 allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen)) 695 write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..." 696 do l = 1,nlayer + 1 697 do ig = 1,ngrid 698 zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) 699 enddo 700 enddo 701 702 ! II.a.3. Surface pressure timeseries 703 write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..." 704 do ig = 1,ngrid 705 ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old 706 enddo 707 708 ! II.a.4. New pressure levels timeseries 709 allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen)) 710 write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..." 711 do l = 1,nlayer + 1 712 do ig = 1,ngrid 713 zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) 714 enddo 715 enddo 716 717 ! II.a.5. Tracers timeseries 718 write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..." 719 720 l = 1 721 do ig = 1,ngrid 722 do t = 1,timelen 723 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))/ & 724 (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) 725 if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30 726 if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1. 727 enddo 728 enddo 729 730 do ig = 1,ngrid 731 do t = 1, timelen 732 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))/ & 733 (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & 734 + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & 735 - (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ & 736 (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) 737 if (q_co2_PEM_phys(ig,t) < 0) then 738 q_co2_PEM_phys(ig,t) = 1.e-30 739 elseif (q_co2_PEM_phys(ig,t) > 1) then 740 q_co2_PEM_phys(ig,t) = 1. 741 endif 742 mmean=1/(A*q_co2_PEM_phys(ig,t) + B) 743 vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 744 enddo 745 enddo 746 747 deallocate(zplev_new_timeseries,zplev_old_timeseries) 748 749 !------------------------ 750 ! II Run 751 ! II_b Evolution of the ice 752 !------------------------ 753 write(*,*) "Evolution of h2o ice" 754 call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope) 755 756 write(*,*) "Evolution of co2 ice" 757 call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope) 758 759 do islope=1, nslope 760 write(str2(1:2),'(i2.2)') islope 761 call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope)) 762 call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope)) 763 call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) 764 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) 765 enddo 766 767 !------------------------ 768 ! II Run 769 ! II_c CO2 & H2O glaciers flows 770 !------------------------ 771 write(*,*) "CO2 glacier flows" 772 773 if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, & 774 global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh) 775 775 776 do ig = 1,ngrid 776 write(*,*) "H2O glacier flows" 777 778 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) 779 780 do islope = 1,nslope 781 write(str2(1:2),'(i2.2)') islope 782 call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope)) 783 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope)) 784 call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) 785 enddo 786 787 !------------------------ 788 ! II Run 789 ! II_d Update surface and soil temperatures 790 !------------------------ 791 ! II_d.1 Update Tsurf 792 write(*,*) "Updating the new Tsurf" 793 bool_sublim = .false. 794 Tsurfave_before_saved(:,:) = tsurf_ave(:,:) 795 do ig = 1,ngrid 777 796 do islope = 1,nslope 778 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.) 779 enddo 780 enddo 781 782 if(adsorption_pem) then 783 totmassco2_adsorbded = 0. 784 totmassh2o_adsorbded = 0. 785 do ig = 1,ngrid 786 do islope =1, nslope 787 do l = 1,nsoilmx_PEM - 1 788 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 789 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & 790 cell_area(ig) 791 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 792 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & 793 cell_area(ig) 794 enddo 795 enddo 796 enddo 797 798 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded 799 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded 800 endif ! adsorption 801 deallocate(tsurf_ave_yr1) 802 803 !------------------------ 804 805 ! I Initialisation 806 ! I_a READ run.def 807 ! I_b READ of start_evol.nc and starfi_evol.nc 808 ! I_c Subslope parametrisation 809 ! I_d READ GCM data and convert to the physical grid 810 ! I_e Initialisation of the PEM variable and soil 811 ! I_f Compute tendencies & Save initial situation 812 ! I_g Save initial PCM situation 813 ! I_h Read the PEMstar 814 ! I_i Compute orbit criterion 815 816 #ifndef CPP_STD 817 CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) 818 #else 819 CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit) 820 #endif 821 822 if (evol_orbit_pem) then 823 call orbit_param_criterion(year_iter_max) 824 else 825 year_iter_max = Max_iter_pem 826 endif 827 828 !--------------------------- END INITIALISATION --------------------- 829 830 !---------------------------- RUN --------------------- 831 832 !------------------------ 833 797 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 798 if (latitude_deg(ig) > 0) then 799 do ig_loop = ig,ngrid 800 do islope_loop = islope,iflat,-1 801 if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then 802 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) 803 bool_sublim = .true. 804 exit 805 endif 806 enddo 807 if (bool_sublim) exit 808 enddo 809 else 810 do ig_loop = ig,1,-1 811 do islope_loop = islope,iflat 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 endif 821 initial_co2_ice(ig,islope) = 0 822 if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then 823 albedo(ig,1,islope) = albedo_h2o_frost 824 albedo(ig,2,islope) = albedo_h2o_frost 825 else 826 albedo(ig,1,islope) = albedodat(ig) 827 albedo(ig,2,islope) = albedodat(ig) 828 emis(ig,islope) = emissiv 829 endif 830 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 831 ave = 0. 832 do t = 1,timelen 833 if (co2_ice_GCM(ig,islope,t) > 1.e-3) then 834 ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.)) 835 else 836 ave = ave + tsurf_GCM_timeseries(ig,islope,t) 837 endif 838 enddo 839 tsurf_ave(ig,islope) = ave/timelen 840 endif 841 enddo 842 enddo 843 844 do t = 1,timelen 845 tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:)) 846 enddo 847 ! for the start 848 do ig = 1,ngrid 849 do islope = 1,nslope 850 tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope)) 851 enddo 852 enddo 853 854 do islope = 1,nslope 855 write(str2(1:2),'(i2.2)') islope 856 call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) 857 enddo 858 859 if (soil_pem) then 860 861 ! II_d.2 Update soil temperature 862 allocate(TI_locslope(ngrid,nsoilmx_PEM)) 863 allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) 864 allocate(Tsurf_locslope(ngrid)) 865 write(*,*)"Updating soil temperature" 866 867 ! Soil averaged 868 do islope = 1,nslope 869 TI_locslope(:,:) = TI_PEM(:,:,islope) 870 do t = 1,timelen 871 Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) 872 Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t) 873 call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) 874 call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope) 875 tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) 876 do ig = 1,ngrid 877 do isoil = 1,nsoilmx_PEM 878 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) 879 if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1) 880 enddo 881 enddo 882 enddo 883 enddo 884 tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen 885 watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen 886 887 write(*,*) "Update of soil temperature done" 888 889 deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope) 890 write(*,*) "Compute ice table" 891 892 ! II_d.3 Update the ice table 893 call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness) 894 write(*,*) "Update soil propreties" 895 896 ! II_d.4 Update the soil thermal properties 897 call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, & 898 porefillingice_depth,porefillingice_thickness,TI_PEM) 899 900 ! II_d.5 Update the mass of the regolith adsorbded 901 if (adsorption_pem) then 902 call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, & 903 qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, & 904 q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) 905 906 totmassco2_adsorbded = 0. 907 totmassh2o_adsorbded = 0. 908 do ig = 1,ngrid 909 do islope =1, nslope 910 do l = 1,nsoilmx_PEM - 1 911 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 912 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & 913 cell_area(ig) 914 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* & 915 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * & 916 cell_area(ig) 917 enddo 918 enddo 919 enddo 920 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded 921 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded 922 endif 923 endif !soil_pem 924 925 !------------------------ 834 926 ! II Run 835 ! II_a update pressure,ice and tracers 836 837 !------------------------ 838 year_iter=0 839 do while (year_iter.LT.year_iter_max) 840 841 ! II.a.1. Compute updated global pressure 842 print *, "Recomputing the new pressure..." 843 844 do i=1,ngrid 845 do islope=1,nslope 846 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 847 enddo 848 enddo 849 print *, 'Global average pressure old time step',global_ave_press_old 850 call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new) 851 852 if(adsorption_pem) then 853 do i=1,ngrid 854 global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface 855 enddo 856 print *, 'Global average pressure old time step',global_ave_press_old 857 print *, 'Global average pressure new time step',global_ave_press_new 858 endif 859 860 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) 861 allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen)) 862 print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..." 863 DO l=1,nlayer+1 864 DO ig=1,ngrid 865 zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) 866 ENDDO 867 ENDDO 868 869 ! II.a.3. Surface pressure timeseries 870 print *, "Recomputing the surface pressure timeserie adapted to the new pressure..." 871 do ig = 1,ngrid 872 ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old 873 enddo 874 875 ! II.a.4. New pressure levels timeseries 876 allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen)) 877 print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..." 878 do l=1,nlayer+1 879 do ig=1,ngrid 880 zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) 881 enddo 882 enddo 883 884 ! II.a.5. Tracers timeseries 885 print *, "Recomputing of tracer VMR timeseries for the new pressure..." 886 887 l=1 888 DO ig=1,ngrid 889 DO t = 1, timelen 890 q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) 891 if(q_h2o_PEM_phys(ig,t).LT.0) then 892 q_h2o_PEM_phys(ig,t)=1E-30 893 endif 894 if(q_h2o_PEM_phys(ig,t).GT.1) then 895 q_h2o_PEM_phys(ig,t)=1. 896 endif 897 enddo 898 enddo 899 900 DO ig=1,ngrid 901 DO t = 1, timelen 902 q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) & 903 + ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - & 904 (zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / & 905 (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) 906 if (q_co2_PEM_phys(ig,t).LT.0) then 907 q_co2_PEM_phys(ig,t)=1E-30 908 elseif (q_co2_PEM_phys(ig,t).GT.1) then 909 q_co2_PEM_phys(ig,t)=1. 910 endif 911 mmean=1/(A*q_co2_PEM_phys(ig,t) +B) 912 vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 913 ENDDO 914 ENDDO 915 916 deallocate(zplev_new_timeseries) 917 deallocate(zplev_old_timeseries) 918 927 ! II_e Update the tendencies 928 !------------------------ 929 write(*,*) "Adaptation of the new co2 tendencies to the current pressure" 930 call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, & 931 global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) 932 933 !------------------------ 919 934 ! II Run 920 ! II_a update pressure, ice and tracers921 ! II_b Evolution of the ice922 923 ! II.b. Evolution of the ice924 print *, "Evolution of h2o ice"925 call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)926 927 print *, "Evolution of co2 ice"928 call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)929 930 DO islope=1, nslope931 write(str2(1:2),'(i2.2)') islope932 call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))933 call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))934 call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))935 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))936 ENDDO937 938 !------------------------939 940 ! II Run941 ! II_a update pressure, ice and tracers942 ! II_b Evolution of the ice943 ! II_c CO2 & H2O glaciers flows944 945 !------------------------946 947 print *, "CO2 glacier flows"948 949 if (co2glaciersflow) then950 call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&951 global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)952 endif953 954 print *, "H2O glacier flows"955 956 if (h2oglaciersflow) then957 call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)958 endif959 DO islope=1, nslope960 write(str2(1:2),'(i2.2)') islope961 call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))962 call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))963 call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))964 ENDDO965 966 !------------------------967 968 ! II Run969 ! II_a update pressure, ice and tracers970 ! II_b Evolution of the ice971 ! II_c CO2 glaciers flows972 ! II_d Update surface and soil temperatures973 974 !------------------------975 976 ! II_d.1 Update Tsurf977 978 print *, "Updating the new Tsurf"979 bool_sublim=.false.980 Tsurfave_before_saved(:,:) = tsurf_ave(:,:)981 DO ig = 1,ngrid982 DO islope = 1,nslope983 if(initial_co2_ice(ig,islope).gt.0.5 .and. qsurf(ig,igcm_co2,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice984 if(latitude_deg(ig).gt.0) then985 do ig_loop=ig,ngrid986 DO islope_loop=islope,iflat,-1987 if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then988 tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)989 bool_sublim=.true.990 exit991 endif992 enddo993 if (bool_sublim.eqv. .true.) then994 exit995 endif996 enddo997 else998 do ig_loop=ig,1,-1999 DO islope_loop=islope,iflat1000 if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then1001 tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)1002 bool_sublim=.true.1003 exit1004 endif1005 enddo1006 if (bool_sublim.eqv. .true.) then1007 exit1008 endif1009 enddo1010 endif1011 initial_co2_ice(ig,islope)=01012 if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then1013 albedo(ig,1,islope) = albedo_h2o_frost1014 albedo(ig,2,islope) = albedo_h2o_frost1015 else1016 albedo(ig,1,islope) = albedodat(ig)1017 albedo(ig,2,islope) = albedodat(ig)1018 emis(ig,islope) = emissiv1019 endif1020 elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co21021 ave=0.1022 do t=1,timelen1023 if(co2_ice_GCM(ig,islope,t).gt.1e-3) then1024 ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))1025 else1026 ave = ave + tsurf_GCM_timeseries(ig,islope,t)1027 endif1028 enddo1029 tsurf_ave(ig,islope)=ave/timelen1030 endif1031 enddo1032 enddo1033 1034 do t = 1,timelen1035 tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:))1036 enddo1037 ! for the start1038 do ig = 1,ngrid1039 do islope = 1,nslope1040 tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))1041 enddo1042 enddo1043 1044 DO islope=1, nslope1045 write(str2(1:2),'(i2.2)') islope1046 call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))1047 ENDDO1048 1049 if(soil_pem) then1050 1051 ! II_d.2 Update soil temperature1052 1053 allocate(TI_locslope(ngrid,nsoilmx_PEM))1054 allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))1055 allocate(Tsurf_locslope(ngrid))1056 print *,"Updating soil temperature"1057 1058 ! Soil averaged1059 do islope = 1,nslope1060 TI_locslope(:,:) = TI_PEM(:,:,islope)1061 do t = 1,timelen1062 Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)1063 Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)1064 call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)1065 call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)1066 tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)1067 do ig = 1,ngrid1068 do isoil = 1,nsoilmx_PEM1069 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)1070 if(isnan(Tsoil_locslope(ig,isoil))) then1071 call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)1072 endif1073 enddo1074 enddo1075 enddo1076 enddo1077 tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen1078 watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen1079 1080 print *, "Update of soil temperature done"1081 1082 deallocate(TI_locslope)1083 deallocate(Tsoil_locslope)1084 deallocate(Tsurf_locslope)1085 write(*,*) "Compute ice table"1086 1087 ! II_d.3 Update the ice table1088 call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)1089 print *, "Update soil propreties"1090 1091 ! II_d.4 Update the soil thermal properties1092 call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &1093 porefillingice_depth,porefillingice_thickness,TI_PEM)1094 1095 ! II_d.5 Update the mass of the regolith adsorbded1096 if(adsorption_pem) then1097 call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), &1098 tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &1099 h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)1100 1101 1102 totmassco2_adsorbded = 0.1103 totmassh2o_adsorbded = 0.1104 do ig = 1,ngrid1105 do islope =1, nslope1106 do l = 1,nsoilmx_PEM - 11107 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &1108 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &1109 cell_area(ig)1110 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &1111 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &1112 cell_area(ig)1113 enddo1114 enddo1115 enddo1116 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded1117 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded1118 endif1119 endif !soil_pem1120 1121 !------------------------1122 1123 ! II Run1124 ! II_a update pressure, ice and tracers1125 ! II_b Evolution of the ice1126 ! II_c CO2 glaciers flows1127 ! II_d Update surface and soil temperatures1128 ! II_e Update the tendencies1129 1130 !------------------------1131 1132 print *, "Adaptation of the new co2 tendencies to the current pressure"1133 call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&1134 global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)1135 1136 !------------------------1137 1138 ! II Run1139 ! II_a update pressure, ice and tracers1140 ! II_b Evolution of the ice1141 ! II_c CO2 glaciers flows1142 ! II_d Update surface and soil temperatures1143 ! II_e Update the tendencies1144 935 ! II_f Checking the stopping criterion 1145 1146 !------------------------ 1147 call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) 1148 1149 call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & 1150 initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope) 1151 1152 year_iter=year_iter+dt_pem 1153 1154 print *, "Checking all the stopping criterion." 1155 if (STOPPING_water) then 1156 print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 1157 criterion_stop=1 1158 endif 1159 if (STOPPING_1_water) then 1160 print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water 1161 criterion_stop=1 1162 endif 1163 if (STOPPING_co2) then 1164 print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 1165 criterion_stop=2 1166 endif 1167 if (STOPPING_pressure) then 1168 print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure 1169 criterion_stop=3 1170 endif 1171 if (year_iter.ge.year_iter_max) then 1172 print *, "STOPPING because maximum number of iterations reached" 1173 criterion_stop=4 1174 endif 1175 1176 if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then 936 !------------------------ 937 call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice) 938 939 call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, & 940 initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope) 941 942 year_iter = year_iter + dt_pem 943 944 write(*,*) "Checking all the stopping criterion." 945 if (STOPPING_water) then 946 write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 947 criterion_stop = 1 948 endif 949 if (STOPPING_1_water) then 950 write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water 951 criterion_stop = 1 952 endif 953 if (STOPPING_co2) then 954 write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 955 criterion_stop = 2 956 endif 957 if (STOPPING_pressure) then 958 write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure 959 criterion_stop = 3 960 endif 961 if (year_iter >= year_iter_max) then 962 write(*,*) "STOPPING because maximum number of iterations reached" 963 criterion_stop = 4 964 endif 965 966 if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then 1177 967 exit 1178 else 1179 print *, "We continue!" 1180 print *, "Number of iteration of the PEM : year_iter=", year_iter 1181 endif 1182 1183 global_ave_press_old=global_ave_press_new 1184 1185 enddo ! big time iteration loop of the pem 1186 1187 1188 !---------------------------- END RUN PEM --------------------- 1189 1190 !---------------------------- OUTPUT --------------------- 1191 1192 !------------------------ 1193 968 else 969 write(*,*) "We continue!" 970 write(*,*) "Number of iteration of the PEM : year_iter=", year_iter 971 endif 972 973 global_ave_press_old=global_ave_press_new 974 975 enddo ! big time iteration loop of the pem 976 !------------------------------ END RUN -------------------------------- 977 978 !------------------------------- OUTPUT -------------------------------- 979 !------------------------ 1194 980 ! III Output 1195 981 ! III_a Update surface value for the PCM start files 1196 1197 !------------------------ 1198 982 !------------------------ 1199 983 ! III_a.1 Ice update (for startfi) 1200 984 1201 985 ! H2O ice 1202 DO ig=1,ngrid 1203 if(watercaptag(ig)) then 1204 watercap_sum=0. 1205 DO islope=1,nslope 1206 if(qsurf(ig,igcm_h2o_ice,islope).GT. (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 1207 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 1208 else 1209 ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap 1210 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.)) 1211 qsurf(ig,igcm_h2o_ice,islope)=0. 986 do ig = 1,ngrid 987 if (watercaptag(ig)) then 988 watercap_sum = 0. 989 do islope = 1,nslope 990 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 991 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 992 else 993 ! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap 994 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.)) 995 qsurf(ig,igcm_h2o_ice,islope)=0. 996 endif 997 watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 998 watercap(ig,islope) = 0. 999 enddo 1000 water_reservoir(ig) = water_reservoir(ig) + watercap_sum 1001 endif 1002 enddo 1003 1004 do ig = 1,ngrid 1005 water_sum = 0. 1006 do islope = 1,nslope 1007 water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1008 enddo 1009 if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming. 1010 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 1011 watercaptag(ig) = .true. 1012 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 1013 do islope = 1,nslope 1014 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.) 1015 enddo 1016 endif 1017 else ! let's check that the infinite source of water disapear 1018 if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then 1019 watercaptag(ig) = .false. 1020 do islope = 1,nslope 1021 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) 1022 enddo 1023 water_reservoir(ig) = 0. 1024 endif 1025 endif 1026 enddo 1027 1028 ! CO2 ice 1029 do ig = 1,ngrid 1030 do islope = 1,nslope 1031 if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then 1032 perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1033 qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1034 endif 1035 enddo 1036 enddo 1037 1038 ! III_a.2 Tsoil update (for startfi) 1039 if (soil_pem) then 1040 call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) 1041 tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) 1042 endif 1043 1044 ! III_a.4 Pressure (for start) 1045 do i = 1,ip1jmp1 1046 ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM 1047 enddo 1048 1049 do i = 1,ngrid 1050 ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM 1051 enddo 1052 1053 ! III_a.5 Tracer (for start) 1054 allocate(zplev_new(ngrid,nlayer + 1)) 1055 1056 do l = 1,nlayer + 1 1057 do ig = 1,ngrid 1058 zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) 1059 enddo 1060 enddo 1061 1062 do nnq = 1,nqtot 1063 if (noms(nnq) /= "co2") then 1064 do l = 1,llm - 1 1065 do ig = 1,ngrid 1066 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)) 1067 enddo 1068 q(:,llm,nnq) = q(:,llm - 1,nnq) 1069 enddo 1070 else 1071 do l = 1,llm - 1 1072 do ig = 1,ngrid 1073 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)) & 1074 + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ & 1075 (zplev_new(ig,l) - zplev_new(ig,l + 1)) 1076 enddo 1077 q(:,llm,nnq) = q(:,llm - 1,nnq) 1078 enddo 1079 endif 1080 enddo 1081 1082 ! Conserving the tracers mass for GCM start files 1083 do nnq = 1,nqtot 1084 do ig = 1,ngrid 1085 do l = 1,llm - 1 1086 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 1087 extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) 1088 q(ig,l,nnq)=1. 1089 q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2)) 1090 write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" 1212 1091 endif 1213 watercap_sum=watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1214 watercap(ig,islope)=0. 1215 enddo 1216 water_reservoir(ig)=water_reservoir(ig)+watercap_sum 1217 endif 1218 enddo 1219 1220 DO ig=1,ngrid 1221 water_sum = 0. 1222 DO islope=1,nslope 1223 water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1224 ENDDO 1225 if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming. 1226 if(water_sum.gt.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 1227 watercaptag(ig)=.true. 1228 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 1229 DO islope = 1,nslope 1230 qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.) 1231 ENDDO 1232 endif 1233 else ! let's check that the infinite source of water disapear 1234 if((water_sum + water_reservoir(ig)).lt.threshold_water_frost2perenial) then 1235 watercaptag(ig)=.false. 1236 DO islope = 1,nslope 1237 qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.) 1238 ENDDO 1239 water_reservoir(ig) = 0. 1240 endif 1241 endif 1242 enddo 1243 1244 ! CO2 ice 1245 1246 DO ig=1,ngrid 1247 DO islope=1,nslope 1248 if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then 1249 perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1250 qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope) 1251 endif 1252 ENDDO 1253 ENDDO 1254 1255 1256 ! III_a.2 Tsoil update (for startfi) 1257 1258 if(soil_pem) then 1259 call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil) 1260 tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen) 1261 endif !soil_pem 1262 1263 ! III_a.4 Pressure (for start) 1264 do i=1,ip1jmp1 1265 ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM 1266 enddo 1267 1268 do i = 1,ngrid 1269 ps_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM 1270 enddo 1271 1272 ! III_a.5 Tracer (for start) 1273 allocate(zplev_new(ngrid,nlayer+1)) 1274 1275 do l=1,nlayer+1 1276 do ig=1,ngrid 1277 zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig) 1278 enddo 1279 enddo 1280 1281 DO nnq=1,nqtot 1282 if (noms(nnq).NE."co2") then 1283 DO l=1,llm-1 1284 DO ig=1,ngrid 1285 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)) 1286 ENDDO 1287 q(:,llm,nnq)=q(:,llm-1,nnq) 1288 ENDDO 1289 else 1290 DO l=1,llm-1 1291 DO ig=1,ngrid 1292 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)) & 1293 + ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - & 1294 (zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / & 1295 (zplev_new(ig,l)-zplev_new(ig,l+1)) 1296 ENDDO 1297 q(:,llm,nnq)=q(:,llm-1,nnq) 1298 ENDDO 1299 endif 1300 ENDDO 1301 1302 ! Conserving the tracers's mass for GCM start files 1303 DO nnq=1,nqtot 1304 DO ig=1,ngrid 1305 DO l=1,llm-1 1306 if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") .and. (noms(nnq).NE."stormdust_number") .and. (noms(nnq).NE."topdust_number")) then 1307 extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1)) 1308 q(ig,l,nnq)=1. 1309 q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2)) 1310 write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number" 1311 endif 1312 if(q(ig,l,nnq).LT.0) then 1313 q(ig,l,nnq)=1E-30 1314 endif 1315 ENDDO 1316 enddo 1317 enddo 1318 1319 !------------------------ 1320 if(evol_orbit_pem) then 1321 call recomp_orb_param(year_iter) 1322 endif 1323 1092 if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 1093 enddo 1094 enddo 1095 enddo 1096 1097 if (evol_orbit_pem) call recomp_orb_param(year_iter) 1098 1099 !------------------------ 1324 1100 ! III Output 1325 ! III_a Update surface value for the PCM start files 1326 ! III_b Write start and starfi.nc 1327 1328 !------------------------ 1329 1330 ! III_b.1 WRITE restart.nc 1331 1332 ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys 1333 pday=day_ini 1334 ztime_fin=0. 1335 1336 allocate(p(ip1jmp1,nlayer+1)) 1101 ! III_b Write restart_evol.nc and restartfi_evol.nc 1102 !------------------------ 1103 ! III_b.1 Write restart_evol.nc 1104 ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys 1105 pday = day_ini 1106 ztime_fin = 0. 1107 1108 allocate(p(ip1jmp1,nlayer + 1)) 1337 1109 #ifndef CPP_1D 1338 CALL pression (ip1jmp1,ap,bp,ps,p) 1339 CALL massdair(p,masse) 1340 1341 CALL dynredem0("restart_evol.nc", day_ini, phis) 1342 1343 CALL dynredem1("restart_evol.nc", & 1344 time_0,vcov,ucov,teta,q,masse,ps) 1345 print *, "restart_evol.nc has been written" 1346 1110 call pression (ip1jmp1,ap,bp,ps,p) 1111 call massdair(p,masse) 1112 call dynredem0("restart_evol.nc", day_ini, phis) 1113 call dynredem1("restart_evol.nc", & 1114 time_0,vcov,ucov,teta,q,masse,ps) 1115 write(*,*) "restart_evol.nc has been written" 1347 1116 #else 1348 DOnnq = 1, nqtot1349 1350 ENDDO1117 do nnq = 1, nqtot 1118 call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf) 1119 enddo 1351 1120 #endif 1352 1121 1353 ! III_b.2 W RITE restartfi.nc1122 ! III_b.2 Write restartfi_evol.nc 1354 1123 #ifndef CPP_STD 1355 call physdem0("restartfi_evol.nc",longitude,latitude, & 1356 nsoilmx,ngrid,nlayer,nq, & 1357 ptimestep,pday,0.,cell_area, & 1358 albedodat,inertiedat,def_slope, & 1359 subslope_dist) 1360 1361 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & 1362 ptimestep,ztime_fin, & 1363 tsurf,tsoil,inertiesoil,albedo, & 1364 emis,q2,qsurf,tauscaling,totcloudfrac,wstar,& 1365 watercap,perenial_co2ice) 1124 call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & 1125 nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & 1126 inertiedat,def_slope,subslope_dist) 1127 1128 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & 1129 ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & 1130 albedo,emis,q2,qsurf,tauscaling,totcloudfrac, & 1131 wstar,watercap,perenial_co2ice) 1366 1132 #else 1367 call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, & 1368 ptimestep,pday,time_phys,cell_area, & 1369 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 1370 1371 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & 1372 ptimestep,ztime_fin, & 1373 tsurf,tsoil,emis,q2,qsurf, & 1374 cloudfrac,totcloudfrac,hice, & 1375 rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) 1133 call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, & 1134 nlayer,nq,ptimestep,pday,time_phys,cell_area, & 1135 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 1136 1137 call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, & 1138 ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, & 1139 cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, & 1140 tslab,tsea_ice,sea_ice) 1376 1141 #endif 1377 1378 print *, "restartfi_evol.nc has been written" 1379 !------------------------ 1380 1142 write(*,*) "restartfi_evol.nc has been written" 1143 1144 !------------------------ 1381 1145 ! III Output 1382 ! III_a Update surface value for the PCM start files 1383 ! III_b Write start and starfi.nc 1384 ! III_c Write start_pem 1385 1386 !------------------------ 1387 call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & 1388 float(day_ini),0.,nslope,def_slope,subslope_dist) 1389 1390 1391 call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , & 1392 tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, & 1393 co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) 1394 1395 call info_run_PEM(year_iter, criterion_stop) 1396 1397 print *, "restartfi_PEM.nc has been written" 1398 print *, "The PEM had run for ", year_iter, " years." 1399 print *, "LL & RV : So far so good" 1400 1401 deallocate(vmr_co2_gcm) 1402 deallocate(ps_timeseries) 1403 deallocate(tsurf_GCM_timeseries) 1404 deallocate(q_co2_PEM_phys) 1405 deallocate(q_h2o_PEM_phys) 1406 deallocate(co2_ice_GCM) 1407 deallocate(watersurf_density_ave) 1408 deallocate(watersoil_density_timeseries) 1409 deallocate(Tsurfave_before_saved) 1410 deallocate(tsoil_phys_PEM_timeseries) 1411 deallocate(watersoil_density_PEM_timeseries) 1412 deallocate(watersoil_density_PEM_ave) 1413 deallocate(delta_co2_adsorbded) 1414 deallocate(delta_h2o_adsorbded) 1415 deallocate(vmr_co2_pem_phys) 1146 ! III_c Write restartfi_PEM.nc 1147 !------------------------ 1148 call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, & 1149 float(day_ini),0.,nslope,def_slope,subslope_dist) 1150 1151 1152 call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , & 1153 tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, & 1154 co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) 1155 1156 call info_run_PEM(year_iter,criterion_stop) 1157 1158 write(*,*) "restartfi_PEM.nc has been written" 1159 write(*,*) "The PEM had run for ", year_iter, " years." 1160 write(*,*) "LL & RV : So far so good" 1161 1162 deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys) 1163 deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved) 1164 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) 1165 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys) 1166 !----------------------------- END OUTPUT ------------------------------ 1416 1167 1417 1168 END PROGRAM pem 1169 -
trunk/LMDZ.MARS/changelog.txt
r3026 r3028 4169 4169 == 09/08/2023 == JBC 4170 4170 Small fixes to be able to run the Mars PCM 1D without "water" + Improvements/addition of scripts in deftank/pem to run the PEM 1D model according to Laskar orbital parameters. 4171 4172 == 14/08/2023 == JBC 4173 Related to commit r3026: improvement of error message in initracer.F (now it gives correctly the only identified tracers) + one small correction to run PCM 1D without water. 4174 Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
Note: See TracChangeset
for help on using the changeset viewer.
