PROGRAM PROFILES_PROCESSING

implicit none

real(8), parameter :: day_sec = 24.*60.*60.

integer(4), parameter :: nfiles = 1 ! Number of input files with profiles
integer(4), parameter :: nlevels = 20 ! Number of levels in profiles

! nlevels = 21 for Sparkling, except LAKEoneD

integer(4) :: year, month, month_old, year_old, day
integer(4) :: year_term, month_term
integer(4) :: i, j, k
integer(4), external :: MONTH_DAYS

real(8) :: hour, month_sec
real(8) :: dt ! Time step of the data

real(8) :: profiles(1:nlevels, 1:nfiles)
real(8) :: profiles_sum(1:nlevels, 1:nfiles)
real(8) :: levels(1:nlevels)

real(8) :: time

integer(4) :: data_gaps(1:nlevels, 1:nfiles)

character(len=40) :: profile_files(1:nfiles)
character(len=20) :: data_names(0:nfiles)
character(len=40) :: file_name ! Working variable

data_names(0) = 'Depth'

!profile_files(1) = 'FX_LAKE_SP_Minlake_F2.dat'
!data_names(1) = 'Minlake'

profile_files(1) = 'KJ_keps_SP_20m_dz1m_Ndt6_F2.dat'
data_names(1) = 'LAKEoneD'

!profile_files(1) = 'AM_FL_SP_20m_NOSNOW_F2.DAT'
!profile_files(2) = 'AM_HS_SP_20m_NOSNOW_F2.DAT'
!profile_files(3) = 'VS_LAKE_SP_20M_F2.DAT'
!profile_files(4) = 'MP_Keps_SP_20m_F2.dat'
!profile_files(5) = 'SP0205WATERTHR.DAT'

!data_names(1) = 'FLake'
!data_names(2) = 'Hostetler'
!data_names(3) = 'LAKE'
!data_names(4) = 'Simstrat'
!data_names(5) = 'Observed'

! The initial time moment
year = 2002 
month = 1
day = 1
hour = 1.

! Termination year and month
year_term = 2005
month_term = 12

! Timestep of data in files with profiles
dt = 3600. ! 1 hour
!dt = 24.*3600. ! 1 day

do i = 1, nfiles
  open(unit = i, file = profile_files(i), status = 'old')
enddo

profiles_sum(:,:) = 0.
month_old = month
year_old = year
k = nfiles + 1
data_gaps(:,:) = 0
do while (.not.(year == year_term .and. month == month_term))

  do i = 1, nfiles
    read(i,*) time, ( (levels(j), profiles(j,i)), j = 1, nlevels)
    do j = 1, nlevels
      profiles_sum(j,i) = profiles_sum(j,i) + profiles(j,i)
      if (profiles(j,i) == 9999) data_gaps(j,i) = 1
    enddo
  enddo
  
  call JULIAN_DATE(year,month,day,hour,dt)
  if (month /= month_old) then
    month_sec = MONTH_DAYS(month_old, year)*day_sec
    profiles_sum(:,:) = profiles_sum(:,:)*dt/month_sec
    write(file_name,'(a7, i5, i3, i3, a1)') 'monthly', year_old, month_old, nfiles, 'M'
    open(unit = k, file = file_name, status = 'unknown')
    write (k, fmt = 20) (data_names(i), i = 0, nfiles)
    do j = 1, nlevels
      write (k, fmt = 10) levels(j), &
      & ( (profiles_sum(j,i)*(1-data_gaps(j,i)) + 9999*data_gaps(j,i)), i = 1, nfiles)
    enddo
    close(k)
    profiles_sum(:,:) = 0.
    month_old = month
    year_old = year
    data_gaps(:,:) = 0
  endif
  
enddo

do i = 1, nfiles
  close(i)
enddo

10 format (f10.2, <nfiles>f10.2)
20 format (<nfiles+1>a20)

END PROGRAM PROFILES_PROCESSING


SUBROUTINE JULIAN_DATE(year,month,day,hour,dt)

! The subroutine JULIAN_DATE updates the julian date
! after timestep dt

implicit none

real(8), parameter:: hour_dur = 3600.d0 ! The duration of hour, s 

!Input variables
real(8)   , intent(in)   :: dt  ! Timestep, s

!Input/output variables
integer(4), intent(inout):: year
integer(4), intent(inout):: month
integer(4), intent(inout):: day

real(8)   , intent(inout):: hour

!Local variables
integer(4) ndaym(12)
data       ndaym /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/

SAVE ndaym

if (mod(year,4)==0) then 
  ndaym(2) = 29 !Checking for leap-year
else
  ndaym(2) = 28
endif

hour = hour + dt/hour_dur
if (hour > 24.d0) then
  hour = hour - 24.d0
  day  = day  + 1
  if (day > ndaym(month)) then
    day   = 1
    month = month + 1
    if (month > 12) then
      month = 1
      year  = year + 1
    endif 
  endif
endif

END SUBROUTINE JULIAN_DATE


FUNCTION MONTH_DAYS(month, year)

integer(4) :: MONTH_DAYS

integer(4), intent(in) :: month, year

!Local variables
integer(4), save :: ndaym(12)
data       ndaym /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/

if (mod(year,4)==0) then 
  ndaym(2) = 29 !Checking for leap-year
else
  ndaym(2) = 28
endif

MONTH_DAYS = ndaym(month)

END FUNCTION MONTH_DAYS