Index: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 (revision 2944) @@ -93,4 +93,7 @@ use comslope_mod, only : subslope_dist,def_slope_mean use vertical_layers_mod, ONLY: ap,bp + use constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, & + m_h2o,m_co2,m_noco2, & + rho_regolith #endif @@ -121,11 +124,6 @@ real :: as_breccia = 1.7e4 ! Specific area of basalt, Zent real :: inertie_thresold = 800. ! TI > 800 means cementation - real :: m_h2o = 18.01528E-3 ! Molecular weight of h2o (kg/mol) - real :: m_co2 = 44.01E-3 ! Molecular weight of co2 (kg/mol) - real :: m_noco2 = 33.37E-3 ! Molecular weight of non co2 (kg/mol) - real :: rho_regolith = 2000. ! density of the regolith (Zent 1995, Buhler & piqueux 2021) - real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005 - real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005 -! local variables + + ! local variables #ifndef CPP_STD REAL :: deltam_reg_complete(ngrid,index_breccia,nslope) ! Difference in the mass per slope and soil layer (kg/m^3) @@ -202,5 +200,5 @@ theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu else - p_sat =exp(alpha_clapeyron/tsoil_PEM(ig,iloop,islope) +beta_clapeyron) ! we assume fixed temperature in the ice ... not really:q good but ... + p_sat =exp(beta_clap_h2o/tsoil_PEM(ig,iloop,islope) +alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ... theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu endif @@ -250,4 +248,5 @@ use comslope_mod, only : subslope_dist,def_slope_mean use vertical_layers_mod, ONLY: ap,bp + use constants_marspem_mod, only: m_co2, m_noco2,rho_regolith #endif @@ -272,7 +271,4 @@ REAL :: beta = -1541.5 ! Zent & Quinn 1995 REAL :: inertie_thresold = 800. ! TI > 800 means cementation - REAL :: rho_regolith = 2000. ! density of the regolith, buhler & piqueux 2021 - real :: m_co2 = 44.01E-3 ! Molecular weight of co2 (kg/mol) - real :: m_noco2 = 33.37E-3 ! Molecular weight of h2o (kg/mol) real :: m_theta = 4.27e-7 ! Mass of co2 per m^2 absorbed ! real :: as = 18.9e3 ! Specific area, Buhler & Piqueux 2021 Index: trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90 (revision 2944) @@ -64,4 +64,5 @@ USE comconst_mod, ONLY: pi + use comcstfi_h, only: g !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -87,5 +88,4 @@ REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m] ! Local - REAL, PARAMETER :: g = 3.71 ! surface gravity [m/s^2] INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000 REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol] @@ -226,5 +226,9 @@ !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -implicit none + +use constants_marspem_mod,only : alpha_clap_co2,beta_clap_co2 + +implicit none + ! arguments: ! ---------- @@ -243,7 +247,4 @@ INTEGER :: ig,it ! for loop REAL :: ave ! intermediate to compute average - REAL :: alpha_clap, beta_clap ! Clapeyron law for CO2 - alpha_clap = 23.3494 ! James et al. 1992 - beta_clap = 3182.48 ! James et al. 1992 !!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -253,5 +254,5 @@ ave = 0 DO it = 1,timelen - ave = ave + beta_clap/(alpha_clap-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) + ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100)) ENDDO Tcond(ig) = ave/timelen Index: trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90 (revision 2939) +++ (revision ) @@ -1,55 +1,0 @@ - SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute the ice table depth knowing the yearly average water -!!! density at the surface and at depth. -!!! Computations are made following the methods in Schorgofer et al., 2005 -!!! -!!! Author: LL -!!! -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -#ifndef CPP_STD - USE comsoil_h_PEM, only: layer_PEM ! Depth of the vertical grid - implicit none - - integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM - real,intent(in) :: rhowatersurf_ave(ngrid,nslope) ! Water density at the surface, yearly averaged [kg/m^3] - real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope) ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3] - real,intent(inout) :: ice_table(ngrid,nslope) ! ice table depth [m] - real :: z1,z2 ! intermediate variables used when doing a linear interpolation between two depths to find the root - integer ig, islope,isoil ! loop variables - real :: diff_rho(nsoil_PEM) ! difference of water vapor density between the surface and at depth [kg/m^3] - - - do ig = 1,ngrid - do islope = 1,nslope - ice_table(ig,islope) = -10. - - do isoil = 1,nsoil_PEM - diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) - - enddo - if (ig.eq.400) then - write(*,*) 'ig,islope,diff',ig,islope,diff_rho(:) - write(*,*) 'rho surf',rhowatersurf_ave(ig,islope) - write(*,*) 'rho soil',rhowatersoil_ave(ig,:,islope) - endif - if(diff_rho(1) > 0) then ! ice is at the surface - ice_table(ig,islope) = 0. - else - do isoil = 1,nsoil_PEM -1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root - if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then - z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1)) - z2 = -layer_PEM(isoil+1)*z1 + diff_rho(isoil+1) - ice_table(ig,islope) = -z2/z1 - exit - endif - enddo - endif - enddo - enddo -!======================================================================= - RETURN -#endif - END Index: trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 (revision 2944) +++ trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90 (revision 2944) @@ -0,0 +1,36 @@ + MODULE constants_marspem_mod + + + IMPLICIT NONE + + +! Molecular masses for CO2,H2O and non condensible gaz, following Franz et al. 2017 + REAL,PARAMETER :: m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) + REAL,PARAMETER :: m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) + REAL,PARAMETER :: m_h2o = 18.01528E-3 ! Molecular weight of h2o (kg/mol) + +! Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992 + REAL,PARAMETER :: alpha_clap_co2 = 23.3494 !Uniteless,James et al. 1992 + REAL,PARAMETER :: beta_clap_co2 = 3182.48 !Kelvin, James et al. 1992 + +! Coefficient for Clapeyron law for psat(psat = exp(beta/Th2o+alpha)),following Murphie and Kood 2005 + REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphie and Kood 2005 + REAL,PARAMETER :: beta_clap_h2o = -6143.7 ! Kelvin, Murphie and Kood 2005 + +! Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021) + REAL,PARAMETER :: rho_regolith = 2000. ! kg/m^3 + +! Average Thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008 + REAL,PARAMETER :: TI_regolith_avg = 250. ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI] + REAL,PARAMETER :: TI_breccia = 750. ! Thermal inertia of Breccia following Wood 2009 [SI] + REAL,PARAMETER :: TI_bedrock = 2300. ! Thermal inertia of Bedrock following Wood 2009 [SI] + +! Stefan Boltzmann constant + REAL,PARAMETER :: sigmaB=5.678E-8 + +! Latent heat of CO2 + REAL,PARAMETER :: Lco2 = 5.71E5 ! Pilorget and Forget 2016 + + + END MODULE constants_marspem_mod + Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2944) @@ -87,4 +87,7 @@ use co2glaciers_mod,only: co2glaciers_evol,co2glaciersflow use criterion_pem_stop_mod,only: criterion_waterice_stop,criterion_co2_stop + use constants_marspem_mod,only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, & + m_co2,m_noco2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, & @@ -179,5 +182,5 @@ INTEGER :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param - real,save :: m_co2, m_noco2, A , B, mmean ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg] + real,save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] real ,allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3] real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM @@ -189,6 +192,4 @@ REAL, ALLOCATABLE :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) REAL :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 - REAL :: beta_clap_co2 = 3182.48 ! clapeyron's law for CO2 - REAL :: alpha_clap_co2 = 23.3494 ! Clapeyron's law for CO2 @@ -238,6 +239,4 @@ REAL :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] REAL :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] - REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K] - REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1] LOGICAL :: bool_sublim ! logical to check if there is sublimation or not @@ -278,7 +277,5 @@ ngrid=ngridmx nlayer=llm - - m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) - m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) + A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 @@ -971,5 +968,5 @@ do ig = 1,ngrid do isoil = 1,nsoilmx_PEM - watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil) + watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil) if(isnan(Tsoil_locslope(ig,isoil))) then call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1) Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 2944) @@ -10,4 +10,6 @@ USE temps_mod_evol, ONLY: year_PEM USE ice_table_mod, only: computeice_table_equilibrium + USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, & + TI_breccia,TI_bedrock #ifndef CPP_STD use surfdat_h, only: watercaptag @@ -53,6 +55,4 @@ LOGICAL :: found2 ! check if variables are found in the start integer :: iloop,ig,islope,it,isoil ! index for loops - REAL :: TI_breccia = 750. ! Thermal inertia of Breccia following Wood 2009 [SI] - REAL :: TI_bedrock = 2300. ! Thermal inertia of Bedrock following Wood 2009 [SI] real :: kcond ! Thermal conductivity, intermediate variable [SI] real :: delta ! Depth of the interface regolith-breccia, breccia -bedrock [m] @@ -64,6 +64,4 @@ real :: beta_tmp(ngrid,nsoil_PEM-1) ! Intermediate for tsoil computatio [] real :: year_PEM_read ! Year of the PEM previous run - real :: alpha_clap_h2o = -6143.7 ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [K^-1] - real :: beta_clap_h2o = 28.9074 ! Intermediate coefficient to compute psat using clapeyron law, Murphie et al. 2005 [1] LOGICAL :: startpem_file ! boolean to check if we read the startfile or not @@ -244,5 +242,5 @@ do isoil = nsoil_GCM+1,nsoil_PEM do ig = 1,ngrid - watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) + watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope) enddo enddo @@ -448,5 +446,5 @@ do isoil = nsoil_GCM+1,nsoil_PEM do ig = 1,ngrid - watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) + watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope) enddo enddo Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2944) @@ -11,5 +11,5 @@ use comsoil_h, only: nsoilmx USE comsoil_h_PEM, ONLY: soil_pem - + use constants_marspem_mod,only: m_co2,m_noco2 IMPLICIT NONE @@ -55,6 +55,6 @@ INTEGER :: edges(4),corner(4) - INTEGER :: i,j,l,t ! loop variables - real,save :: m_co2, m_noco2, A , B, mmean ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer + INTEGER :: i,j,l,t ! loop variables + real :: A , B, mmean ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer INTEGER :: islope ! loop for variables @@ -80,6 +80,4 @@ modname="read_data_gcm" - m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) - m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90 (revision 2944) @@ -3,4 +3,6 @@ ! SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope) + + use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2 IMPLICIT NONE @@ -32,13 +34,8 @@ INTEGER :: i,t,islope - REAL :: eps, sigma, L, beta, alpha, coef, ave + REAL :: eps, coef, ave eps=0.95 - sigma=5.678E-8 - L=5.71*10**5 - beta=3182.48 - alpha=23.3494 - - coef=669*88875*eps*sigma/L + coef=669*88875*eps*sigmaB/Lco2 ! Evolution of the water ice for each physical point @@ -49,6 +46,6 @@ if(co2ice_slope(i,islope).gt.1e-4 .and. abs(tendencies_co2_ice_phys(i,islope)).gt.1e-5) then do t=1,timelen - ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 & - -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4 + ave=ave+(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 & + -(beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4 enddo endif Index: trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 (revision 2939) +++ trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 (revision 2944) @@ -5,4 +5,5 @@ USE vertical_layers_mod, ONLY: ap,bp USE comsoil_h_PEM, only: n_1km + USE constants_marspem_mod,only: TI_breccia,TI_bedrock, TI_regolith_avg implicit none ! Input: @@ -19,9 +20,6 @@ REAL :: inertie_thresold = 800. ! look for ice - REAL :: inertie_averaged = 250 ! Mellon et al. 2000 REAL :: ice_inertia = 1200 ! Inertia of ice REAL :: P610 = 610.0 ! current average pressure on Mars [Pa] - REAL :: TI_breccia = 750. ! THermal inertia of Breccia following Wood 2009 [SI] - REAL :: TI_bedrock = 2300. ! Thermal inertia of Bedrock following Wood 2009 [SI] ! Local variables: @@ -45,5 +43,5 @@ do ig = 1,ngrid if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then - regolith_inertia(ig,islope) = inertie_averaged + regolith_inertia(ig,islope) = TI_regolith_avg endif if(reg_thprop_dependp) then