out_mod.f90

MODULE OUT_MOD

use LAKE_DATATYPES, only : ireals, iintegers
use INOUT, only : CHECK_UNIT
use DRIVING_PARAMS, only : missing_value
use GRID_OPERATIONS_LAKE, only : LININTERPOL
  use TIME_LAKE_MOD, only : DATEMINUS, TIMESTR

character, save :: outpath*60 

contains

SUBROUTINE MON_OUT(ix,iy,nx,ny,year,month,day,hour,time,z_out,nout)

use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB, only : E1, E2, eps1, S, Gen, Gen_seiches, row, TKE_turb_trans, TF, KT
use ARRAYS, only : u1, v1
use ARRAYS_GRID, only : dzeta_int, z_full, z_half
use ARRAYS_BATHYM, only : h1, l1, voldef, vol
use DRIVING_PARAMS, only : M
use PHYS_CONSTANTS, only : cw_m_row0

use INOUT_PARAMETERS, only : &
& lake_mon_out_unit_min, &
& lake_mon_out_unit_max

implicit none

! Input variables
! Reals
real(kind=ireals), intent(in) :: hour
real(kind=ireals), intent(in) :: time
real(kind=ireals), intent(in) :: z_out(1:nout) ! The output grid

! Integers
integer(kind=iintegers), intent(in) :: year, month, day
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: nout

! Local variables
! Reals
real(kind=ireals), allocatable :: accum_var (:,:,:,:)
real(kind=ireals), allocatable :: var       (:,:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar (:,:,:)
real(kind=ireals), allocatable :: var_scalar       (:,:,:)
real(kind=ireals), allocatable :: tsteps    (:,:)
real(kind=ireals), allocatable :: Profile_out(:,:)
real(kind=ireals), allocatable :: work(:)
real(kind=ireals) :: voldef0, voldef0y
real(kind=ireals) :: work1, work2
real(kind=ireals), allocatable :: valmax(:)

real(kind=ireals), external:: DZETA

! Integers
integer(kind=iintegers), parameter :: n_var = 11, n_var_scalar = 5, n_var_max = 3
integer(kind=iintegers), parameter :: n_var1 = 4 ! Number of vector variables defined at cell interfaces
integer(kind=iintegers), allocatable :: month_old(:,:)
integer(kind=iintegers), allocatable :: out_unita(:,:)
integer(kind=iintegers) :: out_unit = lake_mon_out_unit_min
integer(kind=iintegers) :: i, j ! Loop indices

! Characters
character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*6
character :: coords_point*6
character :: formatline*50

! Logicals
logical :: firstcall, flag
logical, allocatable :: firstcallixiy(:,:)

data firstcall /.true./

SAVE

if (firstcall) then
  allocate (firstcallixiy(1:nx,1:ny))
  allocate (out_unita(1:nx,1:ny))
  allocate (month_old(1:nx, 1:ny))
  allocate (tsteps    (1:nx, 1:ny))  
  allocate (var      (1:n_var, 1:nx, 1:ny, 1:M+1) )
  allocate (accum_var(1:n_var, 1:nx, 1:ny, 1:M+1) )
  allocate (var_scalar      (1:n_var_scalar, 1:nx, 1:ny) )
  allocate (accum_var_scalar(1:n_var_scalar, 1:nx, 1:ny) )
  allocate (valmax(1:n_var_max))
  allocate (work(1:M+1))
  firstcallixiy(:,:) = .true.
  out_unita(:,:) = lake_mon_out_unit_min
  month_old(:,:) = month
  tsteps   (:,:) = 0.d0  
  accum_var(:,:,:,:) = 0.d0
  accum_var_scalar(:,:,:) = 0.d0
  voldef0 = 0.
  voldef0y = 0.
  valmax(:) = 0.
endif

if (firstcallixiy(ix,iy)) then
  call CHECK_UNIT(lake_mon_out_unit_min,lake_mon_out_unit_max,out_unita(ix,iy))
  write (coords_point,'(2i3)') ix, iy
  open(out_unita(ix,iy),file=outpath(1:len_trim(outpath))//'monthly/'// &
  & 'depvol'//coords_point//'.dat', status='unknown')
  write (out_unita(ix,iy),*) '1 - year' 
  write (out_unita(ix,iy),*) '2 - month'
  write (out_unita(ix,iy),*) '3 - depth, m'
  write (out_unita(ix,iy),*) '4 - ice thickness, m' 
  write (out_unita(ix,iy),*) '5 - volume, m**3'
  write (out_unita(ix,iy),*) '6 - volume deficit, m**3'
  write (out_unita(ix,iy),*) '7 - accumulated volume deficit, m**3'
  write (out_unita(ix,iy),*) '8 - maximal monthly depth, m'
  write (out_unita(ix,iy),*) '9 - maximal monthly ice thickness, m'
  write (out_unita(ix,iy),*) '10 - maximal monthly volume, m**3'
endif

! Variables located at cell interfaces
var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M+1) = row (1:M+1)
var(4,ix,iy,1:M+1) = sqrt(u1 (1:M+1)*u1 (1:M+1) + v1 (1:M+1)*v1 (1:M+1) )
! Variables located at cell centers
var(5,ix,iy,1:M)   = E1  (1:M)
var(6,ix,iy,1:M)   = TF(1:M)*E2(1:M) !eps1(1:M)
var(7,ix,iy,1:M)   = S   (1:M)
var(8,ix,iy,1:M)   = Gen (1:M)
var(9,ix,iy,1:M)   = Gen_seiches (1:M)
var(10,ix,iy,1:M)  = TKE_turb_trans (1:M)
var(11,ix,iy,1:M)  = KT (1:M)/cw_m_row0 !Kinematic heat conductance, m**2/s

var_scalar(1,ix,iy) = h1
var_scalar(2,ix,iy) = l1
var_scalar(3,ix,iy) = vol
var_scalar(4,ix,iy) = voldef - voldef0
var_scalar(5,ix,iy) = voldef - voldef0
voldef0 = voldef

valmax(1) = max(h1,valmax(1)) ! Calculating monthly maximum
valmax(2) = max(l1,valmax(2)) 
valmax(3) = max(vol,valmax(3))

accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + &
& var_scalar(:,ix,iy)
tsteps(ix,iy) = tsteps(ix,iy) + 1._ireals

if (month_old(ix,iy) /= month) then
  !print*, time/(365./12.*24.*3600.)
  !read*
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
  accum_var_scalar(1:3,ix,iy) = accum_var_scalar(1:3,ix,iy)/tsteps(ix,iy) ! excepting volume deficit

  call DATEMINUS(1,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR(6,year1,month1,day1,hour1,timestring)

  call CHECK_UNIT(lake_mon_out_unit_min,lake_mon_out_unit_max,out_unit)
  write (coords_point,'(2i3)') ix, iy
  open(out_unit,file=outpath(1:len_trim(outpath))//'monthly/'// &
  & 'Profiles'//coords_point//timestring//'.dat', status='unknown')
  write (out_unit,*) '1 - depth, m' 
  write (out_unit,*) '2 - temperature, C'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - density, kg/m**3'
  write (out_unit,*) '5 - absolute value of velocity, m/s'
  write (out_unit,*) '6 - turbulent kinetic energy, m**2/s**2'
  write (out_unit,*) '7 - disspation rate, m**2/s**3'
  write (out_unit,*) '8 - TKE production by buoyancy, m**2/s**3'
  write (out_unit,*) '9 - TKE production by shear   , m**2/s**3'
  write (out_unit,*) '10 - TKE production by seiches (Goudsmit paramaterization)  , m**2/s**3'
  write (out_unit,*) '11 - TKE turbulent transport  , m**2/s**3'
  write (out_unit,*) '12 - kinematic heat conductance, m**2/s'
  formatline = '(f12.6,f9.3,e12.4,f9.3,8e12.4)'
  if (nout > 0) then
    allocate (Profile_out(1:nout,1:n_var))
    ! Interpolating to given output levels
    do j = 1, n_var1
      work(1:M+1) = accum_var(j,ix,iy,1:M+1)
      call LININTERPOL (z_full,work,M+1,z_out,Profile_out(1,j),nout,flag)
    enddo
    do j = n_var1+1, n_var
      work(1:M) = accum_var(j,ix,iy,1:M)
      call LININTERPOL (z_half,work,M,z_out,Profile_out(1,j),nout,flag)
    enddo
    j = nout
    if (z_out(nout) > z_full(M+1)) then
      do while (z_out(j) > z_full(M+1))
        j = j - 1
      enddo
    endif
    ! Writing to file
    do i = 1, j 
      write (unit = out_unit, fmt = formatline) &
      & -z_out(i), Profile_out(i,1:n_var)
    enddo
  else
    allocate (Profile_out(1:M+1,1:n_var))
    forall (i = 1:M+1, j = 1:n_var1) Profile_out(i,j) = accum_var(j,ix,iy,i)
    ! Interpolating from cell centers to cell interfaces
    do j = n_var1+1, n_var
      work(1:M) = accum_var(j,ix,iy,1:M)
      call LININTERPOL (z_half,work,M,z_full,Profile_out(1,j),M+1,flag)
    enddo
    ! Writing to file
    do i = 1, M+1 
      write (unit = out_unit, fmt = formatline) &
      & -dzeta_int(i)*h1, Profile_out(i,1:n_var)
    enddo
  endif
  deallocate(Profile_out)
  close(out_unit)


  read(year1,*) work1
  read(month1,*) work2
  write (out_unita(ix,iy),fmt=*) & !, fmt = formatline) &
  & int(work1), int(work2), accum_var_scalar(1:n_var_scalar,ix,iy), &
  & valmax(:)

  
  tsteps(ix,iy) = 0.d0
  accum_var(:,ix,iy,:) = 0.d0
  accum_var_scalar(1:4,ix,iy) = 0.d0
  if (month_old(ix,iy) == 12) accum_var_scalar(5,ix,iy) = 0.d0
  valmax(:) = 0.d0

  month_old(ix,iy) = month
endif

if (firstcall) firstcall = .false.
if (firstcallixiy(ix,iy)) firstcallixiy(ix,iy) = .false.
END SUBROUTINE MON_OUT


SUBROUTINE DAY_OUT(ix,iy,nx,ny,year,month,day,hour)

use ARRAYS_GRID, only : dzeta_int
use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB, only : E1, eps1, KT, Ri
use ARRAYS_BATHYM, only : h1
use ARRAYS, only : u1, v1
use ATMOS , only : discharge
use DRIVING_PARAMS, only : M
use PHYS_CONSTANTS, only : cw_m_row0

use INOUT_PARAMETERS, only : &
& lake_day_out_unit_min, &
& lake_day_out_unit_max

implicit none

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

! Integers
integer(kind=iintegers), intent(in) :: year, month, day
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny

! Local variables
! Reals
real(kind=ireals), allocatable :: tsteps   (:,:)
real(kind=ireals), allocatable :: accum_var(:,:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar (:,:,:)
real(kind=ireals), allocatable :: var_scalar (:,:,:)
real(kind=ireals), allocatable :: var      (:,:,:,:)

real(kind=ireals), external :: DZETA

real(kind=ireals) :: work1, work2, work3

! Integers
integer(kind=iintegers), allocatable :: day_old(:,:)
integer(kind=iintegers), allocatable :: out_unita(:,:)
integer(kind=iintegers) :: i ! Loop index
integer(kind=iintegers) :: out_unit = lake_day_out_unit_min
integer(kind=iintegers) :: nvars = 7
integer(kind=iintegers) :: nvars_scalar = 2

character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*8
character :: coords_point*6

logical :: firstcall
logical, allocatable :: firstcallixiy(:,:)

data firstcall /.true./
      
SAVE

if (firstcall) then
  allocate (firstcallixiy(1:nx,1:ny))
  allocate (out_unita(1:nx,1:ny))
  allocate (day_old   (1:nx, 1:ny))
  allocate (tsteps    (1:nx, 1:ny))  
  allocate (var       (1:nvars, 1:nx, 1:ny, 1:M+1) )
  allocate (var_scalar(1:nvars_scalar, 1:nx, 1:ny ) )
  allocate (accum_var (1:nvars, 1:nx, 1:ny, 1:M+1) )
  allocate (accum_var_scalar (1:nvars_scalar, 1:nx, 1:ny) )
  firstcallixiy(:,:) = .true.
  out_unita(:,:) = lake_day_out_unit_min
  day_old  (:,:) = day
  tsteps   (:,:) = 0.d0  
  accum_var(:,:,:,:) = 0.d0
  accum_var_scalar(:,:,:) = 0.d0
  var       (:,:,:,:) = 0.d0
  var_scalar(:,:,:) = 0.d0
endif

if (firstcallixiy(ix,iy)) then
  call CHECK_UNIT(lake_day_out_unit_min,lake_day_out_unit_max,out_unita(ix,iy))
  write (coords_point,'(2i3)') ix, iy
  open(out_unita(ix,iy),file=outpath(1:len_trim(outpath))//'daily/'// &
  & 'fluxes'//coords_point//'.dat', status='unknown')
  write (out_unita(ix,iy),*) '1 - year' 
  write (out_unita(ix,iy),*) '2 - month'
  write (out_unita(ix,iy),*) '3 - discharge in x direction, m**3/s'
  write (out_unita(ix,iy),*) '4 - discharge in y direction, m**3/s' 
endif

!Vector variables
var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M)   = E1  (1:M)
var(4,ix,iy,1:M)   = eps1(1:M)
var(5,ix,iy,1:M)   = KT(1:M)/cw_m_row0
var(6,ix,iy,1:M+1) = sqrt(u1(1:M+1)**2 + v1(1:M+1)**2)
var(7,ix,iy,1:M)   = Ri(1:M)

!Scalar variables
var_scalar(1,ix,iy) = discharge(1)
var_scalar(2,ix,iy) = discharge(2)

accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
tsteps(ix,iy) = tsteps(ix,iy) + 1.d0

if (day_old(ix,iy)/=day) then
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
  call DATEMINUS (2,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR (8,year1,month1,day1,hour1,timestring)
  call CHECK_UNIT(lake_day_out_unit_min,lake_day_out_unit_max,out_unit)
  write (coords_point, '(2i3)') ix, iy
  open(out_unit, file=outpath(1:len_trim(outpath))//'daily/'// &
  & 'Profiles'//coords_point//timestring//'.dat',status='unknown')
  write (out_unit,*) '1 - depth, m' 
  write (out_unit,*) '2 - temperature, C'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - turbulent kinetic energy, m**2/s**2'
  write (out_unit,*) '5 - dissipation rate, m**2/s**3'
  write (out_unit,*) '6 - kinematic heat conductance, m**2/s'
  write (out_unit,*) '7 - modulus of horizontal velocity, m/s'
  write (out_unit,*) '8 - Richardson number, n/d'
  do i = 1, M+1 
    write (out_unit,'(f12.6,f11.5,6e12.4)') &
    &  -dzeta_int(i)*h1, accum_var(1:nvars,ix,iy,i)
  enddo
  close(out_unit)

  accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy)/tsteps(ix,iy)
  read(year1,*) work1
  read(month1,*) work2
  read(day1,*) work3
  write (out_unita(ix,iy),fmt=*) & !, fmt = formatline) &
  & int(work1), int(work2), int(work3), accum_var_scalar(1:nvars_scalar,ix,iy)

  day_old(ix,iy) = day
  accum_var(:,ix,iy,:) = 0.d0
  accum_var_scalar(:,ix,iy) = 0.d0
  tsteps(ix,iy) = 0.d0
endif

if (firstcall) firstcall=.false.
if (firstcallixiy(ix,iy)) firstcallixiy(ix,iy) = .false.
END SUBROUTINE DAY_OUT


SUBROUTINE HOUR_OUT(ix,iy,nx,ny,year,month,day,hour,time)

use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB!, only : E1, eps1, row, H_mixed_layer, H_entrainment, Wedderburn
use ARRAYS_BATHYM, only : h1
use ARRAYS_GRID, only : dzeta_int, dzeta_05int
use ARRAYS, only : u1, v1
use TIMEVAR, only : hour_sec
use ARRAYS_OXYGEN, only: oxyg, DIC
use ARRAYS_METHANE, only : qwater
use PHYS_FUNC, only : HC_CORR_CARBEQUIL
use PHYS_CONSTANTS, only : Kelvin0
use METH_OXYG_CONSTANTS, only : &
& molm3tomgl_o2, molm3tomcM, molm3tomcM, pH

use DRIVING_PARAMS, only : M, scale_output

use INOUT_PARAMETERS, only : &
& lake_hour_out_unit_min, &
& lake_hour_out_unit_max

implicit none

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

! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day

! Local variables
! Reals
real(kind=ireals), allocatable :: accum_var (:,:,:,:)
real(kind=ireals), allocatable :: var       (:,:,:,:)
real(kind=ireals), allocatable :: var_scalar(:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar(:,:,:)

real(kind=ireals), external :: DZETA

! Integers
integer(kind=iintegers), parameter :: n_var = 30
integer(kind=iintegers), parameter :: n_var_scalar = 10
integer(kind=iintegers), allocatable :: hour_old(:,:)
integer(kind=iintegers), allocatable :: tsteps(:,:)
integer(kind=iintegers), allocatable :: n_unit(:,:)
integer(kind=iintegers) :: i
integer(kind=iintegers) :: out_unit = lake_hour_out_unit_min

character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*10
character :: coords_point*6
character :: format_char*100

logical :: firstcall = .true.
logical, allocatable :: firstcall_ixiy(:,:)

      
SAVE

if (firstcall) then
  allocate (accum_var (n_var,1:nx,1:ny,1:M+1) )
  allocate (var       (n_var,1:nx,1:ny,1:M+1) )
  allocate (accum_var_scalar (n_var_scalar,1:nx,1:ny) )
  allocate (var_scalar       (n_var_scalar,1:nx,1:ny) )  
  allocate (tsteps    (1:nx,1:ny) )
  allocate (hour_old  (1:nx,1:ny) )
  allocate (n_unit(1:nx,1:ny))
  allocate (firstcall_ixiy(1:nx,1:ny))
  accum_var(:,:,:,:) = 0.d0
  accum_var_scalar(:,:,:) = 0.d0
  var(:,:,:,:) = 0.d0
  var_scalar(:,:,:) = 0.d0
  tsteps(:,:) = 0
  hour_old(:,:) = int(hour)
  firstcall_ixiy(:,:) = .true.
  n_unit(:,:) = out_unit + 1
endif

var(1,ix,iy,1:M)   = PEMF          (1:M)
var(2,ix,iy,1:M+1) = PDENS_DOWN    (1:M+1)
var(3,ix,iy,1:M+1) = PT_DOWN       (1:M+1)
var(4,ix,iy,1:M+1) = PSAL_DOWN     (1:M+1)
var(5,ix,iy,1:M)   = k_turb_T_flux (1:M)
var(6,ix,iy,1:M)   = T_massflux    (1:M)
var(7,ix,iy,1:M+1) = row           (1:M+1)
var(8,ix,iy,1:M+1) = Tw1           (1:M+1)
var(9,ix,iy,1:M+1) = Sal1          (1:M+1)
var(10,ix,iy,1:M)  = E1            (1:M)
var(11,ix,iy,1:M)  = eps1          (1:M)

if (scale_output%par == 1) then
  var(12,ix,iy,1)    = H_mixed_layer ! Scale
  var(13,ix,iy,1)    = Buoyancy0     ! Scale
  var(14,ix,iy,1)    = w_conv_scale  ! Scale
  var(15,ix,iy,1)    = T_conv_scale  ! Scale
elseif (scale_output%par == 0) then
  var(12,ix,iy,1)    = 1.d0 ! No scaling
  var(13,ix,iy,1)    = 1.d0 ! No scaling
  var(14,ix,iy,1)    = 1.d0 ! No scaling
  var(15,ix,iy,1)    = 1.d0 ! No scaling
else
  print*, 'Scale_output is ', scale_output%par, &
  & ', must be 0 or 1: STOP'  
  STOP
endif
      
var(16,ix,iy,1:M)  = S              (1:M)
var(17,ix,iy,1:M)  = Gen            (1:M)
var(18,ix,iy,1:M)  = Gen_seiches    (1:M)
var(19,ix,iy,1:M)  = TKE_turb_trans (1:M)
      
var(20,ix,iy,2:M)  = k3_mid         (2:M)

var(21,ix,iy,1:M+1)  = u1(1:M+1)
var(22,ix,iy,1:M+1)  = v1(1:M+1)
var(23,ix,iy,1:M+1)  = sqrt(u1(1:M+1)*u1(1:M+1) + v1(1:M+1)*v1(1:M+1))

!Biogeochemical substances
var(24,ix,iy,1:M+1)  = oxyg(1:M+1,1)*molm3tomgl_o2
do i = 1, M+1 
  var(25,ix,iy,i)  = DIC(i,1) / HC_CORR_CARBEQUIL(Tw1(i)+Kelvin0,pH) * molm3tomcM
enddo
var(26,ix,iy,1:M+1)  = qwater(1:M+1,1)*molm3tomcM

!Coefficients of turbulent viscosity
var(27,ix,iy,1:M)   = k2             (1:M  )
var(28,ix,iy,1:M)   = viscturb_ziriv (1:M  )
var(29,ix,iy,1:M)   = eps_ziriv      (1:M  )
var(30,ix,iy,1:M+1) = speed_ziriv    (1:M+1)

var_scalar(1,ix,iy) = S_integr_positive
var_scalar(2,ix,iy) = S_integr_negative
var_scalar(3,ix,iy) = Gen_integr
var_scalar(4,ix,iy) = eps_integr
var_scalar(5,ix,iy) = TKE_balance
var_scalar(6,ix,iy) = H_entrainment
var_scalar(7,ix,iy) = TKE_turb_trans_integr
var_scalar(8,ix,iy) = Wedderburn
var_scalar(9,ix,iy) = LakeNumber
var_scalar(10,ix,iy) = Rossby_rad

accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
tsteps(ix,iy) = tsteps(ix,iy) + 1

if (int(hour)/=hour_old(ix,iy)) then
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) / real(tsteps(ix,iy))
  accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) / real(tsteps(ix,iy))
  call DATEMINUS(3,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR(10,year1,month1,day1,hour1,timestring)
  call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,out_unit)
  write (coords_point, '(2i3)') ix, iy
  
! Writing to the file Profiles<yyyymmddhh>.dat  
  open(out_unit, file=outpath(1:len_trim(outpath))// &
  & 'hourly/'//'Profiles'//coords_point//timestring//'.dat', &
  &  status='unknown')
  write (out_unit,*) '1 - depth, normalized' 
  write (out_unit,*) '2 - temperature, normalized'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - water density, kg/m**3'  
  write (out_unit,*) '5 - turbulent kinetic energy, normalized'
  write (out_unit,*) '6 - dissipation rate, normalized'
  write (out_unit,*) '7 - eddy viscosity for TKE, m**2/s'
  write (out_unit,*) '8  - mass flux, m/s'
  write (out_unit,*) '9  - downdraft temperature, C' 
  write (out_unit,*) '10 - downdraft salinity, kg/kg' 
  write (out_unit,*) '11 - downdraft density, kg/m**3'
  write (out_unit,*) '12 - x-component of speed, m/s'
  write (out_unit,*) '13 - y-component of speed, m/s'
  write (out_unit,*) '14 - modulus of speed, m/s'
  write (out_unit,*) '15 - oxygen, mg/l'
  write (out_unit,*) '16 - co2, mcmol/l'
  write (out_unit,*) '17 - ch4, mcmol/l'
  write (out_unit,*) '18 - eddy viscosity from k-epsilon model, m**2/s'
  write (out_unit,*) '19 - eddy viscosity from Ziriv model, m**2/s'
  write (out_unit,*) '20 - dissipation rate from Ziriv model, m**2/s**3'
  write (out_unit,*) '21 - longitudinal speed from Ziriv model, m/s'
  do i=1,M 
    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e16.7,3f11.5,3e16.7,3e10.3,f11.5)')      &
    & -dzeta_int(i)*h1                         /accum_var(12,ix,iy,1),  &
    & (accum_var(8,ix,iy,i) - scale_output%par*maxval(accum_var(8,ix,iy,:)) ) &
    &  /accum_var(15,ix,iy,1), &
    & accum_var(9,ix,iy,i),                                               &
    & accum_var(7,ix,iy,i),                                               &
    & accum_var(10,ix,iy,i)    / (accum_var(13,ix,iy,1)*accum_var(12,ix,iy,1) /       &
    &  accum_var(14,ix,iy,1) ),                                           &
    & accum_var(11,ix,iy,i)    / accum_var(13,ix,iy,1),                   &
    & accum_var(20,ix,iy,i),                                              & !accum_var(10,i)**2 / accum_var(11,i),
    & accum_var(1,ix,iy,i),                                               &
    & accum_var(3,ix,iy,i),                                               &
    & accum_var(4,ix,iy,i),                                               &
    & accum_var(2,ix,iy,i),                                               &
    & accum_var(21,ix,iy,i),                                              &
    & accum_var(22,ix,iy,i),                                              &
    & accum_var(23,ix,iy,i),                                              &
    & accum_var(24,ix,iy,i),                                              &
    & accum_var(25,ix,iy,i),                                              &
    & accum_var(26,ix,iy,i),                                              &
    & accum_var(27,ix,iy,i),                                              &
    & accum_var(28,ix,iy,i),                                              &
    & accum_var(29,ix,iy,i),                                              &
    & accum_var(30,ix,iy,i)
  enddo 
  close (out_unit)
  
! Writing to the file EDMF_profiles<yyyymmddhh>.dat
  open(out_unit, file=outpath(1:len_trim(outpath))//  &
  & 'hourly/'//'EDMF_profiles'//coords_point//timestring//'.dat',   &
  & status='unknown')
  write (out_unit,*) '1 - depth, normalized by mixed layer depth' 
  write (out_unit,*) '2 - "k - flux" of temperature, normalized'
  write (out_unit,*) '3 - mass flux of temperature,  normalized'
  write (out_unit,*) '4 - total flux of temperature, normalized'
  write (out_unit,'( 4(i10,6x) )') 1,2,3,4
  do i=1, M
    write (out_unit,'(f18.6,3e16.7)')                   &
    & -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1),          &
    & accum_var(5,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
    & accum_var(6,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
    & (accum_var(5,ix,iy,i) + accum_var(6,ix,iy,i) ) /              &
    & (accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1))
  enddo
  close (out_unit)
  
! Writing to the file TKE_budget<yyyymmddhh>.dat
  open(out_unit, file=outpath(1:len_trim(outpath))//    &
  & 'hourly/'//'TKE_budget'//coords_point//timestring//'.dat',        &
  &  status='unknown')
  write (out_unit,*) '1 - depth, normalised with mixed layer depth'
  write (out_unit,*) '2 - shear production, normalised with surface buoyancy flux'
  write (out_unit,*) '3 - seiches production, normalised with surface buoyancy flux'
  write (out_unit,*) '4 - buoyancy source,  normalised with surface buoyancy flux'
  write (out_unit,*) '5 - dissipation rate, normalised with surface buoyancy flux'
  write (out_unit,*) '6 - turbulent transport, normalised with surface buoyancy flux'
  write (out_unit,'( 6(i10,6x) )') 1,2,3,4,5,6
  do i = 1, M
    write (out_unit,'(f18.6,5e16.7)')          &
    & -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1), &
    & accum_var(17,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(18,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(16,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(11,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(19,ix,iy,i) / accum_var(13,ix,iy,1)
  enddo
!  write (out_unit, '(a)')
  close (out_unit)
  
  if (firstcall_ixiy(ix,iy)) then
    call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,n_unit(ix,iy))
    open(n_unit(ix,iy), file=outpath(1:len_trim(outpath))//           &
    & 'hourly/'//'TKE_integr_conv'//coords_point//timestring//'.dat', &
    &  status='unknown')
    write (n_unit(ix,iy),*) 'Col. 1 - year'
    write (n_unit(ix,iy),*) 'Col. 2 - month'
    write (n_unit(ix,iy),*) 'Col. 3 - day'
    write (n_unit(ix,iy),*) 'Col. 4 - hour'
    write (n_unit(ix,iy),*) 'Col. 5 - the time from the start of integration, hours'
    write (n_unit(ix,iy),*) 'Col. 6 - H_entrainment, m    '
    write (n_unit(ix,iy),*) 'Col. 7 - B_integr+,   m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 8 - B_integr-,   m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 9 - S_integr,    m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 10 - eps_integr,  m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 11 - TKE turb trans integr,  m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 12 - TKE_balance, m**3/s**3'
    write (n_unit(ix,iy),*) 'Col. 13 - Wedderburn number, n/d'
    write (n_unit(ix,iy),*) 'Col. 14 - Lake number, n/d'
    write (n_unit(ix,iy),*) 'Col. 15 - Internal Rossby radius, m'
    firstcall_ixiy(ix,iy) = .false.
  endif
  
  format_char = '(3i7,f8.2,f13.2,f10.2,9e15.6)'
  write (n_unit(ix,iy), format_char) &
  &                    year,month,day,hour,time/hour_sec, &
  &                    accum_var_scalar(6,ix,iy), & 
  &                    accum_var_scalar(1,ix,iy), &
  &                    accum_var_scalar(2,ix,iy), &
  &                    accum_var_scalar(3,ix,iy), &
  &                    accum_var_scalar(4,ix,iy), &
  &                    accum_var_scalar(7,ix,iy), &
  &                    accum_var_scalar(5,ix,iy), &
  &                    accum_var_scalar(8,ix,iy), &
  &                    accum_var_scalar(9,ix,iy), &
  &                    accum_var_scalar(10,ix,iy)
  
  hour_old(ix,iy) = int(hour)
  accum_var(:,ix,iy,:) = 0.d0
  accum_var_scalar(:,ix,iy) = 0.d0
  tsteps(ix,iy)  = 0
   
endif

if (firstcall) firstcall = .false.
!if (firstcall_ixiy(ix,iy)) firstcall_ixiy(ix,iy) = .false.
END SUBROUTINE HOUR_OUT          
      


SUBROUTINE EVERYSTEP_OUT(ix,iy,nx,ny)

use ARRAYS!, only : time, nstep

use ARRAYS_WATERSTATE!, only : Tw1, Sal1 
use ARRAYS_BATHYM!, only : h1
use ARRAYS_GRID!, only : dzeta_int, M

use ARRAYS_TURB!, only : &
! & S_integr_positive, &
! & S_integr_negative, &
! & Gen_integr, &
! & eps_integr, &
! & TKE_balance, &
! & E_integr, &
! & H_entrainment, &
! & E1, eps1, &
! & k_turb_T_flux

use DRIVING_PARAMS, only : everystep, M

use INOUT_PARAMETERS, only : &
& lake_everystep_out_unit_min, &
& lake_everystep_out_unit_max

implicit none

! Input variables
! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny

! Local variables
! Reals
real(kind=ireals), external :: DZETA
real(kind=ireals), parameter :: ACC = 1.d-20
real(kind=ireals), parameter :: hour_sec = 60.*60.

! Integers
integer(kind=iintegers), allocatable :: n_unit(:,:,:)
integer(kind=iintegers) :: i ! Loop index

! Characters
character :: coords_point*6
character :: format_char*100

! Logicals
logical, allocatable :: firstcall(:,:)

      
SAVE

if (.not.allocated(firstcall)) then
  allocate (firstcall(1:nx, 1:ny) )
  allocate (n_unit   (1:2, 1:nx, 1:ny) )
  firstcall = .true.
  n_unit = lake_everystep_out_unit_min
endif  
      
if (firstcall(ix,iy)) then
  write (coords_point, '(2i3)') ix, iy
  
  if (everystep%par /= 2) then
    call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
    & n_unit(1,ix,iy))
    open (n_unit(1,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
    & 'Profiles'//coords_point//'.dat',  status='unknown')
    write (n_unit(1,ix,iy),*) '1 - depth, m' 
    write (n_unit(1,ix,iy),*) '2 - temperature, C'
    write (n_unit(1,ix,iy),*) '3 - salinity, kg/kg' 
    write (n_unit(1,ix,iy),*) '4 - turbulent kinetic energy, m**2/s**2'
    write (n_unit(1,ix,iy),*) '5 - disspation rate, m**2/s**3'
    write (n_unit(1,ix,iy),*) '6 - eddy diffusivity (TKE**2/dissipation), m**2/s'
    write (n_unit(1,ix,iy),*) '7 - k-flux of temperature, K*m/s'
  endif
    
  call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
  & n_unit(2,ix,iy))
  open (n_unit(2,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
  & 'TKE_integr'//coords_point//'.dat',  status='unknown')
  write (n_unit(2,ix,iy),*) '1 - timestep              '
  write (n_unit(2,ix,iy),*) '2 - time,        hours    '
  write (n_unit(2,ix,iy),*) '3 - B_integr+,   m**3/s**3'
  write (n_unit(2,ix,iy),*) '4 - B_integr-,   m**3/s**3'
  write (n_unit(2,ix,iy),*) '5 - S_integr,    m**3/s**3'
  write (n_unit(2,ix,iy),*) '6 - eps_integr,  m**3/s**3'
  write (n_unit(2,ix,iy),*) '7 - TKE_balance, m**3/s**3'
  write (n_unit(2,ix,iy),*) '8 - E_integr,    m**3/s**2'
endif

if (everystep%par /= 2) then      
  write (n_unit(1,ix,iy),*) 'nstep = ', nstep
  format_char = '(f10.6,f9.3,5e12.4)'
  do i=1,M 
    write (n_unit(1,ix,iy), format_char) &
    & -dzeta_int(i)*h1,Tw1(i),Sal1(i),E1(i),eps1(i),  &
    & E1(i)**2/(eps1(i)+ACC), k_turb_T_flux(i)
  enddo
endif

format_char = '(i10, 2f10.2, 6e15.6)'
write(n_unit(2,ix,iy), format_char) &
& nstep, time/hour_sec, H_entrainment, S_integr_positive, S_integr_negative, &
& Gen_integr, eps_integr, TKE_balance, E_integr
      
if (firstcall(ix,iy)) firstcall(ix,iy)=.false.
END SUBROUTINE EVERYSTEP_OUT
      
      
SUBROUTINE SERIES_OUT(ix,iy,nx,ny,year,month,day,hour,tsw)

use DRIVING_PARAMS, only : M, ns, Mice, dt_out
use NUMERIC_PARAMS , only : &
& ms, small_value

use ARRAYS_TURB, only : row, H_mixed_layer, H_entrainment, &
& signwaveheight, w_conv_scale, u_star, Ri_bulk, maxN, i_maxN, &
& ThermThick, ReTherm, RiTherm, TKE_budget_terms

use ARRAYS_WATERSTATE, only : Tw1, Ti1, lamw, salice

use ARRAYS_GRID, only : nsoilcols

use ARRAYS_BATHYM, only : &
& h1,hx1,hx2,hy1,hy2, &
& hx1t,hx2t,hy1t,hy2t, &
& l1,hs1,ls1,  &
& voldef, vol,     &
& area_int, bathymwater

use ARRAYS_METHANE, only : &
& fbbleflx_ch4_sum, &
& fbbleflx_ch4,    &
& febul0, fplant,   &
& fdiff_lake_surf, &
& fdiffbot, foutflow, &
& qwater, qsoil, &
& rprod_total_newC, &
& rprod_total_oldC, &
& h_talik, lammeth, &
& qwateroxidtot, qwateroxidML

use ARRAYS_OXYGEN, only : &
& fbbleflx_co2_sum, &
& fbbleflx_co2,    &
& fbbleflx_o2_sum, &
& fbbleflx_o2, &
& oxyg, &
& fco2, fo2

use ARRAYS_SOIL, only : Tsoil1, lsm

use ARRAYS, only:  &
& Tskin,       &
& zinv,            &
& time,            &
& T_0dim,           &
& roughness,emissivity,albedo,aM,bM,relhums

use ATMOS,  only:    &
& eflux0_kinem,      &
& eflux0,            &
& turb_density_flux, &
& hw, xlew, botflux, &
& tau, windwork, surfrad, &
& discharge, velfrict_bot, &
& shortwave

use INOUT_PARAMETERS, only : &
& lake_series_out_unit_min, &
& lake_series_out_unit_max

use METH_OXYG_CONSTANTS, only : &
& molmass_ch4, mfs

use TIMEVAR, only : &
hour_sec, day_sec, hour_min

use PHYS_CONSTANTS, only : &
Kelvin0, cw_m_row0, roa0, sabs, row0

use DZETA_MOD, only : VARMEAN

use ARRAYS_GRID, only : gsp

use SNOWSOIL_MOD, only : Tsn, itop

implicit none

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

! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day

! Local variables
! Reals

!real(kind=ireals), parameter :: mfs = molmass_ch4*8.64d+7 ! the transform multiplier for methane flux
                                                          ! from mol/(m**2*s) to mg/(m**2*day)
real(kind=ireals), parameter :: small_number = 1.d-5

real(kind=ireals) :: T_mean, voldef0 = 0.
real(kind=ireals) :: methfluxML, methfluxshalsed

real(kind=ireals), allocatable :: tsteps(:,:)
real(kind=ireals), allocatable :: var_scalar (:,:,:) 
real(kind=ireals), allocatable :: accum_var_scalar(:,:,:)

! Integers
integer(kind=iintegers), parameter :: nfiles = 6
integer(kind=iintegers), parameter :: numaccum = 5 ! The number of scalar variables being accumulated (averaged)
integer(kind=iintegers), allocatable :: n_unit(:,:,:)
integer(kind=iintegers), allocatable :: count_out(:,:)

! Characters
character :: coords_point*6
character :: format_char*100, workchar*10

! Logicals
logical, allocatable :: firstcall(:,:)


!real(kind=ireals) :: dz(ms)
integer(kind=iintegers) :: i
!common /SOILDAT/ dz,itop
!real(kind=ireals) :: Tsn(1:ms), cs(1:ms)
!common /SNOW_CHAR/ Tsn,cs

      
SAVE

if (.not.allocated(firstcall)) then
  allocate (firstcall(1:nx, 1:ny) )
  allocate (count_out(1:nx, 1:ny) )
  allocate (n_unit   (1:nfiles, 1:nx, 1:ny) )
  allocate (tsteps           (1:nx, 1:ny))  
  allocate (var_scalar       (1:numaccum, 1:nx, 1:ny) )
  allocate (accum_var_scalar (1:numaccum, 1:nx, 1:ny) )

  tsteps(:,:) = 0.d0
  var_scalar(:,:,:) = 0.d0
  accum_var_scalar(:,:,:) = 0.d0

  firstcall(:,:) = .true.
  n_unit = lake_series_out_unit_min
endif  
            
if (firstcall(ix,iy)) then
  write (coords_point, '(2i3)') ix, iy

  i = 1
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(i,ix,iy))
  open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'layers'//coords_point//'.dat',  status='unknown')
  write (n_unit(i,ix,iy),*)'Col. 1 - year'
  write (n_unit(i,ix,iy),*)'Col. 2 - month'
  write (n_unit(i,ix,iy),*)'Col. 3 - day'
  write (n_unit(i,ix,iy),*)'Col. 4 - hour'
  write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(i,ix,iy),*)'Col. 6 - water layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 7 - W mixed layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 8 - E mixed layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 9 - S mixed layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 10 - N mixed layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 11 - W lower layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 12 - E lower layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 13 - S lower layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 14 - N lower layer thickness, m'
  write (n_unit(i,ix,iy),*)'Col. 15 - ice layer thickness,   m'
  write (n_unit(i,ix,iy),*)'Col. 16 - snow layer thickness,  m'
  write (n_unit(i,ix,iy),*)'Col. 17 - bottom ice thickness,  m'
  write (n_unit(i,ix,iy),*)'Col. 18 - reservoir volume,  m**3'
  write (n_unit(i,ix,iy),*)'Col. 19 - volume deficit (accumulated),  m**3'
    
  i = i + 1
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(i,ix,iy))
  open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'T_fluxes'//coords_point//'.dat',status='unknown')
  write (n_unit(i,ix,iy),*)'Col. 1 - year'
  write (n_unit(i,ix,iy),*)'Col. 2 - month'
  write (n_unit(i,ix,iy),*)'Col. 3 - day'
  write (n_unit(i,ix,iy),*)'Col. 4 - hour'
  write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(i,ix,iy),*)'Col. 6 - surface temperature, C'
  write (n_unit(i,ix,iy),*)'Col. 7 - water skin temperature, C'
  write (n_unit(i,ix,iy),*)'Col. 8 - water surface temperature, C'
  write (n_unit(i,ix,iy),*)'Col. 9 - mean temperature of water coloumn, C'
  write (n_unit(i,ix,iy),*)'Col. 10 - maximal temperature in the water coloumn, C'
  write (n_unit(i,ix,iy),*)'Col. 11 - zero-dimensional model temperature, C'  
  write (n_unit(i,ix,iy),*)'Col. 12 - upper ice surface temperature, C'
  write (n_unit(i,ix,iy),*)'Col. 13 - upper snow surface temperature, C'
  write (n_unit(i,ix,iy),*)'Col. 14 - sensible heat flux,    W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 15 - latent heat flux,      W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 16 - downward heat flux at the upper lake surface, W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 17 - downward heat flux at the lake bottom, W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 18 - friction velocity at the surface (waterside), m/s'
  write (n_unit(i,ix,iy),*)'Col. 19 - friction velocity at the bottom, m/s'
  write (n_unit(i,ix,iy),*)'Col. 20 - wind work at the water surface, W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 21 - albedo of the lake-atmosphere interface, n/d'
  write (n_unit(i,ix,iy),*)'Col. 22 - shortwave radiation penetrated below surface, W/m**2'
  write (n_unit(i,ix,iy),*)'Col. 23 - significant wave height, m'
  write (n_unit(i,ix,iy),*)'Col. 24 - bottom ice salinity, kg/kg'
  write (n_unit(i,ix,iy),*)'Col. 25 - discharge in x direction, m**3/s'
  write (n_unit(i,ix,iy),*)'Col. 26 - discharge in y direction, m**3/s'
  write (n_unit(i,ix,iy),*)'Cols. 27... - top soil columns temperature, m'