Merge pull request #3584 from CliMA/dy/face_density

Change density on cell faces to match dycore paper

Author: dennisYatunin

NEWS.md

@@ -12,6 +12,13 @@ If `netcdf_interpolation_num_points` is not provided, `ClimaAtmos` will
 determine it automatically by matching approximately the same number of points
 as the model grid.
 
+### Change reconstruction of density on cell faces for stretched grids
+
+PR [3584](https://github.com/CliMA/ClimaAtmos.jl/pull/3584) changes the weighted
+interpolation of density from centers to faces so that it uses `ᶜJ` and `ᶠJ`,
+rather than `ᶜJ` and `ᶠint(ᶜJ)`. As of ClimaCore v0.14.25, `ᶠJ` is no longer
+equivalent to `ᶠint(ᶜJ)` for stretched grids.
+
 v0.28.5
 -------
 

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-216
+217
 
 # **README**
 #
@@ -21,6 +21,9 @@
 
 #=
 
+217
+- Change reconstruction of density on cell faces
+
 216
 - Change prescribed aerosols in `aquaplanet_nonequil_allsky_gw_res.yml`
 

src/parameterized_tendencies/microphysics/precipitation.jl

@@ -325,30 +325,31 @@ function compute_precipitation_surface_fluxes!(
     (; surface_rain_flux, surface_snow_flux) = p.precipitation
     (; col_integrated_precip_energy_tendency,) = p.conservation_check
     (; ᶜwᵣ, ᶜwₛ, ᶜspecific) = p.precomputed
-
-    (; ᶠtemp_scalar) = p.scratch
-    slg = Fields.level(Fields.local_geometry_field(ᶠtemp_scalar), Fields.half)
-
-    # Constant extrapolation: - put values from bottom cell center to bottom cell face
-    ˢρ = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), axes(slg))
-    # For density this is equivalent with ᶠwinterp(ᶜJ, Y.c.ρ) and therefore
-    # consistent with the way we do vertical advection
-    ˢqᵣ = Fields.Field(
+    ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
+    sfc_J = Fields.level(ᶠJ, Fields.half)
+    sfc_space = axes(sfc_J)
+
+    # Jacobian-weighted extrapolation from interior to surface, consistent with
+    # the reconstruction of density on cell faces, ᶠρ = ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ
+    int_J = Fields.Field(Fields.field_values(Fields.level(ᶜJ, 1)), sfc_space)
+    int_ρ = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), sfc_space)
+    sfc_ρ = @. lazy(int_ρ * int_J / sfc_J)
+
+    # Constant extrapolation to surface, consistent with simple downwinding
+    sfc_qᵣ = Fields.Field(
         Fields.field_values(Fields.level(ᶜspecific.q_rai, 1)),
-        axes(slg),
+        sfc_space,
     )
-    ˢqₛ = Fields.Field(
+    sfc_qₛ = Fields.Field(
         Fields.field_values(Fields.level(ᶜspecific.q_sno, 1)),
-        axes(slg),
+        sfc_space,
     )
-    ˢwᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵣ, 1)), axes(slg))
-    ˢwₛ = Fields.Field(Fields.field_values(Fields.level(ᶜwₛ, 1)), axes(slg))
-
-    # Project the flux to CT3 vector and convert to physical units.
-    @. surface_rain_flux =
-        -projected_vector_data(CT3, ˢρ * ˢqᵣ * Geometry.WVector(ˢwᵣ), slg)
-    @. surface_snow_flux =
-        -projected_vector_data(CT3, ˢρ * ˢqₛ * Geometry.WVector(ˢwₛ), slg)
+    sfc_wᵣ = Fields.Field(Fields.field_values(Fields.level(ᶜwᵣ, 1)), sfc_space)
+    sfc_wₛ = Fields.Field(Fields.field_values(Fields.level(ᶜwₛ, 1)), sfc_space)
+
+    @. surface_rain_flux = sfc_ρ * sfc_qᵣ * (-sfc_wᵣ)
+    @. surface_snow_flux = sfc_ρ * sfc_qₛ * (-sfc_wₛ)
 end
 
 function precipitation_tendency!(

src/prognostic_equations/implicit/implicit_solver.jl

@@ -560,6 +560,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     ᶜuₕ = Y.c.uₕ
     ᶠu₃ = Y.f.u₃
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ
     ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ
 
@@ -581,7 +582,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix()
 
     @. ᶜadvection_matrix =
-        -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ))
+        -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ)
 
     if use_derivative(topography_flag)
         ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)]
@@ -663,7 +664,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
         #TODO: tetsing explicit vs implicit
         #@. ∂ᶜρe_tot_err_∂ᶜρe_tot +=
         #    dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅
-        #    DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) ⋅ ᶠright_bias_matrix() ⋅
+        #    DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅
         #    DiagonalMatrixRow(
         #        -(1 + ᶜkappa_m) / ᶜρ * ifelse(
         #            ᶜh_tot == 0,
@@ -677,7 +678,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
         #TODO: tetsing explicit vs implicit
         #@. ∂ᶜρe_tot_err_∂ᶜρq_tot =
         #    dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅
-        #    DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) ⋅ ᶠright_bias_matrix() ⋅
+        #    DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅
         #    DiagonalMatrixRow(
         #        -(ᶜkappa_m) * ∂e_int_∂q_tot / ᶜρ * ifelse(
         #            ᶜh_tot == 0,
@@ -691,7 +692,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
         #TODO: tetsing explicit vs implicit
         #@. ∂ᶜρq_tot_err_∂ᶜρq_tot =
         #    dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅
-        #    DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) ⋅ ᶠright_bias_matrix() ⋅
+        #    DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅
         #    DiagonalMatrixRow(
         #        -1 / ᶜρ * ifelse(
         #            ᶜspecific.q_tot == 0,
@@ -706,7 +707,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
             ᶜwₚ = MatrixFields.get_field(p, wₚ_name)
             @. ∂ᶜρqₚ_err_∂ᶜρqₚ =
                 dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρ)) ⋅ ᶠright_bias_matrix() ⋅
+                DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅
+                ᶠright_bias_matrix() ⋅
                 DiagonalMatrixRow(-Geometry.WVector(ᶜwₚ) / ᶜρ) - (I,)
         end
 
@@ -881,13 +883,13 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
                             ᶠu³ʲs.:(1),
                             draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
                         ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
+                    ) / ᶠJ,
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜJ * ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp) *
                     ∂e_int_∂q_tot,
                 )
             @. ᶠbidiagonal_matrix_ct3_2 =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
+                DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
                 ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
                 DiagonalMatrixRow(
                     Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp) *
@@ -906,12 +908,12 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
                             ᶠu³ʲs.:(1),
                             draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)),
                         ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
+                    ) / ᶠJ,
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜJ * ᶜkappa_mʲ * (ᶜρʲs.:(1))^2 / ((ᶜkappa_mʲ + 1) * ᶜp),
                 )
             @. ᶠbidiagonal_matrix_ct3_2 =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1))) ⋅
+                DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
                 ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅
                 DiagonalMatrixRow(
                     Y.c.sgsʲs.:(1).ρa * ᶜkappa_mʲ / ((ᶜkappa_mʲ + 1) * ᶜp),
@@ -923,7 +925,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
                 matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)]
             @. ᶜadvection_matrix =
                 -(ᶜadvdivᵥ_matrix()) ⋅
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ, ᶜρʲs.:(1)))
+                DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ)
             @. ∂ᶜρaʲ_err_∂ᶜρaʲ =
                 dtγ * ᶜadvection_matrix ⋅
                 ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1))) ⋅

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -54,48 +54,54 @@ end
 
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none})
     ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ))))
+    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠu³ * ᶠinterp(ᶜχ))))
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order})
     ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ))))
+    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ))))
 end
 @static if pkgversion(ClimaCore) ≥ v"0.14.22"
     function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:vanleer_limiter})
         ᶜJ = Fields.local_geometry_field(ᶜρ).J
+        ᶠJ = Fields.local_geometry_field(ᶠu³).J
         return @. lazy(
-            -(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠlin_vanleer(ᶠu³, ᶜχ, dt))),
+            -(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠlin_vanleer(ᶠu³, ᶜχ, dt))),
         )
     end
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order})
     ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ))))
+    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind3(ᶠu³, ᶜχ))))
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:boris_book})
     ᶜJ = Fields.local_geometry_field(ᶜρ).J
+    ᶠJ = Fields.local_geometry_field(ᶠu³).J
     return @. lazy(
         -(ᶜadvdivᵥ(
-            ᶠwinterp(ᶜJ, ᶜρ) * (
+            ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * (
                 ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_boris_book(
                     ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
                     ᶜχ / dt -
-                    ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
+                    ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
                 )
             ),
         )),
     )
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak})
     ᶜJ = Fields.local_geometry_field(ᶜρ).J
+    ᶠJ = Fields.local_geometry_field(ᶠu³).J
     return @. lazy(
         -(ᶜadvdivᵥ(
-            ᶠwinterp(ᶜJ, ᶜρ) * (
+            ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * (
                 ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak(
                     ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
                     ᶜχ / dt,
                     ᶜχ / dt -
-                    ᶜadvdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
+                    ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
                 )
             ),
         )),
@@ -114,10 +120,11 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; dt) = p
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     (; ᶠgradᵥ_ᶜΦ) = p.core
     (; ᶜh_tot, ᶜspecific, ᶠu³, ᶜp) = p.precomputed
 
-    @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠu³)
+    @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠu³)
 
     # Central advection of active tracers (e_tot and q_tot)
     vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜh_tot, dt, Val(:none))
@@ -133,11 +140,11 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     if moisture_model isa NonEquilMoistModel
         (; ᶜwₗ, ᶜwᵢ) = p.precomputed
         @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ(
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ *
             ᶠright_bias(Geometry.WVector(-(ᶜwₗ)) * ᶜspecific.q_liq),
         )
         @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ(
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ *
             ᶠright_bias(Geometry.WVector(-(ᶜwᵢ)) * ᶜspecific.q_ice),
         )
     end
@@ -149,11 +156,11 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         # using downward biasing and free outflow bottom boundary condition
         (; ᶜwᵣ, ᶜwₛ) = p.precomputed
         @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ *
             ᶠright_bias(Geometry.WVector(-(ᶜwᵣ)) * ᶜspecific.q_rai),
         )
         @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(
-            ᶠwinterp(ᶜJ, Y.c.ρ) *
+            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ *
             ᶠright_bias(Geometry.WVector(-(ᶜwₛ)) * ᶜspecific.q_sno),
         )
     end

src/prognostic_equations/water_advection.jl

@@ -3,14 +3,16 @@ import ClimaCore: Fields
 function vertical_advection_of_water_tendency!(Yₜ, Y, p, t)
 
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     (; ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed
 
     if !(p.atmos.moisture_model isa DryModel)
-        @. Yₜ.c.ρ -= ᶜprecipdivᵥ(ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠright_bias(-(ᶜwₜqₜ)))
+        @. Yₜ.c.ρ -=
+            ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₜqₜ)))
         @. Yₜ.c.ρe_tot -=
-            ᶜprecipdivᵥ(ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠright_bias(-(ᶜwₕhₜ)))
+            ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₕhₜ)))
         @. Yₜ.c.ρq_tot -=
-            ᶜprecipdivᵥ(ᶠwinterp(ᶜJ, Y.c.ρ) * ᶠright_bias(-(ᶜwₜqₜ)))
+            ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₜqₜ)))
     end
     return nothing
 end