Changeset 3320 for trunk/LMDZ.COMMON/libf/evolution
- Timestamp:
- Apr 30, 2024, 12:08:09 PM (8 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 4 edited
-
changelog.txt (modified) (1 diff)
-
deftank/run_PEM.def (modified) (1 diff)
-
ice_table_mod.F90 (modified) (6 diffs)
-
layering_mod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3317 r3320 286 286 == 25/04/2024 == JBC 287 287 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". 288 289 == 26/04/2024 == JBC 290 - Small corrections to make the PEM work in 3D. 291 - Addition of the flag "layering_algo" (Default = .false.) defined in the "run_PEM.def" to choose to use the layering algorithm. 292 293 == 30/04/2024 == JBC 294 Small correction of a bug in the reading of "run_PEM.def" + some cleanings. -
trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def
r3319 r3320 92 92 # co2ice_flow=.true. 93 93 94 !#---------- Layering parameters ----------#94 #---------- Layering parameters ----------# 95 95 # Do you want to run with the layering algorithm? Default = .false. 96 96 # layering=.true. -
trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
r3264 r3320 1 moduleice_table_mod2 3 1 MODULE ice_table_mod 2 3 implicit none 4 4 5 5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! … … 10 10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 11 11 12 LOGICAL,SAVE :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium 13 LOGICAL,SAVE :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) 14 real,save,allocatable :: porefillingice_depth(:,:) ! ngrid x nslope: Depth of the ice table [m] 15 real,save,allocatable :: porefillingice_thickness(:,:) ! ngrid x nslope: Thickness of the ice table [m] 16 12 logical, save :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium 13 logical, save :: icetable_dynamic ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method) 14 real, allocatable, dimension(:,:), save :: porefillingice_depth ! ngrid x nslope: Depth of the ice table [m] 15 real, allocatable, dimension(:,:), save :: porefillingice_thickness ! ngrid x nslope: Thickness of the ice table [m] 16 17 !---------------------------------------- 17 18 contains 18 19 20 subroutine ini_ice_table_porefilling(ngrid,nslope) 21 22 implicit none 23 integer,intent(in) :: ngrid ! number of atmospheric columns 24 integer,intent(in) :: nslope ! number of slope within a mesh 25 26 allocate(porefillingice_depth(ngrid,nslope)) 27 allocate(porefillingice_thickness(ngrid,nslope)) 28 29 end subroutine ini_ice_table_porefilling 30 31 subroutine end_ice_table_porefilling 32 33 implicit none 34 if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) 35 if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) 36 37 end subroutine end_ice_table_porefilling 38 39 !!! -------------------------------------- 40 41 SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) 19 !---------------------------------------- 20 SUBROUTINE ini_ice_table_porefilling(ngrid,nslope) 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 27 allocate(porefillingice_depth(ngrid,nslope)) 28 allocate(porefillingice_thickness(ngrid,nslope)) 29 30 END SUBROUTINE ini_ice_table_porefilling 31 32 !---------------------------------------- 33 SUBROUTINE end_ice_table_porefilling 34 35 implicit none 36 37 if (allocated(porefillingice_depth)) deallocate(porefillingice_depth) 38 if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness) 39 40 END SUBROUTINE end_ice_table_porefilling 41 42 !---------------------------------------- 43 SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness) 42 44 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 43 45 !!! … … 45 47 !!! density at the surface and at depth. 46 48 !!! Computations are made following the methods in Schorgofer et al., 2005 47 !!! This subroutineonly gives the ice table at equilibrium and does not consider exchange with the atmosphere48 !!! 49 !!! This SUBROUTINE only gives the ice table at equilibrium and does not consider exchange with the atmosphere 50 !!! 49 51 !!! Author: LL 50 52 !!! 51 53 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 52 use math_mod,only: findroot 53 USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM ! Depth of the vertical grid 54 USE soil_thermalproperties_mod, only: ice_thermal_properties 55 implicit none 56 ! inputs 57 ! ----------- 58 integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 59 logical,intent(in) :: watercaptag(ngrid) ! Boolean to check the presence of a perennial glacier 60 real,intent(in) :: rhowatersurf_ave(ngrid,nslope) ! Water density at the surface, yearly averaged [kg/m^3] 61 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] 62 real,intent(in) :: regolith_inertia(ngrid,nslope) ! Thermal inertia of the regolith layer [SI] 54 use math_mod, only: findroot 55 use comsoil_h_PEM, only: mlayer_PEM, layer_PEM ! Depth of the vertical grid 56 use soil_thermalproperties_mod, only: ice_thermal_properties 57 58 implicit none 59 60 ! Inputs 61 ! -------- 62 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 63 logical, dimension(ngrid), intent(in) :: watercaptag ! Boolean to check the presence of a perennial glacier 64 real, dimension(ngrid,nslope), intent(in) :: rhowatersurf_ave ! Water density at the surface, yearly averaged [kg/m^3] 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 real, dimension(ngrid,nslope), intent(in) :: regolith_inertia ! Thermal inertia of the regolith layer [SI] 63 67 ! Ouputs 64 ! -------- ---65 real,intent(out) :: ice_table_beg(ngrid,nslope)! ice table depth [m]66 real,intent(out) :: ice_table_thickness(ngrid,nslope)! ice table thickness [m]67 ! Local 68 ! -------- ---69 70 71 72 73 74 real :: wice_inertia ! Water Ice thermal Inertia [USI] 68 ! -------- 69 real, dimension(ngrid,nslope), intent(out) :: ice_table_beg ! ice table depth [m] 70 real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m] 71 ! Locals 72 ! -------- 73 integer ig, islope,isoil,isoilend ! loop variables 74 real :: diff_rho(nsoil_PEM) ! difference of water vapor density between the surface and at depth [kg/m^3] 75 real :: ice_table_end ! depth of the end of the ice table [m] 76 real :: previous_icetable_depth(ngrid,nslope) ! Ice table computed at previous ice depth [m] 77 real :: stretch ! stretch factor to improve the convergence of the ice table 78 real :: wice_inertia ! Water Ice thermal Inertia [USI] 75 79 ! Code 76 ! ----------- 77 80 ! ------ 78 81 previous_icetable_depth(:,:) = ice_table_beg(:,:) 79 80 if(watercaptag(ig)) then81 82 83 82 do ig = 1,ngrid 83 if (watercaptag(ig)) then 84 ice_table_beg(ig,:) = 0. 85 ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) 86 else 84 87 do islope = 1,nslope 85 ice_table_beg(ig,islope) = -1.86 ice_table_thickness(ig,islope) = 0.87 do isoil = 1,nsoil_PEM88 diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)89 enddo90 if(diff_rho(1) > 0) then! ice is at the surface91 ice_table_beg(ig,islope) = 0.92 do isoilend = 2,nsoil_PEM-193 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then94 call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)95 ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)96 exit97 endif98 enddo99 else100 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 root101 if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then102 call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope))103 ! Now let's find the end of the ice table104 ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)105 do isoilend = isoil+1,nsoil_PEM-1106 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then107 call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)108 ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)109 exit110 endif111 enddo112 endif !diff_rho(z) <0 & diff_rho(z+1) > 0113 enddo114 endif ! diff_rho(1) > 088 ice_table_beg(ig,islope) = -1. 89 ice_table_thickness(ig,islope) = 0. 90 do isoil = 1,nsoil_PEM 91 diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) 92 enddo 93 if (diff_rho(1) > 0) then ! ice is at the surface 94 ice_table_beg(ig,islope) = 0. 95 do isoilend = 2,nsoil_PEM-1 96 if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend+1) < 0.) then 97 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) 98 ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) 99 exit 100 endif 101 enddo 102 else 103 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 104 if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then 105 call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope)) 106 ! Now let's find the end of the ice table 107 ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.) 108 do isoilend = isoil+1,nsoil_PEM-1 109 if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then 110 call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end) 111 ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope) 112 exit 113 endif 114 enddo 115 endif !diff_rho(z) <0 & diff_rho(z+1) > 0 116 enddo 117 endif ! diff_rho(1) > 0 115 118 enddo 116 117 119 endif ! watercaptag 120 enddo 118 121 119 122 ! Small trick to speed up the convergence, Oded's idea. 120 121 122 if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then123 do islope = 1,nslope 124 do ig = 1,ngrid 125 if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then 123 126 call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia) 124 127 stretch = (regolith_inertia(ig,islope)/wice_inertia)**2 125 126 128 ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - & 127 previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) 128 ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch 129 endif 130 enddo 131 enddo 132 133 RETURN 134 END 135 136 !!! -------------------------------------- 137 138 139 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) 140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 141 !!! 142 !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed 143 !!! 129 previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch) 130 ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch 131 endif 132 enddo 133 enddo 134 135 END SUBROUTINE 136 137 !---------------------------------------- 138 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) 139 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 140 !!! 141 !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed 142 !!! 144 143 !!! Author: LL 145 144 !!! 146 145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 147 use comsoil_h_PEM,only: mlayer_PEM148 use comslope_mod, only: subslope_dist,def_slope_mean149 150 use vdifc_mod,only: compute_Tice146 use comsoil_h_PEM, only: mlayer_PEM 147 use comslope_mod, only: subslope_dist, def_slope_mean 148 use constants_marspem_mod, only: porosity 149 use vdifc_mod, only: compute_Tice 151 150 #ifndef CPP_STD 152 151 use comcstfi_h, only: pi … … 155 154 #endif 156 155 157 implicit none 158 ! inputs 159 ! ----------- 160 integer,intent(in) :: ngrid,nslope,nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 161 real,intent(in) :: former_ice_table_thickness(ngrid,nslope) ! ice table thickness at the former iteration [m] 162 real,intent(in) :: new_ice_table_thickness(ngrid,nslope) ! ice table thickness at the current iteration [m] 163 real,intent(in) :: ice_table_depth(ngrid,nslope) ! ice table depth [m] 164 real,intent(in) :: tsurf(ngrid,nslope) ! Surface temperature [K] 165 real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope) ! Soil temperature [K] 166 ! outputs 167 ! ----------- 168 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] 169 170 ! local 171 real :: rho(ngrid,nslope) ! density of water ice [kg/m^3] 172 integer :: ig,islope,ilay,iref ! loop index 173 real :: Tice(ngrid,nslope) ! ice temperature [k] 156 implicit none 157 158 ! Inputs 159 ! -------- 160 integer, intent(in) :: ngrid, nslope, nsoil_PEM ! Size of the physical grid, number of subslope, number of soil layer in the PEM 161 real, dimension(ngrid,nslope), intent(in) :: former_ice_table_thickness ! ice table thickness at the former iteration [m] 162 real, dimension(ngrid,nslope), intent(in) :: new_ice_table_thickness ! ice table thickness at the current iteration [m] 163 real, dimension(ngrid,nslope), intent(in) :: ice_table_depth ! ice table depth [m] 164 real, dimension(ngrid,nslope), intent(in) :: tsurf ! Surface temperature [K] 165 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil ! Soil temperature [K] 166 ! Outputs 167 ! --------- 168 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 ! Locals 170 !--------- 171 integer :: ig, islope, ilay, iref ! loop index 172 real, dimension(ngrid,nslope) :: rho ! density of water ice [kg/m^3] 173 real, dimension(ngrid,nslope) :: Tice ! ice temperature [k] 174 174 ! Code 175 ! ------ -----176 rho(:,:)= 0.177 Tice(:,:)= 0.175 ! ------ 176 rho = 0. 177 Tice = 0. 178 178 !1. First let's compute Tice using a linear interpolation between the mlayer level 179 do ig = 1,ngrid 180 do islope = 1,nslope 181 call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope), ice_table_depth(ig,islope), Tice(ig,islope)) 182 rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+ 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 183 enddo 184 enddo 185 186 179 do ig = 1,ngrid 180 do islope = 1,nslope 181 call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope),ice_table_depth(ig,islope),Tice(ig,islope)) 182 rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2 + 0.0351* Tice(ig,islope) + 933.5030 ! Rottgers, 2012 183 enddo 184 enddo 185 187 186 !2. Let's compute the amount of ice that has sublimated in each subslope 188 do ig = 1,ngrid 189 do islope = 1,nslope 190 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 191 *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) 192 enddo 193 enddo 194 195 return 196 end subroutine 197 198 !!! -------------------------------------- 199 200 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) 201 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 202 !!! 203 !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM 204 !!! 187 do ig = 1,ngrid 188 do islope = 1,nslope 189 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 190 *subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.) 191 enddo 192 enddo 193 194 END SUBROUTINE 195 196 !---------------------------------------- 197 198 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) 199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 200 !!! 201 !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM 202 !!! 205 203 !!! Author: LL 206 204 !!! 207 205 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 208 use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM 209 use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple 210 implicit none 211 ! inputs 206 use comsoil_h_PEM, only: nsoilmx_PEM, mlayer_PEM 207 use math_mod, only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple 208 209 implicit none 210 211 ! Inputs 212 212 ! ------ 213 214 real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice 215 real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa] 216 real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] 217 real,intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] 218 real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] 219 real,intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) 220 integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] 213 real, dimension(nsoilmx_PEM), intent(in) :: porefill ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice 214 real, dimension(nsoilmx_PEM), intent(in) :: psat_soil ! Soil water pressure at saturation, yearly averaged [Pa] 215 real, intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa] 216 real, intent(in) :: pwat_surf ! Water vapor pressure at the surface, not necesseraly at saturation, yearly averaged [Pa] 217 real, intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa] 218 real, intent(in) :: B ! constant (Eq. 8 from Schorgofer, Icarus (2010).) 219 integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1] 221 220 real, intent(inout) :: depth_filling ! depth where pore filling begins [m] 222 223 ! outputs 221 ! Outputs 224 222 ! ------- 225 integer,intent(out) :: index_filling ! index where the pore filling begins [1] 226 integer, intent(out) :: index_geothermal ! index where the ice table stops [1] 227 real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] 228 real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase 229 230 ! local 231 232 real :: eta(nsoilmx_PEM) ! constriction 233 integer :: ilay ! index for loop 234 real :: old_depth_filling ! depth_filling saved 235 real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn 236 integer :: index_tmp ! for loop 237 real :: Jdry ! flux trought the dry layer 238 real :: Jsat ! flux trought the ice layer 239 real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer 240 integer :: index_firstice ! first index where ice appears (i.e., f > 0) 241 real :: dz_eta(nsoilmx_PEM) ! \partial z \eta 242 real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat 243 real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal 244 245 ! constant 246 real :: pvap2rho = 18.e-3/8.314 247 ! 248 223 integer, intent(out) :: index_filling ! index where the pore filling begins [1] 224 integer, intent(out) :: index_geothermal ! index where the ice table stops [1] 225 real, intent(out) :: depth_geothermal ! depth where the ice table stops [m] 226 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 227 ! Locals 228 !--------- 229 real, dimension(nsoilmx_PEM) :: eta ! constriction 230 real, dimension(nsoilmx_PEM) :: dz_psat ! first derivative of the vapor pressure at saturation 231 real, dimension(nsoilmx_PEM) :: dz_eta ! \partial z \eta 232 real, dimension(nsoilmx_PEM) :: dzz_psat ! \partial \partial psat 233 integer :: ilay, index_tmp ! index for loop 234 real :: old_depth_filling ! depth_filling saved 235 real :: Jdry ! flux trought the dry layer 236 real :: Jsat ! flux trought the ice layer 237 real :: Jdry_prevlay, Jsat_prevlay ! same but for the previous ice layer 238 integer :: index_firstice ! first index where ice appears (i.e., f > 0) 239 real :: massfillabove, massfillafter ! h2O mass above and after index_geothermal 240 ! Constant 241 !----------- 242 real, parameter :: pvap2rho = 18.e-3/8.314 243 ! Code 244 !------- 249 245 ! 0. Compute constriction over the layer 250 246 ! Within the ice sheet, constriction is set to 0. Elsewhere, constriction = (1-porefilling)**2 251 247 if (index_IS.lt.0) then 252 index_tmp = nsoilmx_PEM253 do ilay = 1,nsoilmx_PEM254 call constriction(porefill(ilay),eta(ilay))255 enddo248 index_tmp = nsoilmx_PEM 249 do ilay = 1,nsoilmx_PEM 250 call constriction(porefill(ilay),eta(ilay)) 251 enddo 256 252 else 257 index_tmp = index_IS258 do ilay = 1,index_IS-1259 call constriction(porefill(ilay),eta(ilay))260 enddo261 do ilay = index_IS,nsoilmx_PEM262 eta(ilay) = 0.263 enddo253 index_tmp = index_IS 254 do ilay = 1,index_IS-1 255 call constriction(porefill(ilay),eta(ilay)) 256 enddo 257 do ilay = index_IS,nsoilmx_PEM 258 eta(ilay) = 0. 259 enddo 264 260 endif 265 261 266 ! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) 267 262 ! 1. Depth at which pore filling occurs. We solve Eq. 9 from Schorgofer, Icarus (2010) 268 263 old_depth_filling = depth_filling 269 264 270 call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 265 call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat) 271 266 272 267 do ilay = 1,index_tmp 273 Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010)274 Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)275 if((Jdry - Jsat).le.0) then276 index_filling = ilay277 exit278 endif279 enddo 280 281 if (index_filling.eq.1) depth_filling = mlayer_PEM(1)282 283 if(index_filling.gt.1) then284 Jdry_prevlay = (psat_soil(index_filling -1) - pwat_surf)/mlayer_PEM(index_filling-1)285 Jsat_prevlay = eta(index_filling -1)*dz_psat(index_filling-1)286 call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)268 Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus (2010) 269 Jsat = eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010) 270 if (Jdry - Jsat <= 0) then 271 index_filling = ilay 272 exit 273 endif 274 enddo 275 276 if (index_filling == 1) then 277 depth_filling = mlayer_PEM(1) 278 else if (index_filling > 1) then 279 Jdry_prevlay = (psat_soil(index_filling - 1) - pwat_surf)/mlayer_PEM(index_filling - 1) 280 Jsat_prevlay = eta(index_filling - 1)*dz_psat(index_filling - 1) 281 call findroot(Jdry - Jsat,Jdry_prevlay - Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling - 1),depth_filling) 287 282 endif 288 283 289 290 284 ! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010)) 291 292 293 285 ! 2.0 preliminary: depth to shallowest ice (discontinuity at interface) 294 295 286 index_firstice = -1 296 287 do ilay = 1,nsoilmx_PEM 297 if (porefill(ilay).le.0.) then288 if (porefill(ilay) <= 0.) then 298 289 index_firstice = ilay ! first point with ice 299 290 exit … … 302 293 303 294 ! 2.1: now we can computeCompute d_z (eta* d_z(rho)) 304 305 call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta) 306 307 if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then 308 call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) 295 call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM - 1),dz_eta) 296 297 if (index_firstice > 0 .and. index_firstice < nsoilmx_PEM - 2) call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice)) 298 299 call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) 300 dz_etadz_rho = pvap2rho*(dz_eta*dz_psat + eta*dzz_psat) 301 302 ! 3. Ice table boundary due to geothermal heating 303 if (index_IS > 0) index_geothermal = -1 304 if (index_geothermal < 0) depth_geothermal = -1. 305 if ((index_geothermal > 0).and.(index_IS < 0)) then ! Eq. 21 from Schorfoger, Icarus (2010) 306 index_geothermal = -1 307 do ilay = 2,nsoilmx_PEM 308 if (dz_psat(ilay) > 0.) then ! first point with reversed flux 309 index_geothermal = ilay 310 call findroot(dz_psat(ilay - 1),dz_psat(ilay),mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) 311 exit 312 endif 313 enddo 314 else 315 index_geothermal = -1 309 316 endif 310 311 call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat) 312 dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:)) 313 314 ! 3. Ice table boundary due to geothermal heating 315 316 if(index_IS.gt.0) index_geothermal = -1 317 if(index_geothermal.lt.0) depth_geothermal = -1. 318 if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010) 319 index_geothermal = -1 320 do ilay=2,nsoilmx_PEM 321 if (dz_psat(ilay).gt.0.) then ! first point with reversed flux 322 index_geothermal=ilay 323 call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) 324 exit 325 endif 326 enddo 327 else 328 index_geothermal = -1 329 endif 330 if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 24 from Schorgofer, Icarus (2010) 331 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) 332 index_tmp = -1 333 do ilay=index_geothermal,nsoilmx_PEM 334 if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full 335 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) 336 if (massfillafter<dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG 337 if (ilay>index_geothermal) then 338 ! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal 339 call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, & 340 dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal) 341 ! if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*) 342 ! '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay) 343 index_tmp=ilay 317 if (index_geothermal > 0 .and. index_IS < 0) then ! Eq. 24 from Schorgofer, Icarus (2010) 318 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) 319 index_tmp = -1 320 do ilay = index_geothermal,nsoilmx_PEM 321 if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle ! eta=0 means completely full 322 call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter) 323 if (massfillafter < dz_psat(ilay)*pvap2rho*B) then ! usually executes on i=typeG 324 if (ilay > index_geothermal) then 325 ! write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal 326 call findroot(dz_psat(ilay - 1)*pvap2rho*B - massfillabove, & 327 dz_psat(ilay)*pvap2rho*B - massfillafter,mlayer_PEM(ilay - 1),mlayer_PEM(ilay),depth_geothermal) 328 ! 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) 329 index_tmp = ilay 344 330 endif 345 331 ! otherwise leave depth_geothermal unchanged … … 347 333 endif 348 334 massfillabove = massfillafter 349 enddo 350 if (index_tmp.gt.0) index_geothermal = index_tmp 351 end if 352 return 353 end subroutine 354 355 !!! -------------------------------------- 356 335 enddo 336 if (index_tmp > 0) index_geothermal = index_tmp 337 end if 338 339 END SUBROUTINE 340 341 !---------------------------------------- 357 342 SUBROUTINE constriction(porefill,eta) 358 343 359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 360 !!! 361 !!! Purpose: Compute the constriction of vapor flux by pore ice 362 !!! 363 !!! Author: LL 364 !!! 365 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 366 implicit none 367 real,intent(in) :: porefill ! pore filling fraction 368 real,intent(out) :: eta ! constriction 369 370 !!! 371 372 if (porefill.le.0.) eta = 1. 373 if ((porefill.gt.0.) .and.(porefill.lt.1.)) then 374 eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 375 endif 376 if (porefill.le.1.) eta = 0. 377 return 378 end subroutine 379 380 end module 344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 345 !!! 346 !!! Purpose: Compute the constriction of vapor flux by pore ice 347 !!! 348 !!! Author: LL 349 !!! 350 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 351 352 implicit none 353 354 real, intent(in) :: porefill ! pore filling fraction 355 real, intent(out) :: eta ! constriction 356 357 if (porefill <= 0.) then 358 eta = 1. 359 else if (0 < porefill .and. porefill < 1.) then 360 eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 361 else 362 eta = 0. 363 endif 364 365 END SUBROUTINE 366 367 END MODULE -
trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
r3319 r3320 30 30 ! Stratum = node of the linked list 31 31 type :: stratum 32 real :: thickness! Layer thickness [m]33 real :: top_elevation! Layer top_elevation (top height from the surface) [m]34 real :: co2ice_volfrac! CO2 ice volumetric fraction35 real :: h2oice_volfrac! H2O ice volumetric fraction36 real :: dust_volfrac! Dust volumetric fraction37 real :: air_volfrac! Air volumetric fraction inside pores32 real :: thickness ! Layer thickness [m] 33 real :: top_elevation ! Layer top_elevation (top height from the surface) [m] 34 real :: co2ice_volfrac ! CO2 ice volumetric fraction 35 real :: h2oice_volfrac ! H2O ice volumetric fraction 36 real :: dust_volfrac ! Dust volumetric fraction 37 real :: air_volfrac ! Air volumetric fraction inside pores 38 38 type(stratum), pointer :: up => null() ! Upper stratum (next node) 39 39 type(stratum), pointer :: down => null() ! Lower stratum (previous node)
Note: See TracChangeset
for help on using the changeset viewer.
