Convert to 1-based indexing

integration_tests/utils/Cases.jl

@@ -1010,8 +1010,8 @@ 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]
@@ -1047,15 +1047,15 @@ 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
@@ -1064,7 +1064,8 @@ function update_forcing(self::CasesBase{TRMM_LBA}, GMV::GridMeanVariables, TS::T
                         (TS.t - self.rad_time[ind1]) + 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

integration_tests/utils/main.jl

@@ -30,13 +30,21 @@ end
 
 function TurbulenceConvection.initialize(self::Simulation1d, namelist)
     Cases.initialize_reference(self.Case, self.Gr, self.Ref, self.Stats)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0")
     Cases.initialize_profiles(self.Case, self.Gr, self.GMV, self.Ref)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.1")
     Cases.initialize_surface(self.Case, self.Gr, self.Ref)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.2")
     Cases.initialize_forcing(self.Case, self.Gr, self.Ref, self.GMV)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.3")
     Cases.initialize_radiation(self.Case, self.Gr, self.Ref, self.GMV)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.4")
     TurbulenceConvection.initialize(self.Turb, self.Case, self.GMV, self.Ref)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.5")
     TurbulenceConvection.initialize_io(self)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.6")
     TurbulenceConvection.io(self)
+    TurbulenceConvection.debug_nans(self.GMV, self.Turb, "0.7")
 
     return
 end
@@ -45,15 +53,24 @@ function run(self::Simulation1d)
     iter = 0
     TurbulenceConvection.open_files(self.Stats) # #removeVarsHack
     while self.TS.t <= self.TS.t_max
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "1, iter=$iter")
         TurbulenceConvection.zero_tendencies(self.GMV)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "2, iter=$iter")
         Cases.update_surface(self.Case, self.GMV, self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "3, iter=$iter")
         Cases.update_forcing(self.Case, self.GMV, self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "4, iter=$iter")
         Cases.update_radiation(self.Case, self.GMV, self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "5, iter=$iter")
         TurbulenceConvection.update(self.Turb, self.GMV, self.Case, self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "6, iter=$iter")
         TurbulenceConvection.update(self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "7, iter=$iter")
         # Apply the tendencies, also update the BCs and diagnostic thermodynamics
         TurbulenceConvection.update(self.GMV, self.TS)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "8, iter=$iter")
         TurbulenceConvection.update_GMV_diagnostics(self.Turb, self.GMV)
+        TurbulenceConvection.debug_nans(self.GMV, self.Turb, "9, iter=$iter")
 
         if mod(iter, 100) == 0
             progress = self.TS.t / self.TS.t_max

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

@@ -10,6 +10,7 @@ struct Grid{A1, FT}
     nzg::Int
     z_half::A1
     z::A1
+    zLL::FT
     function Grid(namelist)
 
         #Get the grid spacing
@@ -24,31 +25,32 @@ 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)
+        zLL = z_half[gw + 1]
+        return new{A1, FT}(zmin, zmax, cinterior, finterior, dz, dzi, gw, nz, nzg, z_half, z, zLL)
     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,41 @@ 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 debug_nans(GMV, Turb, loc)
+    # B = (
+    #     any(isnan.(Turb.mixing_length)),
+    #     any(isnan.(Turb.EnvVar.TKE.values)),
+    #     any(isnan.(GMV.TKE.values)),
+    #     any(isnan.(diffusivity_h(Turb).values)),
+    #     any(isnan.(diffusivity_m(Turb).values)),
+    #     any(isnan.(GMV.H.values)),
+    #     any(isnan.(GMV.QT.values)),
+    #     any(isnan.(GMV.QL.values)),
+    #     any(isnan.(GMV.T.values)),
+    #     any(isnan.(GMV.B.values)),
+    #     any(isnan.(GMV.U.values)),
+    #     any(isnan.(GMV.V.values)),
+    #     any(isnan.(GMV.H.new)),
+    #     any(isnan.(GMV.QT.new)),
+    #     any(isnan.(GMV.U.new)),
+    #     any(isnan.(GMV.V.new)),
+    #     any(isnan.(GMV.H.mf_update)),
+    #     any(isnan.(GMV.QT.mf_update)),
+    #     any(isnan.(GMV.U.mf_update)),
+    #     any(isnan.(GMV.V.mf_update)),
+    #     any(isnan.(GMV.H.tendencies)),
+    #     any(isnan.(GMV.QT.tendencies)),
+    #     any(isnan.(GMV.U.tendencies)),
+    #     any(isnan.(GMV.V.tendencies)),
+    # )
+    # if any(B)
+    #     @show loc
+    #     @show B
+    #     error("Got nans")
+    # end
+end
 
 function parse_namelist(namelist, keys...; default = nothing, valid_options = nothing)
     param = namelist