Index: /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3562)
+++ /trunk/LMDZ.COMMON/libf/evolution/NS_dyn_ss_ice_m.F90	(revision 3563)
@@ -37,5 +37,5 @@
   real(8)  ssi_depth_in, ssi_depth, T1
   real(8), dimension(NP) :: zdepthF, zdepthE, zdepthT, zdepthG
-  real(8), dimension(NMAX,NP) :: porefill, porefill_in
+  real(8), dimension(nz,NP) :: porefill, porefill_in
   real(8), dimension(nz) :: Tb
   real(8), dimension(NP) :: Tmean1, Tmean3, avrho1
@@ -152,5 +152,5 @@
   !print *,'Zt0=  ',ZdepthT
      call icelayer_mars(timestep,nz,NP,thIn,rhoc,z,porosity,pfrost,Tb,zdepthF, &
-          & zdepthE,porefill(1:nz,:),Tmean1,Tmean3,zdepthG, &
+          & zdepthE,porefill,Tmean1,Tmean3,zdepthG, &
           & albedo,p0,icefrac,zdepthT,avrho1, &
           & avrho1prescribed)
Index: /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3562)
+++ /trunk/LMDZ.COMMON/libf/evolution/NS_fast_subs_mars.F90	(revision 3563)
@@ -49,4 +49,6 @@
   real(8), SAVE :: avdrho_old(100), zdepth_old(100)  ! NP<=100
   logical mode2
+
+  avdrho_old = 1. ! initialization
 
   !$omp parallel &
Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3562)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3563)
@@ -525,2 +525,5 @@
 == 17/12/2024 == JBC
 As intended, years computed by the PEM runs are now the only ones to be counted for the duration of the PEM simulation. The possibility to count in addition years computed by the PCM runs is left as an option of "launchPEM.sh" with the variable 'counting' (0 = "only PEM runs count"; any other values = "PCM runs are taken into account") + several small corrections/improvements in the launching scripts.
+
+== 17/12/2024 == JBC
+Correction of Norbert Schorghofer's code due to missing initialization and bad shape array as subroutine argument + some cleanings.
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3562)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3563)
@@ -264,19 +264,4 @@
                 write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
                 write(*,*)'will reconstruct the values of Tsoil'
-!                do ig = 1,ngrid
-!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
-!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
-!                    do iloop=index_breccia+2,index_bedrock
-!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
-!                    enddo
-!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
-!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
-!
-!                    do iloop=index_bedrock+2,nsoil_PEM
-!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
-!                    enddo
-!                enddo
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -452,20 +437,4 @@
 !b) Soil temperature
         do islope = 1,nslope
-!            do ig = 1,ngrid
-!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
-!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
-!
-!                do iloop=index_breccia+2,index_bedrock
-!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
-!                enddo
-!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
-!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
-!
-!                do iloop=index_bedrock+2,nsoil_PEM
-!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
-!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
-!                enddo
-!            enddo
             call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))