Changeset 3327 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- May 15, 2024, 4:45:14 PM (6 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 6 edited
-
changelog.txt (modified) (1 diff)
-
deftank/run_PEM.def (modified) (3 diffs)
-
ice_table_mod.F90 (modified) (11 diffs)
-
pem.F90 (modified) (9 diffs)
-
soil_thermalproperties_mod.F90 (modified) (5 diffs)
-
stopping_crit_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3324 r3327 300 300 Correction of a bug when setting the albedo / emissivity: it should be the PCM 301 301 which sets the emissivity / albedo of the CO2 ice (because of CO2 snowfall) and 302 not the PEM. I thus removed the lines in the PEM accordingly 302 not the PEM. I thus removed the lines in the PEM accordingly 303 304 == 15/05/2024 == JBC 305 - Update of some default parameters; 306 - Correction of a bug when the ice table was too low for the PEM subsurface discretization; 307 - Correction of the computation for the change of sublimating ice necessary to the stopping criteria (the mistake was introduced in r3149) + simplification of the algorithm + renaming of variables with more explicit names; 308 - Cleaning of "soil_thermalproperties_mod.F90". -
trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def
r3320 r3327 12 12 #---------- Orbital parameters ----------# 13 13 # Do you want to follow an orbital forcing read in "obl_ecc_lsp.asc"? Default = .false. 14 evol_orbit_pem=.true.14 # evol_orbit_pem=.false. 15 15 16 16 # If evol_orbit_pem=.true., number of Earth years before present to start the PEM run? Default = 0 17 year_earth_bp_ini=-100000017 # year_earth_bp_ini=0 18 18 19 19 # Do you want to vary the obliquity when following "obl_ecc_lsp.asc"? Default = .true. 20 var_obl=.true.20 # var_obl=.false. 21 21 22 22 # Do you want to vary the eccentricity when following "obl_ecc_lsp.asc"? Default = .true. 23 var_ecc=.true.23 # var_ecc=.false. 24 24 25 25 # Do you want to vary the ls perihelie when following "obl_ecc_lsp.asc"? Default = .true. 26 var_lsp=.true.26 # var_lsp=.false. 27 27 28 28 #---------- Stopping criteria parameters ----------# 29 29 # If evol_orbit_pem=.false., maximal number of iterations if no stopping criterion is reached? Default=100000000 30 # Max_iter_pem=10000000030 # Max_iter_pem=10 31 31 32 32 # Acceptance rate of sublimating H2O ice surface change? Default = 0.2 … … 46 46 # soil_pem=.true. 47 47 48 # Do you want to run with adsoprtion in the PEM? Default = . true.48 # Do you want to run with adsoprtion in the PEM? Default = .false. 49 49 # adsorption_pem=.true. 50 50 … … 90 90 91 91 # Do you want the CO2 ice to flow along subslope inside a cell? Default = .true. 92 # co2ice_flow=. true.92 # co2ice_flow=.false. 93 93 94 94 #---------- Layering parameters ----------# -
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r3321 r3327 16 16 17 17 !----------------------------------------------------------------------- 18 18 contains 19 19 !----------------------------------------------------------------------- 20 20 SUBROUTINE ini_ice_table_porefilling(ngrid,nslope) … … 58 58 implicit none 59 59 60 ! 61 ! ------ --60 ! Inputs 61 ! ------ 62 62 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 63 63 logical, dimension(ngrid), intent(in) :: watercaptag ! Boolean to check the presence of a perennial glacier … … 65 65 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: rhowatersoil_ave ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3] 66 66 real, dimension(ngrid,nslope), intent(in) :: regolith_inertia ! Thermal inertia of the regolith layer [SI] 67 ! 68 ! ------ --67 ! Ouputs 68 ! ------ 69 69 real, dimension(ngrid,nslope), intent(out) :: ice_table_beg ! ice table depth [m] 70 70 real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m] 71 ! 72 ! ------ --71 ! Locals 72 ! ------ 73 73 integer :: ig, islope, isoil, isoilend ! loop variables 74 74 real, dimension(nsoil_PEM) :: diff_rho ! difference of water vapor density between the surface and at depth [kg/m^3] … … 77 77 real :: stretch ! stretch factor to improve the convergence of the ice table 78 78 real :: wice_inertia ! Water Ice thermal Inertia [USI] 79 ! 80 ! ---- --79 ! Code 80 ! ---- 81 81 previous_icetable_depth = ice_table_beg 82 82 do ig = 1,ngrid … … 147 147 use comslope_mod, only: subslope_dist, def_slope_mean 148 148 use constants_marspem_mod, only: porosity 149 use vdifc_mod, only: compute_Tice150 149 #ifndef CPP_STD 151 150 use comcstfi_h, only: pi … … 156 155 implicit none 157 156 158 ! 159 ! ------ --157 ! Inputs 158 ! ------ 160 159 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 161 160 real, dimension(ngrid,nslope), intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m] … … 164 163 real, dimension(ngrid,nslope), intent(in) :: tsurf ! Surface temperature [K] 165 164 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil ! Soil temperature [K] 166 ! 167 ! ------- --165 ! Outputs 166 ! ------- 168 167 real, dimension(ngrid), intent(out) :: delta_m_h2o ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2] 169 ! 170 !------- --168 ! Locals 169 !------- 171 170 integer :: ig, islope, ilay, iref ! loop index 172 171 real, dimension(ngrid,nslope) :: rho ! density of water ice [kg/m^3] 173 172 real, dimension(ngrid,nslope) :: Tice ! ice temperature [k] 174 ! 175 ! ---- --173 ! Code 174 ! ---- 176 175 rho = 0. 177 176 Tice = 0. … … 179 178 do ig = 1,ngrid 180 179 do islope = 1,nslope 181 call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope)) 180 print*, ig, islope, ice_table_depth(ig,islope) 181 call compute_Tice_pem(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope)) 182 182 rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 183 183 enddo … … 224 224 real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] 225 225 real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase 226 ! 227 !------- --226 ! Locals 227 !------- 228 228 real, dimension(nsoilmx_PEM) :: eta ! constriction 229 229 real, dimension(nsoilmx_PEM) :: dz_psat ! first derivative of the vapor pressure at saturation … … 237 237 integer :: index_firstice ! first index where ice appears (i.e., f > 0) 238 238 real :: massfillabove, massfillafter ! h2O mass above and after index_geothermal 239 ! 240 !--------- --239 ! Constant 240 !--------- 241 241 real, parameter :: pvap2rho = 18.e-3/8.314 242 ! 243 !----- --242 ! Code 243 !----- 244 244 ! 0. Compute constriction over the layer 245 245 ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 … … 364 364 END SUBROUTINE constriction 365 365 366 !----------------------------------------------------------------------- 367 SUBROUTINE compute_Tice_pem(nsoil, ptsoil, ptsurf, ice_depth, Tice) 368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 369 !!! 370 !!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells. 371 !!! 372 !!! Author: LL 373 !!! 374 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 375 376 use comsoil_h_PEM, only: layer_PEM, mlayer_PEM 377 378 implicit none 379 380 ! Inputs 381 ! ------ 382 integer, intent(in) :: nsoil ! Number of soil layers 383 real, dimension(nsoil), intent(in) :: ptsoil ! Soil temperature (K) 384 real, intent(in) :: ptsurf ! Soil temperature (K) 385 real, intent(in) :: ice_depth ! Ice depth (m) 386 387 ! Outputs 388 ! ------ 389 real, intent(out) :: Tice ! Ice temperatures (K) 390 391 ! Locals 392 ! ------ 393 integer :: ik ! Loop variables 394 integer :: indexice ! Index of the ice 395 396 ! Code 397 !----- 398 indexice = -1 399 if (ice_depth >= mlayer_PEM(nsoil - 1)) then 400 Tice = ptsoil(nsoil) 401 else 402 if(ice_depth < mlayer_PEM(0)) then 403 indexice = 0. 404 else 405 do ik = 0,nsoil - 2 ! go through all the layers to find the ice locations 406 if (mlayer_PEM(ik) <= ice_depth .and. mlayer_PEM(ik + 1) > ice_depth) then 407 indexice = ik + 1 408 exit 409 endif 410 enddo 411 endif 412 if (indexice < 0) then 413 call abort_physic("compute_Tice_pem","subsurface ice is below the last soil layer",1) 414 else 415 if(indexice >= 1) then ! Linear inteprolation between soil temperature 416 Tice = (ptsoil(indexice) - ptsoil(indexice + 1))/(mlayer_PEM(indexice - 1) - mlayer_PEM(indexice))*(ice_depth - mlayer_PEM(indexice)) + ptsoil(indexice + 1) 417 else ! Linear inteprolation between the 1st soil temperature and the surface temperature 418 Tice = (ptsoil(1) - ptsurf)/mlayer_PEM(0)*(ice_depth - mlayer_PEM(0)) + ptsoil(1) 419 endif ! index ice >= 1 420 endif !indexice < 0 421 endif ! icedepth > mlayer_PEM(nsoil - 1) 422 423 END SUBROUTINE compute_Tice_pem 424 366 425 END MODULE ice_table_mod -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3324 r3327 158 158 159 159 ! Variables for h2o_ice evolution 160 real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM 161 real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2] 162 real :: h2o_surf_ini ! Initial surface of sublimating h2o ice 163 real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa] 164 real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step 165 real :: global_avg_press_new ! constant: Global average pressure of current time step 166 real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] 167 real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2] 168 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] 169 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step 170 integer :: stopPEM ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu 171 real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 172 real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] 173 integer :: timelen ! # time samples 174 real :: ave ! intermediate varibale to compute average 175 real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 176 real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 160 real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM 161 real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2] 162 real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice 163 logical, dimension(:,:), allocatable :: ini_h2oice_sublim ! Logical array to know if h2o ice is sublimating 164 real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa] 165 real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step 166 real :: global_avg_press_new ! constant: Global average pressure of current time step 167 real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] 168 real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2] 169 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] 170 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step 171 integer :: stopPEM ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu 172 real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] 173 real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] 174 integer :: timelen ! # time samples 175 real :: ave ! intermediate varibale to compute average 176 real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) 177 real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 177 178 178 179 ! Variables for co2_ice evolution 179 real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM 180 real, dimension(:,:), allocatable :: co2_ice_ini ! Initial amount of co2 ice in the PEM 181 real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2] 182 real :: co2_surf_ini ! Initial surface of sublimating co2 ice 183 real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] 184 real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM 185 real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] 180 real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM 181 real, dimension(:,:), allocatable :: co2_ice_ini ! Initial amount of co2 ice in the PEM 182 real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2] 183 real :: co2ice_ini_surf ! Initial surface of sublimating co2 ice 184 logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating 185 real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] 186 real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM 187 real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] 186 188 187 189 ! Variables for stratification (layering) evolution … … 556 558 ! We save the places where h2o ice is sublimating 557 559 ! We compute the surface of h2o ice sublimating 558 co2_surf_ini = 0. 559 h2o_surf_ini = 0. 560 allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope)) 561 co2ice_ini_surf = 0. 562 h2oice_ini_surf = 0. 560 563 Total_surface = 0. 564 ini_co2ice_sublim = .false. 565 ini_h2oice_sublim = .false. 561 566 do i = 1,ngrid 562 567 Total_surface = Total_surface + cell_area(i) 563 568 do islope = 1,nslope 564 if (tend_co2_ice(i,islope) < 0.) co2_surf_ini = co2_surf_ini + cell_area(i)*subslope_dist(i,islope) 565 if (tend_h2o_ice(i,islope) < 0.) h2o_surf_ini = h2o_surf_ini + cell_area(i)*subslope_dist(i,islope) 569 if (tend_co2_ice(i,islope) < 0.) then 570 ini_co2ice_sublim(i,islope) = .true. 571 co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope) 572 endif 573 if (tend_h2o_ice(i,islope) < 0.) then 574 ini_h2oice_sublim(i,islope) = .true. 575 h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) 576 endif 566 577 enddo 567 578 enddo 568 579 569 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2 _surf_ini570 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2o _surf_ini580 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf 581 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf 571 582 write(*,*) "Total surface of the planet =", Total_surface 572 583 allocate(zplev_PCM(ngrid,nlayer + 1)) … … 778 789 do ig = 1,ngrid 779 790 do islope = 1,nslope 780 if (co2_ice_ini(ig,islope) > 0. 5.and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice791 if (co2_ice_ini(ig,islope) > 0. .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice 781 792 if (latitude_deg(ig) > 0) then 782 793 do ig_loop = ig,ngrid 783 794 do islope_loop = islope,iflat,-1 784 if (co2_ice_ini(ig_loop,islope_loop) < 0.5.and. co2_ice(ig_loop,islope_loop) < 1.e-10) then795 if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then 785 796 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) 786 797 bool_sublim = .true. … … 793 804 do ig_loop = ig,1,-1 794 805 do islope_loop = islope,iflat 795 if (co2_ice_ini(ig_loop,islope_loop) < 0.5.and. co2_ice(ig_loop,islope_loop) < 1.e-10) then806 if (co2_ice_ini(ig_loop,islope_loop) < 1.e-10 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then 796 807 tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop) 797 808 bool_sublim = .true. … … 802 813 enddo 803 814 endif 804 co2_ice_ini(ig,islope) = 0 815 co2_ice_ini(ig,islope) = 0. 805 816 if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then 806 817 albedo(ig,1,islope) = albedo_h2o_frost … … 954 965 ! II_g Checking the stopping criterion 955 966 !------------------------ 956 call stopping_crit_h2o_ice(cell_area,h2o _surf_ini,h2o_ice,stopPEM,ngrid)957 call stopping_crit_co2(cell_area,co2 _surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)967 call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid) 968 call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope) 958 969 959 970 year_iter = year_iter + dt_pem … … 1107 1118 #ifndef CPP_STD 1108 1119 call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & 1109 nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, 1120 nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & 1110 1121 inertiedat,def_slope,subslope_dist) 1111 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &1122 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & 1112 1123 ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & 1113 1124 albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, & … … 1115 1126 #else 1116 1127 call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & 1117 nlayer,nq,ptimestep,pday,time_phys,cell_area, 1128 nlayer,nq,ptimestep,pday,time_phys,cell_area, & 1118 1129 albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) 1119 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &1130 call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & 1120 1131 ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, & 1121 1132 cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, & … … 1154 1165 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave) 1155 1166 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter) 1156 deallocate(co2_ice_ini, stratif)1167 deallocate(co2_ice_ini,ini_co2ice_sublim,ini_h2oice_sublim,stratif) 1157 1168 !----------------------------- END OUTPUT ------------------------------ 1158 1169 -
trunk/LMDZ.COMMON/libf/evolution/soil_thermalproperties_mod.F90
r3206 r3327 5 5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 6 6 !!! 7 !!! Purpose: Compute the soil thermal properties 7 !!! Purpose: Compute the soil thermal properties 8 8 !!! Author: LL, 04/2023 9 9 !!! … … 19 19 ! -------- 20 20 ! 21 ! author: LL, 04/2023 21 ! author: LL, 04/2023 22 22 ! ------- 23 ! 23 ! 24 24 !======================================================================= 25 25 … … 37 37 logical, intent(in) :: ispureice ! Boolean to check if ice is massive or just pore filling 38 38 real, intent(in) :: pore_filling ! ice pore filling in each layer (1) 39 real, intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2) 39 real, intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2) 40 40 real, intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2) 41 41 42 42 ! Local Variables 43 43 ! -------------- 44 REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004) 44 REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004) 45 45 !======================================================================= 46 46 ! Beginning of the code … … 53 53 endif 54 54 55 END SUBROUTINE 55 END SUBROUTINE 56 56 !======================================================================= 57 57 … … 64 64 implicit none 65 65 66 ! Input: 67 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil 68 real, intent(in) :: p_avg_new ! Global average surface pressure [Pa] 69 real, intent(in) :: tendencies_waterice(ngrid,nslope) ! Tendencies on the water ice [kg/m^2/year] 70 real, intent(in) :: waterice(ngrid,nslope) ! Surface Water ice [kg/m^2] 71 real, intent(in) :: ice_depth(ngrid,nslope) ! Ice table depth [m] 72 real, intent(in) :: ice_thickness(ngrid,nslope) ! Ice table thickness [m] 66 ! Input: 67 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil 68 real, intent(in) :: p_avg_new ! Global average surface pressure [Pa] 69 real, dimension(ngrid,nslope), intent(in) :: tendencies_waterice ! Tendencies on the water ice [kg/m^2/year] 70 real, dimension(ngrid,nslope), intent(in) :: waterice ! Surface Water ice [kg/m^2] 71 real, dimension(ngrid,nslope), intent(in) :: ice_depth ! Ice table depth [m] 72 real, dimension(ngrid,nslope), intent(in) :: ice_thickness ! Ice table thickness [m] 73 73 74 ! Outputs: 74 75 real, intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal Inertia [J/m^2/K/s^1/2] 76 75 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2] 76 77 77 ! Constants: 78 79 REAL :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2] 80 REAL :: reg_inertie_minvalue = 50. ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] 81 REAL :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] 82 REAL :: ice_inertia ! Inertia of water ice [SI] 83 REAL :: P610 = 610.0 ! current average pressure on Mars [Pa] 84 REAL :: C = 0.0015 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless] 85 REAL :: K = 8.1*1e4 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa] 86 REAL :: Pa2Tor = 1./133.3 ! Conversion from Pa to tor [Pa/tor] 87 78 real :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2] 79 real :: reg_inertie_minvalue = 50. ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] 80 real :: reg_inertie_maxvalue = 370. ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2] 81 real :: ice_inertia ! Inertia of water ice [SI] 82 real :: P610 = 610.0 ! current average pressure on Mars [Pa] 83 real :: C = 0.0015 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [uniteless] 84 real :: K = 8.1*1e4 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [tor, or 133.3Pa] 85 real :: Pa2Tor = 1./133.3 ! Conversion from Pa to tor [Pa/tor] 88 86 89 87 ! Local variables: 90 91 INTEGER :: ig,islope,isoil,iref,iend ! Loop variables 92 REAL :: regolith_inertia(ngrid,nslope) ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2] 93 REAL :: delta ! Difference of depth to compute the thermal inertia in a mixed layer [m] 94 REAL :: ice_bottom_depth ! Bottom depth of the subsurface ice [m] 95 REAL :: d_part ! Regolith particle size [micrometer] 96 97 98 write(*,*) "Update soil propreties" 99 88 integer :: ig, islope, isoil, iref, iend ! Loop variables 89 real, dimension(ngrid,nslope) :: regolith_inertia ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2] 90 real :: delta ! Difference of depth to compute the thermal inertia in a mixed layer [m] 91 real :: ice_bottom_depth ! Bottom depth of the subsurface ice [m] 92 real :: d_part ! Regolith particle size [micrometer] 93 94 write(*,*) "Update soil properties" 95 100 96 ! 1. Modification of the regolith thermal inertia. 101 do islope = 1,nslope 102 regolith_inertia(:,islope) = inertiedat_PEM(:,1) 103 do ig = 1,ngrid 104 if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then 105 regolith_inertia(ig,islope) = TI_regolith_avg 106 endif 107 if(reg_thprop_dependp) then 108 if(TI_PEM(ig,1,islope).lt.reg_inertie_thresold) then 109 d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer 110 regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K))) 111 if(regolith_inertia(ig,islope).gt.reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue 112 if(regolith_inertia(ig,islope).lt.reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue 113 endif 114 endif 115 enddo 116 enddo 117 118 ! 2. Build new Thermal Inertia 119 do islope=1,nslope 120 do ig = 1,ngrid 121 do isoil = 1,index_breccia 122 TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope) 123 enddo 124 if(regolith_inertia(ig,islope).lt.TI_breccia) then 97 do islope = 1,nslope 98 regolith_inertia(:,islope) = inertiedat_PEM(:,1) 99 do ig = 1,ngrid 100 if (tendencies_waterice(ig,islope) < -1.e-5 .and. waterice(ig,islope) == 0) regolith_inertia(ig,islope) = TI_regolith_avg 101 if (reg_thprop_dependp) then 102 if (TI_PEM(ig,1,islope) < reg_inertie_thresold) then 103 d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Tor)**(0.6)))**(-1./(0.11*log10(P610*Pa2Tor/K))) ! compute particle size, in micrometer 104 regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Tor)**(0.6)*d_part**(-0.11*log10(p_avg_new*Pa2Tor/K))) 105 if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue 106 if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue 107 endif 108 endif 109 enddo 110 enddo 111 112 ! 2. Build new Thermal Inertia 113 do islope = 1,nslope 114 do ig = 1,ngrid 115 do isoil = 1,index_breccia 116 TI_PEM(ig,isoil,islope) = regolith_inertia(ig,islope) 117 enddo 118 if (regolith_inertia(ig,islope) < TI_breccia) then 125 119 !!! transition 126 127 TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/&128 (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &129 ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))130 do isoil =index_breccia+2,index_bedrock131 TI_PEM(ig,isoil,islope) = TI_breccia132 enddo 133 else ! we keep the high ti values134 do isoil =index_breccia+1,index_bedrock135 TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope)136 enddo 137 endif ! TI PEM and breccia comparison138 !! transition139 delta = depth_bedrock140 TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/&141 (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &142 ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))143 do isoil=index_bedrock+2,nsoil_PEM120 delta = depth_breccia 121 TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ & 122 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + & 123 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2)))) 124 do isoil = index_breccia + 2,index_bedrock 125 TI_PEM(ig,isoil,islope) = TI_breccia 126 enddo 127 else ! we keep the high ti values 128 do isoil = index_breccia + 1,index_bedrock 129 TI_PEM(ig,isoil,islope) = TI_PEM(ig,index_breccia,islope) 130 enddo 131 endif ! TI PEM and breccia comparison 132 !!! transition 133 delta = depth_bedrock 134 TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ & 135 (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + & 136 ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2)))) 137 do isoil = index_bedrock + 2,nsoil_PEM 144 138 TI_PEM(ig,isoil,islope) = TI_bedrock 145 enddo146 147 ENDDO! islope139 enddo 140 enddo ! ig 141 enddo ! islope 148 142 149 143 ! 3. Build new TI in case of ice table 150 do ig=1,ngrid 151 do islope=1,nslope 152 if (ice_depth(ig,islope).gt.-1e-6) then 153 ! 3.0 Case where it is perennial ice 154 if (ice_depth(ig,islope).lt. 1e-10) then 155 call ice_thermal_properties(.true.,1.,0.,ice_inertia) 156 do isoil = 1,nsoil_PEM 157 TI_PEM(ig,isoil,islope)=ice_inertia 158 enddo 159 else 160 161 call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia) 144 do ig = 1,ngrid 145 do islope = 1,nslope 146 if (ice_depth(ig,islope) > -1.e-6) then 147 ! 3.0 Case where it is perennial ice 148 if (ice_depth(ig,islope) < 1.e-10) then 149 call ice_thermal_properties(.true.,1.,0.,ice_inertia) 150 do isoil = 1,nsoil_PEM 151 TI_PEM(ig,isoil,islope) = ice_inertia 152 enddo 153 else 154 call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),ice_inertia) 162 155 ! 3.1.1 find the index of the mixed layer 163 iref=0 ! initialize iref 164 do isoil=1,nsoil_PEM ! loop on layers to find the beginning of the ice table 165 if (ice_depth(ig,islope).ge.layer_PEM(isoil)) then 166 iref=isoil ! pure regolith layer up to here 167 else 168 ! correct iref was obtained in previous cycle 169 exit 170 endif 171 enddo 156 iref = 0 ! initialize iref 157 do isoil = 1,nsoil_PEM ! loop on layers to find the beginning of the ice table 158 if (ice_depth(ig,islope) >= layer_PEM(isoil)) then 159 iref = isoil ! pure regolith layer up to here 160 else 161 exit ! correct iref was obtained in previous cycle 162 endif 163 enddo 172 164 ! 3.1.2 find the index of the end of the ice table 173 iend=0 ! initialize iend 174 ice_bottom_depth = ice_depth(ig,islope)+ice_thickness(ig,islope) 175 do isoil=1,nsoil_PEM ! loop on layers to find the end of the ice table 176 if (ice_bottom_depth.ge.layer_PEM(isoil)) then 177 iend=isoil ! pure regolith layer up to here 178 else 179 ! correct iref was obtained in previous cycle 180 exit 181 endif 182 enddo 183 ! 3.2 Build the new ti 184 do isoil=1,iref 185 TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope) 186 enddo 187 if (iref.lt.nsoil_PEM) then 188 if(iref.eq.iend) then 189 ! Ice table begins and end in the same mixture Mixtures with three components 190 if (iref.ne.0) then ! 191 ! mixed layer 192 TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ & 193 (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ & 194 ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ & 195 ((layer_PEM(iref+1)-ice_bottom_depth)/(TI_PEM(ig,iref+1,islope)**2)))) 196 else ! first layer is already a mixed layer 197 ! (ie: take layer(iref=0)=0) 198 TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ & 199 (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ & 200 ((ice_bottom_depth-ice_depth(ig,islope))/(ice_inertia**2))+ & 201 ((layer_PEM(2)-ice_bottom_depth)/(TI_PEM(ig,2,islope)**2)))) 202 endif ! of if (iref.ne.0) 203 else 204 if (iref.ne.0) then 205 ! mixed layer 206 TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ & 207 (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ & 208 ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2)))) 209 else ! first layer is already a mixed layer 210 ! (ie: take layer(iref=0)=0) 211 TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ & 212 (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ & 213 ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2)))) 214 endif ! of if (iref.ne.0) 215 endif ! iref.eq.iend 216 ! 3.3 Build the new ti 217 do isoil=iref+2,iend 218 TI_PEM(ig,isoil,islope)=ice_inertia 219 enddo 220 221 222 if(iend.lt.nsoil_PEM) then 223 TI_PEM(ig,iend+1,islope)=sqrt((layer_PEM(iend+1)-layer_PEM(iend))/ & 224 (((ice_bottom_depth-layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + & 225 ((layer_PEM(iend+1)-ice_bottom_depth)/(ice_inertia**2)))) 226 endif 227 228 229 endif ! of if (iref.lt.(nsoilmx)) 230 endif ! permanent glaciers 231 endif !ice_depth(ig,islope) > 0. 232 ! write(*,*) 'TI=', TI_PEM(ig,:,islope) 165 iend = 0 ! initialize iend 166 ice_bottom_depth = ice_depth(ig,islope) + ice_thickness(ig,islope) 167 do isoil = 1,nsoil_PEM ! loop on layers to find the end of the ice table 168 if (ice_bottom_depth >= layer_PEM(isoil)) then 169 iend = isoil ! pure regolith layer up to here 170 else 171 exit ! correct iref was obtained in previous cycle 172 endif 173 enddo 174 ! 3.2 Build the new ti 175 do isoil = 1,iref 176 TI_PEM(ig,isoil,islope) = TI_PEM(ig,1,islope) 177 enddo 178 if (iref < nsoil_PEM) then 179 if (iref == iend) then 180 ! Ice table begins and end in the same mixture Mixtures with three components 181 if (iref /= 0) then ! mixed layer 182 TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/ & 183 (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + & 184 ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + & 185 ((layer_PEM(iref + 1) - ice_bottom_depth)/(TI_PEM(ig,iref + 1,islope)**2)))) 186 else ! first layer is already a mixed layer 187 ! (ie: take layer(iref=0)=0) 188 TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ & 189 (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + & 190 ((ice_bottom_depth - ice_depth(ig,islope))/(ice_inertia**2)) + & 191 ((layer_PEM(2) - ice_bottom_depth)/(TI_PEM(ig,2,islope)**2)))) 192 endif ! of if (iref /= 0) 193 else 194 if (iref /= 0) then ! mixed layer 195 TI_PEM(ig,iref + 1,islope) = sqrt((layer_PEM(iref + 1) - layer_PEM(iref))/ & 196 (((ice_depth(ig,islope) - layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2)) + & 197 ((layer_PEM(iref + 1) - ice_depth(ig,islope))/(ice_inertia**2)))) 198 else ! first layer is already a mixed layer 199 ! (ie: take layer(iref=0)=0) 200 TI_PEM(ig,1,islope) = sqrt((layer_PEM(1))/ & 201 (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2)) + & 202 ((layer_PEM(1) - ice_depth(ig,islope))/(ice_inertia**2)))) 203 endif ! of if (iref /= 0) 204 endif ! iref == iend 205 ! 3.3 Build the new ti 206 do isoil = iref + 2,iend 207 TI_PEM(ig,isoil,islope) = ice_inertia 208 enddo 209 if (iend < nsoil_PEM) then 210 TI_PEM(ig,iend + 1,islope) = sqrt((layer_PEM(iend + 1) - layer_PEM(iend))/ & 211 (((ice_bottom_depth - layer_PEM(iend))/(TI_PEM(ig,iend,islope)**2)) + & 212 ((layer_PEM(iend + 1) - ice_bottom_depth)/(ice_inertia**2)))) 213 endif 214 endif ! of if (iref < nsoilmx) 215 endif ! permanent glaciers 216 endif ! ice_depth(ig,islope) > 0. 217 ! write(*,*) 'TI=', TI_PEM(ig,:,islope) 233 218 enddo !islope 234 219 enddo !ig 235 220 236 221 END SUBROUTINE update_soil_thermalproperties -
trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90
r3206 r3327 14 14 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 15 15 16 SUBROUTINE stopping_crit_h2o_ice(cell_area, surf_ini,h2o_ice,stopPEM,ngrid)16 SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid) 17 17 18 18 use time_evol_mod, only: h2o_ice_crit … … 26 26 ! 27 27 !======================================================================= 28 29 ! arguments:30 ! ---------- 31 ! INPUT 32 integer, intent(in) :: ngrid ! # of physical grid points 33 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells34 real, dimension(ngrid,nslope), intent(in) :: h2o_ice ! Actual density of h2o ice 35 real, intent(in) :: surf_ini ! Initial surface of h2o ice that was sublimating 36 ! OUTPUT28 ! Inputs 29 !------- 30 integer, intent(in) :: ngrid ! # of physical grid points 31 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells 32 real, dimension(ngrid,nslope), intent(in) :: h2o_ice ! Actual density of h2o ice 33 real, intent(in) :: h2oice_ini_surf ! Initial surface of sublimating h2o ice 34 logical, dimension(ngrid,nslope), intent(in) :: ini_h2oice_sublim ! Grid points where h2o ice was initially sublimating 35 ! Outputs 36 !-------- 37 37 integer, intent(inout) :: stopPEM ! Stopping criterion code 38 ! local:39 ! 40 integer :: i, islope ! Loop41 real :: surf_now! Current surface of h2o ice38 ! Locals 39 ! ------ 40 integer :: i, islope ! Loop 41 real :: h2oice_now_surf ! Current surface of h2o ice 42 42 43 43 !======================================================================= 44 ! Computation of the present surface of h2o ice 45 surf_now= 0.44 ! Computation of the present surface of h2o ice still sublimating 45 h2oice_now_surf = 0. 46 46 do i = 1,ngrid 47 47 do islope = 1,nslope 48 if ( h2o_ice(i,islope) > 0.) surf_now = surf_now+ cell_area(i)*subslope_dist(i,islope)48 if (ini_h2oice_sublim(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope) 49 49 enddo 50 50 enddo 51 51 52 52 ! Check of the criterion 53 if ( surf_now < surf_ini*(1. - h2o_ice_crit)) then53 if (h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)) then 54 54 stopPEM = 1 55 55 write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold" 56 write(*,*) " surf_now < surf_ini*(1. - h2o_ice_crit)", surf_now < surf_ini*(1. - h2o_ice_crit)57 write(*,*) "Current surface of h2o ice sublimating =", surf_now58 write(*,*) "Initial surface of h2o ice sublimating =", surf_ini56 write(*,*) "h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)", h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit) 57 write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf 58 write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf 59 59 write(*,*) "Percentage of change accepted =", h2o_ice_crit*100 60 else if ( surf_now > surf_ini*(1. + h2o_ice_crit)) then60 else if (h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)) then 61 61 stopPEM = 1 62 62 write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold" 63 write(*,*) " surf_now > surf_ini*(1. + h2o_ice_crit)", surf_now > surf_ini*(1. + h2o_ice_crit)64 write(*,*) "Current surface of h2o ice sublimating =", surf_now65 write(*,*) "Initial surface of h2o ice sublimating =", surf_ini63 write(*,*) "h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)", h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit) 64 write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf 65 write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf 66 66 write(*,*) "Percentage of change accepted =", h2o_ice_crit*100 67 67 endif 68 68 69 if (abs( surf_ini) < 1.e-5) stopPEM = 069 if (abs(h2oice_ini_surf) < 1.e-5) stopPEM = 0 70 70 71 71 END SUBROUTINE stopping_crit_h2o_ice … … 73 73 !======================================================================= 74 74 75 SUBROUTINE stopping_crit_co2(cell_area, surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)75 SUBROUTINE stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope) 76 76 77 77 use time_evol_mod, only: co2_ice_crit, ps_criterion … … 85 85 ! 86 86 !======================================================================= 87 88 ! arguments:89 ! ---------- 90 ! INPUT 91 integer, intent(in) :: ngrid, nslope ! # of grid physical grid points 92 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells93 real, dimension(ngrid,nslope), intent(in) :: co2_ice ! Actual density of h2o ice 94 real, intent(in) :: surf_ini ! Initial surface of co2 ice that was sublimating95 real, intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files96 real, intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations97 ! OUTPUT87 ! Inputs 88 !------- 89 integer, intent(in) :: ngrid, nslope ! # of grid physical grid points 90 real, dimension(ngrid), intent(in) :: cell_area ! Area of the cells 91 real, dimension(ngrid,nslope), intent(in) :: co2_ice ! Actual density of co2 ice 92 real, intent(in) :: co2ice_ini_surf ! Initial surface of sublimatingco2 ice 93 logical, dimension(ngrid,nslope), intent(in) :: ini_co2ice_sublim ! Grid points where co2 ice was initially sublimating 94 real, intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files 95 real, intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations 96 ! Outputs 97 !-------- 98 98 integer, intent(inout) :: stopPEM ! Stopping criterion code 99 99 100 ! local:101 ! 102 integer :: i, islope ! Loop103 real :: surf_now! Current surface of co2 ice100 ! Locals 101 ! ------ 102 integer :: i, islope ! Loop 103 real :: co2ice_now_surf ! Current surface of co2 ice 104 104 105 105 !======================================================================= 106 ! Computation of the present surface of co2 ice 107 surf_now= 0.106 ! Computation of the present surface of co2 ice still sublimating 107 co2ice_now_surf = 0. 108 108 do i = 1,ngrid 109 109 do islope = 1,nslope 110 if ( co2_ice(i,islope) > 0.) surf_now = surf_now+ cell_area(i)*subslope_dist(i,islope)110 if (ini_co2ice_sublim(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope) 111 111 enddo 112 112 enddo 113 113 114 114 ! Check of the criterion 115 if ( surf_now < surf_ini*(1. - co2_ice_crit)) then115 if (co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)) then 116 116 stopPEM = 3 117 117 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold" 118 write(*,*) " surf_now < surf_ini*(1. - co2_ice_crit)", surf_now < surf_ini*(1. - co2_ice_crit)119 write(*,*) "Current surface of co2 ice sublimating =", surf_now120 write(*,*) "Initial surface of co2 ice sublimating =", surf_ini118 write(*,*) "co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit) 119 write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf 120 write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf 121 121 write(*,*) "Percentage of change accepted =", co2_ice_crit*100. 122 else if ( surf_now > surf_ini*(1. + co2_ice_crit)) then122 else if (co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)) then 123 123 stopPEM = 3 124 124 write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold" 125 write(*,*) " surf_now > surf_ini*(1. + co2_ice_crit)", surf_now > surf_ini*(1. + co2_ice_crit)126 write(*,*) "Current surface of co2 ice sublimating =", surf_now127 write(*,*) "Initial surface of co2 ice sublimating =", surf_ini125 write(*,*) "co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit) 126 write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf 127 write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf 128 128 write(*,*) "Percentage of change accepted =", co2_ice_crit*100. 129 129 endif 130 130 131 if (abs( surf_ini) < 1.e-5) stopPEM = 0131 if (abs(co2ice_ini_surf) < 1.e-5) stopPEM = 0 132 132 133 133 if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
Note: See TracChangeset
for help on using the changeset viewer.
