ZonalRidge orography added, Geometry with lond,latd

src/boundaries.jl

@@ -81,7 +81,47 @@ function initialize_orography!( orog::AbstractOrography,
     geopot_surf .*= gravity                 # turn orography into surface geopotential
     spectral_truncation!(geopot_surf,lmax,mmax) # set the lmax+1 harmonics to zero
 
-    # SCALE OROGRAPHY (or disable)
+    scale_orography!(orog,P)                # make whatever mountains bigger/smaller
+end
+
+function initialize_orography!( orog::AbstractOrography,
+                                ::Type{<:ZonalRidge},
+                                P::Parameters{M},
+                                S::SpectralTransform,
+                                G::Geometry) where {M<:Union{ShallowWater,PrimitiveEquation}}
+
+    @unpack gravity, rotation_earth, radius_earth = P
+    @unpack lmax, mmax = S
+
+    @unpack orography, geopot_surf = orog
+
+    ηₛ = 1
+    η₀ = 0.252
+    ηᵥ = (ηₛ-η₀)*π/2
+    u₀ = 35
+    s = u₀*cos(ηᵥ)^(3/2)
+    aΩ = radius_earth*rotation_earth
+    g⁻¹ = inv(gravity)
+    ϕ = G.latds
+
+    for ij in eachindex(ϕ,orography.data)
+        sinϕ = sind(ϕ[ij])
+        cosϕ = cosd(ϕ[ij])
+        orography[ij] = g⁻¹*s*(s*(-2*sinϕ^6*(cosϕ^2 + 1/3) + 10/63) + (8/5*cosϕ^3*(sinϕ^2 + 2/3) - π/4)*aΩ)
+    end
+
+    spectral!(geopot_surf,orography,S)      # to grid-point space
+    geopot_surf .*= gravity                 # turn orography into surface geopotential
+    spectral_truncation!(geopot_surf,lmax,mmax) # set the lmax+1 harmonics to zero
+
+    scale_orography!(orog,P)                # make whatever mountains bigger/smaller
+end
+
+function scale_orography!(  orog::AbstractOrography,
+                            P::Parameters)
+
+    @unpack orography, geopot_surf = orog
     orography .*= P.orography_scale
     geopot_surf .*= P.orography_scale
+    return nothing
 end

src/column_variables.jl

@@ -9,8 +9,8 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV
 
     # COORDINATES
     nlev::Int = 0                        # number of vertical levels
-    lon::NF = NaN                        # longitude
-    lat::NF = NaN                        # latitude, needed for shortwave radiation
+    lond::NF = NaN                       # longitude
+    latd::NF = NaN                       # latitude, needed for shortwave radiation
     nband::Int = 0                       # number of radiation bands, needed for radiation
     n_stratosphere_levels::Int = 0       # number of stratospheric levels, needed for radiation
 
@@ -120,8 +120,8 @@ function get_column!(   C::ColumnVariables,
 
     @boundscheck C.nlev == D.nlev || throw(BoundsError)
 
-    C.lat = G.lats[ij]        # pull latitude, longitude [˚N,˚E] for gridpoint ij from Geometry
-    C.lon = G.lons[ij]
+    C.latd = G.latds[ij]        # pull latitude, longitude [˚N,˚E] for gridpoint ij from Geometry
+    C.lond = G.londs[ij]
 
     # surface pressure (logarithm used in dynamics, convert back here)
     C.log_pres = D.surface.pres_grid[ij]

src/geometry.jl

@@ -30,8 +30,8 @@ struct Geometry{NF<:AbstractFloat}      # NF: Number format
     lond::Vector{Float64}   # array of longitudes in degrees (0...360˚)
 
     # COORDINATES
-    lons::Vector{NF}        # longitude (0...2π) for each grid point in ring order
-    lats::Vector{NF}        # latitude (π/2...-π/2) for each grid point in ring order
+    londs::Vector{NF}       # longitude (-180˚...180˚) for each grid point in ring order
+    latds::Vector{NF}       # latitude (-90˚...˚90) for each grid point in ring order
 
     # SINES AND COSINES OF LATITUDE
     sinlat::Vector{NF}              # sin of latitudes
@@ -98,7 +98,11 @@ function Geometry(P::Parameters,Grid::Type{<:AbstractGrid})
     lond = lon*360/2π                               # array of longitudes in degrees 0...360˚
 
     # COORDINATES for every grid point in ring order
-    lats,lons = get_colatlons(Grid,nlat_half)       # in radians
+    latds,londs = get_colatlons(Grid,nlat_half)     # in radians
+    latds .*= (360/2π)                              # in degrees
+    latds .= 90 .- latds                            # in -90˚...90˚N
+    londs .*= (360/2π)                              # in degrees
+    londs .-= 180                                   # -180˚...180˚
 
     # SINES AND COSINES OF LATITUDE
     sinlat = sin.(lat)
@@ -155,7 +159,7 @@ function Geometry(P::Parameters,Grid::Type{<:AbstractGrid})
     # conversion to number format NF happens here
     Geometry{P.NF}( Grid,nresolution,
                     nlon_max,nlon,nlat,nlev,nlat_half,npoints,radius_earth,
-                    lat,latd,colat,colatd,lon,lond,lons,lats,
+                    lat,latd,colat,colatd,lon,lond,londs,latds,
                     sinlat,coslat,coslat⁻¹,coslat²,coslat⁻²,f_coriolis,
                     n_stratosphere_levels,
                     σ_levels_half,σ_levels_full,σ_levels_thick,σ_levels_thick⁻¹_half,σ_f,

src/shortwave_radiation.jl

@@ -39,11 +39,11 @@ Compute average daily flux of solar radiation for an atmospheric column,
 from Hartmann (1994).
 """
 function solar!(column::ColumnVariables{NF}) where {NF<:AbstractFloat}
-    @unpack tyear, csol, lat = column
+    @unpack tyear, csol, latd = column
 
     # Compute cosine and sine of latitude
-    clat = cos(lat * π / 180)
-    slat = sin(lat * π / 180)
+    clat = cos(latd)
+    slat = sin(latd)
 
     # 1.0 Compute declination angle and Earth - Sun distance factor
     pigr = 2 * asin(1)
@@ -86,13 +86,13 @@ Compute solar radiation parametres for an atmospheric column.
 function sol_oz!(
     column::ColumnVariables{NF}, model::PrimitiveEquation
 ) where {NF<:AbstractFloat}
-    @unpack tyear, lat = column
+    @unpack tyear, latd = column
     @unpack tropic_cancer = model.parameters
     @unpack solc, epssw = model.parameters.radiation_coefs
 
     # Compute cosine and sine of latitude
-    clat = cos(lat * π / 180)
-    slat = sin(lat * π / 180)
+    clat = cos(latd)
+    slat = sin(latd)
 
     # α = year phase ( 0 - 2pi, 0 = winter solstice = 22dec.h00 )
     α = 4 * asin(1) * (tyear + 10 / 365)

test/orography.jl

@@ -0,0 +1,5 @@
+@testset "Orographies" begin
+    initialize_speedy(PrimitiveEquation,orography=NoOrography)
+    initialize_speedy(PrimitiveEquation,orography=EarthOrography)
+    initialize_speedy(PrimitiveEquation,orography=ZonalRidge)
+end