diff --git a/input_data/vertical_coordinate/L137_phalf_levels.txt b/input_data/vertical_coordinate/L137_phalf_levels.txt deleted file mode 100644 index a1d30a50b..000000000 --- a/input_data/vertical_coordinate/L137_phalf_levels.txt +++ /dev/null @@ -1,139 +0,0 @@ -phalf -0.0000 -0.0200 -0.0310 -0.0467 -0.0683 -0.0975 -0.1361 -0.1861 -0.2499 -0.3299 -0.4288 -0.5496 -0.6952 -0.8690 -1.0742 -1.3143 -1.5928 -1.9134 -2.2797 -2.6954 -3.1642 -3.6898 -4.2759 -4.9262 -5.6441 -6.4334 -7.2974 -8.2397 -9.2634 -10.3720 -11.5685 -12.8561 -14.2377 -15.7162 -17.2945 -18.9752 -20.7610 -22.6543 -24.6577 -26.7735 -29.0039 -31.3512 -33.8174 -36.4047 -39.1149 -41.9493 -44.9082 -47.9915 -51.1990 -54.5299 -57.9834 -61.5607 -65.2695 -69.1187 -73.1187 -77.2810 -81.6182 -86.1450 -90.8774 -95.8280 -101.0047 -106.4153 -112.0681 -117.9714 -124.1337 -130.5637 -137.2703 -144.2624 -151.5493 -159.1403 -167.0450 -175.2731 -183.8344 -192.7389 -201.9969 -211.6186 -221.6146 -231.9954 -242.7719 -253.9549 -265.5556 -277.5852 -290.0548 -302.9762 -316.3607 -330.2202 -344.5663 -359.4111 -374.7666 -390.6450 -407.0583 -424.0190 -441.5395 -459.6321 -478.3096 -497.5845 -517.4198 -537.7195 -558.3430 -579.1926 -600.1668 -621.1624 -642.0764 -662.8084 -683.2620 -703.3467 -722.9795 -742.0855 -760.5996 -778.4661 -795.6396 -812.0847 -827.7756 -842.6959 -856.8376 -870.2004 -882.7910 -894.6222 -905.7116 -916.0815 -925.7571 -934.7666 -943.1399 -950.9082 -958.1037 -964.7584 -970.9046 -976.5737 -981.7968 -986.6036 -991.0230 -995.0824 -998.8081 -1002.2250 -1005.3562 -1008.2239 -1010.8487 -1013.2500 \ No newline at end of file diff --git a/input_data/vertical_coordinate/L31_phalf_levels.txt b/input_data/vertical_coordinate/L31_phalf_levels.txt deleted file mode 100644 index 863bae771..000000000 --- a/input_data/vertical_coordinate/L31_phalf_levels.txt +++ /dev/null @@ -1,33 +0,0 @@ -phalf -0.00000 -20.0000 -40.0000 -60.0000 -80.0000 -100.1574 -121.1638 -143.6299 -167.9520 -194.3569 -222.9420 -253.7096 -286.5963 -321.4966 -358.2815 -396.8118 -436.9459 -478.5425 -521.4575 -565.5364 -610.6005 -656.4283 -702.7315 -749.1255 -795.0949 -839.9533 -882.7980 -922.4593 -957.4447 -985.8772 -1005.4292 -1013.2500 \ No newline at end of file diff --git a/input_data/vertical_coordinate/L40_phalf_levels.txt b/input_data/vertical_coordinate/L40_phalf_levels.txt deleted file mode 100644 index 09a4b1c0e..000000000 --- a/input_data/vertical_coordinate/L40_phalf_levels.txt +++ /dev/null @@ -1,42 +0,0 @@ -phalf -0.0000 -20.0000 -40.0000 -60.0000 -80.0000 -100.0886 -120.6766 -142.1783 -164.9202 -189.1466 -215.0251 -242.6522 -272.0592 -303.2174 -336.0439 -370.4072 -406.1328 -443.0085 -480.7907 -519.2093 -557.9734 -596.7774 -635.3061 -673.2403 -710.2627 -746.0635 -780.3455 -812.8302 -843.2634 -871.4203 -897.1117 -920.1893 -940.5512 -958.1476 -972.9868 -985.1399 -994.7472 -1002.0236 -1007.2639 -1010.8487 -1013.2500 \ No newline at end of file diff --git a/input_data/vertical_coordinate/L60_phalf_levels.txt b/input_data/vertical_coordinate/L60_phalf_levels.txt deleted file mode 100644 index 9e19aff11..000000000 --- a/input_data/vertical_coordinate/L60_phalf_levels.txt +++ /dev/null @@ -1,62 +0,0 @@ -phalf -0.0000 -0.2000 -0.3843 -0.6365 -0.9564 -1.3448 -1.8058 -2.3478 -2.9850 -3.7397 -4.6462 -5.7565 -7.1322 -8.8366 -10.9483 -13.5647 -16.8064 -20.8227 -25.7989 -31.9642 -39.6029 -49.0671 -60.1802 -73.0663 -87.7274 -104.2288 -122.6137 -142.9017 -165.0886 -189.1466 -215.0251 -242.6523 -272.0593 -303.2174 -336.0439 -370.4072 -406.1328 -443.0086 -480.7907 -519.2093 -557.9734 -596.7774 -635.3060 -673.2403 -710.2627 -746.0635 -780.3455 -812.8303 -843.2634 -871.4203 -897.1118 -920.1893 -940.5511 -958.1477 -972.9868 -985.1399 -994.7472 -1002.0236 -1007.2639 -1010.8487 -1013.2500 \ No newline at end of file diff --git a/input_data/vertical_coordinate/L91_phalf_levels.txt b/input_data/vertical_coordinate/L91_phalf_levels.txt deleted file mode 100644 index 8906f6210..000000000 --- a/input_data/vertical_coordinate/L91_phalf_levels.txt +++ /dev/null @@ -1,93 +0,0 @@ -phalf -0.0000 -0.0200 -0.0398 -0.0739 -0.1291 -0.2141 -0.3395 -0.5175 -0.7617 -1.0872 -1.5099 -2.0464 -2.7136 -3.5282 -4.5069 -5.6652 -7.0181 -8.5795 -10.3617 -12.3759 -14.6316 -17.1371 -19.8987 -22.9216 -26.2090 -29.7630 -33.5843 -37.6720 -42.0242 -46.6378 -51.5086 -56.6316 -61.9984 -67.5973 -73.4150 -79.4434 -85.7016 -92.2162 -99.0182 -106.1445 -113.6382 -121.5502 -129.9403 -138.8558 -148.3260 -158.3816 -169.0545 -180.3786 -192.3889 -205.1222 -218.6172 -232.9140 -248.0547 -264.0833 -281.0456 -298.9895 -317.9651 -338.0245 -359.2221 -381.6144 -405.2606 -430.2069 -456.4813 -483.8505 -512.0662 -540.8577 -569.9401 -599.0310 -627.9668 -656.6129 -684.8491 -712.5573 -739.5739 -765.7697 -791.0376 -815.2774 -838.3507 -860.1516 -880.6080 -899.6602 -917.2205 -933.2247 -947.6584 -960.5245 -971.8169 -981.5301 -989.7322 -996.8732 -1002.8013 -1007.4431 -1010.8487 -1013.2500 \ No newline at end of file diff --git a/input_data/vertical_coordinate/vertical_levels.png b/input_data/vertical_coordinate/vertical_levels.png deleted file mode 100644 index cd5f63e03..000000000 Binary files a/input_data/vertical_coordinate/vertical_levels.png and /dev/null differ diff --git a/input_data/vertical_coordinate/vertical_resolution.ipynb b/input_data/vertical_coordinate/vertical_resolution.ipynb deleted file mode 100644 index cc54c4c4c..000000000 --- a/input_data/vertical_coordinate/vertical_resolution.ipynb +++ /dev/null @@ -1,179 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Modelling the vertical resolution of IFS for n levels" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "using CSV\n", - "using DataFrames\n", - "using PyPlot\n", - "using LsqFit" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "L31 = CSV.read(\"L31_phalf_levels.txt\",DataFrame)\n", - "L40 = CSV.read(\"L40_phalf_levels.txt\",DataFrame)\n", - "L60 = CSV.read(\"L60_phalf_levels.txt\",DataFrame)\n", - "L91 = CSV.read(\"L91_phalf_levels.txt\",DataFrame)\n", - "L137 = CSV.read(\"L137_phalf_levels.txt\",DataFrame);" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "p_half31 = Vector(Array(L31)[:,1])\n", - "p_half40 = Vector(Array(L40)[:,1])\n", - "p_half60 = Vector(Array(L60)[:,1])\n", - "p_half91 = Vector(Array(L91)[:,1])\n", - "p_half137 = Vector(Array(L137)[:,1]);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Fit a generalised logistic function" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "$Y(x) = A + \\frac{K - A}{(C + Q \\exp(-B(x-M)))^{1/\\nu}}$" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "7-element Vector{Float16}:\n", - " -0.2827\n", - " 0.871\n", - " 0.4136\n", - " 6.695\n", - " 10.336\n", - " 0.602\n", - " 5.812" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# fit based on L?\n", - "p0 = 1013.25 # surface pressure\n", - "y = p_half31/p0 # dependent data\n", - "x = Vector(range(0,1,length(y))) # independent data\n", - "\n", - "# generalised logistic\n", - "p0 = Float64[0,1,1,1,1,0,1]\n", - "\n", - "function model(x,p)\n", - " A,K,C,Q,B,M,ν = p\n", - " return @. A + (K-A)/(C+Q*exp(-B*(x-M)))^(1/ν)\n", - "end\n", - "\n", - "fit = curve_fit(model,x,y,p0)\n", - "pfit = coef(fit)\n", - "x = range(0,1,9)\n", - "yfit = model(x,Float16.(pfit))\n", - "Float16.(pfit) # round for shorter inclusion in the default parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "for (p,t) in zip([p_half31,p_half40,p_half60,p_half91,p_half137],\n", - " [\"L31\",\"L40\",\"L60\",\"L91\",\"L137\"])\n", - " plot(range(0,1,length(p)),p/p[end],label=t,alpha=.7)\n", - "end\n", - "\n", - "yfit[1] = 0\n", - "yfit[end] = 1\n", - "\n", - "plot(x,yfit,\"ko--\",label=\"GL fit\",alpha=.3)\n", - "\n", - "xlabel(\"level/nlev\")\n", - "ylabel(\"σ\")\n", - "grid(alpha=.3)\n", - "\n", - "legend()\n", - "tight_layout()\n", - "savefig(\"vertical_levels.png\")" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "9×2 Matrix{Float64}:\n", - " 0.0 0.0\n", - " 0.073686 0.05\n", - " 0.162284 0.14\n", - " 0.272612 0.26\n", - " 0.408969 0.42\n", - " 0.573089 0.6\n", - " 0.75429 0.77\n", - " 0.913631 0.9\n", - " 1.0 1.0" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# new default levels become\n", - "hcat(yfit,[0.0, 0.05, 0.14, 0.26, 0.42, 0.6, 0.77, 0.9, 1.0])" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Julia 1.7.0", - "language": "julia", - "name": "julia-1.7" - }, - "language_info": { - "file_extension": ".jl", - "mimetype": "application/julia", - "name": "julia", - "version": "1.7.0" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/src/abstract_types.jl b/src/abstract_types.jl index c4a2ad8a9..cd5ba289a 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -32,7 +32,6 @@ abstract type AbstractDrag{NF} end # PARAMETERIZATIONS abstract type AbstractParameterization{NF} end -abstract type BoundaryLayerDrag{NF} <: AbstractParameterization{NF} end abstract type TemperatureRelaxation{NF} <: AbstractParameterization{NF} end abstract type VerticalDiffusion{NF} <: AbstractParameterization{NF} end abstract type AbstractCondensation{NF} <: AbstractParameterization{NF} end diff --git a/src/dynamics/geopotential.jl b/src/dynamics/geopotential.jl index fb89c8a07..d4746f73b 100644 --- a/src/dynamics/geopotential.jl +++ b/src/dynamics/geopotential.jl @@ -31,24 +31,6 @@ function initialize_geopotential( σ_levels_full::Vector, return Δp_geopot_half, Δp_geopot_full end -# function lapserate_correction( σ_levels_full::Vector, -# σ_levels_half::Vector, -# Δp_geopot_full::Vector) - -# nlev = length(σ_levels_full) -# @assert nlev+1 == length(σ_levels_half) "σ half levels must have length nlev+1" -# @assert nlev == length(Δp_geopot_full) "σ half levels must have length nlev" - -# lapserate_corr = zeros(nlev) -# for k in 2:nlev-1 # only in the free troposphere -# # TODO reference -# lapserate_corr[k] = 0.5*Δp_geopot_full[k]* -# log(σ_levels_half[k+1]/σ_levels_full[k]) / log(σ_levels_full[k+1]/σ_levels_full[k-1]) -# end - -# return lapserate_corr -# end - """ $(TYPEDSIGNATURES) Compute spectral geopotential `geopot` from spectral temperature `temp` @@ -89,19 +71,6 @@ function geopotential!( geopot_k[lm] = geopot_k½ + temp_k[lm]*Δp_geopot_full[k] # 2nd onto full layer end end - - # # LAPSERATE CORRECTION IN THE FREE TROPOSPHERE (>nlev) - # # TODO only for spectral coefficients 1,: ? - # lmax,mmax = size(geopot) - # for k in 2:nlev-1 - # temp_k_above = progn.layers[k-1].timesteps[lf].temp - # temp_k_below = progn.layers[k+1].timesteps[lf].temp - # geopot = diagn.layers[k].dynamics_variables.geopot - - # for l in 1:lmax-1 - # geopot[l,l] += lapserate_corr[k]*(temp_k_below[l,l] - temp_k_above[l,l]) - # end - # end end """ @@ -109,9 +78,11 @@ $(TYPEDSIGNATURES) Calculate the geopotential based on `temp` in a single column. This exclues the surface geopotential that would need to be added to the returned vector. Function not used in the dynamical core but for post-processing and analysis.""" -function geopotential!( geopot::Vector, - temp::Vector, - C::DynamicsConstants) +function geopotential!( geopot::AbstractVector, + temp::AbstractVector, + C::DynamicsConstants, + geopot_surf::Real = 0) + nlev = length(geopot) (;Δp_geopot_half, Δp_geopot_full) = C # = R*Δlnp either on half or full levels @@ -120,7 +91,7 @@ function geopotential!( geopot::Vector, @boundscheck length(Δp_geopot_half) >= nlev || throw(BoundsError) # bottom layer - geopot[nlev] = temp[nlev]*Δp_geopot_full[end] + geopot[nlev] = geopot_surf + temp[nlev]*Δp_geopot_full[end] # OTHER FULL LAYERS, integrate two half-layers from bottom to top @inbounds for k in nlev-1:-1:1 diff --git a/src/physics/boundary_layer.jl b/src/physics/boundary_layer.jl index 5839396dd..39cfdd057 100644 --- a/src/physics/boundary_layer.jl +++ b/src/physics/boundary_layer.jl @@ -1,3 +1,5 @@ +abstract type BoundaryLayerDrag{NF} <: AbstractParameterization{NF} end + """Concrete type that disables the boundary layer drag scheme.""" struct NoBoundaryLayerDrag{NF} <: BoundaryLayerDrag{NF} end @@ -9,15 +11,16 @@ function initialize!( scheme::NoBoundaryLayerDrag, return nothing end -# function barrier +# function barrier for all boundary layer drags function boundary_layer_drag!( column::ColumnVariables, model::PrimitiveEquation) - boundary_layer_drag!(column,model.boundary_layer_drag) + boundary_layer_drag!(column,model.boundary_layer_drag,model) end """NoBoundaryLayer scheme just passes.""" function boundary_layer_drag!( column::ColumnVariables, - scheme::NoBoundaryLayerDrag) + scheme::NoBoundaryLayerDrag, + model::PrimitiveEquation) return nothing end @@ -57,6 +60,13 @@ function initialize!( scheme::LinearDrag, end end +# function barrier +function boundary_layer_drag!( column::ColumnVariables, + scheme::LinearDrag, + model::PrimitiveEquation) + boundary_layer_drag!(column,scheme) +end + """ $(TYPEDSIGNATURES) Compute tendency for boundary layer drag of a `column` and add to its tendencies fields""" @@ -73,4 +83,81 @@ function boundary_layer_drag!( column::ColumnVariables, v_tend[k] += kᵥ*v[k] end end -end \ No newline at end of file +end + +Base.@kwdef struct ConstantDrag{NF} <: BoundaryLayerDrag{NF} + drag::NF = 1e-3 +end + +ConstantDrag(SG::SpectralGrid;kwargs...) = ConstantDrag{SG.NF}(;kwargs...) + +initialize!(::ConstantDrag,::PrimitiveEquation) = nothing + +function boundary_layer_drag!( column::ColumnVariables, + scheme::ConstantDrag, + model::PrimitiveEquation) + column.boundary_layer_drag = scheme.drag +end + +"""Boundary layer drag coefficient from the bulk Richardson number, +following Frierson, 2006, Journal of the Atmospheric Sciences. +$(TYPEDFIELDS)""" +Base.@kwdef struct BulkRichardsonDrag{NF} <: BoundaryLayerDrag{NF} + "von Kármán constant [1]" + κ::NF = 0.4 + + "roughness length [m]" + z₀::NF = 3.21e-5 + + "Critical Richardson number for stable mixing cutoff [1]" + Ri_c::NF = 1 + + "Maximum drag coefficient, κ²/log(zₐ/z₀) for zₐ from reference temperature" + drag_max::Base.RefValue{NF} = Ref(zero(NF)) +end + +BulkRichardsonDrag(SG::SpectralGrid,kwargs...) = BulkRichardsonDrag{SG.NF}(;kwargs...) + +function initialize!(scheme::BulkRichardsonDrag, model::PrimitiveEquation) + + # Typical height Z of lowermost layer from geopotential of reference surface temperature + # minus surface geopotential (orography * gravity) + (;temp_ref) = model.atmosphere + (;Δp_geopot_full, gravity) = model.constants + Z = temp_ref*Δp_geopot_full[end] / gravity + + # maximum drag Cmax from that height, stable conditions would decrease Cmax towards 0 + # Frierson 2006, eq (12) + (;κ, z₀) = scheme + scheme.drag_max[] = (κ/log(Z/z₀))^2 +end + +function boundary_layer_drag!( column::ColumnVariables, + scheme::BulkRichardsonDrag, + model::PrimitiveEquation) + boundary_layer_drag!(column,scheme) +end + +function boundary_layer_drag!( + column::ColumnVariables, + scheme::BulkRichardsonDrag, +) + + (;Ri_c) = scheme + drag_max = scheme.drag_max[] + + # bulk Richardson number at lowermost layer from Frierson, 2006, eq. (15) + Ri_a = column.bulk_richardson[column.nlev] + + # clamp to get the cases, eq (12-14) + # if Ri_a > Ri_c then C = 0 + # if Ri_c > Ri_a > 0 then = κ^2/log(zₐ/z₀)^2 * (1-Ri_a/Ri_c)^2 + # if Ri_c < 0 then κ^2/log(zₐ/z₀)^2 + # but Z ≈ zₐ given that the lowermost layer height is proportional to temperature + # which doesn't change much with instantaneous temperature variations but with + # vertical resolution, hence κ^2/log(Z/z₀)^2 is precalculated in initialize! + Ri_a = clamp(Ri_a, 0, Ri_c) + column.boundary_layer_drag = drag_max*(1-Ri_a/Ri_c)^2 +end + + diff --git a/src/physics/column_variables.jl b/src/physics/column_variables.jl index 182ab5cc6..bce6e3f96 100644 --- a/src/physics/column_variables.jl +++ b/src/physics/column_variables.jl @@ -1,3 +1,15 @@ +# function barrier +function get_column!( + C::ColumnVariables, + D::DiagnosticVariables, + P::PrognosticVariables, + ij::Integer, # grid point index + jring::Integer, # ring index 1 around North Pole to J around South Pole + model::PrimitiveEquation, +) + get_column!(C, D, P, ij, jring, model.geometry, model.constants, model.orography, model.land_sea_mask) +end + """ $(TYPEDSIGNATURES) Update `C::ColumnVariables` by copying the prognostic variables from `D::DiagnosticVariables` @@ -8,6 +20,8 @@ function get_column!( C::ColumnVariables, ij::Integer, # grid point index jring::Integer, # ring index 1 around North Pole to J around South Pole G::Geometry, + constants::DynamicsConstants, + O::AbstractOrography, L::AbstractLandSeaMask) (;σ_levels_full,ln_σ_levels_full) = G @@ -19,6 +33,8 @@ function get_column!( C::ColumnVariables, C.ij = ij # grid-point index C.jring = jring # ring index j of column, used to index latitude vectors C.land_fraction = L.land_sea_mask[ij] + C.orography = O.orography[ij] + C.surface_geopotential = C.orography*constants.gravity # pressure [Pa] or [log(Pa)] lnpₛ = D.surface.pres_grid[ij] # logarithm of surf pressure used in dynamics @@ -31,6 +47,7 @@ function get_column!( C::ColumnVariables, C.u[k] = layer.grid_variables.u_grid[ij] C.v[k] = layer.grid_variables.v_grid[ij] C.temp[k] = layer.grid_variables.temp_grid[ij] + C.temp_virt[k] = layer.grid_variables.temp_virt_grid[ij] # actually diagnostic C.humid[k] = layer.grid_variables.humid_grid[ij] end @@ -45,18 +62,18 @@ function get_column!( C::ColumnVariables, D::DiagnosticVariables, P::PrognosticVariables, ij::Int, # grid point index - G::Geometry, - L::LandSeaMask) + model::PrimitiveEquation) - rings = eachring(G.Grid,G.nlat_half) + SG = model.spectral_grid + rings = eachring(SG.Grid,SG.nlat_half) jring = whichring(ij,rings) - get_column!(C,D,P,ij,jring,G,L) + get_column!(C,D,P,ij,jring,model) end function get_column( S::AbstractSimulation, - ij::Integer) + ij::Integer, + verbose::Bool = true) (;prognostic_variables, diagnostic_variables) = S - (;geometry, land_sea_mask) = S.model column = deepcopy(S.diagnostic_variables.columns[1]) reset_column!(column) @@ -65,13 +82,12 @@ function get_column( S::AbstractSimulation, diagnostic_variables, prognostic_variables, ij, - geometry, - land_sea_mask) + model) # execute all parameterizations for this column to return a consistent state parameterization_tendencies!(column,S.model) - @info "Receiving column at $(column.latd)˚N, $(column.lond)˚E." + verbose && @info "Receiving column at $(column.latd)˚N, $(column.lond)˚E." return column end diff --git a/src/physics/define_column.jl b/src/physics/define_column.jl index ff10b7bcc..01c97e68a 100644 --- a/src/physics/define_column.jl +++ b/src/physics/define_column.jl @@ -17,12 +17,13 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV lond::NF = 0 # longitude latd::NF = 0 # latitude, needed for shortwave radiation land_fraction::NF = 0 # fraction of the column that is over land + orography::NF = 0 # orography height [m] # PROGNOSTIC VARIABLES - const u::Vector{NF} = zeros(NF,nlev) # zonal velocity [m/s] - const v::Vector{NF} = zeros(NF,nlev) # meridional velocity [m/s] - const temp::Vector{NF} = zeros(NF,nlev) # temperature [K] - const humid::Vector{NF} = zeros(NF,nlev) # specific humidity [kg/kg] + const u::Vector{NF} = zeros(NF,nlev) # zonal velocity [m/s] + const v::Vector{NF} = zeros(NF,nlev) # meridional velocity [m/s] + const temp::Vector{NF} = zeros(NF,nlev) # absolute temperature [K] + const humid::Vector{NF} = zeros(NF,nlev) # specific humidity [kg/kg] # (log) pressure per layer, surface is prognostic, last element here, but precompute other layers too const ln_pres::Vector{NF} = zeros(NF,nlev+1) # logarithm of pressure [log(Pa)] @@ -34,9 +35,6 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV const temp_tend::Vector{NF} = zeros(NF,nlev) # absolute temperature [K/s] const humid_tend::Vector{NF} = zeros(NF,nlev) # specific humidity [kg/kg/s] - # DIAGNOSTIC VARIABLES - const geopot::Vector{NF} = zeros(NF,nlev) # gepotential height [m] - # FLUXES, arrays to be used for various parameterizations, on half levels incl top and bottom const flux_u_upward::Vector{NF} = zeros(NF,nlev+1) const flux_u_downward::Vector{NF} = zeros(NF,nlev+1) @@ -50,6 +48,10 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV const flux_humid_upward::Vector{NF} = zeros(NF,nlev+1) const flux_humid_downward::Vector{NF} = zeros(NF,nlev+1) + # boundary layer + boundary_layer_depth::Int = 0 + boundary_layer_drag::NF = 0 + surface_geopotential::NF = 0 surface_u::NF = 0 surface_v::NF = 0 surface_temp::NF = 0 @@ -66,7 +68,10 @@ Base.@kwdef mutable struct ColumnVariables{NF<:AbstractFloat} <: AbstractColumnV const dry_static_energy::Vector{NF} = zeros(NF,nlev) # Dry static energy const moist_static_energy::Vector{NF} = zeros(NF,nlev) # Moist static energy const sat_moist_static_energy::Vector{NF} = zeros(NF,nlev) # Saturation moist static energy - + const temp_virt::Vector{NF} = zeros(NF,nlev) # virtual temperature [K] + const bulk_richardson::Vector{NF} = zeros(NF,nlev) # bulk richardson number [1] + const geopot::Vector{NF} = zeros(NF,nlev) # gepotential height [m] + # and interpolated to half levels const humid_half::Vector{NF} = zeros(NF,nlev) # Specific humidity interpolated to half-levels const sat_humid_half::Vector{NF} = zeros(NF,nlev) # Saturation specific humidity interpolated to half-levels diff --git a/src/physics/surface_fluxes.jl b/src/physics/surface_fluxes.jl index 5e8cd3cb5..12165fb02 100644 --- a/src/physics/surface_fluxes.jl +++ b/src/physics/surface_fluxes.jl @@ -51,14 +51,16 @@ Base.@kwdef struct SurfaceWind{NF<:AbstractFloat} <: AbstractSurfaceWind{NF} "Wind speed of sub-grid scale gusts [m/s]" V_gust::NF = 5 - # TODO make this orography dependent - "Drag coefficient over land (orography = 0) [1]" + "Use (possibly) flow-dependent column.boundary_layer_drag coefficient" + use_boundary_layer_drag::Bool = true + + "Otherwise, drag coefficient over land (orography = 0) [1]" drag_land::NF = 2.4e-3 - "Drag coefficient over sea [1]" + "Otherwise, Drag coefficient over sea [1]" drag_sea::NF = 1.8e-3 - "Flux limiter [kg/m/s²]" + "Flux limiter to cap the max of surface momentum fluxes [kg/m/s²]" max_flux::NF = 0.1 end @@ -79,6 +81,11 @@ function surface_wind_stress!( column::ColumnVariables{NF}, # SPEEDY documentation eq. 50 column.surface_wind_speed = sqrt(surface_u^2 + surface_v^2 + V_gust^2) + # drag coefficient either from SurfaceEvaporation or from a central drag coefficient + drag_sea, drag_land = surface_wind.use_boundary_layer_drag ? + (column.boundary_layer_drag, column.boundary_layer_drag) : + (drag_sea, drag_land) + # surface wind stress: quadratic drag, fractional land-sea mask ρ = column.surface_air_density V₀ = column.surface_wind_speed @@ -98,8 +105,17 @@ function surface_wind_stress!( column::ColumnVariables{NF}, end Base.@kwdef struct SurfaceSensibleHeat{NF<:AbstractFloat} <: AbstractSurfaceHeat{NF} + + "Use (possibly) flow-dependent column.boundary_layer_drag coefficient" + use_boundary_layer_drag::Bool = true + + "Otherwise, use the following drag coefficient for heat fluxes over land" heat_exchange_land::NF = 1.2e-3 # for neutral stability + + "Otherwise, use the following drag coefficient for heat fluxes over sea" heat_exchange_sea::NF = 0.9e-3 + + "Flux limiter for surface heat fluxes [W/m²]" max_flux::NF = 100 end @@ -118,9 +134,14 @@ function sensible_heat_flux!( column::ColumnVariables{NF}, T = column.surface_temp land_fraction = column.land_fraction + # drag coefficient + drag_sea, drag_land = heat_flux.use_boundary_layer_drag ? + (column.boundary_layer_drag, column.boundary_layer_drag) : + (heat_exchange_sea, heat_exchange_land) + # SPEEDY documentation Eq. 54 and 56 - flux_land = ρ*heat_exchange_land*V₀*cₚ*(T_skin_land - T) - flux_sea = ρ*heat_exchange_sea*V₀*cₚ*(T_skin_sea - T) + flux_land = ρ*drag_land*V₀*cₚ*(T_skin_land - T) + flux_sea = ρ*drag_sea*V₀*cₚ*(T_skin_sea - T) # flux limiter flux_land = clamp(flux_land,-max_flux,max_flux) @@ -147,8 +168,18 @@ function sensible_heat_flux!( column::ColumnVariables{NF}, return nothing end +""" +Surface evaporation following a bulk formula with wind from model.surface_wind +$(TYPEDFIELDS)""" Base.@kwdef struct SurfaceEvaporation{NF<:AbstractFloat} <: AbstractEvaporation{NF} - moisture_exchange_land::NF = 1.2e-3 # for neutral stability + + "Use column.boundary_layer_drag coefficient" + use_boundary_layer_drag::Bool = true + + "Otherwise, use the following drag coefficient for evaporation over land" + moisture_exchange_land::NF = 1.2e-3 + + "Or drag coefficient for evaporation over sea" moisture_exchange_sea::NF = 0.9e-3 end @@ -183,8 +214,14 @@ function evaporation!( column::ColumnVariables{NF}, V₀ = column.surface_wind_speed land_fraction = column.land_fraction (;surface_humid) = column - flux_sea = ρ*moisture_exchange_sea*V₀*max(sat_humid_sea - surface_humid,0) - flux_land = ρ*moisture_exchange_land*V₀*α*max(sat_humid_land - surface_humid,0) + + # drag coefficient either from SurfaceEvaporation or from a central drag coefficient + drag_sea, drag_land = evaporation.use_boundary_layer_drag ? + (column.boundary_layer_drag, column.boundary_layer_drag) : + (moisture_exchange_sea, moisture_exchange_land) + + flux_sea = ρ*drag_sea*V₀*max(sat_humid_sea - surface_humid,0) + flux_land = ρ*drag_land*V₀*α*max(sat_humid_land - surface_humid,0) # mix fluxes for fractional land-sea mask land_available = isfinite(skin_temperature_land) && isfinite(α) diff --git a/src/physics/tendencies.jl b/src/physics/tendencies.jl index b003eff93..abfddcdcc 100644 --- a/src/physics/tendencies.jl +++ b/src/physics/tendencies.jl @@ -16,7 +16,6 @@ function parameterization_tendencies!( cos_zenith!(time,model) G = model.geometry - L = model.land_sea_mask rings = eachring(G.Grid,G.nlat_half) @floop for ij in eachgridpoint(diagn) # loop over all horizontal grid points @@ -27,7 +26,7 @@ function parameterization_tendencies!( # extract current column for contiguous memory access reset_column!(column) # set accumulators back to zero for next grid point - get_column!(column,diagn,progn,ij,jring,G,L) + get_column!(column,diagn,progn,ij,jring,model) # execute all parameterizations parameterization_tendencies!(column,model) @@ -45,14 +44,15 @@ function parameterization_tendencies!( # Pre-compute thermodynamic quantities get_thermodynamics!(column,model) + temperature_relaxation!(column,model) + boundary_layer_drag!(column,model) + # VERTICAL DIFFUSION + # diffusion_coefficient!(column,model) + # momentum_diffusion!(column,model) static_energy_diffusion!(column,model) humidity_diffusion!(column,model) - # HELD-SUAREZ - temperature_relaxation!(column,model) - boundary_layer_drag!(column,model) - # Calculate parametrizations (order of execution is important!) convection!(column,model) large_scale_condensation!(column,model) diff --git a/src/physics/thermodynamics.jl b/src/physics/thermodynamics.jl index 56ea895b3..49dea74c8 100644 --- a/src/physics/thermodynamics.jl +++ b/src/physics/thermodynamics.jl @@ -113,6 +113,7 @@ end $(TYPEDSIGNATURES) Calculate the dry static energy for the primitive dry model.""" function get_thermodynamics!(column::ColumnVariables,model::PrimitiveDry) + # use temp here as temp = temp_virt geopotential!(column.geopot,column.temp,model.constants) dry_static_energy!(column, model.constants) end @@ -123,10 +124,11 @@ Calculate thermodynamic quantities like saturation vapour pressure, saturation specific humidity, dry static energy, moist static energy and saturation moist static energy from the prognostic column variables.""" function get_thermodynamics!(column::ColumnVariables,model::PrimitiveWet) - geopotential!(column.geopot,column.temp,model.constants) + geopotential!(column.geopot, column.temp_virt, model.constants, column.surface_geopotential) dry_static_energy!(column, model.constants) saturation_humidity!(column, model.clausis_clapeyron) moist_static_energy!(column, model.clausis_clapeyron) + bulk_richardson!(column, model.constants) # Interpolate certain variables to half-levels vertical_interpolate!(column, model) @@ -170,6 +172,25 @@ function dry_static_energy!(column::ColumnVariables,constants::DynamicsConstants return nothing end +function bulk_richardson!(column::ColumnVariables,constants::DynamicsConstants) + + (;cₚ) = constants + (;u, v, geopot, temp_virt, nlev, bulk_richardson) = column + + V² = u[nlev]^2 + v[nlev]^2 + Θ₀ = cₚ*temp_virt[nlev] + Θ₁ = Θ₀ + geopot[nlev] + bulk_richardson[nlev] = geopot[nlev]*(Θ₁ - Θ₀)/Θ₀/V² + + @inbounds for k in nlev-1:-1:1 + V² = u[k]^2 + v[k]^2 + virtual_dry_static_energy = cₚ * temp_virt[k] + geopot[k] + bulk_richardson[k] = geopot[k]*(virtual_dry_static_energy - Θ₁)/Θ₁/V² + end + + return nothing +end + """$(TYPEDSIGNATURES) Compute the saturation water vapour pressure [Pa], the saturation humidity [kg/kg] and the relative humidity following `clausius_clapeyron`.""" diff --git a/test/column_variables.jl b/test/column_variables.jl index 1d10ef78e..030da6a88 100644 --- a/test/column_variables.jl +++ b/test/column_variables.jl @@ -36,7 +36,7 @@ end column = ColumnVariables{NF}(;nlev) SpeedyWeather.reset_column!(column) - SpeedyWeather.get_column!(column,diagn,progn,1,model.geometry,model.land_sea_mask) + SpeedyWeather.get_column!(column,diagn,progn,1,model) # set a tendency to something humid_tend = rand(NF,nlev)