!-----------------------------------------------------------------------
      module module_timeseries
!-----------------------------------------------------------------------

      use module_kinds
      use module_constants
      use module_my_domain_specs
      use module_vars

!-----------------------------------------------------------------------

      implicit none

!-----------------------------------------------------------------------

      private

      public :: timeseries_initialize, timeseries_run, timeseries_finalize

      integer, parameter :: max_point=20
      logical, save :: initialized=.false.
      integer, save :: npoints
      real, dimension(max_point), save :: points_lon, points_lat
      integer, dimension(max_point), save :: ipnt, jpnt
      real, dimension(max_point), save :: pnt_lon, pnt_lat

      real, parameter :: rtd=180.0/pi

      real, allocatable, dimension(:,:) :: zint   ! height of full levels
      real, allocatable, dimension(:,:) :: zmid   ! height of half levels

      integer, dimension(max_point), save :: tsunit

      integer :: var2d_number
      character(len=8), dimension(1000) :: var2d_name
      character(len=24), dimension(1000) :: var2d_units
      character(len=64), dimension(1000) :: var2d_description
      integer :: var3d_number
      character(len=8), dimension(1000) :: var3d_name
      character(len=24), dimension(1000) :: var3d_units
      character(len=64), dimension(1000) :: var3d_description
      character(len=4), dimension(1000) :: var3d_lvlind

      logical :: nml_exist

      integer, parameter :: max_fulllevel_vars = 10
      character(len=8), dimension(max_fulllevel_vars) :: fulllevel_vars &
           = reshape ( (/                                               &
                              'PINT    '                                &
                       /)                                               &
                      ,(/max_fulllevel_vars/)                           &
                      ,(/'********'/)                                   &
                     )
!-----------------------------------------------------------------------

      contains

!-----------------------------------------------------------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------

      subroutine timeseries_initialize(solver_state,ntsd,ierr)

      use module_solver_internal_state, only : solver_internal_state

      implicit none

      type(solver_internal_state),intent(in) :: solver_state
      integer, intent(in) :: ntsd
      integer, intent(out) :: ierr

      integer :: ios

      integer :: im, jm, lm
      real :: tlm0d,tph0d,dlmd,dphd,wbd,sbd
      real :: tlmd,tphd
      integer :: i,k,j, np, var, inrs,jnrs, nvar, n, l
      logical :: inside
      integer, dimension(8) :: modelstarttime
      character(len=3) :: trnum

      integer year, month, day, hour, minute, second, ten_thousandth
      integer :: ihr
      integer :: nf_hours,nf_minutes
      real    :: nf_seconds
      real    :: secfcst
      integer :: idatv(3),ihrv,idaywk

      character*128  :: filename

      namelist /ts_locations/ npoints, points_lon, points_lat


      ierr = 0

      if (initialized) return

      inquire ( file='ts_locations.nml', exist=nml_exist)

      if (.not.nml_exist) then
        if (mype == 0) then
          write(0,*) ' ts_locations.nml does not exist. will skip timeseries output'
        end if
        return
      end if

      open ( unit=80, file='ts_locations.nml', status='old', iostat = ios)
      if (ios /= 0) then
         write(0,*) 'error missing ts_locations.nml'
         ierr = -1
         return
      end if
      read  ( unit=80, nml=ts_locations,  iostat = ios)
      if (ios /= 0) then
         write(0,*) 'error reading ts_locations.nml'
         ierr = -1
         return
      end if
      close(unit=80,iostat=ios)
      if (ios /= 0) then
         write(0,*) 'error closing ts_locations.nml'
         ierr = -1
         return
      end if

      im = solver_state%im
      jm = solver_state%jm
      lm = solver_state%lm
      tlm0d = solver_state%tlm0d
      tph0d = solver_state%tph0d
      dlmd = solver_state%dlmd
      dphd = solver_state%dphd
      wbd = solver_state%wbd
      sbd = solver_state%sbd

!! calculate forecast time for this time step

      secfcst = solver_state%dt * ntsd
      nf_hours = int(secfcst/3600)
      ihr = nf_hours
      nf_minutes = int( mod(secfcst,3600.0)/60.0 )
      nf_seconds = (secfcst - nf_hours*3600.0) - nf_minutes*60.0
      if (nf_seconds< 0.000001) nf_seconds=0.0
      ten_thousandth = nint((secfcst-int(secfcst))*10000)

      call valid(solver_state%idat,solver_state%ihrst,nf_hours,idatv,ihrv,idaywk)

      year = idatv(3)
      month = idatv(2)
      day = idatv(1)
      hour = ihrv
      minute = nf_minutes
      second = nf_seconds

      modelstarttime(1) = year
      modelstarttime(2) = month
      modelstarttime(3) = day
      modelstarttime(4) = hour
      modelstarttime(5) = minute
      modelstarttime(6) = second
      modelstarttime(7) = 0
      modelstarttime(8) = ntsd

      allocate(zint(npoints,lm+1))
      allocate(zmid(npoints,lm))

      var2d_number = 0
      var3d_number = 0


!! loop over all solver state variables and find which ones are selected for timeseries output.
!! currently only 2d and 3d real variables are supported. this is due to limitations of used
!! timeseries binary file format

      do n=1,solver_state%num_vars
        if (solver_state%vars(n)%tseries) then
        select case(solver_state%vars(n)%tkr)
          case(tkr_r2d)
            var2d_number = var2d_number + 1
            var2d_name(var2d_number) = trim(solver_state%vars(n)%vbl_name)
            var2d_units(var2d_number) = ""
            var2d_description(var2d_number) = trim(solver_state%vars(n)%description)
          case(tkr_r3d)
            var3d_number = var3d_number + 1
            var3d_name(var3d_number) = trim(solver_state%vars(n)%vbl_name)
            var3d_units(var3d_number) = ""
            var3d_description(var3d_number) = trim(solver_state%vars(n)%description)
            var3d_lvlind(var3d_number) = 'H   '

            !! for 3d variables, level indicator is set to 'H' by default which means 'half' level
            !! variable, or variable located at the middle of the layer in vertical. if the variable name
            !! is included in the list of "full" level (or interface) variables then
            !! level indicator is reset to 'F'

            do l=1,max_fulllevel_vars
               if (trim(fulllevel_vars(l))==trim(var3d_name(var3d_number))) then
                  var3d_lvlind(var3d_number) = 'F   '
                  exit
               end if
            end do
          case default
            write(0,*)' unknown tkr = ', solver_state%vars(n)%tkr, trim(solver_state%vars(n)%vbl_name)
            ierr = -1
            return
        end select
        end if
      end do

!! loop over number of points specified in ts_locations namelist and calculate i,j
!! indexes of the nearest H point

      np_loop: do np=1,npoints

         call tll(points_lon(np),points_lat(np),tlmd,tphd,tph0d,tlm0d)
         call ijnrs(tlmd,tphd,dlmd,dphd,wbd,sbd,im,jm,inrs,jnrs,inside)
         if (inside) then
            ipnt(np)=inrs
            jpnt(np)=jnrs
         end if

!! if this point, with i,j indexes ipnt(np),jpnt(np) is inside the tile (its:ite,jts:jte)
!! located on this PE then proceed and create output file and write out file header

         inside_if: if (its <= ipnt(np) .and. ipnt(np) <= ite .and. jts <= jpnt(np) .and. jpnt(np) <= jte ) then

            i = ipnt(np)
            j = jpnt(np)
            zint(np,lm+1)=solver_state%fis(i,j)/g

            do k=lm,1,-1
                  zint(np,k)=zint(np,k+1)+solver_state%t(i,j,k)*(0.608*max(solver_state%q(i,j,k),epsq)+1.)*r_d        &
                            *(solver_state%pint(i,j,k+1)-solver_state%pint(i,j,k))                                     &
                            /((solver_state%sgml2(k)*solver_state%pd(i,j)+solver_state%psgml1(k))*g)
                  zmid(np,k)=0.5*(zint(np,k+1)+zint(np,k))
            end do

!! open the timeseries output file
            tsunit(np) = 90 + np
            write(filename,fmt='(a,i2.2,a,i2.2,a)') "ts_p",np,"_d",my_domain_id,".bin"
            write(0,*) ' open file:', np,tsunit(np),filename
            open(unit=tsunit(np), file=filename, form='unformatted', iostat=ios)
            if (ios /= 0) then
               write(0,*) 'error opening file '//trim(filename)//' in timeseries_initialize'
               ierr = -1
               return
            end if

!! header
            write(tsunit(np)) solver_state%dt
            write(tsunit(np)) modelstarttime
            write(tsunit(np)) 0  ! avgyn
            write(tsunit(np)) 0  ! avgfrq
            write(tsunit(np)) 0  ! avglen
            write(tsunit(np)) 0  ! avgfirst
            write(tsunit(np)) var2d_number
            do nvar=1,var2d_number
               write(tsunit(np)) nvar,var2d_name(nvar),var2d_units(nvar),var2d_description(nvar)
               write(0,*)        nvar,var2d_name(nvar),var2d_units(nvar),var2d_description(nvar)
            end do
            write(tsunit(np)) var3d_number
            do nvar=1,var3d_number
               write(tsunit(np)) nvar,var3d_name(nvar),var3d_lvlind(nvar),var3d_units(nvar),var3d_description(nvar)
               write(0,*)        nvar,var3d_name(nvar),var3d_lvlind(nvar),var3d_units(nvar),var3d_description(nvar)
            end do
            write(tsunit(np)) 1
            write(tsunit(np)) ipnt(np),jpnt(np), points_lat(np),points_lon(np)
            write(tsunit(np)) lm+1            ! number of full levels
            do k=1,lm+1
               write(tsunit(np)) zint(np,k)
            end do
            write(tsunit(np)) lm              ! number of half levels
            do k=1,lm
             write(tsunit(np)) zmid(np,k)
            end do
!! end of header

!! close the file
            close(unit=tsunit(np), iostat=ios)
            if (ios /= 0) then
               write(0,*) 'error closing '//filename
            end if

         end if inside_if

      end do np_loop

      initialized=.true.

      end subroutine timeseries_initialize

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine timeseries_run(solver_state,ntsd,ierr)

      use module_solver_internal_state, only : solver_internal_state

      implicit none

      type(solver_internal_state),intent(in) :: solver_state
      integer,intent(in) :: ntsd
      integer, intent(out) :: ierr

      integer i,j,k, lm, lmax, l
      integer n, np, var
      integer year, month, day, hour, minute, second, ten_thousandth
      real :: tlm0d,tph0d
      real :: pus,pvs

      integer :: ihr
      integer :: nf_hours,nf_minutes
      real    :: nf_seconds
      real    :: secfcst

      integer :: idatv(3),ihrv,idaywk

      integer :: ios
      character*128  :: filename

!!
      if (.not.nml_exist) then
        ierr = 0
        return
      end if

      tlm0d = solver_state%tlm0d
      tph0d = solver_state%tph0d
      lm = solver_state%lm

!! calculate forecast time for this time step

      secfcst = solver_state%dt * ntsd
      nf_hours = int(secfcst/3600)
      ihr = nf_hours
      nf_minutes = int( mod(secfcst,3600.0)/60.0 )
      nf_seconds = (secfcst - nf_hours*3600.0) - nf_minutes*60.0
      if (nf_seconds< 0.000001) nf_seconds=0.0
      ten_thousandth = nint((secfcst-int(secfcst))*10000)

      call valid(solver_state%idat,solver_state%ihrst,nf_hours,idatv,ihrv,idaywk)

      year = idatv(3)
      month = idatv(2)
      day = idatv(1)
      hour = ihrv
      minute = nf_minutes
      second = nf_seconds


      j_loop: do j = jts, jte
      i_loop: do i = its, ite

        np_loop: do np = 1, npoints

          if (i.eq.ipnt(np) .and. j.eq.jpnt(np)) then

            tsunit(np) = 90 + np
            write(filename,fmt='(a,i2.2,a,i2.2,a)') "ts_p",np,"_d",my_domain_id,".bin"
            write(0,*) ' open file:', np,tsunit(np),filename
            open(unit=tsunit(np), file=filename, form='unformatted', position='append', iostat=ios)
            if (ios /= 0) then
               write(0,*) 'error opening file '//trim(filename)//' in timeseries_run'
               ierr = -1
               return
            end if

            write(tsunit(np)) year,month,day,hour,minute,second,ten_thousandth,ntsd

! 2d
            do n=1,solver_state%num_vars
              if (solver_state%vars(n)%tseries .and. solver_state%vars(n)%tkr==tkr_r2d) then
                write(tsunit(np)) solver_state%vars(n)%r2d(i,j)
              end if
            end do

!3d
            do n=1,solver_state%num_vars
              if (solver_state%vars(n)%tseries .and. solver_state%vars(n)%tkr==tkr_r3d) then
                lmax = lm
                do l=1,max_fulllevel_vars
                  if (trim(fulllevel_vars(l))==trim(solver_state%vars(n)%vbl_name)) then
                    lmax = lm+1
                    exit
                  end if
                end do
                do k=1,lmax
                  write(tsunit(np)) solver_state%vars(n)%r3d(i,j,k)
                end do
                if(lmax==lm) write(tsunit(np)) 0.0
              end if
            end do

            close(unit=tsunit(np), iostat=ios)
            if (ios /= 0) then
               write(0,*) 'error closing '//filename
            end if

          end if

        end do np_loop

      end do i_loop
      end do j_loop

      end subroutine timeseries_run

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine timeseries_finalize
      end subroutine timeseries_finalize

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine tll(almd,aphd,tlmd,tphd,tph0d,tlm0d)
!-----------------------------------------------------------------------
      real, intent(in) :: almd, aphd
      real, intent(out) :: tlmd, tphd
      real, intent(in) :: tph0d, tlm0d
!-----------------------------------------------------------------------
      real, parameter :: pi=3.141592654
      real, parameter :: dtr=pi/180.0
!
      real :: tph0, ctph0, stph0, relm, srlm, crlm
      real :: aph, sph, cph, cc, anum, denom
!-----------------------------------------------------------------------

      tph0=tph0d*dtr
      ctph0=cos(tph0)
      stph0=sin(tph0)

      relm=(almd-tlm0d)*dtr
      srlm=sin(relm)
      crlm=cos(relm)
      aph=aphd*dtr
      sph=sin(aph)
      cph=cos(aph)
      cc=cph*crlm
      anum=cph*srlm
      denom=ctph0*cc+stph0*sph

      tlmd=atan2(anum,denom)/dtr
      tphd=asin(ctph0*sph-stph0*cc)/dtr

      return

      end subroutine tll

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine ijnrs(tlmpt,tphpt,dlmd,dphd,wbd,sbd,im,jm,inrs,jnrs,inside)

!     ******************************************************************
!     *                                                                *
!     *  code to:                                                      *
!     *  - find the i and j indices of nearest h point of the b        *
!     *    grid box containing point tlmpt,tphpt                       *
!     *                                                                *
!     *  ictp version - d.jovic                                        *
!     ******************************************************************

      real, intent(in) :: tlmpt,tphpt,dlmd,dphd,wbd,sbd
      integer, intent(in) :: im,jm
      integer, intent(out) :: inrs,jnrs
      logical, intent(out) :: inside

!-----------------------------------------------------------------------

      inrs=nint((tlmpt-wbd)/dlmd)+1
      jnrs=nint((tphpt-sbd)/dphd)+1

      if (jnrs >= 1 .and. jnrs <= jm .and. inrs >= 1 .and. inrs <= im ) then
         inside=.true.
      else
         inside=.false.
      endif

      return

      end subroutine ijnrs

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine tllwin(almd,aphd,tph0d,tlm0d,pus,pvs)

!     ******************************************************************
!     *                                                                *
!     *  ll to tll transformation of velocity                          *
!     *                                                                *
!     *  programer: z.janjic, yugoslav fed. hydromet. inst.,           *
!     *                 beograd, 1982                                  *
!     *                                                                *
!     ******************************************************************

      real, intent (in) :: almd, aphd
      real, intent (in) :: tph0d, tlm0d
      real, intent (inout) :: pus,pvs

      real, parameter :: pi=3.141592654
      real, parameter :: dtr=pi/180.0

      real :: relm,ctph0,stph0,tph0,srlm,crlm,ph,sph,cph,cc,tph
      real :: rctph,cray,dray
      real :: tpus,tpvs
      real :: rdenom

!-----------------------------------------------------------------------

      tph0=tph0d*dtr

      ctph0=cos(tph0)
      stph0=sin(tph0)

         relm=(almd-tlm0d)*dtr
         srlm=sin(relm)
         crlm=cos(relm)

         ph=aphd*dtr
         sph=sin(ph)
         cph=cos(ph)

         cc=cph*crlm
         tph=asin(ctph0*sph-stph0*cc)
         rctph=1.0/cos(tph)

         cray=stph0*srlm*rctph
         dray=(ctph0*cph+stph0*sph*crlm)*rctph

         tpus=pus
         tpvs=pvs
         rdenom=1.0/(dray*dray + cray*cray)
         pus=(dray*tpus+cray*tpvs)*rdenom
         pvs=(-cray*tpus+dray*tpvs)*rdenom

      return

      end subroutine tllwin

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine valid(idat,ihrst,ihr,idatv,ihrv,idaywk)
      integer idat(3),ihrst,ihr,idatv(3),ihrv,idaywk

      integer :: ijulian,iadd,iadday,idayyr

      ijulian=iw3jdn(idat(3),idat(2),idat(1))
      iadd=ihrst+ihr
      iadday=int((ihrst+ihr)/24)
      ihrv=iadd-24*iadday
      ijulian=ijulian+iadday
      call w3fs26(ijulian,idatv(3),idatv(2),idatv(1),idaywk,idayyr)
      return
      end subroutine valid

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      integer  function iw3jdn(iyear,month,iday)

      integer, intent(in) :: iyear,month,iday

      iw3jdn  = iday - 32075                                           &
              + 1461 * (iyear + 4800 + (month - 14) / 12) / 4          &
              + 367 * (month - 2 - (month -14) / 12 * 12) / 12         &
              - 3 * ((iyear + 4900 + (month - 14) / 12) / 100) / 4
      return
      end function iw3jdn

!-----------------------------------------------------------------------
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!-----------------------------------------------------------------------

      subroutine w3fs26(jldayn,iyear,month,iday,idaywk,idayyr)

      integer, intent(in) :: jldayn
      integer, intent(out) :: iyear, month, iday,idaywk,idayyr

      integer :: l,n,i,j

      l      = jldayn + 68569
      n      = 4 * l / 146097
      l      = l - (146097 * n + 3) / 4
      i      = 4000 * (l + 1) / 1461001
      l      = l - 1461 * i / 4 + 31
      j      = 80 * l / 2447
      iday   = l - 2447 * j / 80
      l      = j / 11
      month  = j + 2 - 12 * l
      iyear  = 100 * (n - 49) + i + l
      idaywk = mod((jldayn + 1),7) + 1
      idayyr = jldayn - (-31739 +1461 * (iyear+4799) / 4 - 3 * ((iyear+4899)/100)/4)
      return
      end subroutine w3fs26

!-----------------------------------------------------------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------

      end module module_timeseries

!-----------------------------------------------------------------------