module gms_math_derivative
  use datatype
  use mem_manager
  implicit none
  interface d_x
    module procedure d_x_x
    module procedure d_x_xy
    module procedure d_x_xz
    module procedure d_x_xyz
  end interface

  interface d_y
    module procedure d_y_y
    module procedure d_y_xy
    module procedure d_y_yz
    module procedure d_y_xyz
  end interface

  interface d_z
    module procedure d_z_z
    module procedure d_z_xz
    module procedure d_z_yz
    module procedure d_z_xyz
  end interface

  interface d_2x
    module procedure d_2x_x
    module procedure d_2x_xy
    module procedure d_2x_xz
    module procedure d_2x_xyz
  end interface

  interface d_2y
    module procedure d_2y_y
    module procedure d_2y_xy
    module procedure d_2y_yz
    module procedure d_2y_xyz
  end interface

  interface d_2z
    module procedure d_2z_z
    module procedure d_2z_xz
    module procedure d_2z_yz
    module procedure d_2z_xyz
  end interface

  interface d_4x
    module procedure d_4x_x
    module procedure d_4x_xy
    module procedure d_4x_xz
    module procedure d_4x_xyz
  end interface

  interface d_4y
    module procedure d_4y_y
    module procedure d_4y_xy
    module procedure d_4y_yz
    module procedure d_4y_xyz
  end interface

  interface d_4z
    module procedure d_4z_z
    module procedure d_4z_xz
    module procedure d_4z_yz
    module procedure d_4z_xyz
  end interface

contains
  function d_x_x(input) result(output)
    type(var_x), intent(in) :: input
    type(var_x) :: output
    integer :: new_id

    call get_new_id_x(new_id)

    output%id = new_id

    output%grid(1) = mod( input%grid(1) + 1, 2 )
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_x(lb_axis1+output%grid(1) : ub_axis1-1+output%grid(1),:,:,new_id) & 
         = (   work_x(lb_axis1+1:ub_axis1,:,:,input%id) &
             - work_x(lb_axis1:ub_axis1-1,:,:,input%id) ) / dx
  end function d_x_x

  function d_2x_x(input) result(output)
    type(var_x), intent(in) :: input
    type(var_x) :: output
    integer :: new_id

    call get_new_id_x(new_id)

    output%id = new_id

    output%grid = input%grid

    work_x(lb_axis1+1 : ub_axis1-1,:,:,new_id) & 
         = (   work_x(lb_axis1+2:ub_axis1,:,:,input%id) &
             - work_x(lb_axis1:ub_axis1-2,:,:,input%id) ) / ( 2.0D0 * dx )
  end function d_2x_x

  function d_4x_x(input) result(output)
    type(var_x), intent(in) :: input
    type(var_x) :: output
    integer :: new_id

    call get_new_id_x(new_id)

    output%id = new_id

    output%grid = input%grid

    work_x(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
         = (   work_x(lb_axis1+4:ub_axis1,:,:,input%id) &
             - work_x(lb_axis1:ub_axis1-4,:,:,input%id) ) / ( 4.0D0 * dx )
  end function d_4x_x

  function d_x_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = mod( input%grid(1) + 1, 2 )
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xy(lb_axis1+output%grid(1) : ub_axis1-1+output%grid(1),:,:,new_id) & 
         = (   work_xy(lb_axis1+1:ub_axis1,:,:,input%id) &
             - work_xy(lb_axis1:ub_axis1-1,:,:,input%id) ) / dx
  end function d_x_xy

  function d_2x_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xy(lb_axis1+1 : ub_axis1-1,:,:,new_id) & 
         = (   work_xy(lb_axis1+2:ub_axis1,:,:,input%id) &
             - work_xy(lb_axis1:ub_axis1-2,:,:,input%id) ) / ( 2.0D0 * dx )
  end function d_2x_xy

  function d_4x_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xy(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
         = (   work_xy(lb_axis1+4:ub_axis1,:,:,input%id) &
             - work_xy(lb_axis1:ub_axis1-4,:,:,input%id) ) / ( 4.0D0 * dx )
  end function d_4x_xy

  function d_x_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = mod( input%grid(1) + 1, 2 )
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xz(lb_axis1+output%grid(1) : ub_axis1-1+output%grid(1),:,:,new_id) & 
         = (   work_xz(lb_axis1+1:ub_axis1,:,:,input%id) &
             - work_xz(lb_axis1:ub_axis1-1,:,:,input%id) ) / dx
  end function d_x_xz

  function d_2x_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xz(lb_axis1+1 : ub_axis1-1,:,:,new_id) & 
         = (   work_xz(lb_axis1+2:ub_axis1,:,:,input%id) &
             - work_xz(lb_axis1:ub_axis1-2,:,:,input%id) ) / ( 2.0D0 * dx )
  end function d_2x_xz

  function d_4x_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xz(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
         = (   work_xz(lb_axis1+4:ub_axis1,:,:,input%id) &
             - work_xz(lb_axis1:ub_axis1-4,:,:,input%id) ) / ( 4.0D0 * dx )
  end function d_4x_xz

  function d_x_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = mod( input%grid(1) + 1, 2 )
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xyz(lb_axis1+output%grid(1) : ub_axis1-1+output%grid(1),:,:,new_id) & 
         = (   work_xyz(lb_axis1+1:ub_axis1,:,:,input%id) &
             - work_xyz(lb_axis1:ub_axis1-1,:,:,input%id) ) / dx
  end function d_x_xyz

  function d_2x_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xyz(lb_axis1+1 : ub_axis1-1,:,:,new_id) & 
         = (   work_xyz(lb_axis1+2:ub_axis1,:,:,input%id) &
             - work_xyz(lb_axis1:ub_axis1-2,:,:,input%id) ) / ( 2.0D0 * dx )
  end function d_2x_xyz

  function d_4x_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid = input%grid

    work_xyz(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
         = (   work_xyz(lb_axis1+4:ub_axis1,:,:,input%id) &
             - work_xyz(lb_axis1:ub_axis1-4,:,:,input%id) ) / ( 4.0D0 * dx )
  end function d_4x_xyz

  function d_y_y(input) result(output)
    type(var_y), intent(in) :: input
    type(var_y) :: output
    integer :: new_id

    call get_new_id_y(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = mod( input%grid(2) + 1, 2 )
    output%grid(3) = input%grid(3)

    work_y(:,lb_axis2+output%grid(2) : ub_axis2-1+output%grid(2),:,new_id) & 
         = (   work_y(:,lb_axis2+1:ub_axis2,:,input%id) &
             - work_y(:,lb_axis2:ub_axis2-1,:,input%id) ) / dy
  end function d_y_y

  function d_2y_y(input) result(output)
    type(var_y), intent(in) :: input
    type(var_y) :: output
    integer :: new_id

    call get_new_id_y(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_y(:,lb_axis2+1 : ub_axis2-1,:,new_id) & 
         = (   work_y(:,lb_axis2+2:ub_axis2,:,input%id) &
             - work_y(:,lb_axis2:ub_axis2-2,:,input%id) ) / (2.0D0 * dy )
  end function d_2y_y

  function d_4y_y(input) result(output)
    type(var_y), intent(in) :: input
    type(var_y) :: output
    integer :: new_id

    call get_new_id_y(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_y(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
         = (   work_y(:,lb_axis2+4:ub_axis2,:,input%id) &
             - work_y(:,lb_axis2:ub_axis2-4,:,input%id) ) / (4.0D0 * dy )
  end function d_4y_y

  function d_y_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = mod( input%grid(2) + 1, 2 )
    output%grid(3) = input%grid(3)

    work_xy(:,lb_axis2+output%grid(2) : ub_axis2-1+output%grid(2),:,new_id) & 
         = (   work_xy(:,lb_axis2+1:ub_axis2,:,input%id) &
             - work_xy(:,lb_axis2:ub_axis2-1,:,input%id) ) / dy
  end function d_y_xy

  function d_2y_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xy(:,lb_axis2+1 : ub_axis2-1,:,new_id) & 
         = (   work_xy(:,lb_axis2+2:ub_axis2,:,input%id) &
             - work_xy(:,lb_axis2:ub_axis2-2,:,input%id) ) / (2.0D0 * dy )
  end function d_2y_xy

  function d_4y_xy(input) result(output)
    type(var_xy), intent(in) :: input
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xy(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
         = (   work_xy(:,lb_axis2+4:ub_axis2,:,input%id) &
             - work_xy(:,lb_axis2:ub_axis2-4,:,input%id) ) / (4.0D0 * dy )
  end function d_4y_xy

  function d_y_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = mod( input%grid(2) + 1, 2 )
    output%grid(3) = input%grid(3)

    work_yz(:,lb_axis2+output%grid(2) : ub_axis2-1+output%grid(2),:,new_id) & 
         = (   work_yz(:,lb_axis2+1:ub_axis2,:,input%id) &
             - work_yz(:,lb_axis2:ub_axis2-1,:,input%id) ) / dy
  end function d_y_yz

  function d_2y_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_yz(:,lb_axis2+1 : ub_axis2-1,:,new_id) & 
         = (   work_yz(:,lb_axis2+2:ub_axis2,:,input%id) &
             - work_yz(:,lb_axis2:ub_axis2-2,:,input%id) ) / (2.0D0 * dy )
  end function d_2y_yz

  function d_4y_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_yz(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
         = (   work_yz(:,lb_axis2+4:ub_axis2,:,input%id) &
             - work_yz(:,lb_axis2:ub_axis2-4,:,input%id) ) / (4.0D0 * dy )
  end function d_4y_yz

  function d_y_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = mod( input%grid(2) + 1, 2 )
    output%grid(3) = input%grid(3)

    work_xyz(:,lb_axis2+output%grid(2) : ub_axis2-1+output%grid(2),:,new_id) & 
         = (   work_xyz(:,lb_axis2+1:ub_axis2,:,input%id) &
             - work_xyz(:,lb_axis2:ub_axis2-1,:,input%id) ) / dy
  end function d_y_xyz

  function d_2y_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xyz(:,lb_axis2+1 : ub_axis2-1,:,new_id) & 
         = (   work_xyz(:,lb_axis2+2:ub_axis2,:,input%id) &
             - work_xyz(:,lb_axis2:ub_axis2-2,:,input%id) ) / (2.0D0 * dy )
  end function d_2y_xyz

  function d_4y_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xyz(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
         = (   work_xyz(:,lb_axis2+4:ub_axis2,:,input%id) &
             - work_xyz(:,lb_axis2:ub_axis2-4,:,input%id) ) / (4.0D0 * dy )
  end function d_4y_xyz

  function d_z_z(input) result(output)
    type(var_z), intent(in) :: input
    type(var_z) :: output
    integer :: new_id

    call get_new_id_z(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = mod( input%grid(3) + 1, 2 )

    work_z(:,:,lb_axis3+output%grid(3) : ub_axis3-1+output%grid(3),new_id) & 
         = (   work_z(:,:,lb_axis3+1:ub_axis3,input%id) &
             - work_z(:,:,lb_axis3:ub_axis3-1,input%id) ) / dz
  end function d_z_z

  function d_2z_z(input) result(output)
    type(var_z), intent(in) :: input
    type(var_z) :: output
    integer :: new_id

    call get_new_id_z(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_z(:,:,lb_axis3+1 : ub_axis3-1,new_id) & 
         = (   work_z(:,:,lb_axis3+2:ub_axis3,input%id) &
             - work_z(:,:,lb_axis3:ub_axis3-2,input%id) ) / ( 2.0D0 * dz)
  end function d_2z_z

  function d_4z_z(input) result(output)
    type(var_z), intent(in) :: input
    type(var_z) :: output
    integer :: new_id

    call get_new_id_z(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_z(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
         = (   work_z(:,:,lb_axis3+4:ub_axis3,input%id) &
             - work_z(:,:,lb_axis3:ub_axis3-4,input%id) ) / ( 4.0D0 * dz)
  end function d_4z_z

  function d_z_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = mod( input%grid(3) + 1, 2 )

    work_xz(:,:,lb_axis3+output%grid(3) : ub_axis3-1+output%grid(3),new_id) & 
         = (   work_xz(:,:,lb_axis3+1:ub_axis3,input%id) &
             - work_xz(:,:,lb_axis3:ub_axis3-1,input%id) ) / dz
  end function d_z_xz

  function d_2z_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xz(:,:,lb_axis3+1 : ub_axis3-1,new_id) & 
         = (   work_xz(:,:,lb_axis3+2:ub_axis3,input%id) &
             - work_xz(:,:,lb_axis3:ub_axis3-2,input%id) ) / ( 2.0D0 * dz)
  end function d_2z_xz

  function d_4z_xz(input) result(output)
    type(var_xz), intent(in) :: input
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
         = (   work_xz(:,:,lb_axis3+4:ub_axis3,input%id) &
             - work_xz(:,:,lb_axis3:ub_axis3-4,input%id) ) / ( 4.0D0 * dz)
  end function d_4z_xz

  function d_z_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = mod( input%grid(3) + 1, 2 )

    work_yz(:,:,lb_axis3+output%grid(3) : ub_axis3-1+output%grid(3),new_id) & 
         = (   work_yz(:,:,lb_axis3+1:ub_axis3,input%id) &
             - work_yz(:,:,lb_axis3:ub_axis3-1,input%id) ) / dz
  end function d_z_yz

  function d_2z_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_yz(:,:,lb_axis3+1 : ub_axis3-1,new_id) & 
         = (   work_yz(:,:,lb_axis3+2:ub_axis3,input%id) &
             - work_yz(:,:,lb_axis3:ub_axis3-2,input%id) ) / ( 2.0D0 * dz)
  end function d_2z_yz

  function d_4z_yz(input) result(output)
    type(var_yz), intent(in) :: input
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_yz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
         = (   work_yz(:,:,lb_axis3+4:ub_axis3,input%id) &
             - work_yz(:,:,lb_axis3:ub_axis3-4,input%id) ) / ( 4.0D0 * dz)
  end function d_4z_yz

  function d_z_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = mod( input%grid(3) + 1, 2 )

    work_xyz(:,:,lb_axis3+output%grid(3) : ub_axis3-1+output%grid(3),new_id) & 
         = (   work_xyz(:,:,lb_axis3+1:ub_axis3,input%id) &
             - work_xyz(:,:,lb_axis3:ub_axis3-1,input%id) ) / dz
  end function d_z_xyz

  function d_2z_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xyz(:,:,lb_axis3+1 : ub_axis3-1,new_id) & 
         = (   work_xyz(:,:,lb_axis3+2:ub_axis3,input%id) &
             - work_xyz(:,:,lb_axis3:ub_axis3-2,input%id) ) / ( 2.0D0 * dz)
  end function d_2z_xyz

  function d_4z_xyz(input) result(output)
    type(var_xyz), intent(in) :: input
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input%grid(1)
    output%grid(2) = input%grid(2)
    output%grid(3) = input%grid(3)

    work_xyz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
         = (   work_xyz(:,:,lb_axis3+4:ub_axis3,input%id) &
             - work_xyz(:,:,lb_axis3:ub_axis3-4,input%id) ) / ( 4.0D0 * dz)
  end function d_4z_xyz

end module gms_math_derivative
