Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3609) +++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt (revision 3610) @@ -572,3 +572,7 @@ == 05/02/2025 == JBC - Simplification of "read_data_PCM_mod.F90" to treat the cases with/without slope. -- Few bug corrections related to 1D. +- Few bug corrections related to 1D. + +== 05/02/2025 == JBC +- The sublimation flux can now be modified by the new function 'recomp_tend_h2o' to account for the growth of a dust lag layer (see Eran Vos's note for the formula). +- Addition of a function 'itp_tsoil' to compute the soil temperature at any point in the soil profile. Index: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90 (revision 3609) +++ trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90 (revision 3610) @@ -268,19 +268,6 @@ ! Position in the old discretization of the depth z = mlayer_PEM(isoil - 1) - zshift_surfloc - ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs - iz = 0 - do while (mlayer_PEM(iz) <= z) - iz = iz + 1 - enddo - if (iz == 0) then - tsoil_minus = tsurf(ig,islope) - mlayer_minus = 0. - else - tsoil_minus = tsoil_old(ig,iz,islope) - mlayer_minus = mlayer_PEM(iz - 1) - endif ! Interpolation of the temperature profile from the old discretization - a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) - tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus + tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) enddo endif @@ -300,19 +287,6 @@ ! Position of the lag bottom in the old discretization of the depth z = zlag(ig,islope) - zshift_surfloc - ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs - iz = 0 - do while (mlayer_PEM(iz) <= z) - iz = iz + 1 - enddo - if (iz == 0) then - tsoil_minus = tsurf(ig,islope) - mlayer_minus = 0. - else - tsoil_minus = tsoil_old(ig,iz,islope) - mlayer_minus = mlayer_PEM(iz - 1) - endif ! The "new lag" temperature is set to the ice temperature just below - a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) - tsoil(ig,:ilag,islope) = a*(z - mlayer_minus) + tsoil_minus + tsoil(ig,:ilag,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) endif @@ -326,19 +300,6 @@ ! Position in the old discretization of the depth z = mlayer_PEM(isoil - 1) - zshift_surfloc - ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs - iz = 0 - do while (mlayer_PEM(iz) <= z) - iz = iz + 1 - enddo - if (iz == 0) then - tsoil_minus = tsurf(ig,islope) - mlayer_minus = 0. - else - tsoil_minus = tsoil_old(ig,iz,islope) - mlayer_minus = mlayer_PEM(iz - 1) - endif ! Interpolation of the temperature profile from the old discretization - a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) - tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus + tsoil(ig,isoil,islope) = itp_tsoil(tsoil_old(ig,:,islope),tsurf(ig,islope),z) enddo @@ -352,3 +313,35 @@ END SUBROUTINE shift_tsoil2surf +!======================================================================= + +FUNCTION itp_tsoil(tsoil,tsurf,z) RESULT(tsoil_z) + +use comsoil_h_PEM, only: mlayer_PEM + +implicit none + +real, dimension(:), intent(in) :: tsoil +real, intent(in) :: z, tsurf + +real :: tsoil_z, tsoil_minus, mlayer_minus, a +integer :: iz + +! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs +iz = 0 +do while (mlayer_PEM(iz) <= z) + iz = iz + 1 +enddo +if (iz == 0) then + tsoil_minus = tsurf + mlayer_minus = 0. +else + tsoil_minus = tsoil(iz) + mlayer_minus = mlayer_PEM(iz - 1) +endif +! Interpolation of the temperature profile from the old discretization +a = (tsoil(iz + 1) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus) +tsoil_z = a*(z - mlayer_minus) + tsoil_minus + +END FUNCTION itp_tsoil + END MODULE compute_soiltemp_mod Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3609) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3610) @@ -65,5 +65,5 @@ use pemetat0_mod, only: pemetat0 use read_data_PCM_mod, only: read_data_PCM -use recomp_tend_co2_mod, only: recomp_tend_co2 +use recomp_tend_mod, only: recomp_tend_co2, recomp_tend_h2o use compute_soiltemp_mod, only: compute_tsoil_pem, shift_tsoil2surf use writediagpem_mod, only: writediagpem, writediagsoilpem @@ -211,4 +211,5 @@ real, dimension(:,:,:,:), allocatable :: tsoil_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K] real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K] +real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries_old ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM at the previous time step [K] real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3] real, dimension(:,:), allocatable :: watersurf_density_avg ! Grid points x Slope: Averaged water surface density [kg/m^3] @@ -227,4 +228,5 @@ real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] +real, dimension(:,:), allocatable :: icetable_depth_old ! Old depth of the ice table ! Some variables for the PEM run @@ -857,4 +859,6 @@ ! II_d.3 Update soil temperature write(*,*)"> Updating soil temperature profile" + allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen)) + tsoil_PEM_timeseries_old = tsoil_PEM_timeseries do islope = 1,nslope tsoil_avg_old = tsoil_PEM(:,:,islope) @@ -878,5 +882,5 @@ ! II_d.4 Update the ice table - allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope)) + allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope)) if (icetable_equilibrium) then write(*,*) "> Updating ice table (equilibrium method)" @@ -887,4 +891,5 @@ write(*,*) "> Updating ice table (dynamic method)" ice_porefilling_old = ice_porefilling + icetable_depth_old = icetable_depth allocate(porefill(nsoilmx_PEM)) do ig = 1,ngrid @@ -971,5 +976,10 @@ !------------------------ call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new) - deallocate(vmr_co2_PEM_phys) + do ig = 1,ngrid + do islope = 1,nslope + call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) + enddo + enddo + deallocate(vmr_co2_PEM_phys,icetable_depth_old,tsoil_PEM_timeseries_old) !------------------------ Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90 (revision 3609) +++ (revision ) @@ -1,97 +1,0 @@ -MODULE recomp_tend_co2_mod - -implicit none - -!======================================================================= -contains -!======================================================================= - -SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, & - vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global) - -use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol - -implicit none - -!======================================================================= -! -! To compute the evolution of the tendencie for co2 ice -! -!======================================================================= - -! Inputs -! ------ -integer, intent(in) :: timelen, ngrid, nslope -real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM ! physical point field: Volume mixing ratio of co2 in the first layer -real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! physical point field: Volume mixing ratio of co2 in the first layer -real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! physical point field: Surface pressure in the PCM -real, intent(in) :: ps_avg_global_ini ! global averaged pressure at previous timestep -real, intent(in) :: ps_avg_global ! global averaged pressure at current timestep -real, dimension(ngrid,nslope), intent(in) :: d_co2ice_ini ! physical point field: Evolution of perennial ice over one year -real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice per mesh and sub-grid slope (kg/m^2) -real, dimension(ngrid,nslope), intent(in) :: emissivity ! Emissivity per mesh and sub-grid slope(1) -! Outputs -! ------- -real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year - -! Local: -! ------ -integer :: i, t, islope -real :: coef, avg - -write(*,*) "> Updating the CO2 ice tendency for the new pressure" - -! Evolution of the water ice for each physical point -do i = 1,ngrid - do islope = 1,nslope - coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2 - avg = 0. - if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then - do t = 1,timelen - avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 & - - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4 - enddo - if (avg < 1.e-4) avg = 0. - d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen - endif - enddo -enddo - -END SUBROUTINE recomp_tend_co2 -!======================================================================= - -SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp) - -implicit none - -!======================================================================= -! -! To compute the evolution of the tendencie for h2o ice -! -!======================================================================= - -! Inputs -! ------ -integer, intent(in) :: timelen, ngrid, nslope -real, dimension(:), intent(in) :: PCM_temp, PEM_temp -! Outputs -! ------- -real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year - -! Local: -! ------ - -write(*,*) "Update of the H2O tendency due to lag layer" - -! Flux correction due to lag layer -!~ Rz_old = h2oice_depth_old*0.0325/4.e-4 ! resistance from PCM -!~ Rz_new = h2oice_depth_new*0.0325/4.e-4 ! new resistance based on new depth -!~ R_dec = (1./Rz_old)/(1./Rz_new) ! decrease because of resistance -!~ soil_psv_old = psv(max(PCM_temp(h2oice_depth_old))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the GCM run temperature at the old ice location -!~ soil_psv_new = psv(max(PEM_temp(h2oice_depth_new))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the PEM run temperature at the new ice location -!~ hum_dec = soil_psv_old/soil_psv_new ! decrease because of lower water vapor pressure at the new depth -!~ d_h2oice = d_h2oice*R_dec*hum_dec ! decrease of flux - -END SUBROUTINE recomp_tend_h2o - -END MODULE recomp_tend_co2_mod Index: trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3610) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3610) @@ -0,0 +1,114 @@ +MODULE recomp_tend_mod + +implicit none + +!======================================================================= +contains +!======================================================================= + +SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, & + vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global) + +use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol + +implicit none + +!======================================================================= +! +! To compute the evolution of the tendency for co2 ice +! +!======================================================================= + +! Inputs +! ------ +integer, intent(in) :: timelen, ngrid, nslope +real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM ! physical point field: Volume mixing ratio of co2 in the first layer +real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! physical point field: Volume mixing ratio of co2 in the first layer +real, dimension(ngrid,timelen), intent(in) :: ps_PCM ! physical point field: Surface pressure in the PCM +real, intent(in) :: ps_avg_global_ini ! global averaged pressure at previous timestep +real, intent(in) :: ps_avg_global ! global averaged pressure at current timestep +real, dimension(ngrid,nslope), intent(in) :: d_co2ice_ini ! physical point field: Evolution of perennial ice over one year +real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice per mesh and sub-grid slope (kg/m^2) +real, dimension(ngrid,nslope), intent(in) :: emissivity ! Emissivity per mesh and sub-grid slope(1) +! Outputs +! ------- +real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year + +! Local: +! ------ +integer :: i, t, islope +real :: coef, avg + +write(*,*) "> Updating the CO2 ice tendency for the new pressure" + +! Evolution of the water ice for each physical point +do i = 1,ngrid + do islope = 1,nslope + coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2 + avg = 0. + if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then + do t = 1,timelen + avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 & + - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4 + enddo + if (avg < 1.e-4) avg = 0. + d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen + endif + enddo +enddo + +END SUBROUTINE recomp_tend_co2 +!======================================================================= + +SUBROUTINE recomp_tend_h2o(icetable_depth_old,icetable_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice) + +use compute_soiltemp_mod, only: itp_tsoil +use fast_subs_mars, only: psv + +implicit none + +!======================================================================= +! +! To compute the evolution of the tendency for h2o ice +! +!======================================================================= + +! Inputs +! ------ +real, intent(in) :: icetable_depth_old, icetable_depth_new, tsurf +real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old, tsoil_PEM_timeseries_new +! Outputs +! ------- +real, intent(inout) :: d_h2oice ! Evolution of perennial ice over one year + +! Local: +! ------ +real :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new +integer :: t +real, parameter :: coef_diff = 4.e-4 ! Diffusion coefficient +real, parameter :: zcdv = 0.0325 ! Drag coefficient + +write(*,*) "> Updating the H2O tendency due to lag layer" + +! 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 +R_dec = Rz_new/Rz_old ! Decrease because of resistance + +! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location +psv_max_old = 0. +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))) +enddo + +! Lower humidity due to growing lag layer (higher depth) +hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth + +! Flux correction (decrease) +d_h2oice = d_h2oice*R_dec*hum_dec + +END SUBROUTINE recomp_tend_h2o + +END MODULE recomp_tend_mod