module gms_math_plus use datatype use mem_manager implicit none integer :: ix, iy, iz interface operator(+) module procedure plus_x_x, plus_x_y, plus_x_z, plus_x_xy, plus_x_xz, plus_x_yz, plus_x_xyz module procedure plus_y_x, plus_y_y, plus_y_z, plus_y_xy, plus_y_xz, plus_y_yz, plus_y_xyz module procedure plus_z_x, plus_z_y, plus_z_z, plus_z_xy, plus_z_xz, plus_z_yz, plus_z_xyz module procedure plus_xy_x, plus_xy_y, plus_xy_z, plus_xy_xy, plus_xy_xz, plus_xy_yz, plus_xy_xyz module procedure plus_xz_x, plus_xz_y, plus_xz_z, plus_xz_xy, plus_xz_xz, plus_xz_yz, plus_xz_xyz module procedure plus_yz_x, plus_yz_y, plus_yz_z, plus_yz_xy, plus_yz_xz, plus_yz_yz, plus_yz_xyz module procedure plus_xyz_x, plus_xyz_y, plus_xyz_z, plus_xyz_xy, plus_xyz_xz, plus_xyz_yz, plus_xyz_xyz module procedure plus_real_x, plus_real_y, plus_real_z, plus_real_xy, plus_real_xz, plus_real_yz, plus_real_xyz module procedure plus_x_real, plus_y_real, plus_z_real, plus_xy_real, plus_xz_real, plus_yz_real, plus_xyz_real module procedure plus_x, plus_y, plus_z, plus_xy, plus_xz, plus_yz, plus_xyz end interface contains function plus_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 plus_x_x function plus_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 plus_x_y function plus_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 plus_x_z function plus_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 plus_x_xy function plus_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 plus_x_xz function plus_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 plus_x_yz function plus_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 plus_x_xyz function plus_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 plus_x_real function plus_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 plus_y_x function plus_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 plus_y_y function plus_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 plus_y_z function plus_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 plus_y_xy function plus_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 plus_y_xz function plus_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 plus_y_yz function plus_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 plus_y_xyz function plus_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 plus_y_real function plus_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 plus_z_x function plus_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 plus_z_y function plus_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 plus_z_z function plus_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 plus_z_xy function plus_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 plus_z_xz function plus_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 plus_z_yz function plus_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 plus_z_xyz function plus_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 plus_z_real function plus_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 plus_xy_x function plus_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 plus_xy_y function plus_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 plus_xy_z function plus_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 plus_xy_xy function plus_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 plus_xy_xz function plus_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 plus_xy_yz function plus_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 plus_xy_xyz function plus_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 plus_xy_real function plus_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 plus_xz_x function plus_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 plus_xz_y function plus_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 plus_xz_z function plus_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 plus_xz_xy function plus_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 plus_xz_xz function plus_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 plus_xz_yz function plus_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 plus_xz_xyz function plus_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 plus_xz_real function plus_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 plus_yz_x function plus_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 plus_yz_y function plus_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 plus_yz_z function plus_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 plus_yz_xy function plus_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 plus_yz_xz function plus_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 plus_yz_yz function plus_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 plus_yz_xyz function plus_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 plus_yz_real function plus_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 plus_xyz_x function plus_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 plus_xyz_y function plus_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 plus_xyz_z function plus_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 plus_xyz_xy function plus_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 plus_xyz_xz function plus_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 plus_xyz_yz function plus_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 plus_xyz_xyz function plus_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 plus_xyz_real function plus_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 plus_real_x function plus_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 plus_real_y function plus_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 plus_real_z function plus_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 plus_real_xy function plus_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 plus_real_xz function plus_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 plus_real_yz function plus_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 plus_real_xyz function plus_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 plus_x function plus_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 plus_y function plus_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 plus_z function plus_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 plus_xy function plus_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 plus_xz function plus_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 plus_yz function plus_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 plus_xyz end module gms_math_plus