module gms_math_minus
  use datatype
  use mem_manager
  implicit none
  integer :: ix, iy, iz

  interface operator(-)
  module procedure minus_x_x, minus_x_y, minus_x_z, minus_x_xy, minus_x_xz, minus_x_yz, minus_x_xyz
  module procedure minus_y_x, minus_y_y, minus_y_z, minus_y_xy, minus_y_xz, minus_y_yz, minus_y_xyz
  module procedure minus_z_x, minus_z_y, minus_z_z, minus_z_xy, minus_z_xz, minus_z_yz, minus_z_xyz
  module procedure minus_xy_x, minus_xy_y, minus_xy_z, minus_xy_xy, minus_xy_xz, minus_xy_yz, minus_xy_xyz
  module procedure minus_xz_x, minus_xz_y, minus_xz_z, minus_xz_xy, minus_xz_xz, minus_xz_yz, minus_xz_xyz
  module procedure minus_yz_x, minus_yz_y, minus_yz_z, minus_yz_xy, minus_yz_xz, minus_yz_yz, minus_yz_xyz
  module procedure minus_xyz_x, minus_xyz_y, minus_xyz_z, minus_xyz_xy, minus_xyz_xz, minus_xyz_yz, minus_xyz_xyz
  module procedure minus_real_x, minus_real_y, minus_real_z, minus_real_xy, minus_real_xz, minus_real_yz, minus_real_xyz
  module procedure minus_x_real, minus_y_real, minus_z_real, minus_xy_real, minus_xz_real, minus_yz_real, minus_xyz_real
  module procedure minus_x, minus_y, minus_z, minus_xy, minus_xz, minus_yz, minus_xyz
  end interface
contains
  function minus_x_x(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_x)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_x(new_id)
    output%id = new_id 

        do ix = lb_axis1, ub_axis1

        work_x(ix, 1, 1, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = -1
  end function minus_x_x

  function minus_x_y(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = -1
  end function minus_x_y

  function minus_x_z(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input2%grid(3)
  end function minus_x_z

  function minus_x_xy(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = -1
  end function minus_x_xy

  function minus_x_xz(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input2%grid(3)
  end function minus_x_xz

  function minus_x_yz(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_x_yz

  function minus_x_xyz(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_x(ix, 1, 1, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_x_xyz

  function minus_x_real(input1, input2) result(output)
    type(var_x), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_x)             :: output 

    integer :: new_id

    call get_new_id_x(new_id)
    output%id = new_id 

        do ix = lb_axis1, ub_axis1

        work_x(ix, 1, 1, new_id) & 
                      = work_x(ix, 1, 1, input1%id) -input2
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = -1
  end function minus_x_real

  function minus_y_x(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_y_x

  function minus_y_y(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_y)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_y(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2

        work_y(1, iy, 1, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_y_y

  function minus_y_z(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_y_z

  function minus_y_xy(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_y_xy

  function minus_y_xz(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_y_xz

  function minus_y_yz(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_y_yz

  function minus_y_xyz(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_y(1, iy, 1, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_y_xyz

  function minus_y_real(input1, input2) result(output)
    type(var_y), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_y)             :: output 

    integer :: new_id

    call get_new_id_y(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2

        work_y(1, iy, 1, new_id) & 
                      = work_y(1, iy, 1, input1%id) -input2
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_y_real

  function minus_z_x(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
        end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_z_x

  function minus_z_y(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_z_y

  function minus_z_z(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_z)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_z(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3

        work_z(1, 1, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do

    output%grid(1) = -1
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_z_z

  function minus_z_xy(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_z_xy

  function minus_z_xz(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_z_xz

  function minus_z_yz(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_z_yz

  function minus_z_xyz(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_z_xyz

  function minus_z_real(input1, input2) result(output)
    type(var_z), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_z)             :: output 

    integer :: new_id

    call get_new_id_z(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3

        work_z(1, 1, iz, new_id) & 
                      = work_z(1, 1, iz, input1%id) -input2
        end do

    output%grid(1) = -1
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_z_real

  function minus_xy_x(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_xy_x

  function minus_xy_y(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_xy_y

  function minus_xy_z(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_xy_z

  function minus_xy_xy(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_xy_xy

  function minus_xy_xz(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_xy_xz

  function minus_xy_yz(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_xy_yz

  function minus_xy_xyz(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_xy_xyz

  function minus_xy_real(input1, input2) result(output)
    type(var_xy), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = work_xy(ix, iy, 1, input1%id) -input2
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = -1
  end function minus_xy_real

  function minus_xz_x(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_xz_x

  function minus_xz_y(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xz_y

  function minus_xz_z(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_xz_z

  function minus_xz_xy(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xz_xy

  function minus_xz_xz(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_xz_xz

  function minus_xz_yz(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xz_yz

  function minus_xz_xyz(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xz_xyz

  function minus_xz_real(input1, input2) result(output)
    type(var_xz), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = work_xz(ix, 1, iz, input1%id) -input2
        end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = -1
    output%grid(3) = input1%grid(3)
  end function minus_xz_real

  function minus_yz_x(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_x

  function minus_yz_y(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_y

  function minus_yz_z(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_z

  function minus_yz_xy(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_xy

  function minus_yz_xz(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_xz

  function minus_yz_yz(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_yz

  function minus_yz_xyz(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_xyz

  function minus_yz_real(input1, input2) result(output)
    type(var_yz), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = work_yz(1, iy, iz, input1%id) -input2
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_yz_real

  function minus_xyz_x(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_x(ix, 1, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_x

  function minus_xyz_y(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_y(1, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_y

  function minus_xyz_z(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_z(1, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_z

  function minus_xyz_xy(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_xy(ix, iy, 1, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_xy

  function minus_xyz_xz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_xz(ix, 1, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_xz

  function minus_xyz_yz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_yz(1, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_yz

  function minus_xyz_xyz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id
    if ( input1%grid(1) /= input2%grid(1) ) stop "(-)grid violation"
    if ( input1%grid(2) /= input2%grid(2) ) stop "(-)grid violation"
    if ( input1%grid(3) /= input2%grid(3) ) stop "(-)grid violation"

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) & 
                      - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_xyz

  function minus_xyz_real(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1 
    real(8), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = work_xyz(ix, iy, iz, input1%id) -input2
        end do
      end do
    end do

    output%grid(1) = input1%grid(1)
    output%grid(2) = input1%grid(2)
    output%grid(3) = input1%grid(3)
  end function minus_xyz_real

  function minus_real_x(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_x), intent(in) :: input2 
    type(var_x)             :: output 

    integer :: new_id

    call get_new_id_x(new_id)
    output%id = new_id 

        do ix = lb_axis1, ub_axis1

        work_x(ix, 1, 1, new_id) & 
                      = input1 - work_x(ix, 1, 1, input2%id)
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = -1
    output%grid(3) = -1
  end function minus_real_x

  function minus_real_y(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_y), intent(in) :: input2 
    type(var_y)             :: output 

    integer :: new_id

    call get_new_id_y(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2

        work_y(1, iy, 1, new_id) & 
                      = input1 - work_y(1, iy, 1, input2%id)
      end do

    output%grid(1) = -1
    output%grid(2) = input2%grid(2)
    output%grid(3) = -1
  end function minus_real_y

  function minus_real_z(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_z), intent(in) :: input2 
    type(var_z)             :: output 

    integer :: new_id

    call get_new_id_z(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3

        work_z(1, 1, iz, new_id) & 
                      = input1 - work_z(1, 1, iz, input2%id)
        end do

    output%grid(1) = -1
    output%grid(2) = -1
    output%grid(3) = input2%grid(3)
  end function minus_real_z

  function minus_real_xy(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_xy), intent(in) :: input2 
    type(var_xy)             :: output 

    integer :: new_id

    call get_new_id_xy(new_id)
    output%id = new_id 

      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xy(ix, iy, 1, new_id) & 
                      = input1 - work_xy(ix, iy, 1, input2%id)
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = -1
  end function minus_real_xy

  function minus_real_xz(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_xz), intent(in) :: input2 
    type(var_xz)             :: output 

    integer :: new_id

    call get_new_id_xz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
        do ix = lb_axis1, ub_axis1

        work_xz(ix, 1, iz, new_id) & 
                      = input1 - work_xz(ix, 1, iz, input2%id)
        end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = -1
    output%grid(3) = input2%grid(3)
  end function minus_real_xz

  function minus_real_yz(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_yz), intent(in) :: input2 
    type(var_yz)             :: output 

    integer :: new_id

    call get_new_id_yz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2

        work_yz(1, iy, iz, new_id) & 
                      = input1 - work_yz(1, iy, iz, input2%id)
        end do
      end do

    output%grid(1) = -1
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_real_yz

  function minus_real_xyz(input1, input2) result(output)
    real(8), intent(in) :: input1 
    type(var_xyz), intent(in) :: input2 
    type(var_xyz)             :: output 

    integer :: new_id

    call get_new_id_xyz(new_id)
    output%id = new_id 

    do iz = lb_axis3, ub_axis3
      do iy = lb_axis2, ub_axis2
        do ix = lb_axis1, ub_axis1

        work_xyz(ix, iy, iz, new_id) & 
                      = input1 - work_xyz(ix, iy, iz, input2%id)
        end do
      end do
    end do

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)
  end function minus_real_xyz

  function minus_x(input) result(output)
    type(var_x), intent(in) :: input 
    type(var_x)             :: output 

    output=input
    work_x(:,:,:,input%id) = - work_x(:,:,:,input%id)

  end function minus_x
  function minus_y(input) result(output)
    type(var_y), intent(in) :: input 
    type(var_y)             :: output 

    output=input
    work_y(:,:,:,input%id) = - work_y(:,:,:,input%id)

  end function minus_y
  function minus_z(input) result(output)
    type(var_z), intent(in) :: input 
    type(var_z)             :: output 

    output=input
    work_z(:,:,:,input%id) = - work_z(:,:,:,input%id)

  end function minus_z
  function minus_xy(input) result(output)
    type(var_xy), intent(in) :: input 
    type(var_xy)             :: output 

    output=input
    work_xy(:,:,:,input%id) = - work_xy(:,:,:,input%id)

  end function minus_xy
  function minus_xz(input) result(output)
    type(var_xz), intent(in) :: input 
    type(var_xz)             :: output 

    output=input
    work_xz(:,:,:,input%id) = - work_xz(:,:,:,input%id)

  end function minus_xz
  function minus_yz(input) result(output)
    type(var_yz), intent(in) :: input 
    type(var_yz)             :: output 

    output=input
    work_yz(:,:,:,input%id) = - work_yz(:,:,:,input%id)

  end function minus_yz
  function minus_xyz(input) result(output)
    type(var_xyz), intent(in) :: input 
    type(var_xyz)             :: output 

    output=input
    work_xyz(:,:,:,input%id) = - work_xyz(:,:,:,input%id)

  end function minus_xyz
end module gms_math_minus
