diff --git a/docs/src/userguide.md b/docs/src/userguide.md
index cf451883d..cdc119ae2 100644
--- a/docs/src/userguide.md
+++ b/docs/src/userguide.md
@@ -20,8 +20,6 @@ Any Mimi model is made up of at least one component, so before you construct a m
 
 A component can have any number of parameters and variables. Parameters are data values that will be provided to the component as input, and variables are values that the component will calculate in the `run_timestep` function when the model is run. The index of a parameter or variable determines the number of dimensions that parameter or variable has. They can be scalar values and have no index, such as parameter 'c' in the example below. They can be one-dimensional, such as the variable 'A' and the parameters 'd' and 'f' below. They can be two dimensional such as variable 'B' and parameter 'e' below. Note that any index other than 'time' must be declared at the top of the component, as shown by `regions = Index()` below.
 
-Also note that if a  `Variable` or `Parameter` has `time` as an index, it must be listed as the first index in the definition, eg. `B = Variable(index = [time, regions])` is allowed, but `B = Variable(index = [regions, time])` is not.
-
 The user must define a `run_timestep` function for each component. 
 
 We define a component in the following way:
diff --git a/src/core/build.jl b/src/core/build.jl
index b781054eb..a1bd6b083 100644
--- a/src/core/build.jl
+++ b/src/core/build.jl
@@ -6,19 +6,21 @@ function _instance_datatype(md::ModelDef, def::DatumDef)
     dims = dimensions(def)
     num_dims = dim_count(def)
 
+    ti = get_time_index_position(def)
+
     if num_dims == 0
         T = ScalarModelParameter{dtype}
 
-    elseif dims[1] != :time
+    elseif ti == nothing     # there's no time dimension
         T = Array{dtype, num_dims}
     
     else   
         if isuniform(md)
             first, stepsize = first_and_step(md)
-            T = TimestepArray{FixedTimestep{first, stepsize}, Union{dtype, Missing}, num_dims}
+            T = TimestepArray{FixedTimestep{first, stepsize}, Union{dtype, Missing}, num_dims, ti}
         else
             times = time_labels(md)
-            T = TimestepArray{VariableTimestep{(times...,)}, Union{dtype, Missing}, num_dims}
+            T = TimestepArray{VariableTimestep{(times...,)}, Union{dtype, Missing}, num_dims, ti}
         end
     end
 
diff --git a/src/core/connections.jl b/src/core/connections.jl
index e2b922ab7..4775b3b81 100644
--- a/src/core/connections.jl
+++ b/src/core/connections.jl
@@ -134,16 +134,17 @@ function connect_param!(md::ModelDef,
         if dim_count == 0
             values = backup
         else
+            ti = get_time_index_position(dst_param)
             
             if isuniform(md)
                 # use the first from the comp_def not the ModelDef
                 _, stepsize = first_and_step(md)
-                values = TimestepArray{FixedTimestep{first, stepsize}, T, dim_count}(backup)
+                values = TimestepArray{FixedTimestep{first, stepsize}, T, dim_count, ti}(backup)
             else
                 times = time_labels(md)
                 # use the first from the comp_def 
                 first_index = findfirst(isequal(first), times) 
-                values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, dim_count}(backup)
+                values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, dim_count, ti}(backup)
             end
             
         end
@@ -293,10 +294,11 @@ end
 
 function set_external_param!(md::ModelDef, name::Symbol, value::Union{AbstractArray, AbstractRange, Tuple}; 
                              param_dims::Union{Nothing,Array{Symbol}} = nothing)
-    if param_dims[1] == :time   
+    ti = get_time_index_position(param_dims)
+    if ti != nothing
         value = convert(Array{md.number_type}, value)
         num_dims = length(param_dims)
-        values = get_timestep_array(md, eltype(value), num_dims, value)      
+        values = get_timestep_array(md, eltype(value), num_dims, ti, value)      
     else
         values = value
     end
@@ -426,7 +428,8 @@ function _update_array_param!(md::ModelDef, name, value, update_timesteps, raise
         if param.values isa TimestepArray 
             T = eltype(value)
             N = length(size(value))
-            new_timestep_array = get_timestep_array(md, T, N, value)
+            ti = get_time_index_position(param)
+            new_timestep_array = get_timestep_array(md, T, N, ti, value)
             md.external_params[name] = ArrayModelParameter(new_timestep_array, param.dimensions)
         elseif raise_error
             error("Cannot update timesteps; parameter $name is not a TimestepArray.")
diff --git a/src/core/defcomp.jl b/src/core/defcomp.jl
index 6828c560c..8e8a9be9f 100644
--- a/src/core/defcomp.jl
+++ b/src/core/defcomp.jl
@@ -221,10 +221,6 @@ macro defcomp(comp_name, ex)
                     if !isempty(filter(x -> !(x isa Union{Int,Symbol}), dims))
                         error("Dimensions ($dims) must be defined by a Symbol placeholder or an Int")
                     end
-
-                    if (:time in dims && dims[1] != :time)
-                        error("$elt_type $name: time must be the first dimension ($dims)")
-                    end
                     
                     append!(dimensions, map(Symbol, dims))  # converts, e.g., 4 into Symbol("4")
 
diff --git a/src/core/defs.jl b/src/core/defs.jl
index 0d7701b07..99f11edd3 100644
--- a/src/core/defs.jl
+++ b/src/core/defs.jl
@@ -367,10 +367,12 @@ function set_param!(md::ModelDef, comp_name::Symbol, param_name::Symbol, value,
             value = convert(Array{dtype, num_dims}, value)
         end
 
-        if comp_param_dims[1] == :time
+        ti = get_time_index_position(md, comp_name, param_name)
+
+        if ti != nothing   # there is a time dimension
             T = eltype(value)
 
-            if num_dims == 0
+            if num_dims == 0    
                 values = value
             else
                 # Want to use the first from the comp_def if it has it, if not use ModelDef
@@ -378,12 +380,12 @@ function set_param!(md::ModelDef, comp_name::Symbol, param_name::Symbol, value,
 
                 if isuniform(md)
                     _, stepsize = first_and_step(md)
-                    values = TimestepArray{FixedTimestep{first, stepsize}, T, num_dims}(value)
+                    values = TimestepArray{FixedTimestep{first, stepsize}, T, num_dims, ti}(value)
                 else
                     times = time_labels(md)  
                     #use the first from the comp_def 
                     first_index = findfirst(isequal(first), times)                
-                    values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, num_dims}(value)
+                    values = TimestepArray{VariableTimestep{(times[first_index:end]...,)}, T, num_dims, ti}(value)
                 end 
             end
         else
@@ -748,12 +750,12 @@ function Base.copy(obj::TimestepVector{T_ts, T}) where {T_ts, T}
     return TimestepVector{T_ts, T}(copy(obj.data))
 end
 
-function Base.copy(obj::TimestepMatrix{T_ts, T}) where {T_ts, T}
-    return TimestepMatrix{T_ts, T}(copy(obj.data))
+function Base.copy(obj::TimestepMatrix{T_ts, T, ti}) where {T_ts, T, ti}
+    return TimestepMatrix{T_ts, T, ti}(copy(obj.data))
 end
 
-function Base.copy(obj::TimestepArray{T_ts, T, N}) where {T_ts, T, N}
-    return TimestepArray{T_ts, T, N}(copy(obj.data))
+function Base.copy(obj::TimestepArray{T_ts, T, N, ti}) where {T_ts, T, N, ti}
+    return TimestepArray{T_ts, T, N, ti}(copy(obj.data))
 end
 
 """
diff --git a/src/core/time.jl b/src/core/time.jl
index 109407c7c..03496c2b6 100644
--- a/src/core/time.jl
+++ b/src/core/time.jl
@@ -164,17 +164,23 @@ end
 #
 
 # Get a timestep array of type T with N dimensions. Time labels will match those from the time dimension in md
-function get_timestep_array(md::ModelDef, T, N, value)
+function get_timestep_array(md::ModelDef, T, N, ti, value)
 
 	if isuniform(md)
         first, stepsize = first_and_step(md)
-        return TimestepArray{FixedTimestep{first, stepsize}, T, N}(value)
+        return TimestepArray{FixedTimestep{first, stepsize}, T, N, ti}(value)
     else
         TIMES = (time_labels(md)...,)
-        return TimestepArray{VariableTimestep{TIMES}, T, N}(value)
+        return TimestepArray{VariableTimestep{TIMES}, T, N, ti}(value)
     end
 end
 
+# Return the index position of the time dimension in the datumdef or parameter. If there is no time dimension, return nothing
+get_time_index_position(dims::Union{Nothing, Array{Symbol}}) = findfirst(isequal(:time), dims)
+get_time_index_position(datumdef::DatumDef) = get_time_index_position(datumdef.dimensions)
+get_time_index_position(param::ArrayModelParameter) = get_time_index_position(param.dimensions)
+get_time_index_position(md::ModelDef, comp_name::Symbol, datum_name::Symbol) = get_time_index_position(dimensions(compdef(md, comp_name), datum_name))
+
 const AnyIndex = Union{Int, Vector{Int}, Tuple, Colon, OrdinalRange}
 
 # TBD: can it be reduced to this?
@@ -227,8 +233,8 @@ function Base.getindex(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, ts::F
 	_missing_data_check(data)
 end
 
-function Base.getindex(v::TimestepVector{VariableTimestep{D_FIRST}, T}, ts::VariableTimestep{T_FIRST}) where {T, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
+function Base.getindex(v::TimestepVector{VariableTimestep{D_TIMES}, T}, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
 	data = v.data[t]
 	_missing_data_check(data)
 end
@@ -257,8 +263,8 @@ function Base.setindex!(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, val,
 	setindex!(v.data, val, t)
 end
 
-function Base.setindex!(v::TimestepVector{VariableTimestep{D_FIRST}, T}, val, ts::VariableTimestep{T_FIRST}) where {T, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
+function Base.setindex!(v::TimestepVector{VariableTimestep{D_TIMES}, T}, val, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
 	setindex!(v.data, val, t)
 end
 
@@ -283,74 +289,112 @@ Base.lastindex(v::TimestepVector) = length(v)
 # 3c. TimestepMatrix
 #
 
-function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, ts::FixedTimestep{FIRST, STEP, LAST}, i::AnyIndex) where {T, FIRST, STEP, LAST} 
-	data = mat.data[ts.t, i]
+function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 1}, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST} 
+	data = mat.data[ts.t, idx]
 	_missing_data_check(data)
 end
 
-function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, ts::VariableTimestep{TIMES}, i::AnyIndex) where {T, TIMES}
-	data = mat.data[ts.t, i]
+function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 1}, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
+	data = mat.data[ts.t, idx]
 	_missing_data_check(data)
 end
 
-function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T}, ts::FixedTimestep{T_FIRST, STEP, LAST}, i::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST} 
+function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 1}, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST} 
 	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
-	data = mat.data[t, i]
+	data = mat.data[t, idx]
 	_missing_data_check(data)
 end
 
-function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_FIRST}, T}, ts::VariableTimestep{T_FIRST}, i::AnyIndex) where {T, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
-	data = mat.data[t, i]
+function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 1}, ts::VariableTimestep{T_TIMES}, idx::AnyIndex) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
+	data = mat.data[t, idx]
 	_missing_data_check(data)
 end
 
-# int indexing version supports old-style components and internal functions, not
-# part of the public API
+function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 2}, idx::AnyIndex, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST} 
+	data = mat.data[idx, ts.t]
+	_missing_data_check(data)
+end
 
-function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP}
-	return mat.data[idx1, idx2]
+function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 2}, idx::AnyIndex, ts::VariableTimestep{TIMES}) where {T, TIMES}
+	data = mat.data[ts.t, idx, ts.t]
+	_missing_data_check(data)
 end
 
-function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES}
-	return mat.data[idx1, idx2]
+function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 2}, idx::AnyIndex, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST} 
+	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
+	data = mat.data[idx, ts.t]
+	_missing_data_check(data)
+end
+
+function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 2}, idx::AnyIndex, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
+	data = mat.data[idx, ts.t]
+	_missing_data_check(data)
 end
 
-function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST} 
+
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 1}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST} 
 	setindex!(mat.data, val, ts.t, idx)
 end
 
-function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 1}, val, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
 	setindex!(mat.data, val, ts.t, idx)
 end
 
-function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST} 
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 1}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST} 
 	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
 	setindex!(mat.data, val, t, idx)
 end
 
-function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_FIRST}, T}, val, ts::VariableTimestep{T_FIRST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 1}, val, ts::VariableTimestep{T_TIMES}, idx::AnyIndex) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
 	setindex!(mat.data, val, t, idx)
 end
 
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, 2}, val, idx::AnyIndex, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST} 
+	setindex!(mat.data, val, idx, ts.t)
+end
+
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, 2}, val, idx::AnyIndex, ts::VariableTimestep{TIMES}) where {T, TIMES}
+	setindex!(mat.data, val, idx, ts.t)
+end
+
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T, 2}, val, idx::AnyIndex, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST} 
+	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
+	setindex!(mat.data, val, idx, t)
+end
+
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_TIMES}, T, 2}, val, idx::AnyIndex, ts::VariableTimestep{T_TIMES}) where {T, D_TIMES, T_TIMES}
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1
+	setindex!(mat.data, val, idx, t)
+end
+
 # int indexing version supports old-style components and internal functions, not
 # part of the public API
 
-function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, idx1::Int, idx2::Int) where {T, FIRST, STEP}
+function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP, ti}
+	return mat.data[idx1, idx2]
+end
+
+function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES, ti}
+	return mat.data[idx1, idx2]
+end
+
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, val, idx1::Int, idx2::Int) where {T, FIRST, STEP, ti}
 	setindex!(mat.data, val, idx1, idx2)
 end
 
-function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP}
-	mat.data[idx1,idx2] .= val
+function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T, ti}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP, ti}
+	mat.data[idx1, idx2] .= val
 end
 
-function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, idx1::Int, idx2::Int) where {T, TIMES}
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, val, idx1::Int, idx2::Int) where {T, TIMES, ti}
 	setindex!(mat.data, val, idx1, idx2)
 end
 
-function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES}
-	mat.data[idx1,idx2] .= val
+function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T, ti}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES, ti}
+	mat.data[idx1, idx2] .= val
 end
 
 #
@@ -368,74 +412,81 @@ Base.size(obj::TimestepArray) = size(obj.data)
 
 Base.size(obj::TimestepArray, i::Int) = size(obj.data, i)
 
-Base.ndims(obj::TimestepArray{T_ts, T, N}) where {T_ts,T, N} = N
-
-Base.eltype(obj::TimestepArray{T_ts, T, N}) where {T_ts,T, N} = T
+Base.ndims(obj::TimestepArray{T_ts, T, N, ti}) where {T_ts, T, N, ti} = N
 
-first_period(obj::TimestepArray{FixedTimestep{FIRST,STEP}, T, N}) where {FIRST, STEP, T, N} = FIRST
-first_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES, T, N} = TIMES[1]
+Base.eltype(obj::TimestepArray{T_ts, T, N, ti}) where {T_ts, T, N, ti} = T
 
-last_period(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}) where {FIRST, STEP,T, N} = (FIRST + (size(obj, 1) - 1) * STEP)
-last_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES,T, N} = TIMES[end]
+first_period(obj::TimestepArray{FixedTimestep{FIRST,STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = FIRST
+first_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = TIMES[1]
 
-time_labels(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}) where {FIRST, STEP, T, N} = collect(FIRST:STEP:(FIRST + (size(obj, 1) - 1) * STEP))
-time_labels(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES, T, N} = collect(TIMES)
+last_period(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = (FIRST + (size(obj, 1) - 1) * STEP)
+last_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = TIMES[end]
 
-function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, ts::FixedTimestep{FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, FIRST, STEP, LAST}
-	return arr.data[ts.t, idxs...]
-end
+time_labels(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}) where {FIRST, STEP, T, N, ti} = collect(FIRST:STEP:(FIRST + (size(obj, 1) - 1) * STEP))
+time_labels(obj::TimestepArray{VariableTimestep{TIMES}, T, N, ti}) where {TIMES, T, N, ti} = collect(TIMES)
 
-function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, ts::VariableTimestep{TIMES}, idxs::AnyIndex...) where {T, N, TIMES}
-	return arr.data[ts.t, idxs...]
-end
+split_indices(idxs, ti) = idxs[1:ti - 1], idxs[ti], idxs[ti + 1:end]
 
-function Base.getindex(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, ts::FixedTimestep{T_FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
-	t = Int(ts.t + (FIRST - TIMES[1]) / STEP)					
-	return arr.data[t, idxs...]
+function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, idxs::Union{FixedTimestep{FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, FIRST, STEP, LAST}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	return arr.data[idxs1..., ts.t, idxs2...]
 end
 
-function Base.getindex(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, ts::VariableTimestep{T_FIRST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1	
-	return arr.data[t, idxs...]
+function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, idxs::Union{VariableTimestep{TIMES}, AnyIndex}...) where {T, N, ti, TIMES}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	return arr.data[idxs1..., ts.t, idxs2...]
 end
 
-# int indexing version supports old-style components and internal functions, not
-# part of the public API; first index is Int or Range, rather than a Timestep
-
-function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, FIRST, STEP}
-	return arr.data[idx1, idx2, idxs...]
+function Base.getindex(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, idxs::Union{FixedTimestep{T_FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	t = Int(ts.t + (FIRST - TIMES[1]) / STEP)					
+	return arr.data[idxs1..., t, idxs2...]
 end
 
-function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, TIMES}
-	return arr.data[idx1, idx2, idxs...]
+function Base.getindex(arr::TimestepArray{VariableTimestep{D_TIMES}, T, N, ti}, idxs::Union{VariableTimestep{T_TIMES}, AnyIndex}...) where {T, N, ti, D_TIMES, T_TIMES}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	t = ts.t + findfirst(isequal(T_TIMES[1]), D_TIMES) - 1	
+	return arr.data[idxs1..., t, idxs2...]
 end
 
-function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, FIRST, STEP, LAST}
-	setindex!(arr.data, val, ts.t, idxs...)
+function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, val, idxs::Union{FixedTimestep{FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, FIRST, STEP, LAST}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	setindex!(arr.data, val, idxs1..., ts.t, idxs2...)
 end
 
-function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, val, ts::VariableTimestep{TIMES}, idxs::AnyIndex...) where {T, N, TIMES}
-	setindex!(arr.data, val, ts.t, idxs...)
+function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, val, idxs::Union{VariableTimestep{TIMES}, AnyIndex}...) where {T, N, ti, TIMES}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	setindex!(arr.data, val, idxs1..., ts.t, idxs2...)
 end
 
-function Base.setindex!(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
+function Base.setindex!(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, val, idxs::Union{FixedTimestep{T_FIRST, STEP, LAST}, AnyIndex}...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
 	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1	
-	setindex!(arr.data, val, t, idxs...)
+	setindex!(arr.data, val, idxs1..., t, idxs2...)
 end
 
-function Base.setindex!(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, val, ts::VariableTimestep{T_FIRST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST}
-	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1	
-	setindex!(arr.data, val, t, idxs...)
+function Base.setindex!(arr::TimestepArray{VariableTimestep{D_TIMES}, T, N, ti}, val, idxs::Union{VariableTimestep{T_TIMES}, AnyIndex}...) where {T, N, ti, D_TIMES, T_TIMES}
+	idxs1, ts, idxs2 = split_indices(idxs, ti)
+	t = ts.t + findfirst(isequal(T_FIRST[1]), T_TIMES) - 1	
+	setindex!(arr.data, val, idxs1..., t, idxs2...)
 end
 
 # int indexing version supports old-style components and internal functions, not
 # part of the public API; first index is Int or Range, rather than a Timestep
 
-function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, FIRST, STEP}
+function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, ti, FIRST, STEP}
+	return arr.data[idx1, idx2, idxs...]
+end
+
+function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, ti, TIMES}
+	return arr.data[idx1, idx2, idxs...]
+end
+
+function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, ti, FIRST, STEP}
 	setindex!(arr.data, val, idx1, idx2, idxs...)
 end
 
-function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, TIMES}
+function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, ti, TIMES}
 	setindex!(arr.data, val, idx1, idx2, idxs...)
 end
 
@@ -444,7 +495,7 @@ end
 
 Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
 """
-function hasvalue(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, N, FIRST, STEP, LAST}
+function hasvalue(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N, ti}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, N, ti, FIRST, STEP, LAST}
 	return 1 <= ts.t <= size(arr, 1)	
 end
 
@@ -453,15 +504,15 @@ end
 
 Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
 """
-function hasvalue(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, ts::VariableTimestep{TIMES}) where {T, N, TIMES}
+function hasvalue(arr::TimestepArray{VariableTimestep{TIMES}, T, N, ti}, ts::VariableTimestep{TIMES}) where {T, N, ti, TIMES}
 	return 1 <= ts.t <= size(arr, 1)	
 end
 
-function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
+function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
 	return D_FIRST <= gettime(ts) <= last_period(arr)
 end
 
-function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, ts::VariableTimestep{T_FIRST}) where {T, N, T_FIRST, D_FIRST}
+function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N, ti}, ts::VariableTimestep{T_FIRST}) where {T, N, ti, T_FIRST, D_FIRST}
 	return D_FIRST[1] <= gettime(ts) <= last_period(arr)	
 end
 
@@ -471,9 +522,9 @@ end
 Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
 indices `idxs`. Used when Array and Timestep have different FIRST, validating all dimensions.
 """
-function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, 
+function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N, ti}, 
 	ts::FixedTimestep{T_FIRST, STEP, LAST}, 
-	idxs::Int...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
+	idxs::Int...) where {T, N, ti, D_FIRST, T_FIRST, STEP, LAST}
 	return D_FIRST <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
 end
 
@@ -483,9 +534,9 @@ end
 Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
 indices `idxs`. Used when Array and Timestep different TIMES, validating all dimensions.
 """
-function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, 
+function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N, ti}, 
 	ts::VariableTimestep{T_FIRST}, 
-	idxs::Int...) where {T, N, D_FIRST, T_FIRST}
+	idxs::Int...) where {T, N, ti, D_FIRST, T_FIRST}
 
 	return D_FIRST[1] <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
 end
diff --git a/src/core/types.jl b/src/core/types.jl
index 8ed4d68c3..28b8832f1 100644
--- a/src/core/types.jl
+++ b/src/core/types.jl
@@ -34,14 +34,14 @@ mutable struct Clock{T <: AbstractTimestep}
     end
 end
 
-mutable struct TimestepArray{T_TS <: AbstractTimestep, T, N}
+mutable struct TimestepArray{T_TS <: AbstractTimestep, T, N, ti}
 	data::Array{T, N}
 
-    function TimestepArray{T_TS, T, N}(d::Array{T, N}) where {T_TS, T, N}
+    function TimestepArray{T_TS, T, N, ti}(d::Array{T, N}) where {T_TS, T, N, ti}
 		return new(d)
 	end
 
-    function TimestepArray{T_TS, T, N}(lengths::Int...) where {T_TS, T, N}
+    function TimestepArray{T_TS, T, N, ti}(lengths::Int...) where {T_TS, T, N, ti}
 		return new(Array{T, N}(undef, lengths...))
 	end
 end
@@ -49,8 +49,8 @@ end
 # Since these are the most common cases, we define methods (in time.jl)
 # specific to these type aliases, avoiding some of the inefficiencies
 # associated with an arbitrary number of dimensions.
-const TimestepMatrix{T_TS, T} = TimestepArray{T_TS, T, 2}
-const TimestepVector{T_TS, T} = TimestepArray{T_TS, T, 1}
+const TimestepMatrix{T_TS, T, ti} = TimestepArray{T_TS, T, 2, ti}
+const TimestepVector{T_TS, T} = TimestepArray{T_TS, T, 1, 1}
 
 #
 # 2. Dimensions
diff --git a/src/explorer/buildspecs.jl b/src/explorer/buildspecs.jl
index fadf21797..5dfe0dbcb 100644
--- a/src/explorer/buildspecs.jl
+++ b/src/explorer/buildspecs.jl
@@ -22,15 +22,23 @@ function _spec_for_item(m::Model, comp_name::Symbol, item_name::Symbol; interact
     else
         name = "$comp_name : $item_name"          # the name is needed for the list label
         df = getdataframe(m, comp_name, item_name)
-
         dffields = map(string, names(df))         # convert to string once before creating specs
 
         # check if there are too many dimensions to map and if so, warn
         if length(dffields) > 3
             error()
-                
-        # a 'time' field necessitates a line plot  
-        elseif dffields[1] == "time"
+               
+        # a 'time' field necessitates a line plot
+        elseif "time" in dffields
+
+            # need to reorder the df to have 'time' as the first dimension
+            ti = findfirst(isequal("time"), dffields)
+            if ti != 1    
+                fields1, fields2 = dffields[1:ti-1], dffields[ti+1:end]
+                dffields = ["time", fields1..., fields2...]
+                df = df[:, [Symbol(name) for name in dffields]]
+            end
+
             if length(dffields) > 2
                 spec = createspec_multilineplot(name, df, dffields, interactive=interactive)
             else
@@ -145,6 +153,7 @@ function createspec_lineplot_interactive(name, df, dffields)
     return spec
 end
 
+
 function createspec_lineplot_static(name, df, dffields)
     datapart = getdatapart(df, dffields, :line) #returns JSONtext type 
     spec = Dict(
diff --git a/test/test_parametertypes.jl b/test/test_parametertypes.jl
index c8a20a112..52a83434e 100644
--- a/test/test_parametertypes.jl
+++ b/test/test_parametertypes.jl
@@ -78,7 +78,7 @@ extpars = external_params(m)
 @test isa(extpars[:e], ArrayModelParameter)
 @test isa(extpars[:f], ScalarModelParameter) # note that :f is stored as a scalar parameter even though its values are an array
 
-@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, numtype}
+@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, numtype, 1}
 @test typeof(extpars[:b].values) == TimestepVector{FixedTimestep{2000, 1}, numtype}
 @test typeof(extpars[:c].values) == Array{numtype, 1}
 @test typeof(extpars[:d].value) == numtype
@@ -101,7 +101,7 @@ update_param!(m, :d, 5) # should work, will convert to float
 update_param!(m, :e, [4,5,6,7])
 
 @test length(extpars) == 8
-@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, numtype}
+@test typeof(extpars[:a].values) == TimestepMatrix{FixedTimestep{2000, 1}, numtype, 1}
 @test typeof(extpars[:d].value) == numtype
 @test typeof(extpars[:e].values) == Array{numtype, 1}
 
diff --git a/test/test_timesteparrays.jl b/test/test_timesteparrays.jl
index fd1e13c00..71d6bc4a7 100644
--- a/test/test_timesteparrays.jl
+++ b/test/test_timesteparrays.jl
@@ -102,7 +102,7 @@ years = Tuple(2000:1:2003)
 #3a.  test constructor (with both matching years 
 # and mismatched years)
 
-y = TimestepMatrix{FixedTimestep{2000, 1}, Int}(a[:,1:2])
+y = TimestepMatrix{FixedTimestep{2000, 1}, Int, 1}(a[:,1:2])
 
 #3b.  test hasvalue, getindex, and setindex! (with both matching years and
 # mismatched years)
@@ -135,7 +135,7 @@ y[:,:] = 11
 @test all([y[1,j] == 11 for j in 1:2])    
 
 #3c.  interval wider than 1
-z = TimestepMatrix{FixedTimestep{2000, 2}, Int}(a[:,3:4])
+z = TimestepMatrix{FixedTimestep{2000, 2}, Int, 1}(a[:,3:4])
 t = FixedTimestep{1980, 2, 3000}(11)
 
 @test z[t,1] == 9
@@ -151,7 +151,7 @@ t2 = next_timestep(t)
 #------------------------------------------------------------------------------
 
 years = Tuple([2000:5:2005; 2015:10:2025])
-y = TimestepMatrix{VariableTimestep{years}, Int}(a[:,1:2])
+y = TimestepMatrix{VariableTimestep{years}, Int, 1}(a[:,1:2])
 
 #4a.  test hasvalue, getindex, setindex!, and lastindex (with both matching years and
 # mismatched years)
@@ -192,9 +192,9 @@ x_years = Tuple(2000:5:2015) #fixed
 y_years = Tuple([2000:5:2005; 2015:10:2025]) #variable
 
 x_vec = TimestepVector{FixedTimestep{2000, 5}, Int}(a[:,3]) 
-x_mat = TimestepMatrix{FixedTimestep{2000, 5}, Int}(a[:,1:2])
+x_mat = TimestepMatrix{FixedTimestep{2000, 5}, Int, 1}(a[:,1:2])
 y_vec = TimestepVector{VariableTimestep{y_years}, Int}(a[:,3]) 
-y_mat = TimestepMatrix{VariableTimestep{y_years}, Int}(a[:,1:2])
+y_mat = TimestepMatrix{VariableTimestep{y_years}, Int, 1}(a[:,1:2])
 
 @test first_period(x_vec) == first_period(x_mat) == x_years[1] 
 @test first_period(y_vec) == first_period(y_mat) == y_years[1]
@@ -250,4 +250,55 @@ x = Mimi.TimestepVector{Mimi.FixedTimestep{2005,10}, Float64}(zeros(10))
 x[:] .= 10
 @test all(x.data .== 10)
 
-end #module
+#------------------------------------------------------------------------------
+# 7. Test TimestepArrays with time not as the first dimension
+#------------------------------------------------------------------------------
+
+@defcomp gdp begin
+    growth = Parameter(index=[regions, foo, time, 2])   # test that time is not first but not last
+    gdp = Variable(index=[regions, foo, time, 2])
+    gdp0 = Parameter(index=[regions, foo, 2])
+
+    pgrowth = Parameter(index=[regions, 3, time])       # test time as last
+    pop = Variable(index=[regions, 3, time])
+
+    mat = Parameter(index=[regions, time])              # test time as last for a matrix
+    mat2 = Variable(index=[regions, time])
+
+    function run_timestep(p, v, d, ts)
+        if is_first(ts)
+            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* p.gdp0
+            v.pop[:, :, ts] = zeros(2, 3)
+        else
+            v.gdp[:, :, ts, :] = (1 .+ p.growth[:, :, ts, :]) .* v.gdp[:, :, ts-1, :]
+            v.pop[:, :, ts] = v.pop[:, :, ts-1] .+ p.pgrowth[:, :, ts]
+        end
+        v.mat2[:, ts] = p.mat[:, ts]
+    end
+end
+
+time_index = 2000:2100
+regions = ["OECD","non-OECD"]
+nsteps=length(time_index)
+
+m = Model()
+set_dimension!(m, :time, time_index)
+set_dimension!(m, :regions, regions)
+set_dimension!(m, :foo, 3)
+add_comp!(m, gdp)
+set_param!(m, :gdp, :gdp0, [3; 7] .* ones(length(regions), 3, 2))
+set_param!(m, :gdp, :growth, [0.02; 0.03] .* ones(length(regions), 3, nsteps, 2))
+set_leftover_params!(m, Dict{String, Any}([
+    "pgrowth" => ones(length(regions), 3, nsteps),
+    "mat" => rand(length(regions), nsteps)
+]))
+run(m)
+explore(m)
+
+@test size(m[:gdp, :gdp]) == (length(regions), 3, length(time_index), 2)
+
+@test all(!ismissing, m[:gdp, :gdp])
+@test all(!ismissing, m[:gdp, :pop])
+@test all(!ismissing, m[:gdp, :mat2])
+
+end #module
\ No newline at end of file
diff --git a/test/test_timesteps.jl b/test/test_timesteps.jl
index 9eaed5d58..ab52dffa7 100644
--- a/test/test_timesteps.jl
+++ b/test/test_timesteps.jl
@@ -149,8 +149,8 @@ set_dimension!(m, :time, 2000:2009)
 vector = ones(5)
 matrix = ones(3,2)
 
-t_vector= get_timestep_array(m.md, Float64, 1, vector)
-t_matrix = get_timestep_array(m.md, Float64, 2, matrix)
+t_vector= get_timestep_array(m.md, Float64, 1, 1, vector)
+t_matrix = get_timestep_array(m.md, Float64, 2, 1, matrix)
 
 @test typeof(t_vector) <: TimestepVector
 @test typeof(t_matrix) <: TimestepMatrix
@@ -162,8 +162,8 @@ t_matrix = get_timestep_array(m.md, Float64, 2, matrix)
     set_dimension!(m, :time, [2000:1:2004; 2005:2:2016])
 )
 
-t_vector= get_timestep_array(m.md, Float64, 1, vector)
-t_matrix = get_timestep_array(m.md, Float64, 2, matrix)
+t_vector= get_timestep_array(m.md, Float64, 1, 1, vector)
+t_matrix = get_timestep_array(m.md, Float64, 2, 2, matrix)
 
 @test typeof(t_vector) <: TimestepVector
 @test typeof(t_matrix) <: TimestepMatrix