Index: LMDZ6/trunk/libf/phylmd/inlandsis/VAR0SV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VAR0SV.F90	(revision 5259)
+++ 	(revision )
@@ -1,33 +1,0 @@
-MODULE VAR0SV
-
-USE VAR_SV, only : nsol, nsno 
-USE VARdSV, only : nsot
-
-IMPLICIT NONE
-                                                  
-      INTEGER    :: islpSV(-nsol:0)                                           
-      INTEGER    :: isnpSV(      nsno)                                        
-      INTEGER    :: islmSV(-nsol:0)    
-                                      
-      INTEGER,PARAMETER :: nkhy=50                                            
-      REAL       :: Implic,Explic                                             
-      REAL       :: dzmiSV(-nsol:0)          ! dz_(i-1/2)                     
-      REAL       :: dzi_SV(-nsol:0)          ! dz_(i-1)/(dz_(i)+dz_(i-1))     
-      REAL       :: dziiSV(-nsol:0)          ! dz_(i)  /(dz_(i)+dz_(i-1))     
-      REAL       :: dtz_SV(-nsol:0)          ! dt / dz             
-      REAL       :: dtz_SV2(-nsol:0)          ! dt / dz             
-      REAL       :: dz78SV(-nsol:0)          ! 7/8 (dz)                       
-      REAL       :: dz34SV(-nsol:0)          ! 3/4 (dz)                       
-      REAL       :: dz_8SV(-nsol:0)          ! 1/8 (dz)                       
-      REAL       :: dzAvSV(-nsol:0)          ! 1/8dz_(-1)+3/4dz+1/8dz_(+1)    
-      REAL       :: OcndSV                   ! Swab Ocean / Soil Ratio        
-      REAL       :: rocsSV( 0:nsot)          ! Soil Contribution to (ro c)_s  
-      REAL       :: etamSV( 0:nsot)          ! Soil Minimum Humidity          
-      REAL       :: s1__SV( 0:nsot)          ! ... X eta**( b+2), DR97(3.36)  
-      REAL       :: s2__SV( 0:nsot)          ! ... X eta**(2b+3), DR97(3.35)  
-      REAL       :: aKdtSV( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: a * dt           
-      REAL       :: bKdtSV( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: b * dt     
-      REAL       :: aKdtSV2( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: a * dt           
-      REAL       :: bKdtSV2( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: b * dt                 
-
-END MODULE VAR0SV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VAR0SV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VAR0SV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VAR0SV.f90	(revision 5268)
@@ -0,0 +1,33 @@
+MODULE VAR0SV
+
+USE VAR_SV, only : nsol, nsno 
+USE VARdSV, only : nsot
+
+IMPLICIT NONE
+                                                  
+      INTEGER    :: islpSV(-nsol:0)                                           
+      INTEGER    :: isnpSV(      nsno)                                        
+      INTEGER    :: islmSV(-nsol:0)    
+                                      
+      INTEGER,PARAMETER :: nkhy=50                                            
+      REAL       :: Implic,Explic                                             
+      REAL       :: dzmiSV(-nsol:0)          ! dz_(i-1/2)                     
+      REAL       :: dzi_SV(-nsol:0)          ! dz_(i-1)/(dz_(i)+dz_(i-1))     
+      REAL       :: dziiSV(-nsol:0)          ! dz_(i)  /(dz_(i)+dz_(i-1))     
+      REAL       :: dtz_SV(-nsol:0)          ! dt / dz             
+      REAL       :: dtz_SV2(-nsol:0)          ! dt / dz             
+      REAL       :: dz78SV(-nsol:0)          ! 7/8 (dz)                       
+      REAL       :: dz34SV(-nsol:0)          ! 3/4 (dz)                       
+      REAL       :: dz_8SV(-nsol:0)          ! 1/8 (dz)                       
+      REAL       :: dzAvSV(-nsol:0)          ! 1/8dz_(-1)+3/4dz+1/8dz_(+1)    
+      REAL       :: OcndSV                   ! Swab Ocean / Soil Ratio        
+      REAL       :: rocsSV( 0:nsot)          ! Soil Contribution to (ro c)_s  
+      REAL       :: etamSV( 0:nsot)          ! Soil Minimum Humidity          
+      REAL       :: s1__SV( 0:nsot)          ! ... X eta**( b+2), DR97(3.36)  
+      REAL       :: s2__SV( 0:nsot)          ! ... X eta**(2b+3), DR97(3.35)  
+      REAL       :: aKdtSV( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: a * dt           
+      REAL       :: bKdtSV( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: b * dt     
+      REAL       :: aKdtSV2( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: a * dt           
+      REAL       :: bKdtSV2( 0:nsot, 0:nkhy)  ! Khyd=a*eta+b: b * dt                 
+
+END MODULE VAR0SV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VAR_SV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VAR_SV.F90	(revision 5259)
+++ 	(revision )
@@ -1,26 +1,0 @@
-MODULE VAR_SV
-
-
-IMPLICIT NONE
-
-INCLUDE "dimsoil.h"
-
-! +
-! +   LMDZ_SV dims
-! +   ===================
-! +
-!      INTEGER,PARAMETER :: nsol=nsoilmx-1
-!      INTEGER,PARAMETER :: nsno=nsnowmx
-     
-      INTEGER, PARAMETER    :: nsol=nsoilmx-1
-      INTEGER, PARAMETER    :: nsot=nsoilmx+1
-      INTEGER, PARAMETER    :: nsno=nsnowmx
-      INTEGER, PARAMETER    :: nb_wri=1
-      INTEGER,SAVE      :: klonv
-!$OMP THREADPRIVATE(klonv)
-      INTEGER,SAVE      :: knonv
-!$OMP THREADPRIVATE(knonv)
-
-
-
-END MODULE VAR_SV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VAR_SV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VAR_SV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VAR_SV.f90	(revision 5268)
@@ -0,0 +1,26 @@
+MODULE VAR_SV
+
+
+IMPLICIT NONE
+
+INCLUDE "dimsoil.h"
+
+! +
+! +   LMDZ_SV dims
+! +   ===================
+! +
+!      INTEGER,PARAMETER :: nsol=nsoilmx-1
+!      INTEGER,PARAMETER :: nsno=nsnowmx
+     
+      INTEGER, PARAMETER    :: nsol=nsoilmx-1
+      INTEGER, PARAMETER    :: nsot=nsoilmx+1
+      INTEGER, PARAMETER    :: nsno=nsnowmx
+      INTEGER, PARAMETER    :: nb_wri=1
+      INTEGER,SAVE      :: klonv
+!$OMP THREADPRIVATE(klonv)
+      INTEGER,SAVE      :: knonv
+!$OMP THREADPRIVATE(knonv)
+
+
+
+END MODULE VAR_SV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARdSV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARdSV.F90	(revision 5259)
+++ 	(revision )
@@ -1,101 +1,0 @@
-MODULE VARdSV
-
-USE VAR_SV
-
-IMPLICIT NONE
-
-
-! +--SISVAT Global Variables
-! +  =======================
-
-      LOGICAL           ::   INI_SV=.false.          ! Initialisation Switch     
-      REAL,PARAMETER    ::   eps_21=1.e-21           ! Arbitrary  very small value
-
-
-! +--Snow
-! +  ----
-
-      INTEGER,PARAMETER,DIMENSION(5) ::   istdSV(1:5)=(/1,2,3,4,5/)  ! Snow History
-
-      REAL,PARAMETER    ::   Cn_dSV= 2105.           ! Snow Heat Capacity          [J/kg/K]
-      REAL,PARAMETER    ::   SMndSV= 1.00            ! Minimum Thickness of new Layers
-      REAL,PARAMETER    ::   G1_dSV= 99.             ! Conversion 0/99-->0/1 
-      REAL,PARAMETER    ::   DDcdSV= 1.,DFcdSV= 4.,DScdSV= 3.   
-                                                     ! Snow Grains Optical Diameter [1e-4m]
-      REAL,PARAMETER    ::   ADSdSV= 4.              ! Snow Grains Actual  Diameter [1e-4m]
-      REAL,PARAMETER    ::   So1dSV= 0.580,So2dSV= 0.320,So3dSV= 0.100
-                                                     ! Total Solar Irradiance Fractions [-]
-                                                     ! Tuning ETH camp 0.3--0.8mim Interval
-                                                     ! Tuning ETH camp 0.8--1.5mim Interval
-                                                     ! Tuning ETH camp 1.5--2.8mim Interval
-                                                     ! So1dSV=0.606,So2dSV=0.301,So3dSV=0.093
-      REAL,PARAMETER    ::   aI1dSV= 0.40,aI2dSV= 0.45,aI3dSV= 0.65
-                                                     ! Bare Ice Albedo                  [-]
-                                                     ! Minimum/Maximum/ICE lense albedo at 
-                                                     ! 800 kg/m3 and minimum pure snow albedo
-      REAL,PARAMETER    ::   ws0dSV= 0.07            ! Irreducible Water Saturation in Snow
-      REAL,PARAMETER    ::   roCdSV= 800.            ! Pore Hole Close OFF Density  [kg/m3]
-      REAL,PARAMETER    ::   roBdSV= 450.            ! Max blowing snow Density [kg/m3]
-      REAL,PARAMETER    ::   roSdSV= 450.            ! Max pure    snow Density [kg/m3]
-      REAL,PARAMETER    ::   ru_dSV= 200.            ! Surficial Water Scale Factor [kg/m2]
-
-!C +--Ice
-!C +  ---
-
-      REAL,PARAMETER    ::   CdidSV= 2.1
-
-
-
-!C +--Soil
-!C +  ----
-
-      REAL,PARAMETER    ::   rcwdSV = 4.180e+6    ! Density * Water Specific Heat
-      REAL              ::   dz_dSV(-nsol:0)      ! Vertical  Discretization MARSV:
-                                                  !/0.72,0.20,0.060,0.019,0.001/ 
-                                                  ! Layer's Thickness
-      REAL              ::   zz_dSV               ! Soil      Thickness
-
-      REAL,PARAMETER    ::   etadSV(0:nsot) = (/ 1.000,0.395,0.410,     &
-     &     0.435,0.485,0.451,0.420,0.477,0.476,0.426,0.492,0.482,0.001 /)      
-                                                  ! Water Content at Saturation  [m3/m3]
-
-      REAL,PARAMETER    ::   psidSV(0:nsot) = (/ 1.000,0.121,0.090,     &
-     &     0.218,0.786,0.478,0.299,0.356,0.630,0.153,0.490,0.405,0.001 /)
-                                                  ! Water Succion at Saturation      [m]
-
-      REAL,PARAMETER    ::   Ks_dSV(0:nsot) = (/ 0.e00, 176.0e-6,       &
-     &          156.3e-6,  34.1e-6,   7.2e-6,   7.0e-6,   6.3e-6,       &
-     &            1.7e-6,   2.5e-6,   2.2e-6,   1.0e-6,   1.3e-6,0.0e0 /)
-                                                  ! Hydraulic Conductivity
-                                                  !               at Saturation    [m/s]
-      REAL,PARAMETER    ::   bCHdSV(0:nsot) = (/ 1.00, 4.05, 4.38,      &
-     &      4.90, 5.30, 5.39, 7.12, 7.75, 8.52,10.40,10.40,11.40, 0.02 /)
-                                                  ! Clapp-Hornberger Coefficient b   [-]
-
- !     etadSV,   psidSV,   Ks_dSV    bCHdSV  
- !     /1.000,    1.000,   0.0e00,     1.00,      !  0 WATER
- !      0.395,    0.121, 176.0e-6,     4.05,      !  1 SAND
- !      0.410,    0.090, 156.3e-6,     4.38,      !  2 LOAMY      SAND
- !      0.435,    0.218,  34.1e-6,     4.90,      !  3 SANDY      LOAM
- !      0.485,    0.786,   7.2e-6,     5.30,      !  4 SILT       LOAM
- !      0.451,    0.478,   7.0e-6,     5.39,      !  5            LOAM
- !      0.420,    0.299,   6.3e-6,     7.12,      !  6 SANDY CLAY LOAM
- !      0.477,    0.356,   1.7e-6,     7.75,      !  7 SILTY CLAY LOAM 
- !      0.476,    0.630,   2.5e-6,     8.52,      !  8       CLAY LOAM
- !      0.426,    0.153,   2.2e-6,    10.40,      !  9 SANDY CLAY
- !      0.492,    0.490,   1.0e-6,    10.40,      ! 10 SILTY CLAY
- !      0.482,    0.405,   1.3e-6,    11.40,      ! 11       CLAY
- !      0.001,    0.001,   0.0e00,     0.02/      ! 12       ICE 
-
-
-!C +--Water Bodies
-!C +  ------------
-
-      REAL,PARAMETER    ::   vK_dSV = 1000.       ! Diffusivity in Water          [m2/s]
-      REAL,PARAMETER    ::   TSIdSV = 0.50        ! Sea-Ice Fraction: SST Scale      [K]
-
-      INTEGER ::  iso1  !rajout hjp for ini - check if ivg,iso is possible
-
-!
-
-END MODULE VARdSV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARdSV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARdSV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VARdSV.f90	(revision 5268)
@@ -0,0 +1,101 @@
+MODULE VARdSV
+
+USE VAR_SV
+
+IMPLICIT NONE
+
+
+! +--SISVAT Global Variables
+! +  =======================
+
+      LOGICAL           ::   INI_SV=.false.          ! Initialisation Switch     
+      REAL,PARAMETER    ::   eps_21=1.e-21           ! Arbitrary  very small value
+
+
+! +--Snow
+! +  ----
+
+      INTEGER,PARAMETER,DIMENSION(5) ::   istdSV(1:5)=(/1,2,3,4,5/)  ! Snow History
+
+      REAL,PARAMETER    ::   Cn_dSV= 2105.           ! Snow Heat Capacity          [J/kg/K]
+      REAL,PARAMETER    ::   SMndSV= 1.00            ! Minimum Thickness of new Layers
+      REAL,PARAMETER    ::   G1_dSV= 99.             ! Conversion 0/99-->0/1 
+      REAL,PARAMETER    ::   DDcdSV= 1.,DFcdSV= 4.,DScdSV= 3.   
+                                                     ! Snow Grains Optical Diameter [1e-4m]
+      REAL,PARAMETER    ::   ADSdSV= 4.              ! Snow Grains Actual  Diameter [1e-4m]
+      REAL,PARAMETER    ::   So1dSV= 0.580,So2dSV= 0.320,So3dSV= 0.100
+                                                     ! Total Solar Irradiance Fractions [-]
+                                                     ! Tuning ETH camp 0.3--0.8mim Interval
+                                                     ! Tuning ETH camp 0.8--1.5mim Interval
+                                                     ! Tuning ETH camp 1.5--2.8mim Interval
+                                                     ! So1dSV=0.606,So2dSV=0.301,So3dSV=0.093
+      REAL,PARAMETER    ::   aI1dSV= 0.40,aI2dSV= 0.45,aI3dSV= 0.65
+                                                     ! Bare Ice Albedo                  [-]
+                                                     ! Minimum/Maximum/ICE lense albedo at 
+                                                     ! 800 kg/m3 and minimum pure snow albedo
+      REAL,PARAMETER    ::   ws0dSV= 0.07            ! Irreducible Water Saturation in Snow
+      REAL,PARAMETER    ::   roCdSV= 800.            ! Pore Hole Close OFF Density  [kg/m3]
+      REAL,PARAMETER    ::   roBdSV= 450.            ! Max blowing snow Density [kg/m3]
+      REAL,PARAMETER    ::   roSdSV= 450.            ! Max pure    snow Density [kg/m3]
+      REAL,PARAMETER    ::   ru_dSV= 200.            ! Surficial Water Scale Factor [kg/m2]
+
+!C +--Ice
+!C +  ---
+
+      REAL,PARAMETER    ::   CdidSV= 2.1
+
+
+
+!C +--Soil
+!C +  ----
+
+      REAL,PARAMETER    ::   rcwdSV = 4.180e+6    ! Density * Water Specific Heat
+      REAL              ::   dz_dSV(-nsol:0)      ! Vertical  Discretization MARSV:
+                                                  !/0.72,0.20,0.060,0.019,0.001/ 
+                                                  ! Layer's Thickness
+      REAL              ::   zz_dSV               ! Soil      Thickness
+
+      REAL,PARAMETER    ::   etadSV(0:nsot) = (/ 1.000,0.395,0.410,     &
+     &     0.435,0.485,0.451,0.420,0.477,0.476,0.426,0.492,0.482,0.001 /)      
+                                                  ! Water Content at Saturation  [m3/m3]
+
+      REAL,PARAMETER    ::   psidSV(0:nsot) = (/ 1.000,0.121,0.090,     &
+     &     0.218,0.786,0.478,0.299,0.356,0.630,0.153,0.490,0.405,0.001 /)
+                                                  ! Water Succion at Saturation      [m]
+
+      REAL,PARAMETER    ::   Ks_dSV(0:nsot) = (/ 0.e00, 176.0e-6,       &
+     &          156.3e-6,  34.1e-6,   7.2e-6,   7.0e-6,   6.3e-6,       &
+     &            1.7e-6,   2.5e-6,   2.2e-6,   1.0e-6,   1.3e-6,0.0e0 /)
+                                                  ! Hydraulic Conductivity
+                                                  !               at Saturation    [m/s]
+      REAL,PARAMETER    ::   bCHdSV(0:nsot) = (/ 1.00, 4.05, 4.38,      &
+     &      4.90, 5.30, 5.39, 7.12, 7.75, 8.52,10.40,10.40,11.40, 0.02 /)
+                                                  ! Clapp-Hornberger Coefficient b   [-]
+
+ !     etadSV,   psidSV,   Ks_dSV    bCHdSV  
+ !     /1.000,    1.000,   0.0e00,     1.00,      !  0 WATER
+ !      0.395,    0.121, 176.0e-6,     4.05,      !  1 SAND
+ !      0.410,    0.090, 156.3e-6,     4.38,      !  2 LOAMY      SAND
+ !      0.435,    0.218,  34.1e-6,     4.90,      !  3 SANDY      LOAM
+ !      0.485,    0.786,   7.2e-6,     5.30,      !  4 SILT       LOAM
+ !      0.451,    0.478,   7.0e-6,     5.39,      !  5            LOAM
+ !      0.420,    0.299,   6.3e-6,     7.12,      !  6 SANDY CLAY LOAM
+ !      0.477,    0.356,   1.7e-6,     7.75,      !  7 SILTY CLAY LOAM 
+ !      0.476,    0.630,   2.5e-6,     8.52,      !  8       CLAY LOAM
+ !      0.426,    0.153,   2.2e-6,    10.40,      !  9 SANDY CLAY
+ !      0.492,    0.490,   1.0e-6,    10.40,      ! 10 SILTY CLAY
+ !      0.482,    0.405,   1.3e-6,    11.40,      ! 11       CLAY
+ !      0.001,    0.001,   0.0e00,     0.02/      ! 12       ICE 
+
+
+!C +--Water Bodies
+!C +  ------------
+
+      REAL,PARAMETER    ::   vK_dSV = 1000.       ! Diffusivity in Water          [m2/s]
+      REAL,PARAMETER    ::   TSIdSV = 0.50        ! Sea-Ice Fraction: SST Scale      [K]
+
+      INTEGER ::  iso1  !rajout hjp for ini - check if ivg,iso is possible
+
+!
+
+END MODULE VARdSV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.F90	(revision 5259)
+++ 	(revision )
@@ -1,103 +1,0 @@
-MODULE VARphy
-
-
-IMPLICIT NONE
-
-
-
-! Etienne, 08112020: merging YOMCST_SISVAT,  PHY_SV and VARphy
-! but redundancies remain
-!
-!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
-!                 veillez à n'utiliser que des ! pour les commentaires
-!                 et à bien positionner les & des lignes de continuation
-!                 (les placer en colonne 6 et en colonne 73)
-!
-!
-! A1.0 Fundamental constants
-!----------------------------
-! Constants
-                                                                              
-      CHARACTER(LEN=1), PARAMETER :: labnum(0:9) =                            &
-              (/'0','1','2','3','4','5','6','7','8','9'/)             
-!C +...                labnum: Alphanumeric Character      
-     
-                                                                               
-      INTEGER, PARAMETER ::  iun=1                                             
-      REAL, PARAMETER    ::  zer0 = 0.0e+0, half = 0.5e+0, un_1 = 1.0e+0,     &
-     &                       eps6 = 1.0e-6, R_1000=1.e3   
-      REAL, PARAMETER    ::  zero = 0.0e+0, demi = 0.5e+0, unun = 1.0e+0,     &
-     &                       epsi = 1.0e-6, eps9 = 1.0e-9         
-      REAL               ::  ea_MAX,ea_MIN                                                                               
-      REAL, PARAMETER    ::  pi = 3.141592653589793238462643e0
-     
-! A1.2 Geoide
-!-------------
-      REAL, PARAMETER    ::  Grav_F = 9.81e0, Grav_I = 1./Grav_F                       
-!C +                   Grav_F:  Gravitational  Force         =  9.81    m/s2           
-!C +                   Grav_I: 1 / Grav_F                    = 1 /9.81    s2/m    
-
-                                                                               
-      REAL, PARAMETER    ::  gravit = 9.81e0,      gravi2 = gravit**2 ,       &
-     &                       grvinv = 1./gravit                       
-!C +                   gravit: Earth Gravity Acceleration    = 9.81    m/s2    
-!C +                   gravi2: idem (squared)                                  
-!C +                   grvinv: idem (inverse)                 
-! A1.3 Radiation
-!----------------
-      REAL,PARAMETER :: StefBo=5.67e-8  
-!C +                   StefBo: Stefan-Boltzman Constant       =    5.67d-8 W/m2/K4
-
-! A1.4 Thermodynamic 
-!-----------------------------
-     
-      REAL, PARAMETER    ::  cp=1004.708845
-                                                       
-      REAL, PARAMETER    ::  RVapor=461.e0, Lv_H2O=2.5008e+6,&
-     &                       Ls_H2O=2.8345e+6, r_LvCp=2490.04,r_LcCp=332.27,  &
-     &                       r_LsCp=2822.31            
-!C +                   cp    : dry air specific heat at constant p (1004 J/kg/K) 
-
-      REAL, PARAMETER    ::  LhfH2O=3.34e+5, LhvH2O=2.5008e+6,&
-     &                       LhsH2O=2.8345e+6     
-      REAL, PARAMETER    ::  rhoWat=1000.00e0    
-!C +                   rhoWat: Water Specific Mass            = 1000.00d+0 kg/m3 
-
-      REAL, PARAMETER    ::  rhoIce=920.e0, BSnoRo=3.30e+2, frsno=300.0                             
-!C +...                rhoIce:  Density      of Pure Ice      =  920.00d+0 kg/m3
-!C +                   BSnoRo : Blowed Snow Density           =  255.00d+0 kg/m3 
-                                                                                
-      REAL, PARAMETER    ::  Tf_Sno=273.16e+0 
-!C +...                Tf_Sno: Snow  Melting  Point           =  273.16    K     
-
-
-      REAL, PARAMETER    ::  ro_Wat=1000.00e0, C__Wat=4186.00e0
-!C +                   ro_Wat: Density       of Water         = 1000.00d+0 kg/m3 
-!C +                   C__Wat: Heat Capacity of Water         = 4186.00d+0 J/kg/K
-                                                                                
-      REAL, PARAMETER    ::  ro_Ice=920.e0, cdice=2.10e0                                 
-!C +...                ro_Ice:  Density      of Pure Ice       =  920.00d+0 kg/m3
-!C +                   cdice :  Conductivity of Pure Ice       =    2.51d+0 W/m/K
-                                                                                
-      REAL, PARAMETER    ::  TfSnow=273.15e+0,csnow=2105.00e+0,r0sno=3.00e+1, &
-     &                       blsno=3.30e+2,   Lf_H2O=3.337e+5                          
-!C +...                TfSnow:        Snow melting Temperature=  273.15d+0 K     
-!C +                    csnow:Heat Capacity of Snow             2105      J/kg/K 
-!C +                         (Loth et al. 1993, JGR 98 D6, 2.2.2 2e para p.10453)
-!C +                   r0sno : Fresh  Snow Density                50.00d+0 kg/m3 
-!C +                   blsno : Blowed Snow Density               330.00d+0 kg/m3 
-!C +                   Lf_H2O: Latent Heat of Fusion of Snow  =    3.34d+5 J/kg  
-
-! A1.6 Turbulent and molecular diffusion
-!----------------------------------------
-      REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0, r_turb=3.0
-      REAL, PARAMETER    ::  A_turb=5.8, akmol=1.35e-5
-!C +...                A_MolV: Air Viscosity                 = 1.35d-5 m2/s   
-!C +                   vonKrm: von Karman constant           = 0.4            
-!C +                   r_turb:   Turbulent Diffusivities Ratio K*/Km        
-!C +                   A_turb:    Stability  Coefficient Moment                                 
-!C +                   Air Viscosity                 = 1.35d-5 m2/s                 
-
-
-
-END MODULE VARphy
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VARphy.f90	(revision 5268)
@@ -0,0 +1,103 @@
+MODULE VARphy
+
+
+IMPLICIT NONE
+
+
+
+! Etienne, 08112020: merging YOMCST_SISVAT,  PHY_SV and VARphy
+! but redundancies remain
+!
+!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
+!                 veillez à n'utiliser que des ! pour les commentaires
+!                 et à bien positionner les & des lignes de continuation
+!                 (les placer en colonne 6 et en colonne 73)
+!
+!
+! A1.0 Fundamental constants
+!----------------------------
+! Constants
+                                                                              
+      CHARACTER(LEN=1), PARAMETER :: labnum(0:9) =                            &
+              (/'0','1','2','3','4','5','6','7','8','9'/)             
+!C +...                labnum: Alphanumeric Character      
+     
+                                                                               
+      INTEGER, PARAMETER ::  iun=1                                             
+      REAL, PARAMETER    ::  zer0 = 0.0e+0, half = 0.5e+0, un_1 = 1.0e+0,     &
+     &                       eps6 = 1.0e-6, R_1000=1.e3   
+      REAL, PARAMETER    ::  zero = 0.0e+0, demi = 0.5e+0, unun = 1.0e+0,     &
+     &                       epsi = 1.0e-6, eps9 = 1.0e-9         
+      REAL               ::  ea_MAX,ea_MIN                                                                               
+      REAL, PARAMETER    ::  pi = 3.141592653589793238462643e0
+     
+! A1.2 Geoide
+!-------------
+      REAL, PARAMETER    ::  Grav_F = 9.81e0, Grav_I = 1./Grav_F                       
+!C +                   Grav_F:  Gravitational  Force         =  9.81    m/s2           
+!C +                   Grav_I: 1 / Grav_F                    = 1 /9.81    s2/m    
+
+                                                                               
+      REAL, PARAMETER    ::  gravit = 9.81e0,      gravi2 = gravit**2 ,       &
+     &                       grvinv = 1./gravit                       
+!C +                   gravit: Earth Gravity Acceleration    = 9.81    m/s2    
+!C +                   gravi2: idem (squared)                                  
+!C +                   grvinv: idem (inverse)                 
+! A1.3 Radiation
+!----------------
+      REAL,PARAMETER :: StefBo=5.67e-8  
+!C +                   StefBo: Stefan-Boltzman Constant       =    5.67d-8 W/m2/K4
+
+! A1.4 Thermodynamic 
+!-----------------------------
+     
+      REAL, PARAMETER    ::  cp=1004.708845
+                                                       
+      REAL, PARAMETER    ::  RVapor=461.e0, Lv_H2O=2.5008e+6,&
+     &                       Ls_H2O=2.8345e+6, r_LvCp=2490.04,r_LcCp=332.27,  &
+     &                       r_LsCp=2822.31            
+!C +                   cp    : dry air specific heat at constant p (1004 J/kg/K) 
+
+      REAL, PARAMETER    ::  LhfH2O=3.34e+5, LhvH2O=2.5008e+6,&
+     &                       LhsH2O=2.8345e+6     
+      REAL, PARAMETER    ::  rhoWat=1000.00e0    
+!C +                   rhoWat: Water Specific Mass            = 1000.00d+0 kg/m3 
+
+      REAL, PARAMETER    ::  rhoIce=920.e0, BSnoRo=3.30e+2, frsno=300.0                             
+!C +...                rhoIce:  Density      of Pure Ice      =  920.00d+0 kg/m3
+!C +                   BSnoRo : Blowed Snow Density           =  255.00d+0 kg/m3 
+                                                                                
+      REAL, PARAMETER    ::  Tf_Sno=273.16e+0 
+!C +...                Tf_Sno: Snow  Melting  Point           =  273.16    K     
+
+
+      REAL, PARAMETER    ::  ro_Wat=1000.00e0, C__Wat=4186.00e0
+!C +                   ro_Wat: Density       of Water         = 1000.00d+0 kg/m3 
+!C +                   C__Wat: Heat Capacity of Water         = 4186.00d+0 J/kg/K
+                                                                                
+      REAL, PARAMETER    ::  ro_Ice=920.e0, cdice=2.10e0                                 
+!C +...                ro_Ice:  Density      of Pure Ice       =  920.00d+0 kg/m3
+!C +                   cdice :  Conductivity of Pure Ice       =    2.51d+0 W/m/K
+                                                                                
+      REAL, PARAMETER    ::  TfSnow=273.15e+0,csnow=2105.00e+0,r0sno=3.00e+1, &
+     &                       blsno=3.30e+2,   Lf_H2O=3.337e+5                          
+!C +...                TfSnow:        Snow melting Temperature=  273.15d+0 K     
+!C +                    csnow:Heat Capacity of Snow             2105      J/kg/K 
+!C +                         (Loth et al. 1993, JGR 98 D6, 2.2.2 2e para p.10453)
+!C +                   r0sno : Fresh  Snow Density                50.00d+0 kg/m3 
+!C +                   blsno : Blowed Snow Density               330.00d+0 kg/m3 
+!C +                   Lf_H2O: Latent Heat of Fusion of Snow  =    3.34d+5 J/kg  
+
+! A1.6 Turbulent and molecular diffusion
+!----------------------------------------
+      REAL, PARAMETER    ::  A_MolV = 1.35e-5, vonKrm = 0.40e0, r_turb=3.0
+      REAL, PARAMETER    ::  A_turb=5.8, akmol=1.35e-5
+!C +...                A_MolV: Air Viscosity                 = 1.35d-5 m2/s   
+!C +                   vonKrm: von Karman constant           = 0.4            
+!C +                   r_turb:   Turbulent Diffusivities Ratio K*/Km        
+!C +                   A_turb:    Stability  Coefficient Moment                                 
+!C +                   Air Viscosity                 = 1.35d-5 m2/s                 
+
+
+
+END MODULE VARphy
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.F90	(revision 5259)
+++ 	(revision )
@@ -1,92 +1,0 @@
-MODULE VARtSV
-
-USE VAR_SV, only : nsol, nsno, klonv
-
-IMPLICIT NONE
-! +
-      INTEGER, PARAMETER      ::   ntaver = 4
-! +
-      REAL,ALLOCATABLE,SAVE                ::   toicSV(:)    ! Snow to ice balance    
-!$OMP THREADPRIVATE(toicSV)
-
-      REAL,ALLOCATABLE,SAVE                ::   dz1_SV(:,:) ! "inverse" layer thicknes
-!$OMP THREADPRIVATE(dz1_SV)
-      REAL,ALLOCATABLE,SAVE                ::   dz2_SV(:,:) ! layer thickness
-!$OMP THREADPRIVATE(dz2_SV)
-
-
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   Tsf_SV  ! surface temperature !hj12032010
-!$OMP THREADPRIVATE(Tsf_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   TsfnSV  ! new surface temperature 
-!$OMP THREADPRIVATE(TsfnSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   AcoHSV, BcoHSV ! coefficients for Enthalpy
-!$OMP THREADPRIVATE(AcoHSV, BcoHSV)                             ! evolution,from atmosphere
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   AcoQSV, BcoQSV ! coefficients for Humidity
-!$OMP THREADPRIVATE( AcoQSV, BcoQSV)                            ! evolution,from atmosphere
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   ps__SV  ! surface pressure 
-!$OMP THREADPRIVATE(ps__SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   p1l_SV  ! 1st layer pressure 
-!$OMP THREADPRIVATE(p1l_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   cdH_SV  ! drag coeff Energy (?)
-!$OMP THREADPRIVATE(cdH_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   cdM_SV  ! drag coeff momentum (?)
-!$OMP THREADPRIVATE(cdM_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   rsolSV  ! Radiation balance surface
-!$OMP THREADPRIVATE(rsolSV)
-      REAL,SAVE                            ::   lambSV  ! coefficient soil layers
-!$OMP THREADPRIVATE(lambSV)
-
- 
-CONTAINS
-
-  SUBROUTINE INIT_VARtSV
-  
-  IMPLICIT NONE
-
-  INTEGER ikl
-
-
-
-
-
-
-
-      ALLOCATE(toicSV(klonv))
-
-      ALLOCATE(dz1_SV(klonv,-nsol:nsno))     ! "inverse" layer thicknes
-      ALLOCATE(dz2_SV(klonv,-nsol:nsno))     ! layer thickness
-      ALLOCATE(Tsf_SV(klonv))                ! surface temperature !hj12032010
-      ALLOCATE(TsfnSV(klonv))                ! new surface temperature 
-      ALLOCATE(AcoHSV(klonv), BcoHSV(klonv)) ! coefficients for Enthalpy
-                                             ! evolution,from atmosphere
-      ALLOCATE(AcoQSV(klonv), BcoQSV(klonv)) ! coefficients for Humidity
-                                             ! evolution,from atmosphere                                                                                   
-      ALLOCATE(ps__SV(klonv))                ! surface pressure 
-      ALLOCATE(p1l_SV(klonv))                ! 1st layer pressure 
-      ALLOCATE(cdH_SV(klonv))                ! drag coeff Energy (?)
-      ALLOCATE(cdM_SV(klonv))                ! drag coeff momentum (?)
-      ALLOCATE(rsolSV(klonv))                ! Radiation balance surface
-
-      DO ikl=1,klonv        
-           
-         toicSV(ikl)   = 0.
-         dz1_SV(ikl,:) = 0.
-         dz2_SV(ikl,:) = 0.
-         Tsf_SV(ikl)   = 0.
-         TsfnSV(ikl)   = 0.
-         AcoHSV(ikl)   = 0.
-         BcoHSV(ikl)   = 0.
-         AcoQSV(ikl)   = 0.
-         ps__SV(ikl)   = 0.
-         p1l_SV(ikl)   = 0.
-         cdH_SV(ikl)   = 0.
-         cdM_SV(ikl)   = 0.
-         rsolSV(ikl)   = 0.
-      END DO
-
-
-
-  END SUBROUTINE INIT_VARtSV
-
-END MODULE VARtSV
-
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VARtSV.f90	(revision 5268)
@@ -0,0 +1,92 @@
+MODULE VARtSV
+
+USE VAR_SV, only : nsol, nsno, klonv
+
+IMPLICIT NONE
+! +
+      INTEGER, PARAMETER      ::   ntaver = 4
+! +
+      REAL,ALLOCATABLE,SAVE                ::   toicSV(:)    ! Snow to ice balance    
+!$OMP THREADPRIVATE(toicSV)
+
+      REAL,ALLOCATABLE,SAVE                ::   dz1_SV(:,:) ! "inverse" layer thicknes
+!$OMP THREADPRIVATE(dz1_SV)
+      REAL,ALLOCATABLE,SAVE                ::   dz2_SV(:,:) ! layer thickness
+!$OMP THREADPRIVATE(dz2_SV)
+
+
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   Tsf_SV  ! surface temperature !hj12032010
+!$OMP THREADPRIVATE(Tsf_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   TsfnSV  ! new surface temperature 
+!$OMP THREADPRIVATE(TsfnSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   AcoHSV, BcoHSV ! coefficients for Enthalpy
+!$OMP THREADPRIVATE(AcoHSV, BcoHSV)                             ! evolution,from atmosphere
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   AcoQSV, BcoQSV ! coefficients for Humidity
+!$OMP THREADPRIVATE( AcoQSV, BcoQSV)                            ! evolution,from atmosphere
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   ps__SV  ! surface pressure 
+!$OMP THREADPRIVATE(ps__SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   p1l_SV  ! 1st layer pressure 
+!$OMP THREADPRIVATE(p1l_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   cdH_SV  ! drag coeff Energy (?)
+!$OMP THREADPRIVATE(cdH_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   cdM_SV  ! drag coeff momentum (?)
+!$OMP THREADPRIVATE(cdM_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE  ::   rsolSV  ! Radiation balance surface
+!$OMP THREADPRIVATE(rsolSV)
+      REAL,SAVE                            ::   lambSV  ! coefficient soil layers
+!$OMP THREADPRIVATE(lambSV)
+
+ 
+CONTAINS
+
+  SUBROUTINE INIT_VARtSV
+  
+  IMPLICIT NONE
+
+  INTEGER ikl
+
+
+
+
+
+
+
+      ALLOCATE(toicSV(klonv))
+
+      ALLOCATE(dz1_SV(klonv,-nsol:nsno))     ! "inverse" layer thicknes
+      ALLOCATE(dz2_SV(klonv,-nsol:nsno))     ! layer thickness
+      ALLOCATE(Tsf_SV(klonv))                ! surface temperature !hj12032010
+      ALLOCATE(TsfnSV(klonv))                ! new surface temperature 
+      ALLOCATE(AcoHSV(klonv), BcoHSV(klonv)) ! coefficients for Enthalpy
+                                             ! evolution,from atmosphere
+      ALLOCATE(AcoQSV(klonv), BcoQSV(klonv)) ! coefficients for Humidity
+                                             ! evolution,from atmosphere                                                                                   
+      ALLOCATE(ps__SV(klonv))                ! surface pressure 
+      ALLOCATE(p1l_SV(klonv))                ! 1st layer pressure 
+      ALLOCATE(cdH_SV(klonv))                ! drag coeff Energy (?)
+      ALLOCATE(cdM_SV(klonv))                ! drag coeff momentum (?)
+      ALLOCATE(rsolSV(klonv))                ! Radiation balance surface
+
+      DO ikl=1,klonv        
+           
+         toicSV(ikl)   = 0.
+         dz1_SV(ikl,:) = 0.
+         dz2_SV(ikl,:) = 0.
+         Tsf_SV(ikl)   = 0.
+         TsfnSV(ikl)   = 0.
+         AcoHSV(ikl)   = 0.
+         BcoHSV(ikl)   = 0.
+         AcoQSV(ikl)   = 0.
+         ps__SV(ikl)   = 0.
+         p1l_SV(ikl)   = 0.
+         cdH_SV(ikl)   = 0.
+         cdM_SV(ikl)   = 0.
+         rsolSV(ikl)   = 0.
+      END DO
+
+
+
+  END SUBROUTINE INIT_VARtSV
+
+END MODULE VARtSV
+
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.F90	(revision 5259)
+++ 	(revision )
@@ -1,417 +1,0 @@
-MODULE VARxSV                      
-
-USE VAR_SV, only : klonv, nsol, nsno, nb_wri
-                                                          
-IMPLICIT NONE                                                                                
-! +--SISVAT INPUT        Variables                                              
-! +  -----------------------------                                              
-  
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   LSmask  ! Land-Sea   Mask                 
-!$OMP THREADPRIVATE(LSmask)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   isotSV  ! Soil       Type                 
-!$OMP THREADPRIVATE(isotSV)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   iWaFSV  ! Soil       Drainage:(1,0)=(y,n) 
-!$OMP THREADPRIVATE(iWaFSV)                                                                              
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   coszSV  ! Cosine of Sun zenithal Angle    
-!$OMP THREADPRIVATE(coszSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   sol_SV  ! Downward  Solar    Radiation    
-!$OMP THREADPRIVATE(sol_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRd_SV  ! Downward  Longwave Radiation    
-!$OMP THREADPRIVATE(IRd_SV)
-                                                                                
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   drr_SV  ! Rain  Intensity       [kg/m2/s] 
-!$OMP THREADPRIVATE(drr_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dsn_SV  ! Snow  Intensity       [kg/m2/s] 
-!$OMP THREADPRIVATE(dsn_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dsnbSV  ! Idem, fraction, from Drift  [-] 
-!$OMP THREADPRIVATE(dsnbSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   esnbSV  ! Idem, fraction, from Drift  [-]
-!$OMP THREADPRIVATE(esnbSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_SV  ! Drift Amount            [kg/m2] 
-!$OMP THREADPRIVATE(dbs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_Er  ! BS Erosion              [kg/m2]
-!$OMP THREADPRIVATE(dbs_Er)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_Ac  !  BS Accumulation         [kg/m2]
-!$OMP THREADPRIVATE(dbs_Ac)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BrosSV  ! Buffer Snow Layer Density
-!$OMP THREADPRIVATE(BrosSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BG1sSV  ! Buffer Snow Layer Dendr/Sphe[-]
-!$OMP THREADPRIVATE(BG1sSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BG2sSV  ! Buffer Snow Layer Spher/Size[-][0.0001m]
-!$OMP THREADPRIVATE(BG2sSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dz0_SV  ! dz0(Sastrugi dh)            [m] 
-!$OMP THREADPRIVATE(dz0_SV)
-                                                                               
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   cld_SV  ! Cloudiness (seen from SBL)      
-!$OMP THREADPRIVATE(cld_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   za__SV  ! SBL Height                      
-!$OMP THREADPRIVATE(za__SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   VV__SV  !(SBL Top)  Wind Velocity         
-!$OMP THREADPRIVATE(VV__SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   VVs_SV  !(Sastr,V)  Relevance
-!$OMP THREADPRIVATE(VVs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   RRs_SV  !(Sastr,V)  Counter
-!$OMP THREADPRIVATE(RRs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   DDs_SV  !(Sastr,V)  Angle
-!$OMP THREADPRIVATE(DDs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   TaT_SV  ! SBL Top   Temperature           
-!$OMP THREADPRIVATE(TaT_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   ExnrSV  ! Exner     Potential             
-!$OMP THREADPRIVATE(ExnrSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dSdTSV  ! Sensible Heat Flux T Derivat.   
-!$OMP THREADPRIVATE(dSdTSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dLdTSV  ! Latent   Heat Flux T Derivat.   
-!$OMP THREADPRIVATE(dLdTSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rhT_SV  ! SBL Top   Air  Density          
-!$OMP THREADPRIVATE(rhT_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QaT_SV  ! SBL Top   Specific Humidity     
-!$OMP THREADPRIVATE(QaT_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QsT_SV  ! SBL Top   Specific Humidity
-!$OMP THREADPRIVATE(QsT_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dQa_SV  ! SBL Flux  Limitation of Qa      
-!$OMP THREADPRIVATE(dQa_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   qsnoSV  ! SBL Mean  Snow       Content    
-!$OMP THREADPRIVATE(qsnoSV)                                                                   
-                                                                               
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   alb0SV  ! Soil    Albedo                  
-!$OMP THREADPRIVATE(alb0SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   slopSV  ! Snow/Ice/Soil-Water Surf. Slope 
-!$OMP THREADPRIVATE(slopSV)
-
-                                                              
-      REAL,SAVE                      ::   dt__SV  ! Time Step                       
-!$OMP THREADPRIVATE(dt__SV)
-      CHARACTER (len=18),SAVE        ::   daHost  ! Date Host Model                 
-!$OMP THREADPRIVATE(daHost)
-                                                                               
-                                                                               
-! +--SISVAT INPUT/OUTPUT Variables                                             
-! +  -----------------------------                                             
-                                                                               
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   isnoSV  ! Nb of Ice/Snow Layers           
-!$OMP THREADPRIVATE(isnoSV)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   ispiSV  ! Uppermost superimposed ice      
-!$OMP THREADPRIVATE(ispiSV)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   iiceSV  ! Nb of Ice      Layers           
-!$OMP THREADPRIVATE(iiceSV)
-      INTEGER ,ALLOCATABLE,SAVE              ::   istoSV(:,:)  ! Snow Layer     History          
-!$OMP THREADPRIVATE(istoSV)                                                                              
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   alb_SV  ! Surface Albedo           
-!$OMP THREADPRIVATE(alb_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   emi_SV  ! Surface Emissivity       
-!$OMP THREADPRIVATE(emi_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRs_SV  ! Soil           IR Flux          
-!$OMP THREADPRIVATE(IRs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   LMO_SV  ! Monin-Obukhov  Scale            
-!$OMP THREADPRIVATE(LMO_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   us__SV  ! Friction       Velocity         
-!$OMP THREADPRIVATE(us__SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uts_SV  ! Temperature Turbulent Scale     
-!$OMP THREADPRIVATE(uts_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   cutsSV  ! Temperature Turbulent Scale C.  
-!$OMP THREADPRIVATE(cutsSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uqs_SV  ! Spec.Humid. Turbulent Scale     
-!$OMP THREADPRIVATE(uqs_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uss_SV  ! Blow.Snow   Turbulent Scale     
-!$OMP THREADPRIVATE(uss_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   usthSV  ! Blowing Snow Erosion Thresh.    
-!$OMP THREADPRIVATE(usthSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rCDmSV  ! Square  Root Contribut. Drag_m  
-!$OMP THREADPRIVATE(rCDmSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rCDhSV  ! Square  Root Contribut. Drag_h  
-!$OMP THREADPRIVATE(rCDhSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0m_SV  ! Momentum     Roughness Length   
-!$OMP THREADPRIVATE(Z0m_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0mmSV  !  z0(Momentum,    Time Mean) [m] 
-!$OMP THREADPRIVATE(Z0mmSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0mnSV  !  z0(Momentum,    instanta.) [m] 
-!$OMP THREADPRIVATE(Z0mnSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0roSV  ! Subgrid Topo Roughness Length   
-!$OMP THREADPRIVATE(Z0roSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0SaSV  !  z0(Sastrugi  h)            [m] 
-!$OMP THREADPRIVATE(Z0SaSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0e_SV  !  z0(Snow eroded)            [m] 
-!$OMP THREADPRIVATE(Z0e_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0emSV  !  z0(Snow eroded, Time Mean) [m] 
-!$OMP THREADPRIVATE(Z0emSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0enSV  !  z0(Snow eroded, instanta.) [m] 
-!$OMP THREADPRIVATE(Z0enSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0h_SV  ! Heat         Roughness Length   
-!$OMP THREADPRIVATE(Z0h_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0hmSV  !  z0(Heat,        Time Mean) [m] 
-!$OMP THREADPRIVATE(Z0hmSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0hnSV  !  z0(Heat,        instanta.) [m] 
-!$OMP THREADPRIVATE(Z0hnSV)
-                                                                         
-
-                                                                        
-      REAL ,ALLOCATABLE,SAVE   ::   TsisSV(:,:)  ! Snow/Ice/Soil-Water Temperature 
-!$OMP THREADPRIVATE(TsisSV)
-      REAL ,ALLOCATABLE,SAVE   ::   ro__SV(:,:)  ! Snow/Ice/Soil-Water VolumicMass 
-!$OMP THREADPRIVATE(ro__SV)
-      REAL,ALLOCATABLE,SAVE    ::   eta_SV(:,:)  ! Snow/Ice/Soil     Water Content 
-!$OMP THREADPRIVATE(eta_SV)
-      REAL,ALLOCATABLE,SAVE    ::   G1snSV(:,:)  ! Snow Dendricity/Sphericity      
-!$OMP THREADPRIVATE(G1snSV)
-      REAL,ALLOCATABLE,SAVE    ::   G2snSV(:,:)  ! Snow Sphericity/Size            
-!$OMP THREADPRIVATE(G2snSV)
-      REAL,ALLOCATABLE,SAVE    ::   dzsnSV(:,:)  ! Snow Layer  Thickness           
-!$OMP THREADPRIVATE(dzsnSV)
-      REAL,ALLOCATABLE,SAVE    ::   agsnSV(:,:)  ! Snow Age                        
-!$OMP THREADPRIVATE(agsnSV)
-      REAL,ALLOCATABLE,SAVE    ::   DOPsnSV(:,:)  ! Snow optical diameter [m]
-!$OMP THREADPRIVATE(DOPsnSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BufsSV  ! Snow Buffer Layer               
-!$OMP THREADPRIVATE(BufsSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rusnSV  ! Surficial   Water               
-!$OMP THREADPRIVATE(rusnSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   SWf_SV  ! Normalized  Decay               
-!$OMP THREADPRIVATE(SWf_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   SWS_SV  ! Surficial Water Status          
-!$OMP THREADPRIVATE(SWS_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   HFraSV  ! Frazil      Thickness           
-!$OMP THREADPRIVATE(HFraSV)
-                                                                            
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zWE_SV  ! Current   Snow Thickness [mmWE] 
-!$OMP THREADPRIVATE(zWE_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zWEcSV  ! Compacted Snow Thickness [mmWE] 
-!$OMP THREADPRIVATE(zWEcSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wem_SV  ! Only Melting             [mmWE] 
-!$OMP THREADPRIVATE(wem_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wer_SV  ! Refreezing               [mmWE] 
-!$OMP THREADPRIVATE(wer_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wes_SV  ! Sublimation              [mmWE] 
-!$OMP THREADPRIVATE(wes_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zn4_SV  ! snowheight change         [mm]
-!$OMP THREADPRIVATE(zn4_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zn5_SV  ! snowheight change         [mm]
-!$OMP THREADPRIVATE(zn5_SV)                                                                            
-                                                                              
-! +--SISVAT OUTPUT       Variables                                            
-! +  -----------------------------                                            
-                                                                               
-      INTEGER,DIMENSION(nb_wri),SAVE ::   no__SV  ! OUTPUT file Unit Number         
-!$OMP THREADPRIVATE(no__SV)
-      INTEGER,DIMENSION(nb_wri),SAVE ::   i___SV  ! OUTPUT point   i Coordinate    
-!$OMP THREADPRIVATE(i___SV)
-      INTEGER,DIMENSION(nb_wri),SAVE ::   j___SV  ! OUTPUT point   j Coordinate     
-!$OMP THREADPRIVATE(j___SV)
-      INTEGER,DIMENSION(nb_wri),SAVE ::   n___SV  ! OUTPUT point   n Coordinate     
-!$OMP THREADPRIVATE(n___SV)
-      INTEGER,DIMENSION(nb_wri),SAVE ::   lwriSV  ! OUTPUT point vec Index          
-!$OMP THREADPRIVATE(lwriSV)
-!                                                                              
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   ii__SV  ! WORK   point   i Coordinate     
-!$OMP THREADPRIVATE(ii__SV)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   jj__SV  ! WORK   point   j Coordinate     
-!$OMP THREADPRIVATE(jj__SV)
-      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   nn__SV  ! WORK   point   n Coordinate     
-!$OMP THREADPRIVATE(nn__SV)
-                                                                               
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRu_SV  ! UPward    IR Flux (effective)   
-!$OMP THREADPRIVATE(IRu_SV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   hSalSV  ! Saltating Layer Height          
-!$OMP THREADPRIVATE(hSalSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   qSalSV  ! Saltating Snow  Concentration   
-!$OMP THREADPRIVATE(qSalSV)
-      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   RnofSV  ! RunOFF    Intensity             
-!$OMP THREADPRIVATE(RnofSV)
-      REAL, DIMENSION(:,:),ALLOCATABLE,SAVE    ::   RuofSV  ! RunOFF    Intensity             
-!$OMP THREADPRIVATE(RuofSV)
-
-
-
-                  
-CONTAINS
-
-
-
-  SUBROUTINE INIT_VARxSV
-  IMPLICIT NONE
-
-      INTEGER                                :: ikl
-
-      ALLOCATE(LSmask(klonv))  ! Land-Sea   Mask                 
-      ALLOCATE(isotSV(klonv))  ! Soil       Type                 
-      ALLOCATE(iWaFSV(klonv))  ! Soil       Drainage:(1,0)=(y,n)                                                                               
-      ALLOCATE(coszSV(klonv))  ! Cosine of Sun zenithal Angle    
-      ALLOCATE(sol_SV(klonv))  ! Downward  Solar    Radiation    
-      ALLOCATE(IRd_SV(klonv))  ! Downward  Longwave Radiation    
-                                                                                
-      ALLOCATE(drr_SV(klonv))  ! Rain  Intensity       [kg/m2/s] 
-      ALLOCATE(dsn_SV(klonv))  ! Snow  Intensity       [kg/m2/s] 
-      ALLOCATE(dsnbSV(klonv))  ! Idem, from Drift        [kg/m2] 
-      ALLOCATE(esnbSV(klonv))  ! Idem, from Drift        [kg/m2] 
-      ALLOCATE(dbs_SV(klonv))  ! Drift Amount            [kg/m2] 
-      ALLOCATE(dbs_er(klonv))  ! Drift erosion           [kg/m2] 
-      ALLOCATE(dbs_ac(klonv))  ! Drift accumulation      [kg/m2] 
-      ALLOCATE(BrosSV(klonv))
-      ALLOCATE(BG1sSV(klonv))          
-      ALLOCATE(BG2sSV(klonv))
-      ALLOCATE(dz0_SV(klonv))  ! dz0(Sastrugi dh)            [m] 
-                                                                               
-      ALLOCATE(cld_SV(klonv))  ! Cloudiness (seen from SBL)      
-      ALLOCATE(za__SV(klonv))  ! SBL Height                      
-      ALLOCATE(VV__SV(klonv))  !(SBL Top)  Wind Velocity         
-      ALLOCATE(VVs_SV(klonv))
-      ALLOCATE(RRs_SV(klonv))
-      ALLOCATE(DDs_SV(klonv))
-      ALLOCATE(TaT_SV(klonv))  ! SBL Top   Temperature           
-      ALLOCATE(ExnrSV(klonv))  ! Exner     Potential             
-      ALLOCATE(dSdTSV(klonv))  ! Sensible Heat Flux T Derivat.   
-      ALLOCATE(dLdTSV(klonv))  ! Latent   Heat Flux T Derivat.   
-      ALLOCATE(rhT_SV(klonv))  ! SBL Top   Air  Density          
-      ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity    
-      ALLOCATE(QsT_SV(klonv))  ! surface   Specific Humidity 
-      ALLOCATE(dQa_SV(klonv))  ! SBL Flux  Limitation of Qa      
-      ALLOCATE(qsnoSV(klonv))  ! SBL Mean  Snow       Content    
-                                                                               
-      ALLOCATE(alb0SV(klonv))  ! Soil    Albedo                  
-      ALLOCATE(slopSV(klonv))  ! Snow/Ice/Soil-Water Surf. Slope 
-
-                                                                               
-                                                                               
-! +--SISVAT INPUT/OUTPUT Variables                                             
-! +  -----------------------------                                             
-                                                                               
-      ALLOCATE(isnoSV(klonv))  ! Nb of Ice/Snow Layers           
-      ALLOCATE(ispiSV(klonv))  ! Uppermost superimposed ice      
-      ALLOCATE(iiceSV(klonv))  ! Nb of Ice      Layers           
-      ALLOCATE(istoSV(klonv,0:nsno))  ! Snow Layer     History          
-                                                                               
-      ALLOCATE(alb_SV(klonv))  ! Surface-Canopy Albedo           
-      ALLOCATE(emi_SV(klonv))  ! Surface-Canopy Emissivity       
-      ALLOCATE(IRs_SV(klonv))  ! Soil           IR Flux          
-      ALLOCATE(LMO_SV(klonv))  ! Monin-Obukhov  Scale            
-      ALLOCATE(us__SV(klonv))  ! Friction       Velocity         
-      ALLOCATE(uts_SV(klonv))  ! Temperature Turbulent Scale     
-      ALLOCATE(cutsSV(klonv))  ! Temperature Turbulent Scale C.  
-      ALLOCATE(uqs_SV(klonv))  ! Spec.Humid. Turbulent Scale     
-      ALLOCATE(uss_SV(klonv))  ! Blow.Snow   Turbulent Scale     
-      ALLOCATE(usthSV(klonv))  ! Blowing Snow Erosion Thresh.    
-      ALLOCATE(rCDmSV(klonv))  ! Square  Root Contribut. Drag_m  
-      ALLOCATE(rCDhSV(klonv))  ! Square  Root Contribut. Drag_h  
-      ALLOCATE(Z0m_SV(klonv))  ! Momentum     Roughness Length   
-      ALLOCATE(Z0mmSV(klonv))  !  z0(Momentum,    Time Mean) [m] 
-      ALLOCATE(Z0mnSV(klonv))  !  z0(Momentum,    instanta.) [m] 
-      ALLOCATE(Z0roSV(klonv))  ! Subgrid Topo Roughness Length   
-      ALLOCATE(Z0SaSV(klonv))  !  z0(Sastrugi  h)            [m] 
-      ALLOCATE(Z0e_SV(klonv))  !  z0(Snow eroded)            [m] 
-      ALLOCATE(Z0emSV(klonv))  !  z0(Snow eroded, Time Mean) [m] 
-      ALLOCATE(Z0enSV(klonv))  !  z0(Snow eroded, instanta.) [m] 
-      ALLOCATE(Z0h_SV(klonv))  !  Heat         Roughness Length   
-      ALLOCATE(Z0hmSV(klonv))  !  z0(Heat,        Time Mean) [m] 
-      ALLOCATE(Z0hnSV(klonv))  !  z0(Heat,        instanta.) [m] 
-                                                                          
-                                                                        
-      ALLOCATE(TsisSV(klonv,-nsol:nsno))  ! Snow/Ice/Soil-Water Temperature 
-      ALLOCATE(ro__SV(klonv,-nsol:nsno))  ! Snow/Ice/Soil-Water VolumicMass 
-      ALLOCATE(eta_SV(klonv,-nsol:nsno))  ! Snow/Ice/Soil     Water Content 
-      ALLOCATE(G1snSV(klonv,    0:nsno))  ! Snow Dendricity/Sphericity      
-      ALLOCATE(G2snSV(klonv,    0:nsno))  ! Snow Sphericity/Size            
-      ALLOCATE(dzsnSV(klonv,    0:nsno))  ! Snow Layer  Thickness           
-      ALLOCATE(agsnSV(klonv,    0:nsno))  ! Snow Age                        
-      ALLOCATE(DOPsnSV(klonv,    0:nsno))  ! Snow Optical diameter                        
-      ALLOCATE(BufsSV(klonv))  ! Snow Buffer Layer               
-      ALLOCATE(rusnSV(klonv))  ! Surficial   Water               
-      ALLOCATE(SWf_SV(klonv))  ! Normalized  Decay               
-      ALLOCATE(SWS_SV(klonv))  ! Surficial Water Status          
-      ALLOCATE(HFraSV(klonv))  ! Frazil      Thickness           
-                                                                            
-      ALLOCATE(zWE_SV(klonv))  ! Current   Snow Thickness [mmWE] 
-      ALLOCATE(zWEcSV(klonv))  ! Compacted Snow Thickness [mmWE] 
-      ALLOCATE(wem_SV(klonv)) ! Only Melting             [mmWE] 
-      ALLOCATE(wer_SV(klonv)) ! Refreezing               [mmWE] 
-      ALLOCATE(wes_SV(klonv))  ! Sublimation              [mmWE] 
-      ALLOCATE(zn4_SV(klonv))  ! snow height change       [mm] 
-      ALLOCATE(zn5_SV(klonv))  ! snow height change       [mm]                                                                             
-                                                                              
-! +--SISVAT OUTPUT       Variables                                            
-! +  -----------------------------                                            
-                                                                               
-!                                                                              
-      ALLOCATE(ii__SV(klonv))  ! WORK   point   i Coordinate     
-      ALLOCATE(jj__SV(klonv))  ! WORK   point   j Coordinate     
-      ALLOCATE(nn__SV(klonv))  ! WORK   point   n Coordinate     
-                                                                               
-      ALLOCATE(IRu_SV(klonv))  ! UPward    IR Flux (effective)   
-      ALLOCATE(hSalSV(klonv))  ! Saltating Layer Height          
-      ALLOCATE(qSalSV(klonv))  ! Saltating Snow  Concentration   
-      ALLOCATE(RnofSV(klonv))  ! RunOFF    Intensity      
-      ALLOCATE(RuofSV(klonv,4))  ! RunOFF    Intensity      
-
-
-      DO ikl=1,klonv        
-
-
-      isnoSV(ikl)  =0.        
-      ispiSV(ikl)  =0.   
-      iiceSV(ikl)  =0.        
-      istoSV(ikl,:)=0.                                                                        
-      alb_SV(ikl)  =0.      
-      emi_SV(ikl)  =0.    
-      IRs_SV(ikl)  =0.      
-      LMO_SV(ikl)  =0.
-      us__SV(ikl)  =0.
-      uts_SV(ikl)  =0.
-      cutsSV(ikl)  =0.
-      uqs_SV(ikl)  =0.
-      uss_SV(ikl)  =0.
-      usthSV(ikl)  =0.
-      rCDmSV(ikl)  =0.
-      rCDhSV(ikl)  =0.
-      Z0m_SV(ikl)  =0.
-      Z0mmSV(ikl)  =0.
-      Z0mnSV(ikl)  =0.
-      Z0roSV(ikl)  =0.
-      Z0SaSV(ikl)  =0.
-      Z0e_SV(ikl)  =0.
-      Z0emSV(ikl)  =0.
-      Z0enSV(ikl)  =0.
-      Z0h_SV(ikl)  =0.
-      Z0hmSV(ikl)  =0.
-      Z0hnSV(ikl)  =0.
-                                                                          
-                                                                        
-      TsisSV(ikl,:)  =0.
-      ro__SV(ikl,:)  =0.
-      eta_SV(ikl,:)  =0.
-      G1snSV(ikl,:)  =0.  
-      G2snSV(ikl,:)  =0.          
-      dzsnSV(ikl,:)  =0.      
-      agsnSV(ikl,:)  =0.              
-      DOPsnSV(ikl,:) =0.                  
-      BufsSV(ikl)  =0.         
-      rusnSV(ikl)  =0.     
-      SWf_SV(ikl)  =0.        
-      SWS_SV(ikl)  =0.
-      HFraSV(ikl)  =0.       
-                                                                            
-      zWE_SV(ikl)  =0.
-      zWEcSV(ikl)  =0.
-      wem_SV(ikl)  =0.
-      wer_SV(ikl)  =0.
-      wes_SV(ikl)  =0.
-      zn4_SV(ikl)  =0.
-      zn5_SV(ikl)  =0.                                                        
-                                       
-                                                                                                                                                          
-      ii__SV(ikl)  =0. 
-      jj__SV(ikl)  =0.
-      nn__SV(ikl)  =0.
-                                                                               
-      IRu_SV(ikl)  =0. 
-      hSalSV(ikl)  =0.     
-      qSalSV(ikl)  =0.
-      RnofSV(ikl)  =0.    
-      RuofSV(ikl,:)  =0.
-
-
-
-
-      END DO
-  END SUBROUTINE INIT_VARxSV
-
-
-                                               
-END MODULE VARxSV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VARxSV.f90	(revision 5268)
@@ -0,0 +1,417 @@
+MODULE VARxSV                      
+
+USE VAR_SV, only : klonv, nsol, nsno, nb_wri
+                                                          
+IMPLICIT NONE                                                                                
+! +--SISVAT INPUT        Variables                                              
+! +  -----------------------------                                              
+  
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   LSmask  ! Land-Sea   Mask                 
+!$OMP THREADPRIVATE(LSmask)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   isotSV  ! Soil       Type                 
+!$OMP THREADPRIVATE(isotSV)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   iWaFSV  ! Soil       Drainage:(1,0)=(y,n) 
+!$OMP THREADPRIVATE(iWaFSV)                                                                              
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   coszSV  ! Cosine of Sun zenithal Angle    
+!$OMP THREADPRIVATE(coszSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   sol_SV  ! Downward  Solar    Radiation    
+!$OMP THREADPRIVATE(sol_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRd_SV  ! Downward  Longwave Radiation    
+!$OMP THREADPRIVATE(IRd_SV)
+                                                                                
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   drr_SV  ! Rain  Intensity       [kg/m2/s] 
+!$OMP THREADPRIVATE(drr_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dsn_SV  ! Snow  Intensity       [kg/m2/s] 
+!$OMP THREADPRIVATE(dsn_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dsnbSV  ! Idem, fraction, from Drift  [-] 
+!$OMP THREADPRIVATE(dsnbSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   esnbSV  ! Idem, fraction, from Drift  [-]
+!$OMP THREADPRIVATE(esnbSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_SV  ! Drift Amount            [kg/m2] 
+!$OMP THREADPRIVATE(dbs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_Er  ! BS Erosion              [kg/m2]
+!$OMP THREADPRIVATE(dbs_Er)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dbs_Ac  !  BS Accumulation         [kg/m2]
+!$OMP THREADPRIVATE(dbs_Ac)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BrosSV  ! Buffer Snow Layer Density
+!$OMP THREADPRIVATE(BrosSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BG1sSV  ! Buffer Snow Layer Dendr/Sphe[-]
+!$OMP THREADPRIVATE(BG1sSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BG2sSV  ! Buffer Snow Layer Spher/Size[-][0.0001m]
+!$OMP THREADPRIVATE(BG2sSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dz0_SV  ! dz0(Sastrugi dh)            [m] 
+!$OMP THREADPRIVATE(dz0_SV)
+                                                                               
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   cld_SV  ! Cloudiness (seen from SBL)      
+!$OMP THREADPRIVATE(cld_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   za__SV  ! SBL Height                      
+!$OMP THREADPRIVATE(za__SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   VV__SV  !(SBL Top)  Wind Velocity         
+!$OMP THREADPRIVATE(VV__SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   VVs_SV  !(Sastr,V)  Relevance
+!$OMP THREADPRIVATE(VVs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   RRs_SV  !(Sastr,V)  Counter
+!$OMP THREADPRIVATE(RRs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   DDs_SV  !(Sastr,V)  Angle
+!$OMP THREADPRIVATE(DDs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   TaT_SV  ! SBL Top   Temperature           
+!$OMP THREADPRIVATE(TaT_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   ExnrSV  ! Exner     Potential             
+!$OMP THREADPRIVATE(ExnrSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dSdTSV  ! Sensible Heat Flux T Derivat.   
+!$OMP THREADPRIVATE(dSdTSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dLdTSV  ! Latent   Heat Flux T Derivat.   
+!$OMP THREADPRIVATE(dLdTSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rhT_SV  ! SBL Top   Air  Density          
+!$OMP THREADPRIVATE(rhT_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QaT_SV  ! SBL Top   Specific Humidity     
+!$OMP THREADPRIVATE(QaT_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   QsT_SV  ! SBL Top   Specific Humidity
+!$OMP THREADPRIVATE(QsT_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   dQa_SV  ! SBL Flux  Limitation of Qa      
+!$OMP THREADPRIVATE(dQa_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   qsnoSV  ! SBL Mean  Snow       Content    
+!$OMP THREADPRIVATE(qsnoSV)                                                                   
+                                                                               
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   alb0SV  ! Soil    Albedo                  
+!$OMP THREADPRIVATE(alb0SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   slopSV  ! Snow/Ice/Soil-Water Surf. Slope 
+!$OMP THREADPRIVATE(slopSV)
+
+                                                              
+      REAL,SAVE                      ::   dt__SV  ! Time Step                       
+!$OMP THREADPRIVATE(dt__SV)
+      CHARACTER (len=18),SAVE        ::   daHost  ! Date Host Model                 
+!$OMP THREADPRIVATE(daHost)
+                                                                               
+                                                                               
+! +--SISVAT INPUT/OUTPUT Variables                                             
+! +  -----------------------------                                             
+                                                                               
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   isnoSV  ! Nb of Ice/Snow Layers           
+!$OMP THREADPRIVATE(isnoSV)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   ispiSV  ! Uppermost superimposed ice      
+!$OMP THREADPRIVATE(ispiSV)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   iiceSV  ! Nb of Ice      Layers           
+!$OMP THREADPRIVATE(iiceSV)
+      INTEGER ,ALLOCATABLE,SAVE              ::   istoSV(:,:)  ! Snow Layer     History          
+!$OMP THREADPRIVATE(istoSV)                                                                              
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   alb_SV  ! Surface Albedo           
+!$OMP THREADPRIVATE(alb_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   emi_SV  ! Surface Emissivity       
+!$OMP THREADPRIVATE(emi_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRs_SV  ! Soil           IR Flux          
+!$OMP THREADPRIVATE(IRs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   LMO_SV  ! Monin-Obukhov  Scale            
+!$OMP THREADPRIVATE(LMO_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   us__SV  ! Friction       Velocity         
+!$OMP THREADPRIVATE(us__SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uts_SV  ! Temperature Turbulent Scale     
+!$OMP THREADPRIVATE(uts_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   cutsSV  ! Temperature Turbulent Scale C.  
+!$OMP THREADPRIVATE(cutsSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uqs_SV  ! Spec.Humid. Turbulent Scale     
+!$OMP THREADPRIVATE(uqs_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   uss_SV  ! Blow.Snow   Turbulent Scale     
+!$OMP THREADPRIVATE(uss_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   usthSV  ! Blowing Snow Erosion Thresh.    
+!$OMP THREADPRIVATE(usthSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rCDmSV  ! Square  Root Contribut. Drag_m  
+!$OMP THREADPRIVATE(rCDmSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rCDhSV  ! Square  Root Contribut. Drag_h  
+!$OMP THREADPRIVATE(rCDhSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0m_SV  ! Momentum     Roughness Length   
+!$OMP THREADPRIVATE(Z0m_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0mmSV  !  z0(Momentum,    Time Mean) [m] 
+!$OMP THREADPRIVATE(Z0mmSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0mnSV  !  z0(Momentum,    instanta.) [m] 
+!$OMP THREADPRIVATE(Z0mnSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0roSV  ! Subgrid Topo Roughness Length   
+!$OMP THREADPRIVATE(Z0roSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0SaSV  !  z0(Sastrugi  h)            [m] 
+!$OMP THREADPRIVATE(Z0SaSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0e_SV  !  z0(Snow eroded)            [m] 
+!$OMP THREADPRIVATE(Z0e_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0emSV  !  z0(Snow eroded, Time Mean) [m] 
+!$OMP THREADPRIVATE(Z0emSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0enSV  !  z0(Snow eroded, instanta.) [m] 
+!$OMP THREADPRIVATE(Z0enSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0h_SV  ! Heat         Roughness Length   
+!$OMP THREADPRIVATE(Z0h_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0hmSV  !  z0(Heat,        Time Mean) [m] 
+!$OMP THREADPRIVATE(Z0hmSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   Z0hnSV  !  z0(Heat,        instanta.) [m] 
+!$OMP THREADPRIVATE(Z0hnSV)
+                                                                         
+
+                                                                        
+      REAL ,ALLOCATABLE,SAVE   ::   TsisSV(:,:)  ! Snow/Ice/Soil-Water Temperature 
+!$OMP THREADPRIVATE(TsisSV)
+      REAL ,ALLOCATABLE,SAVE   ::   ro__SV(:,:)  ! Snow/Ice/Soil-Water VolumicMass 
+!$OMP THREADPRIVATE(ro__SV)
+      REAL,ALLOCATABLE,SAVE    ::   eta_SV(:,:)  ! Snow/Ice/Soil     Water Content 
+!$OMP THREADPRIVATE(eta_SV)
+      REAL,ALLOCATABLE,SAVE    ::   G1snSV(:,:)  ! Snow Dendricity/Sphericity      
+!$OMP THREADPRIVATE(G1snSV)
+      REAL,ALLOCATABLE,SAVE    ::   G2snSV(:,:)  ! Snow Sphericity/Size            
+!$OMP THREADPRIVATE(G2snSV)
+      REAL,ALLOCATABLE,SAVE    ::   dzsnSV(:,:)  ! Snow Layer  Thickness           
+!$OMP THREADPRIVATE(dzsnSV)
+      REAL,ALLOCATABLE,SAVE    ::   agsnSV(:,:)  ! Snow Age                        
+!$OMP THREADPRIVATE(agsnSV)
+      REAL,ALLOCATABLE,SAVE    ::   DOPsnSV(:,:)  ! Snow optical diameter [m]
+!$OMP THREADPRIVATE(DOPsnSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   BufsSV  ! Snow Buffer Layer               
+!$OMP THREADPRIVATE(BufsSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   rusnSV  ! Surficial   Water               
+!$OMP THREADPRIVATE(rusnSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   SWf_SV  ! Normalized  Decay               
+!$OMP THREADPRIVATE(SWf_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   SWS_SV  ! Surficial Water Status          
+!$OMP THREADPRIVATE(SWS_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   HFraSV  ! Frazil      Thickness           
+!$OMP THREADPRIVATE(HFraSV)
+                                                                            
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zWE_SV  ! Current   Snow Thickness [mmWE] 
+!$OMP THREADPRIVATE(zWE_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zWEcSV  ! Compacted Snow Thickness [mmWE] 
+!$OMP THREADPRIVATE(zWEcSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wem_SV  ! Only Melting             [mmWE] 
+!$OMP THREADPRIVATE(wem_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wer_SV  ! Refreezing               [mmWE] 
+!$OMP THREADPRIVATE(wer_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   wes_SV  ! Sublimation              [mmWE] 
+!$OMP THREADPRIVATE(wes_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zn4_SV  ! snowheight change         [mm]
+!$OMP THREADPRIVATE(zn4_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   zn5_SV  ! snowheight change         [mm]
+!$OMP THREADPRIVATE(zn5_SV)                                                                            
+                                                                              
+! +--SISVAT OUTPUT       Variables                                            
+! +  -----------------------------                                            
+                                                                               
+      INTEGER,DIMENSION(nb_wri),SAVE ::   no__SV  ! OUTPUT file Unit Number         
+!$OMP THREADPRIVATE(no__SV)
+      INTEGER,DIMENSION(nb_wri),SAVE ::   i___SV  ! OUTPUT point   i Coordinate    
+!$OMP THREADPRIVATE(i___SV)
+      INTEGER,DIMENSION(nb_wri),SAVE ::   j___SV  ! OUTPUT point   j Coordinate     
+!$OMP THREADPRIVATE(j___SV)
+      INTEGER,DIMENSION(nb_wri),SAVE ::   n___SV  ! OUTPUT point   n Coordinate     
+!$OMP THREADPRIVATE(n___SV)
+      INTEGER,DIMENSION(nb_wri),SAVE ::   lwriSV  ! OUTPUT point vec Index          
+!$OMP THREADPRIVATE(lwriSV)
+!                                                                              
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   ii__SV  ! WORK   point   i Coordinate     
+!$OMP THREADPRIVATE(ii__SV)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   jj__SV  ! WORK   point   j Coordinate     
+!$OMP THREADPRIVATE(jj__SV)
+      INTEGER, DIMENSION(:),ALLOCATABLE,SAVE ::   nn__SV  ! WORK   point   n Coordinate     
+!$OMP THREADPRIVATE(nn__SV)
+                                                                               
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   IRu_SV  ! UPward    IR Flux (effective)   
+!$OMP THREADPRIVATE(IRu_SV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   hSalSV  ! Saltating Layer Height          
+!$OMP THREADPRIVATE(hSalSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   qSalSV  ! Saltating Snow  Concentration   
+!$OMP THREADPRIVATE(qSalSV)
+      REAL, DIMENSION(:),ALLOCATABLE,SAVE    ::   RnofSV  ! RunOFF    Intensity             
+!$OMP THREADPRIVATE(RnofSV)
+      REAL, DIMENSION(:,:),ALLOCATABLE,SAVE    ::   RuofSV  ! RunOFF    Intensity             
+!$OMP THREADPRIVATE(RuofSV)
+
+
+
+                  
+CONTAINS
+
+
+
+  SUBROUTINE INIT_VARxSV
+  IMPLICIT NONE
+
+      INTEGER                                :: ikl
+
+      ALLOCATE(LSmask(klonv))  ! Land-Sea   Mask                 
+      ALLOCATE(isotSV(klonv))  ! Soil       Type                 
+      ALLOCATE(iWaFSV(klonv))  ! Soil       Drainage:(1,0)=(y,n)                                                                               
+      ALLOCATE(coszSV(klonv))  ! Cosine of Sun zenithal Angle    
+      ALLOCATE(sol_SV(klonv))  ! Downward  Solar    Radiation    
+      ALLOCATE(IRd_SV(klonv))  ! Downward  Longwave Radiation    
+                                                                                
+      ALLOCATE(drr_SV(klonv))  ! Rain  Intensity       [kg/m2/s] 
+      ALLOCATE(dsn_SV(klonv))  ! Snow  Intensity       [kg/m2/s] 
+      ALLOCATE(dsnbSV(klonv))  ! Idem, from Drift        [kg/m2] 
+      ALLOCATE(esnbSV(klonv))  ! Idem, from Drift        [kg/m2] 
+      ALLOCATE(dbs_SV(klonv))  ! Drift Amount            [kg/m2] 
+      ALLOCATE(dbs_er(klonv))  ! Drift erosion           [kg/m2] 
+      ALLOCATE(dbs_ac(klonv))  ! Drift accumulation      [kg/m2] 
+      ALLOCATE(BrosSV(klonv))
+      ALLOCATE(BG1sSV(klonv))          
+      ALLOCATE(BG2sSV(klonv))
+      ALLOCATE(dz0_SV(klonv))  ! dz0(Sastrugi dh)            [m] 
+                                                                               
+      ALLOCATE(cld_SV(klonv))  ! Cloudiness (seen from SBL)      
+      ALLOCATE(za__SV(klonv))  ! SBL Height                      
+      ALLOCATE(VV__SV(klonv))  !(SBL Top)  Wind Velocity         
+      ALLOCATE(VVs_SV(klonv))
+      ALLOCATE(RRs_SV(klonv))
+      ALLOCATE(DDs_SV(klonv))
+      ALLOCATE(TaT_SV(klonv))  ! SBL Top   Temperature           
+      ALLOCATE(ExnrSV(klonv))  ! Exner     Potential             
+      ALLOCATE(dSdTSV(klonv))  ! Sensible Heat Flux T Derivat.   
+      ALLOCATE(dLdTSV(klonv))  ! Latent   Heat Flux T Derivat.   
+      ALLOCATE(rhT_SV(klonv))  ! SBL Top   Air  Density          
+      ALLOCATE(QaT_SV(klonv))  ! SBL Top   Specific Humidity    
+      ALLOCATE(QsT_SV(klonv))  ! surface   Specific Humidity 
+      ALLOCATE(dQa_SV(klonv))  ! SBL Flux  Limitation of Qa      
+      ALLOCATE(qsnoSV(klonv))  ! SBL Mean  Snow       Content    
+                                                                               
+      ALLOCATE(alb0SV(klonv))  ! Soil    Albedo                  
+      ALLOCATE(slopSV(klonv))  ! Snow/Ice/Soil-Water Surf. Slope 
+
+                                                                               
+                                                                               
+! +--SISVAT INPUT/OUTPUT Variables                                             
+! +  -----------------------------                                             
+                                                                               
+      ALLOCATE(isnoSV(klonv))  ! Nb of Ice/Snow Layers           
+      ALLOCATE(ispiSV(klonv))  ! Uppermost superimposed ice      
+      ALLOCATE(iiceSV(klonv))  ! Nb of Ice      Layers           
+      ALLOCATE(istoSV(klonv,0:nsno))  ! Snow Layer     History          
+                                                                               
+      ALLOCATE(alb_SV(klonv))  ! Surface-Canopy Albedo           
+      ALLOCATE(emi_SV(klonv))  ! Surface-Canopy Emissivity       
+      ALLOCATE(IRs_SV(klonv))  ! Soil           IR Flux          
+      ALLOCATE(LMO_SV(klonv))  ! Monin-Obukhov  Scale            
+      ALLOCATE(us__SV(klonv))  ! Friction       Velocity         
+      ALLOCATE(uts_SV(klonv))  ! Temperature Turbulent Scale     
+      ALLOCATE(cutsSV(klonv))  ! Temperature Turbulent Scale C.  
+      ALLOCATE(uqs_SV(klonv))  ! Spec.Humid. Turbulent Scale     
+      ALLOCATE(uss_SV(klonv))  ! Blow.Snow   Turbulent Scale     
+      ALLOCATE(usthSV(klonv))  ! Blowing Snow Erosion Thresh.    
+      ALLOCATE(rCDmSV(klonv))  ! Square  Root Contribut. Drag_m  
+      ALLOCATE(rCDhSV(klonv))  ! Square  Root Contribut. Drag_h  
+      ALLOCATE(Z0m_SV(klonv))  ! Momentum     Roughness Length   
+      ALLOCATE(Z0mmSV(klonv))  !  z0(Momentum,    Time Mean) [m] 
+      ALLOCATE(Z0mnSV(klonv))  !  z0(Momentum,    instanta.) [m] 
+      ALLOCATE(Z0roSV(klonv))  ! Subgrid Topo Roughness Length   
+      ALLOCATE(Z0SaSV(klonv))  !  z0(Sastrugi  h)            [m] 
+      ALLOCATE(Z0e_SV(klonv))  !  z0(Snow eroded)            [m] 
+      ALLOCATE(Z0emSV(klonv))  !  z0(Snow eroded, Time Mean) [m] 
+      ALLOCATE(Z0enSV(klonv))  !  z0(Snow eroded, instanta.) [m] 
+      ALLOCATE(Z0h_SV(klonv))  !  Heat         Roughness Length   
+      ALLOCATE(Z0hmSV(klonv))  !  z0(Heat,        Time Mean) [m] 
+      ALLOCATE(Z0hnSV(klonv))  !  z0(Heat,        instanta.) [m] 
+                                                                          
+                                                                        
+      ALLOCATE(TsisSV(klonv,-nsol:nsno))  ! Snow/Ice/Soil-Water Temperature 
+      ALLOCATE(ro__SV(klonv,-nsol:nsno))  ! Snow/Ice/Soil-Water VolumicMass 
+      ALLOCATE(eta_SV(klonv,-nsol:nsno))  ! Snow/Ice/Soil     Water Content 
+      ALLOCATE(G1snSV(klonv,    0:nsno))  ! Snow Dendricity/Sphericity      
+      ALLOCATE(G2snSV(klonv,    0:nsno))  ! Snow Sphericity/Size            
+      ALLOCATE(dzsnSV(klonv,    0:nsno))  ! Snow Layer  Thickness           
+      ALLOCATE(agsnSV(klonv,    0:nsno))  ! Snow Age                        
+      ALLOCATE(DOPsnSV(klonv,    0:nsno))  ! Snow Optical diameter                        
+      ALLOCATE(BufsSV(klonv))  ! Snow Buffer Layer               
+      ALLOCATE(rusnSV(klonv))  ! Surficial   Water               
+      ALLOCATE(SWf_SV(klonv))  ! Normalized  Decay               
+      ALLOCATE(SWS_SV(klonv))  ! Surficial Water Status          
+      ALLOCATE(HFraSV(klonv))  ! Frazil      Thickness           
+                                                                            
+      ALLOCATE(zWE_SV(klonv))  ! Current   Snow Thickness [mmWE] 
+      ALLOCATE(zWEcSV(klonv))  ! Compacted Snow Thickness [mmWE] 
+      ALLOCATE(wem_SV(klonv)) ! Only Melting             [mmWE] 
+      ALLOCATE(wer_SV(klonv)) ! Refreezing               [mmWE] 
+      ALLOCATE(wes_SV(klonv))  ! Sublimation              [mmWE] 
+      ALLOCATE(zn4_SV(klonv))  ! snow height change       [mm] 
+      ALLOCATE(zn5_SV(klonv))  ! snow height change       [mm]                                                                             
+                                                                              
+! +--SISVAT OUTPUT       Variables                                            
+! +  -----------------------------                                            
+                                                                               
+!                                                                              
+      ALLOCATE(ii__SV(klonv))  ! WORK   point   i Coordinate     
+      ALLOCATE(jj__SV(klonv))  ! WORK   point   j Coordinate     
+      ALLOCATE(nn__SV(klonv))  ! WORK   point   n Coordinate     
+                                                                               
+      ALLOCATE(IRu_SV(klonv))  ! UPward    IR Flux (effective)   
+      ALLOCATE(hSalSV(klonv))  ! Saltating Layer Height          
+      ALLOCATE(qSalSV(klonv))  ! Saltating Snow  Concentration   
+      ALLOCATE(RnofSV(klonv))  ! RunOFF    Intensity      
+      ALLOCATE(RuofSV(klonv,4))  ! RunOFF    Intensity      
+
+
+      DO ikl=1,klonv        
+
+
+      isnoSV(ikl)  =0.        
+      ispiSV(ikl)  =0.   
+      iiceSV(ikl)  =0.        
+      istoSV(ikl,:)=0.                                                                        
+      alb_SV(ikl)  =0.      
+      emi_SV(ikl)  =0.    
+      IRs_SV(ikl)  =0.      
+      LMO_SV(ikl)  =0.
+      us__SV(ikl)  =0.
+      uts_SV(ikl)  =0.
+      cutsSV(ikl)  =0.
+      uqs_SV(ikl)  =0.
+      uss_SV(ikl)  =0.
+      usthSV(ikl)  =0.
+      rCDmSV(ikl)  =0.
+      rCDhSV(ikl)  =0.
+      Z0m_SV(ikl)  =0.
+      Z0mmSV(ikl)  =0.
+      Z0mnSV(ikl)  =0.
+      Z0roSV(ikl)  =0.
+      Z0SaSV(ikl)  =0.
+      Z0e_SV(ikl)  =0.
+      Z0emSV(ikl)  =0.
+      Z0enSV(ikl)  =0.
+      Z0h_SV(ikl)  =0.
+      Z0hmSV(ikl)  =0.
+      Z0hnSV(ikl)  =0.
+                                                                          
+                                                                        
+      TsisSV(ikl,:)  =0.
+      ro__SV(ikl,:)  =0.
+      eta_SV(ikl,:)  =0.
+      G1snSV(ikl,:)  =0.  
+      G2snSV(ikl,:)  =0.          
+      dzsnSV(ikl,:)  =0.      
+      agsnSV(ikl,:)  =0.              
+      DOPsnSV(ikl,:) =0.                  
+      BufsSV(ikl)  =0.         
+      rusnSV(ikl)  =0.     
+      SWf_SV(ikl)  =0.        
+      SWS_SV(ikl)  =0.
+      HFraSV(ikl)  =0.       
+                                                                            
+      zWE_SV(ikl)  =0.
+      zWEcSV(ikl)  =0.
+      wem_SV(ikl)  =0.
+      wer_SV(ikl)  =0.
+      wes_SV(ikl)  =0.
+      zn4_SV(ikl)  =0.
+      zn5_SV(ikl)  =0.                                                        
+                                       
+                                                                                                                                                          
+      ii__SV(ikl)  =0. 
+      jj__SV(ikl)  =0.
+      nn__SV(ikl)  =0.
+                                                                               
+      IRu_SV(ikl)  =0. 
+      hSalSV(ikl)  =0.     
+      qSalSV(ikl)  =0.
+      RnofSV(ikl)  =0.    
+      RuofSV(ikl,:)  =0.
+
+
+
+
+      END DO
+  END SUBROUTINE INIT_VARxSV
+
+
+                                               
+END MODULE VARxSV
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.F90	(revision 5259)
+++ 	(revision )
@@ -1,153 +1,0 @@
-MODULE VARySV
-
-USE VAR_SV, only : klonv, nsol, nsno
-
-
-IMPLICIT NONE
-
-! +
-      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   NLaysv  ! New   Snow     Layer   Switch   
-!$OMP THREADPRIVATE(NLaysv)
-      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   i_thin  ! Index of the thinest Layer      
-!$OMP THREADPRIVATE(i_thin)
-      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   LIndsv  ! Contiguous Layer relative Index
-!$OMP THREADPRIVATE(LIndsv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albisv  ! Integrated Surface Albedo
-!$OMP THREADPRIVATE(albisv)
-      !hj prelim                                            
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb1sv  ! Surface Albedo VIS      
-!$OMP THREADPRIVATE(alb1sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb2sv  ! Surface Albedo NIR      
-!$OMP THREADPRIVATE(alb2sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb3sv  ! Surface Albedo FIR      
-!$OMP THREADPRIVATE(alb3sv)
-      REAL, DIMENSION(:,:),SAVE,ALLOCATABLE  ::   alb6sv  ! 6 band-albedo    
-!$OMP THREADPRIVATE(alb6sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albssv  ! Soil               Albedo [-]   
-!$OMP THREADPRIVATE(albssv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   SoSosv  !? Surface Absorbed Solar Radiat.  
-!$OMP THREADPRIVATE(SoSosv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Eso_sv  !? Emissivity of       Soil+Snow   
-!$OMP THREADPRIVATE(Eso_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HSv_sv  ! Sensible Heat Flux  [W/m2]      
-!$OMP THREADPRIVATE(HSv_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HLv_sv  ! Latent   Heat Flux  [W/m2]      
-!$OMP THREADPRIVATE(HLv_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HSs_sv  !? Sensible Heat Flux (t)          
-!$OMP THREADPRIVATE(HSs_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HLs_sv  !? Latent   Heat Flux (t)          
-!$OMP THREADPRIVATE(HLs_sv) 
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   sqrCm0  ! in Neutral Drag Coef.Moment.    
-!$OMP THREADPRIVATE(sqrCm0)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   sqrCh0  ! in Neutral Drag Coef.Heat       
-!$OMP THREADPRIVATE(sqrCh0)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Lx_H2O  ! Latent Heat of Vaporiz./Sublim. 
-!$OMP THREADPRIVATE(Lx_H2O)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   ram_sv  ! Aerodyn.Resistance (Moment.)    
-!$OMP THREADPRIVATE(ram_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   rah_sv  ! Aerodyn.Resistance (Heat)       
-!$OMP THREADPRIVATE(rah_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Fh__sv  ! Stability Function              
-!$OMP THREADPRIVATE(Fh__sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   dFh_sv  ! Stability Function (Deriv.)     
-!$OMP THREADPRIVATE(dFh_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Evp_sv  !x Evaporation        [kg/m2]      
-!$OMP THREADPRIVATE(Evp_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   EvT_sv  !x Evapotranspiration [kg/m2]      
-!$OMP THREADPRIVATE(EvT_sv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   LSdzsv  ! Land/Sea Vert. Discretiz. Fact. 
-!$OMP THREADPRIVATE(LSdzsv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Tsrfsv  ! Surface    Temperature          
-!$OMP THREADPRIVATE(Tsrfsv)
-      REAL,SAVE,ALLOCATABLE  ::   sEX_sv(:,:)  ! Verticaly Integr.Extinct.Coef.  
-!$OMP THREADPRIVATE(sEX_sv)
-      REAL,SAVE,ALLOCATABLE  ::   zzsnsv(:,:)    ! Snow  Pack Thickness      [m]   
-!$OMP THREADPRIVATE(zzsnsv)
-      REAL,SAVE,ALLOCATABLE  ::   psi_sv(:,:)    ! Soil   Water        Potential   
-!$OMP THREADPRIVATE(psi_sv)
-      REAL,SAVE,ALLOCATABLE  ::   Khydsv(:,:)    ! Soil   Hydraulic    Conductiv.  
-!$OMP THREADPRIVATE(Khydsv)
-      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   EExcsv  ! Energy in Excess, current       
-!$OMP THREADPRIVATE(EExcsv)
-
-CONTAINS
-  SUBROUTINE INIT_VARySV
-  IMPLICIT NONE
-      INTEGER                ::   ikl
-
-      ALLOCATE(NLaysv(klonv))  ! New   Snow     Layer   Switch   
-      ALLOCATE(i_thin(klonv))  ! Index of the thinest Layer      
-      ALLOCATE(LIndsv(klonv))  ! Contiguous Layer relative Index
-      ALLOCATE(albisv(klonv))  ! Integrated Surface Albedo
-      !hj prelim                                            
-      ALLOCATE(alb1sv(klonv))  ! Surface Albedo VIS      
-      ALLOCATE(alb2sv(klonv))  ! Surface Albedo NIR      
-      ALLOCATE(alb3sv(klonv))  ! Surface Albedo FIR       
-      ALLOCATE(alb6sv(klonv,6))! 6-band  Albedo      
-
-      !
-      ALLOCATE(albssv(klonv))  ! Soil               Albedo [-]   
-      ALLOCATE(SoSosv(klonv))  !? Surface Absorbed Solar Radiat.  
-      ALLOCATE(Eso_sv(klonv))  !? Emissivity of       Soil+Snow   
-      ALLOCATE(HSv_sv(klonv))  ! Sensible Heat Flux  [W/m2]      
-      ALLOCATE(HLv_sv(klonv))  ! Latent   Heat Flux  [W/m2]      
-      ALLOCATE(HSs_sv(klonv))  !? Sensible Heat Flux (t)          
-      ALLOCATE(HLs_sv(klonv))  !? Latent   Heat Flux (t)     
-      ALLOCATE(sqrCm0(klonv))  ! in Neutral Drag Coef.Moment.    
-      ALLOCATE(sqrCh0(klonv))  ! in Neutral Drag Coef.Heat       
-      ALLOCATE(Lx_H2O(klonv))  ! Latent Heat of Vaporiz./Sublim. 
-      ALLOCATE(ram_sv(klonv))  ! Aerodyn.Resistance (Moment.)    
-      ALLOCATE(rah_sv(klonv))  ! Aerodyn.Resistance (Heat)       
-      ALLOCATE(Fh__sv(klonv))  ! Stability Function              
-      ALLOCATE(dFh_sv(klonv))  ! Stability Function (Deriv.)     
-      ALLOCATE(Evp_sv(klonv))  !x Evaporation        [kg/m2]      
-      ALLOCATE(EvT_sv(klonv))  !x Evapotranspiration [kg/m2]      
-      ALLOCATE(LSdzsv(klonv))  ! Land/Sea Vert. Discretiz. Fact. 
-      ALLOCATE(Tsrfsv(klonv))  ! Surface    Temperature          
-      ALLOCATE(sEX_sv(klonv,-nsol:nsno+1))  ! Verticaly Integr.Extinct.Coef.  
-      ALLOCATE(zzsnsv(klonv,    0:nsno))    ! Snow  Pack Thickness      [m]   
-      ALLOCATE(psi_sv(klonv,-nsol:0   ))    ! Soil   Water        Potential   
-      ALLOCATE(Khydsv(klonv,-nsol:0   ))    ! Soil   Hydraulic    Conductiv.  
-      ALLOCATE(EExcsv(klonv))  ! Energy in Excess, current         
-
-      DO ikl=1,klonv
-
-      NLaysv(ikl) =0.
-      i_thin(ikl) =0.
-      LIndsv(ikl) =0.
-      albisv(ikl) =0.
-      alb1sv(ikl) =0. 
-      alb2sv(ikl) =0.   
-      alb3sv(ikl) =0.
-      alb6sv(ikl,:)=0.
-      albssv(ikl) =0.
-      SoSosv(ikl) =0.
-      Eso_sv(ikl) =0.
-      HSv_sv(ikl) =0.   
-      HLv_sv(ikl) =0.
-      HSs_sv(ikl) =0.       
-      HLs_sv(ikl) =0.
-      sqrCm0(ikl) =0.   
-      sqrCh0(ikl) =0.    
-      Lx_H2O(ikl) =0.
-      ram_sv(ikl) =0.
-      rah_sv(ikl) =0.
-      Fh__sv(ikl) =0.      
-      dFh_sv(ikl) =0.
-      Evp_sv(ikl) =0.
-      EvT_sv(ikl) =0.
-      LSdzsv(ikl) =0.
-      Tsrfsv(ikl) =0.
-      sEX_sv(ikl,:)  =0.
-      zzsnsv(ikl,:)  =0.
-      psi_sv(ikl,:)  =0. 
-      Khydsv(ikl,:)  =0. 
-      EExcsv(ikl)  =0.    
-
-
-      END DO
-
-  END SUBROUTINE INIT_VARySV
-
-END MODULE VARySV
-
Index: LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/VARySV.f90	(revision 5268)
@@ -0,0 +1,153 @@
+MODULE VARySV
+
+USE VAR_SV, only : klonv, nsol, nsno
+
+
+IMPLICIT NONE
+
+! +
+      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   NLaysv  ! New   Snow     Layer   Switch   
+!$OMP THREADPRIVATE(NLaysv)
+      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   i_thin  ! Index of the thinest Layer      
+!$OMP THREADPRIVATE(i_thin)
+      INTEGER, DIMENSION(:),SAVE,ALLOCATABLE ::   LIndsv  ! Contiguous Layer relative Index
+!$OMP THREADPRIVATE(LIndsv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albisv  ! Integrated Surface Albedo
+!$OMP THREADPRIVATE(albisv)
+      !hj prelim                                            
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb1sv  ! Surface Albedo VIS      
+!$OMP THREADPRIVATE(alb1sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb2sv  ! Surface Albedo NIR      
+!$OMP THREADPRIVATE(alb2sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   alb3sv  ! Surface Albedo FIR      
+!$OMP THREADPRIVATE(alb3sv)
+      REAL, DIMENSION(:,:),SAVE,ALLOCATABLE  ::   alb6sv  ! 6 band-albedo    
+!$OMP THREADPRIVATE(alb6sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   albssv  ! Soil               Albedo [-]   
+!$OMP THREADPRIVATE(albssv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   SoSosv  !? Surface Absorbed Solar Radiat.  
+!$OMP THREADPRIVATE(SoSosv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Eso_sv  !? Emissivity of       Soil+Snow   
+!$OMP THREADPRIVATE(Eso_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HSv_sv  ! Sensible Heat Flux  [W/m2]      
+!$OMP THREADPRIVATE(HSv_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HLv_sv  ! Latent   Heat Flux  [W/m2]      
+!$OMP THREADPRIVATE(HLv_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HSs_sv  !? Sensible Heat Flux (t)          
+!$OMP THREADPRIVATE(HSs_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   HLs_sv  !? Latent   Heat Flux (t)          
+!$OMP THREADPRIVATE(HLs_sv) 
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   sqrCm0  ! in Neutral Drag Coef.Moment.    
+!$OMP THREADPRIVATE(sqrCm0)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   sqrCh0  ! in Neutral Drag Coef.Heat       
+!$OMP THREADPRIVATE(sqrCh0)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Lx_H2O  ! Latent Heat of Vaporiz./Sublim. 
+!$OMP THREADPRIVATE(Lx_H2O)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   ram_sv  ! Aerodyn.Resistance (Moment.)    
+!$OMP THREADPRIVATE(ram_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   rah_sv  ! Aerodyn.Resistance (Heat)       
+!$OMP THREADPRIVATE(rah_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Fh__sv  ! Stability Function              
+!$OMP THREADPRIVATE(Fh__sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   dFh_sv  ! Stability Function (Deriv.)     
+!$OMP THREADPRIVATE(dFh_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Evp_sv  !x Evaporation        [kg/m2]      
+!$OMP THREADPRIVATE(Evp_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   EvT_sv  !x Evapotranspiration [kg/m2]      
+!$OMP THREADPRIVATE(EvT_sv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   LSdzsv  ! Land/Sea Vert. Discretiz. Fact. 
+!$OMP THREADPRIVATE(LSdzsv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   Tsrfsv  ! Surface    Temperature          
+!$OMP THREADPRIVATE(Tsrfsv)
+      REAL,SAVE,ALLOCATABLE  ::   sEX_sv(:,:)  ! Verticaly Integr.Extinct.Coef.  
+!$OMP THREADPRIVATE(sEX_sv)
+      REAL,SAVE,ALLOCATABLE  ::   zzsnsv(:,:)    ! Snow  Pack Thickness      [m]   
+!$OMP THREADPRIVATE(zzsnsv)
+      REAL,SAVE,ALLOCATABLE  ::   psi_sv(:,:)    ! Soil   Water        Potential   
+!$OMP THREADPRIVATE(psi_sv)
+      REAL,SAVE,ALLOCATABLE  ::   Khydsv(:,:)    ! Soil   Hydraulic    Conductiv.  
+!$OMP THREADPRIVATE(Khydsv)
+      REAL, DIMENSION(:),SAVE,ALLOCATABLE    ::   EExcsv  ! Energy in Excess, current       
+!$OMP THREADPRIVATE(EExcsv)
+
+CONTAINS
+  SUBROUTINE INIT_VARySV
+  IMPLICIT NONE
+      INTEGER                ::   ikl
+
+      ALLOCATE(NLaysv(klonv))  ! New   Snow     Layer   Switch   
+      ALLOCATE(i_thin(klonv))  ! Index of the thinest Layer      
+      ALLOCATE(LIndsv(klonv))  ! Contiguous Layer relative Index
+      ALLOCATE(albisv(klonv))  ! Integrated Surface Albedo
+      !hj prelim                                            
+      ALLOCATE(alb1sv(klonv))  ! Surface Albedo VIS      
+      ALLOCATE(alb2sv(klonv))  ! Surface Albedo NIR      
+      ALLOCATE(alb3sv(klonv))  ! Surface Albedo FIR       
+      ALLOCATE(alb6sv(klonv,6))! 6-band  Albedo      
+
+      !
+      ALLOCATE(albssv(klonv))  ! Soil               Albedo [-]   
+      ALLOCATE(SoSosv(klonv))  !? Surface Absorbed Solar Radiat.  
+      ALLOCATE(Eso_sv(klonv))  !? Emissivity of       Soil+Snow   
+      ALLOCATE(HSv_sv(klonv))  ! Sensible Heat Flux  [W/m2]      
+      ALLOCATE(HLv_sv(klonv))  ! Latent   Heat Flux  [W/m2]      
+      ALLOCATE(HSs_sv(klonv))  !? Sensible Heat Flux (t)          
+      ALLOCATE(HLs_sv(klonv))  !? Latent   Heat Flux (t)     
+      ALLOCATE(sqrCm0(klonv))  ! in Neutral Drag Coef.Moment.    
+      ALLOCATE(sqrCh0(klonv))  ! in Neutral Drag Coef.Heat       
+      ALLOCATE(Lx_H2O(klonv))  ! Latent Heat of Vaporiz./Sublim. 
+      ALLOCATE(ram_sv(klonv))  ! Aerodyn.Resistance (Moment.)    
+      ALLOCATE(rah_sv(klonv))  ! Aerodyn.Resistance (Heat)       
+      ALLOCATE(Fh__sv(klonv))  ! Stability Function              
+      ALLOCATE(dFh_sv(klonv))  ! Stability Function (Deriv.)     
+      ALLOCATE(Evp_sv(klonv))  !x Evaporation        [kg/m2]      
+      ALLOCATE(EvT_sv(klonv))  !x Evapotranspiration [kg/m2]      
+      ALLOCATE(LSdzsv(klonv))  ! Land/Sea Vert. Discretiz. Fact. 
+      ALLOCATE(Tsrfsv(klonv))  ! Surface    Temperature          
+      ALLOCATE(sEX_sv(klonv,-nsol:nsno+1))  ! Verticaly Integr.Extinct.Coef.  
+      ALLOCATE(zzsnsv(klonv,    0:nsno))    ! Snow  Pack Thickness      [m]   
+      ALLOCATE(psi_sv(klonv,-nsol:0   ))    ! Soil   Water        Potential   
+      ALLOCATE(Khydsv(klonv,-nsol:0   ))    ! Soil   Hydraulic    Conductiv.  
+      ALLOCATE(EExcsv(klonv))  ! Energy in Excess, current         
+
+      DO ikl=1,klonv
+
+      NLaysv(ikl) =0.
+      i_thin(ikl) =0.
+      LIndsv(ikl) =0.
+      albisv(ikl) =0.
+      alb1sv(ikl) =0. 
+      alb2sv(ikl) =0.   
+      alb3sv(ikl) =0.
+      alb6sv(ikl,:)=0.
+      albssv(ikl) =0.
+      SoSosv(ikl) =0.
+      Eso_sv(ikl) =0.
+      HSv_sv(ikl) =0.   
+      HLv_sv(ikl) =0.
+      HSs_sv(ikl) =0.       
+      HLs_sv(ikl) =0.
+      sqrCm0(ikl) =0.   
+      sqrCh0(ikl) =0.    
+      Lx_H2O(ikl) =0.
+      ram_sv(ikl) =0.
+      rah_sv(ikl) =0.
+      Fh__sv(ikl) =0.      
+      dFh_sv(ikl) =0.
+      Evp_sv(ikl) =0.
+      EvT_sv(ikl) =0.
+      LSdzsv(ikl) =0.
+      Tsrfsv(ikl) =0.
+      sEX_sv(ikl,:)  =0.
+      zzsnsv(ikl,:)  =0.
+      psi_sv(ikl,:)  =0. 
+      Khydsv(ikl,:)  =0. 
+      EExcsv(ikl)  =0.    
+
+
+      END DO
+
+  END SUBROUTINE INIT_VARySV
+
+END MODULE VARySV
+
Index: LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.F90	(revision 5259)
+++ 	(revision )
@@ -1,1495 +1,0 @@
-MODULE surf_inlandsis_mod
-
-    IMPLICIT NONE
-
-CONTAINS
-
-
-SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
-            rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
-            precip_rain, precip_snow, &
-            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
-            rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
-            radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
-            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
-            runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
-            tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)
-
-        ! |                                                                        |
-        ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
-        ! |                              (INLANDSIS)                               |
-        ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |
-        ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
-        ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
-        ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
-        ! |   Updated by Etienne Vignon, Cecile Agosta                             |
-        ! |                                                                        |
-        ! +------------------------------------------------------------------------+
-        ! |
-        ! |   In the current setup, SISVAT is used only to model the land ice      |
-        ! |   part of the surface; hence it is called with the compressed variables|
-        ! |   from pbl_surface, and only by the surf_landice routine.              |
-        ! |                                                                        |
-        ! |   In this interface it is assumed that the partitioning of the soil,   |
-        ! |   and hence the number of grid points is constant during a simulation, |
-        ! |   hence eg. snow properties remain stored in the global SISVAT         |
-        ! |   variables between the calls and don't need to be handed over as      |
-        ! |   arguments. When the partitioning is supposed to change, make sure to |
-        ! |   update the variables.                                                |
-        ! |                                                                        |
-        ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV ...)                    |
-        ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |
-        ! |                                                                        |
-        ! +------------------------------------------------------------------------+
-
-        USE dimphy
-        USE VAR_SV
-        USE VARdSV
-        USE VARxSV
-        USE VARySV
-        USE VARtSV
-        USE VARphy
-        USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor
-
-        IMPLICIT NONE
-
-        ! +--INTERFACE Variables
-        ! +  ===================
-        !    include  "dimsoil.h"
-
-        ! +--Global Variables
-        ! +  ================
-        ! Input Variables for SISVAT
-        INTEGER, INTENT(IN) :: knon
-        INTEGER, INTENT(IN) :: itime
-        REAL, INTENT(IN) :: dtime
-        LOGICAL, INTENT(IN) :: debut     ! true if first step
-        LOGICAL, INTENT(IN) :: lafin     ! true if last step
-
-        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i     ! Index Decompression
-        REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
-        REAL, DIMENSION(klon), INTENT(IN) :: rmu0      ! cos sol. zenith angle
-        REAL, DIMENSION(klon), INTENT(IN) :: swdown    !
-        REAL, DIMENSION(klon), INTENT(IN) :: lwdown    !
-        REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
-        REAL, DIMENSION(klon), INTENT(IN) :: pexner    ! Exner potential
-        REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
-        REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
-        REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
-        REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
-        REAL, DIMENSION(klon), INTENT(IN) :: rugos
-        REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
-        REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
-        REAL, DIMENSION(klon), INTENT(IN) :: alt       ! surface elevation
-        REAL, DIMENSION(klon), INTENT(IN) :: cloudf
-        REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
-        REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
-        REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
-        REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! friction velocity
-
-        ! Variables exchanged between LMDZ and SISVAT
-        REAL, DIMENSION(klon), INTENT(IN) :: radsol    ! Surface absorbed rad.
-        REAL, DIMENSION(klon), INTENT(INOUT) :: snow      ! Tot snow mass [kg/m2]
-        REAL, DIMENSION(klon), INTENT(INOUT) :: zfra      ! snwo surface fraction [0-1]
-        REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
-        REAL, DIMENSION(klon), INTENT(OUT) :: qsol      ! Soil Water Content
-        REAL, DIMENSION(klon), INTENT(INOUT) :: z0m    ! Momentum Roughn Lgt
-        REAL, DIMENSION(klon), INTENT(INOUT) :: z0h    ! Momentum Roughn Lgt
-
-        ! Output Variables for LMDZ
-        REAL, DIMENSION(klon), INTENT(OUT) :: alb1      ! Albedo SW
-        REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
-        REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
-        REAL, DIMENSION(klon), INTENT(OUT) :: emis_new  ! Surface Emissivity
-        REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
-        REAL, DIMENSION(klon), INTENT(OUT) :: ffonte    ! enthalpy flux due to surface melting
-        REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte   ! water flux due to surface melting
-        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s   ! d/dT sens. ht flux
-        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l   ! d/dT latent ht flux
-        REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens  ! Sensible ht flux
-        REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat   ! Latent heat flux
-        REAL, DIMENSION(klon), INTENT(OUT) :: evap      ! Evaporation
-        REAL, DIMENSION(klon), INTENT(OUT) :: erod      ! Erosion of surface snow (flux)
-        REAL, DIMENSION(klon), INTENT(OUT) :: agesno    ! Snow age (top layer)
-        REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
-        REAL, DIMENSION(klon), INTENT(OUT) :: qsurf     ! Surface Humidity
-
-        ! Specific INLANDIS outputs
-        REAL, DIMENSION(klon), INTENT(OUT) :: qsnow     ! Total H2O snow[kg/m2]
-        REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt   ! Snow height (m)
-        REAL, DIMENSION(klon), INTENT(OUT) :: to_ice    ! Snow passed to ice
-        REAL, DIMENSION(klon), INTENT(OUT) :: sissnow   ! Snow in model (kg/m2)
-
-        ! +--Internal  Variables
-        ! +  ===================
-
-        CHARACTER(len = 20) :: fn_outfor ! Name for output file
-        CHARACTER (len = 80)              :: abort_message
-        CHARACTER (len = 20)              :: modname = 'surf_inlandsis_mod'
-
-        INTEGER :: i, ig, ikl, isl, isn, nt
-        INTEGER :: gp_outfor, un_outfor
-        REAL, PARAMETER :: f1 = 0.5
-        REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
-        REAL, PARAMETER :: sn_add = 400., sn_div = 2.
-        ! snow mass upper,lower limit,
-        ! added mass/division lowest layer
-        REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
-        REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
-        ! Parameters for drainage
-        ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
-        ! +...        Run Off Parameters
-        ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)
-        ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
-
-        REAL, DIMENSION(klon) :: eps0SL          ! surface Emissivity
-        REAL :: zsigma, Ua_min, Us_min, lati
-        REAL, PARAMETER :: cdmax=0.05
-        REAL :: lambda          ! Par. soil discret.
-        REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2         ! Soil layer thicknesses
-        !$OMP THREADPRIVATE(dz1,dz2)
-        LOGICAL, SAVE :: firstcall
-        !$OMP THREADPRIVATE(firstcall)
-
-        INTEGER :: iso
-        LOGICAL :: file_exists
-        CHARACTER(len = 20) :: fichnom
-        LOGICAL :: is_init_domec
-        ! CA initialization
-        ! dz_profil_15 : 1 m in 15 layers [m]
-        real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
-                                                0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
-        ! mean_temp : mean annual surface temperature [K]
-        real, dimension(klon) :: mean_temp
-        ! mean_dens : mean surface density [kg/m3]
-        real, dimension(klon) :: mean_dens
-        ! lat_scale : temperature lapse rate against latitude [K degree-1]
-        real :: lat_scale
-        ! sh_scale : temperature lapse rate against altitude [K km-1]
-        real :: sh_scale
-        ! variables for density profile
-        ! E0, E1 : exponent
-        real :: E0, E1
-        ! depth at which 550 kg m-3 is reached [m]
-        real :: z550
-        ! depths of snow layers
-        real :: depth, snow_depth, distup
-        ! number of initial snow layers
-        integer :: nb_snow_layer
-        ! For density calc.
-        real :: alpha0, alpha1, ln_smb
-        ! theoritical densities [kg m-3]
-        real :: rho0, rho1, rho1_550
-        ! constants for density profile
-        ! C0, C1 : constant, 0.07 for z <= 550 kg m-3
-        real, parameter :: C0 = 0.07
-        real, parameter :: C1 = 0.03
-        ! rho_i : ice density [kg m-3]
-        real, parameter :: rho_ice = 917.
-        ! E_c : activation energy [J mol-1]
-        real, parameter :: E_c = 60000.
-        ! E_g : activation energy [J mol-1]
-        real, parameter :: E_g = 42400.
-        ! R : gas constant [J mol-1 K-1]
-        real, parameter :: R = 8.3144621
-
-      
-     
-
-
-        ! + PROGRAM START
-        ! + -----------------------------------------
-
-        zsigma = 1000.
-        dt__SV = dtime
-
-        IF (debut) THEN
-            firstcall = .TRUE.
-            INI_SV = .false.
-        ELSE
-            firstcall = .false.
-            INI_SV = .true.
-        END IF
-
-        IF (ok_outfor) THEN
-            un_outfor = 51    ! unit number for point output file
-            gp_outfor = 1    ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
-            fn_outfor = 'outfor_SV.dat'
-        END IF
-
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! + INITIALISATION: BEGIN +++
-        ! + -----------------------------------------
-        IF (firstcall) THEN
-
-            ! +--Array size
-            ! +  -----------------------
-
-            klonv = klon
-            knonv = knon
-                write(*, *) 'ikl, lon and lat in INLANDSIS'
-
-            DO ikl = 1, knon
-                i=ikl2i(ikl)
-                write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
-            END DO
-
-            ! +--Variables initizialisation
-            ! +  ---------------------------
-
-            CALL INIT_VARtSV
-            CALL INIT_VARxSV
-            CALL INIT_VARySV
-
-
-
-            ! +--Surface Fall Line Slope
-            ! +  -----------------------
-            IF (SnoMod)  THEN
-                DO ikl = 1, knon
-                    slopSV(ikl) = slope(ikl)
-                    SWf_SV(ikl) = &   ! Normalized Decay of the
-                            exp(-dt__SV             &   ! Surficial Water Content
-                                    / (c1_zuo                &   !(Zuo and Oerlemans 1996,
-                                            + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
-                END DO
-            END IF
-
-
-
-            ! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
-            ! +  ----------------------------------------------------------------
-            !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
-            CALL get_soil_levels(dz1, dz2, lambda)
-
-            lambSV = lambda
-            dz1_SV(1:knon, 1:) = 0.
-            dz2_SV(1:knon, 1:) = 0.
-
-            DO isl = -nsol, 0
-                dz_dSV(isl) = 0.5e-3 * dz2(1 - isl)           ! Soil layer thickness
-                DO ikl = 1, knon
-                    dz1_SV(ikl, isl) = dz1(1 - isl)    !1.e-3*
-                    dz2_SV(ikl, isl) = dz2(1 - isl)    !1.e-3*
-                END DO
-            END DO
-
-
-            ! Set variables
-            ! =============
-            DO ikl = 1, knon
-                ! LSmask : Land/Sea Mask
-                LSmask(ikl) = 1
-                ! isotSV : Soil Type -> 12 = ice
-                isotSV(ikl) = 12
-                ! iWaFSV : Soil Drainage (1,0)=(y,n)
-                iWaFSV(ikl) = 1
-                ! eps0SL : Surface Emissivity
-                eps0SL(ikl) = 1.
-                ! alb0SV : Soil Albedo
-                alb0SV(ikl) = alb_soil(ikl)
-                ! Tsf_SV : Surface Temperature, must be bellow freezing
-                Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
-            END DO
-
-            ! +--Initialization of soil and snow variables in case startsis is not read
-            ! +  ----------------------------------------------------------------------
-
-            is_init_domec=.FALSE.
-
-
-            IF (is_init_domec) THEN
-            ! Coarse initilization inspired from vertcical profiles at Dome C,
-            ! Antarctic Plateaui (10m of snow, 19 levels)
-
-                 DO ikl = 1,knon
-! + Soil
-                 DO isl =   -nsol,0    
-                   TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)
-                   !tsoil(ikl,1-isl)   Soil Temperature
-                   !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
-                   eta_SV(ikl,isl) = epsi           !etasoil(ikl,1-isl) Soil Water[m3/m3]
-                   ro__SV(ikl,isl) = rhoIce         !rosoil(ikl,1-isl) volumic mass
-                 END DO   
-
-
-           ! Snow
-                 isnoSV(ikl) = 19
-                 istoSV(ikl, 1:isnoSV(ikl)) = 100
-                 ro__SV(ikl, 1:isnoSV(ikl)) = 350.
-                 eta_SV(ikl, 1:isnoSV(ikl)) = epsi
-                 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
-                 G1snSV(ikl, 1:isnoSV(ikl)) = 99.
-                 G2snSV(ikl, 1:isnoSV(ikl)) = 2.
-                 agsnSV(ikl, 1:isnoSV(ikl)) = 50.
-                 dzsnSV(ikl, 19) = 0.015
-                 dzsnSV(ikl, 18) = 0.015
-                 dzsnSV(ikl, 17) = 0.020
-                 dzsnSV(ikl, 16) = 0.030
-                 dzsnSV(ikl, 15) = 0.040
-                 dzsnSV(ikl, 14) = 0.060
-                 dzsnSV(ikl, 13) = 0.080
-                 dzsnSV(ikl, 12) = 0.110
-                 dzsnSV(ikl, 11) = 0.150
-                 dzsnSV(ikl, 10) = 0.200
-                 dzsnSV(ikl, 9) = 0.300
-                 dzsnSV(ikl, 8) = 0.420
-                 dzsnSV(ikl, 7) = 0.780
-                 dzsnSV(ikl, 6) = 1.020
-                 dzsnSV(ikl, 5) = 0.980
-                 dzsnSV(ikl, 4) = 1.020
-                 dzsnSV(ikl, 3) = 3.970
-                 dzsnSV(ikl, 2) = 1.020
-                 dzsnSV(ikl, 1) = 1.020
-
-                 END DO
-            ELSE
-
-            ! Initilialisation with climatological temperature and density
-            ! profiles as in MAR. Methodology developed by Cecile Agosta
-  
-            ! initialize with 0., for unused snow layers
-            dzsnSV = 0.
-            G1snSV = 0.
-            G2snSV = 0.
-            istoSV = 0
-            TsisSV = 0.
-
-
-            ! initialize mean variables (unrealistic)
-            mean_temp = TfSnow
-            mean_dens = 300.
-            ! loop on grid cells
-            DO ikl = 1, knon
-            	lati=rlat(ikl2i(ikl))
-                ! approximations for mean_temp and mean_dens
-                ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
-                ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
-                ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
-                ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
-                if (lati > 60.) then
-                    ! CA todo : add longitude bounds
-                    ! Greenland mean temperature : function of altitude and latitude
-                    ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
-                    lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
-                    lat_scale = max(min(lat_scale, 0.9), 0.75)
-                    ! sh_scale equals the environmental lapse rate : 6.5 °C km-1
-                    sh_scale = 6.5
-                    mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
-                    ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
-                    mean_dens(ikl) = 315.
-                else if (lati < -60.) then
-                    ! Antarctica mean temperature : function of altitude and latitude
-                    ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
-                    lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
-                    lat_scale = max(min(lat_scale, 1.3), 0.6)
-                    ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
-                    sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
-                    sh_scale = max(min(sh_scale, 9.8), 6.5)
-                    mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
-                    ! Antarctica surface density : function of mean annual temperature
-                    ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
-                    ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
-                    ! Weinhart et al 2020  https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
-                    ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
-                    ! Dome C : st_ant_param(3233, -75.1) = -47.7
-                    ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
-                    mean_dens(ikl) =  (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
-                    mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
-                else
-
-                !    write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica'
-
-                     mean_dens(ikl) =350.
-                     mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)
-
-                !abort_message='temperature initialization is only defined for Greenland and Antarctica'
-                !CALL abort_physic(modname,abort_message,1)
-
-                end if
- 
-                ! mean_temp is defined for ice ground only
-                mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)
-
-                ! Soil layers
-                ! ===========
-                DO isl = -nsol, 0
-                    ! TsisSV : Temperature [K]
-                    TsisSV(ikl, isl) = mean_temp(ikl)
-                    ! eta_SV : Soil Water [m3/m3]
-                    eta_SV(ikl, isl) = epsi
-                    ! ro__SV : Volumic Mass [kg/m3]
-                    ro__SV(ikl, isl) = rhoIce
-                END DO
-
-                ! Snow layers
-                ! ===========
-                ! snow_depth : initial snow depth
-                snow_depth = 20.
-                ! nb_snow_layer : initial nb of snow layers
-                nb_snow_layer = 15
-                ! isnoSV : total nb of snow layers
-                isnoSV(ikl) = nb_snow_layer
-                ! depth : depth of each layer
-                depth = snow_depth
-                do isl = 1, nb_snow_layer
-                    ! dzsnSV : snow layer thickness
-                    dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
-                    ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
-                    G1snSV(ikl, isl) = 99.
-                    ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
-                    G2snSV(ikl, isl) = 3.
-                    agsnSV(ikl, isl) = 0.
-                    istoSV(ikl, isl) = 0
-                    ! eta_SV : Liquid Water Content [m3/m3]
-                    eta_SV(ikl, isl) = 0.
-                    ! distance to surface
-                    depth = depth - dzsnSV(ikl,isl) / 2.
-                    distup = min(1., max(0., depth / snow_depth))
-                    ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
-                    TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2
-                    ! firn density : densification formulas from :
-                    ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
-                    ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
-                    ! Integration of the steady state equation
-                    ! ln_smb approximated as a function of temperature
-                    ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
-                    ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
-                    alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
-                    alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
-                    E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
-                    E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
-                    z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
-                    rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
-                    rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
-                    if (depth <= z550) then
-                        ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
-                    else
-                        ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
-                    end if
-                    depth = depth - dzsnSV(ikl,isl) / 2.
-                    
-                end do
-
-            END DO
-
-            END IF
-
-
-            ! + Numerics paramaters, SISVAT_ini 
-            ! +  ----------------------
-            CALL SISVAT_ini(knon)
-
-
-            ! +--Read restart file
-            ! +  =================================================
-
-            INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
-            IF (file_exists) THEN
-                CALL sisvatetat0("startsis.nc", ikl2i)
-            END IF
-
-
-
-            ! +--Output ascii file
-            ! +  =================================================
-
-            ! open output file
-            IF (ok_outfor) THEN
-                open(unit = un_outfor, status = 'replace', file = fn_outfor)
-                ikl = gp_outfor     ! index sur la grille land ice
-                write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
-                write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
-                        & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
-            END IF
-
-        END IF  ! firstcall
-        ! +
-        ! +  +++  INITIALISATION:  END  +++
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-
-
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! + READ FORCINGS
-        ! + ------------------------
-
-        ! + Update Forcings for SISVAT given by the LMDZ model.
-        ! +
-        DO ikl = 1, knon
-
-            ! +--Atmospheric Forcing                                    (INPUT)
-            ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
-            za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
-            Ua_min = 0.2 * sqrt(za__SV(ikl))            !
-            VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
-            TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
-            ExnrSV(ikl) = pexner(ikl)                   ! Exner potential
-            rhT_SV(ikl) = dens_air(ikl)                 ! Air density
-            QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
-            ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
-            p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
-
-            ! +--Surface properties
-            ! +  ^^^^^^^^^^^^^^^^^^
-
-            Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
-            Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
-
-            ! +--Energy Fluxes                                          (INPUT)
-            ! +  ^^^^^^^^^^^^^                                           ^^^^^
-            coszSV(ikl) = max(czemin, rmu0(ikl))         ! cos(zenith.Dist.)
-            sol_SV(ikl) = swdown(ikl)                   ! downward Solar
-            IRd_SV(ikl) = lwdown(ikl)                   ! downward IR
-            rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.
-
-            ! +--Water  Fluxes                                          (INPUT)
-            ! +  ^^^^^^^^^^^^^                                           ^^^^^
-            drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
-            dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
-
-            ! #BS    dbs_SV(ikl) = blowSN(i,j,n)
-            ! dbs_SV = Maximum potential erosion amount [kg/m2]
-            ! => Upper bound for eroded snow mass
-            !        uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f)
-            ! #BS  if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
-            ! #BS    dsnbSV(ikl) =1.0-min(qsHY(i,j,kB)     !BS neglib. at kb ~100 magl)
-            ! #BS.                        /max(qshy(i,j,mz),eps9),unun)
-            ! #BS    dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
-            ! #BS    dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
-            ! #BS  else
-            ! #BS    dsnbSV(ikl) = 0.
-            ! #BS  endif
-            !      dsnbSV is the drift fraction of deposited snow updated in sisvat.f
-            !      will be used for characterizing the Buffer Layer
-            !      (see update of  Bros_N, G1same, G2same, zroOLD, zroNEW)
-            ! #BS  if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
-            qsnoSV(ikl) = snow_cont_air(ikl)
-
-
-
-            ! +--Soil/BL                                      (INPUT)
-            ! +  ^^^^^^^                                       ^^^^^
-            alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
-            AcoHSV(ikl) = AcoefH(ikl)
-            BcoHSV(ikl) = BcoefH(ikl)
-            AcoQSV(ikl) = AcoefQ(ikl)
-            BcoQSV(ikl) = BcoefQ(ikl)
-            cdH_SV(ikl) = min(cdragh(ikl),cdmax)
-            cdM_SV(ikl) = min(cdragm(ikl),cdmax)
-            rcdmSV(ikl) = sqrt(cdM_SV(ikl))
-            Us_min = 0.01
-            us__SV(ikl) = max(Us_min, ustar(ikl))
-            ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
-            rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))
-
-            ! +--Energy Fluxes                                          (INPUT/OUTPUT)
-            ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^
-            !IF (.not.firstcall) THEN
-            Tsrfsv(ikl)  = tsurf(ikl)                     !hj 12 03 2010
-            cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness
-            !END IF
-
-         END DO
-
-        !
-        ! +  +++  READ FORCINGS:  END  +++
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- 
-
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! +--SISVAT EXECUTION
-        ! +  ----------------
-
-        call  INLANDSIS(SnoMod, BloMod, 1)
-
-
-       
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        ! + RETURN RESULTS
-        ! + --------------
-        ! + Return (compressed) SISVAT variables to LMDZ
-        ! +
-        DO  ikl = 1, knon                  ! use only 1:knon (actual ice sheet..)
-            dflux_s(ikl) = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
-            dflux_l(ikl) = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
-            fluxsens(ikl) = HSs_sv(ikl)         ! HS
-            fluxlat(ikl) = HLs_sv(ikl)         ! HL
-            evap(ikl) = -1*HLs_sv(ikl) / LHvH2O  ! Evaporation
-            erod(ikl) = 0.
-
-            IF (BloMod) THEN
-                ! + Blowing snow
-
-                !       SLussl(i,j,n)= 0.
-                ! #BS   SLussl(i,j,n)=                     !Effective erosion
-                ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl))    !~u*qs* from previous time step
-                ! #BS   blowSN(i,j,n)=  dt*uss_SV(ikl)     !New max. pot. Erosion [kg/m2]
-                ! #BS.                    *rhT_SV(ikl)     !(further bounded in sisvat_bsn.f)
-                ! #BS  erprev(i,j,n) =     dbs_Er(ikl)/dt__SV
-                erod(ikl) = dbs_Er(ikl) / dt__SV
-            ENDIF
-
-            ! +   Check snow thickness,  substract if too thick, add if too thin
-
-            sissnow(ikl) = 0.  !()
-            DO  isn = 1, isnoSV(ikl)
-                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
-            END DO
-
-            IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
-                IF (isnoSV(ikl).GE.1) THEN
-                    dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
-                    toicSV(ikl) = toicSV(ikl) - sn_add
-                ELSE
-                    write(*, *) 'Attention, bare ice... point ', ikl
-                    isnoSV(ikl) = 1
-                    istoSV(ikl, 1) = 0
-                    ro__SV(ikl, 1) = 350.
-                    dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi)  ! 1.
-                    eta_SV(ikl, 1) = epsi
-                    TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
-                    G1snSV(ikl, 1) = 0.
-                    G2snSV(ikl, 1) = 0.3
-                    agsnSV(ikl, 1) = 10.
-                    toicSV(ikl) = toicSV(ikl) - sn_add
-                END IF
-            END IF
-
-            IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
-                dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
-                toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
-            END IF
-
-            sissnow(ikl) = 0.
-            qsnow(ikl) = 0.
-            snow(ikl) = 0.
-            snowhgt(ikl) = 0.
-
-            DO  isn = 1, isnoSV(ikl)
-                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
-                snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
-                ! Etienne: check calc qsnow
-                qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
-            END DO
-
-            zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
-            ! Etienne: comment following line
-            ! snow(ikl)    = sissnow(ikl)+toicSV(ikl)
-            snow(ikl) = sissnow(ikl)
-
-            to_ice(ikl) = toicSV(ikl)
-            runoff_lic(ikl) = RnofSV(ikl)    ! RunOFF: intensity (flux due to melting + liquid precip)
-            fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
-            ffonte(ikl)=fqfonte(ikl)*Lf_H2O
-
-            qsol(ikl) = 0.
-            DO  isl = -nsol, 0
-                tsoil(ikl, 1 - isl) = TsisSV(ikl, isl)       ! Soil Temperature
-                ! Etienne: check calc qsol
-                qsol(ikl) = qsol(ikl)                      &
-                        + eta_SV(ikl, isl) * dz_dSV(isl)
-            END DO
-            agesno(ikl) = agsnSV(ikl, isnoSV(ikl))        !          [day]
-
-            alb1(ikl) = alb1sv(ikl)             ! Albedo VIS
-!            alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl)                   &
-!                    & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
-            alb2(ikl)=alb2sv(ikl)
-            ! Albedo NIR
-            alb3(ikl) = alb3sv(ikl)             ! Albedo FIR
-            ! 6 band Albedo
-            alb6(ikl,:)=alb6sv(ikl,:)
-
-            tsurf_new(ikl) = Tsrfsv(ikl)
-
-            qsurf(ikl) = QsT_SV(ikl)
-            emis_new(ikl) = eps0SL(ikl)
-            z0m(ikl) = Z0m_SV(ikl)
-            z0h(ikl) = Z0h_SV(ikl)
-
-
-        END DO
-
-            IF (ok_outfor) THEN
-             ikl= gp_outfor
-            write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
-            write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
-            write(un_outfor, *) dzsnSV(ikl, :)
-            write(un_outfor, *) TsisSV(ikl, :)
-            write(un_outfor, *) ro__SV(ikl, :)
-            write(un_outfor, *) eta_SV(ikl, :)
-            write(un_outfor, *) G1snSV(ikl, :)
-            write(un_outfor, *) G2snSV(ikl, :)
-            write(un_outfor, *) agsnSV(ikl, :)
-            write(un_outfor, *) istoSV(ikl, :)
-            write(un_outfor, *) DOPsnSV(ikl, :)
-        ENDIF
-
-
-
-        ! +  -----------------------------
-        ! +  END --- RETURN RESULTS
-        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        IF (lafin) THEN
-            fichnom = "restartsis.nc"
-            CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
-
-            IF (ok_outfor) THEN
-                close(unit = un_outfor)
-            END IF
-        END IF
-
-    END SUBROUTINE surf_inlandsis
-
-
-    !=======================================================================
-
-    SUBROUTINE get_soil_levels(dz1, dz2, lambda)
-        ! ======================================================================
-        ! Routine to compute the vertical discretization of the soil in analogy
-        ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
-        ! of SISVAT, therefore it's needed here.
-        !
-        USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
-        USE mod_phys_lmdz_para
-        USE VAR_SV
-
-
-        !    INCLUDE "dimsoil.h"
-
-        REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
-        REAL, INTENT(OUT) :: lambda
-
-
-        !-----------------------------------------------------------------------
-        !   Depthts:
-        !   --------
-        REAL fz, rk, fz1, rk1, rk2
-        REAL min_period, dalph_soil
-        INTEGER ierr, jk
-
-        fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
-
-        !    write(*,*)'Start soil level computation'
-        !-----------------------------------------------------------------------
-        ! Calculation of some constants
-        ! NB! These constants do not depend on the sub-surfaces
-        !-----------------------------------------------------------------------
-        !-----------------------------------------------------------------------
-        !   ground levels
-        !   grnd=z/l where l is the skin depth of the diurnal cycle:
-        !-----------------------------------------------------------------------
-
-        min_period = 1800. ! en secondes
-        dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
-        ! !$OMP MASTER
-        !     IF (is_mpi_root) THEN
-        !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
-        !        IF (ierr == 0) THEN ! Read file only if it exists
-        !           READ(99,*) min_period
-        !           READ(99,*) dalph_soil
-        !           PRINT*,'Discretization for the soil model'
-        !           PRINT*,'First level e-folding depth',min_period, &
-        !                '   dalph',dalph_soil
-        !           CLOSE(99)
-        !        END IF
-        !     ENDIF
-        ! !$OMP END MASTER
-        !     CALL bcast(min_period)
-        !     CALL bcast(dalph_soil)
-
-        !   la premiere couche represente un dixieme de cycle diurne
-        fz1 = SQRT(min_period / 3.14)
-
-        DO jk = 1, nsoilmx
-            rk1 = jk
-            rk2 = jk - 1
-            dz2(jk) = fz(rk1) - fz(rk2)
-        ENDDO
-        DO jk = 1, nsoilmx - 1
-            rk1 = jk + .5
-            rk2 = jk - .5
-            dz1(jk) = 1. / (fz(rk1) - fz(rk2))
-        ENDDO
-        lambda = fz(.5) * dz1(1)
-        DO jk = 1, nsoilmx
-            rk = jk
-            rk1 = jk + .5
-            rk2 = jk - .5
-        ENDDO
-
-    END SUBROUTINE get_soil_levels
-
-
-    !===========================================================================
-
-    SUBROUTINE SISVAT_ini(knon)
-
-        !C +------------------------------------------------------------------------+
-        !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
-        !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |
-        !C +------------------------------------------------------------------------+
-        !C |   PARAMETERS:  klonv: Total Number of columns =                        |
-        !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
-        !C |                     X       Number of Mosaic Cell per grid box         |
-        !C |                                                                        |
-        !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
-        !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |
-        !C |                                                                        |
-        !C |   OUTPUT:             [-] |
-        !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |
-        !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |
-        !C |                      (based on a prescribed Soil Relative Humidity)    |
-        !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |
-        !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |
-        !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |
-        !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |
-        !C |            dzsnSV(0): Soil first Layer Thickness                   [m] |
-        !C |            dzmiSV   : Distance between two contiguous levels       [m] |
-        !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
-        !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
-        !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
-        !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
-        !C |            dtz_SV   : dt/dz                                      [s/m] |
-        !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
-        !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |
-        !C |            Explic   : Explicit Parameter = 1.0 - Implic                |
-        !C |                                                                        |
-        !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
-        !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
-        !C | #          #SH: Hapex-Sahel Values                                     !
-        !C |                                                                        |
-        !C +------------------------------------------------------------------------+
-        !
-        !
-
-        !C +--Global Variables
-        !C +  ================
-
-        USE dimphy
-        USE VARphy
-        USE VAR_SV
-        USE VARdSV
-        USE VAR0SV
-        USE VARxSV
-        USE VARtSV
-        USE VARxSV
-        USE VARySV
-        IMPLICIT NONE
-
-
-
-        !C +--Arguments
-        !C +  ==================
-        INTEGER, INTENT(IN) :: knon
-
-        !C +--Internal Variables
-        !C +  ==================
-
-        INTEGER :: ivt, ist, ikl, isl, isn, ikh
-        INTEGER :: misl_2, nisl_2
-        REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
-        REAL, PARAMETER :: RHsMin = 0.001        ! Min.Soil Relative Humidity
-        REAL :: PsiMax                        ! Max.Soil Water    Potential
-        REAL :: a_Khyd, b_Khyd                 ! Water conductivity
-
-
-        !c #WR REAL    ::  Khyd_x,Khyd_y
-
-
-
-        !C +--Non Time Dependant SISVAT parameters
-        !C +  ====================================
-
-        !C +--Soil Discretization
-        !C +  -------------------
-
-        !C +--Numerical Scheme Parameters
-        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        Implic = 0.75                           ! 0.5  <==> Crank-Nicholson  
-        Explic = 1.00 - Implic                  !                            
-
-        !C +--Soil/Snow Layers Indices
-        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^
-        DO  isl = -nsol, 0
-            islpSV(isl) = isl + 1
-            islpSV(isl) = min(islpSV(isl), 0)
-            islmSV(isl) = isl - 1
-            islmSV(isl) = max(-nsol, islmSV(isl))
-        END DO
-
-        DO  isn = 1, nsno
-            isnpSV(isn) = isn + 1
-            isnpSV(isn) = min(isnpSV(isn), nsno)
-        END DO
-
-        !C +--Soil      Layers Thicknesses
-        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
-        !c #kd IF (nsol.gt.4)                                              THEN
-        !c #kd   DO isl=-5,-nsol,-1
-        !c #kd     dz_dSV(isl)=   1.
-        !c #kd   END DO
-        !c #kd END IF
-        !
-        !      IF (nsol.ne.4)                                              THEN
-        !        DO isl= 0,-nsol,-1
-        !          misl_2 =     -mod(isl,2)
-        !          nisl_2 =         -isl/2
-        !          dz_dSV(isl)=(((1-misl_2) * 0.001
-        !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.
-        !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm
-        !
-        !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001
-        !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.
-        !
-        !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001
-        !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5
-        !        END DO
-        !          dz_dSV(0)  =               0.001
-        !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004
-        !      END IF
-
-        zz_dSV = 0.
-        DO  isl = -nsol, 0
-            dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
-            dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
-            dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
-            dtz_SV(isl) = dt__SV / dz_dSV(isl)
-            dtz_SV2(isl) = 1. / dz_dSV(isl)
-            dz78SV(isl) = 0.875 * dz_dSV(isl)
-            dz34SV(isl) = 0.750 * dz_dSV(isl)
-            dz_8SV(isl) = 0.125 * dz_dSV(isl)
-            dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl))                        &
-                    & + 0.750 * dz_dSV(isl)                                &
-                    & + 0.125 * dz_dSV(islpSV(isl))
-            zz_dSV = zz_dSV + dz_dSV(isl)
-        END DO
-        DO ikl = 1, knon !v
-            dzsnSV(ikl, 0) = dz_dSV(0)
-        END DO
-
-        !C +--Conversion to a 50 m Swab Ocean Discretization
-        !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        OcndSV = 0.
-        DO isl = -nsol, 0
-            OcndSV = OcndSV + dz_dSV(isl)
-        END DO
-        OcndSV = 50. / OcndSV
-
-
-        !C +--Secondary Soil       Parameters
-        !C +  -------------------------------
-
-        DO  ist = 0, nsot
-            rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6   ! Soil Contrib. to (ro c)_s
-            s1__SV(ist) = bCHdSV(ist)          & ! Factor of (eta)**(b+2)
-                    & * psidSV(ist) * Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)
-                    & / (etadSV(ist)**(bCHdSV(ist) + 3.))     !
-            s2__SV(ist) = Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)
-                    & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.))     !    in DR97, Eqn.(3.35)
-
-            !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
-            !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-            Psimax = -(log(RHsMin)) / 7.2E-5        ! DR97, Eqn 3.15 Inversion
-            etamSV(ist) = etadSV(ist)                                      &
-                    & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
-        END DO
-        etamSV(12) = 0.
-
-        !C +--Piecewise Hydraulic Conductivity Profiles
-        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        DO   ist = 0, nsot
-
-            d__eta = etadSV(ist) / nkhy
-            eta__1 = 0.
-            eta__2 = d__eta
-            DO ikh = 0, nkhy
-                Khyd_1 = s2__SV(ist)             & ! DR97, Eqn.(3.35)
-                        & * (eta__1      **(2. * bCHdSV(ist) + 3.))        !
-                Khyd_2 = s2__SV(ist)             &!
-                        & * (eta__2      **(2. * bCHdSV(ist) + 3.))        !
-
-                a_Khyd = (Khyd_2 - Khyd_1) / d__eta   !
-                b_Khyd = Khyd_1 - a_Khyd * eta__1   !
-                !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !
-                !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !
-                aKdtSV(ist, ikh) = a_Khyd * dt__SV   !
-                bKdtSV(ist, ikh) = b_Khyd * dt__SV   !
-
-                eta__1 = eta__1 + d__eta
-                eta__2 = eta__2 + d__eta
-            END DO
-        END DO
-
-        return
-
-    END SUBROUTINE SISVAT_ini
-
-
-    !***************************************************************************
-
-    SUBROUTINE sisvatetat0 (fichnom, ikl2i)
-
-        USE dimphy
-        USE mod_grid_phy_lmdz
-        USE mod_phys_lmdz_para
-
-        USE iostart
-        USE VAR_SV
-        USE VARdSV
-        USE VARxSV
-        USE VARtSV
-        USE indice_sol_mod
-
-        IMPLICIT none
-        !======================================================================
-        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
-        ! Objet: Lecture du fichier de conditions initiales pour SISVAT
-        !======================================================================
-        include "netcdf.inc"
-        !    include "indicesol.h"
-
-        !    include "dimsoil.h"
-        include "clesphys.h"
-        include "compbl.h"
-
-        !======================================================================
-        CHARACTER(LEN = *) :: fichnom
-
-        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
-        REAL, DIMENSION(klon) :: rlon
-        REAL, DIMENSION(klon) :: rlat
-
-        ! les variables globales ecrites dans le fichier restart
-        REAL, DIMENSION(klon) :: isno
-        REAL, DIMENSION(klon) :: ispi
-        REAL, DIMENSION(klon) :: iice
-        REAL, DIMENSION(klon) :: rusn
-        REAL, DIMENSION(klon, nsno) :: isto
-
-        REAL, DIMENSION(klon, nsismx) :: Tsis
-        REAL, DIMENSION(klon, nsismx) :: eta
-        REAL, DIMENSION(klon, nsismx) :: ro
-
-        REAL, DIMENSION(klon, nsno) :: dzsn
-        REAL, DIMENSION(klon, nsno) :: G1sn
-        REAL, DIMENSION(klon, nsno) :: G2sn
-        REAL, DIMENSION(klon, nsno) :: agsn
-
-        REAL, DIMENSION(klon) :: toic
-
-        INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
-        CHARACTER (len = 2) :: str2
-        LOGICAL :: found
-
-        errT = 0
-        errro = 0
-        erreta = 0
-        errdz = 0
-        snopts = 0
-        ! Ouvrir le fichier contenant l'etat initial:
-
-        CALL open_startphy(fichnom)
-
-        ! Lecture des latitudes, longitudes (coordonnees):
-
-        CALL get_field("latitude", rlat, found)
-        CALL get_field("longitude", rlon, found)
-
-        CALL get_field("n_snows", isno, found)
-        IF (.NOT. found) THEN
-            PRINT*, 'phyetat0: Le champ <n_snows> est absent'
-            PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
-        ENDIF
-
-        CALL get_field("n_ice_top", ispi, found)
-        CALL get_field("n_ice", iice, found)
-        CALL get_field("surf_water", rusn, found)
-
-
-        CALL get_field("to_ice", toic, found)
-        IF (.NOT. found) THEN
-            PRINT*, 'phyetat0: Le champ <to_ice> est absent'
-            toic(:) = 0.
-        ENDIF
-
-        DO isn = 1, nsno
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("AGESNOW" // str2, &
-                        agsn(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsno
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("DZSNOW" // str2, &
-                        dzsn(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsno
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("G2SNOW" // str2, &
-                        G2sn(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsno
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("G1SNOW" // str2, &
-                        G1sn(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsismx
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("ETA" // str2, &
-                        eta(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsismx
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("RO" // str2, &
-                        ro(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsismx
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("TSS" // str2, &
-                        Tsis(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        DO isn = 1, nsno
-            IF (isn.LE.99) THEN
-                WRITE(str2, '(i2.2)') isn
-                CALL get_field("HISTORY" // str2, &
-                        isto(:, isn), found)
-            ELSE
-                PRINT*, "Trop de couches"
-                CALL abort
-            ENDIF
-        ENDDO
-        write(*, *)'Read ', fichnom, ' finished!!'
-
-        !*********************************************************************************
-        ! Compress restart file variables for SISVAT
-
-        DO  ikl = 1, klon
-            i = ikl2i(ikl)
-            IF (i > 0) THEN
-                isnoSV(ikl) = INT(isno(i))          ! Nb Snow/Ice Lay.
-                ispiSV(ikl) = INT(ispi(i))          ! Nb Supr.Ice Lay.
-                iiceSV(ikl) = INT(iice(i))          ! Nb      Ice Lay.
-
-                DO isl = -nsol, 0
-                    ro__SV(ikl, isl) = ro(i, nsno + 1 - isl)       !
-                    eta_SV(ikl, isl) = eta(i, nsno + 1 - isl)         ! Soil Humidity
-                    !hjp 15/10/2010
-                    IF (eta_SV(ikl, isl) <= 1.e-6) THEN          !hj check
-                        eta_SV(ikl, isl) = 1.e-6
-                    ENDIF
-                    TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl)        ! Soil Temperature
-                    IF (TsisSV(ikl, isl) <= 1.) THEN             !hj check
-                        !                errT=errT+1
-                        TsisSV(ikl, isl) = 273.15 - 0.2              ! Etienne: negative temperature since soil is ice
-                    ENDIF
-
-                END DO
-                write(*, *)'Copy histo', ikl
-
-                DO  isn = 1, isnoSV(ikl) !nsno
-                    snopts = snopts + 1
-                    IF (isto(i, isn) > 10.) THEN          !hj check
-                        write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn)
-                        isto(i, isn) = 1.
-                    ENDIF
-
-                    istoSV(ikl, isn) = INT(isto(i, isn))     ! Snow     History
-                    ro__SV(ikl, isn) = ro(i, isn)            !        [kg/m3]
-                    eta_SV(ikl, isn) = eta(i, isn)           !        [m3/m3]
-                    TsisSV(ikl, isn) = Tsis(i, isn)          !            [K]
-
-                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
-                        errT = errT + 1
-                        TsisSV(ikl, isn) = TsisSV(ikl, 0)
-                    ENDIF
-                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
-                        TsisSV(ikl, isn) = 263.15
-                    ENDIF
-                    IF (eta_SV(ikl, isn) < 1.e-9) THEN          !hj check
-                        eta_SV(ikl, isn) = 1.e-6
-                        erreta = erreta + 1
-                    ENDIF
-                    IF (ro__SV(ikl, isn) <= 10.) THEN          !hj check
-                        ro__SV(ikl, isn) = 11.
-                        errro = errro + 1
-                    ENDIF
-                    write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
-                    G1snSV(ikl, isn) = G1sn(i, isn)          ! [-]        [-]
-                    G2snSV(ikl, isn) = G2sn(i, isn)          ! [-] [0.0001 m]
-                    dzsnSV(ikl, isn) = dzsn(i, isn)          !            [m]
-                    agsnSV(ikl, isn) = agsn(i, isn)          !          [day]
-                END DO
-                rusnSV(ikl) = rusn(i)              ! Surficial Water
-                toicSV(ikl) = toic(i)              ! bilan snow to ice
-            END IF
-        END DO
-
-    END SUBROUTINE sisvatetat0
-
-
-    !======================================================================
-    SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)
-
-
-
-        !======================================================================
-        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
-        ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
-        !======================================================================
-        USE mod_grid_phy_lmdz
-        USE mod_phys_lmdz_para
-        USE iostart
-        USE VAR_SV
-        USE VARxSV
-        USE VARySV !hj tmp 12 03 2010
-        USE VARtSV
-        USE indice_sol_mod
-        USE dimphy
-
-        IMPLICIT none
-
-        include "netcdf.inc"
-        include "clesphys.h"
-        include "compbl.h"
-
-        !======================================================================
-
-        CHARACTER(LEN = *) :: fichnom
-        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
-        REAL, DIMENSION(klon), INTENT(IN) :: rlon
-        REAL, DIMENSION(klon), INTENT(IN) :: rlat
-
-        ! les variables globales ecrites dans le fichier restart
-        REAL, DIMENSION(klon) :: isno
-        REAL, DIMENSION(klon) :: ispi
-        REAL, DIMENSION(klon) :: iice
-        REAL, DIMENSION(klon, nsnowmx) :: isto
-
-        REAL, DIMENSION(klon, nsismx) :: Tsis
-        REAL, DIMENSION(klon, nsismx) :: eta
-        REAL, DIMENSION(klon, nsnowmx) :: dzsn
-        REAL, DIMENSION(klon, nsismx) :: ro
-        REAL, DIMENSION(klon, nsnowmx) :: G1sn
-        REAL, DIMENSION(klon, nsnowmx) :: G2sn
-        REAL, DIMENSION(klon, nsnowmx) :: agsn
-        REAL, DIMENSION(klon) :: IRs
-        REAL, DIMENSION(klon) :: LMO
-        REAL, DIMENSION(klon) :: rusn
-        REAL, DIMENSION(klon) :: toic
-        REAL, DIMENSION(klon) :: Bufs
-        REAL, DIMENSION(klon) :: alb1, alb2, alb3
-
-        INTEGER isl, ikl, i, isn, ierr
-        CHARACTER (len = 2) :: str2
-        INTEGER :: pass
-
-        isno(:) = 0
-        ispi(:) = 0
-        iice(:) = 0
-        IRs(:) = 0.
-        LMO(:) = 0.
-        eta(:, :) = 0.
-        Tsis(:, :) = 0.
-        isto(:, :) = 0
-        ro(:, :) = 0.
-        G1sn(:, :) = 0.
-        G2sn(:, :) = 0.
-        dzsn(:, :) = 0.
-        agsn(:, :) = 0.
-        rusn(:) = 0.
-        toic(:) = 0.
-        Bufs(:) = 0.
-        alb1(:) = 0.
-        alb2(:) = 0.
-        alb3(:) = 0.
-
-        !***************************************************************************
-        ! Uncompress SISVAT output variables for storage
-
-        DO  ikl = 1, klon
-            i = ikl2i(ikl)
-            IF (i > 0) THEN
-                isno(i) = 1. * isnoSV(ikl)               ! Nb Snow/Ice Lay.
-                ispi(i) = 1. * ispiSV(ikl)               ! Nb Supr.Ice Lay.
-                iice(i) = 1. * iiceSV(ikl)               ! Nb      Ice Lay.
-
-                !        IRs(i)        = IRs_SV(ikl)
-                !        LMO(i)        = LMO_SV(ikl)
-
-                DO isl = -nsol, 0                           !
-                    eta(i, nsno + 1 - isl) = eta_SV(ikl, isl)            ! Soil Humidity
-                    Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl)            ! Soil Temperature
-                    ro(i, nsno + 1 - isl) = ro__SV(ikl, isl)            !        [kg/m3]
-                END DO
-
-                DO  isn = 1, nsno
-                    isto(i, isn) = 1. * istoSV(ikl, isn)         ! Snow     History
-                    ro(i, isn) = ro__SV(ikl, isn)            !        [kg/m3]
-                    eta(i, isn) = eta_SV(ikl, isn)            !        [m3/m3]
-                    Tsis(i, isn) = TsisSV(ikl, isn)            !            [K]
-                    G1sn(i, isn) = G1snSV(ikl, isn)            ! [-]        [-]
-                    G2sn(i, isn) = G2snSV(ikl, isn)            ! [-] [0.0001 m]
-                    dzsn(i, isn) = dzsnSV(ikl, isn)            !            [m]
-                    agsn(i, isn) = agsnSV(ikl, isn)            !          [day]
-                END DO
-                rusn(i) = rusnSV(ikl)                  ! Surficial Water
-                toic(i) = toicSV(ikl)                  ! to ice
-                alb1(i) = alb1sv(ikl)
-                alb2(i) = alb2sv(ikl)
-                alb3(i) = alb3sv(ikl)
-                !        Bufs(i)       = BufsSV(ikl)
-            END IF
-        END DO
-
-        CALL open_restartphy(fichnom)
-
-        DO pass = 1, 2
-            CALL put_field(pass, "longitude", &
-                    "Longitudes de la grille physique", rlon)
-            CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
-
-            CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
-            CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
-            CALL put_field(pass, "n_ice", "number of ice layers", iice)
-            CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
-            CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
-            CALL put_field(pass, "surf_water", "Surficial water", rusn)
-            CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
-            CALL put_field(pass, "alb_1", "albedo sw", alb1)
-            CALL put_field(pass, "alb_2", "albedo nIR", alb2)
-            CALL put_field(pass, "alb_3", "albedo fIR", alb3)
-            CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
-
-            DO isn = 1, nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "AGESNOW" // str2, &
-                            "Age de la neige layer No." // str2, &
-                            agsn(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "DZSNOW" // str2, &
-                            "Snow/ice thickness layer No." // str2, &
-                            dzsn(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "G2SNOW" // str2, &
-                            "Snow Property 2, layer No." // str2, &
-                            G2sn(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "G1SNOW" // str2, &
-                            "Snow Property 1, layer No." // str2, &
-                            G1sn(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsismx
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "ETA" // str2, &
-                            "Soil/snow water content layer No." // str2, &
-                            eta(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsismx   !nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "RO" // str2, &
-                            "Snow density layer No." // str2, &
-                            ro(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsismx
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "TSS" // str2, &
-                            "Soil/snow temperature layer No." // str2, &
-                            Tsis(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-            DO isn = 1, nsno
-                IF (isn.LE.99) THEN
-                    WRITE(str2, '(i2.2)') isn
-                    CALL put_field(pass, "HISTORY" // str2, &
-                            "Snow history layer No." // str2, &
-                            isto(:, isn))
-                ELSE
-                    PRINT*, "Trop de couches"
-                    CALL abort
-                ENDIF
-            ENDDO
-
-            CALL enddef_restartphy
-        ENDDO
-        CALL close_restartphy
-
-    END SUBROUTINE sisvatredem
-
-END MODULE surf_inlandsis_mod
Index: LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.f90
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.f90	(revision 5268)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/surf_inlandsis_mod.f90	(revision 5268)
@@ -0,0 +1,1495 @@
+MODULE surf_inlandsis_mod
+
+    IMPLICIT NONE
+
+CONTAINS
+
+
+SUBROUTINE surf_inlandsis(knon, rlon, rlat, ikl2i, itime, dtime, debut, lafin, &
+            rmu0, swdown, lwdown, albedo_old, pexner, ps, p1lay, &
+            precip_rain, precip_snow, &
+            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
+            rugos, snow_cont_air, alb_soil, alt, slope, cloudf, &
+            radsol, qsol, tsoil, snow, zfra, snowhgt, qsnow, to_ice, sissnow, agesno, &
+            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
+            runoff_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat, dflux_s,dflux_l, &
+            tsurf_new, alb1, alb2, alb3, alb6, emis_new, z0m, z0h, qsurf)
+
+        ! |                                                                        |
+        ! |   SubRoutine surf_inlandsis: Interfacing Lmdz AND Sisvat's Ice and Snow|
+        ! |                              (INLANDSIS)                               |
+        ! |   SISVAT (Soil/Ice Snow Vegetation Atmosphere Transfer Scheme)         |
+        ! |   surface scheme of the Modele Atmospherique Regional (MAR)            |
+        ! |   Author: Heinz Juergen Punge, LSCE                June 2009           |
+        ! |     based on the MAR-SISVAT interface by Hubert Gallee                 |
+        ! |   Updated by Etienne Vignon, Cecile Agosta                             |
+        ! |                                                                        |
+        ! +------------------------------------------------------------------------+
+        ! |
+        ! |   In the current setup, SISVAT is used only to model the land ice      |
+        ! |   part of the surface; hence it is called with the compressed variables|
+        ! |   from pbl_surface, and only by the surf_landice routine.              |
+        ! |                                                                        |
+        ! |   In this interface it is assumed that the partitioning of the soil,   |
+        ! |   and hence the number of grid points is constant during a simulation, |
+        ! |   hence eg. snow properties remain stored in the global SISVAT         |
+        ! |   variables between the calls and don't need to be handed over as      |
+        ! |   arguments. When the partitioning is supposed to change, make sure to |
+        ! |   update the variables.                                                |
+        ! |                                                                        |
+        ! |   INPUT    (via MODULES VARxSV, VARySV, VARtSV ...)                    |
+        ! |   ^^^^^     xxxxSV: SISVAT/LMDZ interfacing variables                  |
+        ! |                                                                        |
+        ! +------------------------------------------------------------------------+
+
+        USE dimphy
+        USE VAR_SV
+        USE VARdSV
+        USE VARxSV
+        USE VARySV
+        USE VARtSV
+        USE VARphy
+        USE surface_data, only : iflag_tsurf_inlandsis, SnoMod, BloMod, ok_outfor
+
+        IMPLICIT NONE
+
+        ! +--INTERFACE Variables
+        ! +  ===================
+        !    include  "dimsoil.h"
+
+        ! +--Global Variables
+        ! +  ================
+        ! Input Variables for SISVAT
+        INTEGER, INTENT(IN) :: knon
+        INTEGER, INTENT(IN) :: itime
+        REAL, INTENT(IN) :: dtime
+        LOGICAL, INTENT(IN) :: debut     ! true if first step
+        LOGICAL, INTENT(IN) :: lafin     ! true if last step
+
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i     ! Index Decompression
+        REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat
+        REAL, DIMENSION(klon), INTENT(IN) :: rmu0      ! cos sol. zenith angle
+        REAL, DIMENSION(klon), INTENT(IN) :: swdown    !
+        REAL, DIMENSION(klon), INTENT(IN) :: lwdown    !
+        REAL, DIMENSION(klon), INTENT(IN) :: albedo_old
+        REAL, DIMENSION(klon), INTENT(IN) :: pexner    ! Exner potential
+        REAL, DIMENSION(klon), INTENT(IN) :: precip_rain, precip_snow
+        REAL, DIMENSION(klon), INTENT(IN) :: zsl_height, wind_velo
+        REAL, DIMENSION(klon), INTENT(IN) :: temp_air, spechum, ps, p1lay
+        REAL, DIMENSION(klon), INTENT(IN) :: dens_air, tsurf
+        REAL, DIMENSION(klon), INTENT(IN) :: rugos
+        REAL, DIMENSION(klon), INTENT(IN) :: snow_cont_air
+        REAL, DIMENSION(klon), INTENT(IN) :: alb_soil, slope
+        REAL, DIMENSION(klon), INTENT(IN) :: alt       ! surface elevation
+        REAL, DIMENSION(klon), INTENT(IN) :: cloudf
+        REAL, DIMENSION(klon), INTENT(IN) :: AcoefH, AcoefQ
+        REAL, DIMENSION(klon), INTENT(IN) :: BcoefH, BcoefQ
+        REAL, DIMENSION(klon), INTENT(IN) :: cdragm, cdragh
+        REAL, DIMENSION(klon), INTENT(IN) :: ustar   ! friction velocity
+
+        ! Variables exchanged between LMDZ and SISVAT
+        REAL, DIMENSION(klon), INTENT(IN) :: radsol    ! Surface absorbed rad.
+        REAL, DIMENSION(klon), INTENT(INOUT) :: snow      ! Tot snow mass [kg/m2]
+        REAL, DIMENSION(klon), INTENT(INOUT) :: zfra      ! snwo surface fraction [0-1]
+        REAL, DIMENSION(klon, nsoilmx), INTENT(OUT) :: tsoil ! Soil Temperature
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsol      ! Soil Water Content
+        REAL, DIMENSION(klon), INTENT(INOUT) :: z0m    ! Momentum Roughn Lgt
+        REAL, DIMENSION(klon), INTENT(INOUT) :: z0h    ! Momentum Roughn Lgt
+
+        ! Output Variables for LMDZ
+        REAL, DIMENSION(klon), INTENT(OUT) :: alb1      ! Albedo SW
+        REAL, DIMENSION(klon), INTENT(OUT) :: alb2, alb3 ! Albedo NIR and LW
+        REAL, DIMENSION(klon,6), INTENT(OUT) :: alb6 ! 6 band Albedo
+        REAL, DIMENSION(klon), INTENT(OUT) :: emis_new  ! Surface Emissivity
+        REAL, DIMENSION(klon), INTENT(OUT) :: runoff_lic ! Runoff
+        REAL, DIMENSION(klon), INTENT(OUT) :: ffonte    ! enthalpy flux due to surface melting
+        REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte   ! water flux due to surface melting
+        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_s   ! d/dT sens. ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: dflux_l   ! d/dT latent ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: fluxsens  ! Sensible ht flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: fluxlat   ! Latent heat flux
+        REAL, DIMENSION(klon), INTENT(OUT) :: evap      ! Evaporation
+        REAL, DIMENSION(klon), INTENT(OUT) :: erod      ! Erosion of surface snow (flux)
+        REAL, DIMENSION(klon), INTENT(OUT) :: agesno    ! Snow age (top layer)
+        REAL, DIMENSION(klon), INTENT(OUT) :: tsurf_new ! Surface Temperature
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsurf     ! Surface Humidity
+
+        ! Specific INLANDIS outputs
+        REAL, DIMENSION(klon), INTENT(OUT) :: qsnow     ! Total H2O snow[kg/m2]
+        REAL, DIMENSION(klon), INTENT(OUT) :: snowhgt   ! Snow height (m)
+        REAL, DIMENSION(klon), INTENT(OUT) :: to_ice    ! Snow passed to ice
+        REAL, DIMENSION(klon), INTENT(OUT) :: sissnow   ! Snow in model (kg/m2)
+
+        ! +--Internal  Variables
+        ! +  ===================
+
+        CHARACTER(len = 20) :: fn_outfor ! Name for output file
+        CHARACTER (len = 80)              :: abort_message
+        CHARACTER (len = 20)              :: modname = 'surf_inlandsis_mod'
+
+        INTEGER :: i, ig, ikl, isl, isn, nt
+        INTEGER :: gp_outfor, un_outfor
+        REAL, PARAMETER :: f1 = 0.5
+        REAL, PARAMETER :: sn_upp = 10000., sn_low = 500.
+        REAL, PARAMETER :: sn_add = 400., sn_div = 2.
+        ! snow mass upper,lower limit,
+        ! added mass/division lowest layer
+        REAL, PARAMETER :: c1_zuo = 12.960e+4, c2_zuo = 2.160e+6
+        REAL, PARAMETER :: c3_zuo = 1.400e+2, czemin = 1.e-3
+        ! Parameters for drainage
+        ! c1_zuo/ 2.796e+4/,c2_zuo/2.160e+6/,c3_zuo/1.400e+2/ !     Tuning
+        ! +...        Run Off Parameters
+        ! +           86400*1.5 day     ...*25 days (Modif. ETH Camp: 86400*0.3day)
+        ! +           (Zuo and Oerlemans 1996, J.Glacio. 42, 305--317)
+
+        REAL, DIMENSION(klon) :: eps0SL          ! surface Emissivity
+        REAL :: zsigma, Ua_min, Us_min, lati
+        REAL, PARAMETER :: cdmax=0.05
+        REAL :: lambda          ! Par. soil discret.
+        REAL, DIMENSION(nsoilmx), SAVE :: dz1, dz2         ! Soil layer thicknesses
+        !$OMP THREADPRIVATE(dz1,dz2)
+        LOGICAL, SAVE :: firstcall
+        !$OMP THREADPRIVATE(firstcall)
+
+        INTEGER :: iso
+        LOGICAL :: file_exists
+        CHARACTER(len = 20) :: fichnom
+        LOGICAL :: is_init_domec
+        ! CA initialization
+        ! dz_profil_15 : 1 m in 15 layers [m]
+        real, parameter :: dz_profil_15(15) = (/0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, &
+                                                0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.17/)
+        ! mean_temp : mean annual surface temperature [K]
+        real, dimension(klon) :: mean_temp
+        ! mean_dens : mean surface density [kg/m3]
+        real, dimension(klon) :: mean_dens
+        ! lat_scale : temperature lapse rate against latitude [K degree-1]
+        real :: lat_scale
+        ! sh_scale : temperature lapse rate against altitude [K km-1]
+        real :: sh_scale
+        ! variables for density profile
+        ! E0, E1 : exponent
+        real :: E0, E1
+        ! depth at which 550 kg m-3 is reached [m]
+        real :: z550
+        ! depths of snow layers
+        real :: depth, snow_depth, distup
+        ! number of initial snow layers
+        integer :: nb_snow_layer
+        ! For density calc.
+        real :: alpha0, alpha1, ln_smb
+        ! theoritical densities [kg m-3]
+        real :: rho0, rho1, rho1_550
+        ! constants for density profile
+        ! C0, C1 : constant, 0.07 for z <= 550 kg m-3
+        real, parameter :: C0 = 0.07
+        real, parameter :: C1 = 0.03
+        ! rho_i : ice density [kg m-3]
+        real, parameter :: rho_ice = 917.
+        ! E_c : activation energy [J mol-1]
+        real, parameter :: E_c = 60000.
+        ! E_g : activation energy [J mol-1]
+        real, parameter :: E_g = 42400.
+        ! R : gas constant [J mol-1 K-1]
+        real, parameter :: R = 8.3144621
+
+      
+     
+
+
+        ! + PROGRAM START
+        ! + -----------------------------------------
+
+        zsigma = 1000.
+        dt__SV = dtime
+
+        IF (debut) THEN
+            firstcall = .TRUE.
+            INI_SV = .false.
+        ELSE
+            firstcall = .false.
+            INI_SV = .true.
+        END IF
+
+        IF (ok_outfor) THEN
+            un_outfor = 51    ! unit number for point output file
+            gp_outfor = 1    ! grid point number for point output 1 for 1D, 273 for zoom-nudg DC
+            fn_outfor = 'outfor_SV.dat'
+        END IF
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + INITIALISATION: BEGIN +++
+        ! + -----------------------------------------
+        IF (firstcall) THEN
+
+            ! +--Array size
+            ! +  -----------------------
+
+            klonv = klon
+            knonv = knon
+                write(*, *) 'ikl, lon and lat in INLANDSIS'
+
+            DO ikl = 1, knon
+                i=ikl2i(ikl)
+                write(*, *) 'ikl=', ikl, 'rlon=', rlon(i), 'rlat=', rlat(i)
+            END DO
+
+            ! +--Variables initizialisation
+            ! +  ---------------------------
+
+            CALL INIT_VARtSV
+            CALL INIT_VARxSV
+            CALL INIT_VARySV
+
+
+
+            ! +--Surface Fall Line Slope
+            ! +  -----------------------
+            IF (SnoMod)  THEN
+                DO ikl = 1, knon
+                    slopSV(ikl) = slope(ikl)
+                    SWf_SV(ikl) = &   ! Normalized Decay of the
+                            exp(-dt__SV             &   ! Surficial Water Content
+                                    / (c1_zuo                &   !(Zuo and Oerlemans 1996,
+                                            + c2_zuo * exp(-c3_zuo * abs(slopSV(ikl)))))  ! J.Glacio. 42, 305--317)
+                END DO
+            END IF
+
+
+
+            ! +--Soil layer thickness . Compute soil discretization (as for LMDZ)
+            ! +  ----------------------------------------------------------------
+            !        write(*,'(/a)') 'Start SISVAT init: soil discretization ', nsoilmx
+            CALL get_soil_levels(dz1, dz2, lambda)
+
+            lambSV = lambda
+            dz1_SV(1:knon, 1:) = 0.
+            dz2_SV(1:knon, 1:) = 0.
+
+            DO isl = -nsol, 0
+                dz_dSV(isl) = 0.5e-3 * dz2(1 - isl)           ! Soil layer thickness
+                DO ikl = 1, knon
+                    dz1_SV(ikl, isl) = dz1(1 - isl)    !1.e-3*
+                    dz2_SV(ikl, isl) = dz2(1 - isl)    !1.e-3*
+                END DO
+            END DO
+
+
+            ! Set variables
+            ! =============
+            DO ikl = 1, knon
+                ! LSmask : Land/Sea Mask
+                LSmask(ikl) = 1
+                ! isotSV : Soil Type -> 12 = ice
+                isotSV(ikl) = 12
+                ! iWaFSV : Soil Drainage (1,0)=(y,n)
+                iWaFSV(ikl) = 1
+                ! eps0SL : Surface Emissivity
+                eps0SL(ikl) = 1.
+                ! alb0SV : Soil Albedo
+                alb0SV(ikl) = alb_soil(ikl)
+                ! Tsf_SV : Surface Temperature, must be bellow freezing
+                Tsf_SV(ikl) = min(temp_air(ikl), TfSnow)
+            END DO
+
+            ! +--Initialization of soil and snow variables in case startsis is not read
+            ! +  ----------------------------------------------------------------------
+
+            is_init_domec=.FALSE.
+
+
+            IF (is_init_domec) THEN
+            ! Coarse initilization inspired from vertcical profiles at Dome C,
+            ! Antarctic Plateaui (10m of snow, 19 levels)
+
+                 DO ikl = 1,knon
+! + Soil
+                 DO isl =   -nsol,0    
+                   TsisSV(ikl,isl) = min(tsoil(ikl,1+nsol),TfSnow-0.2)   !temp_air(ikl)
+                   !tsoil(ikl,1-isl)   Soil Temperature
+                   !TsisSV(ikl,isl) = min(temp_air(ikl),TfSnow-0.2)
+                   eta_SV(ikl,isl) = epsi           !etasoil(ikl,1-isl) Soil Water[m3/m3]
+                   ro__SV(ikl,isl) = rhoIce         !rosoil(ikl,1-isl) volumic mass
+                 END DO   
+
+
+           ! Snow
+                 isnoSV(ikl) = 19
+                 istoSV(ikl, 1:isnoSV(ikl)) = 100
+                 ro__SV(ikl, 1:isnoSV(ikl)) = 350.
+                 eta_SV(ikl, 1:isnoSV(ikl)) = epsi
+                 TsisSV(ikl, 1:isnoSV(ikl)) = min(tsoil(ikl, 1), TfSnow - 0.2)
+                 G1snSV(ikl, 1:isnoSV(ikl)) = 99.
+                 G2snSV(ikl, 1:isnoSV(ikl)) = 2.
+                 agsnSV(ikl, 1:isnoSV(ikl)) = 50.
+                 dzsnSV(ikl, 19) = 0.015
+                 dzsnSV(ikl, 18) = 0.015
+                 dzsnSV(ikl, 17) = 0.020
+                 dzsnSV(ikl, 16) = 0.030
+                 dzsnSV(ikl, 15) = 0.040
+                 dzsnSV(ikl, 14) = 0.060
+                 dzsnSV(ikl, 13) = 0.080
+                 dzsnSV(ikl, 12) = 0.110
+                 dzsnSV(ikl, 11) = 0.150
+                 dzsnSV(ikl, 10) = 0.200
+                 dzsnSV(ikl, 9) = 0.300
+                 dzsnSV(ikl, 8) = 0.420
+                 dzsnSV(ikl, 7) = 0.780
+                 dzsnSV(ikl, 6) = 1.020
+                 dzsnSV(ikl, 5) = 0.980
+                 dzsnSV(ikl, 4) = 1.020
+                 dzsnSV(ikl, 3) = 3.970
+                 dzsnSV(ikl, 2) = 1.020
+                 dzsnSV(ikl, 1) = 1.020
+
+                 END DO
+            ELSE
+
+            ! Initilialisation with climatological temperature and density
+            ! profiles as in MAR. Methodology developed by Cecile Agosta
+  
+            ! initialize with 0., for unused snow layers
+            dzsnSV = 0.
+            G1snSV = 0.
+            G2snSV = 0.
+            istoSV = 0
+            TsisSV = 0.
+
+
+            ! initialize mean variables (unrealistic)
+            mean_temp = TfSnow
+            mean_dens = 300.
+            ! loop on grid cells
+            DO ikl = 1, knon
+            	lati=rlat(ikl2i(ikl))
+                ! approximations for mean_temp and mean_dens
+                ! from Feulner et al., 2013 (DOI: 10.1175/JCLI-D-12-00636.1)
+                ! Fig. 3 and 5 : the lapse rate vs. latitude at high latitude is about 0.55 °C °lat-1
+                ! with a moist-adiabatic lapse rate of 5 °C km-1 everywhere except for Antarctica,
+                ! for Antarctica, a dry-adiabatic lapse rate of 9.8 °C km-1 is assumed.
+                if (lati > 60.) then
+                    ! CA todo : add longitude bounds
+                    ! Greenland mean temperature : function of altitude and latitude
+                    ! for altitudes 0. to 1000. m, lat_scale varies from 0.9 to 0.75 °C °lat-1
+                    lat_scale = (0.75 - 0.9) / 1000. * alt(ikl) + 0.9
+                    lat_scale = max(min(lat_scale, 0.9), 0.75)
+                    ! sh_scale equals the environmental lapse rate : 6.5 °C km-1
+                    sh_scale = 6.5
+                    mean_temp(ikl) = TfSnow + 1.5 - sh_scale * alt(ikl) / 1000. - lat_scale * (lati - 60.)
+                    ! surface density: Fausto et al. 2018, https://doi.org/10.3389/feart.2018.00051
+                    mean_dens(ikl) = 315.
+                else if (lati < -60.) then
+                    ! Antarctica mean temperature : function of altitude and latitude
+                    ! for altitudes 0. to 500. m, lat_scale varies from 1.3 to 0.6 °C °lat-1
+                    lat_scale = (0.6 - 1.3) / 500. * alt(ikl) + 1.3
+                    lat_scale = max(min(lat_scale, 1.3), 0.6)
+                    ! for altitudes 0. to 500. m, sh_scale varies from 6.5 to 9.8 °C km-1
+                    sh_scale = (9.8 - 6.5) / 500. * alt(ikl) + 6.5
+                    sh_scale = max(min(sh_scale, 9.8), 6.5)
+                    mean_temp(ikl) = TfSnow - 7. - sh_scale * alt(ikl) / 1000. + lat_scale * (lati + 60.)
+                    ! Antarctica surface density : function of mean annual temperature
+                    ! surface density of 350. kg m-3 at Dome C and 450. kg m-3 at Prud'homme (Agosta et al. 2013)
+                    ! 350 kg m-3 is a typical value for the Antarctic plateau around 3200 m.
+                    ! Weinhart et al 2020  https://doi.org/10.5194/tc-14-3663-2020 and Sugiyama et al. 2011 oi: 10.3189/2012JoG11J201
+                    ! 320 kg m-3 is reached at Dome A, 4100 m a.s.l.
+                    ! Dome C : st_ant_param(3233, -75.1) = -47.7
+                    ! Dumont d'Urville : st_ant_param(0, -66.66) = -15.7
+                    mean_dens(ikl) =  (450. - 320.) / (-15.7 + 47.7) * (mean_temp(ikl) - TfSnow + 15.7) + 450.
+                    mean_dens(ikl) = min(450., max(320., mean_dens(ikl)))
+                else
+
+                !    write(*, *) 'Attention: temperature initialization is only defined for Greenland and Antarctica'
+
+                     mean_dens(ikl) =350.
+                     mean_temp(ikl) = min(tsoil(ikl,1),TfSnow-0.2)
+
+                !abort_message='temperature initialization is only defined for Greenland and Antarctica'
+                !CALL abort_physic(modname,abort_message,1)
+
+                end if
+ 
+                ! mean_temp is defined for ice ground only
+                mean_temp(ikl) = min(mean_temp(ikl), TfSnow - 0.2)
+
+                ! Soil layers
+                ! ===========
+                DO isl = -nsol, 0
+                    ! TsisSV : Temperature [K]
+                    TsisSV(ikl, isl) = mean_temp(ikl)
+                    ! eta_SV : Soil Water [m3/m3]
+                    eta_SV(ikl, isl) = epsi
+                    ! ro__SV : Volumic Mass [kg/m3]
+                    ro__SV(ikl, isl) = rhoIce
+                END DO
+
+                ! Snow layers
+                ! ===========
+                ! snow_depth : initial snow depth
+                snow_depth = 20.
+                ! nb_snow_layer : initial nb of snow layers
+                nb_snow_layer = 15
+                ! isnoSV : total nb of snow layers
+                isnoSV(ikl) = nb_snow_layer
+                ! depth : depth of each layer
+                depth = snow_depth
+                do isl = 1, nb_snow_layer
+                    ! dzsnSV : snow layer thickness
+                    dzsnSV(ikl, isl) = max(0.01, snow_depth * dz_profil_15(nb_snow_layer - isl + 1))
+                    ! G1snSV : dendricity (<0) or sphericity (>0) : 99. = sperical
+                    G1snSV(ikl, isl) = 99.
+                    ! G2snSV : Sphericity (>0) or Size [1/10 mm] : 2. = small grain size
+                    G2snSV(ikl, isl) = 3.
+                    agsnSV(ikl, isl) = 0.
+                    istoSV(ikl, isl) = 0
+                    ! eta_SV : Liquid Water Content [m3/m3]
+                    eta_SV(ikl, isl) = 0.
+                    ! distance to surface
+                    depth = depth - dzsnSV(ikl,isl) / 2.
+                    distup = min(1., max(0., depth / snow_depth))
+                    ! TsisSV : Temperature [K], square interpolation between Tsf_SV (surface) and mean_temp (bottom)
+                    TsisSV(ikl, isl) = Tsf_SV(ikl) * (1. - distup**2) + mean_temp(ikl) * distup**2
+                    ! firn density : densification formulas from :
+                    ! Ligtenberg et al 2011 eq. (6) (www.the-cryosphere.net/5/809/2011/)
+                    ! equivalent to Arthern et al. 2010 eq. (4) "Nabarro-Herring" (doi:10.1029/2009JF001306)
+                    ! Integration of the steady state equation
+                    ! ln_smb approximated as a function of temperature
+                    ln_smb = max((mean_temp(ikl) - TfSnow) * 5. / 60. + 8., 3.)
+                    ! alpha0, alpha1 : correction coefficient as a function of ln_SMB from Ligtenberg 2011, adjusted for alpha1
+                    alpha0 = max(1.435 - 0.151 * ln_smb, 0.25)
+                    alpha1 = max(2.0111 - 0.2051 * ln_smb, 0.25)
+                    E0 = C0 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha0
+                    E1 = C1 * gravit * exp((E_g - E_c)/(R * mean_temp(ikl))) * rho_ice * alpha1
+                    z550 = log((rho_ice/mean_dens(ikl) - 1.)/(rho_ice/550. - 1.)) / E0
+                    rho0 = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
+                    rho1 = exp(E1 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E1 * depth)) * rho_ice
+                    if (depth <= z550) then
+                        ro__SV(ikl, isl) = exp(E0 * depth) / (rho_ice / mean_dens(ikl) - 1 + exp(E0 * depth)) * rho_ice
+                    else
+                        ro__SV(ikl, isl) = exp(E1 * (depth - z550)) / (rho_ice / 550. - 1 + exp(E1 * (depth - z550))) * rho_ice
+                    end if
+                    depth = depth - dzsnSV(ikl,isl) / 2.
+                    
+                end do
+
+            END DO
+
+            END IF
+
+
+            ! + Numerics paramaters, SISVAT_ini 
+            ! +  ----------------------
+            CALL SISVAT_ini(knon)
+
+
+            ! +--Read restart file
+            ! +  =================================================
+
+            INQUIRE(FILE = "startsis.nc", EXIST = file_exists)
+            IF (file_exists) THEN
+                CALL sisvatetat0("startsis.nc", ikl2i)
+            END IF
+
+
+
+            ! +--Output ascii file
+            ! +  =================================================
+
+            ! open output file
+            IF (ok_outfor) THEN
+                open(unit = un_outfor, status = 'replace', file = fn_outfor)
+                ikl = gp_outfor     ! index sur la grille land ice
+                write(un_outfor, *) fn_outfor, ikl, dt__SV, rlon(ikl2i(ikl)), rlat(ikl2i(ikl))
+                write(un_outfor, *) 'nsnow - albedo - z0m - z0h , dz [m,30], temp [K,41], rho [kg/m3,41], eta [kg/kg,41] &
+                        & G1 [-,30], G2 [-,30], agesnow [d,30], history [-,30], DOP [m,30]'
+            END IF
+
+        END IF  ! firstcall
+        ! +
+        ! +  +++  INITIALISATION:  END  +++
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + READ FORCINGS
+        ! + ------------------------
+
+        ! + Update Forcings for SISVAT given by the LMDZ model.
+        ! +
+        DO ikl = 1, knon
+
+            ! +--Atmospheric Forcing                                    (INPUT)
+            ! +  ^^^^^^^^^^^^^^^^^^^                                     ^^^^^
+            za__SV(ikl) = zsl_height(ikl)               ! surface layer height (fisr model level) [m]
+            Ua_min = 0.2 * sqrt(za__SV(ikl))            !
+            VV__SV(ikl) = max(Ua_min, wind_velo(ikl))   ! Wind velocity       [m/s]
+            TaT_SV(ikl) = temp_air(ikl)                 ! BL top Temperature    [K]
+            ExnrSV(ikl) = pexner(ikl)                   ! Exner potential
+            rhT_SV(ikl) = dens_air(ikl)                 ! Air density
+            QaT_SV(ikl) = spechum(ikl)                  ! Specific humidity
+            ps__SV(ikl) = ps(ikl)                       ! surface pressure     [Pa]
+            p1l_SV(ikl) = p1lay(ikl)                    ! lowest atm. layer press[Pa]
+
+            ! +--Surface properties
+            ! +  ^^^^^^^^^^^^^^^^^^
+
+            Z0m_SV(ikl) = z0m(ikl)                      ! Moment.Roughn.L.
+            Z0h_SV(ikl) = z0h(ikl)                      ! Moment.Roughn.L.
+
+            ! +--Energy Fluxes                                          (INPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^
+            coszSV(ikl) = max(czemin, rmu0(ikl))         ! cos(zenith.Dist.)
+            sol_SV(ikl) = swdown(ikl)                   ! downward Solar
+            IRd_SV(ikl) = lwdown(ikl)                   ! downward IR
+            rsolSV(ikl) = radsol(ikl)                   ! surface absorbed rad.
+
+            ! +--Water  Fluxes                                          (INPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^
+            drr_SV(ikl) = precip_rain(ikl)              ! Rain fall rate  [kg/m2/s]
+            dsn_SV(ikl) = precip_snow(ikl)              ! Snow fall rate  [kg/m2/s]
+
+            ! #BS    dbs_SV(ikl) = blowSN(i,j,n)
+            ! dbs_SV = Maximum potential erosion amount [kg/m2]
+            ! => Upper bound for eroded snow mass
+            !        uss_SV(ikl) = SLussl(i,j,n) ! u*qs* (only for Tv in sisvatesbl.f)
+            ! #BS  if(dsn_SV(ikl)>eps12.and.erprev(i,j,n).gt.eps9) then
+            ! #BS    dsnbSV(ikl) =1.0-min(qsHY(i,j,kB)     !BS neglib. at kb ~100 magl)
+            ! #BS.                        /max(qshy(i,j,mz),eps9),unun)
+            ! #BS    dsnbSV(ikl) = max(dsnbSV(ikl),erprev(i,j,n)/dsn_SV(ikl))
+            ! #BS    dsnbSV(ikl) = max(0.,min(1.,dsnbSV(ikl)))
+            ! #BS  else
+            ! #BS    dsnbSV(ikl) = 0.
+            ! #BS  endif
+            !      dsnbSV is the drift fraction of deposited snow updated in sisvat.f
+            !      will be used for characterizing the Buffer Layer
+            !      (see update of  Bros_N, G1same, G2same, zroOLD, zroNEW)
+            ! #BS  if(n==1) qbs_HY(i,j) = dsnbSV(ikl)
+            qsnoSV(ikl) = snow_cont_air(ikl)
+
+
+
+            ! +--Soil/BL                                      (INPUT)
+            ! +  ^^^^^^^                                       ^^^^^
+            alb0SV(ikl) = alb_soil(ikl)                 ! Soil background Albedo
+            AcoHSV(ikl) = AcoefH(ikl)
+            BcoHSV(ikl) = BcoefH(ikl)
+            AcoQSV(ikl) = AcoefQ(ikl)
+            BcoQSV(ikl) = BcoefQ(ikl)
+            cdH_SV(ikl) = min(cdragh(ikl),cdmax)
+            cdM_SV(ikl) = min(cdragm(ikl),cdmax)
+            rcdmSV(ikl) = sqrt(cdM_SV(ikl))
+            Us_min = 0.01
+            us__SV(ikl) = max(Us_min, ustar(ikl))
+            ram_sv(ikl) = 1. / (cdM_SV(ikl) * max(VV__SV(ikl), eps6))
+            rah_sv(ikl) = 1. / (cdH_SV(ikl) * max(VV__SV(ikl), eps6))
+
+            ! +--Energy Fluxes                                          (INPUT/OUTPUT)
+            ! +  ^^^^^^^^^^^^^                                           ^^^^^^^^^^^^
+            !IF (.not.firstcall) THEN
+            Tsrfsv(ikl)  = tsurf(ikl)                     !hj 12 03 2010
+            cld_SV(ikl) = cloudf(ikl)                    ! Cloudiness
+            !END IF
+
+         END DO
+
+        !
+        ! +  +++  READ FORCINGS:  END  +++
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ 
+
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! +--SISVAT EXECUTION
+        ! +  ----------------
+
+        call  INLANDSIS(SnoMod, BloMod, 1)
+
+
+       
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        ! + RETURN RESULTS
+        ! + --------------
+        ! + Return (compressed) SISVAT variables to LMDZ
+        ! +
+        DO  ikl = 1, knon                  ! use only 1:knon (actual ice sheet..)
+            dflux_s(ikl) = dSdTSV(ikl)         ! Sens.H.Flux T-Der.
+            dflux_l(ikl) = dLdTSV(ikl)         ! Latn.H.Flux T-Der.
+            fluxsens(ikl) = HSs_sv(ikl)         ! HS
+            fluxlat(ikl) = HLs_sv(ikl)         ! HL
+            evap(ikl) = -1*HLs_sv(ikl) / LHvH2O  ! Evaporation
+            erod(ikl) = 0.
+
+            IF (BloMod) THEN
+                ! + Blowing snow
+
+                !       SLussl(i,j,n)= 0.
+                ! #BS   SLussl(i,j,n)=                     !Effective erosion
+                ! #BS. (- dbs_ER(ikl))/(dt*rhT_SV(ikl))    !~u*qs* from previous time step
+                ! #BS   blowSN(i,j,n)=  dt*uss_SV(ikl)     !New max. pot. Erosion [kg/m2]
+                ! #BS.                    *rhT_SV(ikl)     !(further bounded in sisvat_bsn.f)
+                ! #BS  erprev(i,j,n) =     dbs_Er(ikl)/dt__SV
+                erod(ikl) = dbs_Er(ikl) / dt__SV
+            ENDIF
+
+            ! +   Check snow thickness,  substract if too thick, add if too thin
+
+            sissnow(ikl) = 0.  !()
+            DO  isn = 1, isnoSV(ikl)
+                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
+            END DO
+
+            IF (sissnow(ikl) .LE. sn_low) THEN  !add snow
+                IF (isnoSV(ikl).GE.1) THEN
+                    dzsnSV(ikl, 1) = dzsnSV(ikl, 1) + sn_add / max(ro__SV(ikl, 1), epsi)
+                    toicSV(ikl) = toicSV(ikl) - sn_add
+                ELSE
+                    write(*, *) 'Attention, bare ice... point ', ikl
+                    isnoSV(ikl) = 1
+                    istoSV(ikl, 1) = 0
+                    ro__SV(ikl, 1) = 350.
+                    dzsnSV(ikl, 1) = sn_add / max(ro__SV(ikl, 1), epsi)  ! 1.
+                    eta_SV(ikl, 1) = epsi
+                    TsisSV(ikl, 1) = min(TsisSV(ikl, 0), TfSnow - 0.2)
+                    G1snSV(ikl, 1) = 0.
+                    G2snSV(ikl, 1) = 0.3
+                    agsnSV(ikl, 1) = 10.
+                    toicSV(ikl) = toicSV(ikl) - sn_add
+                END IF
+            END IF
+
+            IF (sissnow(ikl) .ge. sn_upp) THEN  !thinnen snow layer below
+                dzsnSV(ikl, 1) = dzsnSV(ikl, 1) / sn_div
+                toicSV(ikl) = toicSV(ikl) + dzsnSV(ikl, 1) * ro__SV(ikl, 1) / sn_div
+            END IF
+
+            sissnow(ikl) = 0.
+            qsnow(ikl) = 0.
+            snow(ikl) = 0.
+            snowhgt(ikl) = 0.
+
+            DO  isn = 1, isnoSV(ikl)
+                sissnow(ikl) = sissnow(ikl) + dzsnSV(ikl, isn) * ro__SV(ikl, isn)
+                snowhgt(ikl) = snowhgt(ikl) + dzsnSV(ikl, isn)
+                ! Etienne: check calc qsnow
+                qsnow(ikl) = qsnow(ikl) + rhoWat * eta_SV(ikl, isn) * dzsnSV(ikl, isn)
+            END DO
+
+            zfra(ikl) = max(min(isnoSV(ikl) - iiceSV(ikl), 1), 0)
+            ! Etienne: comment following line
+            ! snow(ikl)    = sissnow(ikl)+toicSV(ikl)
+            snow(ikl) = sissnow(ikl)
+
+            to_ice(ikl) = toicSV(ikl)
+            runoff_lic(ikl) = RnofSV(ikl)    ! RunOFF: intensity (flux due to melting + liquid precip)
+            fqfonte(ikl)= max(0., (wem_SV(ikl)-wer_SV(ikl))/dtime) ! net melting = melting - refreezing
+            ffonte(ikl)=fqfonte(ikl)*Lf_H2O
+
+            qsol(ikl) = 0.
+            DO  isl = -nsol, 0
+                tsoil(ikl, 1 - isl) = TsisSV(ikl, isl)       ! Soil Temperature
+                ! Etienne: check calc qsol
+                qsol(ikl) = qsol(ikl)                      &
+                        + eta_SV(ikl, isl) * dz_dSV(isl)
+            END DO
+            agesno(ikl) = agsnSV(ikl, isnoSV(ikl))        !          [day]
+
+            alb1(ikl) = alb1sv(ikl)             ! Albedo VIS
+!            alb2(ikl) = ((So1dSV - f1) * alb1sv(ikl)                   &
+!                    & + So2dSV * alb2sv(ikl) + So3dSV * alb3sv(ikl)) / f1
+            alb2(ikl)=alb2sv(ikl)
+            ! Albedo NIR
+            alb3(ikl) = alb3sv(ikl)             ! Albedo FIR
+            ! 6 band Albedo
+            alb6(ikl,:)=alb6sv(ikl,:)
+
+            tsurf_new(ikl) = Tsrfsv(ikl)
+
+            qsurf(ikl) = QsT_SV(ikl)
+            emis_new(ikl) = eps0SL(ikl)
+            z0m(ikl) = Z0m_SV(ikl)
+            z0h(ikl) = Z0h_SV(ikl)
+
+
+        END DO
+
+            IF (ok_outfor) THEN
+             ikl= gp_outfor
+            write(un_outfor, *) '+++++++++++', rlon(ikl2i(ikl)), rlat(ikl2i(ikl)),alt(ikl),'+++++++++++'
+            write(un_outfor, *) isnoSV(ikl), alb_SV(ikl), Z0m_SV(ikl), Z0h_SV(ikl),HSs_sv(ikl),HLs_sv(ikl),alb1(ikl),alb2(ikl)
+            write(un_outfor, *) dzsnSV(ikl, :)
+            write(un_outfor, *) TsisSV(ikl, :)
+            write(un_outfor, *) ro__SV(ikl, :)
+            write(un_outfor, *) eta_SV(ikl, :)
+            write(un_outfor, *) G1snSV(ikl, :)
+            write(un_outfor, *) G2snSV(ikl, :)
+            write(un_outfor, *) agsnSV(ikl, :)
+            write(un_outfor, *) istoSV(ikl, :)
+            write(un_outfor, *) DOPsnSV(ikl, :)
+        ENDIF
+
+
+
+        ! +  -----------------------------
+        ! +  END --- RETURN RESULTS
+        ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+        IF (lafin) THEN
+            fichnom = "restartsis.nc"
+            CALL sisvatredem("restartsis.nc", ikl2i, rlon, rlat)
+
+            IF (ok_outfor) THEN
+                close(unit = un_outfor)
+            END IF
+        END IF
+
+    END SUBROUTINE surf_inlandsis
+
+
+    !=======================================================================
+
+    SUBROUTINE get_soil_levels(dz1, dz2, lambda)
+        ! ======================================================================
+        ! Routine to compute the vertical discretization of the soil in analogy
+        ! to LMDZ. In LMDZ it is done in soil.F, which is not used in the case
+        ! of SISVAT, therefore it's needed here.
+        !
+        USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
+        USE mod_phys_lmdz_para
+        USE VAR_SV
+
+
+        !    INCLUDE "dimsoil.h"
+
+        REAL, DIMENSION(nsoilmx), INTENT(OUT) :: dz2, dz1
+        REAL, INTENT(OUT) :: lambda
+
+
+        !-----------------------------------------------------------------------
+        !   Depthts:
+        !   --------
+        REAL fz, rk, fz1, rk1, rk2
+        REAL min_period, dalph_soil
+        INTEGER ierr, jk
+
+        fz(rk) = fz1 * (dalph_soil**rk - 1.) / (dalph_soil - 1.)
+
+        !    write(*,*)'Start soil level computation'
+        !-----------------------------------------------------------------------
+        ! Calculation of some constants
+        ! NB! These constants do not depend on the sub-surfaces
+        !-----------------------------------------------------------------------
+        !-----------------------------------------------------------------------
+        !   ground levels
+        !   grnd=z/l where l is the skin depth of the diurnal cycle:
+        !-----------------------------------------------------------------------
+
+        min_period = 1800. ! en secondes
+        dalph_soil = 2.    ! rapport entre les epaisseurs de 2 couches succ.
+        ! !$OMP MASTER
+        !     IF (is_mpi_root) THEN
+        !        OPEN(99,file='soil.def',status='old',form='formatted',iostat=ierr)
+        !        IF (ierr == 0) THEN ! Read file only if it exists
+        !           READ(99,*) min_period
+        !           READ(99,*) dalph_soil
+        !           PRINT*,'Discretization for the soil model'
+        !           PRINT*,'First level e-folding depth',min_period, &
+        !                '   dalph',dalph_soil
+        !           CLOSE(99)
+        !        END IF
+        !     ENDIF
+        ! !$OMP END MASTER
+        !     CALL bcast(min_period)
+        !     CALL bcast(dalph_soil)
+
+        !   la premiere couche represente un dixieme de cycle diurne
+        fz1 = SQRT(min_period / 3.14)
+
+        DO jk = 1, nsoilmx
+            rk1 = jk
+            rk2 = jk - 1
+            dz2(jk) = fz(rk1) - fz(rk2)
+        ENDDO
+        DO jk = 1, nsoilmx - 1
+            rk1 = jk + .5
+            rk2 = jk - .5
+            dz1(jk) = 1. / (fz(rk1) - fz(rk2))
+        ENDDO
+        lambda = fz(.5) * dz1(1)
+        DO jk = 1, nsoilmx
+            rk = jk
+            rk1 = jk + .5
+            rk2 = jk - .5
+        ENDDO
+
+    END SUBROUTINE get_soil_levels
+
+
+    !===========================================================================
+
+    SUBROUTINE SISVAT_ini(knon)
+
+        !C +------------------------------------------------------------------------+
+        !C | MAR          SISVAT_ini                             Jd 11-10-2007  MAR |
+        !C |   SubRoutine SISVAT_ini generates non time dependant SISVAT parameters |
+        !C +------------------------------------------------------------------------+
+        !C |   PARAMETERS:  klonv: Total Number of columns =                        |
+        !C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
+        !C |                     X       Number of Mosaic Cell per grid box         |
+        !C |                                                                        |
+        !C |   INPUT:   dt__SV   : Time  Step                                   [s] |
+        !C |   ^^^^^    dz_dSV   : Layer Thickness                              [m] |
+        !C |                                                                        |
+        !C |   OUTPUT:             [-] |
+        !C |   ^^^^^^   rocsSV   : Soil Contrib. to (ro c)_s exclud.Water  [J/kg/K] |
+        !C |            etamSV   : Soil Minimum Humidity                    [m3/m3] |
+        !C |                      (based on a prescribed Soil Relative Humidity)    |
+        !C |            s1__SV   : Factor of eta**( b+2) in Hydraul.Diffusiv.       |
+        !C |            s2__SV   : Factor of eta**( b+2) in Hydraul.Conduct.        |
+        !C |            aKdtSV   : KHyd: Piecewise Linear Profile:  a * dt    [m]   |
+        !C |            bKdtSV   : KHyd: Piecewise Linear Profile:  b * dt    [m/s] |
+        !C |            dzsnSV(0): Soil first Layer Thickness                   [m] |
+        !C |            dzmiSV   : Distance between two contiguous levels       [m] |
+        !C |            dz78SV   : 7/8 (Layer Thickness)                        [m] |
+        !C |            dz34SV   : 3/4 (Layer Thickness)                        [m] |
+        !C |            dz_8SV   : 1/8 (Layer Thickness)                        [m] |
+        !C |            dzAvSV   : 1/8  dz_(i-1) + 3/4 dz_(i) + 1/8 dz_(i+1)    [m] |
+        !C |            dtz_SV   : dt/dz                                      [s/m] |
+        !C |            OcndSV   : Swab Ocean / Soil Ratio                      [-] |
+        !C |            Implic   : Implicit Parameter  (0.5:  Crank-Nicholson)      |
+        !C |            Explic   : Explicit Parameter = 1.0 - Implic                |
+        !C |                                                                        |
+        !C | # OPTIONS: #ER: Richards Equation is not smoothed                      |
+        !C | # ^^^^^^^  #kd: De Ridder   Discretization                             |
+        !C | #          #SH: Hapex-Sahel Values                                     !
+        !C |                                                                        |
+        !C +------------------------------------------------------------------------+
+        !
+        !
+
+        !C +--Global Variables
+        !C +  ================
+
+        USE dimphy
+        USE VARphy
+        USE VAR_SV
+        USE VARdSV
+        USE VAR0SV
+        USE VARxSV
+        USE VARtSV
+        USE VARxSV
+        USE VARySV
+        IMPLICIT NONE
+
+
+
+        !C +--Arguments
+        !C +  ==================
+        INTEGER, INTENT(IN) :: knon
+
+        !C +--Internal Variables
+        !C +  ==================
+
+        INTEGER :: ivt, ist, ikl, isl, isn, ikh
+        INTEGER :: misl_2, nisl_2
+        REAL :: d__eta, eta__1, eta__2, Khyd_1, Khyd_2
+        REAL, PARAMETER :: RHsMin = 0.001        ! Min.Soil Relative Humidity
+        REAL :: PsiMax                        ! Max.Soil Water    Potential
+        REAL :: a_Khyd, b_Khyd                 ! Water conductivity
+
+
+        !c #WR REAL    ::  Khyd_x,Khyd_y
+
+
+
+        !C +--Non Time Dependant SISVAT parameters
+        !C +  ====================================
+
+        !C +--Soil Discretization
+        !C +  -------------------
+
+        !C +--Numerical Scheme Parameters
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        Implic = 0.75                           ! 0.5  <==> Crank-Nicholson  
+        Explic = 1.00 - Implic                  !                            
+
+        !C +--Soil/Snow Layers Indices
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^
+        DO  isl = -nsol, 0
+            islpSV(isl) = isl + 1
+            islpSV(isl) = min(islpSV(isl), 0)
+            islmSV(isl) = isl - 1
+            islmSV(isl) = max(-nsol, islmSV(isl))
+        END DO
+
+        DO  isn = 1, nsno
+            isnpSV(isn) = isn + 1
+            isnpSV(isn) = min(isnpSV(isn), nsno)
+        END DO
+
+        !C +--Soil      Layers Thicknesses
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        ! Not used here as LMDZ method is applied, see SUBROUTINE get_soil_levels!
+        !c #kd IF (nsol.gt.4)                                              THEN
+        !c #kd   DO isl=-5,-nsol,-1
+        !c #kd     dz_dSV(isl)=   1.
+        !c #kd   END DO
+        !c #kd END IF
+        !
+        !      IF (nsol.ne.4)                                              THEN
+        !        DO isl= 0,-nsol,-1
+        !          misl_2 =     -mod(isl,2)
+        !          nisl_2 =         -isl/2
+        !          dz_dSV(isl)=(((1-misl_2) * 0.001
+        !     .                  +  misl_2  * 0.003) * 10**(nisl_2)) * 4.
+        !C +...    dz_dSV(0)  =         Hapex-Sahel Calibration:       4 mm
+        !
+        !c +SH     dz_dSV(isl)=(((1-misl_2) * 0.001
+        !c +SH.                  +  misl_2  * 0.003) * 10**(nisl_2)) * 1.
+        !
+        !c #05     dz_dSV(isl)=(((1-misl_2) * 0.001
+        !c #05.                  +  misl_2  * 0.008) * 10**(nisl_2)) * 0.5
+        !        END DO
+        !          dz_dSV(0)  =               0.001
+        !          dz_dSV(-1) = dz_dSV(-1)  - dz_dSV(0)              + 0.004
+        !      END IF
+
+        zz_dSV = 0.
+        DO  isl = -nsol, 0
+            dzmiSV(isl) = 0.500 * (dz_dSV(isl) + dz_dSV(islmSV(isl)))
+            dziiSV(isl) = 0.500 * dz_dSV(isl) / dzmiSV(isl)
+            dzi_SV(isl) = 0.500 * dz_dSV(islmSV(isl)) / dzmiSV(isl)
+            dtz_SV(isl) = dt__SV / dz_dSV(isl)
+            dtz_SV2(isl) = 1. / dz_dSV(isl)
+            dz78SV(isl) = 0.875 * dz_dSV(isl)
+            dz34SV(isl) = 0.750 * dz_dSV(isl)
+            dz_8SV(isl) = 0.125 * dz_dSV(isl)
+            dzAvSV(isl) = 0.125 * dz_dSV(islmSV(isl))                        &
+                    & + 0.750 * dz_dSV(isl)                                &
+                    & + 0.125 * dz_dSV(islpSV(isl))
+            zz_dSV = zz_dSV + dz_dSV(isl)
+        END DO
+        DO ikl = 1, knon !v
+            dzsnSV(ikl, 0) = dz_dSV(0)
+        END DO
+
+        !C +--Conversion to a 50 m Swab Ocean Discretization
+        !C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        OcndSV = 0.
+        DO isl = -nsol, 0
+            OcndSV = OcndSV + dz_dSV(isl)
+        END DO
+        OcndSV = 50. / OcndSV
+
+
+        !C +--Secondary Soil       Parameters
+        !C +  -------------------------------
+
+        DO  ist = 0, nsot
+            rocsSV(ist) = (1.0 - etadSV(ist)) * 1.2E+6   ! Soil Contrib. to (ro c)_s
+            s1__SV(ist) = bCHdSV(ist)          & ! Factor of (eta)**(b+2)
+                    & * psidSV(ist) * Ks_dSV(ist)          & !    in DR97, Eqn.(3.36)
+                    & / (etadSV(ist)**(bCHdSV(ist) + 3.))     !
+            s2__SV(ist) = Ks_dSV(ist)          & ! Factor of (eta)**(2b+3)
+                    & / (etadSV(ist)**(2. * bCHdSV(ist) + 3.))     !    in DR97, Eqn.(3.35)
+
+            !C +--Soil Minimum Humidity (from a prescribed minimum relative Humidity)
+            !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+            Psimax = -(log(RHsMin)) / 7.2E-5        ! DR97, Eqn 3.15 Inversion
+            etamSV(ist) = etadSV(ist)                                      &
+                    & * (PsiMax / psidSV(ist))**(-min(10., 1. / bCHdSV(ist)))
+        END DO
+        etamSV(12) = 0.
+
+        !C +--Piecewise Hydraulic Conductivity Profiles
+        !C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+        DO   ist = 0, nsot
+
+            d__eta = etadSV(ist) / nkhy
+            eta__1 = 0.
+            eta__2 = d__eta
+            DO ikh = 0, nkhy
+                Khyd_1 = s2__SV(ist)             & ! DR97, Eqn.(3.35)
+                        & * (eta__1      **(2. * bCHdSV(ist) + 3.))        !
+                Khyd_2 = s2__SV(ist)             &!
+                        & * (eta__2      **(2. * bCHdSV(ist) + 3.))        !
+
+                a_Khyd = (Khyd_2 - Khyd_1) / d__eta   !
+                b_Khyd = Khyd_1 - a_Khyd * eta__1   !
+                !c #WR     Khyd_x          =  a_Khyd*eta__1 +b_Khyd   !
+                !c #WR     Khyd_y          =  a_Khyd*eta__2 +b_Khyd   !
+                aKdtSV(ist, ikh) = a_Khyd * dt__SV   !
+                bKdtSV(ist, ikh) = b_Khyd * dt__SV   !
+
+                eta__1 = eta__1 + d__eta
+                eta__2 = eta__2 + d__eta
+            END DO
+        END DO
+
+        return
+
+    END SUBROUTINE SISVAT_ini
+
+
+    !***************************************************************************
+
+    SUBROUTINE sisvatetat0 (fichnom, ikl2i)
+
+        USE dimphy
+        USE mod_grid_phy_lmdz
+        USE mod_phys_lmdz_para
+
+        USE iostart
+        USE VAR_SV
+        USE VARdSV
+        USE VARxSV
+        USE VARtSV
+        USE indice_sol_mod
+
+        IMPLICIT none
+        !======================================================================
+        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
+        ! Objet: Lecture du fichier de conditions initiales pour SISVAT
+        !======================================================================
+        include "netcdf.inc"
+        !    include "indicesol.h"
+
+        !    include "dimsoil.h"
+        include "clesphys.h"
+        include "compbl.h"
+
+        !======================================================================
+        CHARACTER(LEN = *) :: fichnom
+
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
+        REAL, DIMENSION(klon) :: rlon
+        REAL, DIMENSION(klon) :: rlat
+
+        ! les variables globales ecrites dans le fichier restart
+        REAL, DIMENSION(klon) :: isno
+        REAL, DIMENSION(klon) :: ispi
+        REAL, DIMENSION(klon) :: iice
+        REAL, DIMENSION(klon) :: rusn
+        REAL, DIMENSION(klon, nsno) :: isto
+
+        REAL, DIMENSION(klon, nsismx) :: Tsis
+        REAL, DIMENSION(klon, nsismx) :: eta
+        REAL, DIMENSION(klon, nsismx) :: ro
+
+        REAL, DIMENSION(klon, nsno) :: dzsn
+        REAL, DIMENSION(klon, nsno) :: G1sn
+        REAL, DIMENSION(klon, nsno) :: G2sn
+        REAL, DIMENSION(klon, nsno) :: agsn
+
+        REAL, DIMENSION(klon) :: toic
+
+        INTEGER :: isl, ikl, i, isn, errT, erreta, errro, errdz, snopts
+        CHARACTER (len = 2) :: str2
+        LOGICAL :: found
+
+        errT = 0
+        errro = 0
+        erreta = 0
+        errdz = 0
+        snopts = 0
+        ! Ouvrir le fichier contenant l'etat initial:
+
+        CALL open_startphy(fichnom)
+
+        ! Lecture des latitudes, longitudes (coordonnees):
+
+        CALL get_field("latitude", rlat, found)
+        CALL get_field("longitude", rlon, found)
+
+        CALL get_field("n_snows", isno, found)
+        IF (.NOT. found) THEN
+            PRINT*, 'phyetat0: Le champ <n_snows> est absent'
+            PRINT *, 'fichier startsisvat non compatible avec sisvatetat0'
+        ENDIF
+
+        CALL get_field("n_ice_top", ispi, found)
+        CALL get_field("n_ice", iice, found)
+        CALL get_field("surf_water", rusn, found)
+
+
+        CALL get_field("to_ice", toic, found)
+        IF (.NOT. found) THEN
+            PRINT*, 'phyetat0: Le champ <to_ice> est absent'
+            toic(:) = 0.
+        ENDIF
+
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("AGESNOW" // str2, &
+                        agsn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("DZSNOW" // str2, &
+                        dzsn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("G2SNOW" // str2, &
+                        G2sn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("G1SNOW" // str2, &
+                        G1sn(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("ETA" // str2, &
+                        eta(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("RO" // str2, &
+                        ro(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsismx
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("TSS" // str2, &
+                        Tsis(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        DO isn = 1, nsno
+            IF (isn.LE.99) THEN
+                WRITE(str2, '(i2.2)') isn
+                CALL get_field("HISTORY" // str2, &
+                        isto(:, isn), found)
+            ELSE
+                PRINT*, "Trop de couches"
+                CALL abort
+            ENDIF
+        ENDDO
+        write(*, *)'Read ', fichnom, ' finished!!'
+
+        !*********************************************************************************
+        ! Compress restart file variables for SISVAT
+
+        DO  ikl = 1, klon
+            i = ikl2i(ikl)
+            IF (i > 0) THEN
+                isnoSV(ikl) = INT(isno(i))          ! Nb Snow/Ice Lay.
+                ispiSV(ikl) = INT(ispi(i))          ! Nb Supr.Ice Lay.
+                iiceSV(ikl) = INT(iice(i))          ! Nb      Ice Lay.
+
+                DO isl = -nsol, 0
+                    ro__SV(ikl, isl) = ro(i, nsno + 1 - isl)       !
+                    eta_SV(ikl, isl) = eta(i, nsno + 1 - isl)         ! Soil Humidity
+                    !hjp 15/10/2010
+                    IF (eta_SV(ikl, isl) <= 1.e-6) THEN          !hj check
+                        eta_SV(ikl, isl) = 1.e-6
+                    ENDIF
+                    TsisSV(ikl, isl) = Tsis(i, nsno + 1 - isl)        ! Soil Temperature
+                    IF (TsisSV(ikl, isl) <= 1.) THEN             !hj check
+                        !                errT=errT+1
+                        TsisSV(ikl, isl) = 273.15 - 0.2              ! Etienne: negative temperature since soil is ice
+                    ENDIF
+
+                END DO
+                write(*, *)'Copy histo', ikl
+
+                DO  isn = 1, isnoSV(ikl) !nsno
+                    snopts = snopts + 1
+                    IF (isto(i, isn) > 10.) THEN          !hj check
+                        write(*, *)'Irregular isto', ikl, i, isn, isto(i, isn)
+                        isto(i, isn) = 1.
+                    ENDIF
+
+                    istoSV(ikl, isn) = INT(isto(i, isn))     ! Snow     History
+                    ro__SV(ikl, isn) = ro(i, isn)            !        [kg/m3]
+                    eta_SV(ikl, isn) = eta(i, isn)           !        [m3/m3]
+                    TsisSV(ikl, isn) = Tsis(i, isn)          !            [K]
+
+                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
+                        errT = errT + 1
+                        TsisSV(ikl, isn) = TsisSV(ikl, 0)
+                    ENDIF
+                    IF (TsisSV(ikl, isn) <= 1.) THEN          !hj check
+                        TsisSV(ikl, isn) = 263.15
+                    ENDIF
+                    IF (eta_SV(ikl, isn) < 1.e-9) THEN          !hj check
+                        eta_SV(ikl, isn) = 1.e-6
+                        erreta = erreta + 1
+                    ENDIF
+                    IF (ro__SV(ikl, isn) <= 10.) THEN          !hj check
+                        ro__SV(ikl, isn) = 11.
+                        errro = errro + 1
+                    ENDIF
+                    write(*, *)ikl, i, isn, Tsis(i, isn), G1sn(i, isn)
+                    G1snSV(ikl, isn) = G1sn(i, isn)          ! [-]        [-]
+                    G2snSV(ikl, isn) = G2sn(i, isn)          ! [-] [0.0001 m]
+                    dzsnSV(ikl, isn) = dzsn(i, isn)          !            [m]
+                    agsnSV(ikl, isn) = agsn(i, isn)          !          [day]
+                END DO
+                rusnSV(ikl) = rusn(i)              ! Surficial Water
+                toicSV(ikl) = toic(i)              ! bilan snow to ice
+            END IF
+        END DO
+
+    END SUBROUTINE sisvatetat0
+
+
+    !======================================================================
+    SUBROUTINE sisvatredem (fichnom, ikl2i, rlon, rlat)
+
+
+
+        !======================================================================
+        ! Auteur(s) HJ PUNGE (LSCE) date: 07/2009
+        ! Objet: Ecriture de l'etat de redemarrage pour SISVAT
+        !======================================================================
+        USE mod_grid_phy_lmdz
+        USE mod_phys_lmdz_para
+        USE iostart
+        USE VAR_SV
+        USE VARxSV
+        USE VARySV !hj tmp 12 03 2010
+        USE VARtSV
+        USE indice_sol_mod
+        USE dimphy
+
+        IMPLICIT none
+
+        include "netcdf.inc"
+        include "clesphys.h"
+        include "compbl.h"
+
+        !======================================================================
+
+        CHARACTER(LEN = *) :: fichnom
+        INTEGER, DIMENSION(klon), INTENT(IN) :: ikl2i
+        REAL, DIMENSION(klon), INTENT(IN) :: rlon
+        REAL, DIMENSION(klon), INTENT(IN) :: rlat
+
+        ! les variables globales ecrites dans le fichier restart
+        REAL, DIMENSION(klon) :: isno
+        REAL, DIMENSION(klon) :: ispi
+        REAL, DIMENSION(klon) :: iice
+        REAL, DIMENSION(klon, nsnowmx) :: isto
+
+        REAL, DIMENSION(klon, nsismx) :: Tsis
+        REAL, DIMENSION(klon, nsismx) :: eta
+        REAL, DIMENSION(klon, nsnowmx) :: dzsn
+        REAL, DIMENSION(klon, nsismx) :: ro
+        REAL, DIMENSION(klon, nsnowmx) :: G1sn
+        REAL, DIMENSION(klon, nsnowmx) :: G2sn
+        REAL, DIMENSION(klon, nsnowmx) :: agsn
+        REAL, DIMENSION(klon) :: IRs
+        REAL, DIMENSION(klon) :: LMO
+        REAL, DIMENSION(klon) :: rusn
+        REAL, DIMENSION(klon) :: toic
+        REAL, DIMENSION(klon) :: Bufs
+        REAL, DIMENSION(klon) :: alb1, alb2, alb3
+
+        INTEGER isl, ikl, i, isn, ierr
+        CHARACTER (len = 2) :: str2
+        INTEGER :: pass
+
+        isno(:) = 0
+        ispi(:) = 0
+        iice(:) = 0
+        IRs(:) = 0.
+        LMO(:) = 0.
+        eta(:, :) = 0.
+        Tsis(:, :) = 0.
+        isto(:, :) = 0
+        ro(:, :) = 0.
+        G1sn(:, :) = 0.
+        G2sn(:, :) = 0.
+        dzsn(:, :) = 0.
+        agsn(:, :) = 0.
+        rusn(:) = 0.
+        toic(:) = 0.
+        Bufs(:) = 0.
+        alb1(:) = 0.
+        alb2(:) = 0.
+        alb3(:) = 0.
+
+        !***************************************************************************
+        ! Uncompress SISVAT output variables for storage
+
+        DO  ikl = 1, klon
+            i = ikl2i(ikl)
+            IF (i > 0) THEN
+                isno(i) = 1. * isnoSV(ikl)               ! Nb Snow/Ice Lay.
+                ispi(i) = 1. * ispiSV(ikl)               ! Nb Supr.Ice Lay.
+                iice(i) = 1. * iiceSV(ikl)               ! Nb      Ice Lay.
+
+                !        IRs(i)        = IRs_SV(ikl)
+                !        LMO(i)        = LMO_SV(ikl)
+
+                DO isl = -nsol, 0                           !
+                    eta(i, nsno + 1 - isl) = eta_SV(ikl, isl)            ! Soil Humidity
+                    Tsis(i, nsno + 1 - isl) = TsisSV(ikl, isl)            ! Soil Temperature
+                    ro(i, nsno + 1 - isl) = ro__SV(ikl, isl)            !        [kg/m3]
+                END DO
+
+                DO  isn = 1, nsno
+                    isto(i, isn) = 1. * istoSV(ikl, isn)         ! Snow     History
+                    ro(i, isn) = ro__SV(ikl, isn)            !        [kg/m3]
+                    eta(i, isn) = eta_SV(ikl, isn)            !        [m3/m3]
+                    Tsis(i, isn) = TsisSV(ikl, isn)            !            [K]
+                    G1sn(i, isn) = G1snSV(ikl, isn)            ! [-]        [-]
+                    G2sn(i, isn) = G2snSV(ikl, isn)            ! [-] [0.0001 m]
+                    dzsn(i, isn) = dzsnSV(ikl, isn)            !            [m]
+                    agsn(i, isn) = agsnSV(ikl, isn)            !          [day]
+                END DO
+                rusn(i) = rusnSV(ikl)                  ! Surficial Water
+                toic(i) = toicSV(ikl)                  ! to ice
+                alb1(i) = alb1sv(ikl)
+                alb2(i) = alb2sv(ikl)
+                alb3(i) = alb3sv(ikl)
+                !        Bufs(i)       = BufsSV(ikl)
+            END IF
+        END DO
+
+        CALL open_restartphy(fichnom)
+
+        DO pass = 1, 2
+            CALL put_field(pass, "longitude", &
+                    "Longitudes de la grille physique", rlon)
+            CALL put_field(pass, "latitude", "Latitudes de la grille physique", rlat)
+
+            CALL put_field(pass, "n_snows", "number of snow/ice layers", isno)
+            CALL put_field(pass, "n_ice_top", "number of top ice layers", ispi)
+            CALL put_field(pass, "n_ice", "number of ice layers", iice)
+            CALL put_field(pass, "IR_soil", "Soil IR flux", IRs)
+            CALL put_field(pass, "LMO", "Monin-Obukhov Scale", LMO)
+            CALL put_field(pass, "surf_water", "Surficial water", rusn)
+            CALL put_field(pass, "snow_buffer", "Snow buffer layer", Bufs)
+            CALL put_field(pass, "alb_1", "albedo sw", alb1)
+            CALL put_field(pass, "alb_2", "albedo nIR", alb2)
+            CALL put_field(pass, "alb_3", "albedo fIR", alb3)
+            CALL put_field(pass, "to_ice", "Snow passed to ice", toic)
+
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "AGESNOW" // str2, &
+                            "Age de la neige layer No." // str2, &
+                            agsn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "DZSNOW" // str2, &
+                            "Snow/ice thickness layer No." // str2, &
+                            dzsn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "G2SNOW" // str2, &
+                            "Snow Property 2, layer No." // str2, &
+                            G2sn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "G1SNOW" // str2, &
+                            "Snow Property 1, layer No." // str2, &
+                            G1sn(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "ETA" // str2, &
+                            "Soil/snow water content layer No." // str2, &
+                            eta(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx   !nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "RO" // str2, &
+                            "Snow density layer No." // str2, &
+                            ro(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsismx
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "TSS" // str2, &
+                            "Soil/snow temperature layer No." // str2, &
+                            Tsis(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+            DO isn = 1, nsno
+                IF (isn.LE.99) THEN
+                    WRITE(str2, '(i2.2)') isn
+                    CALL put_field(pass, "HISTORY" // str2, &
+                            "Snow history layer No." // str2, &
+                            isto(:, isn))
+                ELSE
+                    PRINT*, "Trop de couches"
+                    CALL abort
+                ENDIF
+            ENDDO
+
+            CALL enddef_restartphy
+        ENDDO
+        CALL close_restartphy
+
+    END SUBROUTINE sisvatredem
+
+END MODULE surf_inlandsis_mod
