Merge pull request #404 from SpeedyWeather/mk/humid

Overhauled initial conditions for humidity

src/SpeedyWeather.jl

@@ -96,7 +96,7 @@ export  NoTemperatureRelaxation,
         JablonowskiRelaxation
 
 # EXPORT BOUNDARY LAYER SCHEMES
-export  NoBoundaryLayer,
+export  NoBoundaryLayerDrag,
         LinearDrag,
         QuadraticDrag
 

src/dynamics/initial_conditions.jl

@@ -441,33 +441,30 @@ function initialize_humidity!(  progn::PrognosticVariables,
                                 pres_surf_grid::AbstractGrid,
                                 model::PrimitiveWet)
 
-    # reference saturation water vapour pressure [Pa]
-    # relative humidity reference [1]
-    (;water_pres_ref, relhumid_ref, R_dry, R_vapour, pres_ref) = model.atmosphere
-    gas_ratio = R_dry/R_vapour
-    humid_ref = relhumid_ref*gas_ratio*water_pres_ref   # reference specific humidity [Pa]
+    (;relhumid_ref) = model.atmosphere      # relative humidity reference [1]
 
     # ratio of scale heights [1], scale height [km], scale height for spec humidity [km]     
     (;scale_height, scale_height_humid, σ_tropopause) = model.atmosphere
     scale_height_ratio = scale_height/scale_height_humid
 
     (;nlev, σ_levels_full) = model.geometry
-    n_stratosphere_levels = findfirst(σ->σ>=σ_tropopause,σ_levels_full)
+    n_stratosphere_levels = findlast(σ->σ<=σ_tropopause,σ_levels_full)
 
     # Specific humidity at the surface (grid space)
+    temp_grid = gridded(progn.layers[end].timesteps[1].temp,model.spectral_transform)
     humid_surf_grid = zero(pres_surf_grid)
-    # @. humid_surf_grid = humid_ref*(exp(pres_surf_grid)/(pres_ref*100))^scale_height_ratio
-    q_ref = 1e-3       # kg/kg at the surface
-    @. humid_surf_grid .= q_ref
-    RingGrids.scale_coslat²!(humid_surf_grid)
+    for ij in eachgridpoint(humid_surf_grid)
+        q_sat = saturation_humidity(temp_grid[ij],exp(pres_surf_grid[ij]),model.thermodynamics)
+        humid_surf_grid[ij] = relhumid_ref*q_sat
+    end
 
     humid_surf = spectral(humid_surf_grid,model.spectral_transform)
     spectral_truncation!(humid_surf)
 
     # Specific humidity at tropospheric levels (stratospheric humidity remains zero)
     a = model.spectral_transform.norm_sphere
     for k in n_stratosphere_levels+1:nlev
-        for lm in eachharmonic(humid_surf,progn.layers[1].timesteps[1].humid)
+        for lm in eachharmonic(humid_surf)
             progn.layers[k].timesteps[1].humid[lm] = humid_surf[lm]*σ_levels_full[k]^scale_height_ratio
         end
     end

src/dynamics/models.jl

@@ -159,7 +159,7 @@ Base.@kwdef struct PrimitiveDryModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
 
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
-    boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid)
+    boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
     surface_thermodynamics::AbstractSurfaceThermodynamics{NF} = SurfaceThermodynamicsConstant(spectral_grid)
@@ -248,7 +248,7 @@ Base.@kwdef struct PrimitiveWetModel{NF<:AbstractFloat, D<:AbstractDevice} <: Pr
     # PHYSICS/PARAMETERIZATIONS
     physics::Bool = true
     thermodynamics::Thermodynamics{NF} = Thermodynamics(spectral_grid,atmosphere)
-    boundary_layer_drag::BoundaryLayerDrag{NF} = LinearDrag(spectral_grid)
+    boundary_layer_drag::BoundaryLayerDrag{NF} = NoBoundaryLayerDrag(spectral_grid)
     temperature_relaxation::TemperatureRelaxation{NF} = HeldSuarez(spectral_grid)
     static_energy_diffusion::VerticalDiffusion{NF} = StaticEnergyDiffusion(spectral_grid)
     large_scale_condensation::AbstractCondensation{NF} = SpeedyCondensation(spectral_grid)

src/physics/large_scale_condensation.jl

@@ -42,7 +42,7 @@ function initialize!(scheme::SpeedyCondensation,model::PrimitiveEquation)
         # the relative humidity threshold above which condensation occurs per layer
         σₖ² = σ_levels_full[k]^2
         relative_threshold[k] = threshold_max + threshold_range * (σₖ² - 1)
-        if k >= σ_boundary_layer
+        if σ_levels_full[k] >= σ_boundary_layer
             relative_threshold[k] = max(relative_threshold[k], threshold_boundary_layer)
         end
 
@@ -87,7 +87,7 @@ function large_scale_condensation!(
     ) where NF
 
     (;relative_threshold,humid_tend_max) = scheme
-    time_scale = convert(NF,3600*scheme.time_scale)
+    time_scale = convert(NF,3600*scheme.time_scale)     # [hrs] -> [s]
     n_stratosphere_levels = scheme.n_stratosphere_levels[]
 
     (;humid, pres) = column               # prognostic variables: specific humidity, surface pressure
@@ -101,7 +101,7 @@ function large_scale_condensation!(
     pₛΔt_gρ = pₛ*Δt_sec/gravity/water_density   # precompute constant
 
     (;σ_levels_thick) = geometry
-    latent_heat = convert(NF, atmosphere.latent_heat_condensation/atmosphere.cₚ)
+    latent_heat = convert(NF, atmosphere.latent_heat_condensation)
 
     # 1. Tendencies of humidity and temperature due to large-scale condensation
     @inbounds for k in n_stratosphere_levels+1:nlev   # top to bottom, skip stratospheric levels
@@ -112,7 +112,8 @@ function large_scale_condensation!(
         if humid[k] > humid_threshold
             # accumulate in tendencies (nothing is added if humidity not above threshold)
             humid_tend_k = -(humid[k] - humid_threshold) / time_scale                   # Formula 22
-            temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23
+            # temp_tend[k] += -latent_heat * min(humid_tend_k, humid_tend_max[k]*pres[k]) # Formula 23
+            temp_tend[k] += -latent_heat * humid_tend_k             # Formula 23, without limiter
 
             # If there is large-scale condensation at a level higher (i.e. smaller k) than
             # the cloud-top previously diagnosed due to convection, then increase the cloud-top
@@ -124,7 +125,7 @@ function large_scale_condensation!(
             column.precip_large_scale += -ΔpₖΔt_gρ * humid_tend_k   # Formula 25, unit [m]
 
             # only write into humid_tend now to allow humid_tend != 0 before this scheme is called
-            humid_tend[k] += humid_tend_k                           
+            humid_tend[k] += humid_tend_k
         end
     end
 end

src/physics/thermodynamics.jl

@@ -88,20 +88,16 @@ function saturation_humidity!(
     thermodynamics::Thermodynamics,
 )
     (;sat_humid, sat_vap_pres, pres, temp) = column
-    # (;e₀, T₀, C₁, C₂, T₁, T₂) = thermodynamics.magnus_coefs
     (;mol_ratio) = thermodynamics      # = mol_mass_vapour/mol_mass_dry_air = 0.622
 
     for k in eachlayer(column)
-        # change coefficients for water (temp > T₀) or ice (else)
-        # C, T = temp[k] > T₀ ? (C₁, T₁) : (C₂, T₂)
-        # sat_vap_pres[k] = e₀ * exp(C * (temp[k] - T₀) / (temp[k] - T))
         sat_vap_pres[k] = saturation_vapour_pressure(temp[k],thermodynamics.magnus_coefs)
-        # sat_humid[k] = mol_ratio*sat_vap_pres[k] / (pres[k] - (1-mol_ratio)*sat_vap_pres[k])
         sat_humid[k] = saturation_humidity(sat_vap_pres[k],pres[k];mol_ratio)
     end
 end
 
 function saturation_vapour_pressure(temp::NF,magnus_coefs::MagnusCoefs{NF}) where NF
+    # change coefficients for water (temp > T₀) or ice (else)
     (;e₀, T₀, C₁, C₂, T₁, T₂) = magnus_coefs
     C, T = temp > T₀ ? (C₁, T₁) : (C₂, T₂)
     return e₀ * exp(C * (temp - T₀) / (temp - T))