diff --git a/setup/testlake_driver.dat b/setup/testlake_driver.dat
index 15a33d6ab3d7c1227ea1f64fb60d76ac3a6bf499..4c8ec4ce9d40414608a40e89069d4cf4102d0ed8 100644
--- a/setup/testlake_driver.dat
+++ b/setup/testlake_driver.dat
@@ -71,6 +71,8 @@ N_LatentFlux    -1
 N_Ustar         -1
 N_surfrad       -1
 N_cloud         -1
+N_NetRad        -1
+N_SurfTemp      -1
 #
 #-----------------------------------------------------------------------------------------
 #                              TIME INTEGRATION PARAMETERS
@@ -95,6 +97,8 @@ hour0       13.33333333
 tinteg      10.8333333
 spinup_times  0
 spinup_period 0
+cp_period     0
+control_point 0
 dt          5
 call_Flake  0
 #
@@ -219,7 +223,7 @@ morphometry    5
 1	4.E+6
 2       3.E+6
 4       1.E+6
-7       100.
+7.1       100.
 #
 #-----------------------------------------------------------------------------------------
 #                               NETCDF OUTPUT PARAMETERS
diff --git a/setup/testlake_setup.dat b/setup/testlake_setup.dat
index 1019ff06a665bc1ce78c7dee64446ee0c7438d19..d96166daaf21952fbb039d519c93a54fe2a21d84 100644
--- a/setup/testlake_setup.dat
+++ b/setup/testlake_setup.dat
@@ -146,6 +146,11 @@ ifbubble 1
 sedim    0
 salsoil  0
 dyn_pgrad 0
+pgrad    0.
+nManning 5.E-2
+horvisc  0.
+backdiff 0
+botfric  1
 zero_model 0
 thermokarst_meth_prod 0.
 soil_meth_prod 1.
@@ -179,12 +184,16 @@ deadvol 0.
 #4.0946210218    11.6325191677   23.7094310675E-3  0.  0.  0.
 #7.              5.              30.E-3            0.  0.  0.
 T_profile 6
-0.2             15.3            1.0967626609E-3  0.  0.  0.
-1.15            15.3            1.2157976331E-3  0.  0.  0.
-2.1             15.3            6.5826901114E-3  0.  0.  0.
-3.05            11.8            20.3271090163E-3  0.  0.  0.
-4.0             10.1             23.7094310675E-3  0.  0.  0.
-7.              5.              30.E-3            0.  0.  0.
+0.2             15.3            1.0967626609E-3  0.  0.  0. 0.
+1.15            15.3            1.2157976331E-3  0.  0.  0. 0.
+2.1             15.3            6.5826901114E-3  0.  0.  0. 0.
+3.05            11.8            20.3271090163E-3  0.  0.  0. 0.
+4.0             10.1            23.7094310675E-3  0.  0.  0. 0.
+7.              5.              30.E-3            0.  0.  0. 0.
+#
+T_soilprofile 2
+0.  4.
+10. 4.
 #.
 #----------------------------------------------------------------------------------------
 #                               BOUNDARY CONDITIONS: TRIBUTARIES AND EFFLUENTS
@@ -206,51 +215,6 @@ iefflloc 1
 fileinflow 'BolshoiVilui20153_inflows.dat'
 fileoutflow 'BolshoiVilui20153_outflow.dat'
 dttribupdate 0.25
-#inflowprof
-#1  100. 0.05 8. 0. 0.05 0.  
-#2  100. 0.05 8. 0. 0.05 0.
-#3  100. 0.05 8. 0. 0.05 0.
-#3  100. 0.0  8. 0. 0.05 0.
-#5  100. 0.0  8. 0. 0.05 0.
-#6  100. 0.0  8. 0. 0.05 0.
-#7  100. 0.0  8. 0. 0.05 0.
-#8  100. 0.0  8. 0. 0.05 0.
-#9  100. 0.0  8. 0. 0.05 0.
-#10 100. 0.0  8. 0. 0.05 0.
-#11 100. 0.0  8. 0. 0.05 0.
-#12 100. 0.0  8. 0. 0.05 0.
-#13 100. 0.0  8. 0. 0.05 0.
-#14 100. 0.0  8. 0. 0.05 0.
-#15 100. 0.0  8. 0. 0.05 0.
-#16 100. 0.0  8. 0. 0.05 0.
-#17 100. 0.0  8. 0. 0.05 0.
-#18 100. 0.0  8. 0. 0.05 0.
-#19 100. 0.0  8. 0. 0.05 0.
-#20 100. 0.0  8. 0. 0.05 0.
-#21 100. 0.0  8. 0. 0.05 0.
-#
-#outflowprof
-#1  100. 0.05 
-#2  100. 0.05
-#3  100. 0.05
-#4  100. 0.
-#5  100. 0.
-#6  100. 0.
-#7  100. 0.
-#8  100. 0.
-#9  100. 0.
-#10 100. 0.
-#11 100. 0.
-#12 100. 0.
-#13 100. 0.
-#14 100. 0.
-#15 100. 0.
-#16 100. 0.
-#17 100. 0.
-#18 100. 0.
-#19 100. 0.
-#20 100. 0.
-#21 100. 0.
 #
 #----------------------------------------------------------------------------------------
 #            DATA ASSIMILATION CONTROLS (NOT OPERATIONAL: PUT EVERYTHING TO 0)
@@ -302,6 +266,7 @@ nscreen     1000
 scale_output 0
 accum_begin 2003060100
 accum_end   2004060100
+zserout -999.
 rtemp 1
 -999. -999. -999.
 #
@@ -315,6 +280,10 @@ ngrid_out 8
 6.
 7.
 #
+ngridice_out 2
+0.
+0.2
+#
 ngridsoil_out 11
 0.
 1.
diff --git a/source/Makefile b/source/Makefile
index 875071eb31ca0bf9330b10c8ededba991833b6fe..d0a1ff6a73d776b5538d341636f08c7a1e59df16 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -18,7 +18,7 @@
  ifeq ($(FC),ifort)
    PREPROCESS_KEY = -fpp
    opt_keys = -qopenmp -O3 #-fp-model source
-   check_keys = -g -traceback -check all -fpe-all=0
+   check_keys = -traceback #-g -check all -fpe-all=0
  endif
  ifeq ($(FC),gfortran)
    PREPROCESS_KEY = -cpp
diff --git a/source/model/bathym_mod.f90 b/source/model/bathym_mod.f90
index 8079f2a80065e9491e77a6276ee2f7afb9173e65..5ca2124dd07d8dad8a3effa25498cd1c0e060a95 100644
--- a/source/model/bathym_mod.f90
+++ b/source/model/bathym_mod.f90
@@ -114,9 +114,9 @@ else
     if (init(ix,iy) == 0) then
       allocate(work5(1:nsoilcols+1), work6(1:nsoilcols))
       call LININTERPOL (depth_area(1,1),depth_area(1,2),ndatamax, &
-                       & zsoilcols(1,ix,iy),work5,nsoilcols+1,flag)
+                       & zsoilcols(1,ix,iy),work5,nsoilcols+1,flag,.true.)
       call LININTERPOL (depth_area(1,1),depth_area(1,2),ndatamax, &
-                       & work3,work6,nsoilcols,flag) 
+                       & work3,work6,nsoilcols,flag,.true.) 
       bathymsoil(1:nsoilcols+1,ix,iy)%area_int  = work5(1:nsoilcols+1)
       bathymsoil(1:nsoilcols,  ix,iy)%area_half = work6(1:nsoilcols)
       deallocate (work5, work6)
@@ -199,8 +199,8 @@ else
       work2(:) = depth_area(:,2)
       work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + &
       & h1 + ls1 ! Relating coordinate systems before interpolation
-      call LININTERPOL (work1,work2,ndatamax,z_full,area_int,M+1,flag)
-      call LININTERPOL (work1,work2,ndatamax,z_half,area_half,M,flag)
+      call LININTERPOL (work1,work2,ndatamax,z_full,area_int,M+1,flag,.true.)
+      call LININTERPOL (work1,work2,ndatamax,z_half,area_half,M,flag,.true.)
     endif
   endif
   !... in ice layer
@@ -222,8 +222,8 @@ else
     elseif (soilcolconjtype == 2) then
       work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + &
       & h1 + l1 + ls1 
-      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag)
-      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag) 
+      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag,.true.)
+      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag,.true.) 
     endif
     bathymice(1:Mice+1)%area_int  = work5(1:Mice+1)
     bathymice(1:Mice)  %area_half = work6(1:Mice)
@@ -249,8 +249,8 @@ else
       work5(Mice+1) = bathymsoil(nsoilcols+1,ix,iy)%area_int
     elseif (soilcolconjtype == 2) then
       work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + ls1 
-      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag)
-      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag) 
+      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag,.true.)
+      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag,.true.) 
     endif
     bathymdice(1:Mice+1)%area_int  = work5(1:Mice+1)
     bathymdice(1:Mice)  %area_half = work6(1:Mice)
diff --git a/source/model/init_var_mod.f90 b/source/model/init_var_mod.f90
index 155236187fdd2805f9c546a9408e8661ef2cbfe6..b29d1b072f3b843041b52400a5afa49e8056a296 100644
--- a/source/model/init_var_mod.f90
+++ b/source/model/init_var_mod.f90
@@ -176,9 +176,9 @@
     ! The initial temperature profile is given from the input file
   elseif (init_T == 3) then
     call LININTERPOL (ip%zTinitprof,ip%Tinitprof,ip%lenprof, &
-    & z_full,Tw1,M+1,flag)
+    & z_full,Tw1,M+1,flag,.true.)
     call LININTERPOL (ip%zTinitprof,ip%Sinitprof,ip%lenprof, &
-    & z_full,Sal1,M+1,flag)
+    & z_full,Sal1,M+1,flag,.true.)
     if (.not.flag) then
       print*, 'The error while interpolating the initial &
       &temperature profile: terminating program'
@@ -246,7 +246,7 @@
     ! The initial soil temperature profile is given from the input file
     elseif (init_T == 3) then
       call LININTERPOL (isp%zTinitprof,isp%Tinitprof,isp%lenprof, &
-      & zsoil,Tsoil1(1,nsoilcols),ns,flag)
+      & zsoil,Tsoil1(1,nsoilcols),ns,flag,.true.)
       if (.not.flag) then
         print*, 'The error while interpolating the initial &
         &soil temperature profile: terminating program'
@@ -309,7 +309,7 @@
     enddo
   else
     call LININTERPOL (ip%zTinitprof,ip%ch4initprof,ip%lenprof, &
-    & z_full,qwater,M+1,flag)
+    & z_full,qwater,M+1,flag,.true.)
   endif
   
 
@@ -348,7 +348,7 @@
     DIC(1:M+1) = 10.*co2_atm0 ! Atmospheric concentration
   else
     call LININTERPOL (ip%zTinitprof,ip%co2initprof,ip%lenprof, &
-    & z_full,DIC,M+1,flag)
+    & z_full,DIC,M+1,flag,.true.)
   endif
 
 ! Inorganic dissolved phosphorus
@@ -356,7 +356,7 @@
     phosph(:) = 0.01/molmass_p
   else
     call LININTERPOL (ip%zTinitprof,ip%phosphinitprof,ip%lenprof, &
-    & z_full,phosph,M+1,flag)
+    & z_full,phosph,M+1,flag,.true.)
   endif
 
 ! Oxygen initialization
@@ -371,7 +371,7 @@
     enddo
   else
     call LININTERPOL (ip%zTinitprof,ip%o2initprof,ip%lenprof, &
-    & z_full,oxyg,M+1,flag)
+    & z_full,oxyg,M+1,flag,.true.)
   endif
   oxygsoil(1:nsoilcols) = 0. !Zero initial oxygen concentration in the aerobic soil layer