lake.f90

MODULE LAKE_MOD

type ISIMIP_Lake_type
  ! Type for ISIMIP-Lake output
  sequence
  integer :: istratstart, istratend, istratdur
  integer :: icestart, icetend, icedur
  real :: thermodepth
  real :: surftemp, bottemp
  real :: lakeicefrac
  real :: icethick
  real :: snowthick
  real :: waterthick
  real :: icetemp
  real :: snowtemp
  real :: sensheatf
  real :: latentheatf
  real :: momf
  real :: lwup
  real :: swup
  real :: lakeheatf
  real :: albedo
  real :: extcoeff
  real :: sedheatf
  real, allocatable :: watertemp(:)
  real, allocatable :: zwatertemp(:)
end type ISIMIP_Lake_type

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 ARRAYS_GRID, only : gs
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_LAKE_MOD, only : COMSOILFORLAKE
use SURF_SCHEME1_MOD, only : PBLDAT_LAKE
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
use ALLOCARRAYS_MOD, only : ALLOCARRAYS

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_LAKE
call PBLDAT_INM
call COMSOILFORLAKE
call TURB_CONST_DERIVED_SET
call SET_DENS_CONSTS(eos%par)
call NUMERIC_PARAMS_SET(gs)

if (error_cov==1) then
  if (runmode%par == 2) then
    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, fdiff_lake_surf1, 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, LakePhysISIMIP)

!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_LAKE_MOD, only : UnDenGrMix
use SALINITY_LAKE_MOD, only : S_DIFF, UPDATE_ICE_SALINITY
use TURB_LAKE_MOD, only : K_EPSILON, ED_TEMP_HOSTETLER, &
&                         ED_TEMP_HOSTETLER2, ZIRIVMODEL
use TURB_CONST, only : min_diff
use DZETA_MOD, only : VARMEAN

use MOMENTUM_LAKE_MOD, only : &
& MOMENTUM_SOLVER

use VERTADV, only : &
& VERTADV_UPDATE

use METHANE_LAKE_MOD, only : &
& METHANE, METHANE_OXIDATION

use SOIL_LAKE_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, DIFFMIN_KPP, &
& LONGWAVE_RADIATION, &
& SHORTWAVE_RADIATION

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_LAKE_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_LAKE_MOD, only : BUBBLE, BUBBLEFLUXAVER

use CONTROL_POINT_LAKE_MOD, only : CONTROL_POINT_OUT, CONTROL_POINT_IN

use TRIBUTARIES, only : TRIBTEMP

use BATHYMSUBR, only : BATHYMETRY

use NUMERICS, only : ACCUMM, VALUE_IS_FINITE

use BL_MOD_LAKE

use SNOWSOIL_MOD
use SNOWTEMP_LAKE, only : SNOWTEMP, SNOW_COND_HEAT_COEF

use RADIATION, only : &
& RadWater, RadIce, RadDeepIce, &
& fracbands, nbands, extwatbands

use ZERODIM_MODEL_LAKE_MOD, only : ZERODIM_MODEL

implicit none

!Input variables
!Reals
real(kind=ireals), intent(in) :: tempair1   !> Air temperature, K
real(kind=ireals), intent(in) :: humair1    !> Specific humidity, kg/kg
real(kind=ireals), intent(in) :: pressure1  !> Air pressure, Pa
real(kind=ireals), intent(in) :: uwind1, vwind1 !> x- and y-components of wind speed, m/s
real(kind=ireals), intent(in) :: longwave1  !> Longwave downward radiation, W/m**2
real(kind=ireals), intent(in) :: netrad1    !> Net radiation, W/m**2
real(kind=ireals), intent(in) :: zref1      !> Height above the lake surface, where the forcing is given, m
real(kind=ireals), intent(in) :: shortwave1 !> Global shortwave radiation, W/m**2
real(kind=ireals), intent(in) :: precip1    !> Precipitation intensity, m/s
real(kind=ireals), intent(in) :: Sflux01    !> Salinity flux at the water surface by rain
real(kind=ireals), intent(in) :: cloud1     !> Cloud sky fraction, 0-1, n/d
real(kind=ireals), intent(in) :: SurfTemp1  !> Surface temperature, K

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

real(kind=ireals), intent(in), target :: 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


real(kind=ireals), intent(in) :: trib_inflow(1:nvartribext) !> The group of tributary data: discharge, ...
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), target :: 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, fdiff_lake_surf1
real(kind=ireals), intent(out) :: h2_out

type(ISIMIP_Lake_type), intent(inout), optional :: LakePhysISIMIP

!------------------------ 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:nveclen) :: a, b, c, d 
real(kind=ireals), allocatable :: 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:nveclen) :: Temp
real(kind=ireals), allocatable :: 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
integer(kind=iintegers) :: nstep_cp_written = 0
integer(kind=iintegers) :: iwork(1:2)

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

!Characters
character(len=60) :: workchar

data uniform_depth /.true./
data nstep_meas /0/
 
SAVE
  

allocate(a(1:nveclen),b(1:nveclen),c(1:nveclen), &
&        d(1:nveclen),Temp(1:nveclen))

!BODY OF PROGRAM

! Getting the number of cores
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)

time_group%year  => year
time_group%month => month
time_group%day   => day
time_group%hour  => hour

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

!Incident global shortwave radiation is computed by empirical formulae, if missing value is provided 
!in the LAKE arguments
if (shortwave1 == missing_value) then
  if (cloud1 /= missing_value) then
    shortwave = SHORTWAVE_RADIATION(time_group,phi,cloud/10.)
  else
    print*, 'Global solar radiaion and cloud fraction are not provided: STOP'
    STOP
  endif
endif

if (longwave1 == missing_value) then
  if (firstcall) write(*,*) 'Atmospheric radiation is missing in input file and will be calculated &
  &from net longwave radiation, net radiation or empirical formulae'
  if (cloud1 == missing_value) then
    if (firstcall) 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
  else
    longwave = LONGWAVE_RADIATION(tempair,humair,pressure,cloud/10.)
  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)
!Control of the input atmospheric forcing
flag = .false.
if (tempair > 60. .or. tempair < -90. .or. (.not. VALUE_IS_FINITE(tempair) )) then
  print*, 'The air temperature ', tempair, 'is illegal: STOP'
  flag = .true.
elseif (abs(humair) > 1. .or. (.not. VALUE_IS_FINITE(humair) )) then
  print*, 'The air humidity ', humair, 'is illegal: STOP'
  flag = .true.
elseif (pressure > 110000 .or. pressure < 20000. .or. (.not. VALUE_IS_FINITE(pressure) )) then
  print*, 'The air pressure ', pressure, 'is illegal: STOP'
  flag = .true.
elseif (abs(uwind) > 200. .or. (.not. VALUE_IS_FINITE(uwind) )) then
  print*, 'The x-component of wind ', uwind, 'is illegal: STOP'
  flag = .true.
elseif (abs(vwind) > 200. .or. (.not. VALUE_IS_FINITE(vwind) )) then
  print*, 'The y-component of wind ', vwind, 'is illegal: STOP'
  flag = .true.
elseif (abs(longwave) > 1000. .or. (.not. VALUE_IS_FINITE(longwave) )) then
  print*, 'The longwave radiation ',longwave,'is illegal: STOP'
  flag = .true.
elseif (abs(shortwave) > 1400. .or. (.not. VALUE_IS_FINITE(shortwave) )) then
  print*, 'The shortwave radiation ', shortwave, 'is illegal: STOP'
  flag = .true.
elseif (abs(precip) > 1. .or. (.not. VALUE_IS_FINITE(precip) )) then
  print*, 'The atmospheric precipitation ',precip, &
  & 'is illegal: STOP'
  flag = .true.
endif
if (flag) then
  write(*,*) 'LAKE: atmospheric forcing is incorrect: 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. &
  & (.not. control_point_read) ) 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'
  control_point_read = .true.
  !call CONTROL_POINT_OUT(year,month,day,hour, &
  !& nx,nx0,nx_max,ny,ny0,ny_max,gs,parparams,time)
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, &
  & 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

!Checking temperature for NaNs and extreme values
do i = 1, M+1
  if (Tw1(i) /= Tw1(i) .or. abs(Tw1(i)) > 100.) then
    print*, 'Tw1(i) = ', Tw1(i), ', nstep = ', nstep, &
    & ', init(ix,iy) = ', init(ix,iy), ', ix = ', ix, ', iy = ', iy
    print*, 'Other coordinatess: ', nx, ny, nx0, ny0, nx_max, ny_max
    print*, 'Terminating LAKE model ...'
    stop
  endif
enddo

!At the next timestep after control point output, the control point may be written again
!Writing the control point for all lakes of the domain
!Note: the lakes state is written for the beginning of the timestep
if (icp == 1_iintegers .and. (nstep /= nstep_cp_written)) then
  !print*, 'Before CP output: ', icp, nstep, nstep_cp_written, rank_comm3d, ix, iy
  call CONTROL_POINT_OUT(year,month,day,hour, &
  & nx,nx0,nx_max,ny,ny0,ny_max,gs,parparams,time)
  nstep_cp_written = nstep
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 (may be overwritten later in BATHYMETRY and/or in TRIBTEMP)
if (trib_inflow(1) /= -9999.) then
  dhwtrib = trib_inflow(1)*dt/area_int(1)
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 .and. .not. (tribheat%par == 4 .and. trib_inflow(1) == -9999.) ) then
  call TRIBTEMP(time,dt,h1,h10,dhwtrib,trib_inflow, &
  &             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 = min(max(dt*(flux1 - flux2)/(row0_m_Lwi),-ddz(M)*h1),ddzi(1)*ls1*roi_d_row0)
!  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