module gms_math_mean
  use datatype
  use mem_manager
  implicit none
  interface mean_x
    module procedure mean_x_x
    module procedure mean_x_xy
    module procedure mean_x_xz
    module procedure mean_x_xyz
end interface

  interface mean_y
    module procedure mean_y_y
    module procedure mean_y_xy
    module procedure mean_y_yz
    module procedure mean_y_xyz
end interface

  interface mean_z
    module procedure mean_z_z
    module procedure mean_z_xz
    module procedure mean_z_yz
    module procedure mean_z_xyz
end interface

contains
  function mean_x_x(input) result(output)
    type(var_x), intent(in) :: input
    real(8) :: output
    if ( input%grid(1) == 0 ) then
      output & 
         = sum( work_x(plb1:pub1,1,1,input%id),1) / (pub1 - plb1 + 1)
    else 
      output & 
         = sum( work_x(plb1+1:pub1,1,1,input%id),1) / (pub1 - plb1 )
    end if
  end function mean_x_x

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

    call get_new_id_y(new_id)

    output%id = new_id

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

    if ( input%grid(1) == 0 ) then
      work_y(1,:,:,new_id) & 
         = sum( work_xy(plb1:pub1,:,:,input%id),1) / (pub1 - plb1 + 1)
    else 
      work_y(1,:,:,new_id) & 
         = sum( work_xy(plb1+1:pub1,:,:,input%id),1) / (pub1 - plb1 )
    end if
     
  end function mean_x_xy

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

    call get_new_id_z(new_id)

    output%id = new_id

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

    if ( input%grid(1) == 0 ) then
      work_z(1,:,:,new_id) & 
         = sum( work_xz(plb1:pub1,:,:,input%id),1) / (pub1 - plb1 + 1)
    else 
      work_z(1,:,:,new_id) & 
         = sum( work_xz(plb1+1:pub1,:,:,input%id),1) / (pub1 - plb1 )
    end if
     
  end function mean_x_xz

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

    call get_new_id_yz(new_id)

    output%id = new_id

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

    if ( input%grid(1) == 0 ) then
      work_yz(1,:,:,new_id) & 
         = sum( work_xyz(plb1:pub1,:,:,input%id),1) / (pub1 - plb1 + 1)
    else 
      work_yz(1,:,:,new_id) & 
         = sum( work_xyz(plb1+1:pub1,:,:,input%id),1) / (pub1 - plb1 )
    end if
     
  end function mean_x_xyz

  function mean_y_y(input) result(output)
    type(var_y), intent(in) :: input
    real(8) :: output
    if ( input%grid(2) == 0 ) then
      output & 
         = sum( work_y(1,plb2:pub2,1,input%id),1) / (pub2 - plb2 + 1)
    else 
      output & 
         = sum( work_y(1,plb2+1:pub2,1,input%id),1) / (pub2 - plb2 )
    end if
  end function mean_y_y

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

    call get_new_id_x(new_id)

    output%id = new_id

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

    if ( input%grid(2) == 0 ) then
      work_x(:,1,:,new_id) & 
         = sum( work_xy(:,plb2:pub2,:,input%id),2) / (pub2 - plb2 + 1)
    else 
      work_x(:,1,:,new_id) & 
         = sum( work_xy(:,plb2+1:pub2,:,input%id),2) / (pub2 - plb2 )
    end if
     
  end function mean_y_xy

  function mean_y_yz(input) result(output)
    type(var_yz), 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) = 1 
    output%grid(3) = input%grid(3)

    if ( input%grid(2) == 0 ) then
      work_z(:,1,:,new_id) & 
         = sum( work_yz(:,plb2:pub2,:,input%id),2) / (pub2 - plb2 + 1)
    else 
      work_z(:,1,:,new_id) & 
         = sum( work_yz(:,plb2+1:pub2,:,input%id),2) / (pub2 - plb2 )
    end if
     
  end function mean_y_yz

  function mean_y_xyz(input) result(output)
    type(var_xyz), 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) = 1 
    output%grid(3) = input%grid(3)

    if ( input%grid(2) == 0 ) then
      work_xz(:,1,:,new_id) & 
         = sum( work_xyz(:,plb2:pub2,:,input%id),2) / (pub2 - plb2 + 1)
    else 
      work_xz(:,1,:,new_id) & 
         = sum( work_xyz(:,plb2+1:pub2,:,input%id),2) / (pub2 - plb2 )
    end if
     
  end function mean_y_xyz

  function mean_z_z(input) result(output)
    type(var_z), intent(in) :: input
    real(8) :: output
    if ( input%grid(3) == 0 ) then
      output & 
         = sum( work_z(1,1,plb3:pub3,input%id),1) / (pub3 - plb3 + 1)
    else 
      output & 
         = sum( work_z(1,1,plb3+1:pub3,input%id),1) / (pub3 - plb3 )
    end if
  end function mean_z_z

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

    call get_new_id_x(new_id)

    output%id = new_id

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

    if ( input%grid(3) == 0 ) then
      work_x(:,:,1,new_id) & 
         = sum( work_xz(:,:,plb3:pub3,input%id),3) / (pub3 - plb3 + 1)
    else 
      work_x(:,:,1,new_id) & 
         = sum( work_xz(:,:,plb3+1:pub3,input%id),3) / (pub3 - plb3 )
    end if
     
  end function mean_z_xz

  function mean_z_yz(input) result(output)
    type(var_yz), 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) = 1 

    if ( input%grid(3) == 0 ) then
      work_y(:,:,1,new_id) & 
         = sum( work_yz(:,:,plb3:pub3,input%id),3) / (pub3 - plb3 + 1)
    else 
      work_y(:,:,1,new_id) & 
         = sum( work_yz(:,:,plb3+1:pub3,input%id),3) / (pub3 - plb3 )
    end if
     
  end function mean_z_yz

  function mean_z_xyz(input) result(output)
    type(var_xyz), 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) = 1 

    if ( input%grid(3) == 0 ) then
      work_xy(:,:,1,new_id) & 
         = sum( work_xyz(:,:,plb3:pub3,input%id),3) / (pub3 - plb3 + 1)
    else 
      work_xy(:,:,1,new_id) & 
         = sum( work_xyz(:,:,plb3+1:pub3,input%id),3) / (pub3 - plb3 )
    end if
     
  end function mean_z_xyz

end module gms_math_mean