MODULE LAKE_MOD

contains
SUBROUTINE INIT_LAKE(nx0,ny0,nx,ny,fnmap1,dt)

!The subroutine INIT_LAKE implements
!initialisation of the model (but not initial conditions)

use LAKE_DATATYPES, only : ireals, iintegers
use PHYS_PARAMETERS
use NUMERIC_PARAMS
use PHYS_CONSTANTS
use DRIVING_PARAMS 
use ATMOS
use ARRAYS
use TURB_CONST, only : &
& TURB_CONST_DERIVED_SET
use INOUT_PARAMETERS, only : &
& lake_subr_unit_min, &
& lake_subr_unit_max
use WATER_DENSITY
use SOIL_MOD, only : COMSOILFORLAKE
use SURF_SCHEME1_MOD, only : PBLDAT
use SURF_SCHEME_INM, only : PBLDAT_INM
use INOUT, only : fext2, readgrd_lake
use BL_MOD_LAKE, only : BL_MOD_LAKE_ALLOC
use SNOWSOIL_MOD, only : SNOWSOIL_MOD_ALLOC

implicit none

!Input variables
!Reals
real(kind=ireals), intent(in) :: dt

!Integers
integer(kind=iintegers), intent(in) :: nx0, ny0, nx, ny

!Characters
character(len=*), intent(in) :: fnmap1

!Output variables

!Local variables
!Reals
real(kind=ireals), parameter :: hour_sec = 60.*60.
real(kind=ireals) :: x1

!Integers
integer(kind=iintegers) :: n_unit = lake_subr_unit_min
integer(kind=iintegers) :: i, j, dnx, dny

!Logicals
logical :: uniform_depth

!Characters
character(len=80) :: path2
character(len=80) :: fnmap

data uniform_depth /.true./
!data n_unit /999/

dnx = nx - nx0 + 1
dny = ny - ny0 + 1

fnmap = fnmap1

call DEFINE_PHYS_PARAMETERS
call DEFINE_PHYS_CONSTANTS
call DEFINE_DRIVING_PARAMS
call ALLOCARRAYS(dnx,dny)
call BL_MOD_LAKE_ALLOC
call SNOWSOIL_MOD_ALLOC

call PBLDAT
call PBLDAT_INM
call COMSOILFORLAKE
call TURB_CONST_DERIVED_SET
call SET_DENS_CONSTS(eos%par)

if (error_cov==1) then
  if (runmode%par == 2) then
    print*, 'The error covariance estimation algorithm works &
    &only in standalone runs: STOP'
    STOP
  endif
  if (assim == 1) then
    print*, 'assim = 1 and error_cov = 1: these regimes could not &
    &be turned on simultaneously: STOP'
    STOP
  endif
  if (dnx> 1 .or. dny > 1) then
    print*, 'Error covariance calculation is currently adjusted &
    &only for one-point stand-alone runs of the lake model: STOP'
    STOP
  endif
endif
if (assim==1.and.(dnx>1.or.dny>1)) then
  print*, 'Data assimilation is currently adjusted &
  &only for one-point stand-alone lake model: STOP'
  STOP
endif

if (dt_out%par<dt/hour_sec) then
  print*, 'dt_out must be larger or equal to timestep: STOP'
  STOP
endif
   
if (runmode%par == 2.and.(.not.uniform_depth)) then
  path2 = fext2(fnmap,'dep')
  call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,n_unit)
  open (n_unit, file=path(1:len_trim(path))// &
  & path2(1:len_trim(path2)), status='old')
  call READGRD_LAKE(n_unit,dep_2d,0,nx-1,0,ny-1)
  close (n_unit)
  dep_av = 0.
  x1 = 0.
  do i = 1, nx-1
    do j = 1, ny-1
      if (dep_2d(i,j) > 0.) then
        dep_av = dep_av + dep_2d(i,j)
        x1 = x1 + 1.
      endif
    enddo
  enddo
  dep_av = dep_av / x1
endif

!print*, 'Lake model is initialized'

END SUBROUTINE INIT_LAKE

!>MAIN SUBROUTINE OF THE MODEL 
!!Calculates all parameters characterizing the state of a lake (temperature, currents,
!!eddy diffusivity coefficients, snow, ice, sediments characteristics, biogeochemistry etc.) 
!!at the next timestep (i.e. implements evolution of variables 
!!from i-th time step to (i+1)-th )
!!In the atmospheric model must be called once each timestep,
!!or once each N timesteps, where N = dt_lake/dt_atmos 
!!(dt_lake - timestep in the lake model, dt_atmos - timestep the atmospheric model) 
SUBROUTINE LAKE &
& (tempair1,humair1,pressure1,uwind1,vwind1, &
& longwave1,shortwave1,precip1,Sflux01, &
& ch4_pres_atm, co2_pres_atm, o2_pres_atm, zref1, dtl, &
& h10, l10, ls10, hs10, Ts0, Tb0, Tbb0, h_ML0, extwat, extice, &
& kor_in, trib_inflow, Sals0, Salb0, fetch, phi, lam, us0, vs0, &
& Tm, alphax, alphay, a_veg, c_veg, h_veg, area_lake, cellipt, depth_area, &
& tsw,hw1,xlew1,cdmw1,surfrad1,cloud1,ftot,h2_out, &
& ix,iy,nx0,ny0,nx,ny,nx_max,ny_max,ndatamax,year,month,day,hour,init_T,flag_assim,flag_print, &
& outpath1,spinup_done,dataread,lakeform,comm3d,rank_comm3d,coords,parallel,icp,step_final)



!OUTPUT VARIABLES:   
!tsw--- surface temperature of a lake, K;
!hw1--- sensible heat flux from lake, W/m**2; 
!xlew1   --- latent heat flux from lake, W/m**2; 
!cdmw    --- exchange coefficient, m/s ( == wind_speed*exchange_nondimensional_coeficient)

use LAKE_DATATYPES, only : ireals, iintegers
use OMP_LIB

use COMPARAMS, only: &
& num_ompthr, parparams

use NUMERIC_PARAMS
use DRIVING_PARAMS
use ATMOS
use PHYS_CONSTANTS
use ARRAYS
use ARRAYS_GRID
use ARRAYS_BATHYM
use ARRAYS_SOIL
use ARRAYS_TURB
use ARRAYS_WATERSTATE
use ARRAYS_METHANE
use ARRAYS_OXYGEN
use MODI_MASSFLUX_CONVECTION
use CONVECTPAR
use SALINITY, only : S_DIFF
use TURB, only : K_EPSILON, ED_TEMP_HOSTETLER, ED_TEMP_HOSTETLER2
use TURB_CONST, only : min_diff
use DZETA_MOD, only : VARMEAN

use MOMENTUM, only : &
& MOMENTUM_SOLVER

use METHANE_MOD, only : &
& METHANE, METHANE_OXIDATION

use SOIL_MOD, only : &
& SOILFORLAKE, &
& SOIL_COND_HEAT_COEF, &
& SOILCOLSTEMP, &
& LATERHEAT

use INIT_VAR_MOD, only : &
& INIT_VAR

use T_SOLVER_MOD, only : &
& T_SOLVER

use WATER_DENSITY, only: &
& DENSITY_W

use PHYS_FUNC, only: &
& TURB_SCALES, &
& MELTPNT, &
& WATER_FREEZE_MELT, &
& TURB_DENS_FLUX, &
& SINH0, &
& WATER_ALBEDO, &
& EXTINCT_SNOW, &
& UNFRWAT, &
& WI_MAX, &
& WL_MAX, &
& MIXED_LAYER_CALC, &
& HC_CORR_CARBEQUIL

use EVOLUTION_VARIABLES, only : &
& UPDATE_CURRENT_TIMESTEP, &
& UPDATE_NEXT_TIMESTEP

use TURB_CONST, only : &
& lamTM, min_diff

use METH_OXYG_CONSTANTS, only : &
& r0_oliglake, r0_oldorg2_star, &
& ngasb, mf, &
& molm3tonM, &
& molm3toppm_co2,molm3toppm_o2,molm3toppm_ch4, &
& molm3tomgl_co2,molm3tomkgl_ch4,molm3tomgl_o2, &
& molm3tomcM, &
& photic_threshold, &
& pH

use OXYGEN_MOD, only : &
& OXYGEN, OXYGEN_PRODCONS, &
& ADDOXPRODCONS, &
& CHLOROPHYLLA

use CARBON_DIOXIDE_MOD, only : &
& CARBON_DIOXIDE

use TIMEVAR, only : &
& hour_sec, &
& day_sec, &
& year_sec, &
& month_sec

use OUT_MOD

use BUBBLE_MOD, only : BUBBLE, BUBBLEFLUXAVER

use CONTROL_POINT, only : CONTROL_POINT_OUT, CONTROL_POINT_IN

use TRIBUTARIES, only : TRIBTEMP

use BATHYMSUBR, only : BATHYMETRY

use NUMERICS, only : ACCUMM

use BL_MOD_LAKE

use SNOWSOIL_MOD

use RADIATION, only : &
& RadWater, RadIce, RadDeepIce

implicit none



!Input variables
!Reals
!> Air temperature, K
real(kind=ireals), intent(in) :: tempair1
!> Specific humidity, kg/kg
real(kind=ireals), intent(in) :: humair1
!> Air pressure, Pa
real(kind=ireals), intent(in) :: pressure1
!> x- and y-components of wind, m/s
real(kind=ireals), intent(in) :: uwind1, vwind1
!> Longwave downward radiation, W/m**2
real(kind=ireals), intent(in) :: longwave1
!> Height of level above the lake surface, where the forcing is given, m
real(kind=ireals), intent(in) :: zref1
!> Global shortwave radiation, W/m**2
real(kind=ireals), intent(in) :: shortwave1
!> Precipitation intensity, m/s
real(kind=ireals), intent(in) :: precip1
real(kind=ireals), intent(in) :: Sflux01
!> Cloudiness, fraction
real(kind=ireals), intent(in) :: cloud1

real(kind=ireals), intent(in) :: ch4_pres_atm, co2_pres_atm, o2_pres_atm
!> Time step for the LAKE model, s
real(kind=ireals), intent(in) :: dtl

real(kind=ireals), intent(in) :: hour

real(kind=ireals), intent(in) :: h10, l10, ls10, hs10
real(kind=ireals), intent(in) :: Ts0, Tb0, Tbb0
real(kind=ireals), intent(in) :: h_ML0
real(kind=ireals), intent(in) :: extwat, extice
real(kind=ireals), intent(in) :: kor_in

!> Tributary inflow discharge, m**3/s
real(kind=ireals), intent(in) :: trib_inflow
real(kind=ireals), intent(in) :: Sals0, Salb0
real(kind=ireals), intent(in) :: fetch
real(kind=ireals), intent(in) :: phi, lam
real(kind=ireals), intent(in) :: us0, vs0
real(kind=ireals), intent(in) :: Tm
real(kind=ireals), intent(in) :: alphax, alphay
real(kind=ireals), intent(in) :: a_veg, c_veg, h_veg
real(kind=ireals), intent(in) :: area_lake, cellipt
real(kind=ireals), intent(in) :: depth_area(1:ndatamax,1:2)

!Integers

!> Indicator of the lake cross-section form
!! 1 - ellipse
!! 2 - rectangle
integer(kind=iintegers), intent(in) :: lakeform

!> Coordinates of the current lake on a horizontal mesh
integer(kind=iintegers), intent(in) :: ix, iy

!> Dimensions of the horizontal mesh
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: nx0, ny0, nx_max, ny_max
integer(kind=iintegers), intent(in) :: ndatamax
integer(kind=iintegers), intent(in) :: year, month, day
integer(kind=iintegers), intent(in) :: init_T
!> MPI parameters
integer(kind=iintegers), intent(in) :: comm3d, rank_comm3d, coords(1:3) 
!> Control point switch: 0 - do nothing
!!                       1 - write CP                             
!!                       2 - read CP instead of initial conditions                                           
integer(kind=iintegers), intent(in) :: icp 

                                           
!Characters
character(len=60), intent(in) :: outpath1

!Logicals
logical, intent(in) :: flag_assim
logical, intent(in) :: flag_print
logical, intent(in) :: spinup_done 
logical, intent(in) :: dataread
logical, intent(in) :: parallel
!> Indicator of the last time step
logical, intent(in) :: step_final

!Output variables
!Reals
real(kind=ireals), intent(inout) :: hw1, xlew1, cdmw1, surfrad1
real(kind=ireals), intent(out) :: tsw
real(kind=ireals), intent(out) :: ftot
real(kind=ireals), intent(out) :: h2_out

!------------------------ MAIN VARIABLES ------------------------
!  arrays:    Tw1 and Tw2 - temperature of water, C
! Ti1 and Ti2 - temperature of ice, C
! T - temperature of snow, C
!  functions of time:h1 and h2 - thickness of water, m
! l1 and l2 - thickness of ice, m
! hs1 - thickness of snow, m  
! flag_ice_surf - shows if ice is melting on the top surface, n/d  
!  constants: cw - specific heat of water, J/(kg*K) 
! ci - specific heat of ice, J/(kg*K)
! row0 - density of water, kg/m**3
! roi - density of ice, kg/m**3
! lamw - thermal conductivity of water, J*sec/(m**3*K)
! lami - thermal conductivity of ice, J*sec/(m**3*K)
! L - specific heat of melting, J/kg
! Lwv - specific heat of evaporation, J/kg
! Liv - specific heat of sublimation, J/kg    
! ddz - size of grid element (space), m
! dt - size of grid element (time), sec
! kstratwater - coefficient for linear initial conditions in water
! kstratice - coefficient for linear initial conditions in ice 
!  boundary conditions:
! eFlux - total heat flux on the upper surface,
!   shortwave(1-A)-T**4+longwave-H-LE, J/(m**2*sec)
! Elatent - latent heat flux, J/(m**2*sec) 
! Erad - shortwave radiation, penetrated below a surface, J/(m**2*sec)
! tempair - air temperature (2 m), C
! precip - precipitation, m/sec
!(M) number of layers in water and snow   

!Local variables
!Reals
real(kind=ireals), allocatable :: work1(:), work2(:), work3(:), work4(:), work5(:), work6(:)
real(kind=ireals) :: l2, h2, ls2
real(kind=ireals) :: dhwhigh, dhwlow, dhwp, dhwe, dhwtrib
real(kind=ireals) :: dhihigh, dhilow, dhip, dhis, dhif
real(kind=ireals) :: totalevap, totalprecip, totalmelt, totalpen, totalwat
real(kind=ireals) :: snowmass, snowmass_init
real(kind=ireals) :: totalerad, totalhflux
real(kind=ireals) :: h_ML0zv
real(kind=ireals) :: dhwls
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) :: dt
real(kind=ireals) :: b0
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) :: flux1, flux2

real(kind=ireals) :: dTisnmelt, dTiimp, dTmax, fracsn, fracimp

real(kind=ireals) :: calibr_par_min(1:2), calibr_par_max(1:2)
!data calibr_par_min /1.7d-8/5., 1.d-10/5./
!data calibr_par_max /1.7d-8*5., 1.d-10*5./
real(kind=ireals) :: step_calibr(1:2)


!Integers    
integer(kind=iintegers) :: i, j, k, nn = 0, i2, j2
integer(kind=iintegers) :: i_ML
integer(kind=iintegers) :: flag_snow
integer(kind=iintegers) :: nstep_meas
integer(kind=iintegers) :: n_1cm, n_5cm
integer(kind=iintegers) :: layer_case
integer(kind=iintegers) :: ndec
integer(kind=iintegers) :: npoint_calibr(1:2) = 12
integer(kind=iintegers) :: dnx, dny

!Logicals
logical :: uniform_depth
logical :: flag
logical :: accum = .false.
logical :: add_to_winter
logical :: firstcall = .true.
logical :: multsoil
logical :: worklog, worklog1, indTMprev = .false.

!Characters
character(len=60) :: workchar

data uniform_depth /.true./
data nstep_meas /0/
 
SAVE
  
!BODY OF PROGRAM

! Getting the number of cores
num_ompthr = OMP_GET_NUM_PROCS()
! MPI parameters
parparams%comm3d = comm3d
parparams%rank_comm3d = rank_comm3d
parparams%coords = coords
parparams%parallel = parallel

!call TIMEC(0)

gs%ix = ix
gs%iy = iy

dnx = nx - nx0 + 1
dny = ny - ny0 + 1


if (firstcall)  then
  if (calibrate_parameter) then
    ! Calibration runs are performed only when ny = 1
    if (npoint_calibr(1)*npoint_calibr(2) /= dnx) then
      write(*,*) 'npoint_calibr(1)*npoint_calibr(2) /= dnx in LAKE: STOP'
      STOP
    else
      ! Calibration of constants r0_oliglake and r0_oldorg2_star
      calibr_par1 => r0_oliglake
      calibr_par2 => r0_oldorg2_star
      cost_function => febultot
      
!      calibr_par_min(1) = calibr_par1/2. ; calibr_par_max(1) = calibr_par1*2.
      calibr_par_min(1) = 2.25d-8 ; calibr_par_max(1) = 3.d-8
!      calibr_par_min(2) = calibr_par2/sqrt(2.) ; calibr_par_max(2) = calibr_par2*sqrt(2.)
      calibr_par_min(2) = 6.d-11 ; calibr_par_max(2) = 8.d-11
      step_calibr(1:2) = (calibr_par_max(1:2) - calibr_par_min(1:2))/ &
      & real(npoint_calibr(1:2)-1)
    endif
  endif
endif

if (calibrate_parameter) then
  if (mod(ix,npoint_calibr(1)) == 0) then
    i = npoint_calibr(1)
    j = ix/npoint_calibr(1)
  else
    i = mod(ix,npoint_calibr(1))
    j = 1 + ix/npoint_calibr(1)
  endif
  calibr_par1 = calibr_par_min(1) + (i-1)*step_calibr(1)
  calibr_par2 = calibr_par_min(2) + (j-1)*step_calibr(2)
endif

!Checking if forcing data is read
if ( (.not.dataread) .and. PBLpar%par /= 0) then
  write(*,*) 'Atmospheric forcing is not available while PBLpar /= 0: STOP'
  STOP
endif

tempair = tempair1 - Kelvin0
humair = humair1
pressure = pressure1
!pressure = 1.d+5
uwind = uwind1
vwind = vwind1
zref  = zref1
precip = precip1
Sflux0 = Sflux01*(roa0/row0)
outpath = outpath1
cloud = cloud1

hw_input = hw1
xlew_input = xlew1
cdmw_input = cdmw1
surfrad_input = surfrad1

if (kor_in == -999.d0) then
  kor = 2.d0*omega*sin(phi*pi/180.d0)
else
  kor = kor_in
endif

!extwatarr(:) = extwat
RadWater  %extinct(:,:) = extwat
RadIce    %extinct(:,:) = extice
RadDeepIce%extinct(:,:) = extice
if (ifrad%par == 1) then
  longwave  = longwave1
  shortwave = shortwave1
elseif (ifrad%par == 0) then
  longwave  = 0.e0_ireals
  shortwave = 0.e0_ireals
endif

if (longwave1 == missing_value .and. firstcall) then
  write(*,*) 'Atmospheric radiation is missing in input file and will be calculated &
  &from net longwave radiation'
  if (cloud1 == missing_value) then
    STOP 'Cloudiness data is required in forcing when atmospheric radiation is missing!: STOP' 
  endif
endif

if (uwind == 0) uwind = 0.1
if (vwind == 0) vwind = 0.1
wind = sqrt(uwind**2+vwind**2)
!wind10 = wind*log(10./roughness0)/log(zref/roughness0)

!Control of the input atmospheric forcing
flag = .false.
if (tempair>60.or.tempair<-90) then
  print*, 'The air temperature ', tempair, 'is illegal: STOP'
  flag = .true.
elseif (abs(humair)>1.) then
  print*, 'The air humidity ', humair, 'is illegal: STOP'
  flag = .true.
elseif (pressure > 110000 .or. pressure < 80000.) then
  print*, 'The air pressure ', pressure, 'is illegal: STOP'
  flag = .true.
elseif (abs(uwind)>200.) then
  print*, 'The x-component of wind ', uwind, 'is illegal: STOP'
  flag = .true.
elseif (abs(vwind)>200.) then
  print*, 'The y-component of wind ', vwind, 'is illegal: STOP'
  flag = .true.
elseif (abs(longwave)>1000.) then
  print*, 'The longwave radiation ',longwave,'is illegal: STOP'
  flag = .true.
elseif (abs(shortwave)>1400.) then
  print*, 'The shortwave radiation ', shortwave, 'is illegal: STOP'
  flag = .true.
elseif (abs(precip)>1.) then
  print*, 'The atmospheric precipitation ',precip, &
  & 'is illegal: STOP'
  flag = .true.
endif
if (flag) then
  write(*,*) 'Temperature = ', tempair, &
  & 'Humidity = ', humair, &
  & 'Pressure = ', pressure, &
  & 'Uwind = ', uwind, &
  & 'Vwind = ', vwind, &
  & 'Longwave = ', longwave, &
  & 'Shortwave = ', shortwave, &
  & 'Precip = ', precip
  STOP
endif

dt = dtl
dt05 = 0.5d0*dt
dt_keps = dt/real(nstep_keps%par)
dt_inv = 1.d0/dt
dt_inv05 = 0.5d0*dt_inv

! Logical, indicating multiple soil columns configuration
multsoil = (nsoilcols > 1 .and. depth_area(1,2) >= 0. .and. soilswitch%par == 1)

if (init(ix,iy) == 0_iintegers) then
  ! Specification of initial profiles   
  
  rosoil(:) = 1200. ! Bulk soil density for initialization
  call INIT_VAR &
  ( M, Mice, ns, ms, ml, nsoilcols, ndatamax, &
  & init_T, skin%par, zero_model%par, &
  & h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
  & Sals0, Salb0, us0, vs0, h_ML0, &
  & rosoil, rosdry, por, depth_area, ip, &
  
  & flag_snow, flag_snow_init, itop, nstep, itherm, &
  & h1, l1, hs1, ls1, &
  & hx1, hx2, hy1, hy2, &
  & hx1t, hx2t, hy1t, hy2t, &
  & hx1ml, hx2ml, hy1ml, hy2ml, &
  & veg, snmelt, snowmass, cdmw2, velfrict_prev, &
  & roughness, eflux0_kinem, Elatent, &
  & totalevap, totalmelt, totalprecip, totalwat, totalpen, &
  & time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
  
  & E1, eps1, zsoilcols(1,ix,iy), Tsoil1, wi1, wl1, Sals1, &
  & rootss, qsoil, TgrAnn, qwater(1,1), &
  & oxyg(1,1), DIC(1,1), oxygsoil, Sal1, &
  & u1, v1, Tw1, Tskin, T_0dim, Eseiches, &
  & Ti1, Tis1, z_full, ddz, &
  & dz, T, wl, dens, lamw, &
  & dzeta_int, zsoil, pressure )

  Sals2(1:ns,1:nsoilcols) = Sals1(1:ns,1:nsoilcols) !Salinity is currently not calculated in lateral columns 
      
  !if (runmode%par == 2 .and. (uniform_depth .eqv. .false.)) then
  !  if (ix <= nx-1 .and. iy <= ny-1) then
  !    h1 = dep_2d(ix,iy)
  !  else
  !    h1 = dep_av
  !  endif  
  !  if (h1 <= 0) then
  !    print*, 'Negative or zero initial depth at ix=',ix, &
  !    & 'iy=',iy,'h1=',h1,':terminating program'
  !    STOP 
  !  endif
  !endif

endif

if (init(ix,iy) == 0_iintegers .and. icp == 2_iintegers &
  & .and. ix == 1 .and. iy == 1) then
  ! Read control point for inital conditions of all lakes of the domain
  call CONTROL_POINT_IN(nx,nx0,nx_max,ny,ny0,ny_max,gs,parparams)
endif


if (init(ix,iy) /= 0_iintegers .or. &
  & (init(ix,iy) == 0_iintegers .and. icp == 2_iintegers )) then
  call UPDATE_CURRENT_TIMESTEP ( &
  & ix, iy, dnx, dny, M, Mice, ns, ms, ml, nsoilcols, &
  & l1, h1, hx1, hx2, hy1, hy2, ls1, hs1, &
  & hx1t, hx2t, hy1t, hy2t, &
  & hx1ml, hx2ml, hy1ml, hy2ml, &
  & u1, v1, &
  & E1, eps1, k_turb_T_flux, &
  & Tsoil1, Sals1, wi1, wl1, &
  & Tw1, Sal1, lamw, &
  & Tskin(1), &
  & Ti1, Tis1, &
  & ueffl, veffl, Tweffl, Saleffl, &
  & dz, T, wl, dens, &
  & qwater(1,1), qsoil, &
  & oxyg(1,1), oxygsoil, &
  & DIC(1,1), DOC, POCL, POCD, &
  & snmelt, snowmass, &
  & cdmw2, &
  & time, &
  & dhwfsoil, &
  & Elatent, &
  & dhw, dhw0, &
  & dhi, dhi0, dls0, &
  & velfrict_prev, &
  & roughness, &
  & eflux0_kinem, &
  & tot_ice_meth_bubbles, &
  & febul0, Eseiches, &
  
  & flag_snow, flag_snow_init, &
  & itop, &
  & nstep, i_maxN, itherm )
endif

! Updating bathymetry 
call BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
               & depth_area,area_lake,cellipt, &
               & multsoil,trib_inflow,dhwtrib,vol,botar)
! Tributary inflow
if (trib_inflow /= -9999.) then
  dhwtrib = trib_inflow*dt/area_int(1)
  ! Will be corrected for water abstraction in BATHYMETRY 
else
  dhwtrib = (h10 - h1)/10. ! 10 - arbitrary value, the "time" of depth damping towards initial depth
endif

time = time + dt
nstep = nstep + 1

layer_case = 1
if (h1 >  0  .and. l1 >  0) layer_case = 2 
if (h1 == 0  .and. l1 >  0) layer_case = 3
if (h1 == 0  .and. l1 == 0) layer_case = 4

if (layer_case == 1 .and. &
& WATER_FREEZE_MELT(Tw1(1), 0.5*ddz(1)*h1, &
& Meltpnt(Sal1(1),preswat(1),nmeltpoint%par), +1) .and. &
& h1 - min_ice_thick * roi/row0 > min_water_thick) then
  layer_case = 2
endif 

if (layer_case == 3 .and. &
& WATER_FREEZE_MELT(Ti1(Mice+1), 0.5*ddzi(Mice)*l1, Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par), -1) .and. &
& l1 - min_water_thick * row0_d_roi > min_ice_thick) then
  layer_case = 2
endif  

if (layer_case == 3 .and. &
& WATER_FREEZE_MELT(Ti1(1), 0.5*ddzi(1)*l1, Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par), -1) .and. &
& l1 - min_water_thick * row0_d_roi > min_ice_thick .and. flag_snow == 0) then
  h1 = min_water_thick ; dhw = 0.; dhw0 = 0.
  Tw1 = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par) + T_phase_threshold 
  Sal1 = 1.d-5
  ls1 = l1 - min_water_thick * row0_d_roi
  Tis1 = Ti1
  l2 = 0
  Ti2 = 0
  if (ls1 < 0.009999) print*, 'Thin layer! - ls'
  layer_case = 1 
endif

  
if1: IF (layer_case == 1) THEN 

!---------------------------- CASE 1: LIQUID WATER LAYER AT THE TOP ("Summer")-----------------------


!Check the bottom layer water temperature
x = Meltpnt(Sal1(M+1),preswat(M+1),nmeltpoint%par)
if (WATER_FREEZE_MELT(Tw1(M+1), 0.5*ddz(M)*h1, x, +1) .and. &
& ls1 == 0 .and. h1 - min_ice_thick * roi/row0 > min_water_thick) then
  ls1 = min_ice_thick ; dls = 0.; dls0 = 0.
  Tis1 = x - T_phase_threshold
  h1 = h1 - min_ice_thick*roi/row0
endif   

! Updating bathymetry 
call BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
                & depth_area,area_lake,cellipt, &
                & multsoil,trib_inflow,dhwtrib,vol,botar)

!Check the liquid temperature to be above melting point
do i = 2, M
  Tw1(i) = max(Tw1(i),Meltpnt(Sal1(i),preswat(i),nmeltpoint%par))
enddo

!Calculation of water density profile at the current timestep
call DENSITY_W(M,eos%par,lindens%par,Tw1,Sal1,preswat,row)
call MIXED_LAYER_CALC(row,ddz,ddz2,dzeta_int,dzeta_05int,h1,M,i_ML,H_mixed_layer,maxN)

!The salinity flux at the surface are passed to
!TURB_DENS_FLUX as zeros since they are currently      ______
!not taken into account in calculation of density flux w'row'	
 
turb_density_flux = TURB_DENS_FLUX(eflux0_kinem,0.e0_ireals,Tw1(1),0.e0_ireals)
Buoyancy0 = g/row0*turb_density_flux

! Updating tributaries data and adding inflow of all substances
if (tribheat%par > 0) call TRIBTEMP(time,dt,h1,z_full,area_int,gsp, &
& gas,wst,ueffl,veffl,Tweffl,Saleffl,spinup_done)

!if (massflux%par == 1) then
!!Updating density profile for massflux convection
!!Updating z-grid properties 
!  do i = 1, M
!    dz_full (i) = ddz(i)  *h1
!    z_half  (i) = dzeta_05int(i)*h1
!  enddo
!  call MASSFLUX_CONVECTION( &
!  & nstep,dt,z_half,z_full,dz_full, &
!  & turb_density_flux,eflux0_kinem, &
!! The latent heat flux and salinity flux at the surface are passed as zeros,
!! as far as currently they are not taken 
!! into account in massflux conection parameterization
!  & 0.e0_ireals,0.e0_ireals, &
!! The surface wind stress components are passed as zeros,
!! as far as currently they are not taken 
!! into account in massflux conection parameterization
!  & 0.e0_ireals,0.e0_ireals, &
!  & u1,v1,w1,E1(1:M), &
!  & u1,v1,w1,row,Tw1,Sal1, &
!  & PEMF,PDENS_DOWN,PT_DOWN,PSAL_DOWN,zinv,1,M)
!! Debugging purposes only
!! PEMF = 0.e0_ireals
!! call TEMP_MASSFLUX(PEMF,PT_DOWN,ddz,h1,dt,Tw2,T_massflux)
!  do i = 1, M
!    pt_down_f(i) = 0.5d0*(PT_DOWN(i) + PT_DOWN(i+1))
!  enddo
!endif

! The solar radiation, absorbed by the water surface
if (varalb%par == 0) then
  Erad = shortwave*(1-albedoofwater)
elseif (varalb%par == 1) then
  Erad = shortwave* &
  & (1 - WATER_ALBEDO( SINH0(year,month,day,hour,phi) ) ) ! Note: the time here must be a local one
endif  ! time, not UTC

!Radiation fluxes in layers
x = 1._ireals - sabspen*skin%par
call RadWater%RAD_UPDATE(ddz05(0:M)*h1,(/Erad*(1 - sabs0)*x/))
! Calculating photic zone depth (<1% of surface irradiance)
if (RadWater%integr(1) > 0) then
  H_photic_zone = h1
  do i = 1, M
    if (RadWater%integr(i)/RadWater%integr(1) < photic_threshold) then
      H_photic_zone = z_full(i)
      exit
    endif
  enddo
else
  H_photic_zone = 0.
endif
if (ls1 /= 0.e0_ireals) then
  ! xx - radiation penetrated through the water - deep ice interface
  xx = RadWater%flux(M+1,1)*(1-albedoofice)
  call RadDeepIce%RAD_UPDATE(ddzi05(0:Mice)*ls1, (/xx/))
endif

!Calculation of the layer's thickness increments
if (ls1 == 0) then
  ice = 0; snow = 0; water = 1; deepice = 0
  dhwp = precip*dt
  dhwe = - Elatent/(Lwv*row0)*dt
  dhw = dhwp + dhwe + dhwfsoil + dhwtrib
  dhw0 = dhwp + dhwe
  dls = 0.
else
  ice = 0; snow = 0; water = 1; deepice = 1
  flux1 = - lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + RadWater%integr(M) + &
  & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.*dt) + &
  & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
  flux2 = - lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + RadDeepIce%integr(1) + &
  & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.*dt)
  dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
!  dhwls = -lamw(M)*dt/(ddz(M)*h1*row0_m_Lwi)*(Tw1(M+1) - Tw1(M))+
!&  lami*dt/(ddz(1)*ls1*row0_m_Lwi)*(Tis1(2) - Tis1(1))+
!&  (1-albedoofice)*Erad*(1-sabs)*exp(-extwat*h1)/(row0_m_Lwi)*dt 
  dls = - dhwls*row0_d_roi
  dls0 = dls
  dhwp = precip*dt
  dhwe = - Elatent/(Lwv*row0)*dt
  dhw = dhwp + dhwe + dhwfsoil + dhwls + dhwtrib
  dhw0 = dhwp + dhwe
endif

!Calculation of water current's velocities and turbulent characteristics
if (Turbpar%par == 2) then
  call MOMENTUM_SOLVER(ix, iy, dnx, dny, ndatamax, &
  & year, month, day, hour, kor, a_veg, c_veg, h_veg, &
  & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr, tau_wav, fetch, depth_area)
endif

if (soilswitch%par == 1) then
  call SOIL_COND_HEAT_COEF(nsoilcols)
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, &
                   & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                   & bathymwater,bathymice, &
                   & bathymdice,bathymsoil(1,ix,iy), &
                   & soilflux(1,ix,iy),fdiffbot,(/1,0/))
  call LATERHEAT(ix,iy,gs,ls, &
  & bathymwater,bathymice,bathymdice,bathymsoil, &
  & gsp,soilflux(1,ix,iy),.false.,lsh)
endif
! Calculation of the whole temperature profile
call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
& RadWater%extinct(1,1), extice, fetch, dt)

! Diagnostic calculations
do i = 1, M
  k_turb_T_flux(i) = - 0.5*lamw(i)/(cw_m_row0) * &
  & ( Tw2(i+1) + Tw1(i+1) - Tw2(i) - Tw1(i)) / (ddz(i)*h1)
  T_massflux(i) = PEMF(i)*(pt_down_f(i) - 0.5d0*(Tw2(i+1) + Tw2(i)))
enddo

lamsal(:) = (lamw(:) - lamw0)/cw_m_row0 * alsal + lamsal0
lammeth(:) = lamw(:)/cw_m_row0 * almeth
lamoxyg(:) = lamw(:)/cw_m_row0 * aloxyg
lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

! Salinity diffusion in water and soil
CALL S_DIFF(dt,dhilow)

! Oxygen model precedes methane module, in order sedimentary 
! oxygen demand to be available for the latter
call CHLOROPHYLLA(z_full, H_mixed_layer, H_photic_zone, extwat, M, Chl_a, itroph)
call OXYGEN_PRODCONS(gs, gsp, wst, bathymsoil, gas, area_int, area_half, ddz05, &
& h1, dt, Tsoil3, Chl_a, dzs, por, oxygsoil, itroph, &
& prodox, resp, bod, sod, sodbot)

if (soilswitch%par == 1) then
  call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
  ! 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, &
                     & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                     & bathymwater,bathymice, &
                     & bathymdice,bathymsoil(1,ix,iy), &
                     & soilflux(1,ix,iy), fdiffbot, (/0,1/))
    methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
  endif
  call BUBBLE_BLOCK
  ! Calculation of methane sources from methane fluxes from mutliple soil columns
  if (multsoil) then 
    !call LATERHEAT(ix,iy,gs,ls, &
    !& bathymwater,bathymice,bathymdice,bathymsoil, &
    !& gsp,soilflux(1,ix,iy),.false.,lsh)
    call LATERHEAT(ix,iy,gs,ls, &
    & bathymwater,bathymice,bathymdice,bathymsoil, &
    & gsp,methsoilflux,.true.,lsm)
    call LATERHEAT(ix,iy,gs,ls, &
    & bathymwater,bathymice,bathymdice,bathymsoil, &
    & gsp,co2soilflux, .true.,lsc)
    call LATERHEAT(ix,iy,gs,ls, &
    & bathymwater,bathymice,bathymdice,bathymsoil, &
    & gsp,o2soilflux,  .true.,lso)
  endif
  gs%isoilcol = nsoilcols
  call METHANE &
  & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
  & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
  & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
  & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
  & lsm, bathymwater, &
  & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
  & plant_sum,bull_sum,oxid_sum,rprod_sum, &
  & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
  & ice_meth_oxid_total, &
  & h_talik,tot_ice_meth_bubbles,add_to_winter)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !call METHANE2 & ! two-meth
  !& (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !& ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !& fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !& plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !& anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !& h_talik,tot_ice_meth_bubbles2) ! two-meth
  !qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
else
  qwater(:,2) = qwater(:,1)
  qmethane(1:M+1) = qwater(1:M+1,2)
  qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols) 
endif
call ADDOXPRODCONS(M,i_maxN,prodox,resp,bod,sod,dt,ls,gas)
call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
& o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
call CARBON_DIOXIDE (gs, Tw2, fbbleflx_co2_sum, fbbleflx_co2(0,nsoilcols), lamcarbdi, gsp, bathymwater, lsc, fracb0, &
& pressure, wind10, co2_pres_atm, febul0(nsoilcols), ls, dt, eps1(1), sodbot, gas, soilswitch%par, fco2)
call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
& qwateroxidtot, qwateroxidML)
!call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth

!Calculation of water density profile at the next timestep
call DENSITY_W(M,eos%par,lindens%par,Tw2,Sal2,preswat,row2)

if (Turbpar%par == 2) then
  !$OMP PARALLEL NUM_THREADS(num_ompthr) IF(omp%par == 1) PRIVATE(i)
  do i = 1, nstep_keps%par
    call K_EPSILON(ix,iy,dnx,dny, year, month, day, hour, &
    & kor, a_veg, c_veg, h_veg, dt_keps, &
    & b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch) ! 0.5*tau_wav 0.5 - arbitrary (!) fraction of momentum loss due to wave breaking
  enddo
  !$OMP END PARALLEL  
else
! Calculation of the of eddy difusivity for temperature
  SELECT CASE (Turbpar%par)
    CASE (8)
      ! Hostetler formulation
      call ED_TEMP_HOSTETLER &
      & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
      KT(:) = cw_m_row0*KT(:)
    CASE (9)
      ! Modified Hostetler formulation for shallow lake
      call ED_TEMP_HOSTETLER2 &
      & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
      KT(:) = cw_m_row0*KT(:)
    CASE (1)
      do i = 1, M
        KT(i) = (k2(i)-niu_wat)*cw_m_row0 
      enddo
    CASE DEFAULT
      KT = lamTM*k2*cw_m_row0 
  ENDSELECT 
endif

! Heat conductance at the next timestep  
do i = 1, M
  lamw(i) = lamw0 + KT(i) !+KC(i) !+KLengT(i)
enddo

!The scales of turbulence
call TURB_SCALES(gsp,ls, bathymwater, wst, RadWater, &
& k_turb_T_flux, T_massflux, row, eflux0_kinem, &
& turb_density_flux, Buoyancy0, tau_air-tau_wav, kor, i_maxN, H_mixed_layer,maxN, w_conv_scale, &
& T_conv_scale, Wedderburn, LakeNumber, Rossby_rad, ThermThick, ReTherm, RiTherm, trb)

if (zero_model%par == 1) then
! Note that zero model is impemented currently only for
! 1. one-point simulation
! 2. open water conditions
  call ZERODIM_MODEL &
  & (h1, dt, &
  & shortwave, longwave, tempair, humair, pressure, &
  & uwind, vwind, zref, hw_input, xlew_input, cdmw_input, surfrad_input, &
  & cloud, botflux, T_0dim)
endif
  
h2 = h1 + dhw
ls2 = ls1 + dls
l2 = 0.

ENDIF if1   

if2: IF (layer_case == 2) THEN

!---------------------------CASE 2: WATER AND ICE ("Winter")------------------------------

!call TIMEC(1)

! Creation of the initial thin ice layer
if (l1 == 0.) then
  l1 = min_ice_thick; dhi = 0.; dhi0 = 0.; dhilow = 0.; dhihigh = 0.
  h1 = h1 - min_ice_thick * roi/row0
  x = Meltpnt(Sal1(1),0.e0_ireals,nmeltpoint%par)
  if (tempair < 0.) then    
    do i = 1, Mice+1
      Ti1(i) = tempair + float(i-1)/float(Mice) * (x - tempair)
    enddo 
  else 
    Ti1(:) = x - T_phase_threshold
  endif
end if 

!Check for water temperature to be higher than melting point
do i = 2, M
  Tw1(i) = max(Tw1(i),Meltpnt(Sal1(i),preswat(i),nmeltpoint%par))
enddo

!Creation of the initial thin water layer
if (h1 == 0.) then
  h1 = min_water_thick; dhw = 0.; dhw0 = 0.
  Sal1(:) = 1.d-5
  Tw1(:) = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par) + T_phase_threshold
  l1 = l1 - min_water_thick * row0_d_roi
end if

!Creation of the initial thin bottom ice layer
x = Meltpnt(Sal1(M+1),preswat(M+1),nmeltpoint%par)
if (WATER_FREEZE_MELT(Tw1(M+1), 0.5*ddz(M)*h1, x, +1) .and. &
& ls1 == 0 .and. h1 - min_ice_thick * roi/row0 > min_water_thick) then
  ls1 = min_ice_thick; dls = 0.; dls0 = 0.
  Tis1 = x - T_phase_threshold
  h1 = h1 - min_ice_thick * roi/row0
endif

! Updating bathymetry 
call BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
                & depth_area,area_lake,cellipt, &
                & multsoil,trib_inflow,dhwtrib,vol,botar)


!Calculation of water density profile at the current timestep
call DENSITY_W(M,eos%par,lindens%par,Tw1,Sal1,preswat,row)

! Note: the salinity effects and partial ice coverage 
! are currently neglected in density flux at the ice-water interface
eflux0_kinem_water = - lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1*cw_m_row0) ! The first-order approximation
turb_density_flux = TURB_DENS_FLUX(eflux0_kinem_water,0.e0_ireals,Tw1(1),0.e0_ireals)
Buoyancy0 = g/row0*turb_density_flux

!Currently the EDMF parameterization of convection is not implemented in the code
!during 'winter' (if the ice layer is present over
!water layer). This seems to be reasonable, since the temperature profile
!is stable in this case - until spring convection under thin ice
!develops.

PEMF = 0.e0_ireals
pt_down_f = 0.e0_ireals

! Updating tributaries data and adding inflow of all substances
if (tribheat%par > 0) call TRIBTEMP(time,dt,h1,z_full,area_int,gsp, &
& gas,wst,ueffl,veffl,Tweffl,Saleffl,spinup_done)
 
!CASE 2.1: WATER, ICE AND SNOW 

if (flag_snow == 1) then

  call MIXED_LAYER_CALC(row,ddz,ddz2,dzeta_int,dzeta_05int,h1,M,i_ML,H_mixed_layer,maxN)

! Radiation fluxes in physical layers
  SR_botsnow = shortwave*(1-albedoofsnow)
  if (flag_snow_init /= 1) then
    do i = itop, ms-1
      SR_botsnow = SR_botsnow*EXTINCT_SNOW(dens(i))**dz(i)
    enddo
  endif
  call RadIce  %RAD_UPDATE(ddzi05(0:Mice)*l1, (/SR_botsnow/))
  call RadWater%RAD_UPDATE(ddz05 (0:Mice)*h1, (/RadIce%flux(Mice+1,1)/))
! Calculating photic zone depth (<1% of surface irradiance)
  if (RadWater%integr(1) > 0) then
    H_photic_zone = h1
    do i = 1, M
      if (RadWater%integr(i)/RadWater%integr(1) < photic_threshold) then
        H_photic_zone = z_full(i)
        exit
      endif
    enddo
  else
    H_photic_zone = 0.
  endif
  if (ls1 /= 0.e0_ireals) then
    call RadDeepIce%RAD_UPDATE(ddzi05(0:Mice)*ls1, (/RadWater%flux(M+1,1)/))
  endif


! Calculation of the layer's thickness increments
  flux1 = - lami*(Ti1(Mice+1)-Ti1(Mice))/(ddzi(Mice)*l1) + RadIce%integr(Mice) + &
  & ci_m_roi*(dhi-dhi0)*(Ti1(Mice+1)-Ti1(Mice))/(2.d0*dt)
  flux2 = - lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1) + RadWater%integr(1) + &
  & cw_m_row0*dhw0*(Tw1(2)-Tw1(1))/(2.*dt) + &
  & cw_m_row0*PEMF(1)*(pt_down_f(1)-0.5d0*(Tw1(2)+Tw1(1)) )
  ! Limiting the freezing/melting rate by the thicknesses of water and ice
  dhwlow = min(max(dt*(flux1 - flux2)/(row0_m_Lwi),-ddz(1)*h1),ddzi(Mice)*l1*roi_d_row0)
  dhilow = - dhwlow*row0_d_roi
  x = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par)
  if (Ti1(1) > x + T_phase_threshold) then
    dhwhigh = (Ti1(1) - x - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi)
    dhihigh = - dhwhigh*row0_d_roi
    Ti1(1) = x + T_phase_threshold
  else
    dhwhigh = 0.e0_ireals
    dhihigh = 0.e0_ireals
  endif
  !Treatment of water freezing at the snow-ice interface
  !Freezing snowmelt
  call SNOW_COND_HEAT_COEF() !Needed for cs update
  y = 0.5*(cs(ms)*dens(ms)*dz(ms) + ci_m_roi*ddzi(1)*l1)! Heat capacity of snow-ice interface
  dTmax = max((1. - snmeltwat)*snmelt*dt*row0*Lwi/y, 0._ireals)
  dTisnmelt = min(dTmax, x + T_phase_threshold - Ti1(1)) ! x -- melting point!
  dTisnmelt = max(dTisnmelt, 0._ireals)
  fracsn = dTisnmelt/(dTmax + small_value) !Fraction of snowmelt, frozen at the ice top
  snmeltice = dTmax*y/(row0*Lwi)*fracsn
  dhihigh = dhihigh + snmeltice*row0_d_roi ! Adding snowmelt, frozen at the top of ice

  !Assessing impoundment of ice from above by water penetrated through cracks
  dTmax = max((snowmass + l1*(roi - row0))*Lwi/y, 0._ireals)
  dTiimp = min(dTmax, x + T_phase_threshold - Ti1(1) - dTisnmelt) ! x -- melting point!
  dTiimp = max(dTiimp, 0._ireals)
  fracimp = dTiimp/(dTmax + small_value)
  dhiimp = dTmax*y/(Lwi*row0)*fracimp
  dhihigh = dhihigh + dhiimp*row0_d_roi
  dhwhigh = dhwhigh - dhiimp

  if (ls1 /= 0) then
    ice = 1; snow = 1; water = 1; deepice = 1
    flux1 = - lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + RadWater%integr(M) + &
    & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.d0*dt) + &
    & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
    flux2 = - lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + RadDeepIce%integr(1) + &
    & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.d0*dt)
    dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
    dls = - dhwls*row0_d_roi
    dhw = dhwhigh + dhwlow  + dhwfsoil + dhwls + dhwtrib
    dhw = dhw + snmeltwat*snmelt*dt + (1. - snmeltwat)*snmelt*dt*(1. - fracsn)
    dhw0 = dhwlow + dhwhigh 
    dhw0 = dhw0 + snmeltwat*snmelt*dt + (1. - snmeltwat)*snmelt*dt*(1. - fracsn)
    dls0 = dls
  else
    ice = 1; snow = 1; water = 1; deepice = 0
    dls = 0.
    dhwls = 0.
    dhw = dhwhigh + dhwlow + dhwfsoil + dhwtrib
    dhw = dhw + snmeltwat*snmelt*dt + (1. - snmeltwat)*snmelt*dt*(1. - fracsn)
    dhw0 = dhwlow + dhwhigh
    dhw0 = dhw0 + snmeltwat*snmelt*dt + (1. - snmeltwat)*snmelt*dt*(1. - fracsn)
  endif
  dhi = dhilow + dhihigh 
  dhi0 = dhihigh

  if (soilswitch%par == 1) then
    call SOIL_COND_HEAT_COEF(nsoilcols)
  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, &
                     & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                     & bathymwater,bathymice, &
                     & bathymdice,bathymsoil(1,ix,iy), &
                     & soilflux(1,ix,iy), fdiffbot, (/1,0/))
    call LATERHEAT(ix,iy,gs,ls, &
    & bathymwater,bathymice,bathymdice,bathymsoil, &
    & gsp,soilflux(1,ix,iy),.false.,lsh)
  endif
  ! Calculation of the whole temperature profile 
  call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
  & RadWater%extinct(1,1), extice, fetch, dt)

  lamsal(:) = (lamw(:) - lamw0)/cw_m_row0 * alsal + lamsal0
  lammeth(:) = lamw(:)/cw_m_row0 * almeth
  lamoxyg(:)= lamw(:)/cw_m_row0 * aloxyg
  lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

  ! Salinity profile in water and soil
  call S_DIFF(dt,dhilow)
  if (soilswitch%par == 1) then
    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
    ! 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, &
                       & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                       & bathymwater,bathymice, &
                       & bathymdice,bathymsoil(1,ix,iy), &
                       & soilflux(1,ix,iy), fdiffbot, (/0,1/))
      methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
    endif
  endif
  call SNOWTEMP(ix,iy,dnx,dny,year,month,day,hour,snowmass, &
  & snowmass_init,a,b,c,d,Temp,phi,extice,fetch,dt)

  ! Oxygen model precedes methane module, in order sedimentary 
  ! oxygen demand to be available for the latter
  call CHLOROPHYLLA(z_full, H_mixed_layer, H_photic_zone, extwat, M, Chl_a, itroph)
  call OXYGEN_PRODCONS(gs, gsp, wst, bathymsoil, gas, area_int, area_half, ddz05, &
  & h1, dt, Tsoil3, Chl_a, dzs, por, oxygsoil, itroph, &
  & prodox, resp, bod, sod, sodbot)

  if (soilswitch%par == 1) then
    call BUBBLE_BLOCK
    ! Calculation of methane sources from methane fluxes from mutliple soil columns
    if (multsoil) then 
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,methsoilflux,.true.,lsm)
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,co2soilflux, .true.,lsc)
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,o2soilflux,  .true.,lso)
    endif 
    gs%isoilcol = nsoilcols
    call METHANE &
    & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
    & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
    & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
    & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false. ,&
    & lsm, bathymwater, &
    & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
    & ice_meth_oxid_total, &
    & h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !  call METHANE2 & ! two-meth
  !  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !  & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !  & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  else
    qwater(:,2) = qwater(:,1)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  endif
  call ADDOXPRODCONS(M,i_maxN,prodox,resp,bod,sod,dt,ls,gas)
  call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
  & o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
  call CARBON_DIOXIDE (gs, Tw2, fbbleflx_co2_sum, fbbleflx_co2(0,nsoilcols), lamcarbdi, gsp, bathymwater, lsc, fracb0, &
  & pressure, wind10, co2_pres_atm, febul0(nsoilcols), ls, dt, eps1(1), sodbot, gas, soilswitch%par, fco2)
  call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
  & qwateroxidtot, qwateroxidML)
!  call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth

!Calculation of water current's velocities and turbulent characteristics
  if (Turbpar%par /= 1) then
  
    if (Turbpar%par == 2) then
      call MOMENTUM_SOLVER(ix, iy, dnx, dny, ndatamax, &
      & year, month, day, hour, kor, a_veg, c_veg, h_veg, &
      & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr, tau_wav, fetch, depth_area)  
    endif

    !Calculation of water density profile at the current timestep
    call DENSITY_W(M,eos%par,lindens%par,Tw2,Sal2,preswat,row2)

    if (Turbpar%par == 2) then
      !$OMP PARALLEL NUM_THREADS(num_ompthr) IF(omp%par == 1)
      do i = 1, nstep_keps%par
        call K_EPSILON(ix,iy,dnx,dny, year, month, day, hour, &
        & kor, a_veg, c_veg, h_veg, dt_keps, &
        & b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch) ! 0.5*tau_wav 0.5 - arbitrary (!) fraction of momentum loss due to wave breaking
      enddo
      !$OMP END PARALLEL
    else
!     Calculation of the of eddy difusivity for temperature
      SELECT CASE (Turbpar%par)
        CASE (8)
          ! Hostetler formulation
          !call ED_TEMP_HOSTETLER &
          !& (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0.
        CASE (9)
          ! Modified Hostetler formulation for shallow lake
          ! call ED_TEMP_HOSTETLER2 &
          ! & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0. !cw_m_row0*KT(:)
        CASE (1)
          do i = 1, M
            KT(i) = (k2(i)-niu_wat)*cw_m_row0 
          enddo
        CASE DEFAULT
          KT = lamTM*k2*cw_m_row0 
      ENDSELECT 
    endif
    
    do i = 1, M
      lamw(i) = lamw0 + KT(i) !+KC(i)
    enddo  

  else
    lamw = lamw0*3.d0
  endif
  
  
  h2 = h1 + dhw
  l2 = l1 + dhi 
  ls2 = ls1 + dls
  
else

!CASE 2.2: WATER AND ICE WITHOUT SNOW

  call MIXED_LAYER_CALC(row,ddz,ddz2,dzeta_int,dzeta_05int,h1,M,i_ML,H_mixed_layer,maxN)

!Radiation fluxes in layers
  Erad = shortwave*(1-albedoofice)
  call RadIce  %RAD_UPDATE(ddzi05(0:Mice)*l1, (/Erad/))
  call RadWater%RAD_UPDATE(ddz05(0:M)*h1, (/RadIce%flux(Mice+1,1)/))
! Calculating photic zone depth (<1% of surface irradiance)
  if (RadWater%integr(1) > 0) then
    H_photic_zone = h1
    do i = 1, M
      if (RadWater%integr(i)/RadWater%integr(1) < photic_threshold) then
        H_photic_zone = z_full(i)
        exit
      endif
    enddo
  else
    H_photic_zone = 0.
  endif
  if (ls1 /= 0.e0_ireals) then
    call RadDeepIce%RAD_UPDATE(ddzi05(0:Mice)*ls1, (/RadWater%flux(M+1,1)/))
  endif

!Calculation of the layer's thickness increments
  x = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par)
  if (Ti1(1) > x + T_phase_threshold) then
    dhwhigh = (Ti1(1) - x - T_phase_threshold) * &
    & ci_m_roi*ddzi(1)*l1*0.5d0/(row0_m_Lwi)
    dhis = 0. 
    dhihigh = - dhwhigh*row0_d_roi
    Ti1(1) = x + T_phase_threshold
  endif
! Creation of the initial thin snow layer - to be considered at the next timestep
  if (precip > 0 .and. tempair < 0 .and. l1 > 0.05) then
    flag_snow = 1
    hs1 = 2.*dzmin !0.02
    itop = ms-2
    dhip = 0
  else
    dhip = precip*row0_d_roi*dt
  end if
  flux1 = - lami*(Ti1(Mice+1)-Ti1(Mice))/(ddzi(Mice)*l1)+RadIce%integr(Mice) + &
  &  ci_m_roi*(dhi-dhi0)*(Ti1(Mice+1)-Ti1(Mice))/(2.d0*dt)
  flux2 = -lamw(1)*(Tw1(2)-Tw1(1))/(ddz(1)*h1) + RadWater%integr(1) + &
  &  cw_m_row0*dhw0*(Tw1(2)-Tw1(1))/(2.*dt) + &
  &  cw_m_row0*PEMF(1)*(pt_down_f(1)-0.5d0*(Tw1(2)+Tw1(1)) )
  ! Limiting the freezing/melting rate by the thicknesses of water and ice
  dhwlow = min(max(dt*(flux1 - flux2)/(row0_m_Lwi),-ddz(1)*h1),ddzi(Mice)*l1*roi_d_row0)
  dhilow = - dhwlow*row0_d_roi

  if (Ti1(1) < Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par) + T_phase_threshold) then
    dhis = - Elatent/(roi*Liv)*dt
    dhwhigh = 0
    dhihigh = 0
  endif
  dhi = dhihigh + dhilow + dhis + dhip
  dhi0 = dhis + dhihigh + dhip
  if (ls1 /=0 ) then
    ice = 1; snow = 0; water = 1; deepice = 1 
    flux1 = - lamw(M)*(Tw1(M+1)-Tw1(M))/(ddz(M)*h1) + RadWater%integr(M) + &
    & cw_m_row0*(dhw-dhw0)*(Tw1(M+1)-Tw1(M))/(2.d0*dt) + &
    & cw_m_row0*PEMF(M)*(pt_down_f(M)-0.5d0*(Tw1(M+1)+Tw1(M)) )
    flux2 = - lami*(Tis1(2)-Tis1(1))/(ddzi(1)*ls1) + RadDeepIce%integr(1) + &
    & ci_m_roi*dls0*(Tis1(2)-Tis1(1))/(2.d0*dt)
    dhwls = dt*(flux1 - flux2)/(row0_m_Lwi)
    dls = - dhwls*row0_d_roi
    dls0 = dls
    dhw = dhwhigh + dhwfsoil + dhwls + dhwlow + dhwtrib
    dhw0 = dhwhigh + dhwlow
  else
    ice = 1; snow = 0; water = 1; deepice = 0
    dhw = dhwhigh + dhwfsoil + dhwlow + dhwtrib
    dhw0 = dhwhigh + dhwlow
    dls = 0.
  endif

  if (soilswitch%par == 1) then
    call SOIL_COND_HEAT_COEF(nsoilcols)
  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, &
                     & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                     & bathymwater,bathymice, &
                     & bathymdice,bathymsoil(1,ix,iy), &
                     & soilflux(1,ix,iy), fdiffbot, (/1,0/))
    call LATERHEAT(ix,iy,gs,ls, &
    & bathymwater,bathymice,bathymdice,bathymsoil, &
    & gsp,soilflux(1,ix,iy),.false.,lsh)
  endif
  ! Calculation of the whole temperature profile
  call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
  & RadWater%extinct(1,1), extice, fetch, dt)


  lamsal(:) = (lamw(:) - lamw0)/cw_m_row0 * alsal + lamsal0
  lammeth(:) = lamw(:)/cw_m_row0 * almeth
  lamoxyg(:) = lamw(:)/cw_m_row0 * aloxyg
  lamcarbdi(:) = lamw(:)/cw_m_row0 * alcarbdi

  ! Salinity profile in water and soil
  call S_DIFF(dt,dhilow)
  ! Oxygen model precedes methane module, in order sedimentary 
  ! oxygen demand to be available for the latter
  call CHLOROPHYLLA(z_full, H_mixed_layer, H_photic_zone, extwat, M, Chl_a, itroph)
  call OXYGEN_PRODCONS(gs, gsp, wst, bathymsoil, gas, area_int, area_half, ddz05, &
  & h1, dt, Tsoil3, Chl_a, dzs, por, oxygsoil, itroph, &
  & prodox, resp, bod, sod, sodbot)

  if (soilswitch%par == 1) then
    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
  ! 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, &
                       & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                       & bathymwater,bathymice, &
                       & bathymdice,bathymsoil(1,ix,iy), &
                       & soilflux(1,ix,iy), fdiffbot, (/0,1/))
      methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
    endif
    ! Calling bubble model
    call BUBBLE_BLOCK
    ! Calculation of methane sources from methane fluxes from mutliple soil columns
    if (multsoil) then 
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,methsoilflux,.true.,lsm)
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,co2soilflux, .true.,lsc)
      call LATERHEAT(ix,iy,gs,ls, &
      & bathymwater,bathymice,bathymdice,bathymsoil, &
      & gsp,o2soilflux,  .true.,lso)
    endif
    gs%isoilcol = nsoilcols
    call METHANE &
    & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
    & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
    & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
    & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
    & lsm, bathymwater, &
    & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
    & ice_meth_oxid_total, &
    & h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !  call METHANE2 & ! two-meth
  !  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !  & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !  & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  else
    qwater(:,2) = qwater(:,1)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  endif
  call ADDOXPRODCONS(M,i_maxN,prodox,resp,bod,sod,dt,ls,gas)
  call OXYGEN (gs, Tw2, fbbleflx_o2_sum, fbbleflx_o2(0,nsoilcols), lamoxyg, gsp, fracb0, pressure, wind10, &
  & o2_pres_atm, ls, bathymwater, lso, dt, febul0(nsoilcols), eps1(1), sodbot, oxyg, soilswitch%par, fo2)
  call CARBON_DIOXIDE (gs, Tw2, fbbleflx_co2_sum, fbbleflx_co2(0,nsoilcols), lamcarbdi, gsp, bathymwater, lsc, fracb0, &
  & pressure, wind10, co2_pres_atm, febul0(nsoilcols), ls, dt, eps1(1), sodbot, gas, soilswitch%par, fco2)
  call METHANE_OXIDATION(M, i_maxN, dt, ddz, h1, bathymwater, oxyg(1,2), qwater(1,2), DIC(1,2), Tw2, &
  & qwateroxidtot, qwateroxidML)
!  call METHANE_OXIDATION2(M, dt, oxyg, qwater2(1,2,1), qwater2(1,2,2), DIC, Tw2) ! two-meth

!Calculation of water current's velocities and turbulent characteristics
  if (Turbpar%par /= 1) then
  
    if (Turbpar%par == 2) then
      call MOMENTUM_SOLVER(ix, iy, dnx, dny, ndatamax, &
      & year, month, day, hour, kor, a_veg, c_veg, h_veg, &
      & alphax, alphay, dt, b0, tau_air, tau_i, tau_gr, tau_wav, fetch, depth_area)  
    endif

    !Calculation of water density profile at the current timestep
    call DENSITY_W(M,eos%par,lindens%par,Tw2,Sal2,preswat,row2)

    if (Turbpar%par == 2) then
      !$OMP PARALLEL NUM_THREADS(num_ompthr) IF(omp%par == 1)
      do i = 1, nstep_keps%par
        call K_EPSILON(ix,iy,dnx,dny, &
        & year, month, day, hour, kor, a_veg, c_veg, h_veg, dt_keps, &
        & b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch) ! 0.5*tau_wav  0.5 - arbitrary (!) fraction of momentum loss due to wave breaking
      enddo
      !$OMP END PARALLEL
    else
!     Calculation of the of eddy difusivity for temperature
      SELECT CASE (Turbpar%par)
        CASE (8)
          ! Hostetler formulation
          ! call ED_TEMP_HOSTETLER &
          ! & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0.
        CASE (9)
          ! Modified Hostetler formulation for shallow lake
          ! call ED_TEMP_HOSTETLER2 &
          ! & (wind,velfrict_prev,zref,phi,z_full,row,KT,M)
          KT(:) = 0.
        CASE (1)
          do i = 1, M
            KT(i) = (k2(i)-niu_wat)*cw_m_row0 
          enddo
        CASE DEFAULT
          KT = lamTM*k2*cw_m_row0 
      ENDSELECT 
    endif
    
    do i = 1, M
      lamw(i) = lamw0 + KT(i) !+KC(i)
    enddo 
    
  else
    lamw = lamw0*3.d0
  endif
  

  h2 = h1 + dhw
  l2 = l1 + dhi
  ls2 = ls1 + dls
 
endif

!Diagnostic calculations
do i = 1, M
  k_turb_T_flux(i) = - 0.5*lamw(i)/(cw_m_row0) * &
  & ( Tw2(i+1) + Tw1(i+1) - Tw2(i) - Tw1(i)) / (ddz(i)*h1)
  T_massflux(i) = PEMF(i)*(pt_down_f(i)-0.5d0*(Tw2(i+1)+Tw2(i)))
enddo


!The scales of turbulence
call TURB_SCALES(gsp,ls, bathymwater, wst, RadWater, &
& k_turb_T_flux, T_massflux, row, eflux0_kinem_water, &
& turb_density_flux, Buoyancy0, tau_air-tau_wav, kor, i_maxN, H_mixed_layer,maxN, w_conv_scale, &
& T_conv_scale, Wedderburn, LakeNumber, Rossby_rad, ThermThick, ReTherm, RiTherm, trb)
 
ENDIF if2

!CASE 3: ICE WITHOUT WATER

if3: IF (layer_case == 3) THEN
    
!CASE 3.1: ICE WITH SNOW

if (flag_snow == 1) then

! Radiation fluxes
  SR_botsnow = shortwave*(1-albedoofsnow)
  if (flag_snow_init /= 1) then
    do i = itop, ms-1
      SR_botsnow = SR_botsnow*EXTINCT_SNOW(dens(i))**dz(i)
    enddo
  endif
  call RadIce%RAD_UPDATE(ddzi05(0:Mice)*l1, (/SR_botsnow/))

  x = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par)
  if (Ti1(1) > x + T_phase_threshold) then
    dhwhigh = (Ti1(1) - x - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi)
    dhif = - dhwhigh*row0_d_roi
    Ti1(1) = x + T_phase_threshold
  else
    dhwhigh = 0
    dhif = 0
  endif
  dhw = snmeltwat*snmelt*dt + dhwhigh + dhwtrib
  dhi = dhif
  dhi0 = dhi
  ice = 1; snow = 1; water = 0; deepice = 0

  call SNOW_COND_HEAT_COEF()
  if (soilswitch%par == 1) then
    call SOIL_COND_HEAT_COEF(nsoilcols)
  endif
  ! Calculation of heat sources from heat fluxes from mutliple soil columns
  if (multsoil) call LATERHEAT(ix,iy,gs,ls, &
  & bathymwater,bathymice,bathymdice,bathymsoil, &
  & gsp,soilflux(1,ix,iy),.false.,lsh)
  ! Temperature profile calculation
  call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
  & RadWater%extinct(1,1), extice, fetch, dt)

  ! Salinity profile
  call S_DIFF(dt,dhilow)
  if (soilswitch%par == 1) then
    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
  endif
  call SNOWTEMP(ix,iy,dnx,dny,year,month,day,hour,snowmass, &
  & snowmass_init,a,b,c,d,Temp,phi,extice,fetch,dt)
  if (soilswitch%par == 1) then
    call BUBBLE_BLOCK
    gs%isoilcol = nsoilcols
    call METHANE &
    & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
    & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
    & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
    & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
    & lsm, bathymwater, &
    & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
    & ice_meth_oxid_total, &
    & h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !  call METHANE2 & ! two-meth
  !  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !  & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !  & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  endif

  l2 = l1 + dhi 
  h2 = h1 + dhw
  ls2 = 0.

else

!CASE 3.2: ICE WITHOUT SNOW !

  ! Radiation fluxes
  Erad = shortwave*(1-albedoofice)
  call RadIce%RAD_UPDATE(ddzi05(0:Mice)*l1, (/Erad/))

  if (precip > 0 .and. tempair < 0) then
    flag_snow = 1
    hs1 = 2.*dzmin !0.02
    dhip = 0
  else
    dhip = precip*row0_d_roi*dt
  end if
  x = Meltpnt(0.e0_ireals,0.e0_ireals,nmeltpoint%par)
  if (Ti1(1) > x + T_phase_threshold) then
    dhw = (Ti1(1) - x - T_phase_threshold) * &
    & ci_m_roi*l1*ddzi(1)*0.5d0/(row0_m_Lwi) + dhwtrib
    dhihigh = - dhw*row0_d_roi
    dhis = 0
    Ti1(1) = x + T_phase_threshold
  else
    dhis = - Elatent/(roi*Liv)*dt
    dhw = 0.
    dhihigh = 0
  endif
  dhi = dhis + dhihigh + dhip
  dhi0 = dhi
  ice = 1; snow = 0; water = 0; deepice = 0 

  ! Temperature profile caluclation
  if (soilswitch%par == 1) then
    call SOIL_COND_HEAT_COEF(nsoilcols)
  endif
  ! Calculation of heat sources from heat fluxes from mutliple soil columns
  if (multsoil) call LATERHEAT(ix,iy,gs,ls, &
  & bathymwater,bathymice,bathymdice,bathymsoil, &
  & gsp,soilflux(1,ix,iy),.false.,lsh)
  call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
  & RadWater%extinct(1,1), extice, fetch, dt)

  ! Salinity profile
  call S_DIFF(dt,dhilow)
  if (soilswitch%par == 1) then
    call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
    call BUBBLE_BLOCK
    gs%isoilcol = nsoilcols
    call METHANE &
    & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
    & fbbleflx_ch4(0,nsoilcols), wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
    & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
    & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
    & lsm, bathymwater, &
    & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
    & plant_sum,bull_sum,oxid_sum,rprod_sum, &
    & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
    & ice_meth_oxid_total, &
    & h_talik,tot_ice_meth_bubbles,add_to_winter)
    qmethane(1:M+1) = qwater(1:M+1,2)
    qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !  call METHANE2 & ! two-meth
  !  & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !  & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !  & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !  & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !  & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !  & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !  qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !  qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
  endif
    
  l2 = l1 + dhi
  h2 = h1 + dhw
  ls2 = 0.
 
endif

ENDIF if3 

! CASE 4: SNOW AND SOIL WITHOUT ICE AND WATER
    
if4: IF (layer_case==4) THEN
  print*, 'The non-operational case: &
  &there is no water and ice layers: STOP'
  STOP
  if (flag_snow == 1) then
    ice = 0; snow = 1; water = 0; deepice = 0

    ! Temperature profile caluclation
    call SNOW_COND_HEAT_COEF()
    if (soilswitch%par == 1) then
      call SOIL_COND_HEAT_COEF(nsoilcols)
    endif
    call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
    & RadWater%extinct(1,1), extice, fetch, dt)

    ! Salinity profile
    call S_DIFF(dt,dhilow)
    if (soilswitch%par == 1) then
      call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
    endif
    call SNOWTEMP(ix,iy,dnx,dny,year,month,day,hour,snowmass, &
    & snowmass_init,a,b,c,d,Temp,phi,extice,fetch,dt)
    if (soilswitch%par == 1) then
      call BUBBLE_BLOCK 
      gs%isoilcol = nsoilcols
      call METHANE &
      & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
      & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
      & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
      & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
      & lsm, bathymwater, &
      & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
      & plant_sum,bull_sum,oxid_sum,rprod_sum, &
      & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
      & ice_meth_oxid_total, &
      & h_talik,tot_ice_meth_bubbles,add_to_winter)
      qmethane(1:M+1) = qwater(1:M+1,2)
      qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !    call METHANE2 & ! two-meth
  !    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !    & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !    & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !    & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !    & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !    & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !    qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !    qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth    
    endif
  else
    ice = 0; snow = 0; water = 0; deepice = 0
    
    ! Temperature profile calculation
    if (soilswitch%par == 1) then
      call SOIL_COND_HEAT_COEF(nsoilcols)
    endif
    call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
    & RadWater%extinct(1,1), extice, fetch, dt)

    ! Salnity profile
    call S_DIFF(dt,dhilow)
    if (soilswitch%par == 1) then
      call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
      call BUBBLE_BLOCK 
      gs%isoilcol = nsoilcols
      call METHANE(gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,nsoilcols),rosoil,fbbleflx_ch4_sum, &
      & fbbleflx_ch4(0,nsoilcols),wl4(1,nsoilcols),wi2(1,nsoilcols),wa,Sals2, &
      & rootss,rosdry,por,veg,0._ireals,TgrAnn, methgenmh, &
      & ddz,ddz05,wst, lammeth, h1, ls, dhw, dhw0, .true., .false., &
      & lsm, bathymwater, &
      & fplant, febul0(nsoilcols), fdiffbot(nsoilcols), ftot, fdiff_lake_surf, &
      & plant_sum,bull_sum,oxid_sum,rprod_sum, &
      & anox,gs,gsp,dt,eps1(1),sodbot,rprod_total_oldC,rprod_total_newC, &
      & ice_meth_oxid_total, &
      & h_talik,tot_ice_meth_bubbles,add_to_winter)
      qmethane(1:M+1) = qwater(1:M+1,2)
      qmethane(M+2:M+ns) = qsoil(2:ns,nsoilcols)
  !    call METHANE2 & ! two-meth
  !    & (pressure,wind10,zsoil,Tsoil3,wl4,wi2,wa,rootss,rosdry,por,veg,qsoil2,TgrAnn, & ! two-meth
  !    & ddz, Tw2, lammeth, qwater2, h1, l1, ls1, dhw, dhw0, & ! two-meth
  !    & fplant, febul2, fdiff2, ftot, fdiff_lake_surf2, & ! two-meth
  !    & plant_sum,bull_sum,oxid_sum,rprod_sum, & ! two-meth
  !    & anox,M,ns,dt,eps1(1),rprod_total_oldC,rprod_total_newC, & ! two-meth
  !    & h_talik,tot_ice_meth_bubbles2) ! two-meth
  !    qmethane(1:M+1) = qwater2(1:M+1,2,1) + qwater2(1:M+1,2,2) ! two-meth
  !    qmethane(M+2:M+ns) = qsoil2(2:ns,1) + qsoil2(2:ns,2) ! two-meth
    endif
  endif 
ENDIF if4


! Calculation of volume deficit (only if dead volume is specified), diagnostic
if (deadvol%par > 0. .and. h2 < deadvol%par) then
  voldef = voldef + area_int(1)*(deadvol%par - h2)
  h2 = deadvol%par
endif
   

if (h2 < min_water_thick .and. ls2 /= 0 .and. l2 == 0) then
  ls2 = ls2 + h2*row0_d_roi
  if (ls2 < min_ice_thick) then 
    ls2 = 0; Tis2 = 0
  endif
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0.
endif

if ((h2 < min_water_thick .and. h2 > 0) .and. l2 /= 0) then
  l2 = l2 + h2*row0_d_roi
  if (l2 < min_ice_thick) then 
    l2 = 0; Ti2 = 0
  endif
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0. 
end if
    
if (h2 < min_water_thick .and. l2 /= 0 .and. ls2 /= 0) then
  Ti2 = Ti2*l2/(l2+ls2) + Tis2*ls2/(l2+ls2)
  l2 = l2 + h2*row0_d_roi + ls2 
  ls2 = 0.
  Tis2 = 0.
  h2 = 0 
  Tw2 = 0.
  Sal2 = 0.
endif
   
if (l2 < min_ice_thick .and. h2 /= 0) then
  h2 = h2 + l2*roi/row0
  if (h2 < min_water_thick) then 
    h2 = 0; Tw2 = 0; Sal2 = 0.
  endif
  l2 = 0.
  Ti2 = 0.
end if

if (h2 < min_water_thick .and. l2 == 0) then
  h2 = 0.
  Tw2 = 0.
  Sal2 = 0. 
end if

if (l2 < min_ice_thick .and. h2 == 0) then
  l2 = 0.
  Ti2 = 0.
end if
   
if ((hs1 < min_ice_thick .and. hs1 > 0) .or. &
& (hs1 > 0 .and. l2 == 0 .and. h2 /= 0)) then
! h2 = h2 + (totalprecips-totalevaps-totalmelts)
  h2 = h2 + snowmass/row0
  if (h2 < min_water_thick) then
    l2 = l2 + row0*h2/roi
    h2 = 0; Tw2 = 0; Sal2 = 0.
    if (l2 < min_ice_thick) then 
      l2 = 0; Ti2 = 0
    endif
  endif
  hs1 = 0.
  flag_snow = 0
  flag_snow_init = 1
end if

if (hs1==0 .and. flag_snow==1) then
  flag_snow = 0
  flag_snow_init = 1
endif
 
if (ls2 < min_ice_thick .and. h2 /= 0) then
  h2 = h2 + ls2*roi/row0
  if (h2 < min_water_thick) then 
    h2 = 0; Tw2 = 0; Sal2 = 0.
  endif
  ls2 = 0
  Tis2 = 0.
endif 

if (ls2 > 0 .and. h2 == 0) then
  l2 = ls2
  Ti2 = Tis2
  if (l2 < min_ice_thick) then 
    l2 = 0; Ti2 = 0
  endif
  ls2 = 0
  Tis2 = 0.
endif

if (outflpar%par == 2 .and. tribheat%par > 0) then
  print*, 'LAGREFFL not operational: STOP'
  STOP
  !call TENDCALC(dt, utend, vtend, Twtend, Saltend)
  !allocate(work1(1:M+1),work2(1:M+1))
  !work1(:) = 0.5*(u1(:) + u2(:))
  !work2(:) = 0.5*(v1(:) + v2(:))
  !call LAGREFFL(Ux_tribin(1,:),Uy_tribin(1,:), &
  !& T_tribin(1,:),Sal_tribin(1,:), &
  !& utend,vtend,Twtend,Saltend, &
  !& work1,work2,area_int,dt,M, &
  !& ueffl,veffl,Tweffl,Saleffl)
  !deallocate(work1,work2)
endif

!! Deepened temperature maximum diagnostics
!! Output of temperature equation terms
!if (firstcall)   then
!  open(1234,file='maxT.dat',status='unknown')
!  write(1234,*), 'nstep dT/dt -dT_{ML}/dt -ddz/dt/dzT flux1/dz -flux2/dz flux1 -flux2 &
!  &S1/dz -S2/dz S1 -S2 botheat z_{ML} dz T_mean  Tmax-T_{ML} locTmax'
!  allocate (work3(1:17)); work3 = 0.; nn = 0
!  allocate (work4(1:17)); work4 = 0.
!endif
!!write(12345,*) tau/(wind*wind*1.273)
!!if (mod(nstep,60*6) == 0_iintegers) then
!! Identifying the location of the heated layer
!worklog1 = .false.
!i = 1
!c1:do
!  i = i + 1
!  !Averaging over the mixed layer
!  vol_ = 0.
!  do k = 1, i-1
!    vol_ = vol_ + ddz05(k-1)*h1*area_int(k)
!  enddo
!  ttt = 0.
!  do k = 1, i-1
!    ttt = ttt + Tw2(k)*ddz05(k-1)*h1*area_int(k)
!  enddo
!  ttt = ttt/vol_
!  if (Tw2(i) > ttt + 0.01 .and. z_full(i) > 1.) then
!  !if (KT(i-1)/cw_m_row0 < min_diff*100.5 .and. Tw2(i+1) > Tw2(i) .and. z_full(i) > 1.) then
!    !Upper boundary
!    x = 0.5*(z_full(i-1) + z_full(i))
!    j = i
!    c2:do 
!      j = j + 1
!      if (j > M+1) exit c2
!      if (Tw2(j) < ttt) then
!        !Lower boundary
!        xx = 0.5*(z_full(j) + z_full(j-1))
!        worklog1 = .true.
!        exit c2
!      endif
!    enddo c2
!    if (worklog1) then
!      ! Averaging over the temperature maximum layer
!      vol_2 = 0.
!      do k = i, j-1
!        vol_2 = vol_2 + ddz05(k-1)*h1*area_int(k)
!      enddo
!      zzz = 0.
!      do k = i, j-1
!        zzz = zzz + Tw2(k)*ddz05(k-1)*h1*area_int(k)
!      enddo
!      zzz = zzz/vol_2
!      zzz = zzz - ttt
!    endif
!    i2 = i
!    j2 = j
!    exit c1
!  endif
!  if (i == M+1) exit c1
!enddo c1
!worklog = .false.
!i = 1
!c3:do
!  i = i + 1
!  !Averaging over the mixed layer
!  vol_ = 0.
!  do k = 1, i-1
!    vol_ = vol_ + ddz05(k-1)*h1*area_int(k)
!  enddo
!  tt = 0.
!  do k = 1, i-1
!    tt = tt + Tw1(k)*ddz05(k-1)*h1*area_int(k)
!  enddo
!  tt = tt/vol_
!  if (Tw1(i) > tt + 0.01 .and. z_full(i) > 1.) then
!  !if (KT(i-1)/cw_m_row0 < min_diff*100.5 .and. Tw1(i+1) > Tw1(i) .and. z_full(i) > 1.) then
!    !Upper boundary
!    y = 0.5*(z_full(i-1) + z_full(i))
!    j = i
!    c4:do 
!      j = j + 1
!      if (j > M+1) exit c4
!      if (Tw1(j) < tt) then
!        !Lower boundary
!        yy = 0.5*(z_full(j) + z_full(j-1))
!        worklog = .true.
!        exit c4
!      endif
!    enddo c4
!    if (worklog) then
!      ! Averaging over the temperature maximum layer
!      vol_ = 0.
!      do k = i, j-1
!        vol_ = vol_ + ddz05(k-1)*h1*area_int(k)
!      enddo
!      zz = 0.
!      do k = i, j-1
!        zz = zz + Tw1(k)*ddz05(k-1)*h1*area_int(k)
!      enddo
!      zz = zz/vol_
!      zz = zz - tt
!    endif
!    exit c3
!  endif
!  if (i == M+1) exit c3
!enddo c3
!
!if (worklog .and. worklog1 .and. (j-1 > i+1) .and. yy < 4.) then
!
!  !Initializing TM mean temperature
!  if (.not. indTMprev) then
!    work3(1) = zz 
!    work3(2:12) = 0.
!    !print*, zz
!    !print*, Tw1(i:j-1) - tt
!    !print*, Tw1
!    !read*
!    indTMprev = .true.
!  endif
!
!  ! Volume change term
!  ww1 = 0.
!  if (vol_2 /= vol_) then
!    !print*, i,i2,j,j2
!    if (i2 > i) then
!      do k = i, i2-1
!        ww1 = ww1 - (Tw2(k) - ttt)*ddz05(k-1)*h1*area_int(k)
!      enddo
!    elseif (i2 < i) then
!      do k = i2, i-1
!        ww1 = ww1 + (Tw2(k) - ttt)*ddz05(k-1)*h1*area_int(k)
!      enddo
!    endif
!    if (j2 > j) then
!      do k = j, j2-1
!        ww1 = ww1 + (Tw2(k) - ttt)*ddz05(k-1)*h1*area_int(k)
!      enddo
!    elseif (j2 < j) then
!      do k = j2, j-1
!        ww1 = ww1 - (Tw2(k) - ttt)*ddz05(k-1)*h1*area_int(k)
!      enddo
!    endif
!  endif
!
!  ! Summing increments of temperature in the maximum
!  nn = nn + 1
!  work3(1) = work3(1) + (zzz - zz)
!  work3(2) = work3(2) - (ttt - tt)
!  work3(3) = work3(3) - (vol_2 - vol_)/(vol_)*zzz + ww1/(vol_)
!  work3(4) = work3(4) + 1./vol_*area_half(i-1)*k_turb_T_flux(i-1)*dt
!  work3(5) = work3(5) - 1./vol_*area_half(j-1)*k_turb_T_flux(j-1)*dt
!  work3(6) = work3(6) + k_turb_T_flux(i-1)*dt
!  work3(7) = work3(7) - k_turb_T_flux(j-1)*dt
!  work3(8) = work3(8) + 1./vol_*area_half(i-1)*SR(i-1)/cw_m_row0*dt
!  work3(9) = work3(9) - 1./vol_*area_half(j-1)*SR(j-1)/cw_m_row0*dt
!  work3(10) = work3(10) + SR(i-1)/cw_m_row0*dt
!  work3(11) = work3(11) - SR(j-1)/cw_m_row0*dt
!  ww = 0.
!  do k = i, j-1
!    ww = ww + lsh%water(k)*area_int(k)*ddz05(k-1)
!  enddo
!  work3(12) = work3(12) + ww*h1/cw_m_row0/vol_*dt
!  work3(13) = work3(13) + y
!  work3(14) = work3(14) + yy-y
!  work3(15) = work3(15) + zz
!  work3(16) = work3(16) + maxval(Tw1(i:j-1))-tt
!  work3(17) = work3(17) + (maxloc(Tw1,1) - i + 0.5)*ddz(i)*h1 !Assuming regular mesh
!
!  ! Time averaging
!  work4(1:17) = work4(1:17) + work3(1:17)
!
!  !www = (zzz - zz)/dt + (ttt - tt)/dt + (vol_2 - vol_)/(dt*vol_)*zzz - ww1/(vol_*dt) &
!  !& - 1./vol_*area_half(i-1)*k_turb_T_flux(i-1) &
!  !& + 1./vol_*area_half(j-1)*k_turb_T_flux(j-1) &
!  !& - 1./vol_*area_half(i-1)*SR(i-1)/cw_m_row0 &
!  !& + 1./vol_*area_half(j-1)*SR(j-1)/cw_m_row0 &
!  !& - ww*h1/cw_m_row0/vol_
!  !if (abs(www/(zzz - zz)*dt) > 1.E-3) then
!  !  print*, (zzz - zz)/dt, www, (vol_2 - vol_)/(dt*vol_)*zzz, - ww1/(vol_*dt) 
!  !  read*
!  !endif
!
!  if (abs(work3(1) - zzz) > 1.e-10) then
!    print*, 'TMdiag', nstep, zz, zzz, work3(1) - zzz
!    read*
!  endif
!
!  if (mod(nstep,60*6) == 0_iintegers) then
!    write(1234,*) nstep, (zzz - zz)/dt, - (ttt - tt)/dt, &
!    !& - (area_half(j-1)*(xx - yy) - area_half(i-1)*(x - y))/(dt*vol_)*zz, &
!    & - (vol_2 - vol_)/(dt*vol_)*zzz + ww1/(vol_*dt), &
!    & + 1./vol_*area_half(i-1)*k_turb_T_flux(i-1), &
!    & - 1./vol_*area_half(j-1)*k_turb_T_flux(j-1), &
!    & + k_turb_T_flux(i-1), - k_turb_T_flux(j-1), &
!    &   1./vol_*area_half(i-1)*SR(i-1)/cw_m_row0, &
!    & - 1./vol_*area_half(j-1)*SR(j-1)/cw_m_row0, &
!    &   SR(i-1)/cw_m_row0, - SR(j-1)/cw_m_row0, &
!    &   ww*h1/cw_m_row0/vol_, &
!    &   y, yy-y, zz, maxval(Tw1(i:j-1))-tt, &
!    &   (maxloc(Tw1(i:j-1),1)-0.5)*ddz(i)*h1 !Assuming regular mesh
!  endif 
!
!else
!
!  indTMprev = .false.
!
!endif
!
!if (step_final) then
!  open(12345,file='maxTav.dat',status='unknown')
!  write(12345,*), 'nstep dT/dt -dT_{ML}/dt -ddz/dt/dzT flux1/dz -flux2/dz flux1 -flux2 &
!  &S1/dz -S2/dz S1 -S2 botheat z_{ML} dz T_mean  Tmax-T_{ML} locTmax'
!  write(12345,*) nstep, work3(1:17)/real(nn)
!  write(12345,*) nstep, work3(1:17)
!  write(12345,*) nstep, work4(1:17)/real(nn)
!  !print*, 'resid = ', work3(1) - work3(2) - work3(3) - work3(4) - work3(5) - work3(8) - work3(9) - work3(12)
!  close(12345)
!  deallocate(work3,work4)
!endif
!!endif

!i = maxloc(Tw2,1)
!if (i > 1 .and. i < M+1 .and. (mod(nstep,60) == 0_iintegers) ) then
!  write(1234,*) cw_m_row0*(Tw2(i) - Tw1(i))/dt,  &
!  & 0.5/h1**2*( lamw(i)*(Tw2(i+1) + Tw1(i+1) - Tw2(i) - Tw1(i) )/ddz(i) - &
!  & lamw(i-1)*( Tw2(i) + Tw1(i) - Tw2(i-1) - Tw1(i-1) )/ddz(i-1) )/ddz05(i-1), &
!  & lsh%water(i), (SR(i) - SR(i+1))/(h1*ddz05(i-1))
!endif

!! On-screen diagnostics of TKE budget
!i = M-40 !Level in thermocline
!if (month == 10 .and. nn == 0) then
!  xx = tau/row0*u1(1)
!  yy = - H_mixed_layer*lamw0/cw_m_row0*g/row0*(row(i+1)-row(i))/(ddz(i)*h1)
!  zz = - sum(eps1(1:i_maxN)*ddz(1:i_maxN))*h1
!  nn = 1
!elseif (nn > 0) then
!  nn = nn + 1
!  xx  = ACCUMM(nn, xx, tau/row0*u1(1))
!  yy  = ACCUMM(nn, yy, - H_mixed_layer*lamw0/cw_m_row0*g/row0*(row(i+1)-row(i))/(ddz(i)*h1))
!  zz  = ACCUMM(nn, zz, - sum(eps1(1:i_maxN)*ddz(1:i_maxN))*h1)
!  !print*, xx,yy,zz
!  !read*
!endif
!if (mod(nstep,100) == 0) then
!  !i = M-40 !Level in thermocline
!  print*, tau/row0*u1(1), -H_mixed_layer*lamw0/cw_m_row0*g/row0*(row(i+1)-row(i))/(ddz(i)*h1), - &
!  & sum(eps1(1:i_maxN)*ddz(1:i_maxN))*h1, &
!  & tau/row0*u1(1) - H_mixed_layer*lamw0/cw_m_row0*g/row0*(row(i+1)-row(i))/(ddz(i)*h1) - &
!  & sum(eps1(1:i_maxN)*ddz(1:i_maxN))*h1
!  print*, xx + yy + zz
!endif

h1    = h2
l1    = l2
ls1   = ls2
Tw1   = Tw2
Ti1   = Ti2
Tis1  = Tis2
Sal1  = Sal2
Sals1 = Sals2
qwater(:,1) = qwater(:,2)
!qwater2(:,1,1) = qwater2(:,2,1) ! two-meth
!qwater2(:,1,2) = qwater2(:,2,2) ! two-meth
oxyg  (:,1) = oxyg  (:,2)
DIC(:,1) = DIC(:,2)
u1    = u2
v1    = v2
Tskin(1) = Tskin(2)
    
! CALCULATION OF SUMMARY FLUXES !
totalpen = totalpen + dhwfsoil 
totalevap = totalevap + Elatent/(row0*Lwv)*dt
totalprecip = totalprecip + precip*dt
totalhflux = totalhflux + hflux*dt
totalerad = totalerad + erad*dt

if (l1 == 0) tsw = Tw1(1) + 273.15
if (l1 /=0 .and. flag_snow == 0) tsw = Ti1(1) + 273.15
if (flag_snow == 1) tsw = T(itop) + 273.15

if (init(ix,iy) == 0) then
  init(ix,iy) = 1  
endif

!VALUES AT NEXT TIME STEP (time + dt) IN CURRENT POINT (ix,iy)

call UPDATE_NEXT_TIMESTEP ( &
& ix, iy, dnx, dny, M, Mice, ns, ms, ml, nsoilcols, &
& l1, h1, hx1, hx2, hy1, hy2, ls1, hs1, &
& hx1t, hx2t, hy1t, hy2t, &
& hx1ml, hx2ml, hy1ml, hy2ml, &
& u1, v1, &
& E1, eps1, k_turb_T_flux, &
& Tsoil1, Sals1, wi1, wl1, &
& Tw1, Sal1, lamw, &
& Tskin(1), &
& Ti1, Tis1, &
& ueffl, veffl, Tweffl, Saleffl, &
& dz, T, wl, dens, &
& qwater(1,1), qsoil, &
& oxyg(1,1), oxygsoil, &
& DIC(1,1), DOC, POCL, POCD, &
& snmelt, snowmass, &
& cdmw, &
& time, &
& dhwfsoil, &
& Elatent, &
& dhw, dhw0, &
& dhi, dhi0, dls0, &
& velfrict, &
& roughness, &
& eflux0_kinem, &
& tot_ice_meth_bubbles, &
& febul0, Eseiches, &

& flag_snow, flag_snow_init, &
& itop, &
& nstep, i_maxN, itherm )

!At the next timestep after control point output, the control point may be written again
!if (ix == nx - nx0 + 1 .and. iy == ny - ny0 + 1 .and. &
!&  cpwrite_done .eqv. .true.) cpwrite_done = .false.
!Writing the control point for all lakes of the domain
if (icp == 1_iintegers .and. ix == nx - nx0 + 1 .and. iy == ny - ny0 + 1) then
  call CONTROL_POINT_OUT(nx,nx0,nx_max,ny,ny0,ny_max,gs,parparams)
endif

hw1 = hw
xlew1 = xlew
cdmw1 = cdmw
surfrad1 = surfrad
h2_out = h2


!Diagnoctics

!call fluxes
!print*,h1,l1,hs1,ls1,Tw2(1),Sal2(1),Sal2(M+1)
!if (ix == 10 .and. iy == 10) print*, 'Lake', &
!& 'T1=',Tw1(1),'T2=',Tw1(2),'h1=',h1,'S=',shortwave, &
!& 'hh=',hflux,'LE=',Elatent, &
!& 'eflux=',eflux, 'Longwave=', longwave, &
!& 'Bal=', (shortwave*sabs + &
!&  shortwave*(1 - sabs)*(1 - albedoofwater)*(1 - exp( - extwat*0.5*ddz(1)*h1)) + &
!&  longwave*(1 - albedoofwater_lw) - &
!&  surfrad - hflux - Elatent - eflux)/(cw_m_row0*ddz(1)*h1/2)*dt

if (accum .eqv. .false. .and. &
&   year  == year_accum_begin  .and. &
&   month == month_accum_begin .and. &
&   day   == day_accum_begin   .and. &
&   hour  >= hour_accum_begin) then
  accum = .true.
  totmeth0 = VARMEAN(qwater(1,1),bathymwater,11_iintegers) * &
  & (h2 - dhw)*dzeta_05int(M)*mf ! Depth at previous timestep
endif

if (accum .eqv. .true. .and. &
&   year  == year_accum_end  .and. &
&   month == month_accum_end .and. &
&   day   == day_accum_end   .and. &
&   hour  >= hour_accum_end) accum = .false.

if (accum) then
! Note: variables accumulation is implemented currently only for one-point simulations!
  call ACCUM_VAR &
  & (dt, l1, fbbleflx_ch4_sum(0), fbbleflx_ch4_sum(M+1), & 
  & fdiffbot(nsoilcols), fdiff_lake_surf, qwateroxidtot, &
  & ice_meth_oxid_total, &
  & rprod_total_newC(nsoilcols), rprod_total_oldC(nsoilcols), &
  & febultot, febulbottot, fdifftot, fdiff_lake_surftot, &
  & rprod_total_newC_integr, rprod_total_oldC_integr, &
  & methoxidwat, ice_meth_oxid, add_to_winter)
  totmethc = VARMEAN(qwater(1,2),bathymwater,11_iintegers) * &
  & h2**dzeta_05int(M)*mf
!  metracc = metracc + qwater(M+1,2)*(dhw - dhw0) + dhw0*qwater(1,2)
!  call ACCUM_VAR & ! two-meth
!  & (dt, l1, febul2, fdiff2, fdiff_lake_surf2, & ! two-meth
!  & rprod_total_newC, rprod_total_oldC, & ! two-meth
!  & febultot2, fdifftot2, fdiff_lake_surftot2, & ! two-meth
!  & rprod_total_newC_integr2, rprod_total_oldC_integr2) ! two-meth
endif

if (flag_print) then ! The output in ASCII files 
  if (monthly_out%par == 1) call MON_OUT(ix,iy,dnx,dny,year,month,day,hour,time,zgrid_out,ngrid_out%par)
  if (daily_out%par   == 1) call DAY_OUT(ix,iy,dnx,dny,year,month,day,hour)
  if (hourly_out%par  == 1) call HOUR_OUT(ix,iy,dnx,dny,year,month,day,hour,time)
  if (everystep%par   >  0) call EVERYSTEP_OUT(ix,iy,dnx,dny)
  if (time_series%par == 1) then
    call SERIES_OUT(ix,iy,dnx,dny,year,month,day,hour,tsw)
    allocate(work1(1:M+ns))
    ndec = -1 ! Number of decimal digits, negative meaning exponential format
    i = 1
    !Interpolating diffusivity to layers' interfaces, result is stored in work1
    call LININTERPOL (z_half,lamsal,M,z_full,work1,M+1,flag) 
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  Tw1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, & ! ngrid_out%par
    &  outpath, 'water_temp', i, ndec, .false.)  !, row, work1 - last optional argument!
    ndec = 4
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  Tsoil3, zsoil, ns, &
    &  zgridsoil_out, ngridsoil_out%par, & !ngridsoil_out%par, &
    &  outpath, 'soil_temp', i, ndec, .false.)
    z_watersoil(1:M+1) = z_full(1:M+1)    
    do j = 2, ns
      z_watersoil(M+j) = h1 + zsoil(j)
    enddo
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  rowc, z_full, M+1, &
    &  zgrid_out, -1, &
    &  outpath, 'water_dens_layers', i, ndec, .false.)
    ndec = -1
    i = i + 2
    work1(1:M+ns) = qmethane(1:M+ns)*molm3tonM
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_watersoil, M+ns, &
    &  zgridsoil_out, -1, &
    &  outpath, 'methane_water_soil', i, ndec, .false.)
    ndec = -1
    i = i + 2
    work1(1:M+1) = oxyg(1:M+1,1)*molm3tomgl_o2 !molm3tonM
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'oxygen_water', i, ndec, .false.)
    ndec = -1
    i = i + 2
    work1(1:M+1) = qwater(1:M+1,1)*molm3tomcM !molm3tomkgl_ch4 ! molm3tonM
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'methane_water', i, ndec, .false.)
    if (nsoilcols > 1) then
      ndec = -1
      i = i + 2
      ! Methane production in soil columns
      work1(1:nsoilcols) = 0.5*(zsoilcols(1:nsoilcols,ix,iy) + zsoilcols(2:nsoilcols+1,ix,iy))
      call PROFILE_OUTPUT &
      & (ix, iy, dnx, dny, &
      &  year, month, day, hour, &
      &  time, dt_out%par, &
      &  rprod_total_newC, work1, nsoilcols, &
      &  work1, nsoilcols, &
      &  outpath, 'rprod_soil', i, ndec, .false.)
    endif 
    ndec = -1
    i = i + 2
    do j = 1, M+1
      work1(j) = DIC(j,1) / HC_CORR_CARBEQUIL(Tw1(j)+Kelvin0,pH) * molm3tomcM !molm3tomgl_co2 !molm3tonM
    enddo
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'co2_water', i, ndec, .false.)
    ndec = -1
    i = i + 2    
    work1(1:M) = fbbleflx_ch4_sum(1:M) 
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_half, M, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'fbblflx_ch4', i, ndec, .false.)
    ndec = -1
    i = i + 2
    work1(1:M) = fbbleflx_co2_sum(1:M)
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_half, M, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'fbblflx_co2', i, ndec, .false.)
    ndec = -1
    i = i + 2
    work1(1:M) = fbbleflx_o2_sum(1:M)
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  work1, z_half, M, &
    &  zgrid_out, ngrid_out%par, &
    &  outpath, 'fbblflx_o2', i, ndec, .false.)
    deallocate(work1)
    ndec = 4
    i = i + 2    
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  wl1, zsoil, ns, &
    &  zgridsoil_out, ngridsoil_out%par, &
    &  outpath, 'wl_soil', i, ndec, .false.)
    ndec = 4
    i = i + 2    
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  wi1, zsoil, ns, &
    &  zgridsoil_out, ngridsoil_out%par, &
    &  outpath, 'wi_soil', i, ndec, .false.)    
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  Sal1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'sal_water', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  lamw/cw_m_row0, z_half, M, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'lamw', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  eps1, z_half, M, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'epsilon', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  Ri, z_half, M, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'Ri', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  trb%Rp, z_half, M, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'Rp', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  trb%Rpdens, z_half, M, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'Rpdens', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  u1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'u', i, ndec, .false.)
    ndec = -1
    i = i + 2
    call PROFILE_OUTPUT &
    & (ix, iy, dnx, dny, &
    &  year, month, day, hour, &
    &  time, dt_out%par, &
    &  v1, z_full, M+1, &
    &  zgrid_out, ngrid_out%par, & !ngrid_out
    &  outpath, 'v', i, ndec, .true.)
  endif

  ! Output time series of scalars in selected points
  if ( rtemp%arr1(1,1) /= missing_value .and. dyn_pgrad%par == 4 ) then 
      call TEMPLOC(M,3_iintegers,1_iintegers,itherm,h1,dt_out%par,time,hx1ml,hy1ml,Tw1,gsp,gs,bathymwater,rtemp, &
    &            'temploc',firstcall)
      call TEMPLOC(M,3_iintegers,2_iintegers,itherm,h1,dt_out%par,time,hx1ml,hy1ml,qwater(1,1),gsp,gs,bathymwater,rtemp, &
    &            'ch4loc',firstcall)
      call TEMPLOC(M,3_iintegers,3_iintegers,itherm,h1,dt_out%par,time,hx1ml,hy1ml,DIC(1,1),gsp,gs,bathymwater,rtemp, &
    &            'co2loc',firstcall)
  endif

endif
 

if (runmode%par == 1)  then
  if (nscreen%par > 0 .and. mod(nstep,nscreen%par)==0) then

    write (*,'(3a25)') 'timestep', 'N of point', &
    & 'Surface temperature'
    write (*,'(2i25,f25.1)') nstep, ix, tsw
    write (*,'(5a8)') 'Year', 'Month', 'Day', 'Hour', 'Months'
    i = year 
    if (year<1000) i = year+1000
    write (*,'(3i8,2f8.2)') i, month, day, hour, time/(month_sec)

    write(*,'(a)')
    write(*,*) 'Methane generation, transport and oxidation rates'
    write(*,*) 'Total bottom diff, Total surface diff, Total bottom ebul, Total surface ebul'
    write(*,'(a,4(2x,e12.5))') 'Summer ', fdifftot(1), fdiff_lake_surftot(1), febulbottot(1), febultot(1)
    write(*,'(a,4(2x,e12.5))') 'Winter ', fdifftot(2), fdiff_lake_surftot(2), febulbottot(2), febultot(2)
    write(*,*) 'Total prod from "new" org, Total prod from "old" org, Total oxidation in water, &
    & Total oxidation in ice (for winter)'
    write(*,'(a,3(2x,e12.5))') 'Summer ', &
    & rprod_total_newC_integr(1), rprod_total_oldC_integr(1), methoxidwat(1)
    write(*,'(a,4(2x,e12.5))') 'Winter ', &
    & rprod_total_newC_integr(2), rprod_total_oldC_integr(2), methoxidwat(2), ice_meth_oxid
    write(*,'(a,(2x,e12.5))') 'Change of integral methane amount in water column and ice cover', &
    & totmethc - totmeth0 + tot_ice_meth_bubbles*mf
    write(*,'(a)')
    write(*,*) 'Checking methane balance residuals:'   
    write(*,'(a,1(2x,e12.5))') 'Total production in soil - total bottom flux ', &
    & sum(rprod_total_newC_integr(1:2)) + sum(rprod_total_oldC_integr(1:2)) - &
    & sum(febulbottot(1:2)) - sum(fdifftot(1:2))
    write(*,'(a,1(2x,e12.5))') 'Total bottom flux - total surface flux - total oxidation - & 
    & methane total column amount change', &
    & sum(febulbottot(1:2)) + sum(fdifftot(1:2)) - &
    & sum(febultot(1:2)) - sum(fdiff_lake_surftot(1:2)) - &
    & sum(methoxidwat(1:2)) - ice_meth_oxid - &
    & (totmethc - totmeth0 + tot_ice_meth_bubbles*mf)


    open(12345,file=outpath(1:len_trim(outpath))//'lastscrout.dat')
    write(12345,'(a)')
    write(12345,*) 'Methane generation, transport and oxidation rates'
    write(12345,*) 'Total bottom diff, Total surface diff, Total bottom ebul, Total surface ebul'
    write(12345,'(a,4(2x,e12.5))') 'Summer ', fdifftot(1), fdiff_lake_surftot(1), febulbottot(1), febultot(1)
    write(12345,'(a,4(2x,e12.5))') 'Winter ', fdifftot(2), fdiff_lake_surftot(2), febulbottot(2), febultot(2)
    write(12345,*) 'Total prod from "new" org, Total prod from "old" org, Total oxidation in water, &
    & Total oxidation in ice (for winter)'    
    write(12345,'(a,3(2x,e12.5))') 'Summer ', &
    & rprod_total_newC_integr(1), rprod_total_oldC_integr(1), methoxidwat(1)
    write(12345,'(a,4(2x,e12.5))') 'Winter ', &
    & rprod_total_newC_integr(2), rprod_total_oldC_integr(2), methoxidwat(2), ice_meth_oxid
    write(12345,'(a,(2x,e12.5))') 'Change of integral methane amount in water column and ice cover', &
    & totmethc - totmeth0 + tot_ice_meth_bubbles*mf
    write(12345,'(a)')
    write(12345,*) 'Checking methane balance residuals:'   
    write(12345,'(a,1(2x,e12.5))') 'Total production in soil - total bottom flux ', &
    & sum(rprod_total_newC_integr(1:2)) + sum(rprod_total_oldC_integr(1:2)) - &
    & sum(febulbottot(1:2)) - sum(fdifftot(1:2))
    write(12345,'(a,1(2x,e12.5))') 'Total bottom flux - total surface flux - total oxidation - & 
    & methane total column amount change', &
    & sum(febulbottot(1:2)) + sum(fdifftot(1:2)) - &
    & sum(febultot(1:2)) - sum(fdiff_lake_surftot(1:2)) - &
    & sum(methoxidwat(1:2)) - ice_meth_oxid - &
    & (totmethc - totmeth0 + tot_ice_meth_bubbles*mf)
    close(12345)

  endif
endif
  

!write(workchar,'(i5)') ncomp
!workchar = '(a6,'//workchar(1:len_trim(workchar))//'(1x,f6.2))'
!write(*,fmt = workchar) 'Timing', comptime(1:ncomp)/60. ! Converting to minutes
!if (l1 > 0.) print*, T(itop), Ti1(1), Ti1(Mice+1)

if (firstcall) firstcall = .false.

!FORMATS!
7   format (f7.3, 4i5,37f7.3) 
60  format (f5.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
! 61  format (i5, f6.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
61  format (f6.2, f5.2, 2e16.5, 2f8.3, e15.5, f11.5, f7.1)
62  format (f6.2, f5.2, 2f12.3, 2f8.3, f11.5, f7.1,3f12.4,2f7.2,3f7.1)
80  format (f7.3, 10f9.2)
90  format (f7.3, 12f9.2)   
100 format (a5, 2a16, 2a8, a15, a11, a7) 


contains
SUBROUTINE BUBBLE_BLOCK
! The block of code invoking procedures, that calculate the bubble flux
implicit none

! Calling bubble model
fbbleflx_ch4_sum(0:M+1) = 0.
fbbleflx_o2_sum (0:M+1) = 0.
fbbleflx_co2_sum(0:M+1) = 0.
if (ifbubble%par == 1) then
  call BUBBLE(M,M+1,ddz,ddz05,dzeta_int,qwater(:,1),Tw1,Sal1,h1,pressure, &
  & DIC(:,1),oxyg(:,1),ngasb,febul0(nsoilcols),dt, &
  & fracb0,fbbleflx_ch4(0,nsoilcols), &
  & fbbleflx_co2(0,nsoilcols),fbbleflx_o2(0,nsoilcols))
endif 
if (multsoil) then
  do i = 1, nsoilcols-1
    if (ifbubble%par == 1) then
      call BUBBLE(M,bathymsoil(i,ix,iy)%ibot,ddz,ddz05,dzeta_int, &
      & qwater(:,1),Tw1,Sal1,z_half(isoilcolc(i)),pressure, &
      & DIC(:,1),oxyg(:,1),ngasb,febul0(i),dt, &
      & fracb0,fbbleflx_ch4(0,i), &
      & fbbleflx_co2(0,i),fbbleflx_o2(0,i)) 
    else
      fbbleflx_ch4(1:bathymsoil(i,ix,iy)%ibot,i) = 1.; fbbleflx_ch4(bathymsoil(i,ix,iy)%ibot+1:M+1,i) = 0.
      fbbleflx_co2(1:bathymsoil(i,ix,iy)%ibot,i) = 1.; fbbleflx_co2(bathymsoil(i,ix,iy)%ibot+1:M+1,i) = 0.
      fbbleflx_o2 (1:bathymsoil(i,ix,iy)%ibot,i) = 1.; fbbleflx_o2 (bathymsoil(i,ix,iy)%ibot+1:M+1,i) = 0. 
    endif
    ! Adding bubble flux from i-th soil column to the horizontally averaged bubble flux
    call BUBBLEFLUXAVER(ix,iy,i)
    ! Adding bubble fluxes to fluxes of CH_4, CO_2 and O_2 at soil columns tops
    ! Sedimentary oxygen demand is added to CO_2 and O_2 as sources in OXYGEN_MOD module
    methsoilflux(i) =    methsoilflux(i) - febul0(i)
    co2soilflux (i) =  - febul0(i)*fracb0(2)/fracb0(1)
    o2soilflux  (i) =  - febul0(i)*fracb0(3)/fracb0(1)
  enddo
endif

END SUBROUTINE BUBBLE_BLOCK

END SUBROUTINE LAKE


END MODULE LAKE_MOD
