! Read data files written from dataset of meteorological and soil data
! of Baltex project in Lindenberg (Germany) 2002-2009.
! Variuos data files contains data about various parameter in soil,
! at the surface and in the surface air layer.
! Program searchs data with the same date and time (date-time coincidence)
! and writes information about this date and time.
!

implicit none

integer, parameter :: ndata0_inp=200000, npar_inp=29, ndata0=500000

integer, parameter :: year_ini=2009, month_ini=8, day_ini=23, hour_ini=06, minute_ini=0, &
                      year_fin=2009, month_fin=9, day_fin=06, hour_fin=21, minute_fin=0

integer, dimension(npar_inp) :: ndata_inp
integer :: unit_inp=21, unit_out=22, ndata_out

character (len=30), dimension(npar_inp) :: file_name_inp=&
(/"atm_ps.dat     ", "atm_t.dat      ", "atm_q.dat      ", "fl_prec.dat    ", &
  "fl_sh.dat      ", "fl_lh.dat      ", &
  "fl_rswd.dat    ", "fl_rswu.dat    ", "fl_rlwd.dat    ", "fl_rlwu.dat    ", &
  "surf_t.dat     ", &
  "soil_t_0p05.dat", "soil_t_0p10.dat", "soil_t_0p15.dat", "soil_t_0p20.dat", &
  "soil_t_0p30.dat", "soil_t_0p45.dat", "soil_t_0p50.dat", "soil_t_0p60.dat", &
  "soil_t_0p90.dat", "soil_t_1p00.dat", "soil_t_1p20.dat", "soil_t_1p50.dat", &
  "soil_q_0p08.dat", "soil_q_0p15.dat", "soil_q_0p30.dat", &
  "soil_q_0p45.dat", "soil_q_0p60.dat", "soil_q_0p90.dat" /), &
 file_name_inp_control=(/ &
 "observ_atm_ps.dat     ", "observ_atm_t.dat      ", "observ_atm_q.dat      ", &
 "observ_fl_prec.dat    ", "observ_fl_sh.dat      ", "observ_fl_lh.dat      ", &
 "observ_fl_rswu.dat    ", "observ_fl_rswd.dat    ", "observ_fl_rlwu.dat    ", &
 "observ_fl_rlwd.dat    ", "observ_surf_t.dat     ", &
 "observ_soil_t_0p05.dat", "observ_soil_t_0p10.dat", "observ_soil_t_0p15.dat", &
 "observ_soil_t_0p20.dat", "observ_soil_t_0p30.dat", "observ_soil_t_0p45.dat", &
 "observ_soil_t_0p50.dat", "observ_soil_t_0p60.dat", "observ_soil_t_0p90.dat", &
 "observ_soil_t_1p00.dat", "observ_soil_t_1p20.dat", "observ_soil_t_1p50.dat", &
 "observ_soil_q_0p08.dat", "observ_soil_q_0p15.dat", "observ_soil_q_0p30.dat", &
 "observ_soil_q_0p45.dat", "observ_soil_q_0p60.dat", "observ_soil_q_0p90.dat"/)

character (len=30) :: file_name_out="observ_atm_forcing.dat2"

real :: datetime_ini, datetime_fin, nminyear, dtime_period
integer*8 :: idatetime_ini, idatetime_fin

real, dimension(ndata0_inp,npar_inp) :: datetime
integer*8, dimension(ndata0_inp,npar_inp) :: idatetime
integer, dimension(ndata0_inp,npar_inp) :: year, month, day, hour, minute
integer, dimension(npar_inp) :: i_ini=-999, i_fin=-999

real, dimension(ndata0_inp,npar_inp) :: value_inp

integer :: ierr_open=0, ierr_read=0, ifile_inp, nnn, dd_feb, ndyear, imon, jday, ipar_inp, i, ii
character (len=103) :: rline
character (len=1) :: a1
integer, dimension(12) :: dd_mon=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)

integer :: nstep, istep
real :: dtstep=1. ! minute
real, dimension(ndata0_inp) :: xinp, finp
real, dimension(ndata0) :: xout
real, dimension(ndata0,npar_inp) :: value_out 
real, dimension(ndata0) :: p, t, q, rl_rain, fl_snow, frtot, frvis, flh_sens, flh_latent, flq
real :: h_lev=2.4 ! FALKENBERG
!real :: h_lev=30.50 ! FOREST
real :: fsnow

 dd_feb=28
 ndyear=365
 if (mod(year_ini, 4) == 0) then
   dd_feb=29
   ndyear=366
 endif
 dd_mon(2)=dd_feb

 nminyear=float(ndyear)*1440.

 jday=0
 if (month_ini >= 2) then
 do imon=2,month_ini
   jday=jday+dd_mon(imon-1)
 enddo
 endif
 jday=jday+day_ini
 jday=jday-1

 datetime_ini=(float(jday)*1440.+float(hour_ini)*60.+float(minute_ini))/(float(ndyear)*1440.)
 idatetime_ini=year_ini*1000000+nint(datetime_ini*1000000)

 dd_feb=28
 ndyear=365
 if (mod(year_fin, 4) == 0) then
   dd_feb=29
   ndyear=366
 endif
 dd_mon(2)=dd_feb

 jday=0
 do imon=2,month_fin
   jday=jday+dd_mon(imon-1)
 enddo
 jday=jday+day_fin
 jday=jday-1

 datetime_fin=(float(jday)*1440.+float(hour_fin)*60.+float(minute_fin))/(float(ndyear)*1440.)
 idatetime_fin=year_fin*1000000+nint(datetime_fin*1000000)

 dtime_period=datetime_fin-datetime_ini

 nstep=int(dtime_period*nminyear)+2 ! minute number in defined period

 101 format (i4.4,f9.6,1x,f8.2,5x,i4.4,2i2.2,1x,i2.2,a1,i2.2)

 do ifile_inp=1,npar_inp

   open (unit_inp, file=trim(file_name_inp(ifile_inp)), form="formatted", status="old", iostat=ierr_open)

   if (ierr_open /= 0) then
     print *
     print *,"Requested file ",file_name_inp(ifile_inp)," not found, stop"
     stop
   endif

   read (unit_inp, *) rline

   nnn=0
   ierr_read=0
   do while (ierr_read == 0)
     nnn=nnn+1

     read (unit_inp, 101, iostat=ierr_read) year(nnn,ifile_inp),datetime(nnn,ifile_inp),value_inp(nnn,ifile_inp), &
 year(nnn,ifile_inp),month(nnn,ifile_inp),day(nnn,ifile_inp),hour(nnn,ifile_inp),a1,minute(nnn,ifile_inp)

     idatetime(nnn,ifile_inp)=year(nnn,ifile_inp)*1000000+nint(datetime(nnn,ifile_inp)*1000000)

     if (idatetime(nnn,ifile_inp) == idatetime_ini) i_ini(ifile_inp)=nnn
     if (idatetime(nnn,ifile_inp) > idatetime_ini.and.i_ini(ifile_inp) == -999) i_ini(ifile_inp)=nnn
     if (idatetime(nnn,ifile_inp) == idatetime_fin) i_fin(ifile_inp)=nnn
     if (idatetime(nnn,ifile_inp) > idatetime_fin.and.i_fin(ifile_inp) == -999) i_fin(ifile_inp)=nnn

     if (ierr_read /= 0) exit
   enddo

   ndata_inp(ifile_inp)=nnn-1

   print *,ndata_inp(ifile_inp),i_ini(ifile_inp),i_fin(ifile_inp)

   close (unit_inp)

 enddo

! Procedure of data interpolation in time

! X-axis is time (in minute) starting from start data-time determed

! X-axis for output interpolation (250920 minutes)

 do istep=1,nstep
   xout(istep)=float(istep-1)*dtstep
 enddo

 do ipar_inp=1,npar_inp

! X-axis and Function values for input interpolation 

   ii=0
   do i=i_ini(ipar_inp),i_fin(ipar_inp)
     ii=ii+1
     xinp(ii)=(datetime(i,ipar_inp)-datetime_ini)*nminyear
     finp(ii)=value_inp(i,ipar_inp)
   enddo

   nnn=ii

   call interp(1., 0.5, 0.5, nnn, xinp(1:nnn), finp(1:nnn), xout(1:nstep), value_out(1:nstep,ipar_inp), nstep)

! Control output

   open (unit_out, file=file_name_inp_control(ipar_inp), form="formatted", status="unknown")
   do ii=1,nnn
     write (unit_out,'(10f15.6)') xinp(ii),finp(ii)
   enddo

   close (unit_out)

!if (ipar_inp==3) then
!do ii=1,nnn
!  write (41,'(10f15.6)') xinp(ii),finp(ii)
!enddo
!
!do istep=1,nstep
!  write (42,'(10f15.6)') xout(istep),value_out(istep,ipar_inp)
!enddo
!endif

 enddo

! Definition of input data for Pochva model

 do istep=1,nstep

   p(istep)=value_out(istep,1)*1.e2 ! hPa -> Pa
   t(istep)=value_out(istep,2)+273.15 ! deg C -> K
   q(istep)=value_out(istep,3)*1.e-3 ! g/kg -> kg/kg
   rl_rain(istep)=value_out(istep,4)/1800. ! accum in kg/m^2 30 min -> rate kg/m^2/s
   fsnow=0.
   if (value_out(istep,2) > 5.) then
     fsnow=0.
   elseif (value_out(istep,2) < -5.) then
     fsnow=1.
   else
     fsnow=1.-(value_out(istep,2)+5.)/10.
   endif
   fl_snow(istep)=rl_rain(istep)*fsnow
   rl_rain(istep)=rl_rain(istep)*(1.-fsnow)
   frtot(istep)=value_out(istep,7)-value_out(istep,8)+value_out(istep,9)-value_out(istep,10)
   frvis(istep)=value_out(istep,7)-value_out(istep,8)
   flh_sens(istep)=-value_out(istep,5)
   flh_latent(istep)=-value_out(istep,6)
   flq(istep)=flh_latent(istep)/2.5008e6 ! Watt/m^2 -> kg/m^2/s 

!  write (43,'(10f15.6)') xout(istep),rl_rain(istep),fl_snow(istep), &
! frtot(istep),frvis(istep),flh_sens(istep),flh_latent(istep),flq(istep),t(istep),p(istep)
   
 enddo

 open (unit_out, file=file_name_out, form="formatted", status="unknown")

 do istep=2,nstep ! 1 step = time 0

   write (unit_out, 102) istep-1, h_lev, p(istep), t(istep), q(istep), &
 rl_rain(istep), fl_snow(istep), frtot(istep), frvis(istep), &
 flh_sens(istep), flh_latent(istep), flq(istep)

 enddo

 close (unit_out)

102 format (i8,100e16.8)

stop
end
!=======================================================================
subroutine interp(alfa, ex1, ex2, npi, xi, g, x, f, nval)

! (Similar to INTERP_SPLINE_1D, but with different input, and does not require
!  call of SUBR. NEAR - moreover, it checks that input coordinates are monotonic).

!  Interpolates with splines with tension in one dimension.
!  The spline is defined imposing that the second derivative is the average
!  of second derivatives computed at the two adjacent points.
!  At interval extremes the second derivative is assumed null.
!  This subroutine also extrapolates out of the interval where the input funtion G is defined

!  INPUT:  function G defined at coordinates XI (CAUTION: can be changed by this subroutine)
!          G(1:NPI) values at irregular but strictly growing coordinates XI(1:NPI)
!  OUTPUT: F(1:NVAL) interpolated values at arbitrary coordinates X(1:NVAL)

!  ALFA: spline tension parameter, comprised between 0 and 1
!  If ALFA=1, pure linear interpolation; if ALFA=0, pure spline

!  EX1: param. determining extrapolation for X < XI(1)
!  EX2: param. determining extrapolation for X > XI(NPI)
!  If EX1=0 OR EX2=0, constant value extrapolation is used at corresponding extreme
!  If EX1=1 OR EX2=1, linear extrapolation is used at corresponding extreme
!  Intermediate values of EX1 AND EX2 give intermediate extrapolation values

  real, dimension(npi )  :: xi, g
  real, dimension(nval)  :: x,  f

  if(alfa.lt..0.or.alfa.gt.1) then
  print*, 'CAUTION: in INTERP, ALFA out of interval 0-1'
  endif
  if(ex1.lt..0.or.ex1.gt.1) then
  print*, 'CAUTION: in INTERP, EX1 out of interval 0-1'
  endif
  if(ex2.lt..0.or.ex2.gt.1) then
  print*, 'CAUTION: in INTERP, EX2 out of interval 0-1'
  endif

! Fix for the case in which coordinates of the input function are not strictly increasing
! Note that this changes the input coordinates in the calling programme

  do k=2,npi
   if(xi(k).le.xi(k-1)) then
   print*, "CAUTION: in INTERP, coordinates of input function changed because not monotonic!"
   exit
   endif
  enddo

  zeps=(xi(npi)-xi(1))*1.e-6   ! Small deviation used to set apart interlaced coordinates
200 do k=2,npi
     if(xi(k).le.xi(k-1)) then
     ximed=0.5*(xi(k)+xi(k-1))
     xi(k-1)=ximed-zeps
     xi(k)=ximed+zeps
     gmed=0.5*(g(k)+g(k-1))
     g(k-1)=gmed
     g(k)=gmed
     endif
    enddo

 do k=2,npi
  if(xi(k).le.xi(k-1)) then
  goto 200
  endif
 enddo

 do 100 jval =1, nval

!  2 cases of extrapolation

 if(x(jval).lt.xi(1)) then
 f(jval) = g(1) + ex1*(g(1)-g(2))/(xi(1)-xi(2)) * (x(jval)-xi(1))
 go to 100
 elseif (x(jval).ge.xi(npi)) then
 f(jval) = g(npi) + ex2*(g(npi)-g(npi-1))/(xi(npi)-xi(npi-1)) * (x(jval)-xi(npi))
 go to 100
 endif

 ir = 0

!  IR is a reference index determining the interpolation interval
!  The interpolation expression is applied also if X=XI(J)

 do j = 1, npi
 if (x(jval).ge.xi(j)) ir = ir + 1
 enddo

 if (ir.eq.1) then
 fmm = 2*g(1) - g(2)
 xmm = 2*xi(1) - xi(2)
 fpp = g(ir+2)
 xpp = xi(ir+2)
 elseif (ir.eq.(npi-1)) then
 fpp = 2*g(npi) - g(npi-1)
 xpp = 2*xi(npi) - xi(npi-1)
 fmm = g(ir-1)
 xmm = xi(ir-1)
 else
 fmm = g(ir-1)
 xmm = xi(ir-1)
 fpp = g(ir+2)
 xpp = xi(ir+2)
 endif

 fm     = g(ir)
 xm     = xi(ir)
 fp     = g(ir+1)
 xp     = xi(ir+1)
 delx   = xp - xm
 delxp  = xpp - xp
 delxm  = xm - xmm
 delx1  = x(jval) - xm
 delx2  = xp - x(jval)
 delxs  = delx**2
 delx1s = delx1**2
 delx2s = delx2**2

!  Spline contribution to interpolation

 spl = fm*(delx2/delx + delx1*delx2s/(delxs*delxm) - delx1s*     &
       delx2/((delx+delxp)*delxs)) + fp*(delx1/delx +            &
       delx1s*delx2/(delxs*delxp) - delx1*delx2s/((delx+delxm)*  &
       delxs)) - fmm * delx1*delx2s/((delx+delxm)*delx*delxm) -  &
       fpp * delx1s*delx2/((delx+delxp)*delx*delxp)

!  Linear interpolation contribution

 clin = (fm*delx2 + fp*delx1)/delx

!  Final interpolation combined using ALFA

 f(jval) = alfa*clin + (1.-alfa)*spl

 100  continue

 return
end subroutine interp
!=======================================================================
