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