module test_windspeed_routines
  use fruit
  use windspeed_routines
  use windspeed_cliper, only : radcal
  use windspeed_model, only : coef, vpcal

  implicit none

  contains

  subroutine test_rcal()
    real, dimension(6) :: rlat1, rlon1
    real, dimension(8) :: lat_delta, lon_delta, angles
    real :: rnm, azdeg
    integer :: i, j

    rlat1 = (/0., 10., 20., 30., 40., 50./)
    rlon1 = (/0., 10., 20., 30., 40., 50./)
    lat_delta = (/ 0., 10., 10.,  10.,   0., -10., -10., -10. /)
    lon_delta = (/10., 10.,  0., -10., -10., -10.,   0.,  10./)
    angles = (/   0.,  45., 90., 135., 180.,  225., 270., 315./)
    do j=1, size(rlat1)
      do i=1, size(angles)
        call rcal(rlat1(j), rlon1(j), rlat1(j)+lat_delta(i), rlon1(j)+lon_delta(i), rnm, azdeg)
        call assert_equals(angles(i), azdeg, 15.5)
      end do
    end do
    end subroutine test_rcal


  subroutine test_shcal
    real, dimension(0:7) :: flon,flat,ftime
    real, dimension(0:7):: fspd,fhead
    integer ::  mt = 7

    flon = (/-77.6, -78.2, -78.7, -79.1, -79.3, -79.2, -78.5, -75.5/)
    flat = (/29.6, 28.9, 28.5, 28.6, 29.0, 31.0, 32.5, 35.0 /)
    ftime = (/  0,12,24,36,48,72,96,120/)

    call shcal(flon, flat, ftime, fspd, fhead, mt, mt)
    call assert_equals(4.4, fspd(0), 0.1)
    call assert_equals(2.18, fspd(4), 0.1)

    call assert_equals(216.8, fhead(0), 0.1)
    call assert_equals(336.3, fhead(4), 0.1)
  end subroutine test_shcal

  subroutine test_vpcal
    real :: rlat, vmax, asym, beta, th0, rmax, xt, c
    real :: e_asym, e_th0, e_rmax, e_xt
    coef( 1) = 17.0000
    coef( 2) = 0.0800
    coef( 3) = -1.0500
    coef( 4) = 1.0600
    coef( 5) = 0.2800
    coef( 6) = -0.0026
    coef( 7) = -0.0800
    coef( 8) = 0.1147
    coef( 9) = 0.0055
    coef(10) = -0.0010
    coef(11) = 36.1000
    coef(12) = -0.0492
    coef(13) = 0.5740

    rlat=31.0
    vmax=50
    c = 20.054245755768374
    e_asym=5.149539602249986
    e_th0=-3.576958043556793
    e_rmax=37.084
    e_xt=0.38369999999999993

    call vpcal(vmax,rlat,c,asym,th0,rmax,xt)
    call assert_equals(e_asym, asym, 0.1)
    call assert_equals(e_th0, th0, 0.1)
    call assert_equals(e_rmax, rmax, 0.1)
    call assert_equals(e_xt, xt, 0.1)
  end subroutine test_vpcal

    subroutine test_radcal
    ! This routine calculates the wind radii from the parameteric vortex
      real, dimension(4) :: r34,r50,r64,r100
      real, dimension(4) :: e34
      integer :: i
      real :: ap, rmwp, thp
      real :: rlat, vmax, asym, beta, th0, rmax, xt
      r34=0
      r50=0
      r64=0
      r100=0
      ap=0
      rmwp=0
      thp=0
      rlat=31.0
      vmax=50
      asym=5.149539602249986
      beta=78.68008707289383
      th0=-3.576958043556793
      rmax=37.084
      xt=0.38369999999999993
      e34(1)= 93.778605
      e34(2)=110.064238
      e34(3)=63.068076
      e34(4)=55.375533
      call radcal(rlat,vmax,asym,beta,th0,rmax,xt,&
        & r34,r50,r64,r100,ap,rmwp,thp)
      do i=1, 4
        call assert_equals(e34(i), r34(i), 0.5)
      end do
    end subroutine test_radcal

    subroutine test_upcase
      character :: charr
      Character(26), Parameter :: ucl = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
      Character(26), Parameter :: lcl = 'abcdefghijklmnopqrstuvwxyz'
      integer i

      do i=1, 26
        charr = lcl(i:i)
        call uppcase(charr, 1)
        call assert_equals(ucl(i:i), charr)
      end do

      do i=1, 26
        charr = ucl(i:i)
        call uppcase(charr, 1)
        call assert_equals(ucl(i:i), charr)
      end do

    end subroutine

    subroutine test_lcase
      character :: charr
      Character(26), Parameter :: ucl = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
      Character(26), Parameter :: lcl = 'abcdefghijklmnopqrstuvwxyz'
      integer i

      do i=1, 26
        charr = ucl(i:i)
        call lcase(charr, 1)
        call assert_equals(lcl(i:i), charr)
      end do

      do i=1, 26
        charr = lcl(i:i)
        call lcase(charr, 1)
        call assert_equals(lcl(i:i), charr)
      end do

    end subroutine test_lcase

    subroutine test_ctorh
      real, dimension(9) :: x, y,exp_r, exp_theta
      real ::  r, theta
      integer :: i

      x = (/1.0, 1.0, 1.0, 0.0, -1.0, -1.0, -1.0, 0., 0./)
      y = (/1.0, 0.0, -1.0, -1.0, -1.0, 0., 1.0, 1.0, 0./)
      exp_theta = (/45., 90., 135., 180., 225., 270., 315., 0., 0./)
      exp_r = (/1.42, 1.0, 1.42, 1., 1.42, 1., 1.42, 1., 0./)


      r=-999.0
      theta =-999.0
      do i=1, 8
        call ctorh(x(i),y(i),r,theta)
        call assert_equals(exp_r(i), r, 0.1)
        call assert_equals(exp_theta(i), theta, 0.1)
      end do
    end subroutine test_ctorh

    subroutine test_rhtoc
      real, dimension(8) :: exp_x, exp_y,r, theta
      real ::  x, y
      integer :: i

      exp_x = (/1.0, 1.0, 1.0, 0.0, -1.0, -1.0, -1.0, 0./)
      exp_y = (/1.0, 0.0, -1.0, -1.0, -1.0, 0., 1.0, 1.0/)
      theta = (/45., 90., 135., 180., 225., 270., 315., 0./)
      r = (/1.42, 1.0, 1.42, 1., 1.42, 1., 1.42, 1./)


      y=-999.0
      x =-999.0
      do i=1, 8
        call rhtoc(r(i),theta(i), x, y)
        call assert_equals(exp_x(i), x,  0.1)
        call assert_equals(exp_y(i), y, 0.1)
      end do
    end subroutine test_rhtoc

    subroutine test_distk
      real, dimension(2) ::  rlon1, rlon2, rlat1, rlat2, exp_rad

      real :: dx, dy, rad, ret_lon, ret_lat
      integer :: i

      rlon1 = (/10., 290./)
      rlon2 = (/15., 295./)
      rlat1 = (/40., 50./)
      rlat2 = (/60., 70./)
      exp_rad = (/2250, 2239/)
      dx = 0.0
      dy = 0.0
      rad = 0.0
      do i=1,2
        call distk(rlon1(i),rlat1(i),rlon2(i),rlat2(i),dx,dy,rad)
        call assert_equals(exp_rad(i), rad, 1.)
        call distki( rlon1(i), rlat1(i), dx, dy, ret_lon, ret_lat)
        call assert_equals( rlon2(i), ret_lon, 1.)
        call assert_equals( rlat2(i), ret_lat, 1.)
      end do
    end subroutine test_distk

    subroutine test_yint
      INTEGER :: n, init_flag, lin_flag, ierr
      PARAMETER (n = 11)
      REAL :: yout, xi
      REAL, DIMENSION(n) :: x, y
      x = (/ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10./)
      y = (/ 0., 0., 1., 2., 1., 0., -1., -2., -1., 0., 0. /)

      init_flag = 1
      lin_flag = 0
      xi=2.5
      yout = 0
      call yint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.5, yout, 0.05)

      lin_flag = 1
      xi=2.5
      yout = 0
      call yint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.75, yout, 0.05)

      lin_flag = 1
      xi=5
      yout = 0
      call yint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(0., yout, 0.05)

      init_flag = 0
      xi=5
      yout = 0
      call yint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(0., yout, 0.05)

    end subroutine test_yint

    subroutine test_xint
      INTEGER :: n, init_flag, lin_flag, ierr
      PARAMETER (n = 11)
      REAL :: yout, xi
      REAL, DIMENSION(n) :: x, y
      x = (/ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10./)
      y = (/ 0., 0., 1., 2., 1., 0., -1., -2., -1., 0., 0. /)

      init_flag = 1
      lin_flag = 0
      xi=2.5
      yout = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.5, yout, 0.05)

      lin_flag = 1
      xi=2.5
      yout = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.75, yout, 0.05)

      lin_flag = 1
      xi=5
      yout = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(0.,yout, 0.05)

      init_flag = 0
      xi=5
      yout = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(0., yout, 0.05)

    end subroutine test_xint

    subroutine test_xyint
      INTEGER :: n, init_flag, lin_flag, ierr
      PARAMETER (n = 11)
      REAL :: yout, xi
      REAL, DIMENSION(n) :: x, y, y2
      x = (/ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10./)
      y = (/ 0., 0., 1., 2., 1., 0., -1., -2., -1., 0., 0. /)
      y2 = (/ 1., 1., 2., 3., 2., 1., 0., -1., 0., 1., 1. /)

      init_flag = 1
      lin_flag = 0
      xi=2.5
      yout = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.5, yout, 0.05)

      call yint(x,y2,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(2.5, yout, 0.05)

      init_flag = 0
      call xint(x,y,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(1.5, yout, 0.05)

      call yint(x,y2,n,init_flag,lin_flag,xi,yout,ierr)
      call assert_equals(2.5, yout, 0.05)



    end subroutine test_xyint

    subroutine test_linfill
      real, dimension(10) :: f_arr
      real :: missvalue, fillvalue
      integer :: f_len, maxgap
      ! Test simple fills between missing sequence
      f_arr = (/ 1., 2., -999., 4., 5., -999., -999., 8., 9., 10./)
      missvalue = -999.
      fillvalue = 0.
      f_len =10

      maxgap = 3

      call linfill(f_arr, f_len, missvalue, fillvalue, maxgap)
      call assert_equals(6., f_arr(6), 0.1)
      call assert_equals(7., f_arr(7), 0.1)
      call assert_equals(3., f_arr(3), 0.1)
      call assert_equals(5., f_arr(5), 0.1)

      ! Test decay back to fill value
      f_arr = (/ -999., 2., -999., 4., 5., -999., 6., -999., -999., -999./)

      maxgap = 2

      call linfill(f_arr, f_len, missvalue, fillvalue, maxgap)
      call assert_equals(1., f_arr(1), 0.1)
      call assert_equals(5.5, f_arr(6), 0.1)
      call assert_equals(6., f_arr(7), 0.1)
      call assert_equals(0., f_arr(10), 0.1)

            ! Test decay back to fill value
      f_arr = (/ -999., 2., -999., 4., 5., -999., 7., -999., -999., -999./)

      maxgap = -2

      call linfill(f_arr, f_len, missvalue, fillvalue, maxgap)
      call assert_equals(1., f_arr(1), 0.1)
      call assert_equals(3., f_arr(3), 0.1)
      call assert_equals(6., f_arr(6), 0.1)
      call assert_equals(7., f_arr(7), 0.1)
      call assert_equals(7., f_arr(8), 0.1)
      call assert_equals(7., f_arr(9), 0.1)

      ! Test initial missing but later values extrapolation
      f_arr = (/ -999., -999., -999., 4., 5., -999., 7., -999., -999., -999./)

      maxgap = -2

      call linfill(f_arr, f_len, missvalue, fillvalue, maxgap)
      call assert_equals(4., f_arr(1), 0.1)
      call assert_equals(4., f_arr(3), 0.1)
      call assert_equals(6., f_arr(6), 0.1)
      call assert_equals(7., f_arr(7), 0.1)
      call assert_equals(7., f_arr(8), 0.1)
      call assert_equals(7., f_arr(9), 0.1)


      ! Test initial missing but later values interpolation
      f_arr = (/ -999., -999., 3., 4., 5., -999., 6., -999., -999., -999./)

      maxgap = 2

      call linfill(f_arr, f_len, missvalue, fillvalue, maxgap)
      call assert_equals(2., f_arr(2), 0.1)
      call assert_equals(3., f_arr(3), 0.1)
      call assert_equals(5.5, f_arr(6), 0.1)
      call assert_equals(6., f_arr(7), 0.1)
      call assert_equals(0., f_arr(10), 0.1)
    end subroutine test_linfill

    subroutine test_randn
      integer :: ret_rand

      call randn(ret_rand, 1 )
      call assert_true( ret_rand .ge. 0 .and. ret_rand .lt. 10)

      call randn(ret_rand, 2 )
      call assert_true( ret_rand .ge. 10 .and. ret_rand .lt. 100)

      call randn(ret_rand, 3 )
      call assert_true( ret_rand .ge. 100 .and. ret_rand .lt. 1000)

      call randn(ret_rand, 4 )
      call assert_true( ret_rand .ge. 1000 .and. ret_rand .lt. 10000)

    end subroutine test_randn

    subroutine test_rand_distribuion
      integer, parameter  :: N=1000, NTEST=4
      integer :: i, j
      integer, dimension(N) :: ret_rand
      real :: mean
      real, dimension(NTEST) :: answers, accuracies

      answers = (/5., 50., 500., 5000./)
      accuracies = (/1., 2., 20., 200./)

      do j=1, NTEST
        do i=1, N
          call randn(ret_rand(i), 1 )
        end do
        mean  = sum(ret_rand) / N
        call assert_equals( 5., mean, 1.)
      end do

      do i=1, N
        call randn(ret_rand(i), 2 )
      end do
      mean  = sum(ret_rand) / N
      call assert_equals( 50., mean, 2.)

      do i=1, N
        call randn(ret_rand(i), 3 )
      end do
      mean  = sum(ret_rand) / N
      call assert_equals( 500., mean, 20.)

      do i=1, N
        call randn(ret_rand(i), 4 )
      end do
      mean  = sum(ret_rand) / N
      call assert_equals( 5000., mean, 200.)

    end subroutine test_rand_distribuion


    subroutine test_sindex
      integer :: i, ni, j, nj

      i = 50
      ni = 100

      nj = 50
      call sindex(i, ni, j, nj)
      call assert_equals(25, j)

      i = 50
      ni = 100
      nj = ni
      call sindex(i, ni, j, nj)
      call assert_equals(i, j)
    end subroutine test_sindex

    subroutine test_tolck
      real :: x1, x2, tol
      integer :: idiff

      ! test idiff if tolerance is ok
      x1 = 25
      x2 = 30
      tol = 5
      call tolck(x1,x2,tol,idiff)
      call assert_equals(0, idiff)

      call tolck(x2,x1,tol,idiff)
      call assert_equals(0, idiff)

      !  test that idiff is set when tolerance is exceeded
      x1 = 25
      x2 = 30
      tol = 2
      call tolck(x1,x2,tol,idiff)
      call assert_equals( 1, idiff)

      call tolck(x2,x1,tol,idiff)
      call assert_equals(-1, idiff)

    end subroutine test_tolck

    subroutine test_ofti
      integer iego,nx, ny, i, j, ij

      iego = 1
      nx = 10
      ny = 5
      i = 7
      j = 2
      call ofti(iego,nx,ny,i,j,ij)
      call assert_equals(17, ij)

      iego = 2
      call ofti(iego,nx,ny,i,j,ij)
      call assert_equals(37, ij)

      iego = 3
      call ofti(iego,nx,ny,i,j,ij)
      call assert_equals(34, ij)

      iego = 4
      call ofti(iego,nx,ny,i,j,ij)
      call assert_equals(14, ij)
    end subroutine test_ofti

    subroutine test_rint
      real :: azdeg(4), rad(4), radt
      integer :: i

      rad = (/ 45, 315, 225, 135/)
      azdeg = (/90., 100., 180., 250./)

      do i =1, 4
        call rint(azdeg(i),rad,radt)
        call assert_equals(azdeg(i), radt, 0.00001)
      end do


    end subroutine test_rint

end module test_windspeed_routines