diff --git a/NEWS.md b/NEWS.md
index 6c8ac53..055aa41 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -120,6 +120,46 @@ contour2D_on_globe!(fig,
 CairoMakie.save("myfigure.pdf", fig)
 ```
 
+## Integrating `OutputVar` with respect to longitude or latitude
+
+You can use the `integrate_lon(var)`, `integrate_lat(var)`, or `integrate_lonlat(var)`
+functions for integrating along longitude, latitude, or both respectively. The bounds of
+integration are determined by the range of the dimensions longitude and latitude in `var`.
+The convention is that latitude ranges from -90° to 90° and longitude ranges from -180° to
+180° and the unit of both longitude and latitude should be degree. See the example of
+integrating over a sphere where the data is all ones to find the surface area of a sphere.
+
+```@julia integrate_lonlat
+julia> lon = collect(range(-179.5, 179.5, 360));
+
+julia> lat = collect(range(-89.5, 89.5, 180));
+
+julia> data = ones(length(lon), length(lat));
+
+julia> dims = OrderedDict(["lon" => lon, "lat" => lat]);
+
+julia> dim_attribs = OrderedDict([
+           "lon" => Dict("units" => "degrees_east"),
+           "lat" => Dict("units" => "degrees_north"),
+       ]);
+
+julia> attribs = Dict("long_name" => "f");
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> integrated_var = integrate_lonlat(var);
+
+julia> integrated_var.dims # no dimensions since longitude and latitude are integrated out
+OrderedDict{String, Vector{Float64}}()
+
+julia> integrated_var.data # approximately 4π (the surface area of a sphere)
+0-dimensional Array{Float64, 0}:
+12.566530113084296
+
+julia> long_name(integrated_var) # updated long name to reflect the data being integrated
+"f integrated over lon (-179.5 to 179.5degrees_east) and integrated over lat (-89.5 to 89.5degrees_north)"
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.
diff --git a/docs/src/api.md b/docs/src/api.md
index d227772..9c7a1c6 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -55,6 +55,9 @@ Var.dim_units
 Var.range_dim
 Var.resampled_as
 Var.convert_units
+Var.integrate_lonlat
+Var.integrate_lat
+Var.integrate_lon
 ```
 
 ## Utilities
diff --git a/docs/src/howdoi.md b/docs/src/howdoi.md
index 3a42088..42cae81 100644
--- a/docs/src/howdoi.md
+++ b/docs/src/howdoi.md
@@ -107,3 +107,43 @@ OrderedDict{String, Vector{Float64}} with 2 entries:
   "lon"      => [0.0, 1.0]
   "latitude" => [0.0, 1.0, 2.0]
 ```
+
+## How do I integrate the data in a `OutputVar` with respect to longitude or latitude?
+
+You can use the `integrate_lon(var)`, `integrate_lat(var)`, or `integrate_lonlat(var)`
+functions for integrating along longitude, latitude, or both respectively. The bounds of
+integration are determined by the range of the dimensions longitude and latitude in `var`.
+The convention is that latitude ranges from -90° to 90° and longitude ranges from -180° to
+180° and the unit of both longitude and latitude should be degree. See the example of
+integrating over a sphere where the data is all ones to find the surface area of a sphere.
+
+```@julia integrate_lonlat
+julia> lon = collect(range(-179.5, 179.5, 360));
+
+julia> lat = collect(range(-89.5, 89.5, 180));
+
+julia> data = ones(length(lon), length(lat));
+
+julia> dims = OrderedDict(["lon" => lon, "lat" => lat]);
+
+julia> dim_attribs = OrderedDict([
+           "lon" => Dict("units" => "degrees_east"),
+           "lat" => Dict("units" => "degrees_north"),
+       ]);
+
+julia> attribs = Dict("long_name" => "f");
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> integrated_var = integrate_lonlat(var);
+
+julia> integrated_var.dims # no dimensions since longitude and latitude are integrated out
+OrderedDict{String, Vector{Float64}}()
+
+julia> integrated_var.data # approximately 4π (the surface area of a sphere)
+0-dimensional Array{Float64, 0}:
+12.566530113084296
+
+julia> long_name(integrated_var) # updated long name to reflect the data being integrated
+"f integrated over lon (-179.5 to 179.5degrees_east) and integrated over lat (-89.5 to 89.5degrees_north)"
+```
\ No newline at end of file
diff --git a/src/ClimaAnalysis.jl b/src/ClimaAnalysis.jl
index aff9deb..229179c 100644
--- a/src/ClimaAnalysis.jl
+++ b/src/ClimaAnalysis.jl
@@ -4,6 +4,9 @@ import Reexport: @reexport
 include("Utils.jl")
 import .Utils
 
+include("Numerics.jl")
+@reexport using .Numerics
+
 include("Var.jl")
 @reexport using .Var
 include("Sim.jl")
diff --git a/src/Numerics.jl b/src/Numerics.jl
new file mode 100644
index 0000000..2caabed
--- /dev/null
+++ b/src/Numerics.jl
@@ -0,0 +1,132 @@
+module Numerics
+
+"""
+    isequispaced(arr::Vector)
+
+Return whether the array is equispaced or not.
+"""
+function _isequispaced(arr::Vector)
+    return all(diff(arr) .≈ arr[begin + 1] - arr[begin])
+end
+
+"""
+    _integrate_lon(data::AbstractArray, lon::AbstractVector; dims)
+
+Integrate out longitude from `data`. `data` has to be discretized on `lon`.
+
+`dims` indicates which axis of `data` is longitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+longitude should be degrees. The convention is that longitude ranges from -180° and 180°.
+"""
+function _integrate_lon(data::AbstractArray, lon::AbstractVector; dims)
+    length(lon) == 1 &&
+        error("Cannot integrate when longitude is a single point")
+    _isequispaced(lon) ?
+    int_weights = _integration_weights_lon_equispaced(lon) :
+    int_weights = _integration_weights_lon_left(lon)
+    return _integrate_over_angle(data, lon, dims, int_weights)
+end
+
+"""
+    _integrate_lat(data::AbstractArray, lat::AbstractVector; dims)
+
+Integrate out latitude from `data`. `data` has to be discretized on `lat`.
+
+`dims` indicates which axis of `data` is latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+latitude should be degrees. The convention is that latitude ranges from -90° to 90°.
+"""
+function _integrate_lat(data::AbstractArray, lat::AbstractVector; dims)
+    length(lat) == 1 &&
+        error("Cannot integrate when latitude is a single point")
+    _isequispaced(lat) ?
+    int_weights = _integration_weights_lat_equispaced(lat) :
+    int_weights = _integration_weights_lat_left(lat)
+    return _integrate_over_angle(data, lat, dims, int_weights)
+end
+
+"""
+    _integrate_over_angle(
+    data::AbstractArray,
+    angle_arr::AbstractVector,
+    angle_idx,
+    int_weights::AbstractVector,
+)
+
+Integrate out angle (latitude or longitude) from `data` using the weights `int_weights`.
+`data` has to be discretized on `angle_arr`.
+
+`angle_idx` indicates which axis of `data` is angle.
+"""
+function _integrate_over_angle(
+    data::AbstractArray,
+    angle_arr::AbstractVector,
+    angle_idx,
+    int_weights::AbstractVector,
+)
+    # Reshape to add extra dimensions for int_weights for broadcasting if needed
+    size_to_reshape =
+        (i == angle_idx ? length(angle_arr) : 1 for i in 1:ndims(data))
+    int_weights = reshape(int_weights, size_to_reshape...)
+    int_on_angle = sum(data .* int_weights, dims = angle_idx)
+    return int_on_angle
+end
+
+"""
+    _integration_weights_lon_left(lon::AbstractVector)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along longitude.
+"""
+function _integration_weights_lon_left(lon::AbstractVector)
+    d_lon = deg2rad.(diff(lon))
+    # We are doing integration using the left endpoints, so we weight the rightmost endpoint
+    # zero so that it make no contribution to the integral
+    push!(d_lon, zero(eltype(d_lon)))
+    return d_lon
+end
+
+"""
+    _integration_weights_lat_left(lat::AbstractVector)
+
+Return integration weights for rectangular integration using left endpoints for integrating
+along latitude.
+"""
+function _integration_weights_lat_left(lat::AbstractVector)
+    d_lat = deg2rad.(diff(lat))
+    # We are doing integration using the left endpoints, so we weight the rightmost endpoint
+    # zero so that it make no contribution to the integral
+    push!(d_lat, zero(eltype(d_lat)))
+    cos_lat = cosd.(lat)
+    return d_lat .* cos_lat
+end
+
+"""
+    _integration_weights_lon_equispaced(lon::AbstractVector)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along longitude.
+"""
+function _integration_weights_lon_equispaced(lon::AbstractVector)
+    return deg2rad.([lon[begin + 1] - lon[begin] for _ in lon])
+end
+
+"""
+    _integration_weights_lat_equispaced(lat::AbstractVector)
+
+Return integration weights for rectangular integration when the points are equispaced for
+integrating along latitude.
+"""
+function _integration_weights_lat_equispaced(lat::AbstractVector)
+    d_lat = deg2rad.([lat[begin + 1] - lat[begin] for _ in lat])
+    cos_lat = cosd.(lat)
+    return d_lat .* cos_lat
+end
+
+end
diff --git a/src/Var.jl b/src/Var.jl
index d3395bb..f2f2965 100644
--- a/src/Var.jl
+++ b/src/Var.jl
@@ -7,6 +7,7 @@ import Interpolations as Intp
 import Statistics: mean
 import NaNStatistics: nanmean
 
+import ..Numerics
 import ..Utils: nearest_index, seconds_to_prettystr, squeeze
 
 export OutputVar,
@@ -30,7 +31,10 @@ export OutputVar,
     range_dim,
     resampled_as,
     has_units,
-    convert_units
+    convert_units,
+    integrate_lonlat,
+    integrate_lon,
+    integrate_lat
 
 """
     Representing an output variable
@@ -785,6 +789,105 @@ function resampled_as(src_var::OutputVar, dest_var::OutputVar)
     )
 end
 
+"""
+    integrate_lonlat(var::OutputVar)
+
+Integrate `data` in `var` on longitude and latitude with a first-order scheme. `data` has to
+be discretized on longitude and latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints for integrating
+longitude and latitude. The units for longitude and latitude should be degrees. The
+convention is that longitude vary from -180° and 180° and latitude vary from -90° to 90°.
+"""
+function integrate_lonlat(var::OutputVar)
+    var_integrate_lon = var |> integrate_lon
+    # Update long name so that we get "...integrated over lat... and integrated over lon..."
+    # instead of "...integrated over lat... integrated over lon..."
+    if haskey(var_integrate_lon.attributes, "long_name")
+        var_integrate_lon.attributes["long_name"] *= " and"
+    end
+    return var_integrate_lon |> integrate_lat
+end
+
+"""
+    integrate_lon(var::OutputVar)
+
+Integrate `data` in `var` on longitude with a first-order scheme. `data` has to be
+discretized on longitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+longitude should be degrees. The convention is that longitude ranges from -180° and 180°.
+"""
+function integrate_lon(var::OutputVar)
+    has_longitude(var) || error("var does not has longitude as a dimension")
+    lon_name = longitude_name(var)
+    return _integrate_over_angle(var, Numerics._integrate_lon, lon_name)
+end
+
+"""
+    integrate_lat(var::OutputVar)
+
+Integrate `data` in `var` on latitude with a first-order scheme. `data` has to be
+discretized on latitude.
+
+If the points are equispaced, it is assumed that each point correspond to the midpoint of a
+cell which results in rectangular integration using the midpoint rule. Otherwise, the
+integration being done is rectangular integration using the left endpoints. The unit for
+latitude should be degrees. The convention is that latitude ranges from -90° to 90°.
+"""
+function integrate_lat(var::OutputVar)
+    has_latitude(var) || error("var does not has latitude as a dimension")
+    lat_name = latitude_name(var)
+    return _integrate_over_angle(var, Numerics._integrate_lat, lat_name)
+end
+
+"""
+    _integrate_over_angle(var::OutputVar, integrate_on, angle_dim_name)
+
+Integrate `data` in `var` on `angle_dim_name` with a first-order scheme. `data` has to be
+discretized on `angle_dim_name`.
+
+`angle_dim_name` is the name of the angle that is being integrated over. `integrate_on` is a
+function that compute the integration for data over `angle_dim_name`.
+"""
+function _integrate_over_angle(var::OutputVar, integrate_on, angle_dim_name)
+    # Enforce constraint that unit is degree because we compute integration weights assuming
+    # degrees (see functions _integration_weights_lon_left and
+    # _integration_weights_lon_equispaced for examples in Numerics.jl)
+    deg_unit_names = [
+        "degrees",
+        "degree",
+        "deg",
+        "degs",
+        "°",
+        "degrees_north",
+        "degrees_east",
+    ]
+    angle_dim_unit = dim_units(var, angle_dim_name)
+    angle_dim_unit in deg_unit_names ||
+        lowercase(angle_dim_unit) |> (x -> occursin("deg", x)) ||
+        error("The unit for $angle_dim_name is missing or is not degree")
+
+    integrated_var = _reduce_over(
+        integrate_on,
+        angle_dim_name,
+        var,
+        var.dims[angle_dim_name],
+    )
+
+    _update_long_name_generic!(
+        integrated_var,
+        var,
+        angle_dim_name,
+        "integrated",
+    )
+    return integrated_var
+end
+
 """
     overload_binary_op(op)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 94f54bf..0e65c89 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,15 +2,16 @@ using SafeTestsets
 using Test
 
 #! format: off
-# @safetestset "Aqua" begin @time include("aqua.jl") end
-# @safetestset "Docstrings" begin @time include("doctest.jl") end
-# @safetestset "Format" begin @time include("format.jl") end
+@safetestset "Aqua" begin @time include("aqua.jl") end
+@safetestset "Docstrings" begin @time include("doctest.jl") end
+@safetestset "Format" begin @time include("format.jl") end
 
-# @safetestset "Utils" begin @time include("test_Utils.jl") end
-# @safetestset "SimDir" begin @time include("test_Sim.jl") end
-# @safetestset "Atmos" begin @time include("test_Atmos.jl") end
-# @safetestset "OutputVar" begin @time include("test_Var.jl") end
-# @safetestset "MakieExt" begin @time include("test_MakieExt.jl") end
+@safetestset "Utils" begin @time include("test_Utils.jl") end
+@safetestset "Numerics" begin @time include("test_Numerics.jl") end
+@safetestset "SimDir" begin @time include("test_Sim.jl") end
+@safetestset "Atmos" begin @time include("test_Atmos.jl") end
+@safetestset "OutputVar" begin @time include("test_Var.jl") end
+@safetestset "MakieExt" begin @time include("test_MakieExt.jl") end
 @safetestset "GeoMakieExt" begin @time include("test_GeoMakieExt.jl") end
 #! format: on
 
diff --git a/test/test_Numerics.jl b/test/test_Numerics.jl
new file mode 100644
index 0000000..ad11d87
--- /dev/null
+++ b/test/test_Numerics.jl
@@ -0,0 +1,113 @@
+using Test
+import ClimaAnalysis
+
+@testset "integration weights for lon and lat" begin
+    # Integration weights for lon (not equispaced)
+    lon = [-180.0, -45.0, 100.0, 180.0]
+    lon_weights = [135.0, 145.0, 80.0, 0.0] .* (π / 180.0)
+    @test all(
+        isapprox.(
+            lon_weights,
+            ClimaAnalysis.Numerics._integration_weights_lon_left(lon),
+        ),
+    )
+
+    # Integration weights for lat (not equispaced)
+    lat = [-90.0, 20.0, 45.0, 90.0]
+    lat_weights = [110.0, 25.0, 45.0, 0.0] .* (π / 180.0) .* cosd.(lat)
+    @test all(
+        isapprox.(
+            lat_weights,
+            ClimaAnalysis.Numerics._integration_weights_lat_left(lat),
+        ),
+    )
+
+    # Integration weights for lon (not equispaced)
+    lon = collect(range(-180.0, 180.0, 5))
+    lon_weights = [90.0 for _ in lon] .* (π / 180.0)
+    @test all(
+        isapprox.(
+            lon_weights,
+            ClimaAnalysis.Numerics._integration_weights_lon_equispaced(lon),
+        ),
+    )
+
+    lat = collect(range(-90.0, 90.0, 5))
+    lat_weights = [45.0 for _ in lat] .* (π / 180.0) .* cosd.(lat)
+    @test all(
+        isapprox.(
+            lat_weights,
+            ClimaAnalysis.Numerics._integration_weights_lat_equispaced(lat),
+        ),
+    )
+end
+
+@testset "Integrating on lon and lat" begin
+    # Integrating only lon (non equispaced)
+    lon = collect(range(-180.0, 179.0, 100))
+    # Force array to be non equispaced for testing _integration_weights_lon
+    push!(lon, 180.0)
+    lon_data = ones(length(lon))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lon(lon_data, lon, dims = 1)[1],
+        2.0π,
+        atol = 0.01,
+    )
+
+    # Integrating only lat (non equispaced)
+    lat = collect(range(-90.0, 89.0, 100))
+    # Force array to be non equispaced for testing _integration_weights_lat
+    push!(lat, 90.0)
+    lat_data = ones(length(lat))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lat(lat_data, lat, dims = 1)[1],
+        2.0,
+        atol = 0.01,
+    )
+
+    # Integrating both lon and lat
+    data = ones(length(lat), length(lon))
+    integrated_lat = ClimaAnalysis.Numerics._integrate_lat(data, lat, dims = 1)
+    integrated_latlon =
+        ClimaAnalysis.Numerics._integrate_lon(integrated_lat, lon, dims = 1)
+
+    integrated_lon = ClimaAnalysis.Numerics._integrate_lon(data, lon, dims = 2)
+    integrated_lonlat =
+        ClimaAnalysis.Numerics._integrate_lat(integrated_lon, lat, dims = 1)
+
+    # Order of integration should not matter
+    @test isapprox(integrated_latlon[1], integrated_lonlat[1])
+    @test isapprox(integrated_latlon[1], 4π, atol = 0.01)
+
+    # Error checking for length of lon and lat and values in lon and lat
+    @test_throws "Cannot integrate when latitude is a single point" ClimaAnalysis.Numerics._integrate_lat(
+        lat_data,
+        [0.0],
+        dims = 1,
+    )
+
+    @test_throws "Cannot integrate when latitude is a single point" ClimaAnalysis.Numerics._integrate_lat(
+        lon_data,
+        [0.0],
+        dims = 1,
+    )
+
+    # Integrating only lon (non equispaced)
+    lon = collect(range(-179.5, 179.5, 360))
+    lon_data = ones(length(lon))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lon(lon_data, lon, dims = 1)[1],
+        2.0π,
+        atol = 0.01,
+    )
+
+    # Integrating only lat (non equispaced)
+    lat = collect(range(-89.5, 89.5, 180))
+    lat_data = ones(length(lat))
+    @test isapprox(
+        ClimaAnalysis.Numerics._integrate_lat(lat_data, lat, dims = 1)[1],
+        2.0,
+        atol = 0.01,
+    )
+
+end
diff --git a/test/test_Var.jl b/test/test_Var.jl
index 7243b96..d256baf 100644
--- a/test/test_Var.jl
+++ b/test/test_Var.jl
@@ -687,3 +687,137 @@ end
     )
 
 end
+
+@testset "Integrating on lat and lon" begin
+    # Tests for integrate_lon
+    lon = collect(range(-179.5, 179.5, 360))
+    lon_data = ones(length(lon))
+    lon_dims = OrderedDict(["lon" => lon])
+    lon_attribs = Dict("long_name" => "hi")
+    lon_dim_attribs = OrderedDict(["lon" => Dict("units" => "deg")])
+    var = ClimaAnalysis.OutputVar(
+        lon_attribs,
+        lon_dims,
+        lon_dim_attribs,
+        lon_data,
+    )
+    var_integrated_lon = ClimaAnalysis.Var.integrate_lon(var)
+
+    @test isapprox(var_integrated_lon.data[1], 2.0 * π, atol = 0.01)
+    @test var_integrated_lon.dims == OrderedDict()
+    @test var_integrated_lon.dim_attributes == OrderedDict()
+    @test "hi integrated over lon (-179.5 to 179.5deg)" ==
+          var_integrated_lon.attributes["long_name"]
+    @test_throws "var does not has latitude as a dimension" ClimaAnalysis.Var.integrate_lat(
+        var,
+    )
+
+    # Tests for integrate_lat
+    lat = collect(range(-89.5, 89.5, 180))
+    lat_data = ones(length(lat))
+    lat_dims = OrderedDict(["lat" => lat])
+    lat_attribs = Dict("long_name" => "hi")
+    lat_dim_attribs = OrderedDict(["lat" => Dict("units" => "deg")])
+    var = ClimaAnalysis.OutputVar(
+        lat_attribs,
+        lat_dims,
+        lat_dim_attribs,
+        lat_data,
+    )
+    var_integrated_lat = ClimaAnalysis.Var.integrate_lat(var)
+
+    @test isapprox(var_integrated_lat.data[1], 2.0, atol = 0.01)
+    @test var_integrated_lat.dims == OrderedDict()
+    @test var_integrated_lat.dim_attributes == OrderedDict()
+    @test "hi integrated over lat (-89.5 to 89.5deg)" ==
+          var_integrated_lat.attributes["long_name"]
+    @test_throws "var does not has longitude as a dimension" ClimaAnalysis.Var.integrate_lon(
+        var,
+    )
+
+    # Unit checking
+    dim_attribs_no_units = OrderedDict([
+        "lon" => Dict("units" => ""),
+        "lat" => Dict("units" => ""),
+    ])
+    var_lon_no_units = ClimaAnalysis.OutputVar(
+        lon_attribs,
+        lon_dims,
+        dim_attribs_no_units,
+        lon_data,
+    )
+    @test_throws ErrorException ClimaAnalysis.Var.integrate_lon(
+        var_lon_no_units,
+    )
+
+    var_lat_no_units = ClimaAnalysis.OutputVar(
+        lat_attribs,
+        lat_dims,
+        dim_attribs_no_units,
+        lat_data,
+    )
+    @test_throws ErrorException ClimaAnalysis.Var.integrate_lat(
+        var_lat_no_units,
+    )
+end
+
+@testset "Integrating on sphere" begin
+    # Integrate out all dimensions (lat and lon) from OutputVar
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    data = ones(length(lon), length(lat))
+    dims = OrderedDict(["lon" => lon, "lat" => lat])
+    attribs = Dict("long_name" => "hi")
+    dim_attribs = OrderedDict([
+        "lon" => Dict("units" => "deg"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+    integrated_var = ClimaAnalysis.Var.integrate_lonlat(var)
+    @test isapprox(integrated_var.data[1], 4 * π, atol = 0.1)
+    @test integrated_var.dims == OrderedDict()
+    @test integrated_var.dim_attributes == OrderedDict()
+
+    # Integrating out lon and lat to get time series data
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    time = collect(range(0.0, 10.0, 10))
+    data = ones(length(lat), length(time), length(lon))
+    dims = OrderedDict(["lat" => lat, "time" => time, "lon" => lon])
+    attribs = Dict("long_name" => "hi")
+    dim_attribs = OrderedDict([
+        "lat" => Dict("units" => "deg"),
+        "time" => Dict("units" => "days"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+    integrated_var = ClimaAnalysis.Var.integrate_lonlat(var)
+
+    @test all(
+        isapprox.(integrated_var.data, [4 * π for _ in 1:10], atol = 0.01),
+    )
+    @test integrated_var.dims == OrderedDict(["time" => time])
+    @test integrated_var.dim_attributes ==
+          OrderedDict(["time" => Dict("units" => "days")])
+    @test "hi integrated over lon (-179.5 to 179.5deg) and integrated over lat (-89.5 to 89.5deg)" ==
+          integrated_var.attributes["long_name"]
+
+    # Unit checking
+    dim_attribs_no_lon = OrderedDict([
+        "time" => Dict("units" => "days"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs_no_lon, data)
+    @test_throws "The unit for lon is missing or is not degrees" ClimaAnalysis.Var.integrate_lonlat(
+        var,
+    )
+
+    dim_attribs_no_lat = OrderedDict([
+        "time" => Dict("units" => "days"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs_no_lat, data)
+    @test_throws "The unit for lat is missing or is not degrees" ClimaAnalysis.Var.integrate_lonlat(
+        var,
+    )
+end