diff --git a/data/BBS_unknow_error/BBS.dat b/data/BBS_unknow_error/BBS.dat
new file mode 100644
index 0000000000000000000000000000000000000000..e65c7f686ca7e8c5954d4acb0f340f45e39b1f5b
--- /dev/null
+++ b/data/BBS_unknow_error/BBS.dat
@@ -0,0 +1,2 @@
+Date ___________Temp. K Pressure. Pa___Humidity.Kg/kg SOx. m\s SOy. m|s Clfrac Prec.m/s Sol.rad W/m^2 Atm.Rad
+25.05.21,0:00,291.75,100140,.0111989814,1.6629831585,0.6888301783,0,0,0,350.3
diff --git a/data/BBS_unknow_error/BBS_driver.dat b/data/BBS_unknow_error/BBS_driver.dat
new file mode 100644
index 0000000000000000000000000000000000000000..641d387ed0fcfa1bf4c6d968e8b13f7324e44c2c
--- /dev/null
+++ b/data/BBS_unknow_error/BBS_driver.dat
@@ -0,0 +1,236 @@
+# The parameters for driver of Lake model
+#-----------------------------------------------------------------------------------------
+# INFORMATION ON THE FILE WITH ATMOSPHERIC DATA
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# dataname* --- name of file with atmospheric data (must be in data directory)
+# height_T_q* --- height of temperature and humidity measurements, m
+# height_u* --- height of wind measurements, m
+# interval* --- time interval of measurements, hours
+# rad* --- defines, atmospheric radiation (1) or net radiation (2) is
+# in the appropriate coloumn of datafile, relevant if input file is ASCII (forc_format = 0)
+# forc_format*--- defines the input file format: 0 - ASCII(text), 1 - netcdf
+# npoints* --- the number of points of the forcing, must be 1, if input file is ASCII (forc_format = 0)
+# select_call --- the length of the set of numbers of forcing points (maximal 20), for which the Lake model
+# will be launched, other points will be omitted
+# form* --- defines the input file format, relevant if input file is ASCII (forc_format = 0):
+# 0 - "free" (adjustable) format (see below)
+# other options are disabled
+#
+# The parameters for adjustable format of input text file, relevant if forc_format = 0, form = 0 :
+#
+# N_header_lines* --- the number of lines, occupied by file header
+# N_coloumns* --- the total number of coloumns in the file
+# N_Year* --- the number of coloumn with the number of year (not used in the model)
+# N_Month* --- the number of coloumn with the number of month (not used in the model)
+# N_Day* --- the number of coloumn with the number of day (not used in the model)
+# N_Hour* --- the number of coloumn with the number of hour (not used in the model)
+# N_Uspeed* --- the number of coloumn with x-component speed values, (m/s)
+# N_Vspeed* --- the number of coloumn with y-component speed values, (m/s)
+# N_Temp* --- the number of coloumn with air temperature values, (K)
+# N_Hum* --- the number of coloumn with air humidity values, (kg/kg)
+# N_Pres* --- the number of coloumn with atmospheric pressure value,s (Pa)
+# N_SWdown* --- the number of coloumn with net solar radiation values, (W/m**2)
+# N_LWdown* --- the number of coloumn with net longwave radiation values, (W/m**2)
+# N_Precip* --- the number of coloumn with precipitation intensity, (m/s)
+#-----------------------------------------------------------------------------------------
+#
+dataname 'BBS'
+forc_format 0
+npoints 1
+#select_call
+lakinterac 1
+form 0
+height_T_q 2.
+height_u 11.
+interval 3.
+rad 1
+#
+N_header_lines 1
+N_coloumns 11
+#
+N_Year -1
+N_Month -1
+N_Day -1
+N_Hour -1
+N_Uspeed 6
+N_Vspeed 7
+N_Temp 3
+N_Hum 5
+N_Pres 4
+N_SWdown 10
+N_LWdown 11
+N_NetRad -1
+N_Precip 9
+N_SensFlux -1
+N_LatentFlux -1
+N_Ustar -1
+N_surfrad -1
+N_cloud 8
+N_SurfTemp -1
+#
+#-----------------------------------------------------------------------------------------
+# TIME INTEGRATION PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# year0* --- julian year of start of integration
+# month0* --- julian month of start of integration
+# day0* --- julian day of start of integration
+# hour0* --- hour of start of integration (is real value)
+# dt* --- timestep, s
+# tinteg* --- integration time (including spinup period!), days
+# spinup_times* --- number for spinup periods
+# spinup_period--- the duration of spinup period, s
+# call_Flake --- the switch for integrating Flake model (1 - on, 0 - off)
+#-----------------------------------------------------------------------------------------
+#
+year0 2021
+month0 5
+day0 25
+hour0 0.
+#
+tinteg 98
+spinup_times 0
+spinup_period 0
+dt 3600
+call_Flake 0
+control_point 0
+cp_period 0.
+#
+#----------------------------------------------------------------------------------------
+# PHYSICAL PARAMETERS
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# extwat --- coefficient of solar radiation extinction in water body, m**(-1)
+# extice --- coefficient of solar radiation extinction in ice layer, m**(-1)
+# alphax --- slope angle of water surface in the x-direction, deg
+# alphay --- slope angle of water surface in the y-direction, deg
+# c_veg --- friction coefficient of vegetation in water, n/d
+# a_veg --- effective cross-section of vegetation, m**2/m**2
+# h_veg --- the height of vegetation in the lake, m
+# kor --- Coriolis parameter, s**(-1)
+# if -999, it is calculated from latitude
+# phi* --- latitude (positive to North), deg (required only for water albedo calculation)
+# lam* --- longitude (positive to East), dag (required only for water albedo calculation)
+# fetch --- the wind fetch, m
+#----------------------------------------------------------------------------------------
+#
+extwat 1.04
+extice 0.
+alphax 0.0
+alphay 0.0
+a_veg 1.
+c_veg 1.e-3
+h_veg 0.
+kor -999.
+phi 66.548333
+lam 33.134719
+fetch 1.E+3
+#
+#----------------------------------------------------------------------------------------
+# INITIAL CONDITIONS
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# l10 --- thickness of ice, m
+# h10 --- thickness of liquid water, m
+# select_h10 --- thickness of liquid water in selected points (maximal 20), m
+# (must be specified AFTER select_call)
+# ls10 --- thickness of ice at the bottom, m
+# hs10 --- thickness of snow cover, m
+# Ts0 --- temperature of mixed layer, C
+# Tb0 --- temperature at the bottom, C
+# Sals0 --- salinity in mixed layer, kg/kg
+# Salb0 --- salinity at the bottom, kg/kg
+# us0 --- x-component of speed at the surface, m/s
+# vs0 --- y-component of speed at the surface, m/s
+# Tbb0 --- temperature at the lower boundary of soil, C
+# h_ML0 --- thickness of mixed layer, m
+# init_T--- the type of temperature profile initialization:
+# 1 - using h_ML0, Ts0 and Tb0
+# 2 - using Tm, Ts0 and Tb0
+# 3 - using the temperature profile, specified in *_setup.dat file after keyword T_profile
+#-----------------------------------------------------------------------------------------
+#
+l10 0.
+h10 4.8
+#select_h10 3
+#
+ls10 0.
+hs10 0.
+Ts0 0.
+Tb0 0.
+Tbb0 0.
+Tm 0.
+h_ML0 0.
+Sals0 0.
+Salb0 0.
+us0 0.
+vs0 0.
+init_T 3
+#-----------------------------------------------------------------------------------------
+# SOME LAKE PARAMETERS
+#-----------------------------------------------------------------------------------------
+# area_lake --- the area of the lake, m**2
+# trib_inflow* --- total tributaries' inflow, m**3/s
+# morphometry --- depth - lake cross-section area table, m**2
+# effl_outflow --- effluent discharge parameters group, specifying polynomial dependence
+# of discharge on water level; the last value is relative altitude of
+# effluent bottom over lake bottom (at deepest points, respectively);
+# the first N values are polynomial coefficients, where N stands after
+# 'effl_outflow' keyword
+area_lake 6000.
+cellipt 2
+lakeform 1
+trib_inflow -9999.
+effl_outflow 1
+0. 0. 0.
+morphometry 12
+0. 6000.
+1. 5400.
+2. 4800.
+3. 4200.
+4. 3600.
+5. 3000.
+6. 2400.
+7. 1800.
+8. 1200.
+9. 600.
+10. 200.
+11. 50.
+#
+#-----------------------------------------------------------------------------------------
+# NETCDF OUTPUT PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# nstep_ncout --- the interval of netcdf output from driver, timesteps (if -1 no netcdf output from driver)
+#-----------------------------------------------------------------------------------------
+nstep_ncout -1
+#
+#-----------------------------------------------------------------------------------------
+# FLAKE MODEL OUTPUT PARAMETERS
+#-----------------------------------------------------------------------------------------
+# DESCRIPTION
+# nstep_out_Flake --- the interval of output of Flake variables from driver, timesteps
+# (if -1 the output of Flake variables from driver is not implemented)
+# relevant if call_Flake = 1
+#-----------------------------------------------------------------------------------------
+#
+nstep_out_Flake -1
+#
+#-----------------------------------------------------------------------------------------
+# POSTPROCESSING OPTIONS
+#-----------------------------------------------------------------------------------------
+# moving_average_window --- the ineterval of moving average, netcdf output steps (intervals)
+# mean_cycle_period --- the length of mean cycle, netcdf output steps (intervals)
+#-----------------------------------------------------------------------------------------
+#
+moving_average_window -1
+mean_cycle_period -1
+#
+end
+
+
+
+
+
+
diff --git a/data/BBS_unknow_error/BBS_setup.dat b/data/BBS_unknow_error/BBS_setup.dat
new file mode 100644
index 0000000000000000000000000000000000000000..3359b80a006965a0f9f7ec1fb32c4a4bc3f8c72d
--- /dev/null
+++ b/data/BBS_unknow_error/BBS_setup.dat
@@ -0,0 +1,318 @@
+#---------------------------------------------------------------------------------------
+# GENERAL CONTROLS
+#---------------------------------------------------------------------------------------
+# DESCRIPTION
+# path --- the directory, in which the model is launched
+# (required to be set in some UNIX-systems)
+# runmode 1 --- stand alone run
+# runmode 2 --- running in atmospheric model as lake parametrization
+# omp 0 --- OpenMP is not used
+# omp 1 --- OpenMP is used
+#---------------------------------------------------------------------------------------
+#
+path ''
+runmode 1
+omp 0
+#
+#----------------------------------------------------------------------------------------
+# SPATIAL RESOLUTION OF THE MODEL
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# nstep_keps --- number of timesteps of k-epsilon parameterization per on model timestep
+# M --- number of layers in water layer
+# Mice --- number of layers in upper and deep layers
+# ns --- number of levels in soil
+# d_surf --- grid zooming parameter at the surface, n/d
+# d_bot --- grid zooming parameter at the bottom, n/d
+#----------------------------------------------------------------------------------------
+#
+nstep_keps 1
+M 20
+ns 1
+nsoilcols 1
+Mice 0
+d_surf 1.E-2
+d_bot 1.E-2
+#
+#
+#----------------------------------------------------------------------------------------
+# CONTROLS FOR PHYSICS OF THE MODEL
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# PBL parameterization
+# PBLpar -1 --- sensible, latent heat and momentum fluxes are given as input for the model
+# PBLpar 0 --- the latent heat flux is set to zero, while sensible heat and momentum fluxes
+# are constant in time, specified by sensflux0 and momflux0
+# PBLpar 1 --- Businger-Dayer formulas (Monin-Obukhov theory) for exchange coefficients
+# PBLpar 2 --- formulation from NH3d
+# PBLpar 3 --- formulation from FLake
+# PBLpar 4 --- formulation implemented by M.Chechin
+# c_d --- the momentum exchange coefficient, n/d
+# (if -999, momentum flux is calculated by surface flux scheme)
+# waveenh 0 --- the shallow water correction of surface fluxes (Panin et al., 1996) is OFF
+# waveenh 1 --- the shallow water correction of surface fluxes (Panin et al., 1996) is ON
+# momflxpart 0 --- all momentum flux from the atmosphere is consumed by currents acceleration
+# momflxpart 1 --- momentum flux from the atmosphere is partitioned between wave developemnt
+# (controlled by fetch) and currents acceleration,
+# following Lin et al. (2002, J. Phys. Ocean.)
+# kwe --- the factor of turbulence enhancement by wave breaking (wave energy factor), n/d
+# Relative to water currents wind
+# relwind 0 --- relative wind is off
+# relwind 1 --- relative wind is on
+# Equation of state
+# eos 1 --- from Hostetler model
+# eos 2 --- from TEOS-2010
+# eos 3 --- for Kivu lake including salinity
+# nmeltpoint 1 --- melting point linearly dependent on salinity
+# nmeltpoint 2 --- TEOS-2010 formula
+# Turbulent mixing parameterization
+# Turbpar 1 --- "Empirical" parametrization: Stepanenko, Lykosov (2005)
+# Turbpar 2 --- "E-epsilon"("K-epsilon") parameterization: k=E**2/eps with
+# prognostic equations for E and eps
+# Turbpar 3 --- Nickuradze (NICK) formulation: Rodi (1993)
+# Turbpar 4 --- Parabolic (PARAB) formulation: Engelund (1976)
+# Turbpar 7 --- RNG (re-normalization group) formulation: Simoes (1998)
+# stabfunc 1 --- constant stability functions (standard k-epsilon model)
+# stabfunc 2 --- stability functions according to (Canuto et al., 2001)
+# stabfunc 3 --- stability functions according to (Galperin et al., 1988)
+# kepsbc 1 --- Neuman boundary conditions for unstratified sheared flow (Burchard, 2002)
+# kepsbc 2 --- Neuman boundary conditions for unstratified non-sheared flow with wave breaking (Burchard, 2002)
+# kepsbc 3 --- Neuman boundary conditions unstratified sheared flow with wave breaking (Burchard, 2002)
+# kepsbc 4 --- Neuman boundary conditions for free convection
+# Water surface albedo: variable or constant
+# varalb 0 --- constant
+# varalb 1 --- sun height dependent
+# soiltype 1 --- the soil type is "sand"
+# soiltype 2 --- the soil type is "loamy sand"
+# soiltype 3 --- the soil type is "sandy loam"
+# soiltype 4 --- the soil type is "loam"
+# soiltype 5 --- the soil type is "silt loam"
+# soiltype 6 --- the soil type is "sandy clay loam"
+# soiltype 7 --- the soil type is "clay loam"
+# soiltype 8 --- the soil type is "silty clay loam"
+# soiltype 9 --- the soil type is "sandy clay"
+# soiltype 10 --- the soil type is "silty clay"
+# soiltype 11 --- the soil type is "clay"
+# soil_depth --- depth of the soil layer, m
+# thermokarst_meth_prod 0. --- switch for old organics methane production under thermokarst lakes is OFF
+# thermokarst_meth_prod 1. --- switch for old organics methane production under thermokarst lakes is ON
+# soil_meth_prod 0. --- switch for new organics methane production under lakes is OFF
+# soil_meth_prod 1. --- switch for new organics methane production under lakes is ON
+# tricemethhydr 0. --- ice in soil pores is treated as pure ice
+# tricemethhydr 1. --- ice in soil pores is treated as methane hydrate
+# skin 0 --- the skin temperature parameterization is off
+# skin 1 --- the skin temperature parameterisation is on
+# sedim 0 --- gravitational sedimentation of tracer is NOT taken into account
+# sedim 1 --- gravitational sedimentation of tracer is taken into account
+# massflux 0 --- the massflux parameterization of convection (Siebesma et al., 2007) if OFF
+# massflux 1 --- the massflux parameterization of convection (Siebesma et al., 2007) if ON
+# sensflux0 --- sensible heat flux upwards, constant in time (relevant if PBLpar = 0), W/m**2
+# momflux0 --- momentum flux downwards (positive), constant in time (relevant if PBLpar = 0), N/m**2
+# ifrad 1 --- all radiation fluxes at the water surface are taken into account
+# ifrad 0 --- all radiation fluxes are set to zero
+# dyn_pgrad 0 --- dynamic pressure gradient is OFF
+# dyn_pgrad 1 --- dynamic pressure gradient is ON
+# zero_model 0 --- zero-dimensional model is ON
+# zero_model 1 --- zero-dimensional model is OFF
+# outflpar 0 --- variables value at the outflow = cross-section mean
+# outflpar 1 --- the cross-section mean = 0.5*(inflow value + outflow value)
+# outflpar 2 --- variables at the outflow are calculated using Lagrangian approach
+# Note: zero-dimensional model is now implemented only for open water season and one-point simulation
+# deadvol --- the depth (m) corresponding to "dead volume" - the minimal allowed reservoir volume
+#----------------------------------------------------------------------------------------
+#
+varalb 1
+PBLpar 3
+waveenh 0
+momflxpart 2
+c_d -999
+kwe 100.
+relwind 0
+eos 5
+lindens 0
+nmeltpoint 0
+Turbpar 2
+stabfunc 1
+kepsbc 1
+soiltype 5
+soil_depth 2.
+soilswitch 1
+saltice 0
+tricemethhydr 0.
+carbon_model 2
+skin 0
+massflux 0
+ifrad 1
+ifbubble 0
+sedim 0
+salsoil 0
+dyn_pgrad 0
+pgrad -1
+botfric 1
+horvisc 0
+backdiff 0
+backdiff0 -999.
+nManning -1
+zero_model 0
+thermokarst_meth_prod 0.
+soil_meth_prod 0.
+VmaxCH4aeroboxid -999.
+khsCH4 -999.
+khsO2 -999.
+r0methprod -999.
+outflpar 0
+#
+sensflux0 0.
+momflux0 0.
+soilbotflx 0.
+cuette 0
+#
+deadvol 0.
+#
+#----------------------------------------------------------------------------------------
+# INITIAL CONDITIONS FOR TEMPERATURE
+#----------------------------------------------------------------------------------------
+#
+T_profile 14
+0.1 10.4 9.3E-3 0. 0. 0. 0.
+0.5 10.0 9.4E-3 0. 0. 0. 0.
+1.0 10.5 13.7E-3 0. 0. 0. 0.
+1.5 9.1 22,9E-3 0. 0. 0. 0.
+2.0 6.8 23E-3 0. 0. 0. 0.
+2.1 6.4 23.1E-3 0. 0. 0. 0.
+2.2 6.2 22.9E-3 0. 0. 0. 0.
+2.3 6.3 23.1E-3 0. 0. 0. 0.
+2.4 6.2 23E-3 0. 0. 0. 0.
+2.5 6.4 22.9E-3 0. 0. 0. 0.
+3.0 5.0 23.4E-3 0. 0. 0. 0.
+3.5 4.9 23.6E-3 0. 0. 0. 0.
+4.0 4.4 23.7E-3 0. 0. 0. 0.
+4.5 4.2 23.9E-3 0. 0. 0.0 0.
+T_soilprofile 1
+0. 4.2
+#----------------------------------------------------------------------------------------
+# DATA ASSIMILATION CONTROLS (NOT OPERATIONAL: PUT EVERYTHING TO 0)
+#----------------------------------------------------------------------------------------
+# assim --- data assimilation technique: 0 - no data assimilation
+# 1 -
+# 2 - Raleigh damping towards observations
+# 3 - Cressman weighting
+# 4 - ?
+# as_window --- assimilation window: 1 - as_window is spinup period
+#
+#----------------------------------------------------------------------------------------
+#
+error_cov 0
+assim 0
+#as_window 1
+#
+#----------------------------------------------------------------------------------------
+# BOUNDARY CONDITIONS: TRIBUTARIES AND EFFLUENTS
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# tribheat --- the switch for thermal effect of tributaries and effluents, 0 - OFF
+# 1 - ON
+#-----------------------------------------------------------------------------------------
+#
+#
+tribheat 0
+N_tribin -1
+N_triblev -1
+fileinflow -1
+fileoutflow -1
+iefflloc -1
+dttribupdate -1
+#
+#----------------------------------------------------------------------------------------
+# OUTPUT CONTROLS (for ASCII files)
+#----------------------------------------------------------------------------------------
+# DESCRIPTION
+# turb_out --- output turbulence characteristics: 1 - on, 0 - off
+# monthly* --- monthly mean profiles output: 1 - on, 0 - off
+# daily* --- daily mean profiles output: 1 - on, 0 - off
+# hourly* --- hourly mean profiles output: 1 - on, 0 - off
+# everystep* --- every time step profiles output: 1 - on, 2 - on but without profiles, 0 - off
+# time_series* --- output of time series of layer thickness and surface values:
+# 1 - on, 0 - off
+# dt_out* --- time interval for time series output, hours
+# nscreen --- the period of screen output, timesteps
+# scale_output --- the switch for scaling of output of turbulent characteristics: 1 - on, 0 - off
+# (set 0 for simulations of ice-covered lakes)
+# ngrid_out --- the number of output levels for vertical water temperature profiles:
+# -1 - use numerical grid levels
+# >0 - use ngrid_out levels (in meters) given below
+# ngridsoil_out--- the number of output levels for vertical soil temperature profiles:
+# -1 - use numerical grid levels
+# >0 - use ngridsoil_out levels (in meters) given below
+#----------------------------------------------------------------------------------------
+#
+turb_out 1
+monthly 1
+daily 1
+hourly 1
+everystep 0
+time_series 1
+dt_out 1
+nscreen 1000
+zserout -999.
+scale_output 0
+accum_begin 2021052500
+accum_end 2021083100
+rtemp 1
+-999. -999. -999.
+#
+ngrid_out 43
+0.
+0.1
+0.2
+0.3
+0.4
+0.5
+0.6
+0.7
+0.8
+0.9
+1.0
+1.1
+1.2
+1.3
+1.4
+1.5
+1.6
+1.7
+1.8
+1.9
+2.0
+2.1
+2.2
+2.3
+2.4
+2.5
+2.6
+2.7
+2.8
+2.9
+3.0
+3.1
+3.2
+3.3
+3.4
+3.5
+3.6
+3.7
+3.8
+3.9
+4.0
+4.5
+4.8
+ngridice_out 0
+ngridsoil_out 0
+#
+#
+#----------------------------------------------------------------------------------------
+# NOTE: VARIABLES, DENOTED BY ASTERISK *, ARE USED ONLY IN STANDALONE RUNS OF THE MODEL
+#----------------------------------------------------------------------------------------
+end
+
+
+