lake.f90

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
use METH_OXYG_CONSTANTS, only : SET_METH_OXYG_CONSTANTS

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 SET_METH_OXYG_CONSTANTS
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,netrad1,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,SurfTemp1,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, only : UnDenGrMix
use SALINITY, only : S_DIFF, UPDATE_ICE_SALINITY
use TURB, only : K_EPSILON, ED_TEMP_HOSTETLER, ED_TEMP_HOSTETLER2, ZIRIVMODEL
use TURB_CONST, only : min_diff
use DZETA_MOD, only : VARMEAN

use MOMENTUM, only : &
& MOMENTUM_SOLVER

use VERTADV, only : &
& VERTADV_UPDATE

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, &
& DIFFMIN_HS

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

use TURB_CONST, only : &
& lamTM, min_diff

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

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

use CARBON_DIOXIDE_MOD, only : &
& CARBON_DIOXIDE

use PHOSPHORUS_MOD, only : &
& PHOSPHORUS

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, &
& fracbands, nbands, extwatbands

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
!> Net radiation, W/m**2
real(kind=ireals), intent(in) :: netrad1
!> 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
!> Surface temperature, K
real(kind=ireals), intent(in) :: SurfTemp1

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(inout) :: 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 (CP) switch: 0 - do nothing
!!                       1 - write CP                             
!!                       2 - read CP as 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), allocatable :: radflux(:)
real(kind=ireals) :: l2, h2, ls2
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) :: 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_methprod_phot and r0_oldorg2_star
!      calibr_par1 => r0_methprod_phot
!      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
uwind = uwind1
vwind = vwind1
zref  = zref1
precip = precip1
Sflux0 = Sflux01*(roa0/row0)
outpath = outpath1
cloud = cloud1
if (SurfTemp1 /= missing_value) then
  SurfTemp = SurfTemp1 - Kelvin0
else
  SurfTemp = missing_value
endif

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
allocate(radflux(1:nbands))
if (nbands > 1) then
  forall (i=1:M+1) RadWater%extinct(i,1:nbands) = extwatbands(1:nbands)
  RadWater  %extinct(:,2) = extwat ! extwat is treated as extinction coefficient for PAR
else
  RadWater  %extinct(:,:) = extwat
endif
RadIce    %extinct(:,:) = extice
RadDeepIce%extinct(:,:) = extice
if (ifrad%par == 1) then
  longwave  = longwave1
  shortwave = max(shortwave1,0._ireals)
elseif (ifrad%par == 0) then
  longwave  = 0.e0_ireals
  shortwave = 0.e0_ireals
endif
! If atmospheric radiation is missing, net radiation to be used to compute this variable
if (longwave1 /= missing_value) then
  netrad = missing_value
else
  netrad = netrad1
endif

!print*, tempair, humair, pressure, uwind, vwind, longwave, shortwave

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

if (uwind == 0) uwind = minwind
if (vwind == 0) vwind = minwind
wind = sqrt(uwind**2+vwind**2)
!wind10 = wind*log(10./roughness0)/log(zref/roughness0)
precip = 0.
!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, isp, &
  
  & 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, cdmw, 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), phosph, oxygsoil, Sal1, &
  & u1, v1, Tw1, Tskin, T_0dim, Eseiches, &
  & Ti1, salice, porice, 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)
  print*, 'Control point is read'
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, &
  & gradpx_tend, gradpy_tend, &
  & 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, phosph, &
  & snmelt, snowmass, &
  & cdmw, &
  & time, &
  & dhwfsoil, &
  & Elatent, &
  & dhw, dhw0, &
  & dhi, dhi0, dls0, &
  & velfrict_prev, &
  & roughness, &
  & eflux0_kinem, &
  & tot_ice_meth_bubbles, &
  & febul0, Eseiches, salice, porice, &
  
  & 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
  porice(:) = poricebot * saltice%par
  salice(:) = porice(:) * Sal1(1) * row0
  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) then
  call TRIBTEMP(time,dt,h1,dhwtrib,z_full,area_int,area_half,gsp,gas,wst,spinup_done)
  ! Adding tendency due to mean vertical velocity
  call VERTADV_UPDATE(M,dt,h1,gsp,bathymwater,wst,gas)
endif

!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
!Partitioning of shortwave energy between wavelength bands
forall(i=1:nbands) radflux(i) = fracbands(i)*Erad*(1 - sabs0)*x
call RadWater%RAD_UPDATE(ddz05(0:M)*h1,radflux) !(/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
  forall (i=1:nbands) radflux(i) = RadWater%flux(M+1,i)*(1-albedoofice)
  call RadDeepIce%RAD_UPDATE(ddzi05(0:Mice)*ls1, radflux) !(/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, RadIce, SurfTemp, fetch, dt)
if (soilswitch%par == 1) then
  call SOILFORLAKE(dt,a,b,c,d,nsoilcols)
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

!Assuming all dissolved species diffuse with the same molecular diffusivity as salinity
lamsal(:)    = (lamw(:) - lamw0)/cw_m_row0 * alsal    + lamsal0
lammeth(:)   = (lamw(:) - lamw0)/cw_m_row0 * almeth   + lamsal0
lamoxyg(:)   = (lamw(:) - lamw0)/cw_m_row0 * aloxyg   + lamsal0
lamcarbdi(:) = (lamw(:) - lamw0)/cw_m_row0 * alcarbdi + lamsal0
lampart  (:) = (lamw(:) - lamw0)/cw_m_row0 * alcarbdi !no molecular diffusion for particles

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

! 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, RadWater, gas, 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
  ! 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,Sals1, &
  & 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)