Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3308)
@@ -277,2 +277,5 @@
 == 15/04/2024 == JBC
 Update of "output_layering.py" in the deftank.
+
+== 19/04/2024 == JBC
+Few small corrections to make the PEM work in 3D, in particular concerning the initialization of the planet type and the evolution of ice with slopes.
Index: trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/deftank/run_PEM.def	(revision 3308)
@@ -1,6 +1,8 @@
-#
-#-------------------------------------
-# Run control parameters for the PEM
-#-------------------------------------
+#------------------------------------#
+# Run control parameters for the PEM #
+#------------------------------------#
+
+#---------- Planet type ----------#
+planet_type=mars
 
 #---------- Output parameters ----------#
Index: trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90	(revision 3308)
@@ -28,5 +28,5 @@
 
 !   local:
-!   ----
+!   ------
 
 !=======================================================================
@@ -105,19 +105,7 @@
         enddo
     enddo
-    ! We adapt the tendencies to conserve h2o and do only exchange between grid points
-    if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
-        where (tend_h2o_ice < 0.) ! We lower the sublimating rate by a coefficient
-            new_tendencies = tend_h2o_ice*pos_tend/neg_tend
-        elsewhere ! We don't change the condensing rate
-            new_tendencies = tend_h2o_ice
-        end where
-    else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
-        where (tend_h2o_ice < 0.) ! We don't change the sublimating rate
-            new_tendencies = tend_h2o_ice
-        elsewhere ! We lower the condensing rate by a coefficient
-            new_tendencies = tend_h2o_ice*neg_tend/pos_tend
-        end where
-    else if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) == 1.e-10) then
-        write(*,*) "Reason of stopping: there is no sublimating or condensing h20 ice!"
+
+    if (abs(pos_tend) < 1.e-10 .or. abs(neg_tend) < 1.e-10) then
+        write(*,*) "Reason of stopping: there is no sublimating or condensing h2o ice!"
         write(*,*) "Tendencies on ice sublimating =", neg_tend
         write(*,*) "Tendencies on ice increasing =", pos_tend
@@ -125,4 +113,19 @@
         stopPEM = 2
         new_tendencies = 0.
+    else
+        ! We adapt the tendencies to conserve h2o and do only exchange between grid points
+        if (neg_tend > pos_tend .and. pos_tend > 0.) then ! More sublimation on the planet than condensation
+            where (tend_h2o_ice < 0.) ! We lower the sublimating rate by a coefficient
+                new_tendencies = tend_h2o_ice*pos_tend/neg_tend
+            elsewhere ! We don't change the condensing rate
+                new_tendencies = tend_h2o_ice
+            end where
+        else if (neg_tend < pos_tend .and. neg_tend > 0.) then ! More condensation on the planet than sublimation
+            where (tend_h2o_ice < 0.) ! We don't change the sublimating rate
+                new_tendencies = tend_h2o_ice
+            elsewhere ! We lower the condensing rate by a coefficient
+                new_tendencies = tend_h2o_ice*neg_tend/pos_tend
+            end where
+        endif
     endif
 
@@ -150,5 +153,9 @@
     endif
     ! In the place of accumulation of ice, we remove a bit of ice in order to conserve h2o
-    where (new_tendencies > 0.) h2o_ice = h2o_ice - new_tendencies*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
+    do islope = 1,nslope
+        do i = 1,ngrid
+             if (new_tendencies(i,islope) > 0.) h2o_ice(i,islope) = h2o_ice(i,islope) - new_tendencies(i,islope)*real_coefficient*dt_pem*cos(def_slope_mean(islope)*pi/180.)
+        enddo
+     enddo
 else ! ngrid == 1, i.e. in 1D
     h2o_ice = h2o_ice + tend_h2o_ice*dt_pem
Index: trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90	(revision 3308)
@@ -5,6 +5,6 @@
 
 ! Flags for ice management
-logical :: h2oice_flow ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
-logical :: co2ice_flow ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
+logical :: h2oice_flow  ! True by default, to compute H2O ice flow. Read in "run_PEM.def"
+logical :: co2ice_flow  ! True by default, to compute CO2 ice flow. Read in "run_PEM.def"
 logical :: metam_h2oice ! False by default, to compute H2O ice metamorphism. Read in "run_PEM.def"
 logical :: metam_co2ice ! False by default, to compute CO2 ice metamorphism. Read in "run_PEM.def"
Index: trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90	(revision 3308)
@@ -536,5 +536,5 @@
 
 !------ Sublimation of CO2 ice + Condensation of H2O ice
-    else if (htend_co2ice < 0. .and. htend_h2oice >= 0. ) then
+    else if (htend_co2ice < 0. .and. htend_h2oice > 0.) then
         write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
@@ -577,5 +577,5 @@
 
 !------ Sublimation of CO2 ice + H2O ice
-    else if (htend_co2ice < 0. .and. htend_h2oice < 0.) then
+    else if ((htend_co2ice <= 0. .and. htend_h2oice < 0.) .or. (htend_co2ice < 0. .and. htend_h2oice <= 0.)) then
         write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice'
         ! CO2 ice sublimation in the considered stratum + New stratum for dust lag
@@ -633,5 +633,5 @@
 
 !------ Condensation of CO2 ice + Sublimation of H2O ice
-    else if (htend_co2ice >= 0. .and. htend_h2oice < 0.) then
+    else if (htend_co2ice > 0. .and. htend_h2oice < 0.) then
         error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!'
     endif
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3308)
@@ -351,5 +351,5 @@
     status = nf90_close(ncid)
 
-    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+    call iniphysiq('startfi_evol.nc',iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
 #else
     ps_start_PCM(1) = ps(1)
@@ -998,6 +998,6 @@
     ! H2O ice metamorphism
     if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
-        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
-        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
+        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
     endif
 
@@ -1015,6 +1015,6 @@
     ! CO2 ice metamorphism
     if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
-        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
-        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
+        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
     endif
 enddo
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3300)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3308)
@@ -158,5 +158,5 @@
         write(*,*)'with default value ''ini_huge_h2oice'''
         do ig = 1,ngrid
-            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+            if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
         enddo
     else
@@ -358,5 +358,5 @@
 
     ! Stratification (layerings)
-    write(*,*)'So the stratification (layerings) is initialized with nothing more than the 3 sub-surface strata.'
+    write(*,*)'So the stratification (layerings) is initialized with only the 3 sub-surface strata.'
     nb_str_max = 3
 
@@ -365,5 +365,5 @@
     write(*,*)'So ''h2o_ice'' is initialized with default value ''ini_huge_h2oice'' where ''watercaptag'' is true.'
     do ig = 1,ngrid
-        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice/subslope_dist(ig,:)*cos(pi*def_slope_mean(:)*180.)
+        if (watercaptag(ig)) h2o_ice(ig,:) = ini_huge_h2oice
     enddo