Merge #85

85: Convert to 1-based indexing r=charleskawczynski a=charleskawczynski

This is work that I started this morning, just wanted to push it up to see if/how things run.

Co-authored-by: Charles Kawczynski <[email protected]>

integration_tests/utils/Cases.jl

@@ -1010,14 +1010,14 @@ function initialize_forcing(self::CasesBase{TRMM_LBA}, Gr::Grid, Ref::ReferenceS
     A = off_arr(A')
 
     self.rad = A # store matrix in self
-    ind1 = Int(trunc(10.0 / 600.0))
-    ind2 = Int(ceil(10.0 / 600.0)) - 1
+    ind1 = Int(trunc(10.0 / 600.0)) + 1
+    ind2 = Int(ceil(10.0 / 600.0))
     @inbounds for k in face_indicies(Gr)
         if 10 % 600.0 == 0
             self.Fo.dTdt[k] = self.rad[ind1, k]
         else
             self.Fo.dTdt[k] =
-                (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2 + 1] - self.rad_time[ind1 + 1]) * (10.0) +
+                (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2] - self.rad_time[ind1]) * (10.0) +
                 self.rad[ind1, k]
         end
     end
@@ -1047,24 +1047,25 @@ end
 function update_forcing(self::CasesBase{TRMM_LBA}, GMV::GridMeanVariables, TS::TimeStepping)
 
     Gr = self.Fo.Gr
-    ind2 = Int(ceil(TS.t / 600.0))
-    ind1 = Int(trunc(TS.t / 600.0))
+    ind2 = Int(ceil(TS.t / 600.0)) + 1
+    ind1 = Int(trunc(TS.t / 600.0)) + 1
     @inbounds for k in face_indicies(Gr)
         if Gr.z_half[k] >= 22699.48
             self.Fo.dTdt[k] = 0.0
         else
             if TS.t < 600.0 # first 10 min use the radiative forcing of t=10min (as in the paper)
-                self.Fo.dTdt[k] = self.rad[0, k]
-            elseif TS.t < 21600.0 && ind2 < 36
+                self.Fo.dTdt[k] = self.rad[1, k]
+            elseif TS.t < 21600.0 && ind2 < 37
                 if TS.t % 600.0 == 0
                     self.Fo.dTdt[k] = self.rad[ind1, k]
                 else
                     self.Fo.dTdt[k] =
-                        (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2] - self.rad_time[ind1]) *
-                        (TS.t - self.rad_time[ind1]) + self.rad[ind1, k]
+                        (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2 - 1] - self.rad_time[ind1 - 1]) *
+                        (TS.t - self.rad_time[ind1 - 1]) + self.rad[ind1, k]
                 end
             else
-                self.Fo.dTdt[k] = self.rad[35, k]
+                # TODO: remove hard-coded index
+                self.Fo.dTdt[k] = self.rad[36, k]
             end
         end
     end
@@ -1168,8 +1169,8 @@ function update_surface(self::CasesBase{ARM_SGP}, GMV::GridMeanVariables, TS::Ti
     t_Sur_in = off_arr([0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5]) .* 3600 #LES time is in sec
     SH = off_arr([-30.0, 90.0, 140.0, 140.0, 100.0, -10, -10]) # W/m^2
     LH = off_arr([5.0, 250.0, 450.0, 500.0, 420.0, 180.0, 0.0]) # W/m^2
-    self.Sur.shf = pyinterp(off_arr([TS.t]), t_Sur_in, SH)[0]
-    self.Sur.lhf = pyinterp(off_arr([TS.t]), t_Sur_in, LH)[0]
+    self.Sur.shf = pyinterp(off_arr([TS.t]), t_Sur_in, SH)[1]
+    self.Sur.lhf = pyinterp(off_arr([TS.t]), t_Sur_in, LH)[1]
     # if fluxes vanish bflux vanish and wstar and obukov length are NaNs
     ## CK +++ I commented out the lines below as I don"t think this is how we want to fix things!
     # if self.Sur.shf < 1.0
@@ -1188,8 +1189,8 @@ function update_forcing(self::CasesBase{ARM_SGP}, GMV::GridMeanVariables, TS::Ti
     AT_in = off_arr([0.0, 0.0, 0.0, -0.08, -0.016, -0.016]) ./ 3600.0 # Advective forcing for theta [K/h] converted to [K/sec]
     RT_in = off_arr([-0.125, 0.0, 0.0, 0.0, 0.0, -0.1]) ./ 3600.0  # Radiative forcing for theta [K/h] converted to [K/sec]
     Rqt_in = off_arr([0.08, 0.02, 0.04, -0.1, -0.16, -0.3]) ./ 1000.0 ./ 3600.0 # Radiative forcing for qt converted to [kg/kg/sec]
-    dTdt = pyinterp(off_arr([TS.t]), t_in, AT_in)[0] + pyinterp(off_arr([TS.t]), t_in, RT_in)[0]
-    dqtdt = pyinterp(off_arr([TS.t]), t_in, Rqt_in)[0]
+    dTdt = pyinterp(off_arr([TS.t]), t_in, AT_in)[1] + pyinterp(off_arr([TS.t]), t_in, RT_in)[1]
+    dqtdt = pyinterp(off_arr([TS.t]), t_in, Rqt_in)[1]
     Gr = self.Fo.Gr
     @inbounds for k in center_indicies(Gr)
         if Gr.z_half[k] <= 1000.0

src/EDMF_Environment.jl

@@ -1,27 +1,30 @@
 function set_bcs(self::EnvironmentVariable, Gr::Grid)
-    start_low = Gr.gw - 1
-    start_high = Gr.nzg - Gr.gw - 1
+    start_low = Gr.gw
+    start_high = Gr.nzg - Gr.gw
 
     if self.name == "w"
         self.values[start_high] = 0.0
         self.values[start_low] = 0.0
-        @inbounds for k in xrange(1, Gr.gw)
+        @inbounds for kk in xrange(1, Gr.gw)
+            k = kk - 1
             self.values[start_high + k] = -self.values[start_high - k]
             self.values[start_low - k] = -self.values[start_low + k]
         end
     else
-        @inbounds for k in xrange(Gr.gw)
+        @inbounds for kk in xrange(Gr.gw)
+            k = kk - 1
             self.values[start_high + k + 1] = self.values[start_high - k]
             self.values[start_low - k] = self.values[start_low + 1 + k]
         end
     end
 end
 
 function set_bcs(self::EnvironmentVariable_2m, Gr::Grid)
-    start_low = Gr.gw - 1
-    start_high = Gr.nzg - Gr.gw - 1
+    start_low = Gr.gw
+    start_high = Gr.nzg - Gr.gw
 
-    @inbounds for k in xrange(Gr.gw)
+    @inbounds for kk in xrange(Gr.gw)
+        k = kk - 1
         self.values[start_high + k + 1] = self.values[start_high - k]
         self.values[start_low - k] = self.values[start_low + 1 + k]
     end
@@ -270,7 +273,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari
                 corr = fmax(fmin(EnvVar.HQTcov.values[k] / fmax(sd_h * sd_q, 1e-13), 1.0), -1.0)
 
                 # limit sd_q to prevent negative qt_hat
-                sd_q_lim = (1e-10 - EnvVar.QT.values[k]) / (sqrt2 * abscissas[0])
+                sd_q_lim = (1e-10 - EnvVar.QT.values[k]) / (sqrt2 * abscissas[1])
                 # walking backwards to assure your q_t will not be smaller than 1e-10
                 # TODO - check
                 # TODO - change 1e-13 and 1e-10 to some epislon

src/EDMF_Updrafts.jl

@@ -1,20 +1,22 @@
 function set_bcs(self::UpdraftVariable, Gr::Grid)
-    start_low = Gr.gw - 1
-    start_high = Gr.nzg - Gr.gw - 1
+    start_low = Gr.gw
+    start_high = Gr.nzg - Gr.gw
 
     n_updrafts = size(self.values)[1]
 
     if self.name == "w"
         @inbounds for i in xrange(n_updrafts)
             self.values[i, start_high] = 0.0
             self.values[i, start_low] = 0.0
-            @inbounds for k in xrange(1, Gr.gw)
+            @inbounds for kk in xrange(1, Gr.gw)
+                k = kk - 1
                 self.values[i, start_high + k] = -self.values[i, start_high - k]
                 self.values[i, start_low - k] = -self.values[i, start_low + k]
             end
         end
     else
-        @inbounds for k in xrange(Gr.gw)
+        @inbounds for kk in xrange(Gr.gw)
+            k = kk - 1
             @inbounds for i in xrange(n_updrafts)
                 self.values[i, start_high + k + 1] = self.values[i, start_high - k]
                 self.values[i, start_low - k] = self.values[i, start_low + 1 + k]

src/Grid.jl

@@ -24,31 +24,31 @@ struct Grid{A1, FT}
         nz = namelist["grid"]["nz"]
         nzg = nz + 2 * gw
 
-        cinterior = gw:((nzg - gw) - 1)
-        finterior = gw:((nzg - gw) - 1)
+        cinterior = (gw + 1):(nzg - gw)
+        finterior = (gw + 1):(nzg - gw)
 
         # TODO: make cell centers and cell faces different sizes
         z_half = pyzeros(nz + 2 * gw)
         z = pyzeros(nz + 2 * gw)
-        count = 0
-        @inbounds for i in xrange(-gw, nz + gw)
-            z[count] = (i + 1) * dz
-            z_half[count] = (i + 0.5) * dz
+        count = 1
+        @inbounds for i in range(1 - gw, stop = nz + gw)
+            z[count] = i * dz
+            z_half[count] = (i - 0.5) * dz
             count += 1
         end
-        zmin = z[gw - 1]
-        zmax = z[nzg - gw - 1]
+        zmin = z[gw]
+        zmax = z[nzg - gw]
         A1 = typeof(z)
         FT = typeof(zmax)
         return new{A1, FT}(zmin, zmax, cinterior, finterior, dz, dzi, gw, nz, nzg, z_half, z)
     end
 end
 
 # Index of the first interior cell above the surface
-kc_surface(grid::Grid) = grid.gw
-kf_surface(grid::Grid) = grid.gw - 1
-kc_top_of_atmos(grid::Grid) = grid.nzg - grid.gw - 1
-kf_top_of_atmos(grid::Grid) = grid.nzg - grid.gw - 1
+kc_surface(grid::Grid) = grid.gw + 1
+kf_surface(grid::Grid) = grid.gw
+kc_top_of_atmos(grid::Grid) = grid.nzg - grid.gw
+kf_top_of_atmos(grid::Grid) = grid.nzg - grid.gw
 
 is_surface_center(grid::Grid, k::Int) = k == kc_surface(grid)
 is_toa_center(grid::Grid, k::Int) = k == kc_top_of_atmos(grid)

src/ReferenceState.jl

@@ -68,7 +68,6 @@ function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats)
     # Perform the integration
     # TODO: replace with OrdinaryDiffEq
 
-    grid
     z_span = (grid.zmin, grid.zmax)
     @show z_span
     prob = ODEProblem(rhs, logp, z_span)
@@ -83,12 +82,12 @@ function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats)
 
     # Set boundary conditions (in log-space) by mirroring log-pressure
     # TODO: re-generalize, is setting the BCs like this correct?
-    p[0] = p[1]
+    p[1] = p[2]
     p[end] = p[end - 1]
 
-    p_half[1] = p_half[2]
+    p_half[2] = p_half[3]
     p_half[end - 1] = p_half[end - 2]
-    p_half[0] = p_half[3]
+    p_half[1] = p_half[4]
     p_half[end] = p_half[end - 3]
 
     p .= exp.(p)

src/TurbulenceConvection.jl

@@ -22,7 +22,7 @@ const CPSGS = CLIMAParameters.Atmos.SubgridScale
 # For dispatching to inherited class
 struct BaseCase end
 
-up_sum(vals::OffsetArray) = off_arr(reshape(sum(vals; dims = 1), size(vals, 2)))
+up_sum(vals::AbstractArray) = off_arr(reshape(sum(vals; dims = 1), size(vals, 2)))
 
 function parse_namelist(namelist, keys...; default = nothing, valid_options = nothing)
     param = namelist

src/Turbulence_PrognosticTKE.jl

@@ -824,7 +824,7 @@ function get_GMV_CoVar(
             # -1, indexing phi_e.values[-1] yields the
             # _last_ value in the array. This is certainly
             # not intended
-            if k ≠ 0
+            if k ≠ 1
                 phi_diff = interp2pt(phi_e.values[k - 1] - gmv_phi[k - 1], phi_e.values[k] - gmv_phi[k])
                 psi_diff = interp2pt(psi_e.values[k - 1] - gmv_psi[k - 1], psi_e.values[k] - gmv_psi[k])
             else # just use 0th order approximation
@@ -838,7 +838,7 @@ function get_GMV_CoVar(
                 # -1, indexing phi_e.values[-1] yields the
                 # _last_ value in the array. This is certainly
                 # not intended
-                if k ≠ 0
+                if k ≠ 1
                     phi_diff = interp2pt(phi_u.values[i, k - 1] - gmv_phi[k - 1], phi_u.values[i, k] - gmv_phi[k])
                     psi_diff = interp2pt(psi_u.values[i, k - 1] - gmv_psi[k - 1], psi_u.values[i, k] - gmv_psi[k])
                 else # just use 0th order approximation
@@ -968,9 +968,9 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl
     @inbounds for k in real_center_indicies(grid)
         @inbounds for i in xrange(self.n_updrafts)
             input_p.updraft_top = self.UpdVar.updraft_top[i]
-            alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])) - 1, 0)
+            alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])), 1)
             avals = off_arr(self.UpdVar.Area.values[i, cinterior])
-            input_p.a_med = Statistics.median(avals[0:alen])
+            input_p.a_med = Statistics.median(avals[1:alen])
             input.zi = self.UpdVar.cloud_base[i]
             # entrainment
             input.buoy_ed_flux = self.EnvVar.TKE.buoy[k]
@@ -2092,7 +2092,7 @@ function compute_covariance_interdomain_src(
                 # -1, indexing phi_e.values[-1] yields the
                 # _last_ value in the array. This is certainly
                 # not intended
-                if k ≠ 0
+                if k ≠ 1
                     phi_diff =
                         interp2pt(phi_u.values[i, k - 1], phi_u.values[i, k]) -
                         interp2pt(phi_e.values[k - 1], phi_e.values[k])

src/Variables.jl

@@ -6,11 +6,12 @@ function zero_tendencies(self::VariablePrognostic, Gr::Grid)
 end
 
 function set_bcs(self::VariablePrognostic, Gr::Grid)
-    start_low = Gr.gw - 1
-    start_high = Gr.nzg - Gr.gw - 1
+    start_low = Gr.gw
+    start_high = Gr.nzg - Gr.gw
 
     if self.bc == "sym"
-        @inbounds for k in xrange(Gr.gw)
+        @inbounds for kk in xrange(Gr.gw)
+            k = kk - 1
             self.values[start_high + k + 1] = self.values[start_high - k]
             self.values[start_low - k] = self.values[start_low + 1 + k]
 
@@ -30,7 +31,8 @@ function set_bcs(self::VariablePrognostic, Gr::Grid)
         self.new[start_high] = 0.0
         self.new[start_low] = 0.0
 
-        @inbounds for k in xrange(1, Gr.gw)
+        @inbounds for kk in xrange(1, Gr.gw)
+            k = kk - 1
             self.values[start_high + k] = -self.values[start_high - k]
             self.values[start_low - k] = -self.values[start_low + k]
 
@@ -46,19 +48,21 @@ function set_bcs(self::VariablePrognostic, Gr::Grid)
 end
 
 function set_bcs(self::VariableDiagnostic, Gr::Grid)
-    start_low = Gr.gw - 1
-    start_high = Gr.nzg - Gr.gw
+    start_low = Gr.gw
+    start_high = Gr.nzg - Gr.gw + 1
 
     if self.bc == "sym"
-        @inbounds for k in xrange(Gr.gw)
+        @inbounds for kk in xrange(Gr.gw)
+            k = kk - 1
             self.values[start_high + k] = self.values[start_high - 1]
             self.values[start_low - k] = self.values[start_low + 1]
         end
 
     else
         self.values[start_high] = 0.0
         self.values[start_low] = 0.0
-        @inbounds for k in xrange(1, Gr.gw)
+        @inbounds for kk in xrange(1, Gr.gw)
+            k = kk - 1
             self.values[start_high + k] = 0.0  #-self.values[start_high - k ]
             self.values[start_low - k] = 0.0 #-self.values[start_low + k ]
         end

src/python_primitives.jl

@@ -7,21 +7,17 @@ fabs(a) = abs(a)
 
 argwhere(a) = findall(a .> 0)
 # Use stop-1 to emulate python ranges
-xrange(start, stop, step = 1) = range(start, stop - 1; step = step)
+xrange(start, stop, step = 1) = range(start + 1, stop; step = step)
 xrange(stop) = xrange(0, stop)
 
-pyzeros(n::Int) = OffsetArray(zeros(n), 0:(n - 1))
-pyzeros(m::Int, n::Int) = OffsetArray(zeros(m, n), 0:(m - 1), 0:(n - 1))
 
-pyones(n::Int) = OffsetArray(ones(n), 0:(n - 1))
-pyones(m::Int, n::Int) = OffsetArray(ones(m, n), 0:(m - 1), 0:(n - 1))
+pyzeros(n::Int) = zeros(n)
+pyzeros(m::Int, n::Int) = zeros(m, n)
 
-function off_arr(a::AbstractArray)
-    dims = ntuple(ndims(a)) do i
-        0:(size(a, i) - 1)
-    end
-    return OffsetArray(a, dims)
-end
+pyones(n::Int) = ones(n)
+pyones(m::Int, n::Int) = ones(m, n)
+
+off_arr(a::AbstractArray) = a
 
 pow(a::Real, b::Real) = a^b
 pow(a::AbstractVector, b::AbstractVector) = a .^ b

src/turbulence_functions.jl

@@ -262,8 +262,7 @@ function construct_tridiag_diffusion(grid, dt, rho_ae_K_m, rho, ae, a, b, c)
 end
 
 function tridiag_solve(b_rhs, a, b, c)
-    # Note that `1:end` is zero-based indexing.
-    A = Tridiagonal(a[1:end], parent(b), c[0:(end - 1)])
+    A = Tridiagonal(a[2:end], parent(b), c[1:(end - 1)])
     return A \ parent(b_rhs)
 end