Changeset 3789 for trunk/LMDZ.COMMON/libf

Timestamp:

Jun 3, 2025, 5:19:45 PM (2 weeks ago)

Author:

jbclement

Message:

PEM:
Few corrections for the layering algorithm, in particular when a dust lag layer is created.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• layering_mod.F90 (modified) (13 diffs)
• pem.F90 (modified) (4 diffs)
• pemetat0.F90 (modified) (1 diff)
• pemredem.F90 (modified) (1 diff)
• recomp_tend_mod.F90 (modified) (4 diffs)

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

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

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

@@ -28 +28 @@
 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 +550 @@
 !=======================================================================
 ! 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 +556 @@
 !---- Arguments
 type(layering), intent(in) :: this
-logical, intent(out) :: is_subsurface_ice
 real,    intent(out) :: h2o_ice, h_toplag
 
@@ -565 +564 @@
 h_toplag = 0.
 h2o_ice = 0.
-is_subsurface_ice = .false.
 current => this%top
 ! Is the considered stratum a dust lag?
@@ -574 +572 @@
 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 +616 @@
 
 !---- 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 +622 @@
 !---- Code
 h_ice = str%h_co2ice + str%h_h2oice
+h_pore = str%h_pore
 h2subl = h2subl_co2ice + h2subl_h2oice
 
@@ -654 +649 @@
 ! 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 +663 @@
     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 +697 @@
 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 +742 @@
         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 +755 @@
                     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 +828 @@
     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

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

@@ -69 +69 @@
 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 +201 @@
 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 +846 @@
                     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 +1066 @@
         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

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

@@ -171 +171 @@
                 exit
             endif
-            call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,6),found)
-            if (.not. found) then
-                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_icepore_volfrac>!'
+            call get_field('stratif_slope'//num//'_poreice_volfrac',stratif_array(:,islope,:,6),found)
+            if (.not. found) then
+                write(*,*) 'Pemetat0: failed loading <stratif_slope'//num//'_poreice_volfrac>!'
                 exit
             endif

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

@@ -104 +104 @@
         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)

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

@@ -61 +61 @@
 !=======================================================================
 
-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 +76 @@
 ! 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 +90 @@
 
 ! 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 +99 @@
 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