diff --git a/crproj b/crproj
index 6214e1690cfc5148545d24abfeef41558ee22817..84662e5ac7b12351d79282db612cf3a74db7a109 100755
--- a/crproj
+++ b/crproj
@@ -1,5 +1,17 @@
 #!/bin/sh
 
+if [[ "$OSTYPE" == "linux-gnu"* || "$OSTYPE" == "cygwin" ]]; then
+        echo "Operating system: linux"
+	OS="linux"
+elif [[ "$OSTYPE" == "darwin"* ]]; then
+        echo "Operating system: OS X"
+	OS="OSX"
+else
+	echo "Uknown operating system"
+	# Unknown OS.
+fi
+
+
 mkdir -p -v results/$1/everystep
 mkdir -p -v results/$1/netcdf
 mkdir -p -v results/$1/time_series
@@ -8,18 +20,24 @@ mkdir -p -v results/$1/monthly
 mkdir -p -v results/$1/daily
 
 # Modifying driver file
-sed -i '2d' driver_file.dat #Linux
-#sed -i '' '2d' driver_file.dat #OS X
-sed -i "\$a setup/$1_driver.dat" driver_file.dat  #Linux
-#sed -i '' '$ a \
-#	setup/'$1'_driver.dat' driver_file.dat #OS X
+if [[ "$OS" == "linux" ]]; then
+	sed -i '2d' driver_file.dat #Linux
+	sed -i "\$a setup/$1_driver.dat" driver_file.dat  #Linux
+elif [[ "$OS" == "OSX" ]]; then
+	sed -i '' '2d' driver_file.dat #OS X
+	sed -i '' '$ a \
+		setup/'$1'_driver.dat' driver_file.dat #OS X
+fi
 
 # Modifying setup file
-sed -i '2d' setup_file.dat #Linux
-#sed -i '' '2d' setup_file.dat #OS X
-sed -i "\$a setup/$1_setup.dat" setup_file.dat #Linux
-#sed -i '' '$ a \
-#	setup/'$1'_setup.dat' setup_file.dat #OS X
+if [[ "$OS" == "linux" ]]; then
+	sed -i '2d' setup_file.dat #Linux
+	sed -i "\$a setup/$1_setup.dat" setup_file.dat #Linux
+elif [[ "$OS" == "OSX" ]]; then
+	sed -i '' '2d' setup_file.dat #OS X
+	sed -i '' '$ a \
+		setup/'$1'_setup.dat' setup_file.dat #OS X
+fi
 
 file=./setup/$1_setup.dat
 if [ ! -f ${file} ];
@@ -39,3 +57,4 @@ then
 	echo "Warning: The file ${file} does not exist"
 fi
 
+echo "Project for LAKE model created successfully"
diff --git a/source/model/carbon_dioxide.f90 b/source/model/carbon_dioxide.f90
index 58733c2a6218fdccd3acf8c4dc4756e6ed7758b0..b98bb9d493fa58c2c2c40a27f188a9f7775834b3 100644
--- a/source/model/carbon_dioxide.f90
+++ b/source/model/carbon_dioxide.f90
@@ -34,7 +34,7 @@ use T_SOLVER_MOD, only : &
 
 use NUMERICS, only : PROGONKA
 
-use NUMERIC_PARAMS, only : vector_length
+use NUMERIC_PARAMS, only : nveclen
 
 use DRIVING_PARAMS, only : carbon_model
 
@@ -100,8 +100,8 @@ real(kind=ireals) :: Flux_atm, Fcarbdi1 ! the flux of carbon dioxide at the bott
 real(kind=ireals) :: x, xx, xxx, xxxx
 integer(kind=iintegers) :: i, nspec, water_indic
 
-allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
-&         f(1:vector_length),y(1:vector_length))
+allocate (a(1:nveclen),b(1:nveclen),c(1:nveclen), &
+&         f(1:nveclen),y(1:nveclen))
 allocate (conc(1:gs%M+1))
 
 Fcarbdi1 = - sodbot*CO2O2_sod
@@ -178,10 +178,10 @@ do i = 1, nspec
   f(1)   = - x*conc(1) + Flux_atm + b(1)*conc(2) - xx*conc(1) - x*xxxx
   ! Specifying coefficients for internal layers
   if (i == 1) then
-    call DIFF_COEF(a,b,c,f,2,gs%M,2,gs%M,water_indic,dt)
+    call DIFF_COEF(nveclen,a,b,c,f,2,gs%M,2,gs%M,water_indic,dt)
   else
     workc => conc
-    call DIFF_COEF(a,b,c,f,2,gs%M,2,gs%M,water_indic,dt)
+    call DIFF_COEF(nveclen,a,b,c,f,2,gs%M,2,gs%M,water_indic,dt)
     nullify(workc)
   endif
   ! Tridiagonal system coefficients for the bottom layer
@@ -197,7 +197,7 @@ do i = 1, nspec
   f(gs%M+1) = - x*conc(gs%M+1) + a(gs%M+1)*conc(gs%M) - xx*conc(gs%M+1)
   ! Bottom diffusive flux only for DIC
   if (ls%ls1 == 0. .and. i == 1) f(gs%M+1) = f(gs%M+1) + Fcarbdi1
-  call PROGONKA (a,b,c,f,y,1,gs%M+1)
+  call PROGONKA (nveclen,a,b,c,f,y,1,gs%M+1)
   y(1:gs%M+1) = max(y(1:gs%M+1),0.d0)
 
   select case(i)
diff --git a/source/model/diffusion_mod.f90 b/source/model/diffusion_mod.f90
index 2a38f8be78b2c0b3f62d9a78167e470f8012acf3..cbc7712ede689be1b7551646b53ce069c1c4ad5c 100644
--- a/source/model/diffusion_mod.f90
+++ b/source/model/diffusion_mod.f90
@@ -1,6 +1,7 @@
 MODULE DIFFUSION_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
+use NUMERIC_PARAMS, only : nveclen
 
 implicit none
 
@@ -14,7 +15,7 @@ SUBROUTINE DIFFUSION(N,bctype,bcs,dt,h,dh,dh0, &
 use T_SOLVER_MOD, only : DIFF_COEF
 use ARRAYS_BATHYM,  only : bathym
 use NUMERICS, only : PROGONKA, STEP
-use NUMERIC_PARAMS, only : vector_length
+use NUMERIC_PARAMS, only : nveclen
 implicit none
 !Input variables
 integer(kind=iintegers), intent(in) :: N           !> Number of levels of finite-difference grid
@@ -42,8 +43,8 @@ real(kind=ireals), allocatable :: a(:), b(:), c(:), f(:), y(:)
 
 integer(kind=iintegers) :: i !Loop index
 
-allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
-&         f(1:vector_length),y(1:vector_length))
+allocate (a(1:nveclen),b(1:nveclen),c(1:nveclen), &
+&         f(1:nveclen),y(1:nveclen))
 a(:) = 0.; b(:) = 0.; c(:) = 0.; f(:) = 0.; y(:) = 0.
 
 dt_inv = 1./dt
@@ -92,7 +93,7 @@ elseif (bctype(2) == 2) then !Neumann condition
   & a(N)*var1(N-1) - xx*var1(N)
 endif
 
-call PROGONKA (a,b,c,f,var2,1_iintegers,N)
+call PROGONKA (nveclen,a,b,c,f,var2,1_iintegers,N)
 
 deallocate (a,b,c,f,y)
 END SUBROUTINE DIFFUSION
diff --git a/source/model/heatbalsurf_mod.f90 b/source/model/heatbalsurf_mod.f90
index 26ad8e54e7627d60c7f51b2fdb04b39a0e19a6a7..3f5eae9619cd48359873808af7021778dc48b5ec 100644
--- a/source/model/heatbalsurf_mod.f90
+++ b/source/model/heatbalsurf_mod.f90
@@ -128,12 +128,6 @@ FUNCTION HEATBALANCE(Tsurf,surftyp,RadWater,RadIce,fetch,dt)
  real(kind=ireals), intent(out) :: dHdt
  real(kind=ireals), intent(out) :: convect_flux
 
- !integer(kind=iintegers) itop 
- !real(kind=ireals) AL,DLT,DVT,ALLL,DL, &
- !& ALV,DV,Z,T,WL,WV,WI,dens,dz
- !common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
- !& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
- !common /SOILDAT/ dz(ms),itop
 
  SELECT CASE (surftyp)   
   case (soil_indic)
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index 8769ccfb10088f219f4f3cbc276be043847b1a5d..22599d797bf19134e5aef0c93b57e1131743a3b7 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -13,6 +13,7 @@ use PHYS_CONSTANTS
 use DRIVING_PARAMS 
 use ATMOS
 use ARRAYS
+use ARRAYS_GRID, only : gs
 use TURB_CONST, only : &
 & TURB_CONST_DERIVED_SET
 use INOUT_PARAMETERS, only : &
@@ -79,6 +80,7 @@ call PBLDAT_INM
 call COMSOILFORLAKE
 call TURB_CONST_DERIVED_SET
 call SET_DENS_CONSTS(eos%par)
+call NUMERIC_PARAMS_SET(gs)
 
 if (error_cov==1) then
   if (runmode%par == 2) then
@@ -418,7 +420,8 @@ real(kind=ireals) :: x, xx, y, yy, zz, zzz, vol_, vol_2, tt, ttt, ww, ww1, www
 real(kind=ireals) :: kor
 real(kind=ireals) :: Ti10 ! Initial temperature of the ice surface (if l10 > 0)
 
-real(kind=ireals), dimension(1:vector_length) :: a, b, c, d 
+!real(kind=ireals), dimension(1:nveclen) :: a, b, c, d 
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:) 
 
 real(kind=ireals) :: dt
 real(kind=ireals) :: b0
@@ -426,7 +429,8 @@ real(kind=ireals) :: tau_air, tau_i, tau_gr
 real(kind=ireals) :: eflux0_kinem_water
 real(kind=ireals) :: totmeth0, totmethc, metracc = 0.
 
-real(kind=ireals), dimension(1:vector_length) :: Temp
+!real(kind=ireals), dimension(1:nveclen) :: Temp
+real(kind=ireals), allocatable :: Temp(:)
 
 real(kind=ireals) :: flux1, flux2
 
@@ -466,6 +470,9 @@ data nstep_meas /0/
  
 SAVE
   
+allocate(a(1:nveclen),b(1:nveclen),c(1:nveclen), &
+&        d(1:nveclen),Temp(1:nveclen))
+
 !BODY OF PROGRAM
 
 ! Getting the number of cores
@@ -925,7 +932,7 @@ if (soilswitch%par == 1) then
 endif
 if (multsoil) then 
 ! Calculation of heat sources from heat fluxes from mutliple soil columns
-  call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+  call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                    & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                    & bathymwater,bathymice, &
                    & bathymdice,bathymsoil(1,ix,iy), &
@@ -938,7 +945,7 @@ endif
 call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
 & RadWater, RadIce, SurfTemp, fetch, dt)
 if (soilswitch%par == 1) then
-  call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+  call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
 endif
 
 ! Diagnostic calculations
@@ -968,7 +975,7 @@ call OXYGEN_PRODCONS(gs, gsp, wst, bathymsoil, gas, area_int, area_half, ddz05,
 if (soilswitch%par == 1) then
   ! Calculation of methane in all soil columns, except for the lowest one
   if (multsoil) then
-    call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
@@ -1298,7 +1305,7 @@ if (flag_snow == 1) then
   endif
   ! Calculation of heat sources from heat fluxes from mutliple soil columns
   if (multsoil) then
-    call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
@@ -1311,7 +1318,7 @@ if (flag_snow == 1) then
   call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
   & RadWater, RadIce, SurfTemp, fetch, dt)
   if (soilswitch%par == 1) then
-    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+    call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
   endif
 
   !Assuming all dissolved species diffuse with the same molecular diffusivity as salinity
@@ -1327,7 +1334,7 @@ if (flag_snow == 1) then
   if (soilswitch%par == 1) then
     ! Calculation of methane in all soil columns, except for the lowest one
     if (multsoil) then
-      call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                        & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                        & bathymwater,bathymice, &
                        & bathymdice,bathymsoil(1,ix,iy), &
@@ -1554,7 +1561,7 @@ else
   endif
   ! Calculation of heat sources from heat fluxes from mutliple soil columns
   if (multsoil) then 
-    call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
@@ -1567,7 +1574,7 @@ else
   call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
   & RadWater, RadIce, SurfTemp, fetch, dt)
   if (soilswitch%par == 1) then
-    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+    call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
   endif
 
   !Assuming all dissolved species diffuse with the same molecular diffusivity as salinity
@@ -1591,7 +1598,7 @@ else
   if (soilswitch%par == 1) then
   ! Calculation of temperature in all soil columns, except for the lowest one
     if (multsoil) then
-      call SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
                        & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                        & bathymwater,bathymice, &
                        & bathymdice,bathymsoil(1,ix,iy), &
@@ -1786,7 +1793,7 @@ if (flag_snow == 1) then
   call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
   & RadWater, RadIce, SurfTemp, fetch, dt)
   if (soilswitch%par == 1) then
-    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+    call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
   endif
 
   ! Salinity profile
@@ -1863,7 +1870,7 @@ else
   call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
   & RadWater, RadIce, SurfTemp, fetch, dt)
   if (soilswitch%par == 1) then
-    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+    call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
   endif
 
   ! Salinity profile
@@ -1912,7 +1919,7 @@ if4: IF (layer_case==4) THEN
     call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
     & RadWater, RadIce, SurfTemp, fetch, dt)
     if (soilswitch%par == 1) then
-      call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+      call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
     endif
 
     ! Salinity profile
@@ -1949,7 +1956,7 @@ if4: IF (layer_case==4) THEN
     ! Salnity profile
     call S_DIFF(dt,ls,salice)
     if (soilswitch%par == 1) then
-      call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
+      call SOILFORLAKE(nveclen,dt,a,b,c,d,nsoilcols)
       call BUBBLE_BLOCK 
     endif
     gs%isoilcol = nsoilcols
@@ -2895,6 +2902,8 @@ endif
 !if (l1 > 0.) print*, T(itop), Ti1(1), Ti1(Mice+1)
 
 deallocate(radflux)
+deallocate(a,b,c,d)
+deallocate(Temp)
 
 if (firstcall) firstcall = .false.
 
diff --git a/source/model/lake_modules.f90 b/source/model/lake_modules.f90
index 1fde2a0ec85e27c5201305d8f8da1fd776b71745..a63a5265b1eabe08c3ff5c1f316a955a2d436863 100644
--- a/source/model/lake_modules.f90
+++ b/source/model/lake_modules.f90
@@ -24,6 +24,45 @@ real(kind=ireals), parameter :: m_to_nm = 1.E+9
 END MODULE UNITS
 
 
+MODULE ARRAYS_GRID
+
+use LAKE_DATATYPES, only : ireals, iintegers
+!use DRIVING_PARAMS
+
+implicit none
+
+! Grid size group 
+integer(kind=iintegers), save, target :: nsoilcols ! Number of soil columns per one lake
+type, public :: gridsize_type
+  sequence
+  integer(kind=iintegers), pointer :: M, Mice, ns, ms, ml, nsoilcols, nx, ny
+  integer(kind=iintegers) :: ix, iy, isoilcol
+end type
+type(gridsize_type) :: gs
+
+! Grid spacing group
+real(kind=ireals), allocatable, target :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
+real(kind=ireals), allocatable, target :: z_full_ice(:), z_half_ice(:)
+real(kind=ireals), allocatable, target :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
+real(kind=ireals), allocatable, target :: ddzi(:), ddzi05(:)
+real(kind=ireals), allocatable, target :: dzeta_int(:), dzeta_05int(:)
+real(kind=ireals), allocatable, target :: dzetai_int(:), dzetai_05int(:)      
+integer(kind=iintegers), allocatable, target :: isoilcolc(:), ksoilcol(:)
+type, public :: gridspacing_type
+  sequence
+  real(kind=ireals), pointer :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
+  real(kind=ireals), pointer :: z_full_ice(:), z_half_ice(:)
+  real(kind=ireals), pointer :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
+  real(kind=ireals), pointer :: ddzi(:), ddzi05(:)
+  real(kind=ireals), pointer :: dzeta_int(:), dzeta_05int(:)
+  real(kind=ireals), pointer :: dzetai_int(:), dzetai_05int(:) 
+  integer(kind=iintegers), pointer :: isoilcolc(:), ksoilcol(:)
+end type
+type(gridspacing_type) :: gsp
+
+END MODULE ARRAYS_GRID
+
+
 MODULE NUMERIC_PARAMS
 
 use LAKE_DATATYPES, only : ireals, iintegers
@@ -34,7 +73,8 @@ integer(kind=iintegers), target, save :: ML_ = ML, MS_ = MS
 real(kind=ireals), parameter :: dznorm = 0.04, dzmin = 0.015, &
 & UpperLayer = 0.02 ! UpperLayer = 0.1, UpperLayer = 0.008 
 
-integer(kind=iintegers), parameter :: vector_length = 250
+!integer(kind=iintegers), parameter :: nveclen = 250
+integer(kind=iintegers), protected :: nveclen
 real(kind=ireals),    parameter :: pi = 4.*datan(1.d0)
 
 real(kind=ireals), parameter :: small_value = 1.d-20
@@ -52,6 +92,15 @@ real(kind=ireals), parameter :: minwind = 1.d-2
  !   depth --- depth of the lowest layer in soil (m)
  !   UpperLayer --- upper soil layer depth (m)
 
+contains
+!> Subroutine NUMERIC_PARAMS_SET assigns values to parameters related to model numerics
+SUBROUTINE NUMERIC_PARAMS_SET(gs)
+use ARRAYS_GRID, only : gridsize_type
+implicit none
+type(gridsize_type), intent(in) :: gs
+nveclen = max(gs%ms+gs%Mice+gs%ns,gs%M+gs%ns)
+END SUBROUTINE NUMERIC_PARAMS_SET
+
 END MODULE NUMERIC_PARAMS
 
 
@@ -433,47 +482,6 @@ implicit none
 END MODULE ARRAYS_WATERSTATE
 
 
-
-MODULE ARRAYS_GRID
-
-use LAKE_DATATYPES, only : ireals, iintegers
-!use DRIVING_PARAMS
-
-implicit none
-
-! Grid size group 
-integer(kind=iintegers), save, target :: nsoilcols ! Number of soil columns per one lake
-type, public :: gridsize_type
-  sequence
-  integer(kind=iintegers), pointer :: M, Mice, ns, ms, ml, nsoilcols, nx, ny
-  integer(kind=iintegers) :: ix, iy, isoilcol
-end type
-type(gridsize_type) :: gs
-
-! Grid spacing group
-real(kind=ireals), allocatable, target :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
-real(kind=ireals), allocatable, target :: z_full_ice(:), z_half_ice(:)
-real(kind=ireals), allocatable, target :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
-real(kind=ireals), allocatable, target :: ddzi(:), ddzi05(:)
-real(kind=ireals), allocatable, target :: dzeta_int(:), dzeta_05int(:)
-real(kind=ireals), allocatable, target :: dzetai_int(:), dzetai_05int(:)      
-integer(kind=iintegers), allocatable, target :: isoilcolc(:), ksoilcol(:)
-type, public :: gridspacing_type
-  sequence
-  real(kind=ireals), pointer :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
-  real(kind=ireals), pointer :: z_full_ice(:), z_half_ice(:)
-  real(kind=ireals), pointer :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
-  real(kind=ireals), pointer :: ddzi(:), ddzi05(:)
-  real(kind=ireals), pointer :: dzeta_int(:), dzeta_05int(:)
-  real(kind=ireals), pointer :: dzetai_int(:), dzetai_05int(:) 
-  integer(kind=iintegers), pointer :: isoilcolc(:), ksoilcol(:)
-end type
-type(gridspacing_type) :: gsp
-
-END MODULE ARRAYS_GRID
-
-
-
 MODULE ARRAYS_METHANE
 
 use LAKE_DATATYPES, only : ireals, iintegers
diff --git a/source/model/methane_mod.f90 b/source/model/methane_mod.f90
index 3f0e9ceda9ec30878c403b8b73a73da5851c1247..69d924fe88e607d0d165b5babd57630c89f1dc14 100644
--- a/source/model/methane_mod.f90
+++ b/source/model/methane_mod.f90
@@ -238,8 +238,8 @@ implicit none
  allocate (z(1:gs%ns))
  allocate (roots(1:gs%ns))
  allocate (pressoil(1:gs%ns))
- allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
- &         f(1:vector_length),y(1:vector_length))
+ allocate (a(1:nveclen),b(1:nveclen),c(1:nveclen), &
+ &         f(1:nveclen),y(1:nveclen))
  a(:) = 0.; b(:) = 0.; c(:) = 0.; f(:) = 0.; y(:) = 0.
 
  if (soilswitch%par == 0 .and. botflux_ch4 == missing_value) then
@@ -638,7 +638,7 @@ if (deepestsoil) then
     b(1)   = xx
     f(1)   = - yy*qwater(1,1) + Flux_atm + xx*(qwater(2,1) - qwater(1,1)) - &
     & lsm%water(1)*yy*dt 
-    call DIFF_COEF(a,b,c,f,2,gs%M,2,gs%M,water_methane_indic,dt)
+    call DIFF_COEF(nveclen,a,b,c,f,2,gs%M,2,gs%M,water_methane_indic,dt)
     water_bottom_bc_if : if (botflux_ch4 /= missing_value) then
       !Methane flux at the deepest point of a lake is set as b.c., coupling to soil is switched off
       xx = 0.5*( - bathymwater(gs%M)%area_half/bathymwater(gs%M+1)%area_int * &
@@ -647,7 +647,7 @@ if (deepestsoil) then
       c(gs%M+1) = xx - yy  ! bug corrected: - (dhw-dhw0)/(2.d0*dt) 
       a(gs%M+1) = xx 
       f(gs%M+1) = - yy*qwater(gs%M+1,1) + xx*(qwater(gs%M,1) - qwater(gs%M+1,1)) + botflux_ch4
-      call PROGONKA (a,b,c,f,y,1,gs%M+1)
+      call PROGONKA (nveclen,a,b,c,f,y,1,gs%M+1)
       y(1:gs%M+1) = max(y(1:gs%M+1),0.e0_ireals)
       qwater(1:gs%M+1,2) = y(1:gs%M+1)
     else
@@ -660,7 +660,7 @@ if (deepestsoil) then
         c(gs%M+1) = xx - yy  ! bug corrected: - (dhw-dhw0)/(2.d0*dt) 
         a(gs%M+1) = xx 
         f(gs%M+1) = - yy*qwater(gs%M+1,1) + xx*(qwater(gs%M,1) - qwater(gs%M+1,1)) ! + qflux1
-        call PROGONKA (a,b,c,f,y,1,gs%M+1)
+        call PROGONKA (nveclen,a,b,c,f,y,1,gs%M+1)
         y(1:gs%M+1) = max(y(1:gs%M+1),0.e0_ireals)
         qwater(1:gs%M+1,2) = y(1:gs%M+1)
         ! Methane evolution in soil under bottom ice
@@ -688,7 +688,7 @@ if (deepestsoil) then
         f(gs%ns) = qsoil(gs%ns) + dt*(methgenmh(gs%ns) + rprod(gs%ns) + &
         & bull(gs%ns)*(ch4_crit(gs%ns) - qsoil(gs%ns)) ) + &
         & a(gs%ns)*qsoil(gs%ns-1) - xx*qsoil(gs%ns)
-        call PROGONKA(a, b, c, f, y, 1, gs%ns)
+        call PROGONKA(nveclen,a, b, c, f, y, 1, gs%ns)
         y(1:gs%ns) = max(y(1:gs%ns),0.e0_ireals)
         qsoil(1:gs%ns) = y(1:gs%ns)     
       else
@@ -730,7 +730,7 @@ if (deepestsoil) then
           f(gs%M+gs%ns) = qsoil(gs%ns) + dt*(methgenmh(gs%ns) + rprod(gs%ns) + &
           & bull(gs%ns)*(ch4_crit(gs%ns) - qsoil(gs%ns))) + &
           & a(gs%M+gs%ns)*qsoil(gs%ns-1) - xx*qsoil(gs%ns)
-          call PROGONKA (a,b,c,f,y,1,gs%M+gs%ns)
+          call PROGONKA (nveclen,a,b,c,f,y,1,gs%M+gs%ns)
           y(1:gs%M+gs%ns) = max(y(1:gs%M+gs%ns),0.e0_ireals)
           qwater(1:gs%M+1,2) = y(1:gs%M+1)
           qsoil(2:gs%ns) = y(gs%M+2:gs%M+gs%ns)
@@ -794,7 +794,7 @@ if (deepestsoil) then
           f(gs%M+1+gs%ns) = qsoil(gs%ns) + dt*(methgenmh(gs%ns) + rprod(gs%ns) + &
           & bull(gs%ns)*(ch4_crit(gs%ns) - qsoil(gs%ns))) + &
           & a(gs%M+1+gs%ns)*qsoil(gs%ns-1) - xx*qsoil(gs%ns)
-          call PROGONKA (a,b,c,f,y,1,gs%M+1+gs%ns)
+          call PROGONKA (nveclen,a,b,c,f,y,1,gs%M+1+gs%ns)
           y(1:gs%M+1+gs%ns) = max(y(1:gs%M+1+gs%ns),0.e0_ireals)
           qwater(1:gs%M+1,2) = y(1:gs%M+1)
           qsoil(1:gs%ns) = y(gs%M+2:gs%M+1+gs%ns)
@@ -835,7 +835,7 @@ if (deepestsoil) then
     f(gs%ns) = qsoil(gs%ns) + dt*(methgenmh(gs%ns) + rprod(gs%ns) + &
     & bull(gs%ns)*(ch4_crit(gs%ns) - qsoil(gs%ns))) + &
     & a(gs%ns)*qsoil(gs%ns-1) - xx*qsoil(gs%ns)
-    call PROGONKA(a, b, c, f, y, 1, gs%ns)
+    call PROGONKA(nveclen,a, b, c, f, y, 1, gs%ns)
     y(1:gs%ns) = max(y(1:gs%ns),0.e0_ireals)
     qsoil(1:gs%ns) = y(1:gs%ns)     
   else ! bare soil
@@ -880,7 +880,7 @@ else
   f(gs%ns) = qsoil(gs%ns) + dt*(methgenmh(gs%ns) + rprod(gs%ns) + &
   & bull(gs%ns)*(ch4_crit(gs%ns) - qsoil(gs%ns))) + &
   & a(gs%ns)*qsoil(gs%ns-1) - xx*qsoil(gs%ns)
-  call PROGONKA(a, b, c, f, y, 1, gs%ns)
+  call PROGONKA(nveclen,a, b, c, f, y, 1, gs%ns)
   y(1:gs%ns) = max(y(1:gs%ns),0.e0_ireals)
   qsoil(1:gs%ns) = y(1:gs%ns)
   if (.not.ifflux) then
diff --git a/source/model/momentum_mod.f90 b/source/model/momentum_mod.f90
index 62b6d26a73585bdbc93c04e476aa4c1a458ea793..f06f8322efcfe46315bfe31bf125cf5077fff980 100644
--- a/source/model/momentum_mod.f90
+++ b/source/model/momentum_mod.f90
@@ -2,7 +2,7 @@ MODULE MOMENTUM_LAKE_MOD
 
 use NUMERICS, only : PROGONKA, KRON, MATRIXPROGONKA
 use TURB_LAKE_MOD, only : CE_CANUTO, SMOMENT_GALPERIN
-use NUMERIC_PARAMS, only : vector_length, pi, small_value, minwind
+use NUMERIC_PARAMS, only : nveclen, pi, small_value, minwind
 
 use LAKE_DATATYPES, only : ireals, iintegers
 
@@ -128,16 +128,21 @@ real(kind=ireals), intent(out) :: tau_air, tau_gr, tau_i, tau_wav
 
 real(kind=ireals) :: CE(M)
 
-real(kind=ireals) :: a(vector_length)
-real(kind=ireals) :: b(vector_length)
-real(kind=ireals) :: c(vector_length)
-real(kind=ireals) :: d(vector_length)
+!real(kind=ireals) :: a(nveclen)
+!real(kind=ireals) :: b(nveclen)
+!real(kind=ireals) :: c(nveclen)
+!real(kind=ireals) :: d(nveclen)
 
-real(kind=ireals) :: am_(vector_length,2,2)
-real(kind=ireals) :: bm_(vector_length,2,2)
-real(kind=ireals) :: cm_(vector_length,2,2)
-real(kind=ireals) :: ym_(vector_length,2)
-real(kind=ireals) :: dm_(vector_length,2)
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:)
+
+!real(kind=ireals) :: am_(nveclen,2,2)
+!real(kind=ireals) :: bm_(nveclen,2,2)
+!real(kind=ireals) :: cm_(nveclen,2,2)
+!real(kind=ireals) :: ym_(nveclen,2)
+!real(kind=ireals) :: dm_(nveclen,2)
+
+real(kind=ireals), allocatable :: am_(:,:,:), bm_(:,:,:), cm_(:,:,:), &
+&                                 ym_(:,:), dm_(:,:)
 
 !real(kind=ireals) :: work(1:M+1,1:2) ! Work array
 
@@ -244,6 +249,11 @@ xbot = botfric%par - 1 !Switch for sloping bottom drag law
 allocate(Force_x(1:M+1), Force_y(1:M+1))
 allocate(um(1:M+1),vm(1:M+1))
 allocate (LxLyCDL(1:M+1,1:2)); LxLyCDL(1:M+1,1:2) = missing_value
+allocate(a(1:nveclen),b(1:nveclen),c(1:nveclen),d(1:nveclen))
+allocate(am_(1:nveclen,1:2,1:2),bm_(1:nveclen,1:2,1:2), &
+&        cm_(1:nveclen,1:2,1:2), &
+&        dm_(1:nveclen,1:2),ym_(1:nveclen,1:2))
+
 
 u = u1(1)
 v = v1(1)
@@ -1336,7 +1346,7 @@ enddo
 ! endif  
 
 
-call MATRIXPROGONKA(am_,bm_,cm_,dm_,ym_,M+1)
+call MATRIXPROGONKA(nveclen,am_,bm_,cm_,dm_,ym_,M+1)
 do k = 1, M+1
   u2(k) = ym_(k,1)
   v2(k) = ym_(k,2)
@@ -1355,6 +1365,8 @@ enddo
 
 deallocate(Force_x, Force_y, um, vm)
 if (allocated(LxLyCDL)) deallocate(LxLyCDL)
+deallocate(a,b,c,d)
+deallocate(am_,bm_,cm_,dm_,ym_)
 
 if (firstcall) firstcall=.false.
 
diff --git a/source/model/numerics_mod.f90 b/source/model/numerics_mod.f90
index 23ea5ae54704a4274e8e42ca6e0cd8412591b51a..8f9f0179d6fcccacb5fe7389432b234d7ed67f02 100644
--- a/source/model/numerics_mod.f90
+++ b/source/model/numerics_mod.f90
@@ -1,15 +1,16 @@
 MODULE NUMERICS
 
 use LAKE_DATATYPES, only : ireals, iintegers
-use NUMERIC_PARAMS, only : vector_length
+use NUMERIC_PARAMS, only : nveclen
 
 contains
-SUBROUTINE MATRIXSUM(a,b,c,k)
+SUBROUTINE MATRIXSUM(nn,a,b,c,k)
  implicit none
  
 !MATRIXES: C=A+B
- real(kind=ireals), dimension (vector_length,2,2)::a,b,c
- integer(kind=iintegers) k,j,i
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension (nn,2,2) :: a,b,c
+ integer(kind=iintegers) :: k,j,i
 
  do j=1,2
   do i=1,2
@@ -20,13 +21,14 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
  END SUBROUTINE MATRIXSUM
 
 
- SUBROUTINE MATRIXMULT(a,b,c,k)
+ SUBROUTINE MATRIXMULT(nn,a,b,c,k)
  implicit none
 
 !MATRIXES: C=A*B      
 
- real(kind=ireals), dimension (vector_length,2,2)::a,b,c
- integer(kind=iintegers) k
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension (nn,2,2) :: a,b,c
+ integer(kind=iintegers) :: k
 
  c(k,1,1)=a(k,1,1)*b(k,1,1)+a(k,1,2)*b(k,2,1)
  c(k,1,2)=a(k,1,1)*b(k,1,2)+a(k,1,2)*b(k,2,2)
@@ -36,14 +38,14 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
  END SUBROUTINE MATRIXMULT
 
 
- SUBROUTINE MATRIXMULTVECTOR(a,g,f,k)
+ SUBROUTINE MATRIXMULTVECTOR(nn,a,g,f,k)
  implicit none
 
 !MATRIX A, VECTORS g, f: Ag=f
 
- real(kind=ireals) a(vector_length,2,2),f(vector_length,2), &
- & g(vector_length,2)
- integer(kind=iintegers) k
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals) :: a(nn,2,2),f(nn,2), g(nn,2)
+ integer(kind=iintegers) :: k
 
  f(k,1)=a(k,1,1)*g(k,1)+a(k,1,2)*g(k,2)
  f(k,2)=a(k,2,1)*g(k,1)+a(k,2,2)*g(k,2)
@@ -52,24 +54,24 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
  END SUBROUTINE MATRIXMULTVECTOR
 
 
- SUBROUTINE VECTORSUM(a,b,c,k)
+ SUBROUTINE VECTORSUM(nn,a,b,c,k)
  implicit none
 
 !VECTORS: C=A+B
-
- real(kind=ireals), dimension(vector_length,2)::a,b,c
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension(nn,2) :: a,b,c
  integer(kind=iintegers) k
  c(k,1)=a(k,1)+b(k,1)
  c(k,2)=a(k,2)+b(k,2)
  END SUBROUTINE VECTORSUM
 
 
- SUBROUTINE INVERSMATRIX(a,a1,k)
+ SUBROUTINE INVERSMATRIX(nn,a,a1,k)
  implicit none 
  
 !MATRIXES: A1=A*(-1)
-
- real(kind=ireals), dimension(vector_length,2,2)::a,a1
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension(nn,2,2) :: a,a1
  integer(kind=iintegers) k
  
  a1(k,1,1)=a(k,2,2)/(a(k,1,1)*a(k,2,2)-a(k,1,2)*a(k,2,1)) 
@@ -80,34 +82,35 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
  END SUBROUTINE INVERSMATRIX
 
  
- SUBROUTINE MATRIXPROGONKA(a,b,c,d,y,N)
+ SUBROUTINE MATRIXPROGONKA(nn,a,b,c,d,y,N)
  
 !MATRIXPROGONKA solves the set of MATRIX three-point diference equations 
  
  implicit none
 
- real(kind=ireals), dimension(vector_length,2,2):: a,b,c,x3,x32,x31,x4,alpha
- real(kind=ireals), dimension(vector_length,2):: y,d,x2,beta,x21
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension(nn,2,2) :: a,b,c,x3,x32,x31,x4,alpha
+ real(kind=ireals), dimension(nn,2) :: y,d,x2,beta,x21
  integer(kind=iintegers) N,i,j,k
 
- call INVERSMATRIX(c,x32,1)
- call MATRIXMULT(x32,b,x3,1)
+ call INVERSMATRIX(nn,c,x32,1)
+ call MATRIXMULT(nn,x32,b,x3,1)
  do j = 1, 2
    do i = 1, 2
      alpha(2,i,j) = x3(1,i,j)
    enddo
  enddo
- call MATRIXMULTVECTOR(x32,d,x2,1)
+ call MATRIXMULTVECTOR(nn,x32,d,x2,1)
  do i = 1, 2
    beta(2,i) = x2(1,i)
  enddo
 
  do k = 3, N
-  CALL MATRIXMULT(A,ALPHA,X3,k-1)
+  CALL MATRIXMULT(nn,A,ALPHA,X3,k-1)
   X4(k-1,1:2,1:2) = - X3(k-1,1:2,1:2)
-  CALL MATRIXSUM(C,X4,X31,k-1)
-  CALL INVERSMATRIX(X31,X32,k-1)
-  CALL MATRIXMULT(X32,B,X3,k-1)
+  CALL MATRIXSUM(nn,C,X4,X31,k-1)
+  CALL INVERSMATRIX(nn,X31,X32,k-1)
+  CALL MATRIXMULT(nn,X32,B,X3,k-1)
   do j = 1, 2
     do i = 1, 2
       alpha(k,i,j) = X3(k-1,i,j)
@@ -115,28 +118,28 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
   enddo
   !call matrixmult(x3,x31,x33,k-1)
   !call matrixsum(-b,x33,x3,k-1)
-  CALL MATRIXMULTVECTOR(A,BETA,X2,K-1)
-  CALL VECTORSUM(D,X2,X21,K-1)
-  CALL MATRIXMULTVECTOR(X32,X21,X2,K-1)
+  CALL MATRIXMULTVECTOR(nn,A,BETA,X2,K-1)
+  CALL VECTORSUM(nn,D,X2,X21,K-1)
+  CALL MATRIXMULTVECTOR(nn,X32,X21,X2,K-1)
   do i = 1, 2
     beta(k,i) = X2(k-1,i)
   enddo
  enddo
 
- CALL MATRIXMULT(A,ALPHA,X3,N)
+ CALL MATRIXMULT(nn,A,ALPHA,X3,N)
  X4(N,1:2,1:2) = - X3(N,1:2,1:2)
- CALL MATRIXSUM(C,X4,X31,N)
- CALL INVERSMATRIX(X31,X32,N)
- CALL MATRIXMULTVECTOR(A,BETA,X2,N)
- CALL VECTORSUM(D,X2,X21,N)
- CALL MATRIXMULTVECTOR(X32,X21,X2,N)
+ CALL MATRIXSUM(nn,C,X4,X31,N)
+ CALL INVERSMATRIX(nn,X31,X32,N)
+ CALL MATRIXMULTVECTOR(nn,A,BETA,X2,N)
+ CALL VECTORSUM(nn,D,X2,X21,N)
+ CALL MATRIXMULTVECTOR(nn,X32,X21,X2,N)
  do i = 1, 2
    Y(N,i) = X2(N,i)
  enddo
 
  do k = N-1, 1, -1
-  CALL MATRIXMULTVECTOR(ALPHA,Y,X2,K+1)
-  CALL VECTORSUM(X2,BETA,X21,K+1)
+  CALL MATRIXMULTVECTOR(nn,ALPHA,Y,X2,K+1)
+  CALL VECTORSUM(nn,X2,BETA,X21,K+1)
   Y(K,1) = X21(K+1,1)
   Y(K,2) = X21(K+1,2)
  enddo
@@ -154,11 +157,12 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
  END FUNCTION 
 
 
- SUBROUTINE IND_STAB_FACT_DB (a,b,c,N,M,ind_stab,ind_bound)
+ SUBROUTINE IND_STAB_FACT_DB (nn,a,b,c,N,M,ind_stab,ind_bound)
 
  implicit none
 
- real(kind=ireals), dimension(1:vector_length):: a,b,c
+ integer(kind=iintegers), intent(in) :: nn
+ real(kind=ireals), dimension(1:nn) :: a,b,c
  integer(kind=iintegers) M,i,N
  logical ind_stab, ind_bound 
 
@@ -213,7 +217,7 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
 
 
 
- SUBROUTINE PROGONKA(a, b, c, f, y, K, N)
+ SUBROUTINE PROGONKA(nn, a, b, c, f, y, K, N)
  implicit none
 !FACTORIZATION METHOD FOR THE FOLLOWING SYSTEM OF LINEAR EQUATIONS:
 !-a(i)*y(i-1)+c(i)*y(i)-b(i)*y(i+1)=f(i) i=K+1,N-1
@@ -221,12 +225,12 @@ SUBROUTINE MATRIXSUM(a,b,c,k)
 !-a(N)*y(N-1)+c(N)*y(N)=f(N)
 !
 
+ integer(kind=iintegers), intent(in) :: nn
  integer(kind=iintegers), intent(in) :: K, N
- real(kind=ireals), intent(in) :: a(vector_length), b(vector_length), &
- & c(vector_length), f(vector_length) 
- real(kind=ireals), intent(out) :: y(vector_length)
+ real(kind=ireals), intent(in) :: a(nn), b(nn), c(nn), f(nn) 
+ real(kind=ireals), intent(out) :: y(nn)
 
- real(kind=ireals) :: alpha(vector_length+2), beta(vector_length+2) 
+ real(kind=ireals) :: alpha(nn+2), beta(nn+2) 
  integer(kind=iintegers) :: i
 
  SAVE
diff --git a/source/model/out_mod.f90 b/source/model/out_mod.f90
index 6d612fa068e33849e41db0c053257448dda6a6c8..65854ac70045656b71667cd56db6f690d63cab8f 100644
--- a/source/model/out_mod.f90
+++ b/source/model/out_mod.f90
@@ -912,11 +912,7 @@ character :: format_char*100, workchar*10
 logical, allocatable :: firstcall(:,:)
 
 
-!real(kind=ireals) :: dz(ms)
 integer(kind=iintegers) :: i
-!common /SOILDAT/ dz,itop
-!real(kind=ireals) :: Tsn(1:ms), cs(1:ms)
-!common /SNOW_CHAR/ Tsn,cs
 
       
 SAVE
diff --git a/source/model/oxygen_mod.f90 b/source/model/oxygen_mod.f90
index 08031111ac94d5520bb739a6b042eb89bbb75bb0..dc1816b950bd01a744cec5bc928158696e2ae249 100644
--- a/source/model/oxygen_mod.f90
+++ b/source/model/oxygen_mod.f90
@@ -1,7 +1,7 @@
 MODULE OXYGEN_MOD
 
 use NUMERICS, only : PROGONKA, STEP
-use NUMERIC_PARAMS, only : vector_length, small_value
+use NUMERIC_PARAMS, only : nveclen, small_value
 use LAKE_DATATYPES, only : ireals, iintegers
 use ARRAYS_GRID, only : gridsize_type
 use ARRAYS_BATHYM, only : layers_type, bathym 
@@ -99,8 +99,8 @@ integer(kind=iintegers) :: i
 integer(kind=iintegers), parameter :: switchbotbc = 1 !1 - Neumann, 2 - Dirichlet
 integer(kind=iintegers) :: bctype(1:2)
 
-allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
-&         f(1:vector_length),y(1:vector_length))
+allocate (a(1:nveclen),b(1:nveclen),c(1:nveclen), &
+&         f(1:nveclen),y(1:nveclen))
 a(:) = 0.; b(:) = 0.; c(:) = 0.; f(:) = 0.; y(:) = 0.
 
 Foxyg1 = sodbot
diff --git a/source/model/phys_func.f90 b/source/model/phys_func.f90
index b64b795429f964cc4b0cf3f3e657143ae2fc1cba..f8035e73c3b54a1dc3bb0274d3a8c68c14c75819 100644
--- a/source/model/phys_func.f90
+++ b/source/model/phys_func.f90
@@ -560,8 +560,6 @@ END FUNCTION GSW_RHO
 ! implicit none
 !  real(kind=ireals) cloud
 !  real(kind=ireals) dirdif0,b,S0,sinh0
-!  common /cloud/ cloud
-! data cloud /0./
 
 !  SAVE
 
diff --git a/source/model/salinity_mod.f90 b/source/model/salinity_mod.f90
index a822cef4dc5889e5c511cf6f704ab4c391d9a798..956fd303e14a36c4aa3ba07b0fffa46a318d9a9c 100644
--- a/source/model/salinity_mod.f90
+++ b/source/model/salinity_mod.f90
@@ -1,7 +1,8 @@
 MODULE SALINITY_LAKE_MOD
 
-use NUMERICS, only : PROGONKA, CHECK_PROGONKA, IND_STAB_FACT_DB
-use NUMERIC_PARAMS, only : vector_length, small_value
+use NUMERICS, only : PROGONKA, CHECK_PROGONKA, IND_STAB_FACT_DB, &
+&                    nveclen        
+use NUMERIC_PARAMS, only : nveclen, small_value
 use LAKE_DATATYPES, only : ireals, iintegers
 use DZETA_MOD, only : VARMEAN
 use ARRAYS_BATHYM, only : layers_type
@@ -52,10 +53,13 @@ real(kind=ireals), intent(in) :: salice(1:Mice+1) !> Ice layer salinity, kg/kg
 
 ! Local variables
 real(kind=ireals) :: Sflux1,Sflux_soil_bot,x,y,z,zz=0.
-real(kind=ireals), dimension(1:vector_length):: a,b,c,d,Sal
+!real(kind=ireals), dimension(1:nveclen):: a,b,c,d,Sal
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:), Sal(:)
 
 integer(kind=iintegers) :: i !loop index
 
+allocate(a(1:nveclen), b(1:nveclen), c(1:nveclen), &
+&        d(1:nveclen), Sal(1:nveclen))
 a(:) = 0.; b(:) = 0.; c(:) = 0.; d(:) = 0
 
 ! Sflux1 --- salinity flux at the bottom boundary       
@@ -105,7 +109,7 @@ if (ls%ls1 > 0.) Sflux1 = Sflux1 - ls%dhwls /dt*( salice0 - Sal1(M+1) )
 ! sedim = 0
 
 if (water == 1) then
-  call DIFF_COEF(a,b,c,d,2,M,2,M,water_salinity_indic,dt)
+  call DIFF_COEF(nveclen,a,b,c,d,2,M,2,M,water_salinity_indic,dt)
   x = 0.5*( - bathymwater(1)%area_half*lamsal(1)/(h1*ddz(1)* &
   & bathymwater(1)%area_int) + dhw0/(2.d0*dt) )
   c(1)   = x - ddz(1)*h1/(2.d0*dt) ! - dhw0/(2*dt) is wrong sign!
@@ -119,12 +123,12 @@ if (water == 1) then
     c(M+1) = x - ddz(M)*h1/(2.d0*dt) 
     a(M+1) = x
     d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1 - x*(Sal1(M+1) - Sal1(M))
-    call PROGONKA (a,b,c,d,Sal,1,M+1)
+    call PROGONKA (nveclen,a,b,c,d,Sal,1,M+1)
     Sal (1:M+1) = max(Sal(1:M+1),0._ireals)
     Sal2(1:M+1) = Sal(1:M+1)
     if (soilswitch%par == 1 .and. salsoil%par == 1) then
 !     Salinity diffusion in soil
-      call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
+      call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
       x = 0.5*dt*wsoil(1)/dzs(1)
       c(1) = - 1.d0 - x
       b(1) = x
@@ -135,7 +139,7 @@ if (water == 1) then
       d(ns) = - Sals1(ns,nsoilcols) + 2.d0*dt*Sflux_soil_bot/dzs(ns-1) - &
       & x*(Sals1(ns,nsoilcols) + Sals1(ns-1,nsoilcols))
       !print*, CHECK_PROGONKA(ns,a,b,c,d,Sal)
-      call PROGONKA (a,b,c,d,Sal,1,ns)
+      call PROGONKA (nveclen,a,b,c,d,Sal,1,ns)
       Sal(1:ns) = max(Sal(1:ns),0._ireals)
       Sals2(1:ns,nsoilcols) = Sal(1:ns)
     else
@@ -154,7 +158,7 @@ if (water == 1) then
       c(M+1) = - ( 0.5*(dzs(1) + ddz(M)*h1)/dt + z - x ) 
       d(M+1) = - Sal1(M+1)*0.5*(dzs(1) + ddz(M)*h1)/dt - &
       & x*(Sal1(M+1) - Sal1(M)) + z*Sal1(1) + z*y*Sals1(2,nsoilcols)
-      call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_salinity_indic,dt)
+      call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,M+2,M+ns-1,soil_salinity_indic,dt)
       x = 0.5*wsoil(ns-1)*dt/(dzs(ns-1))
       c(M+ns) = - 1.d0 + x
       a(M+ns) = - x
@@ -169,7 +173,7 @@ if (water == 1) then
       !read*
       !print*, 'd', d(1:M+ns)
       !read*
-      call PROGONKA (a,b,c,d,Sal,1,M+ns)
+      call PROGONKA (nveclen,a,b,c,d,Sal,1,M+ns)
       Sal(1:M+ns) = max(Sal(1:M+ns),0._ireals)
       Sals2(2:ns,nsoilcols) = Sal(M+2:M+ns)
       Sals2(1,nsoilcols) = Sal(M+1)/y !mind y!
@@ -180,7 +184,7 @@ if (water == 1) then
       c(M+1) = x - ddz(M)*h1/(2.d0*dt) 
       a(M+1) = x
       d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1 - x*(Sal1(M+1) - Sal1(M))
-      call PROGONKA (a,b,c,d,Sal,1,M+1)
+      call PROGONKA (nveclen,a,b,c,d,Sal,1,M+1)
       Sal (1:M+1) = max(Sal(1:M+1),0._ireals)
       Sal2(1:M+1) = Sal(1:M+1)
     endif
@@ -189,7 +193,7 @@ if (water == 1) then
 else
 ! Case of the soil and the ice above     
   if (soilswitch%par == 1 .and. salsoil%par == 1) then
-    call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
+    call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
     x = 0.5*dt*wsoil(1)/dzs(1)
     c(1) = - 1.d0 - x
     b(1) = x
@@ -200,7 +204,7 @@ else
     d(ns) = - Sals1(ns,nsoilcols) + 2.d0*dt*Sflux_soil_bot/dzs(ns-1) - &
     & x*(Sals1(ns,nsoilcols) + Sals1(ns-1,nsoilcols))
     !print*, CHECK_PROGONKA(ns,a,b,c,d,Sal)
-    call PROGONKA (a,b,c,d,Sal,1,ns)
+    call PROGONKA (nveclen,a,b,c,d,Sal,1,ns)
     Sal(1:ns) = max(Sal(1:ns),0._ireals)
     Sals2(1:ns,nsoilcols) = Sal(1:ns)
   else
@@ -208,6 +212,8 @@ else
   endif
 endif 
 
+deallocate(a,b,c,d,Sal)
+
 !print*, VARMEAN(Sal2,bathymwater,1_iintegers)
 
 !print*, 'Sals2', Sals2(:,nsoilcols)
@@ -327,8 +333,11 @@ real   (kind=ireals)   , allocatable :: salprev(:)
 integer(kind=iintegers) :: i, j !Loop indices
 integer(kind=iintegers) :: iter
 
-real(kind=ireals), dimension(1:vector_length):: a,b,c,d
+!real(kind=ireals), dimension(1:nveclen):: a,b,c,d
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:)
       
+allocate(a(1:nveclen), b(1:nveclen), &
+&        c(1:nveclen), d(1:nveclen))
 a(:) = 0.; b(:) = 0.; c(:) = 0.; d(:) = 0
 
 !dh = ls%dhihigh - row0_d_roi*(ls%snmeltice + ls%dhiimp) ! Melting of ice at the top, negative
@@ -347,7 +356,7 @@ allocate(salprev(1:Mice+1))
 forall(i=1:Mice+1) salprev(i) = wst%salice(i)
 
 
-call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,ice_sal_indic,dt)
+call DIFF_COEF(nveclen,a,b,c,d,2,Mice,2,Mice,ice_sal_indic,dt)
 !1. Zero ice salinity at top of ice
 !c(1) = 1.
 !b(1) = 0.
@@ -371,7 +380,7 @@ else
   c(Mice+1) = 1.
   d(Mice+1) = 0.
 endif
-call PROGONKA (a,b,c,d,wst%salice,1,Mice+1)
+call PROGONKA (nveclen,a,b,c,d,wst%salice,1,Mice+1)
 forall(i=1:Mice) wst%salice(i) = max(wst%salice(i), sal_min)
 
 ! Brine salt concentration at the ice bottom equals to the water salinity below
@@ -487,7 +496,7 @@ do j = 2, Mice
   endif ifadjust
 enddo
 
-
+deallocate(a, b, c, d)
 deallocate(salprev)
 
 contains
diff --git a/source/model/snowcalc.f90 b/source/model/snowcalc.f90
index fa1298caa96201e8d170b2952bba0f29ffa43485..f2ad0002398e15e1d3c7a35f2fe33292c28cecdc 100644
--- a/source/model/snowcalc.f90
+++ b/source/model/snowcalc.f90
@@ -53,28 +53,8 @@ contains
 !C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       
       integer(kind=iintegers) j,i,k,iyear,imonth,iday
-       
-      integer(kind=iintegers) nmonth,nyear,nday,nhour
-      
-      common /time/ nyear,nmonth,nday,nhour
-      !COMMON /SOILDAT/ dz(ms),itop
-      !COMMON /BL/      ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-      !&                ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW, &
-      !&                HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF, &
-      !&                ElatOld,HFold,PRSold,extinct(ms)
-      !COMMON /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
-      !&                ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML), &
-      !&                dens(ms)
-      COMMON /spin/ Cond,Phase
-
-      !common /watericesnowarr/ lams, q
+
       real(kind=ireals) dt
-      !real(kind=ireals) dz,q(110)
-      !real(kind=ireals) lams(ms)
-      !real(kind=ireals) ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-      !& ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW,HS,ES,TGRNEW,QGRNEW, &
-      !& WSNEW,SNNEW,RUNOFF,ElatOld,HFold,PRSold
-      !real(kind=ireals) AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
       real(kind=ireals) Erad_sc,ts,snow_sc, &
       & evap,erain,countrad,qbase,dmass,densev,counter,qin,rf, &
       & smelti,ein,dzi,ziz,porosityOfSnow,w0,fukt,dztop,smelt,eout,p1,q0, &
@@ -537,22 +517,8 @@ contains
 !C      snmelt, m/sec  --- mean snow melting speed during dt
 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
-      !COMMON /SOILDAT/ dz(ms),itop
-      !COMMON /BL/       ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-      !&                  ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW, &
-      !&                  HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF, &
-      !&                  ElatOld,HFold,PRSold,extinct(ms)
-      !COMMON /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
-      !&                  ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML), &
-      !&                  dens(ms)
-      COMMON /spin/ Cond,Phase
-      !real(kind=ireals) dz
-      !real(kind=ireals) ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-      !& ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW,HS,ES,TGRNEW,QGRNEW, &
-      !& WSNEW,SNNEW,RUNOFF,ElatOld,HFold,PRSold
-      !real(kind=ireals) AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
       real(kind=ireals) ts,snow,dz0,pmelt,pheat, &
-      & z0z,dmass, swe,densny,Cond,Phase
+      & z0z,dmass, swe,densny
       integer(kind=iintegers) j,i,k,mn,iyear,imonth,iday
       real(kind=ireals) T2,CCT(ML)
       real(kind=ireals) dt 
diff --git a/source/model/snowtemp.f90 b/source/model/snowtemp.f90
index 936e8f9a32c7ed1b17dd60c31b679efd8040cd61..39d4c7297f43df3353e22c1bb3bb709836848b66 100644
--- a/source/model/snowtemp.f90
+++ b/source/model/snowtemp.f90
@@ -52,31 +52,18 @@ real(kind=ireals) :: snowmass_init
 real(kind=ireals) :: CCT(ML)
 real(kind=ireals) :: dt
 
-!real(kind=ireals) :: ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-!& ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW, &
-!& HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF, &
-!& ElatOld,HFold,PRSold,extinct
-!real(kind=ireals) :: AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
-real(kind=ireals), dimension(1:vector_length) :: a, b, c, d, Temp
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:), Temp(:)
 
 integer(kind=iintegers) :: i
 integer(kind=iintegers) :: iyear
 integer(kind=iintegers) :: imonth
 integer(kind=iintegers) :: iday  
 
-
-!common /watericesnowarr/ lams, q
-!common /BL/ ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD,&
-!     & ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW,&
-!     & HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF,&
-!     & ElatOld,HFold,PRSold,extinct(ms)
-!common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML),&
-!     & ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
-!common /SOILDAT/ dz(ms),itop
-!common /snow_char/ Tsn,cs
-
 SAVE
 
+allocate(a(1:nveclen), b(1:nveclen), c(1:nveclen), &
+&        d(1:nveclen), Temp(1:nveclen))
+
 do i = itop,ms
   T(i) = Tsn(i)    
 enddo
@@ -104,6 +91,8 @@ totalmelts = totalmelts + snmelt*dt
 
 !if (hs1 == 0) hs1=0.00001
 
+deallocate(a, b, c, d, Temp)
+
 END SUBROUTINE SNOWTEMP
 
 
@@ -143,24 +132,8 @@ real(kind=ireals) :: dzsn(ms)
 real(kind=ireals) :: rofresh
 real(kind=ireals) :: CCT(ML)
 
-!real(kind=ireals) :: ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD, &
-!& ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW, &
-!& HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF, &
-!& ElatOld,HFold,PRSold,extinct
-!real(kind=ireals) :: AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
-
-!integer(kind=iintegers) :: itop
 integer(kind=iintegers) :: i
 
-!common /watericesnowarr/ lams, q
-!common /BL/ ST,PGR,TGROLD,QGROLD,RADIAT,WSOLD,SNOLD,&
-!     & ZS,thsoil,whsoil,BOLD,RF1,RF2,SNMELT,HSNOW,&
-!     & HS,ES,TGRNEW,QGRNEW,WSNEW,SNNEW,RUNOFF,&
-!     & ElatOld,HFold,PRSold,extinct(ms)
-!common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML),&
-!     & ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
-!common /SOILDAT/ dz(ms),itop
-!common /snow_char/ Tsn,cs
 
 if (flag_snow_init == 1) then
   itop = ms-2 
diff --git a/source/model/soil_mod.f90 b/source/model/soil_mod.f90
index e16063cb5f64621c1a92371d24a3430d0afae364..f434f90d7da0b4dffccd359d5ebd387d16a8306f 100644
--- a/source/model/soil_mod.f90
+++ b/source/model/soil_mod.f90
@@ -2,7 +2,7 @@ MODULE SOIL_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
 use NUMERICS, only : PROGONKA, STEP
-use NUMERIC_PARAMS, only : vector_length
+use NUMERIC_PARAMS, only : nveclen
 
 contains
 SUBROUTINE COMSOILFORLAKE
@@ -207,7 +207,7 @@ RETURN
 END SUBROUTINE COMSOILFORLAKE
 
 
-SUBROUTINE SOILFORLAKE(dt,a,b,c,d,is)
+SUBROUTINE SOILFORLAKE(nn,dt,a,b,c,d,is)
 
 !SOILFORLAKE calculates profiles of temperature, liquid and solid water 
 !content in the soil under a lake
@@ -234,6 +234,8 @@ use ATMOS!, only : &
 
 implicit none
 
+integer(kind=iintegers), intent(in) :: nn 
+
 real(kind=ireals), intent(in) :: dt
 
 !integer(kind=iintegers), intent(in) :: ix, iy
@@ -247,11 +249,11 @@ integer(kind=iintegers), intent(in) :: is ! Number of soil column
 !real(kind=ireals)   , intent(in) :: extwat, extice
 !real(kind=ireals)   , intent(in) :: fetch
     
-real(kind=ireals), intent(inout) :: a(1:vector_length)
-real(kind=ireals), intent(inout) :: b(1:vector_length)
-real(kind=ireals), intent(inout) :: c(1:vector_length)
-real(kind=ireals), intent(inout) :: d(1:vector_length)
-!real(kind=ireals) :: Temp(1:vector_length)
+real(kind=ireals), intent(inout) :: a(1:nn)
+real(kind=ireals), intent(inout) :: b(1:nn)
+real(kind=ireals), intent(inout) :: c(1:nn)
+real(kind=ireals), intent(inout) :: d(1:nn)
+!real(kind=ireals) :: Temp(1:nveclen)
 
 real(kind=ireals) :: dzmean
 real(kind=ireals) :: ma, mi
@@ -397,7 +399,7 @@ do i=2,ns-1
 !& lammoist1/dzs(i-1)*(wl1(i)-wl1(i-1)))/dzmean
 enddo   
     
-call PROGONKA(a,b,c,d,wl2,1,ns)  
+call PROGONKA(nveclen,a,b,c,d,wl2,1,ns)  
     
 if (h1 /= 0.and.ls1 == 0) then 
   lammoist1 = (lammoist(1)+lammoist(2))/2.
@@ -777,7 +779,7 @@ if (firstcall) firstcall = .false.
 END SUBROUTINE SOIL_COND_HEAT_COEF
 
 
-SUBROUTINE SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm, &
+SUBROUTINE SOILCOLSTEMP(nn,gs,gsp,dt,ls,ftot,ch4_pres_atm, &
                        & ddz,ddzi,zsoilcols, &
                        & wst,SR,a,b,c,d,add_to_winter, &
                        & bathymwater,bathymice, &
@@ -843,7 +845,8 @@ type(bathym),           intent(in) :: bathymdice(1:gs%Mice+1), bathymsoil(1:gs%n
 
 integer(kind=iintegers), intent(in) :: contr(1:2)
 
-real(kind=ireals), intent(inout), dimension(1:vector_length) :: a, b, c, d
+integer(kind=iintegers), intent(in) :: nn
+real(kind=ireals), intent(inout), dimension(1:nn) :: a, b, c, d
 
 logical, intent(inout) :: add_to_winter
 
@@ -904,7 +907,7 @@ if (contr(1) == 1) then
       & csoil(1)*rosoil(1)*dt_inv*0.5*dzs(1)*( Tsoil2(1) - Tsoil1(1,is) ) - &
       & 0.25* ( lamsoil(1) + lamsoil(2) ) * &
       & ( Tsoil2(2) + Tsoil1(2,is) - Tsoil2(1) - Tsoil1(1,is) )/dzs(1)
-      call SOILFORLAKE(dt,a,b,c,d,is)
+      call SOILFORLAKE(nveclen,dt,a,b,c,d,is)
     enddo
   elseif (bcswitch == 2) then ! continuity of flux, calculated from surface layer theory above soil surface
     do is = 1, gs%nsoilcols-1
@@ -913,7 +916,7 @@ if (contr(1) == 1) then
       soilflux(is) = LOGFLUX(sqrt(wst%u1(i)*wst%u1(i) + wst%v1(i)*wst%v1(i)), &
       & Tsoil1(1,is) - wst%Tw1(i), work1(is), z0_bot, z0_bot, cw_m_row0, 1) + SR(i)
       call T_DIFF(1,soilflux(is),dt,0,0,0,0,is,.true.)
-      call SOILFORLAKE(dt,a,b,c,d,is)
+      call SOILFORLAKE(nveclen,dt,a,b,c,d,is)
     enddo
   endif
 endif
diff --git a/source/model/surf_scheme_inm_mod.f90 b/source/model/surf_scheme_inm_mod.f90
index 830ff2c124ef6a5f8bed4bcea62fe25450f0ed36..cd0df5446cbbe6721086cef9c8128e54dfb320e3 100644
--- a/source/model/surf_scheme_inm_mod.f90
+++ b/source/model/surf_scheme_inm_mod.f90
@@ -43,7 +43,6 @@
       COMMON /HYDR/ SNCR,WLMMX,CEFF,CA,DZL,DZG,BMIN,HR,ZRM,ZRMM,TRM, &
      & FLXMIN,TOMIN
       COMMON /VEGSW/ CCT,TBEST,FFTMIN,FFQMIN,CK,SWW
-!C      COMMON /DRAIN/ WSDL,WSDG,DMIN,DMAX,D,WIINF
       !FBL(1)=SIG(NLEV)
       !FBL(2)=1.0E0
       !FBL(3)=1.0E0
diff --git a/source/model/t_solver_mod.f90 b/source/model/t_solver_mod.f90
index 2a2bf6e9e5cf35715a9e4d28e538ba586e9df9f9..b74d58efe21930d66ac25c8fb1de97502bcc682c 100644
--- a/source/model/t_solver_mod.f90
+++ b/source/model/t_solver_mod.f90
@@ -235,21 +235,16 @@ integer(kind=iintegers), intent(in) :: isoilcol
 logical, intent(in) :: ifflux
 
 ! Local variables
-!real(kind=ireals), dimension(1:ms) :: Tsn,lams,q,cs
-real(kind=ireals), dimension(1:vector_length) :: a,b,c,d,Temp
-!real(kind=ireals) :: AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
+!real(kind=ireals), dimension(1:nveclen) :: a,b,c,d,Temp
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:), Temp(:)
 real(kind=ireals) :: cc, xx
 integer(kind=iintegers) :: i,j
-   
-!common /snow_char/ Tsn,cs
-!common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
-!& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
-!common /SOILDAT/ dz(ms),itop
-!common /watericesnowarr/ lams,q
-!data num /1/
 
 SAVE
 
+allocate(a(1:nveclen), b(1:nveclen), c(1:nveclen), &
+&        d(1:nveclen), Temp(1:nveclen))
+
 !Meltpnt - Melting point temperature, C degrees 
 !Meltpnt=0.
 !ddz=1./float(M)VS,06.2007
@@ -260,7 +255,7 @@ if (surf == 1) then
     b(itop) = 0.
     d(itop) = Tsurf
     if (itop <= ms-2) then
-      call DIFF_COEF(a,b,c,d,itop+1,ms-1,itop+1,ms-1,snow_indic,dt) 
+      call DIFF_COEF(nveclen,a,b,c,d,itop+1,ms-1,itop+1,ms-1,snow_indic,dt) 
     endif
     if (ice == 1) then ! Snow and ice exist
 !----------------SNOW-ICE INTERFACE
@@ -273,12 +268,12 @@ if (surf == 1) then
       & a(ms)*T(ms-1) - (a(ms) + b(ms))*Ti1(1) + b(ms)*Ti1(2) - &
       & (snmeltice + dhiimp)*Lwi*row0*dt_inv
 !-----------------------------------
-      call DIFF_COEF(a,b,c,d,2,Mice,ms+1,ms+Mice-1,ice_indic,dt) 
+      call DIFF_COEF(nveclen,a,b,c,d,2,Mice,ms+1,ms+Mice-1,ice_indic,dt) 
       if (water == 1) then ! Snow, ice and water exist
         a(ms+Mice) = 0.
         c(ms+Mice) = 1.
         d(ms+Mice) = Meltpnt(Sal2(1),preswat(1),nmeltpoint%par)
-        call PROGONKA (a,b,c,d,Temp,itop,ms+Mice)
+        call PROGONKA (nveclen,a,b,c,d,Temp,itop,ms+Mice)
         do i = itop, ms
           Tsn(i) = Temp(i)
         enddo
@@ -298,11 +293,11 @@ if (surf == 1) then
           & a(ms+Mice)*Ti1(Mice) - (a(ms+Mice) + b(ms+Mice))*Ti1(Mice+1) + &
           & b(ms+Mice)*Tsoil1(2,isoilcol)
 !----------------------------------------------------------- 
-          call DIFF_COEF(a,b,c,d,2,ns-1,ms+Mice+1,ms+Mice+ns-2,soil_indic,dt,isoilcol)
+          call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,ms+Mice+1,ms+Mice+ns-2,soil_indic,dt,isoilcol)
           c(ms+Mice+ns-1) = 1.
           a(ms+Mice+ns-1) = 1.
           d(ms+Mice+ns-1) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1) + lamsoil(ns))
-          call PROGONKA (a,b,c,d,Temp,itop,ms+Mice+ns-1)
+          call PROGONKA (nveclen,a,b,c,d,Temp,itop,ms+Mice+ns-1)
           do i = ms+Mice, ms+Mice+ns-1
             Tsoil2(i-ms-Mice+1) = Temp(i)
           enddo
@@ -318,7 +313,7 @@ if (surf == 1) then
           ! soilbotflx is interpreted as top soil boundary flux in this case
           & (xx*0.5*lami_v(Mice)/(ddzi(Mice)*l1) - &
           & 0.25*ci_m_roi_v(Mice+1)*(dhi-dhi0)*dt_inv)*Ti1(Mice+1) + a(ms+Mice)*Ti1(Mice)
-          call PROGONKA (a,b,c,d,Temp,itop,ms+Mice)
+          call PROGONKA (nveclen,a,b,c,d,Temp,itop,ms+Mice)
           Tsoil2(:) = Tsoil1(:,isoilcol)
         endif
         do i = itop, ms
@@ -339,11 +334,11 @@ if (surf == 1) then
         & csoil(1)*rosoil(1)*dzs(1)) - SR_botsnow + &
         & a(ms)*T(ms-1) - (a(ms) + b(ms))*Tsoil1(1,isoilcol) + b(ms)*Tsoil1(2,isoilcol)
 !---------------------------------------------------------------
-        call DIFF_COEF(a,b,c,d,2,ns-1,ms+1,ms+ns-2,soil_indic,dt,isoilcol) 
+        call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,ms+1,ms+ns-2,soil_indic,dt,isoilcol) 
         c(ms+ns-1) = 1.
         a(ms+ns-1) = 1.
         d(ms+ns-1) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
-        call PROGONKA (a,b,c,d,Temp,itop,ms+ns-1)
+        call PROGONKA (nveclen,a,b,c,d,Temp,itop,ms+ns-1)
         do i = ms, ms+ns-1
           Tsoil2(i-ms+1) = Temp(i)
         enddo
@@ -353,7 +348,7 @@ if (surf == 1) then
         d(ms) = - T(ms)*0.5*cs(ms)*dens(ms)*dz(ms)*dt_inv - &
         & SR_botsnow + soilbotflx%par + & ! soilbotflx is interpreted as top soil boundary flux in this case
         & a(ms)*(T(ms-1) - T(ms))
-        call PROGONKA (a,b,c,d,Temp,itop,ms)
+        call PROGONKA (nveclen,a,b,c,d,Temp,itop,ms)
         Tsoil2(:) = Tsoil1(:,isoilcol)
       endif
       do i = itop, ms
@@ -364,12 +359,12 @@ if (surf == 1) then
     b(1) = 0.
     c(1) = 1.
     d(1) = Tsurf
-    call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,ice_indic,dt)  
+    call DIFF_COEF(nveclen,a,b,c,d,2,Mice,2,Mice,ice_indic,dt)  
     if (water == 1) then ! Ice and water exist
       a(Mice+1) = 0.
       c(Mice+1) = 1.
       d(Mice+1) = Meltpnt(Sal2(1),preswat(1),nmeltpoint%par)
-      call PROGONKA (a,b,c,d,Temp,1,Mice+1)
+      call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+1)
       do i = 1, Mice+1
         Ti2(i) = Temp(i)
       enddo
@@ -386,11 +381,11 @@ if (surf == 1) then
         & a(Mice+1)*Ti1(Mice) - (a(Mice+1) + b(Mice+1))*Ti1(Mice+1) + &
         & b(Mice+1)*Tsoil1(2,isoilcol)
 !-----------------------------------------------------------
-        call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol) 
+        call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol) 
         c(Mice+ns) = 1.
         a(Mice+ns) = 1.
         d(Mice+ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
-        call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+ns)
         do i = Mice+1, Mice+ns
           Tsoil2(i-Mice) = Temp(i)
         enddo
@@ -406,7 +401,7 @@ if (surf == 1) then
         & (xx*0.5*lami_v(Mice)/(ddzi(Mice)*l1) - &
         & 0.25*dt_inv*ci_m_roi_v(Mice+1)*(dhi - dhi0))*Ti1(Mice+1) + &
         & a(Mice+1)*Ti1(Mice)
-        call PROGONKA (a,b,c,d,Temp,1,Mice+1)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+1)
         Tsoil2(:) = Tsoil1(:,isoilcol)
       endif
       do i = 1, Mice+1
@@ -417,12 +412,12 @@ if (surf == 1) then
     b(1) = 0.
     c(1) = 1.
     d(1) = Tsurf
-    call DIFF_COEF(a,b,c,d,2,M,2,M,water_indic,dt) 
+    call DIFF_COEF(nveclen,a,b,c,d,2,M,2,M,water_indic,dt) 
     if (deepice == 1) then ! Water and deep ice exist
       a(M+1) = 0.
       c(M+1) = 1.
       d(M+1) = Meltpnt(Sal2(M+1),preswat(M+1),nmeltpoint%par)
-      call PROGONKA(a,b,c,d,Temp,1,M+1) 
+      call PROGONKA(nveclen,a,b,c,d,Temp,1,M+1) 
       do i = 1, M+1
         Tw2(i) = Temp(i)
       enddo
@@ -440,11 +435,11 @@ if (surf == 1) then
         & xx*cw_m_row0*PEMF(M)*pt_down_f(M) + &
         & a(M+1)*Tw1(M) - (a(M+1) + b(M+1))*Tw1(M+1) + b(M+1)*Tsoil1(2,isoilcol)
 !-----------------------------------------------------------
-        call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt,isoilcol)
+        call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt,isoilcol)
         c(M+ns) = 1.
         a(M+ns) = 1.
         d(M+ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
-        call PROGONKA (a,b,c,d,Temp,1,M+ns)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,M+ns)
         do i = M+1, M+ns
           Tsoil2(i-M) = Temp(i)
         enddo
@@ -459,14 +454,14 @@ if (surf == 1) then
           ! soilbotflx is interpreted as top soil boundary flux in this case
           & 0.5*( - xx*lamw(M)/(ddz(M)*h1) + 0.5*dt_inv*cw_m_row0*(dhw-dhw0)) * &
           & (Tw1(M) - Tw1(M+1))
-          call PROGONKA (a,b,c,d,Temp,1,M+1)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,M+1)
           Tsoil2(:) = Tsoil1(:,isoilcol)
         elseif (cuette%par == 1 .or. cuette%par == 11) then
           !B.c.s for Cuette flow
           a(M+1) = 0.
           c(M+1) = 1.
           d(M+1) = Tw1(M+1)
-          call PROGONKA (a,b,c,d,Temp,1,M+1)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,M+1)
           Tsoil2(:) = Tsoil1(:,isoilcol)
         endif
       endif
@@ -486,12 +481,12 @@ if (surf == 1) then
       c(1) = 1.
       d(1) = Tsurf
     endif
-    call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_indic,dt,isoilcol) 
+    call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,2,ns-1,soil_indic,dt,isoilcol) 
     c(ns) = 0.5
     a(ns) = 0.5
     d(ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns)) + &
     & a(ns)*Tsoil1(ns-1,isoilcol) - c(ns)*Tsoil1(ns,isoilcol)
-    call PROGONKA (a,b,c,d,Temp,1,ns)
+    call PROGONKA (nveclen,a,b,c,d,Temp,1,ns)
     do i = 1, ns
       Tsoil2(i) = Temp(i)
     enddo
@@ -504,19 +499,19 @@ else
       b(1) = 0.
       c(1) = 1.
       d(1) = Meltpnt(Sal2(1),preswat(1),nmeltpoint%par)
-      call DIFF_COEF(a,b,c,d,2,M,2,M,water_indic,dt)  
+      call DIFF_COEF(nveclen,a,b,c,d,2,M,2,M,water_indic,dt)  
       if (deepice == 1) then ! Ice, water and deep ice exist
         a(M+1) = 0.
         c(M+1) = 1.
         d(M+1) = Meltpnt(Sal2(M+1),preswat(M+1),nmeltpoint%par)
-        call PROGONKA (a,b,c,d,Temp,1,M+1)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,M+1)
         do i = 1, M+1
           Tw2(i) = Temp(i)
         enddo
         b(1) = 0.
         c(1) = 1.
         d(1) = Meltpnt(Sal2(M+1),preswat(M+1),nmeltpoint%par)
-        call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
+        call DIFF_COEF(nveclen,a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
         if (soilswitch%par == 1) then ! Ice, water, deep ice and soil exist
 !----------------DEEPICE-SOIL INTERFACE-------------------------
           xx = bathymdice(Mice+1)%area_int/bathymdice(Mice)%area_half
@@ -529,11 +524,11 @@ else
           & a(Mice+1)*Tis1(Mice) - (a(Mice+1) + b(Mice+1))*Tis1(Mice+1) &
           & + b(Mice+1)*Tsoil1(2,isoilcol)
 !-----------------------------------------------------------
-          call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol)  
+          call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol)  
           c(Mice+ns) = 1.
           a(Mice+ns) = 1.
           d(Mice+ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1) + lamsoil(ns))
-          call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+ns)
           do i = Mice+1, Mice+ns
             Tsoil2(i-Mice) = Temp(i)
           enddo
@@ -547,7 +542,7 @@ else
           ! soilbotflx is interpreted as top soil boundary flux in this case
           & (xx*0.5*lami/(ddzi(Mice)*ls1) - 0.25*dt_inv*ci_m_roi*(dls-dls0)) * &
           & (Tis1(Mice+1) - Tis1(Mice))
-          call PROGONKA (a,b,c,d,Temp,1,Mice+1)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+1)
           Tsoil2(:) = Tsoil1(:,isoilcol)
         endif
         do i = 1, Mice+1
@@ -567,11 +562,11 @@ else
           & xx*cw_m_row0*PEMF(M)*pt_down_f(M) + &
           & a(M+1)*Tw1(M) - (a(M+1) + b(M+1))*Tw1(M+1) + b(M+1)*Tsoil1(2,isoilcol)
 !-----------------------------------------------------------
-          call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt,isoilcol)
+          call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,M+2,M+ns-1,soil_indic,dt,isoilcol)
           c(M+ns) = 1.
           a(M+ns) = 1.
           d(M+ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
-          call PROGONKA (a,b,c,d,Temp,1,M+ns)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,M+ns)
           do i = M+1, M+ns
             Tsoil2(i-M) = Temp(i)
           enddo
@@ -585,7 +580,7 @@ else
           ! soilbotflx is interpreted as top soil boundary flux in this case
           & 0.5*( - xx*lamw(M)/(ddz(M)*h1) + 0.5*dt_inv*cw_m_row0*(dhw-dhw0)) * &
           & (Tw1(M) - Tw1(M+1))
-          call PROGONKA (a,b,c,d,Temp,1,M+1)
+          call PROGONKA (nveclen,a,b,c,d,Temp,1,M+1)
           Tsoil2(:) = Tsoil1(:,isoilcol)
         endif
         do i = 1, M+1
@@ -596,7 +591,7 @@ else
       b(1) = 0.
       c(1) = 1.
       d(1) = Meltpnt(Sal2(M+1),preswat(M+1),nmeltpoint%par)
-      call DIFF_COEF(a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
+      call DIFF_COEF(nveclen,a,b,c,d,2,Mice,2,Mice,deepice_indic,dt)
       if (soilswitch%par == 1) then ! Water, deep ice and soil, no top ice
 !----------------DEEPICE-SOIL INTERFACE-------------------------
         xx = bathymdice(Mice+1)%area_int/bathymdice(Mice)%area_half
@@ -609,11 +604,11 @@ else
         & a(Mice+1)*Tis1(Mice) - (a(Mice+1) + b(Mice+1))*Tis1(Mice+1) &
         & + b(Mice+1)*Tsoil1(2,isoilcol)
 !-----------------------------------------------------------
-        call DIFF_COEF(a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol)  
+        call DIFF_COEF(nveclen,a,b,c,d,2,ns-1,Mice+2,Mice+ns-1,soil_indic,dt,isoilcol)  
         c(Mice+ns) = 1.
         a(Mice+ns) = 1.
         d(Mice+ns) = - soilbotflx%par*dzs(ns-1)/(lamsoil(ns-1)+lamsoil(ns))
-        call PROGONKA (a,b,c,d,Temp,1,Mice+ns)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+ns)
         do i = Mice+1, Mice+ns
           Tsoil2(i-Mice) = Temp(i)
         enddo
@@ -627,7 +622,7 @@ else
         ! soilbotflx is interpreted as top soil boundary flux in this case
         & ( - xx*0.5*lami/(ddzi(Mice)*ls1) + 0.25*dt_inv*ci_m_roi*(dls-dls0))* &
         & (Tis1(Mice) - Tis1(Mice+1))
-        call PROGONKA (a,b,c,d,Temp,1,Mice+1)
+        call PROGONKA (nveclen,a,b,c,d,Temp,1,Mice+1)
         Tsoil2(:) = Tsoil1(:,isoilcol)
       endif
       do i = 1, Mice+1
@@ -637,13 +632,15 @@ else
   endif
 endif
 
+deallocate(a, b, c, d, Temp)
+
 END SUBROUTINE T_DIFF
     
 
 !---------------------------------------------------------------------------------
 
 
-SUBROUTINE DIFF_COEF(a,b,c,d,n0,n1,m0,m1,substr,dt,isoilcol)
+SUBROUTINE DIFF_COEF(nn,a,b,c,d,n0,n1,m0,m1,substr,dt,isoilcol)
 
 !-------------------DEFINES COEFFICIENTS FOR SOLVING A THERMAL DIFFUSIVITY-----------------
 !-------------------EQUATION BY FACTORIZATION METHOD----------------------------------------
@@ -674,20 +671,8 @@ integer(kind=iintegers), intent(in) :: substr
 
 integer(kind=iintegers), intent(in), optional :: isoilcol
 
-real(kind=ireals), intent(out), dimension(1:vector_length) :: a,b,c,d
-
-!Common variables
-!integer(kind=iintegers) :: itop
-!real(kind=ireals), dimension(1:ms) :: Tsn, cs
-!real(kind=ireals), dimension(1:ms) :: lams, q
-!real(kind=ireals) :: dz
-!real(kind=ireals) :: AL,DLT,DVT,ALLL,DL,ALV,DV,Z,T,WL,WV,WI,dens
-
-!common /snow_char/ Tsn,cs
-!common /watericesnowarr/ lams,q
-!common /SOILSOL/ AL(ML),DLT(ML),DVT(ML),ALLL(ML),DL(ML), &
-!& ALV(ML),DV(ML),Z(ML),T(ML),WL(ML),WV(ML),WI(ML),dens(ms)
-!common /SOILDAT/ dz(ms),itop
+integer(kind=iintegers), intent(in) :: nn
+real(kind=ireals), intent(out), dimension(1:nn) :: a, b, c, d
 
 !External functions
 real(kind=ireals), external :: DZETA
diff --git a/source/model/turb_mod.f90 b/source/model/turb_mod.f90
index 500c4d6d63d17dc0302fde58e16f6aa7c851d1a7..1f9ccb7095dd4d6d6314d988c38f52a684da4b26 100644
--- a/source/model/turb_mod.f90
+++ b/source/model/turb_mod.f90
@@ -1,7 +1,7 @@
 MODULE TURB_LAKE_MOD
 
 use NUMERICS, only : PROGONKA, IND_STAB_FACT_DB
-use NUMERIC_PARAMS, only : vector_length, small_value
+use NUMERIC_PARAMS, only : nveclen, small_value
 use INOUT, only : CHECK_UNIT
 use LAKE_DATATYPES, only : ireals, iintegers
 use DRIVING_PARAMS, only : missing_value
@@ -157,10 +157,12 @@ real(kind=ireals) :: C_eps1(M), C_eps2(M), C_eps3(M)
 real(kind=ireals) :: CE(M), CEt(M)
 real(kind=ireals) :: lam_E(M+1), lam_eps(M+1)
 
-real(kind=ireals) :: a(vector_length)
-real(kind=ireals) :: b(vector_length)
-real(kind=ireals) :: c(vector_length)
-real(kind=ireals) :: d(vector_length)
+!real(kind=ireals) :: a(nveclen)
+!real(kind=ireals) :: b(nveclen)
+!real(kind=ireals) :: c(nveclen)
+!real(kind=ireals) :: d(nveclen)
+
+real(kind=ireals), allocatable :: a(:), b(:), c(:), d(:) 
 
 real(kind=ireals) :: AG(5)
 
@@ -320,6 +322,9 @@ firstcall_if : if (firstcall) then
 
 endif firstcall_if
 
+allocate(a(1:nveclen), b(1:nveclen), &
+&        c(1:nveclen), d(1:nveclen))
+
 if (l1 > 0._ireals) then
   !The minimal values for E and epsilon are determined from two conditions:
   !1.corresponding turbulent diffusion coefficient is much less than molecular coefficient
@@ -684,7 +689,7 @@ allocate (rhotemp(1:M), rhosal(1:M))
       & (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1)) 
     enddo
     ind_bound = .false.
-    call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
+    call IND_STAB_FACT_DB(nveclen,a,b,c,1,M,indstab,ind_bound)
     if (indstab .eqv. .false.) then
       do i = 2, M-1
         a(i) = - k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
@@ -692,7 +697,7 @@ allocate (rhotemp(1:M), rhosal(1:M))
       enddo
     endif
      
-    call PROGONKA (a,b,c,d,E_it2,1,M)
+    call PROGONKA (nveclen,a,b,c,d,E_it2,1,M)
     E_it2(M+1) = 0.
                 
     call CHECK_MIN_VALUE(E_it2, M, E_min) 
@@ -882,7 +887,7 @@ allocate (rhotemp(1:M), rhosal(1:M))
       & (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
     enddo
     ind_bound = .false.
-    call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
+    call IND_STAB_FACT_DB(nveclen,a,b,c,1,M,indstab,ind_bound)
     if (indstab .eqv. .false.) then
       do i = 2, M-1
         a(i) = - k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1))
@@ -890,7 +895,7 @@ allocate (rhotemp(1:M), rhosal(1:M))
       enddo
     endif
       
-    call PROGONKA (a,b,c,d,eps_it2,1,M)
+    call PROGONKA (nveclen,a,b,c,d,eps_it2,1,M)
     eps_it2(M+1) = 0.
  
     call CHECK_MIN_VALUE(eps_it2, M, eps_min)
@@ -1124,6 +1129,7 @@ allocate (rhotemp(1:M), rhosal(1:M))
 
   deallocate (work)
   deallocate (rhotemp, rhosal)
+  deallocate (a, b, c, d)
     
   do i = 1, M+1     
     E1(i) = E2(i)