More test cases on the sphere with topography #188

filename = "src/meshes/ETOPO_2022_v1_60s_N90W180_surface.nc"
#=
mesh_file = Trixi.download(
    "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf/ETOPO_2022_v1_60s_N90W180_surface.nc",
    joinpath(@__DIR__, "ETOPO_2022_v1_60s_N90W180_surface.nc")
)
=#
earth_topography = ncread(filename, "z")

function makepositive(x)
    if x >= 0
        return x
    else
        return zero(typeof(x))
    end
end

earth_topography = makepositive.(earth_topography)
# Add extra column and row for periodicity
earth_topography = [earth_topography; earth_topography[1:1, :]]
earth_topography = [earth_topography earth_topography[:, 1:1]]

function initial_topography_earth(lat, lon, earth_topography)
    RealT = eltype(x_cart)
    delta_topography = (60.0 / 3600) * (2 * pi) / 360 # 60 seconds in radians

    i = round(Int, (lon + pi) / delta_topography + 1)
    j = round(Int, (lat + 0.5 * pi) / delta_topography + 1)

    i = clamp(i, 1, size(earth_topography, 1))
    j = clamp(j, 1, size(earth_topography, 2))

    return RealT(earth_topography[i, j])
end

mesh = TrixiAtmo.P4estMeshCubedSphereTopography(trees_per_cube_face..., EARTH_RADIUS, 30000,
                                                polydeg = polydeg,
                                                initial_refinement_level = 0,
                                                initial_topography = (lat, lon) -> initial_topography_earth(lat,
                                                                                                            lon,
                                                                                                            earth_topography),
                                                adapt_vertical_grid = TrixiAtmo.Sleve(0.7,
                                                                                      0.8))