Changeset 3784 for trunk/LMDZ.COMMON/libf/evolution/pem.F90

Timestamp:

May 30, 2025, 11:45:56 AM (9 days ago)

Author:

jbclement

Message:

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

File:

• trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (5 diffs)

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