Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3319) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3320) @@ -286,2 +286,9 @@ == 25/04/2024 == JBC The "start" and "startfi" file names are renamed to match those from the Mars PCM. This is necessary to initialize correctly the 3D PEM with "phys_state_var_init_mod.F90". + +== 26/04/2024 == JBC +- Small corrections to make the PEM work in 3D. +- Addition of the flag "layering_algo" (Default = .false.) defined in the "run_PEM.def" to choose to use the layering algorithm. + +== 30/04/2024 == JBC +Small correction of a bug in the reading of "run_PEM.def" + some cleanings. Index: trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def (revision 3319) +++ trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def (revision 3320) @@ -92,5 +92,5 @@ # co2ice_flow=.true. -!#---------- Layering parameters ----------# +#---------- Layering parameters ----------# # Do you want to run with the layering algorithm? Default = .false. # layering=.true. Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 (revision 3319) +++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90 (revision 3320) @@ -1,5 +1,5 @@ -module ice_table_mod - - implicit none +MODULE ice_table_mod + +implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -10,34 +10,36 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - LOGICAL,SAVE :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium - LOGICAL,SAVE :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) - real,save,allocatable :: porefillingice_depth(:,:) ! ngrid x nslope: Depth of the ice table [m] - real,save,allocatable :: porefillingice_thickness(:,:) ! ngrid x nslope: Thickness of the ice table [m] - +logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium +logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) +real, allocatable, dimension(:,:), save :: porefillingice_depth ! ngrid x nslope: Depth of the ice table [m] +real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m] + +!---------------------------------------- contains - - - subroutine ini_ice_table_porefilling(ngrid,nslope) - - implicit none - integer,intent(in) :: ngrid ! number of atmospheric columns - integer,intent(in) :: nslope ! number of slope within a mesh - - allocate(porefillingice_depth(ngrid,nslope)) - allocate(porefillingice_thickness(ngrid,nslope)) - - end subroutine ini_ice_table_porefilling - - subroutine end_ice_table_porefilling - - implicit none - if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) - if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) - - end subroutine end_ice_table_porefilling - -!!! -------------------------------------- - - SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) +!---------------------------------------- +SUBROUTINE ini_ice_table_porefilling(ngrid,nslope) + +implicit none + +integer, intent(in) :: ngrid ! number of atmospheric columns +integer, intent(in) :: nslope ! number of slope within a mesh + +allocate(porefillingice_depth(ngrid,nslope)) +allocate(porefillingice_thickness(ngrid,nslope)) + +END SUBROUTINE ini_ice_table_porefilling + +!---------------------------------------- +SUBROUTINE end_ice_table_porefilling + +implicit none + +if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) +if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) + + END SUBROUTINE end_ice_table_porefilling + +!---------------------------------------- +SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! @@ -45,108 +47,105 @@ !!! density at the surface and at depth. !!! Computations are made following the methods in Schorgofer et al., 2005 -!!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere -!!! +!!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere +!!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use math_mod,only: findroot - USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM ! Depth of the vertical grid - USE soil_thermalproperties_mod, only: ice_thermal_properties - implicit none -! inputs -! ----------- - integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM - logical,intent(in) :: watercaptag(ngrid) ! Boolean to check the presence of a perennial glacier - 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(in) :: regolith_inertia(ngrid,nslope) ! Thermal inertia of the regolith layer [SI] +use math_mod, only: findroot +use comsoil_h_PEM, only: mlayer_PEM, layer_PEM ! Depth of the vertical grid +use soil_thermalproperties_mod, only: ice_thermal_properties + +implicit none + +! Inputs +! -------- +integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM +logical, dimension(ngrid), intent(in) :: watercaptag ! Boolean to check the presence of a perennial glacier +real, dimension(ngrid,nslope), intent(in) :: rhowatersurf_ave ! Water density at the surface, yearly averaged [kg/m^3] +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] +real, dimension(ngrid,nslope), intent(in) :: regolith_inertia ! Thermal inertia of the regolith layer [SI] ! Ouputs -! ----------- - real,intent(out) :: ice_table_beg(ngrid,nslope) ! ice table depth [m] - real,intent(out) :: ice_table_thickness(ngrid,nslope) ! ice table thickness [m] -! Local -! ----------- - integer ig, islope,isoil,isoilend ! loop variables - real :: diff_rho(nsoil_PEM) ! difference of water vapor density between the surface and at depth [kg/m^3] - real :: ice_table_end ! depth of the end of the ice table [m] - real :: previous_icetable_depth(ngrid,nslope) ! Ice table computed at previous ice depth [m] - real :: stretch ! stretch factor to improve the convergence of the ice table - real :: wice_inertia ! Water Ice thermal Inertia [USI] +! -------- +real, dimension(ngrid,nslope), intent(out) :: ice_table_beg ! ice table depth [m] +real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m] +! Locals +! -------- +integer ig, islope,isoil,isoilend ! loop variables +real :: diff_rho(nsoil_PEM) ! difference of water vapor density between the surface and at depth [kg/m^3] +real :: ice_table_end ! depth of the end of the ice table [m] +real :: previous_icetable_depth(ngrid,nslope) ! Ice table computed at previous ice depth [m] +real :: stretch ! stretch factor to improve the convergence of the ice table +real :: wice_inertia ! Water Ice thermal Inertia [USI] ! Code -! ----------- - +! ------ previous_icetable_depth(:,:) = ice_table_beg(:,:) - do ig = 1,ngrid - if(watercaptag(ig)) then - ice_table_beg(ig,:) = 0. - ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) - else +do ig = 1,ngrid + if (watercaptag(ig)) then + ice_table_beg(ig,:) = 0. + ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) + else do islope = 1,nslope - ice_table_beg(ig,islope) = -1. - ice_table_thickness(ig,islope) = 0. - do isoil = 1,nsoil_PEM - diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) - enddo - if(diff_rho(1) > 0) then ! ice is at the surface - ice_table_beg(ig,islope) = 0. - do isoilend = 2,nsoil_PEM-1 - if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then - call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end) - ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) - exit - endif - enddo - 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 - call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope)) - ! Now let's find the end of the ice table - ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) - do isoilend = isoil+1,nsoil_PEM-1 - if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then - call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end) - ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) - exit - endif - enddo - endif !diff_rho(z) <0 & diff_rho(z+1) > 0 - enddo - endif ! diff_rho(1) > 0 + ice_table_beg(ig,islope) = -1. + ice_table_thickness(ig,islope) = 0. + do isoil = 1,nsoil_PEM + diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) + enddo + if (diff_rho(1) > 0) then ! ice is at the surface + ice_table_beg(ig,islope) = 0. + do isoilend = 2,nsoil_PEM-1 + if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend+1) < 0.) then + call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) + ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) + exit + endif + enddo + 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) < 0 .and. diff_rho(isoil + 1) > 0.) then + call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope)) + ! Now let's find the end of the ice table + ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) + do isoilend = isoil+1,nsoil_PEM-1 + if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then + call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) + ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) + exit + endif + enddo + endif !diff_rho(z) <0 & diff_rho(z+1) > 0 + enddo + endif ! diff_rho(1) > 0 enddo - endif ! watercaptag - enddo + endif ! watercaptag +enddo ! Small trick to speed up the convergence, Oded's idea. - do islope = 1,nslope - do ig = 1,ngrid - if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then +do islope = 1,nslope + do ig = 1,ngrid + if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia) stretch = (regolith_inertia(ig,islope)/wice_inertia)**2 - ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - & - previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) - ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch - endif - enddo - enddo - - RETURN - END - -!!! -------------------------------------- - - - SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed -!!! + previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) + ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch + endif + enddo +enddo + +END SUBROUTINE + +!---------------------------------------- +SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed +!!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use comsoil_h_PEM, only: mlayer_PEM - use comslope_mod, only: subslope_dist,def_slope_mean - use constants_marspem_mod, only: porosity - use vdifc_mod, only: compute_Tice +use comsoil_h_PEM, only: mlayer_PEM +use comslope_mod, only: subslope_dist, def_slope_mean +use constants_marspem_mod, only: porosity +use vdifc_mod, only: compute_Tice #ifndef CPP_STD use comcstfi_h, only: pi @@ -155,145 +154,137 @@ #endif - implicit none -! inputs -! ----------- - 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) :: former_ice_table_thickness(ngrid,nslope) ! ice table thickness at the former iteration [m] - real,intent(in) :: new_ice_table_thickness(ngrid,nslope) ! ice table thickness at the current iteration [m] - real,intent(in) :: ice_table_depth(ngrid,nslope) ! ice table depth [m] - real,intent(in) :: tsurf(ngrid,nslope) ! Surface temperature [K] - real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope) ! Soil temperature [K] -! outputs -! ----------- - real,intent(out) :: delta_m_h2o(ngrid) ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2] - -! local - real :: rho(ngrid,nslope) ! density of water ice [kg/m^3] - integer :: ig,islope,ilay,iref ! loop index - real :: Tice(ngrid,nslope) ! ice temperature [k] +implicit none + +! Inputs +! -------- +integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM +real, dimension(ngrid,nslope), intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m] +real, dimension(ngrid,nslope), intent(in) :: new_ice_table_thickness ! ice table thickness at the current iteration [m] +real, dimension(ngrid,nslope), intent(in) :: ice_table_depth ! ice table depth [m] +real, dimension(ngrid,nslope), intent(in) :: tsurf ! Surface temperature [K] +real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil ! Soil temperature [K] +! Outputs +! --------- +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] +! Locals +!--------- +integer :: ig, islope, ilay, iref ! loop index +real, dimension(ngrid,nslope) :: rho ! density of water ice [kg/m^3] +real, dimension(ngrid,nslope) :: Tice ! ice temperature [k] ! Code -! ----------- - rho(:,:) = 0. - Tice(:,:) = 0. +! ------ +rho = 0. +Tice = 0. !1. First let's compute Tice using a linear interpolation between the mlayer level - do ig = 1,ngrid - do islope = 1,nslope - call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope), ice_table_depth(ig,islope), Tice(ig,islope)) - rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 - enddo - enddo - - +do ig = 1,ngrid + do islope = 1,nslope + call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope)) + rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 + enddo +enddo + !2. Let's compute the amount of ice that has sublimated in each subslope - do ig = 1,ngrid - do islope = 1,nslope - delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses - *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) - enddo - enddo - - return -end subroutine - -!!! -------------------------------------- - - SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho) -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!!! -!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM -!!! +do ig = 1,ngrid + do islope = 1,nslope + delta_m_h2o(ig) = delta_m_h2o(ig) + porosity*rho(ig,islope)*(new_ice_table_thickness(ig,islope) - former_ice_table_thickness(ig,islope)) & ! convention > 0. <=> it condenses + *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) + enddo +enddo + +END SUBROUTINE + +!---------------------------------------- + +SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM +!!! !!! Author: LL !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM -use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple -implicit none -! inputs +use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM +use math_mod, only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple + +implicit none + +! Inputs ! ------ - -real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice -real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa] -real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] -real,intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] -real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] -real,intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) -integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] +real, dimension(nsoilmx_PEM), intent(in) :: porefill ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice +real, dimension(nsoilmx_PEM), intent(in) :: psat_soil ! Soil water pressure at saturation, yearly averaged [Pa] +real, intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] +real, intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] +real, intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] +real, intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) +integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] real, intent(inout) :: depth_filling ! depth where pore filling begins [m] - -! outputs +! Outputs ! ------- -integer,intent(out) :: index_filling ! index where the pore filling begins [1] -integer, intent(out) :: index_geothermal ! index where the ice table stops [1] -real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] -real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase - -! local - -real :: eta(nsoilmx_PEM) ! constriction -integer :: ilay ! index for loop -real :: old_depth_filling ! depth_filling saved -real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn -integer :: index_tmp ! for loop -real :: Jdry ! flux trought the dry layer -real :: Jsat ! flux trought the ice layer -real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer -integer :: index_firstice ! first index where ice appears (i.e., f > 0) -real :: dz_eta(nsoilmx_PEM) ! \partial z \eta -real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat -real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal - -! constant -real :: pvap2rho = 18.e-3/8.314 -! - +integer, intent(out) :: index_filling ! index where the pore filling begins [1] +integer, intent(out) :: index_geothermal ! index where the ice table stops [1] +real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] +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 +! Locals +!--------- +real, dimension(nsoilmx_PEM) :: eta ! constriction +real, dimension(nsoilmx_PEM) :: dz_psat ! first derivative of the vapor pressure at saturation +real, dimension(nsoilmx_PEM) :: dz_eta ! \partial z \eta +real, dimension(nsoilmx_PEM) :: dzz_psat ! \partial \partial psat +integer :: ilay, index_tmp ! index for loop +real :: old_depth_filling ! depth_filling saved +real :: Jdry ! flux trought the dry layer +real :: Jsat ! flux trought the ice layer +real :: Jdry_prevlay, Jsat_prevlay ! same but for the previous ice layer +integer :: index_firstice ! first index where ice appears (i.e., f > 0) +real :: massfillabove, massfillafter ! h2O mass above and after index_geothermal +! Constant +!----------- +real, parameter :: pvap2rho = 18.e-3/8.314 +! Code +!------- ! 0. Compute constriction over the layer ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 if (index_IS.lt.0) then - index_tmp = nsoilmx_PEM - do ilay = 1,nsoilmx_PEM - call constriction(porefill(ilay),eta(ilay)) - enddo + index_tmp = nsoilmx_PEM + do ilay = 1,nsoilmx_PEM + call constriction(porefill(ilay),eta(ilay)) + enddo else - index_tmp = index_IS - do ilay = 1,index_IS-1 - call constriction(porefill(ilay),eta(ilay)) - enddo - do ilay = index_IS,nsoilmx_PEM - eta(ilay) = 0. - enddo + index_tmp = index_IS + do ilay = 1,index_IS-1 + call constriction(porefill(ilay),eta(ilay)) + enddo + do ilay = index_IS,nsoilmx_PEM + eta(ilay) = 0. + enddo endif -! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) - +! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) old_depth_filling = depth_filling -call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) +call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) do ilay = 1,index_tmp - Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) - Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) - if((Jdry - Jsat).le.0) then - index_filling = ilay - exit - endif -enddo - -if(index_filling.eq.1) depth_filling = mlayer_PEM(1) - -if(index_filling.gt.1) then - Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1) - Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1) - call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling) + Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) + Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) + if (Jdry - Jsat <= 0) then + index_filling = ilay + exit + endif +enddo + +if (index_filling == 1) then + depth_filling = mlayer_PEM(1) +else if (index_filling > 1) then + Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1) + Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1) + call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling) endif - ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010)) - - ! 2.0 preliminary: depth to shallowest ice (discontinuity at interface) - index_firstice = -1 do ilay = 1,nsoilmx_PEM - if (porefill(ilay).le.0.) then + if (porefill(ilay) <= 0.) then index_firstice = ilay ! first point with ice exit @@ -302,44 +293,39 @@ ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) - -call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta) - -if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then - call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) +call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta) + +if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) + +call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) +dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat) + +! 3. Ice table boundary due to geothermal heating +if (index_IS > 0) index_geothermal = -1 +if (index_geothermal < 0) depth_geothermal = -1. +if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010) + index_geothermal = -1 + do ilay = 2,nsoilmx_PEM + if (dz_psat(ilay) > 0.) then ! first point with reversed flux + index_geothermal = ilay + call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) + exit + endif + enddo +else + index_geothermal = -1 endif - -call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) -dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:)) - -! 3. Ice table boundary due to geothermal heating - -if(index_IS.gt.0) index_geothermal = -1 -if(index_geothermal.lt.0) depth_geothermal = -1. -if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010) - index_geothermal = -1 - do ilay=2,nsoilmx_PEM - if (dz_psat(ilay).gt.0.) then ! first point with reversed flux - index_geothermal=ilay - call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) - exit - endif - enddo - else - index_geothermal = -1 - endif - if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 24 from Schorgofer, Icarus (2010) - call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) - index_tmp = -1 - do ilay=index_geothermal,nsoilmx_PEM - if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full - call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) - if (massfillafterindex_geothermal) then -! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal - call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & - dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) -! if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt. 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010) + call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) + index_tmp = -1 + do ilay = index_geothermal,nsoilmx_PEM + if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full + call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) + if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG + if (ilay > index_geothermal) then +! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal + call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, & + dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) +! if (depth_geothermal > mlayer_PEM(ilay) .or. depth_geothermal < mlayer_PEM(ilay - 1)) write(*,*) '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay - 1),depth_geothermal,mlayer_PEM(ilay) + index_tmp = ilay endif ! otherwise leave depth_geothermal unchanged @@ -347,34 +333,35 @@ endif massfillabove = massfillafter - enddo - if (index_tmp.gt.0) index_geothermal = index_tmp - end if - return -end subroutine - -!!! -------------------------------------- - + enddo + if (index_tmp > 0) index_geothermal = index_tmp +end if + +END SUBROUTINE + +!---------------------------------------- SUBROUTINE constriction(porefill,eta) - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!! - !!! Purpose: Compute the constriction of vapor flux by pore ice - !!! - !!! Author: LL - !!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - implicit none - real,intent(in) :: porefill ! pore filling fraction - real,intent(out) :: eta ! constriction - - !!! - - if (porefill.le.0.) eta = 1. - if ((porefill.gt.0.) .and.(porefill.lt.1.)) then - eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 - endif - if (porefill.le.1.) eta = 0. - return - end subroutine - -end module +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! +!!! Purpose: Compute the constriction of vapor flux by pore ice +!!! +!!! Author: LL +!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +implicit none + +real, intent(in) :: porefill ! pore filling fraction +real, intent(out) :: eta ! constriction + +if (porefill <= 0.) then + eta = 1. +else if (0 < porefill .and. porefill < 1.) then + eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 +else + eta = 0. +endif + +END SUBROUTINE + +END MODULE Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3319) +++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3320) @@ -30,10 +30,10 @@ ! Stratum = node of the linked list type :: stratum - real :: thickness ! Layer thickness [m] - real :: top_elevation ! Layer top_elevation (top height from the surface) [m] - real :: co2ice_volfrac ! CO2 ice volumetric fraction - real :: h2oice_volfrac ! H2O ice volumetric fraction - real :: dust_volfrac ! Dust volumetric fraction - real :: air_volfrac ! Air volumetric fraction inside pores + real :: thickness ! Layer thickness [m] + real :: top_elevation ! Layer top_elevation (top height from the surface) [m] + real :: co2ice_volfrac ! CO2 ice volumetric fraction + real :: h2oice_volfrac ! H2O ice volumetric fraction + real :: dust_volfrac ! Dust volumetric fraction + real :: air_volfrac ! Air volumetric fraction inside pores type(stratum), pointer :: up => null() ! Upper stratum (next node) type(stratum), pointer :: down => null() ! Lower stratum (previous node)