Update y-axis plot ticks

Author: nichollsh

src/atmosphere.jl

@@ -1283,6 +1283,8 @@ module atmosphere
     Performs hydrostatic integration from the ground upwards.
     Requires density, temperature, pressure to have already been set.
 
+    Does not account for ocean thickness.
+
     Arguments:
         - `atmos::Atmos_t`          the atmosphere struct instance to be used.
 

src/chemistry.jl

@@ -152,13 +152,13 @@ module chemistry
 
             # Super-saturated at the surface...
             elseif dp < 0
-                @debug "Surface $c is supersaturated, dp is negative"
+                @debug @sprintf("            %s super-saturated, partial pressure += %+.3f bar", c, dp/1e5)
 
             # Sub-saturated at the surface...
             #     work out change in partial pressure based on initial reservoir amount
             else
-                @debug "Surface $c is subsaturated, dp is positive"
                 dp = min(dp, atmos.ocean_ini[c] * atmos.grav_surf / (atmos.gas_dat[c].mmw/atmos.layer_μ[end]))
+                @debug @sprintf("            %s sub-saturated, partial pressure += %+.3f bar", c, dp/1e5)
             end
 
             # Record that at least one component has as changed
@@ -167,7 +167,6 @@ module chemistry
             # Reduce or increase total pressure and partial pressure
             atmos.p_boa += dp
             p_gas[c] += dp
-            @debug @sprintf("    adjusted partial prssure by %+.3f bar", dp/1e5)
 
             # Change in surface reservoir, with opposite sign to change in pressure
             atmos.ocean_tot[c] -= (dp / atmos.grav_surf) * (atmos.gas_dat[c].mmw/atmos.layer_μ[end])
@@ -185,7 +184,7 @@ module chemistry
         end
 
         # Generate new pressure grid with updated p_boa
-        @debug @sprintf("New surface pressure: %+.3f bar",atmos.p_boa/1e5)
+        @debug @sprintf("            new p_boa = %.3f bar",atmos.p_boa/1e5)
         atmosphere.generate_pgrid!(atmos)
 
         # Calculate new values for layer properties

src/plotting.jl

@@ -40,6 +40,27 @@ module plotting
         return @sprintf("%d",v)
     end
 
+    # Get y-axis limits
+    function _get_ylims(atmos::atmosphere.Atmos_t)::Tuple
+        return (1e-5*atmos.pl[1]/1.5, 1e-5*max(atmos.p_oboa, atmos.p_boa)*1.5 )
+    end
+
+    # Gey y-axis ticks for log10 axis scale
+    function _get_yticks(atmos::atmosphere.Atmos_t)::Array
+        ylims = _get_ylims(atmos)
+        return 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
+    end
+
+    # Reusable plotting snippets
+    macro _plt_pboa()
+        return esc(:(hline!(plt, [atmos.p_boa/1e5 ], color="black", label=L"p_s",   ls=:solid)))
+    end
+    macro _plt_poboa()
+        return esc(:(hline!(plt, [atmos.p_oboa/1e5], color="black", label=L"p_s^o", ls=:dot)))
+    end
+
+
+
     """
     Plot the temperature-pressure profile.
     """
@@ -48,13 +69,14 @@ module plotting
                             incl_magma::Bool=false,
                             title::String="")
 
-        ylims  = (1e-5*atmos.pl[1]/1.5, 1e-5*atmos.pl[end]*1.5)
-        yticks = 10.0 .^ round.(Int,range(log10(ylims[1]), stop=log10(ylims[2]), step=1))
-
         # Create plot
-        plt = plot(ylims=ylims, yticks=yticks, legend=:outertopright,
+        plt = plot(ylims=_get_ylims(atmos), yticks=_get_yticks(atmos), legend=:outertopright,
                         size=(size_x,size_y); plt_default...)
 
+        # Plot current surface pressure and original
+        @_plt_pboa
+        @_plt_poboa
+
         # Plot phase boundary
         if atmos.condense_any
             sat_t::Array{Float64,1} = zeros(Float64, atmos.nlev_c)
@@ -80,15 +102,15 @@ module plotting
 
         # Plot tmp_magma
         if incl_magma
-            scatter!(plt, [atmos.tmp_magma], [atmos.pl[end]*1e-5],
+            scatter!(plt, [atmos.tmp_magma], [atmos.p_boa/1e5],
                         color="cornflowerblue", label=L"T_m")
         end
 
         # Plot tmp_surf
-        scatter!(plt, [atmos.tmp_surf], [atmos.pl[end]*1e-5], color="brown3", label=L"T_s")
+        scatter!(plt, [atmos.tmp_surf], [atmos.p_boa/1e5], color="brown3", label=L"T_s")
 
         # Plot profile
-        plot!(plt, atmos.tmpl, atmos.pl*1e-5, lc="black", lw=2, label=L"T(p)")
+        plot!(plt, atmos.tmpl, atmos.pl/1e5, lc="black", lw=2, label=L"T(p)")
 
         # Decorate
         xlabel!(plt, "Temperature [K]")
@@ -112,13 +134,15 @@ module plotting
                                 size_x::Int=500, size_y::Int=400,
                                 title::String="")
 
-        ylims  = (1e-5*atmos.pl[1]/1.5, 1e-5*atmos.pl[end]*1.5)
-        yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
-
         # Create plot
-        plt = plot(ylims=ylims, yticks=yticks, legend=:outertopright,
+        plt = plot(ylims=_get_ylims(atmos), yticks=_get_yticks(atmos),
+                        legend=:outertopright,
                         size=(size_x,size_y); plt_default...)
 
+        # Plot current surface pressure and original
+        @_plt_pboa
+        @_plt_poboa
+
         # Plot surface
         scatter!(plt, [atmos.rp*1e-3], [atmos.pl[end]*1e-5], color="brown3", label=L"P_s")
 
@@ -151,12 +175,9 @@ module plotting
         xlims = (-1, 101)
         xticks = collect(range(start=0.0, stop=100.0, step=10.0))
 
-        ylims  = (1e-5*atmos.pl[1]/1.5, 1e-5*atmos.pl[end]*1.5)
-        yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
-
         # Create plot
         plt = plot( xlims=xlims, xticks=xticks,
-                    ylims=ylims, yticks=yticks,
+                    ylims=_get_ylims(atmos), yticks=_get_yticks(atmos),
                     legend=:outertopright, size=(size_x,size_y); plt_default...)
 
         # Temperature profile for reference
@@ -195,13 +216,17 @@ module plotting
         la = 0.7
 
         arr_P = atmos.p .* 1.0e-5 # Convert Pa to bar
-        ylims  = (arr_P[1]/1.5, arr_P[end]*1.5)
-        yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
 
         # Create plot
-        plt = plot(ylims=ylims, yticks=yticks, legend=:outertopright,
+        plt = plot(ylims=_get_ylims(atmos), yticks=_get_yticks(atmos),
+                        legend=:outertopright,
                         size=(size_x,size_y); plt_default...)
 
+
+        # Plot current surface pressure and original
+        @_plt_pboa
+        @_plt_poboa
+
         # Plot log10 VMRs for each gas
         gas_xsurf::Array = zeros(Float64, atmos.gas_num)
         gas::String = ""
@@ -287,19 +312,23 @@ module plotting
                         )
 
         arr_P = atmos.pl .* 1.0e-5 # Convert Pa to bar
-        ylims  = (arr_P[1]/1.5, arr_P[end]*1.5)
-        yticks = 10.0 .^ round.(Int,range( log10(ylims[1]), stop=log10(ylims[2]), step=1))
 
         max_fl = log10(max(100.0, maximum(abs.(atmos.flux_tot)), maximum(atmos.flux_u), maximum(atmos.flux_d)))
         xticks_pos = unique(ceil.(Int, range( start=1, stop=max_fl+1, step=1)))
         xticks = unique(vcat(-1.0.*reverse(xticks_pos), 0.0, xticks_pos))
         xlims = (-xticks_pos[end], xticks_pos[end])
         xticklabels = _intstr.(round.(Int, abs.(xticks)))
 
-        plt = plot(legend=:outertopright, ylims=ylims, yticks=yticks,
+        plt = plot(legend=:outertopright,
+                    ylims=_get_ylims(atmos), yticks=_get_yticks(atmos),
                     xticks=(xticks, xticklabels), xlims=xlims,
                     size=(size_x,size_y); plt_default...)
 
+        # Plot current surface pressure and original
+        @_plt_pboa
+        @_plt_poboa
+
+
         col_r::String = "#c0c0c0"
         col_n::String = "#000000"
         col_c::String = "#6495ed"

src/solver.jl

@@ -332,11 +332,11 @@ module solver
             # Do saturation aloft here, only
             # chemistry.calc_composition!(atmos, rainout, chem, rainout)
             if rainout
-                # reset back to post-chemistry value
+                # reset back to post-chemistry mixing ratios
                 for g in atmos.gas_names
-                    @. atmos.gas_vmr[g] = atmos.gas_cvmr
+                    @. atmos.gas_vmr[g] = atmos.gas_cvmr[g]
                 end
-                _sat_aloft!(atmos)
+                chemistry._sat_aloft!(atmos)
             end
 
             # Calculate fluxes
@@ -621,7 +621,7 @@ module solver
 
                         easy_sf = min(1.0, easy_sf*easy_incr)
                         easy_step = true
-                        @debug "    updated easy_sf = $easy_sf"
+                        @debug "        updated easy_sf = $easy_sf"
 
                         # done modulating
                         if easy_sf > 0.99
@@ -803,7 +803,7 @@ module solver
 
                     # Apply best scale from linesearch
                     ls_alpha = min(max(ls_alpha, ls_min_scale), ls_max_scale)
-                    @debug "    linesearch scale ls_alpha = $ls_alpha"
+                    @debug "            ls_alpha = $ls_alpha"
                     x_dif *= ls_alpha
                 end
 
@@ -930,6 +930,7 @@ module solver
         end
 
         # perform one last evaluation to set `atmos` given the final `x_cur`
+        chemistry.calc_composition!(atmos, rainout, chem, rainout)
         _fev!(x_cur, zeros(Float64, arr_len))
 
         # calc kzz profile
@@ -966,15 +967,20 @@ module solver
         # ----------------------------------------------------------
         loss = maximum(abs.(atmos.flux_tot)) - minimum(abs.(atmos.flux_tot))
         loss_pct = 100.0*loss/maximum(abs.(atmos.flux_tot))
-        @info @sprintf("    outgoing LW flux   = %+.2e W m-2     ", atmos.flux_u_lw[1])
-        if (sol_type == 2)
+        if sol_type == 4
+            @info @sprintf("    outgoing LW flux   = %+.2e W m-2     ", atmos.flux_u_lw[1])
+        end
+        if sol_type == 2
             @info @sprintf("    conduct. skin flux = %+.2e W m-2 ", energy.skin_flux(atmos))
         end
         @info @sprintf("    total flux at TOA  = %+.2e W m-2     ", atmos.flux_tot[1])
         @info @sprintf("    total flux at BOA  = %+.2e W m-2     ", atmos.flux_tot[end])
-        @info @sprintf("    total flux lost    = %+.2e W m-2  (%+.2e %%) ", loss, loss_pct)
-        @info @sprintf("    final cost value   = %+.2e W m-2     ", c_cur)
+        @info @sprintf("    global flux loss   = %+.2e W m-2  (%+.2e %%) ", loss, loss_pct)
+        # @info @sprintf("    final cost value   = %+.2e W m-2     ", c_cur)
         @info @sprintf("    surf temperature   = %-9.3f K        ", atmos.tmp_surf)
+        if rainout
+            @info @sprintf("    surf pressure      = %-9.3e bar      ", atmos.p_boa/1e5)
+        end
 
 
         return atmos.is_converged