From cdd7a5771063364df7d4d9c39f08c2092e3ec749 Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Fri, 6 Sep 2024 14:36:19 -0700
Subject: [PATCH 1/6] Fix formatting issue

---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index 9130e0b..439c225 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -120,7 +120,7 @@ contour2D_on_globe!(fig,
 CairoMakie.save("myfigure.pdf", fig)
 ```
 
-## Integrating `OutputVar` with respect to longitude or latitude
+### 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

From b406f63382bcefbbb032180e23d7ae1c76b8e940 Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Mon, 9 Sep 2024 10:21:36 -0700
Subject: [PATCH 2/6] Add Utils._data_at_dim_vals

---
 src/Utils.jl       | 30 ++++++++++++++++++++++++++++++
 test/test_Utils.jl | 13 +++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/src/Utils.jl b/src/Utils.jl
index 326ae5a..b927e4c 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -316,4 +316,34 @@ function _isequispaced(arr::Vector)
     return all(diff(arr) .≈ arr[begin + 1] - arr[begin])
 end
 
+"""
+    _data_at_dim_vals(data, dim_arr, dim_idx, vals)
+
+Return a view of `data` by slicing along `dim_idx`. The slices are indexed by the indices
+corresponding to values in `dim_arr` closest to the values in `vals`.
+
+Examples
+=========
+
+```jldoctest
+julia> data = [[1, 4, 7]  [2, 5, 8]  [3, 6, 9]];
+
+julia> dim_arr = [1.0, 2.0, 4.0];
+
+julia> dim_idx = 2;
+
+julia> vals = [1.1, 4.0];
+
+julia> Utils._data_at_dim_vals(data, dim_arr, dim_idx, vals)
+3×2 view(::Matrix{Int64}, :, [1, 3]) with eltype Int64:
+ 1  3
+ 4  6
+ 7  9
+```
+"""
+function _data_at_dim_vals(data, dim_arr, dim_idx, vals)
+    nearest_indices = map(val -> nearest_index(dim_arr, val), vals)
+    return selectdim(data, dim_idx, nearest_indices)
+end
+
 end
diff --git a/test/test_Utils.jl b/test/test_Utils.jl
index 4b1ecfb..a2783a1 100644
--- a/test/test_Utils.jl
+++ b/test/test_Utils.jl
@@ -103,3 +103,16 @@ end
     equispaced = Utils._isequispaced([0.0, 2.0, 3.0])
     @test equispaced == false
 end
+
+@testset "data_at_dim_vals" begin
+    data = [[1, 2, 3] [4, 5, 6] [7, 8, 9]]
+    dim_arr = [2.0, 3.0, 4.0]
+    dim_idx = 2
+
+    @test Utils._data_at_dim_vals(data, dim_arr, dim_idx, []) ==
+          reshape([], 3, 0)
+    @test Utils._data_at_dim_vals(data, dim_arr, dim_idx, [2.1]) ==
+          reshape([1; 2; 3], 3, 1)
+    @test Utils._data_at_dim_vals(data, dim_arr, dim_idx, [2.1, 2.9, 4.0]) ==
+          data
+end

From 7e62598b7abc7013a947cd1373a3f3212f864965 Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Mon, 9 Sep 2024 17:22:50 -0700
Subject: [PATCH 3/6] Add Base.is_empty(var::OutputVar)

---
 docs/src/api.md  |  1 +
 src/Var.jl       | 17 ++++++++++++++++-
 test/test_Var.jl | 20 ++++++++++++++++++++
 3 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index 9c7a1c6..c11219a 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -20,6 +20,7 @@ Sim.available_periods
 Var.OutputVar
 Var.read_var
 Var.is_z_1D
+Base.isempty(var::OutputVar)
 Var.short_name
 Var.long_name
 Var.units
diff --git a/src/Var.jl b/src/Var.jl
index 4259767..7175019 100644
--- a/src/Var.jl
+++ b/src/Var.jl
@@ -34,7 +34,8 @@ export OutputVar,
     convert_units,
     integrate_lonlat,
     integrate_lon,
-    integrate_lat
+    integrate_lat,
+    isempty
 
 """
     Representing an output variable
@@ -302,6 +303,20 @@ function is_z_1D(var::OutputVar)
     return length(size(altitudes(var))) == 1
 end
 
+"""
+    isempty(var::OutputVar)
+
+Determine whether an OutputVar is empty.
+"""
+function Base.isempty(var::OutputVar)
+    # Do not include :interpolant because var.interpolant is Nothing if data is
+    # zero dimensional and empty and isempty(Nothing) throws an error
+    return map(
+        fieldname -> isempty(getproperty(var, fieldname)),
+        filter(x -> x != :interpolant, fieldnames(OutputVar)),
+    ) |> all
+end
+
 function Base.copy(var::OutputVar)
     fields = fieldnames(OutputVar)
     # We have nested containers and we have to make sure we hand ownership off,
diff --git a/test/test_Var.jl b/test/test_Var.jl
index d451d92..7fe224c 100644
--- a/test/test_Var.jl
+++ b/test/test_Var.jl
@@ -85,6 +85,26 @@ import Unitful: @u_str
     )
 end
 
+@testset "empty" begin
+    dims = OrderedDict{String, Vector{Float64}}()
+    data = Float64[]
+    empty_var = ClimaAnalysis.OutputVar(dims, data)
+    @test ClimaAnalysis.isempty(empty_var)
+
+    dims = OrderedDict{String, Vector{Float64}}()
+    data = fill(1.0)
+    empty_var = ClimaAnalysis.OutputVar(dims, data)
+    @test !ClimaAnalysis.isempty(empty_var)
+
+    long = 0.0:180.0 |> collect
+    dims = OrderedDict(["long" => long])
+    data = ones(size(long))
+    dim_attributes = OrderedDict(["lon" => Dict("b" => 2)])
+    attribs = Dict("short_name" => "bob", "long_name" => "hi")
+    not_empty_var = ClimaAnalysis.OutputVar(attribs, dims, dim_attributes, data)
+    @test !ClimaAnalysis.isempty(not_empty_var)
+end
+
 @testset "Arithmetic operations" begin
     long = 0.0:180.0 |> collect
     lat = 0.0:90.0 |> collect

From 25f6f3304457ff337c89a83dba5d3889ce765e53 Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Mon, 9 Sep 2024 12:05:43 -0700
Subject: [PATCH 4/6] Add Utils.time_to_date and Utils.date_to_time

Also, add Utils.period_to_seconds_float. These functions are
functionally identical to the equivalent functions in ClimaUtilities. We
do not export from ClimaUtilities because we do not want to add
ClimaUtilities as a dependency.
---
 src/Utils.jl       | 105 +++++++++++++++++++++++++++++++++++++++++++++
 test/test_Utils.jl |  30 +++++++++++++
 2 files changed, 135 insertions(+)

diff --git a/src/Utils.jl b/src/Utils.jl
index b927e4c..9f380c1 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -316,6 +316,111 @@ function _isequispaced(arr::Vector)
     return all(diff(arr) .≈ arr[begin + 1] - arr[begin])
 end
 
+"""
+    time_to_date(start_date::Dates.DateTime, time::AbstractFloat)
+
+Convert the given time to a calendar date starting from `start_date`.
+
+Examples
+=========
+
+```jldoctest
+julia> import Dates
+
+julia> Utils.time_to_date(Dates.DateTime(2013, 7, 1, 12), 86400.0)
+2013-07-02T12:00:00
+
+julia> Utils.time_to_date(Dates.DateTime(2013, 7, 1, 12), 3600.0)
+2013-07-01T13:00:00
+
+julia> Utils.time_to_date(Dates.DateTime(2013, 7, 1, 12), 60.0)
+2013-07-01T12:01:00
+
+julia> Utils.time_to_date(Dates.DateTime(2013, 7, 1, 12), 1.0)
+2013-07-01T12:00:01
+```
+"""
+function time_to_date(start_date::Dates.DateTime, time::AbstractFloat)
+    # We go through milliseconds to allow fractions of a second (otherwise, Second(0.8)
+    # would fail). Milliseconds is the level of resolution that one gets when taking the
+    # difference between two DateTimes. In addition to this, we add a round to account for
+    # floating point errors. If the floating point error is small enough, round will correct
+    # it.
+    milliseconds = Dates.Millisecond.(round.(1_000 * time))
+    return start_date + milliseconds
+end
+
+"""
+    date_to_time(reference_date::Dates.DateTime, date::Dates.DateTime)
+
+Convert the given calendar date to a time (in seconds) where t=0 is `reference_date`.
+
+Examples
+=========
+
+```jldoctest
+julia> import Dates
+
+julia> Utils.date_to_time(Dates.DateTime(2013, 7, 1, 12), Dates.DateTime(2013, 7, 2, 12))
+86400.0
+
+julia> Utils.date_to_time(Dates.DateTime(2013, 7, 1, 12), Dates.DateTime(2013, 7, 1, 13))
+3600.0
+
+julia> Utils.date_to_time(Dates.DateTime(2013, 7, 1, 12), Dates.DateTime(2013, 7, 1, 12, 1))
+60.0
+
+julia> Utils.date_to_time(Dates.DateTime(2013, 7, 1, 12), Dates.DateTime(2013, 7, 1, 12, 0, 1))
+1.0
+```
+"""
+function date_to_time(reference_date::Dates.DateTime, date::Dates.DateTime)
+    return period_to_seconds_float(date - reference_date)
+end
+
+"""
+    period_to_seconds_float(period::Dates.Period)
+
+Convert the given `period` to seconds in Float64.
+
+Examples
+=========
+
+```jldoctest
+julia> import Dates
+
+julia> Utils.period_to_seconds_float(Dates.Millisecond(1))
+0.001
+
+julia> Utils.period_to_seconds_float(Dates.Second(1))
+1.0
+
+julia> Utils.period_to_seconds_float(Dates.Minute(1))
+60.0
+
+julia> Utils.period_to_seconds_float(Dates.Hour(1))
+3600.0
+
+julia> Utils.period_to_seconds_float(Dates.Day(1))
+86400.0
+
+julia> Utils.period_to_seconds_float(Dates.Week(1))
+604800.0
+
+julia> Utils.period_to_seconds_float(Dates.Month(1))
+2.629746e6
+
+julia> Utils.period_to_seconds_float(Dates.Year(1))
+3.1556952e7
+```
+"""
+function period_to_seconds_float(period::Dates.Period)
+    # See https://github.com/JuliaLang/julia/issues/55406
+    period isa Dates.OtherPeriod &&
+        (period = Dates.Second(Dates.Day(1)) * Dates.days(period))
+    return period / Dates.Second(1)
+end
+
 """
     _data_at_dim_vals(data, dim_arr, dim_idx, vals)
 
diff --git a/test/test_Utils.jl b/test/test_Utils.jl
index a2783a1..e59a1b8 100644
--- a/test/test_Utils.jl
+++ b/test/test_Utils.jl
@@ -104,6 +104,36 @@ end
     @test equispaced == false
 end
 
+@testset "Date and time conversion" begin
+    reference_date = Dates.DateTime(2013, 7, 1, 12)
+    date_one_day = Dates.DateTime(2013, 7, 2, 12)
+    date_one_hour = Dates.DateTime(2013, 7, 1, 13)
+    date_one_min = Dates.DateTime(2013, 7, 1, 12, 1)
+    date_one_sec = Dates.DateTime(2013, 7, 1, 12, 0, 1)
+
+    # Test time_to_date
+    @test Utils.time_to_date(reference_date, 86400.0) == date_one_day
+    @test Utils.time_to_date(reference_date, 3600.0) == date_one_hour
+    @test Utils.time_to_date(reference_date, 60.0) == date_one_min
+    @test Utils.time_to_date(reference_date, 1.0) == date_one_sec
+
+    # Test date_to_time
+    @test Utils.date_to_time(reference_date, date_one_day) == 86400.0
+    @test Utils.date_to_time(reference_date, date_one_hour) == 3600.0
+    @test Utils.date_to_time(reference_date, date_one_min) == 60.0
+    @test Utils.date_to_time(reference_date, date_one_sec) == 1.0
+
+    # Test period_to_seconds_float
+    @test Utils.period_to_seconds_float(Dates.Millisecond(1)) == 0.001
+    @test Utils.period_to_seconds_float(Dates.Second(1)) == 1.0
+    @test Utils.period_to_seconds_float(Dates.Minute(1)) == 60.0
+    @test Utils.period_to_seconds_float(Dates.Hour(1)) == 3600.0
+    @test Utils.period_to_seconds_float(Dates.Day(1)) == 86400.0
+    @test Utils.period_to_seconds_float(Dates.Week(1)) == 604800.0
+    @test Utils.period_to_seconds_float(Dates.Month(1)) == 2.629746e6
+    @test Utils.period_to_seconds_float(Dates.Year(1)) == 3.1556952e7
+end
+
 @testset "data_at_dim_vals" begin
     data = [[1, 2, 3] [4, 5, 6] [7, 8, 9]]
     dim_arr = [2.0, 3.0, 4.0]

From dd0835f9c9b2739b6d9b3322dbdf4a537d8fdbbe Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Mon, 9 Sep 2024 15:03:02 -0700
Subject: [PATCH 5/6] Remove Utils.split_by_season from public API

We do this because we want to show the documentation for
Var.split_by_season and not Utils.split_by_season.
---
 docs/src/api.md | 1 -
 src/Utils.jl    | 7 +------
 2 files changed, 1 insertion(+), 7 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index c11219a..890f1c9 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -72,7 +72,6 @@ Utils.nearest_index
 Utils.kwargs
 Utils.seconds_to_prettystr
 Utils.warp_string
-Utils.split_by_season
 ```
 
 ## Atmos
diff --git a/src/Utils.jl b/src/Utils.jl
index 9f380c1..fd3659e 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -1,12 +1,7 @@
 module Utils
 
 export match_nc_filename,
-    squeeze,
-    nearest_index,
-    kwargs,
-    seconds_to_prettystr,
-    warp_string,
-    split_by_season
+    squeeze, nearest_index, kwargs, seconds_to_prettystr, warp_string
 
 import Dates
 

From c9c0eecd6667796b49505e0ce2cdfc3874492fd9 Mon Sep 17 00:00:00 2001
From: Kevin Phan <98072684+ph-kev@users.noreply.github.com>
Date: Tue, 10 Sep 2024 12:54:01 -0700
Subject: [PATCH 6/6] Add Var.split_by_season

Bug introduced: The documentation show the wrong module name for
Var.split_by_season. It is ClimaAnalysis.Utils.split_by_season instead
of ClimaAnalysis.Var.split_by_season.
---
 NEWS.md          | 40 ++++++++++++++++++++++++++++
 docs/src/api.md  |  1 +
 docs/src/var.md  | 40 ++++++++++++++++++++++++++++
 src/Var.jl       | 69 ++++++++++++++++++++++++++++++++++++++++++++++--
 test/test_Var.jl | 63 +++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 211 insertions(+), 2 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 439c225..e508913 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -164,6 +164,46 @@ julia> long_name(integrated_var) # updated long name to reflect the data being i
 "f integrated over lon (-179.5 to 179.5degrees_east) and integrated over lat (-89.5 to 89.5degrees_north)"
 ```
 
+### Split by season
+`OutputVar`s can be split by seasons using `split_by_season(var)` provided that a start date
+can be found in `var.attributes["start_date"]` and time is a dimension in the `OutputVar`.
+The unit of time is expected to be second. The function `split_by_season(var)` returns a
+vector of four `OutputVar`s with each `OutputVar` corresponding to a season. The months of
+the seasons are March to May, June to August, September to November, and December to
+February. The order of the vector is MAM, JJA, SON, and DJF. If there are no dates found for
+a season, then the `OutputVar` for that season will be an empty `OutputVar`.
+
+```@julia split_by_season
+julia> attribs = Dict("start_date" => "2024-1-1");
+
+julia> time = [0.0, 5_184_000.0, 13_132_800.0]; # correspond to dates 2024-1-1, 2024-3-1, 2024-6-1
+
+julia> dims = OrderedDict(["time" => time]);
+
+julia> dim_attribs = OrderedDict(["time" => Dict("units" => "s")]); # unit is second
+
+julia> data = [1.0, 2.0, 3.0];
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> MAM, JJA, SON, DJF = ClimaAnalysis.split_by_season(var);
+
+julia> ClimaAnalysis.isempty(SON) # empty OutputVar because no dates between September to November
+true
+
+julia> [MAM.dims["time"], JJA.dims["time"], DJF.dims["time"]]
+3-element Vector{Vector{Float64}}:
+ [5.184e6]
+ [1.31328e7]
+ [0.0]
+
+julia> [MAM.data, JJA.data, DJF.data]
+3-element Vector{Vector{Float64}}:
+ [2.0]
+ [3.0]
+ [1.0]
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.
diff --git a/docs/src/api.md b/docs/src/api.md
index 890f1c9..819ba1a 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -59,6 +59,7 @@ Var.convert_units
 Var.integrate_lonlat
 Var.integrate_lat
 Var.integrate_lon
+Var.split_by_season(var::OutputVar)
 ```
 
 ## Utilities
diff --git a/docs/src/var.md b/docs/src/var.md
index 81d2517..64eedc6 100644
--- a/docs/src/var.md
+++ b/docs/src/var.md
@@ -100,3 +100,43 @@ julia> integrated_var.data # approximately 4π (the surface area of a sphere)
 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)"
 ```
+
+## Split by season
+`OutputVar`s can be split by seasons using `split_by_season(var)` provided that a start date
+can be found in `var.attributes["start_date"]` and time is a dimension in the `OutputVar`.
+The unit of time is expected to be second. The function `split_by_season(var)` returns a
+vector of four `OutputVar`s with each `OutputVar` corresponding to a season. The months of
+the seasons are March to May, June to August, September to November, and December to
+February. The order of the vector is MAM, JJA, SON, and DJF. If there are no dates found for
+a season, then the `OutputVar` for that season will be an empty `OutputVar`.
+
+```@julia split_by_season
+julia> attribs = Dict("start_date" => "2024-1-1");
+
+julia> time = [0.0, 5_184_000.0, 13_132_800.0]; # correspond to dates 2024-1-1, 2024-3-1, 2024-6-1
+
+julia> dims = OrderedDict(["time" => time]);
+
+julia> dim_attribs = OrderedDict(["time" => Dict("units" => "s")]); # unit is second
+
+julia> data = [1.0, 2.0, 3.0];
+
+julia> var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data);
+
+julia> MAM, JJA, SON, DJF = ClimaAnalysis.split_by_season(var);
+
+julia> ClimaAnalysis.isempty(SON) # empty OutputVar because no dates between September to November
+true
+
+julia> [MAM.dims["time"], JJA.dims["time"], DJF.dims["time"]]
+3-element Vector{Vector{Float64}}:
+ [5.184e6]
+ [1.31328e7]
+ [0.0]
+
+julia> [MAM.data, JJA.data, DJF.data]
+3-element Vector{Vector{Float64}}:
+ [2.0]
+ [3.0]
+ [1.0]
+```
\ No newline at end of file
diff --git a/src/Var.jl b/src/Var.jl
index 7175019..9546a11 100644
--- a/src/Var.jl
+++ b/src/Var.jl
@@ -1,5 +1,6 @@
 module Var
 
+import Dates
 import NCDatasets
 import OrderedCollections: OrderedDict
 
@@ -8,7 +9,14 @@ import Statistics: mean
 import NaNStatistics: nanmean
 
 import ..Numerics
-import ..Utils: nearest_index, seconds_to_prettystr, squeeze
+import ..Utils:
+    nearest_index,
+    seconds_to_prettystr,
+    squeeze,
+    split_by_season,
+    time_to_date,
+    date_to_time,
+    _data_at_dim_vals
 
 export OutputVar,
     read_var,
@@ -35,7 +43,8 @@ export OutputVar,
     integrate_lonlat,
     integrate_lon,
     integrate_lat,
-    isempty
+    isempty,
+    split_by_season
 
 """
     Representing an output variable
@@ -901,6 +910,62 @@ function _integrate_over_angle(var::OutputVar, integrate_on, angle_dim_name)
     return integrated_var
 end
 
+
+"""
+    split_by_season(var::OutputVar)
+
+Return a vector of four `OutputVar`s split by season.
+
+The months of the seasons are March to May, June to August, September to November, and
+December to February. The order of the vector is MAM, JJA, SON, and DJF. If there are no
+dates found for a season, then the `OutputVar` for that season will be an empty `OutputVar`.
+
+The function will use the start date in `var.attributes["start_date"]`. The unit of time is
+expected to be second. Also, the interpolations will be inaccurate in time intervals
+outside of their respective season for the returned `OutputVar`s.
+"""
+function split_by_season(var::OutputVar)
+    # Check time exists and unit is second
+    has_time(var) || error("Time is not a dimension in var")
+    dim_units(var, time_name(var)) == "s" ||
+        error("Unit for time is not second")
+
+    # Check start date exists
+    haskey(var.attributes, "start_date") ?
+    start_date = Dates.DateTime(var.attributes["start_date"]) :
+    error("Start date is not found in var")
+
+    season_dates = split_by_season(time_to_date.(start_date, times(var)))
+    season_times =
+        (date_to_time.(start_date, season) for season in season_dates)
+
+    # Split data according to seasons
+    season_data = (
+        collect(
+            _data_at_dim_vals(
+                var.data,
+                times(var),
+                var.dim2index[time_name(var)],
+                season_time,
+            ),
+        ) for season_time in season_times
+    )
+
+    # Construct an OutputVar for each season
+    return map(season_times, season_data) do time, data
+        if isempty(time)
+            dims = empty(var.dims)
+            data = similar(var.data, 0)
+            return OutputVar(dims, data)
+        end
+        ret_dims = deepcopy(var.dims)
+        ret_attribs = deepcopy(var.attributes)
+        ret_dim_attribs = deepcopy(var.dim_attributes)
+        ret_dims[time_name(var)] = time
+        OutputVar(ret_attribs, ret_dims, ret_dim_attribs, data)
+    end
+end
+
 """
     overload_binary_op(op)
 
diff --git a/test/test_Var.jl b/test/test_Var.jl
index 7fe224c..0c4acdb 100644
--- a/test/test_Var.jl
+++ b/test/test_Var.jl
@@ -841,3 +841,66 @@ end
         var,
     )
 end
+
+@testset "split_by_season" begin
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    time = [0.0]
+    push!(time, 5_184_000.0) # correspond to 2024-3-1
+    push!(time, 5_184_001.0)
+    push!(time, 13_132_800.0) # correspond to 2024-6-1
+    push!(time, 13_132_802.0)
+    push!(time, 13_132_803.0)
+    data = ones(length(lat), length(time), length(lon))
+    dims = OrderedDict(["lat" => lat, "time" => time, "lon" => lon])
+    attribs = Dict("long_name" => "hi", "start_date" => "2024-1-1")
+    dim_attribs = OrderedDict([
+        "lat" => Dict("units" => "deg"),
+        "time" => Dict("units" => "s"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+
+    MAM, JJA, SON, DJF = ClimaAnalysis.split_by_season(var)
+
+    # Check size of data
+    @test size(MAM.data) == (length(lat), 2, length(lon))
+    @test size(JJA.data) == (length(lat), 3, length(lon))
+    @test size(SON.data) == (0,)
+    @test size(DJF.data) == (length(lat), 1, length(lon))
+
+    # Check times are correct in OutputVars
+    @test MAM.dims["time"] == [5_184_000.0, 5_184_001.0]
+    @test JJA.dims["time"] == [13_132_800.0, 13_132_802.0, 13_132_803.0]
+    @test DJF.dims["time"] == [0.0]
+
+    # Check start date
+    MAM.attributes["start_date"] == "2024-1-1"
+    JJA.attributes["start_date"] == "2024-1-1"
+    DJF.attributes["start_date"] == "2024-1-1"
+
+    # Check empty OutputVar
+    @test isempty(SON)
+
+    # Check error handling
+    attribs_no_start_date = Dict("long_name" => "hi")
+    var =
+        ClimaAnalysis.OutputVar(attribs_no_start_date, dims, dim_attribs, data)
+    @test_throws ErrorException ClimaAnalysis.split_by_season(var)
+
+    dim_attribs_no_sec = OrderedDict([
+        "lat" => Dict("units" => "deg"),
+        "time" => Dict("units" => "min"),
+        "lon" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs_no_sec, data)
+    @test_throws ErrorException ClimaAnalysis.split_by_season(var)
+
+    lon = collect(range(-179.5, 179.5, 360))
+    data = ones(length(lon))
+    dims = OrderedDict(["lon" => lon])
+    attribs = Dict("long_name" => "hi", "start_date" => "2024-1-1")
+    dim_attribs = OrderedDict(["lon" => Dict("units" => "deg")])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+    @test_throws ErrorException ClimaAnalysis.split_by_season(var)
+end