Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3173)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3178)
@@ -199,2 +199,7 @@
 - Correction of script "launch_orb_1Dchained.sh" which read orbital parameters missing one interval of years in two.
 - Addition of the Martian date in "info_PEM.txt" for post-processing.
+
+== 17/01/2023 == LL
+Cleaning of the several subroutine regarding soil temperatures: they are now
+gathered in an unique Tsoil module.
+(cosmetic commit)
Index: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3178)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3178)
@@ -0,0 +1,225 @@
+MODULE compute_soiltemp_mod
+
+implicit none
+!-----------------------------------------------------------------------
+!  Author: LL
+!  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation. 
+!  
+!  Note: depths of layers and mid-layers, soil thermal inertia and 
+!        heat capacity are commons in comsoil_PEM.h
+!-----------------------------------------------------------------------
+contains
+!=======================================================================
+
+
+
+SUBROUTINE compute_tsoil_pem(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
+
+use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
+use comsoil_h,     only: volcapa
+
+implicit none
+
+!-----------------------------------------------------------------------
+!  Author: LL
+!  Purpose: Compute soil temperature using an implict 1st order scheme
+!  
+!  Note: depths of layers and mid-layers, soil thermal inertia and 
+!        heat capacity are commons in comsoil_PEM.h
+!-----------------------------------------------------------------------
+
+#include "dimensions.h"
+
+!-----------------------------------------------------------------------
+!  arguments
+!  ---------
+!  inputs:
+integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points 
+integer,                      intent(in) :: nsoil     ! number of soil layers  
+logical,                      intent(in) :: firstcall ! identifier for initialization call 
+real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
+real,                         intent(in) :: timestep  ! time step [s]
+real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
+ 
+! outputs:
+real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
+! local variables:
+integer :: ig, ik    
+
+! 0. Initialisations and preprocessing step
+ if (firstcall) then
+! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
+    do ig = 1,ngrid
+        do ik = 0,nsoil - 1
+            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa    
+        enddo
+    enddo
+
+! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
+    do ig = 1,ngrid
+        do ik = 1,nsoil - 1
+            thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
+                                   + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
+        enddo
+    enddo
+
+! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
+    ! mu_PEM
+    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
+
+    ! q_{1/2}
+    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
+    ! q_{k+1/2}
+    do ik = 1,nsoil - 1
+        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
+    enddo
+
+    do ig = 1,ngrid
+        ! d_k
+        do ik = 1,nsoil - 1
+            coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1))
+        enddo
+
+        ! alph_PEM_{N-1}
+        alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
+        ! alph_PEM_k
+        do ik = nsoil - 2,1,-1
+            alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+        enddo
+      enddo ! of do ig=1,ngrid
+endif ! of if (firstcall)
+
+IF (.not. firstcall) THEN
+! 2. Compute soil temperatures
+! First layer:
+    do ig = 1,ngrid
+        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
+                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
+
+! Other layers:
+        do ik = 1,nsoil - 1
+                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) 
+        enddo
+    enddo
+endif
+
+!  2. Compute beta_PEM coefficients (preprocessing for next time step)
+! Bottom layer, beta_PEM_{N-1}
+do ig = 1,ngrid
+    beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
+                             + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
+enddo
+! Other layers
+do ik = nsoil-2,1,-1
+    do ig = 1,ngrid
+        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ &
+                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+    enddo
+enddo
+
+END SUBROUTINE compute_tsoil_pem
+
+
+
+SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
+
+use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
+use comsoil_h,     only: volcapa
+
+implicit none
+
+!-----------------------------------------------------------------------
+!  Author: LL
+!  Purpose: Initialize the soil with the solution of the stationnary problem of Heat Conduction. Boundarry conditions: Tsurf averaged from the PCM; Geothermal flux at the bottom layer
+!
+!  Note: depths of layers and mid-layers, soil thermal inertia and
+!        heat capacity are commons in comsoil_PEM.h
+!-----------------------------------------------------------------------
+
+#include "dimensions.h"
+
+!-----------------------------------------------------------------------
+!  arguments
+!  ---------
+!  inputs:
+integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
+integer,                      intent(in) :: nsoil   ! number of soil layers
+real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
+real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
+
+! outputs:
+real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
+! local variables:
+integer :: ig, ik, iloop
+
+! 0. Initialisations and preprocessing step
+! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
+do ig = 1,ngrid
+    do ik = 0,nsoil - 1
+        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
+    enddo
+enddo
+
+! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
+do ig = 1,ngrid
+    do ik = 1,nsoil - 1
+        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
+                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
+    enddo
+enddo
+
+! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
+! mu_PEM
+mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
+
+! q_{1/2}
+coefq_PEM(:) = 0.
+! q_{k+1/2}
+
+do ig = 1,ngrid
+    ! d_k
+    do ik = 1,nsoil-1
+        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
+    enddo
+
+    ! alph_PEM_{N-1}
+    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
+    ! alph_PEM_k
+    do ik = nsoil - 2,1,-1
+        alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+    enddo
+enddo ! of do ig=1,ngrid
+
+!  1. Compute beta_PEM coefficients
+! Bottom layer, beta_PEM_{N-1}
+do ig = 1,ngrid
+        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
+                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
+enddo
+! Other layers
+do ik = nsoil - 2,1,-1
+    do ig = 1,ngrid
+        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
+                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
+    enddo
+enddo
+
+!  2. Compute soil temperatures
+do iloop = 1,10 !just convergence
+    ! First layer:
+    do ig = 1,ngrid
+        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
+                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
+    ! Other layers:
+    do ik = 1,nsoil - 1
+        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
+    enddo
+enddo
+
+enddo ! iloop
+
+END SUBROUTINE ini_tsoil_pem
+
+
+
+END MODULE compute_soiltemp_mod
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3173)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3178)
@@ -67,5 +67,5 @@
 use read_data_PCM_mod,          only: read_data_PCM
 use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
-use soil_pem_compute_mod,       only: soil_pem_compute
+use compute_soiltemp_mod,       only: compute_tsoil_pem
 use writediagpem_mod,           only: writediagfipem, writediagsoilpem
 
@@ -819,6 +819,6 @@
                 Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t)
                 Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t)
-                call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
-                call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
+                call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
                 tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope
                 do ig = 1,ngrid
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3173)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3178)
@@ -21,6 +21,5 @@
 use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
 use abort_pem_mod,              only: abort_pem
-use soil_pem_compute_mod,       only: soil_pem_compute
-use soil_pem_ini_mod,           only: soil_pem_ini
+use compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
 use comslope_mod,               only: def_slope_mean, subslope_dist
 
@@ -244,14 +243,14 @@
 !                    enddo
 !                enddo
-                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
-                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+                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))
             else
 ! predictor corrector: restart from year 1 of the PCM and build the evolution of
 ! tsoil at depth
                 tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-                call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
                 tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
-                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
 
                 do iloop = nsoil_PCM+1,nsoil_PEM
@@ -283,5 +282,5 @@
                 call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
                 do islope = 1,nslope
-                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
                 enddo
             endif
@@ -418,6 +417,6 @@
 !                enddo
 !            enddo
-            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
-            call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+            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))
 
 ! First raw initialization
@@ -430,5 +429,5 @@
             do it = 1,timelen
                 do isoil = nsoil_PCM + 1,nsoil_PEM
-                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
+                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
                 enddo
             enddo
@@ -448,5 +447,5 @@
             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
             do islope = 1,nslope
-                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
+                call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             enddo
             write(*,*) 'PEMETAT0: Ice table done'
Index: trunk/LMDZ.COMMON/libf/evolution/soil_pem_compute_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_pem_compute_mod.F90	(revision 3173)
+++ 	(revision )
@@ -1,115 +1,0 @@
-MODULE soil_pem_compute_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE soil_pem_compute(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
-
-use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
-use comsoil_h,     only: volcapa
-
-implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme
-!  
-!  Note: depths of layers and mid-layers, soil thermal inertia and 
-!        heat capacity are commons in comsoil_PEM.h
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points 
-integer,                      intent(in) :: nsoil     ! number of soil layers  
-logical,                      intent(in) :: firstcall ! identifier for initialization call 
-real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
-real,                         intent(in) :: timestep  ! time step [s]
-real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
- 
-! outputs:
-real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
-! local variables:
-integer :: ig, ik    
-
-! 0. Initialisations and preprocessing step
- if (firstcall) then
-! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
-    do ig = 1,ngrid
-        do ik = 0,nsoil - 1
-            mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa    
-        enddo
-    enddo
-
-! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
-    do ig = 1,ngrid
-        do ik = 1,nsoil - 1
-            thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
-                                   + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
-        enddo
-    enddo
-
-! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
-    ! mu_PEM
-    mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
-
-    ! q_{1/2}
-    coefq_PEM(0) = volcapa*layer_PEM(1)/timestep
-    ! q_{k+1/2}
-    do ik = 1,nsoil - 1
-        coefq_PEM(ik) = volcapa*(layer_PEM(ik + 1) - layer_PEM(ik))/timestep
-    enddo
-
-    do ig = 1,ngrid
-        ! d_k
-        do ik = 1,nsoil - 1
-            coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik)-mlayer_PEM(ik - 1))
-        enddo
-
-        ! alph_PEM_{N-1}
-        alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil-1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
-        ! alph_PEM_k
-        do ik = nsoil - 2,1,-1
-            alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
-        enddo
-      enddo ! of do ig=1,ngrid
-endif ! of if (firstcall)
-
-IF (.not. firstcall) THEN
-! 2. Compute soil temperatures
-! First layer:
-    do ig = 1,ngrid
-        tsoil(ig,1) = (tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
-                      (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
-
-! Other layers:
-        do ik = 1,nsoil - 1
-                tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik) 
-        enddo
-    enddo
-endif
-
-!  2. Compute beta_PEM coefficients (preprocessing for next time step)
-! Bottom layer, beta_PEM_{N-1}
-do ig = 1,ngrid
-    beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
-                             + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
-enddo
-! Other layers
-do ik = nsoil-2,1,-1
-    do ig = 1,ngrid
-        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik + 1)*beta_PEM(ig,ik + 1))/ &
-                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
-    enddo
-enddo
-
-END SUBROUTINE soil_pem_compute
-
-END MODULE soil_pem_compute_mod
Index: trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/soil_pem_ini_mod.F90	(revision 3173)
+++ 	(revision )
@@ -1,108 +1,0 @@
-MODULE soil_pem_ini_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE soil_pem_ini(ngrid,nsoil,therm_i,tsurf,tsoil)
-
-use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
-use comsoil_h,     only: volcapa
-
-implicit none
-
-!-----------------------------------------------------------------------
-!  Author: LL
-!  Purpose: Compute soil temperature using an implict 1st order scheme, stationnary solution
-!
-!  Note: depths of layers and mid-layers, soil thermal inertia and
-!        heat capacity are commons in comsoil_PEM.h
-!-----------------------------------------------------------------------
-
-#include "dimensions.h"
-
-!-----------------------------------------------------------------------
-!  arguments
-!  ---------
-!  inputs:
-integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
-integer,                      intent(in) :: nsoil   ! number of soil layers
-real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
-real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
-
-! outputs:
-real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
-! local variables:
-integer :: ig, ik, iloop
-
-! 0. Initialisations and preprocessing step
-! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
-do ig = 1,ngrid
-    do ik = 0,nsoil - 1
-        mthermdiff_PEM(ig,ik) = therm_i(ig,ik + 1)*therm_i(ig,ik + 1)/volcapa
-    enddo
-enddo
-
-! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
-do ig = 1,ngrid
-    do ik = 1,nsoil - 1
-        thermdiff_PEM(ig,ik) = ((layer_PEM(ik) - mlayer_PEM(ik - 1))*mthermdiff_PEM(ig,ik) &
-                             + (mlayer_PEM(ik) - layer_PEM(ik))*mthermdiff_PEM(ig,ik - 1))/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
-    enddo
-enddo
-
-! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alph_PEMa_k and capcal
-! mu_PEM
-mu_PEM = mlayer_PEM(0)/(mlayer_PEM(1) - mlayer_PEM(0))
-
-! q_{1/2}
-coefq_PEM(:) = 0.
-! q_{k+1/2}
-
-do ig = 1,ngrid
-    ! d_k
-    do ik = 1,nsoil-1
-        coefd_PEM(ig,ik) = thermdiff_PEM(ig,ik)/(mlayer_PEM(ik) - mlayer_PEM(ik - 1))
-    enddo
-
-    ! alph_PEM_{N-1}
-    alph_PEM(ig,nsoil - 1) = coefd_PEM(ig,nsoil - 1)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
-    ! alph_PEM_k
-    do ik = nsoil - 2,1,-1
-        alph_PEM(ig,ik) = coefd_PEM(ig,ik)/(coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
-    enddo
-enddo ! of do ig=1,ngrid
-
-!  1. Compute beta_PEM coefficients
-! Bottom layer, beta_PEM_{N-1}
-do ig = 1,ngrid
-        beta_PEM(ig,nsoil - 1) = coefq_PEM(nsoil - 1)*tsoil(ig,nsoil)/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1)) &
-                                 + fluxgeo/(coefq_PEM(nsoil - 1) + coefd_PEM(ig,nsoil - 1))
-enddo
-! Other layers
-do ik = nsoil - 2,1,-1
-    do ig = 1,ngrid
-        beta_PEM(ig,ik) = (coefq_PEM(ik)*tsoil(ig,ik + 1) + coefd_PEM(ig,ik+1)*beta_PEM(ig,ik + 1))/ &
-                          (coefq_PEM(ik) + coefd_PEM(ig,ik + 1)*(1. - alph_PEM(ig,ik + 1)) + coefd_PEM(ig,ik))
-    enddo
-enddo
-
-!  2. Compute soil temperatures
-do iloop = 1,10 !just convergence
-    ! First layer:
-    do ig = 1,ngrid
-        tsoil(ig,1)=(tsurf(ig) + mu_PEM*beta_PEM(ig,1)*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))/ &
-                   (1. + mu_PEM*(1. - alph_PEM(ig,1))*thermdiff_PEM(ig,1)/mthermdiff_PEM(ig,0))
-    ! Other layers:
-    do ik = 1,nsoil - 1
-        tsoil(ig,ik + 1) = alph_PEM(ig,ik)*tsoil(ig,ik) + beta_PEM(ig,ik)
-    enddo
-enddo
-
-enddo ! iloop
-
-END SUBROUTINE soil_pem_ini
-
-END MODULE soil_pem_ini_mod
