Changeset 3610 for trunk

Timestamp:

Feb 5, 2025, 3:06:50 PM (6 hours ago)

Author:

jbclement

Message:

PEM:

• 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.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• compute_soiltemp_mod.F90 (modified) (4 diffs)
• pem.F90 (modified) (7 diffs)
• recomp_tend_mod.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_mod.F90) (5 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -572 +572 @@
 == 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.

trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90

@@ -268 +268 @@
                     ! 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 +287 @@
                     ! 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 +300 @@
                     ! 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 +313 @@
 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

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -65 +65 @@
 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 +211 @@
 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 +228 @@
 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 +859 @@
 ! 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 +882 @@
 
 ! 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 +891 @@
             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 +976 @@
 !------------------------
     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)
 
 !------------------------

trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90

@@ -1 @@
-MODULE recomp_tend_co2_mod
+MODULE recomp_tend_mod
 
 implicit none
@@ -16 +16 @@
 !=======================================================================
 !
-!  To compute the evolution of the tendencie for co2 ice
+!  To compute the evolution of the tendency for co2 ice
 !
 !=======================================================================
@@ -61 +61 @@
 !=======================================================================
 
-SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp)
+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
@@ -67 +70 @@
 !=======================================================================
 !
-!  To compute the evolution of the tendencie for h2o ice
+!  To compute the evolution of the tendency for h2o ice
 !
 !=======================================================================
@@ -73 +76 @@
 ! Inputs
 ! ------
-integer,            intent(in) :: timelen, ngrid, nslope
-real, dimension(:), intent(in) :: PCM_temp, PEM_temp
+real,                 intent(in) :: icetable_depth_old, icetable_depth_new, tsurf
+real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old, tsoil_PEM_timeseries_new
 ! Outputs
 ! -------
-real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year
+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(*,*) "Update of the H2O tendency due to lag layer"
+write(*,*) "> Updating 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
+! 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_co2_mod
+END MODULE recomp_tend_mod