Changeset 3456 for trunk/LMDZ.COMMON
- Timestamp:
- Oct 11, 2024, 10:37:48 AM (6 weeks ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
r3264 r3456 1 module adsorption_mod 2 3 implicit none 4 5 LOGICAL adsorption_pem ! True by default, to compute adsorption/desorption. Read in pem.def 6 real, save, allocatable :: co2_adsorbded_phys(:,:,:) ! co2 that is in the regolith [kg/m^2] 7 real, save, allocatable :: h2o_adsorbded_phys(:,:,:) ! h2o that is in the regolith [kg/m^2] 8 9 contains 1 MODULE adsorption_mod 2 3 implicit none 4 5 logical :: adsorption_pem ! True by default, to compute adsorption/desorption. Read in pem.def 6 real, save, allocatable, dimension(:,:,:) :: co2_adsorbded_phys ! co2 that is in the regolith [kg/m^2] 7 real, save, allocatable, dimension(:,:,:) :: h2o_adsorbded_phys ! h2o that is in the regolith [kg/m^2] 8 9 !======================================================================= 10 contains 11 !======================================================================= 10 12 11 13 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 16 18 !!! 17 19 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 18 subroutine ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) 19 20 implicit none 21 integer,intent(in) :: ngrid ! number of atmospheric columns 22 integer,intent(in) :: nslope ! number of slope within a mesh 23 integer,intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM 24 allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope)) 25 allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope)) 26 end subroutine ini_adsorption_h_PEM 27 28 !!! ----------------------------------------------- 29 30 subroutine end_adsorption_h_PEM 31 32 implicit none 33 if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys) 34 if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys) 35 end subroutine end_adsorption_h_PEM 36 37 !!! ----------------------------------------------- 38 39 subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, & 40 m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg) 41 42 ! inputs 43 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension: physics x subslope x soil x timeseries 44 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers [1] 45 REAL,INTENT(IN) :: waterice(ngrid,nslope) ! water ice at the surface [kg/m^2] 46 REAL,INTENT(IN) :: co2ice(ngrid,nslope) ! co2 ice at the surface [kg/m^2] 47 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal inertia (J/K/^2/s^1/2) 48 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) ! Soil temperature (K) 49 REAL,INTENT(IN) :: ps(ngrid,timelen) ! Average surface pressure [Pa] 50 REAL,INTENT(IN) :: q_co2(ngrid,timelen) ! Mass mixing ratio of co2 in the first layer (kg/kg) 51 REAL,INTENT(IN) :: q_h2o(ngrid,timelen) ! Mass mixing ratio of H2o in the first layer (kg/kg) 52 53 ! outputs 54 REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope) 55 REAL,INTENT(OUT) :: delta_mh2oreg(ngrid) ! Difference density of h2o adsorbed (kg/m^3) 56 57 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope) ! Density of co2 adsorbed (kg/m^3) 58 REAL,INTENT(OUT) :: delta_mco2reg(ngrid) ! Difference density of co2 adsorbed (kg/m^3) 59 60 ! local variables 61 REAL :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules 20 SUBROUTINE ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) 21 22 implicit none 23 24 integer, intent(in) :: ngrid ! number of atmospheric columns 25 integer, intent(in) :: nslope ! number of slope within a mesh 26 integer, intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM 27 28 allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope)) 29 allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope)) 30 31 END SUBROUTINE ini_adsorption_h_PEM 32 33 !======================================================================= 34 SUBROUTINE end_adsorption_h_PEM 35 36 if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys) 37 if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys) 38 39 END SUBROUTINE end_adsorption_h_PEM 40 41 !======================================================================= 42 SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, & 43 m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg) 44 45 implicit none 46 47 ! Inputs 48 integer, intent(in) :: ngrid, nslope, nsoil_PEM, timelen ! size dimension: physics x subslope x soil x timeseries 49 real, dimension(ngrid,nslope), intent(in) :: tend_h2oglaciers, tend_co2glaciers ! tendancies on the glaciers [1] 50 real, dimension(ngrid,nslope), intent(in) :: waterice ! water ice at the surface [kg/m^2] 51 real, dimension(ngrid,nslope), intent(in) :: co2ice ! co2 ice at the surface [kg/m^2] 52 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM ! Soil Thermal inertia (J/K/^2/s^1/2) 53 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM ! Soil temperature (K) 54 real, dimension(ngrid,timelen), intent(in) :: ps ! Average surface pressure [Pa] 55 real, dimension(ngrid,timelen), intent(in) :: q_co2 ! Mass mixing ratio of co2 in the first layer (kg/kg) 56 real, dimension(ngrid,timelen), intent(in) :: q_h2o ! Mass mixing ratio of H2o in the first layer (kg/kg) 57 58 ! Outputs 59 real, dimension(ngrid), intent(out) :: delta_mh2oreg ! Difference density of h2o adsorbed (kg/m^3) 60 real, dimension(ngrid), intent(out) :: delta_mco2reg ! Difference density of co2 adsorbed (kg/m^3) 61 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3) 62 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3) 63 64 ! Local variables 65 real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbded ! Fraction of the pores occupied by H2O molecules 62 66 ! ------------- 63 64 65 67 ! Compute H2O adsorption, then CO2 adsorption 66 67 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 68 theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg) 69 70 71 call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, & 72 tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg) 73 74 RETURN 75 end subroutine 76 77 !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 78 79 subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 80 theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg) 68 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 69 theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg) 70 call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, & 71 tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg) 72 73 END SUBROUTINE 74 75 !======================================================================= 76 SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 77 theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg) 81 78 82 79 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 90 87 use comsoil_h_PEM, only: layer_PEM, index_breccia 91 88 use comslope_mod, only: subslope_dist, def_slope_mean 92 use vertical_layers_mod, only: ap, bp89 use vertical_layers_mod, only: ap, bp 93 90 use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith 94 91 … … 99 96 #endif 100 97 101 implicit none 102 ! inputs 103 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension 104 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) ! Tendencies on the glaciers () 105 REAL,INTENT(IN) :: waterice(ngrid,nslope) ! Water ice at the surface [kg/m^2] 106 REAL,INTENT(IN) :: co2ice(ngrid,nslope) ! CO2 ice at the surface [kg/m^2] 107 REAL,INTENT(IN) :: ps(ngrid,timelen) ! Surface pressure (Pa) 108 REAL,INTENT(IN) :: q_co2(ngrid,timelen) ! Mass mixing ratio of co2 in the first layer (kg/kg) 109 REAL,INTENT(IN) :: q_h2o(ngrid,timelen) ! Mass mixing ratio of H2o in the first layer (kg/kg) 110 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Soil Thermal inertia (J/K/^2/s^1/2) 111 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) ! Soil temperature (K) 112 113 ! outputs 114 REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope) 115 REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules 116 REAL,INTENT(OUT) :: delta_mreg(ngrid) ! Difference density of h2o adsorbed (kg/m^3) 117 118 ! constants 119 REAL :: Ko = 1.57e-8 ! Jackosky et al. 1997 120 REAL :: e = 2573.9 ! Jackosky et al. 1997 121 REAL :: mu = 0.48 ! Jackosky et al. 1997 122 real :: m_theta = 2.84e-7 ! Mass of h2o per m^2 absorbed Jackosky et al. 1997 123 ! real :: as = 18.9e3 ! Specific area, Buhler & Piqueux 2021 124 real :: as = 9.48e4 ! Specific area, Zent 125 real :: inertie_thresold = 800. ! TI > 800 means cementation 126 127 ! local variables 128 REAL :: deltam_reg_complete(ngrid,index_breccia,nslope) ! Difference in the mass per slope and soil layer (kg/m^3) 129 real :: K ! Used to compute theta 130 integer ig, iloop, islope, it ! For loops 131 INTEGER :: ispermanent_co2glaciers(ngrid,nslope) ! Check if the co2 glacier is permanent 132 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope) ! Check if the h2o glacier is permanent 133 REAL :: deltam_reg_slope(ngrid,nslope) ! Difference density of h2o adsorbed per slope (kg/m^3) 134 REAL :: dm_h2o_regolith_slope(ngrid,nsoil_PEM,nslope) ! Elementary h2o mass adsorded per mesh per slope 135 real :: A,B ! Used to compute the mean mass above the surface 136 real :: p_sat ! Saturated vapor pressure of ice 137 real,allocatable :: mass_mean(:,:) ! Mean mass above the surface 138 real,allocatable :: zplev_mean(:,:) ! Pressure above the surface 139 real,allocatable :: pvapor(:,:) ! Partial pressure above the surface 140 real, allocatable :: pvapor_avg(:) ! Yearly averaged 98 implicit none 99 100 ! Inputs 101 integer, intent(in) :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension 102 real, dimension(ngrid,nslope), intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers () 103 real, dimension(ngrid,nslope), intent(in) :: waterice ! Water ice at the surface [kg/m^2] 104 real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice at the surface [kg/m^2] 105 real, dimension(ngrid,timelen), intent(in) :: ps ! Surface pressure (Pa) 106 real, dimension(ngrid,timelen), intent(in) :: q_co2 ! Mass mixing ratio of co2 in the first layer (kg/kg) 107 real, dimension(ngrid,timelen), intent(in) :: q_h2o ! Mass mixing ratio of H2o in the first layer (kg/kg) 108 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM ! Soil Thermal inertia (J/K/^2/s^1/2) 109 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM ! Soil temperature (K) 110 111 ! Outputs 112 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope) 113 real, dimension(ngrid,nsoil_PEM,nslope), intent(out) :: theta_h2o_adsorbded ! Fraction of the pores occupied by H2O molecules 114 real, dimension(ngrid), intent(out) :: delta_mreg ! Difference density of h2o adsorbed (kg/m^3) 115 116 ! Constants 117 real :: Ko = 1.57e-8 ! Jackosky et al. 1997 118 real :: e = 2573.9 ! Jackosky et al. 1997 119 real :: mu = 0.48 ! Jackosky et al. 1997 120 real :: m_theta = 2.84e-7 ! Mass of h2o per m^2 absorbed Jackosky et al. 1997 121 ! real :: as = 18.9e3 ! Specific area, Buhler & Piqueux 2021 122 real :: as = 9.48e4 ! Specific area, Zent 123 real :: inertie_thresold = 800. ! TI > 800 means cementation 124 125 ! Local variables 126 real, dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete ! Difference in the mass per slope and soil layer (kg/m^3) 127 real :: K ! Used to compute theta 128 integer :: ig, iloop, islope, it ! For loops 129 logical, dimension(ngrid,nslope) :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent 130 logical, dimension(ngrid,nslope) :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent 131 real, dimension(ngrid,nslope) :: deltam_reg_slope ! Difference density of h2o adsorbed per slope (kg/m^3) 132 real, dimension(ngrid,nsoil_PEM,nslope) :: dm_h2o_regolith_slope ! Elementary h2o mass adsorded per mesh per slope 133 real :: A, B ! Used to compute the mean mass above the surface 134 real :: p_sat ! Saturated vapor pressure of ice 135 real, dimension(:,:), allocatable :: mass_mean ! Mean mass above the surface 136 real, dimension(:,:), allocatable :: zplev_mean ! Pressure above the surface 137 real, dimension(:,:), allocatable :: pvapor ! Partial pressure above the surface 138 real, dimension(:) , allocatable :: pvapor_avg ! Yearly averaged 141 139 142 140 ! 0. Some initializations 143 144 145 allocate(mass_mean(ngrid,timelen)) 146 allocate(zplev_mean(ngrid,timelen)) 147 allocate(pvapor(ngrid,timelen)) 148 allocate(pvapor_avg(ngrid)) 149 A =(1/m_co2 - 1/m_noco2) 150 B=1/m_noco2 151 theta_h2o_adsorbded(:,:,:) = 0. 152 dm_h2o_regolith_slope(:,:,:) = 0. 153 154 #ifndef CPP_STD 155 156 !0.1 Look at perennial ice 157 do ig = 1,ngrid 158 do islope = 1,nslope 159 if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then 160 ispermanent_h2oglaciers(ig,islope) = 1 161 else 162 ispermanent_h2oglaciers(ig,islope) = 0 163 endif 164 165 if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then 166 ispermanent_co2glaciers(ig,islope) = 1 167 else 168 ispermanent_co2glaciers(ig,islope) = 0 169 endif 170 enddo 171 enddo 172 173 ! 0.2 Compute the partial pressure of vapor 174 !a. the molecular mass into the column 175 do ig = 1,ngrid 176 mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B) 177 enddo 178 179 180 ! b. pressure level 181 do it = 1,timelen 182 do ig = 1,ngrid 183 zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) 184 enddo 185 enddo 186 ! c. Vapor pressure 187 pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:) 188 pvapor_avg(:) = sum(pvapor(:,:),2)/timelen 189 #endif 190 deallocate(pvapor) 191 deallocate(zplev_mean) 192 deallocate(mass_mean) 193 #ifndef CPP_STD 194 195 ! 1. we compute the mass of H2O adsorded in each layer of the meshes 196 197 do ig = 1,ngrid 198 do islope = 1,nslope 199 do iloop = 1,index_breccia 200 K = Ko*exp(e/tsoil_PEM(ig,iloop,islope)) 201 if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then 202 theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu 203 else 204 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 ... 205 theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu 206 endif 207 dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith 208 enddo 209 enddo 210 enddo 211 212 ! 2. Check the exchange between the atmosphere and the regolith 213 214 do ig = 1,ngrid 215 delta_mreg(ig) = 0. 216 do islope = 1,nslope 217 deltam_reg_slope(ig,islope) = 0. 218 do iloop = 1,index_breccia 219 if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then 220 if(iloop==1) then 221 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) & 222 *(layer_PEM(iloop)) 223 else 224 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) & 225 *(layer_PEM(iloop) - layer_PEM(iloop-1)) 226 endif 227 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! 228 deltam_reg_complete(ig,iloop,islope) = 0. 229 endif 230 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) 231 enddo 232 delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 233 enddo 234 enddo 235 m_h2o_completesoil(:,:,:) = dm_h2o_regolith_slope(:,:,:) 236 237 RETURN 238 #endif 239 end subroutine 240 241 !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 141 allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pvapor(ngrid,timelen),pvapor_avg(ngrid)) 142 A = 1./m_co2 - 1./m_noco2 143 B = 1./m_noco2 144 theta_h2o_adsorbded = 0. 145 dm_h2o_regolith_slope = 0. 146 ispermanent_h2oglaciers = .false. 147 ispermanent_co2glaciers = .false. 148 149 #ifndef CPP_STD 150 ! 0.1 Look at perennial ice 151 do ig = 1,ngrid 152 do islope = 1,nslope 153 if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true. 154 if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true. 155 enddo 156 enddo 157 158 ! 0.2 Compute the partial pressure of vapor 159 ! a. the molecular mass into the column 160 do ig = 1,ngrid 161 mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B) 162 enddo 163 164 ! b. pressure level 165 do it = 1,timelen 166 do ig = 1,ngrid 167 zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) 168 enddo 169 enddo 170 171 ! c. Vapor pressure 172 pvapor = mass_mean/m_h2o*q_h2o*zplev_mean 173 pvapor_avg = sum(pvapor,2)/timelen 174 #endif 175 deallocate(pvapor,zplev_mean,mass_mean) 176 177 #ifndef CPP_STD 178 ! 1. we compute the mass of H2O adsorded in each layer of the meshes 179 do ig = 1,ngrid 180 do islope = 1,nslope 181 do iloop = 1,index_breccia 182 K = Ko*exp(e/tsoil_PEM(ig,iloop,islope)) 183 if (TI_PEM(ig,iloop,islope) < inertie_thresold) then 184 theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1. + K*pvapor_avg(ig)))**mu 185 else 186 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 ... 187 theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1. + K*p_sat))**mu 188 endif 189 dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith 190 enddo 191 enddo 192 enddo 193 194 ! 2. Check the exchange between the atmosphere and the regolith 195 do ig = 1,ngrid 196 delta_mreg(ig) = 0. 197 do islope = 1,nslope 198 deltam_reg_slope(ig,islope) = 0. 199 do iloop = 1,index_breccia 200 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then 201 if (iloop == 1) then 202 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop)) 203 else 204 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1)) 205 endif 206 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! 207 deltam_reg_complete(ig,iloop,islope) = 0. 208 endif 209 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) 210 enddo 211 delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 212 enddo 213 enddo 214 m_h2o_completesoil = dm_h2o_regolith_slope 215 #endif 216 END SUBROUTINE 217 218 !======================================================================= 242 219 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 243 220 !!! … … 247 224 !!! 248 225 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 249 250 SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,& 251 tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg) 226 SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg) 252 227 253 228 use comsoil_h_PEM, only: layer_PEM, index_breccia, index_breccia … … 262 237 #endif 263 238 264 IMPLICIT NONE 265 ! Inputs: 266 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension 267 REAL,INTENT(IN) :: ps(ngrid,timelen) ! Average surface pressure [Pa] 268 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) ! Mean Soil Temperature [K] 269 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope) ! Mean Thermal Inertia [USI] 270 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) ! Tendencies on the glaciers () 271 REAL,INTENT(IN) :: q_co2(ngrid,timelen),q_h2o(ngrid,timelen) ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg) 272 REAL,INTENT(IN) :: waterice(ngrid,nslope) ! Water ice at the surface [kg/m^2] 273 REAL,INTENT(IN) :: co2ice(ngrid,nslope) ! CO2 ice at the surface [kg/m^2] 239 implicit none 240 241 ! Inputs: 242 integer, intent(in) :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension 243 real, dimension(ngrid,timelen), intent(in) :: ps ! Average surface pressure [Pa] 244 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM ! Mean Soil Temperature [K] 245 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM ! Mean Thermal Inertia [USI] 246 real, dimension(ngrid,nslope), intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers () 247 real, dimension(ngrid,timelen), intent(in) :: q_co2, q_h2o ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg) 248 real, dimension(ngrid,nslope), intent(in) :: waterice ! Water ice at the surface [kg/m^2] 249 real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice at the surface [kg/m^2] 274 250 275 251 ! Outputs: 276 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)! Density of co2 adsorbed (kg/m^3)277 REAL,INTENT(OUT) :: delta_mreg(ngrid)! Difference density of co2 adsorbed (kg/m^3)278 252 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3) 253 real, dimension(ngrid), intent(out) :: delta_mreg ! Difference density of co2 adsorbed (kg/m^3) 254 279 255 ! Constants: 280 281 REAL :: alpha = 7.512e-6 ! Zent & Quinn 1995 282 REAL :: beta = -1541.5 ! Zent & Quinn 1995 283 REAL :: inertie_thresold = 800. ! TI > 800 means cementation 284 real :: m_theta = 4.27e-7 ! Mass of co2 per m^2 absorbed 256 real :: alpha = 7.512e-6 ! Zent & Quinn 1995 257 real :: beta = -1541.5 ! Zent & Quinn 1995 258 real :: inertie_thresold = 800. ! TI > 800 means cementation 259 real :: m_theta = 4.27e-7 ! Mass of co2 per m^2 absorbed 285 260 ! real :: as = 18.9e3 ! Specific area, Buhler & Piqueux 2021 286 real :: as = 9.48e4 ! Same as previous but from zent 287 ! Local 288 real :: A,B ! Used to compute the mean mass above the surface 289 INTEGER :: ig,islope,iloop,it ! For loops 290 REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope) ! Elementary mass adsorded per mesh per slope 291 INTEGER :: ispermanent_co2glaciers(ngrid,nslope) ! Check if the co2 glacier is permanent 292 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope) ! Check if the h2o glacier is permanent 293 REAL :: deltam_reg_complete(ngrid,index_breccia,nslope) ! Difference in the mass per slope and soil layer (kg/m^3) 294 REAL :: deltam_reg_slope(ngrid,nslope) ! Difference in the mass per slope (kg/m^3) 295 REAL :: m_h2o_adsorbed(ngrid,nsoil_PEM,nslope) ! Density of CO2 adsorbed (kg/m^3) 296 REAL :: theta_h2o_adsorbed(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules 297 REAL :: delta_mh2o(ngrid) ! Difference density of h2o adsorbed (kg/m^3) 298 !timelen array are allocated because heavy ... 299 real,allocatable :: mass_mean(:,:) ! Mean mass above the surface 300 real,allocatable :: zplev_mean(:,:) ! Pressure above the surface 301 real,allocatable :: pco2(:,:) ! Partial pressure above the surface 302 real, allocatable :: pco2_avg(:) ! Yearly averaged 261 real :: as = 9.48e4 ! Same as previous but from zent 262 263 ! Local 264 real :: A, B ! Used to compute the mean mass above the surface 265 integer :: ig, islope, iloop, it ! For loops 266 real, dimension(ngrid,nsoil_PEM,nslope) :: dm_co2_regolith_slope ! Elementary mass adsorded per mesh per slope 267 logical, dimension(ngrid,nslope) :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent 268 logical, dimension(ngrid,nslope) :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent 269 real, dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete ! Difference in the mass per slope and soil layer (kg/m^3) 270 real, dimension(ngrid,nslope) :: deltam_reg_slope ! Difference in the mass per slope (kg/m^3) 271 real, dimension(ngrid,nsoil_PEM,nslope) :: m_h2o_adsorbed ! Density of CO2 adsorbed (kg/m^3) 272 real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules 273 real, dimension(ngrid) :: delta_mh2o ! Difference density of h2o adsorbed (kg/m^3) 274 ! timelen array are allocated because heavy... 275 real, dimension(:,:), allocatable :: mass_mean ! Mean mass above the surface 276 real, dimension(:,:), allocatable :: zplev_mean ! Pressure above the surface 277 real, dimension(:,:), allocatable :: pco2 ! Partial pressure above the surface 278 real, dimension(:), allocatable :: pco2_avg ! Yearly averaged 303 279 304 280 ! 0. Some initializations 305 306 allocate(mass_mean(ngrid,timelen)) 307 allocate(zplev_mean(ngrid,timelen)) 308 allocate(pco2(ngrid,timelen)) 309 allocate(pco2_avg(ngrid)) 310 311 312 313 m_h2o_adsorbed(:,:,:) = 0. 314 A =(1/m_co2 - 1/m_noco2) 315 B=1/m_noco2 316 317 dm_co2_regolith_slope(:,:,:) = 0 318 delta_mreg(:) = 0. 319 320 #ifndef CPP_STD 321 322 !0.1 Look at perennial ice 323 do ig = 1,ngrid 324 do islope = 1,nslope 325 if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then 326 ispermanent_h2oglaciers(ig,islope) = 1 327 else 328 ispermanent_h2oglaciers(ig,islope) = 0 329 endif 330 331 if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then 332 ispermanent_co2glaciers(ig,islope) = 1 333 else 334 ispermanent_co2glaciers(ig,islope) = 0 335 endif 336 enddo 337 enddo 338 339 ! 0.2 Compute the partial pressure of CO2 340 !a. the molecular mass into the column 341 do ig = 1,ngrid 342 mass_mean(ig,:) = 1./(A*q_co2(ig,:) +B) 343 enddo 344 345 ! b. pressure level 346 do it = 1,timelen 347 do ig = 1,ngrid 348 zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) 349 enddo 350 enddo 351 352 ! c. Vapor pressure 353 pco2(:,:) = mass_mean(:,:)/m_co2*q_co2(:,:)*zplev_mean(:,:) 354 pco2_avg(:) = sum(pco2(:,:),2)/timelen 355 356 deallocate(zplev_mean) 357 deallocate(mass_mean) 358 deallocate(pco2) 359 360 361 ! 1. Compute the fraction of the pores occupied by H2O 362 363 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 364 theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o) 365 366 ! 2. we compute the mass of co2 adsorded in each layer of the meshes 367 368 do ig = 1,ngrid 369 do islope = 1,nslope 370 do iloop = 1,index_breccia 371 if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then 372 dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ & 373 (alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 374 else 375 if(abs(m_co2_completesoil(ig,iloop,islope)).lt.(1e-10)) then !!! we are at first call 376 dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) & 377 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 378 else ! no change: permanent ice stick the atoms of CO2 379 dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope) 380 endif 381 endif 382 enddo 383 enddo 384 enddo 385 386 ! 3. Check the exchange between the atmosphere and the regolith 387 388 do ig = 1,ngrid 389 delta_mreg(ig) = 0. 390 do islope = 1,nslope 391 deltam_reg_slope(ig,islope) = 0. 392 do iloop = 1,index_breccia 393 if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then 394 if(iloop == 1) then 395 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) & 396 *(layer_PEM(iloop)) 397 else 398 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) & 399 *(layer_PEM(iloop) - layer_PEM(iloop-1)) 400 endif 401 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! 402 deltam_reg_complete(ig,iloop,islope) = 0. 403 endif 404 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) 405 enddo 406 delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 407 enddo 408 enddo 409 m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:) 410 411 !======================================================================= 412 RETURN 413 #endif 414 END 415 416 417 end module 281 allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pco2(ngrid,timelen),pco2_avg(ngrid)) 282 m_h2o_adsorbed = 0. 283 A = 1./m_co2 - 1./m_noco2 284 B = 1./m_noco2 285 dm_co2_regolith_slope = 0. 286 delta_mreg = 0. 287 ispermanent_h2oglaciers = .false. 288 ispermanent_co2glaciers = .false. 289 290 #ifndef CPP_STD 291 ! 0.1 Look at perennial ice 292 do ig = 1,ngrid 293 do islope = 1,nslope 294 if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true. 295 if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true. 296 enddo 297 enddo 298 299 ! 0.2 Compute the partial pressure of CO2 300 ! a. the molecular mass into the column 301 do ig = 1,ngrid 302 mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B) 303 enddo 304 305 ! b. pressure level 306 do it = 1,timelen 307 do ig = 1,ngrid 308 zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) 309 enddo 310 enddo 311 312 ! c. Vapor pressure 313 pco2 = mass_mean/m_co2*q_co2*zplev_mean 314 pco2_avg(:) = sum(pco2(:,:),2)/timelen 315 316 deallocate(zplev_mean,mass_mean,pco2) 317 318 ! 1. Compute the fraction of the pores occupied by H2O 319 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, & 320 theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o) 321 322 ! 2. we compute the mass of co2 adsorded in each layer of the meshes 323 do ig = 1,ngrid 324 do islope = 1,nslope 325 do iloop = 1,index_breccia 326 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then 327 dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ & 328 (alpha*pco2_avg(ig) + sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 329 else 330 if (abs(m_co2_completesoil(ig,iloop,islope)) < 1.e-10) then !!! we are at first call 331 dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) & 332 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 333 else ! no change: permanent ice stick the atoms of CO2 334 dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope) 335 endif 336 endif 337 enddo 338 enddo 339 enddo 340 341 ! 3. Check the exchange between the atmosphere and the regolith 342 do ig = 1,ngrid 343 delta_mreg(ig) = 0. 344 do islope = 1,nslope 345 deltam_reg_slope(ig,islope) = 0. 346 do iloop = 1,index_breccia 347 if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then 348 if (iloop == 1) then 349 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop)) 350 else 351 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1)) 352 endif 353 else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC! 354 deltam_reg_complete(ig,iloop,islope) = 0. 355 endif 356 deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope) 357 enddo 358 delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 359 enddo 360 enddo 361 m_co2_completesoil = dm_co2_regolith_slope 362 #endif 363 364 END SUBROUTINE regolith_co2adsorption 365 366 END MODULE adsorption_mod 367 -
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3446 r3456 436 436 == 30/09/2024 == JBC 437 437 Correction of the launching script for a relaunch situation. 438 439 == 11/10/2024 == JBC 440 Cleaning of the adsorption module to make the debugging easier. -
trunk/LMDZ.COMMON/libf/evolution/deftank/multiple_exec.sh
r3354 r3456 1 1 #!/bin/bash 2 ############################################################ #3 ### Script to exec tute mutliple scripts in subdirectories ###4 ############################################################ #2 ############################################################ 3 ### Script to execute multiple scripts in subdirectories ### 4 ############################################################ 5 5 6 6 # Name of the file to execute in all subdirectories -
trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
r3430 r3456 217 217 ! b. Special case for inertiedat, inertiedat_PEM 218 218 call get_field("inertiedat_PEM",inertiedat_PEM,found) 219 if (.not. found) then219 if (.not. found) then 220 220 do iloop = 1,nsoil_PCM 221 221 inertiedat_PEM(:,iloop) = inertiedat(:,iloop) … … 259 259 write(num,fmt='(i2.2)') islope 260 260 call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found) 261 if (.not. found) then261 if (.not. found) then 262 262 write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>' 263 263 write(*,*)'will reconstruct the values of Tsoil' … … 346 346 347 347 if (.not. found) deltam_co2_regolith_phys = 0. 348 if (.not. found2) deltam_h2o_regolith_phys = 0.348 if (.not. found2) deltam_h2o_regolith_phys = 0. 349 349 350 350 write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
Note: See TracChangeset
for help on using the changeset viewer.