Changeset 3784

Timestamp:

May 30, 2025, 11:45:56 AM (4 weeks ago)

Author:

jbclement

Message:

PEM:
Correction of 'h2o_ice_depth' update in the main loop.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (5 diffs)
• recomp_tend_mod.F90 (modified) (2 diffs)

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

@@ -678 +678 @@
 - Gathering the processes of "ice condensation" and "dust sedimentation" in the same framework of the layering algorithm
 - Reordering the workflow to manage the different situations of the layering algorithm
+
+== 30/05/2025 == JBC
+Correction of 'h2o_ice_depth' update in the main loop.

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

@@ -200 +200 @@
 type(ptrarray), dimension(:,:), allocatable :: current           ! Current active stratum in the layering
 logical,        dimension(:,:), allocatable :: new_str, new_lag  ! Flags for the layering algorithm
-real                                        :: h2o_ice_depth_old ! Old depth of subsurface ice layer
+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
 
@@ -751 +751 @@
 stopPEM = 0
 if (layering_algo) then
-    allocate(new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
+    allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
     new_str = .true.
     new_lag = .true.
@@ -834 +834 @@
     allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
     if (layering_algo) then
+        h2o_ice_depth_old = h2o_ice_depth
         do islope = 1,nslope
             do ig = 1,ngrid
                 call make_layering(layerings_map(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p)
-                co2_ice = 0.
-                h2o_ice = 0.
-                h2o_ice_depth = 0.
+                co2_ice(ig,islope) = 0.
+                h2o_ice(ig,islope) = 0.
+                h2o_ice_depth(ig,islope) = 0.
                 if (is_co2ice_str(layerings_map(ig,islope)%top)) then
                     co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
@@ -1064 +1065 @@
         do ig = 1,ngrid
             do islope = 1,nslope
-                if (is_h2oice_sublim_ini(ig,islope)) then
-                    h2o_ice_depth_old = h2o_ice_depth(ig,islope)
-                    h2o_ice_depth(ig,islope) = thickness_toplag(layerings_map(ig,islope))
-                    call recomp_tend_h2o(h2o_ice_depth_old,h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope))
-                endif
+                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))
             enddo
         enddo
@@ -1127 +1124 @@
 deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
 deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
-if (layering_algo) deallocate(new_str,new_lag,current)
+if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
 !------------------------------ END RUN --------------------------------
 

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

@@ -92 +92 @@
 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
-if (abs(Rz_old) < 1.e-10) then
-    R_dec = 1.
-else
-    R_dec = Rz_new/Rz_old ! Decrease because of resistance
-endif
+R_dec = 1.
+if (Rz_new >= Rz_old) 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
@@ -107 +104 @@
 
 ! Lower humidity due to growing lag layer (higher depth)
-print*, psv_max_old,psv_max_new
 hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth
 
 ! Flux correction (decrease)
-print*, 'aaaaah', R_dec,hum_dec
 d_h2oice = d_h2oice/R_dec/hum_dec