Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3789) @@ -692,2 +692,5 @@ == 02/06/2025 == JBC Optimization (computing time is devided by 2.2) of the program "reshape_XIOS_output" to convert XIOS output onto the PCM grid. + +== 03/06/2025 == JBC +Few corrections for the layering algorithm, in particular when a dust lag layer is created. Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3789) @@ -28,6 +28,6 @@ real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] -real, parameter :: h_patchy_dust = 0.05 ! Height under which a dust lag is considered patchy [m] -real, parameter :: h_patchy_ice = 0.05 ! Height under which a H2O ice lag is considered patchy [m] +real, parameter :: h_patchy_dust = 5.e-4 ! Height under which a dust lag is considered patchy [m] +real, parameter :: h_patchy_ice = 5.e-4 ! Height under which a H2O ice lag is considered patchy [m] ! Parameters for CO2 flux correction (inspired by Levrard et al. 2007) @@ -550,5 +550,5 @@ !======================================================================= ! To get data about possible subsurface ice in a layering -SUBROUTINE subsurface_ice_layering(this,is_subsurface_ice,h_toplag,h2o_ice) +SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice) implicit none @@ -556,5 +556,4 @@ !---- Arguments type(layering), intent(in) :: this -logical, intent(out) :: is_subsurface_ice real, intent(out) :: h2o_ice, h_toplag @@ -565,5 +564,4 @@ h_toplag = 0. h2o_ice = 0. -is_subsurface_ice = .false. current => this%top ! Is the considered stratum a dust lag? @@ -574,11 +572,7 @@ enddo ! Is the stratum below the dust lag made of ice? -if (is_h2oice_str(current)) then - if (h_toplag >= h_patchy_dust) then - is_subsurface_ice = .true. - else - h_toplag = 0. - h2o_ice = current%h_h2oice - endif +if (.not. is_h2oice_str(current) .or. h_toplag < h_patchy_dust) then + h_toplag = 0. + h2o_ice = current%h_h2oice endif @@ -622,5 +616,5 @@ !---- Local variables -real :: h_ice, h2subl, h_pore, h_tot +real :: h2subl, h_ice, h_pore, h_pore_new, h_tot real :: hlag_dust, hlag_h2oice type(stratum), pointer :: current @@ -628,4 +622,5 @@ !---- Code h_ice = str%h_co2ice + str%h_h2oice +h_pore = str%h_pore h2subl = h2subl_co2ice + h2subl_h2oice @@ -654,13 +649,13 @@ ! Which porosity is considered? if (hlag_dust >= hlag_h2oice) then - h_pore = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity) + h_pore_new = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity) else - h_pore = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity) -endif -h_tot = hlag_dust + hlag_h2oice + h_pore + h_pore_new = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity) +endif +h_tot = hlag_dust + hlag_h2oice + h_pore_new ! New stratum above the current stratum as a lag layer if (new_lag) then - call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore,0.) + call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore_new,0.) new_lag = .false. else @@ -668,5 +663,5 @@ str%up%h_h2oice = str%up%h_h2oice + hlag_h2oice str%up%h_dust = str%up%h_dust + hlag_dust - str%up%h_pore = str%up%h_pore + h_pore + str%up%h_pore = str%up%h_pore + h_pore_new endif @@ -702,5 +697,10 @@ dh_co2ice = d_co2ice/rho_co2ice dh_h2oice = d_h2oice/rho_h2oice -dh_dust = d_dust/rho_dust +! To disable dust sedimenation when there is ice sublimation (because dust tendency is fixed) +if (dh_co2ice > 0. .or. dh_h2oice > 0.) then + dh_dust = d_dust/rho_dust +else + dh_dust = 0. +endif zshift_surf = this%top%top_elevation zlag = 0. @@ -742,5 +742,5 @@ h2subl_h2oice_tot = abs(dh_h2oice) h2lift_tot = abs(dh_dust) - do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. h2lift_tot > 0. .and. associated(current)) + do while ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. associated(current)) h_co2ice_old = current%h_co2ice h_h2oice_old = current%h_h2oice @@ -755,7 +755,7 @@ currentloc => currentloc%up enddo - if (currentloc%h_h2oice > h_patchy_ice) then + if (currentloc%h_h2oice >= h_patchy_ice) then R_subl = 0. - else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice <= h_patchy_ice) then + else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then h_toplag = h_toplag + thickness(currentloc) R_subl = fred_subl**(h_toplag/hmin_lag) @@ -828,5 +828,5 @@ write(*,'(a,es12.3)') ' > dh_co2ice [m/y] = ', dh_co2ice write(*,'(a,es12.3)') ' > dh_h2oice [m/y] = ', dh_h2oice - write(*,'(a,es12.3)') ' > dh_dust [m/y] = ', dh_dust, tol + write(*,'(a,es12.3)') ' > dh_dust [m/y] = ', dh_dust error stop endif Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3789) @@ -69,6 +69,7 @@ use writediagpem_mod, only: writediagpem, writediagsoilpem use co2condens_mod, only: CO2cond_ps -use layering_mod, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, & - ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str +use layering_mod, only: layering, del_layering, make_layering, layering_algo, subsurface_ice_layering, & + ptrarray, stratum, get_nb_str_max, nb_str_max, is_dust_lag, is_co2ice_str, is_h2oice_str, & + print_layering use dyn_ss_ice_m_mod, only: dyn_ss_ice_m use version_info_mod, only: print_version_info @@ -200,5 +201,4 @@ logical, dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm real, dimension(:,:), allocatable :: h2o_ice_depth_old ! Old depth of subsurface ice layer -logical :: is_subsurface_ice ! Boolean to know if there is subsurface ice ! Variables for slopes @@ -846,8 +846,9 @@ h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice else - call subsurface_ice_layering(layerings_map(ig,islope),is_subsurface_ice,h2o_ice_depth(ig,islope),h2o_ice(ig,islope)) + call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope)) endif enddo enddo + call print_layering(layerings_map(1,1)) else zlag = 0. @@ -1065,5 +1066,5 @@ do ig = 1,ngrid do islope = 1,nslope - if (is_h2oice_sublim_ini(ig,islope)) call recomp_tend_h2o(h2o_ice_depth_old(ig,islope),h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) + if (is_h2oice_sublim_ini(ig,islope) .and. h2o_ice_depth(ig,islope) > 0.) call recomp_tend_h2o(h2o_ice_depth_old(ig,islope),h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) enddo enddo Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3789) @@ -171,7 +171,7 @@ exit endif - call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,6),found) - if (.not. found) then - write(*,*) 'Pemetat0: failed loading !' + call get_field('stratif_slope'//num//'_poreice_volfrac',stratif_array(:,islope,:,6),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' exit endif Index: trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3789) @@ -104,5 +104,5 @@ call put_field('stratif_slope'//num//'_h_dust','Layering dust height',stratif_array(:,islope,:,4),Year) call put_field('stratif_slope'//num//'_h_pore','Layering pore height',stratif_array(:,islope,:,5),Year) - call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,6),Year) + call put_field('stratif_slope'//num//'_poreice_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,6),Year) enddo deallocate(stratif_array) Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3786) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3789) @@ -61,5 +61,5 @@ !======================================================================= -SUBROUTINE recomp_tend_h2o(icetable_depth_old,icetable_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice) +SUBROUTINE recomp_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice) use compute_soiltemp_mod, only: itp_tsoil @@ -76,5 +76,5 @@ ! Inputs ! ------ -real, intent(in) :: icetable_depth_old, icetable_depth_new, tsurf +real, intent(in) :: h2oice_depth_old, h2oice_depth_new, tsurf real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old, tsoil_PEM_timeseries_new ! Outputs @@ -90,6 +90,6 @@ ! Higher resistance due to growing lag layer (higher depth) -Rz_old = icetable_depth_old*zcdv/coef_diff ! Old resistance from PCM -Rz_new = icetable_depth_new*zcdv/coef_diff ! New resistance based on new depth +Rz_old = h2oice_depth_old*zcdv/coef_diff ! Old resistance from PCM +Rz_new = h2oice_depth_new*zcdv/coef_diff ! New resistance based on new depth R_dec = 1. if (Rz_new >= Rz_old) R_dec = Rz_new/Rz_old ! Decrease because of resistance @@ -99,6 +99,6 @@ psv_max_new = 0. do t = 1,size(tsoil_PEM_timeseries_old,2) - psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,icetable_depth_old))) - psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,icetable_depth_new))) + psv_max_old = max(psv_max_old,psv(itp_tsoil(tsoil_PEM_timeseries_old(:,t),tsurf,h2oice_depth_old))) + psv_max_new = max(psv_max_new,psv(itp_tsoil(tsoil_PEM_timeseries_new(:,t),tsurf,h2oice_depth_new))) enddo