Merge #1329

1329: Add SLEVE coordinate option r=akshaysridhar a=akshaysridhar

Adds exponential decay of vertical levels, with specified decay-scale, and relaxation strength parameters. Follows Clima Design docs eq. (3.5). 
- [x] Code
- [x] Unit-test
- [x] Docs
Linked issue: CliMA/ClimaAtmos.jl#1774

**Example:** 

- Plots show cell-center values of `z` coordinate in a test space with uniform vertical stretching.

**_Default Linear_**
<img width="300" alt="GCS_Linear" src="https://github.com/CliMA/ClimaCore.jl/assets/43710045/6927fb01-8c81-448b-a223-2b99c74fcf9c">

**_Optional SLEVE_**
<img width="300" alt="SLEVE_Schar_et_al" src="https://github.com/CliMA/ClimaCore.jl/assets/43710045/d53a72b0-7fcb-4294-aa0a-fb21042b4485">



Co-authored-by: Akshay Sridhar <[email protected]>

docs/refs.bib

@@ -146,6 +146,19 @@ @article{Ronchi1996
   doi = {10.1006/jcph.1996.0047}
 }
 
+@article{Schar2002,
+      author = {Christoph Schär and Daniel Leuenberger and Oliver Fuhrer and Daniel Lüthi and Claude Girard"},
+      title = {A New Terrain-Following Vertical Coordinate Formulation for Atmospheric Prediction Models"},
+      journal = {Monthly Weather Review},
+      year = {2002},
+      publisher = {American Meteorological Society},
+      address = {Boston MA, USA},
+      volume = {130},
+      number = {10},
+      doi = {https://doi.org/10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2},
+      pages= {2459 - 2480},
+}
+
 @article{Taylor2010,
   title={A compatible and conservative spectral element method on unstructured grids},
   author={Taylor, Mark A and Fournier, Aim\'{e}},

docs/src/api.md

@@ -250,6 +250,8 @@ Fields.Δz_field
 ## Hypsography
 
 ```@docs
+Hypsography.LinearAdaption
+Hypsography.SLEVEAdaption
 Hypsography.diffuse_surface_elevation!
 ```
 

src/Hypsography/Hypsography.jl

@@ -18,6 +18,23 @@ struct LinearAdaption{F <: Union{Fields.Field, Nothing}} <: HypsographyAdaption
     surface::F
 end
 
+"""
+    SLEVEAdaption(surface::Field, ηₕ::FT, s::FT)
+
+Locate vertical levels using an exponential function between the surface field and the top
+of the domain, using the method of [Schar2002](@cite). This method is modified
+such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper 
+bounds represent the domain bottom and top respectively. `s` governs the decay rate. 
+If the decay-scale is poorly specified (i.e., `s * zₜ` is lower than the maximum 
+surface elevation), a warning is thrown and `s` is adjusted such that it `szₜ > maximum(z_surface)`. 
+"""
+struct SLEVEAdaption{F <: Fields.Field, FT <: Real} <: HypsographyAdaption
+    # Union can be removed once deprecation removed.
+    surface::F
+    ηₕ::FT
+    s::FT
+end
+
 # deprecated in 0.10.12
 LinearAdaption() = LinearAdaption(nothing)
 
@@ -66,11 +83,31 @@ function ExtrudedFiniteDifferenceSpace(
 
     grad = Operators.Gradient()
     wdiv = Operators.WeakDivergence()
+    z_surface = Fields.field_values(adaption.surface)
+
+    FT = eltype(z_surface)
 
+    # TODO: Function dispatch
     if adaption isa LinearAdaption
-        z_surface = Fields.field_values(adaption.surface)
         fZ_data = @. z_ref + (1 - z_ref / z_top) * z_surface
         fZ = Fields.Field(fZ_data, face_space)
+    elseif adaption isa SLEVEAdaption
+        ηₕ = adaption.ηₕ
+        s = adaption.s
+        @assert FT(0) <= ηₕ <= FT(1)
+        @assert s >= FT(0)
+        η = @. z_ref ./ z_top
+        if s * z_top <= maximum(z_surface)
+            @warn "Decay scale (s*z_top = $(s*z_top)) must be higher than max surface elevation (max(z_surface) = $(maximum(z_surface))). 
+                   \n Returning (5/3)*max(zₛ)"
+            s = @. maximum(z_surface) / z_top * eltype(z_surface) * (5 // 3)
+        end
+        fZ_data = @. ifelse(
+            η <= ηₕ,
+            η * z_top + z_surface * (sinh((ηₕ - η) / s / ηₕ)) / (sinh(1 / s)),
+            η * z_top,
+        )
+        fZ = Fields.Field(fZ_data, face_space)
     end
 
     # Take the horizontal gradient for the Z surface field