Merge pull request #64 from exoclime/picket_fence

Picket fence

ifile/picket_fence_wasp43b_hydro.thr

@@ -0,0 +1,348 @@
+# config file for THOR
+# config format version
+config_version = 1
+
+# WASP-43b with grey RT and sponge layer (Mendonca+ 2018)
+
+
+#-- Time stepping and output options -----------------------------------------#
+# number of steps
+num_steps = 288000
+
+# length of timesteps in seconds
+timestep = 300
+
+# output
+# output every n steps
+n_out = 28800
+
+# write to log file at a different cadence from n_out
+custom_log_n_out = false
+# how often to write to log file (if above = true)
+log_n_out = 2880
+
+# output directory (relative to current working directory)
+# defaults to 'results'
+results_path = picket_fence_wasp43b_hydro
+
+# output global diagnostics (energy, mass, momentum, entropy)
+globdiag = true
+
+# output global diagnostics (to text file) at a different cadence from n_out
+custom_global_n_out = false
+# how often to output global diag to text (if above = true)
+global_n_out = 1000
+
+# output time mean quantities (averaged over interval between outputs)
+output_mean = true
+
+
+#-- Planetary parameters -----------------------------------------------------#
+# name of simulation for output files
+simulation_ID = wasp43b_hydro
+
+# Radius [m]
+radius = 72427000.0
+
+# Rotation rate [rad s^-1]
+rotation_rate = 9.09E-5
+
+# Gravitational acceleration [m/s^2]
+gravitation = 47.39
+
+# Gas constant [J/(Kg K)]
+Rd = 3714
+
+# Specific heat capacities [J/(Kg K)]
+Cp = 13000
+
+# Mean atmospheric temperature [K]
+Tmean = 1800.0
+
+# Reference surface pressure [Pa]
+P_ref = 1e8 #1e8
+
+
+#-- Grid options -------------------------------------------------------------#
+# Altitude of the top of the model domain [m]
+Top_altitude = 2.45e6 
+
+# Horizontal resolution level.
+glevel = 5
+
+# Number of vertical layers
+vlevel = 40
+
+# Spring dynamics
+spring_dynamics = true
+
+# Parameter beta for spring dynamics
+spring_beta = 1.15
+
+# use grid refinement in vertical for lower atmos
+vert_refined = false
+
+#lowest layer thickness (used with vert_refined = true) (meters)
+lowest_layer_thickness = 2
+
+#altitude to transition to linear spacing (used with vert_refined=true) (meters)
+transition_altitude = 1000
+
+
+#-- Model options ------------------------------------------------------------#
+# Non-hydrostatic parameter
+NonHydro = false
+
+# Deep atmosphere
+DeepModel = true
+
+# Initial conditions
+rest = true
+
+# initial conditions file, used if rest is set to false
+# (path relative to current working directory)
+# defaults to 'ifile/esp_initial.h5'
+# 'planet' file must be present with name a la 'ifile/esp_initial_planet.h5'
+initial = ifile/esp_initial.h5
+
+# initial temperature-pressure profile (isothermal or guillot)
+# overridden if rest = false !
+# if isothermal: T = Tmean
+# if guillot: equations 27 (Guillot 2010)
+#    with Teq = Tmean and mu* = 0.5
+#    also uses the gray RT params Tint, kappa_lw, kappa_sw, f_lw
+#    includes the collision-induced absorption approx. from Heng+ 2011
+#    this can be omitted by setting f_lw = 1
+
+# Lee comment - generally try guillot first then if that fails try isothermal
+init_PT_profile = parmentier
+
+# Core benchmark tests
+# Held-Suarez test for Earth == HeldSuarez
+# Benchmark test for shallow hot Jupiter == ShallowHotJupiter
+# Benchmark test for deep hot Jupiter == DeepHotJupiter
+# Benchmark test for tidally locked Earth == TidallyLockedEarth
+# Acoustic wave test = AcousticTest
+# Gravity wave test = GWaveTest
+# No benchmark test == NoBenchmark (model is then forced with grey RT by default)
+core_benchmark = NoBenchmark
+
+# Switch off dynamical core (gcm)
+# This is useful for testing/debugging physics modules
+gcm_off = false
+
+# use (dry) convective adjustment scheme
+conv_adj = true
+
+# type of thermodynamic equation = "entropy" (ready) or "energy" (in development!!)
+thermo_equation = entropy
+
+# ultrahot atmosphere options (hydrogen only)
+# under development!!
+# heating due to H/H2 chemisty ('none', 'quasi_eql', or 'relax_chem')
+ultrahot_heating = none
+# variable thermodynamics due to H/H2 chem ('none', 'vary_R_CP', or 'full' )
+ultrahot_thermo = none
+
+## diffusion ############################################
+# Hyper-diffusion
+HyDiff  =    true
+
+# Order of hyperdiffusion operator (4,6,8,...)
+HyDiffOrder = 6
+
+# Divergence-damping
+DivDampP =    true
+
+# Strength of diffusion
+Diffc = 0.0025 #0.015
+
+# Strength of divergence damping
+DivDampc = 0.01 #0.015
+
+# Add vertical terms to hyperdiffusion (6th order)
+VertHyDiff = true
+
+# Strength of vertical hyperdiffusion
+Diffc_v = 0.001
+
+#########################################################
+
+
+#-- Sponge layer (Rayleigh drag or diffusive sponge) ----------------------#
+# use sponge layer (Rayleigh drag) at top of atmosphere?
+RayleighSponge = true
+
+# use temperature sponge layer (Rayleigh drag) at top of atmosphere?
+# (not well tested!)
+RayleighSpongeT = false
+
+# if true, damp to zonal mean (i.e., damp eddy component) (rayleigh sponge)
+# if false, damp to zero
+damp_uv_to_mean = true
+damp_w_to_mean = true
+
+# latitude bins for rayleigh sponge (zonal mean is calculated over these)
+nlat_bins = 20
+
+# bottom of rayleigh sponge layer (fractional height)
+ns_ray_sponge = 0.8
+
+# strength of rayleigh sponge layer (1/damping time)
+# horizontal component
+Ruv_sponge = 1e-3
+# vertical component
+Rw_sponge = 1e-4
+
+# Technical setting: change how rayleigh sponge is applied
+# imp = implicitly update momentum in profx, averages computed once per dt
+# exp1 = explicitly add damping to slow modes, averages computed once per dt
+# exp3 = explicitly add damping to slow modes, averages computed 3x per dt
+raysp_calc_mode = imp
+
+# use diffusive sponge layer at top of atmosphere? (in development!!)
+DiffSponge = true
+
+# strength of diffusive sponge layer (unitless)
+Dv_sponge = 0.01 #0.005
+
+# bottom of diffusive sponge layer (fractional height)
+ns_diff_sponge = 0.8
+
+# order of diffusion operator in diff sponge (2 or 4)
+order_diff_sponge = 4
+
+
+#-- Radiative transfer (gray) options (core_benchmark = NoBenchmark) -------------------#
+## RT parameters #######################################
+radiative_transfer = true
+
+# stellar temperature (K)
+Tstar = 4300 #4520
+
+# stellar metalicity [Fe/H] (for picket-fence scheme essential)
+MetStar =  	0.0
+
+# orbital distance or semi-major axis (au)
+planet_star_dist = 0.01525
+
+# radius of host star (R_sun)
+radius_star = 0.667
+
+# bond albedo of planet
+albedo = 0.18
+
+# diff_ang = 1/(diffusivity factor) for lw radiation
+diff_ang = 0.5
+
+# power law index of unmixed absorbers (lw and sw)
+# optical depth at P_ref is: tau_lw = (kappa_lw/(f_lw*g)) * P_ref
+# n_lw = 2 approximations collision-induced absorption
+# n_lw = 4 approximations water-vapor in Earth troposph.
+n_lw = 2
+
+# optical depth at P_ref is: tau_sw = (kappa_sw/g) * (P_ref)
+n_sw = 1
+
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+# These params are also used if init_PT_profile = guillot
+# grey opacity of thermal wavelengths
+kappa_lw = 0.0025
+
+# grey opacity of incoming stellar flux
+kappa_sw = 0.00125
+
+# strength of unmixed absorbers in lw
+f_lw = 0.5
+
+# temperature of internal heat flux (bottom boundary) (K)
+Tint = 550.87
+#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+# add sin(lat)^2 dependence to tau lw (advanced)
+latf_lw = false
+# opacity at poles (lw)
+kappa_lw_pole = 0.0025
+
+# include surface heating
+surface = false
+# heat capacity of surface
+Csurf = 1e7
+
+# run gray RT code without GCM, every column identical zenith angle (experimental!)
+rt1Dmode = false
+
+#########################################################
+
+## insolation (orbit + spin-state) parameters ###########
+# synchronous rotation (tidally-locking at 1:1)
+sync_rot = true
+
+# mean motion of orbit (if sync_rot=false and ecc>0) (rad/s)
+#mean_motion = 1.98e-7
+
+# initial substellar longitude (deg)
+#alpha_i = 0
+
+# initial orbital position (deg)
+#true_long_i = 0
+
+# eccentricity of orbit
+#ecc = 0
+
+# obliquity (axial-tilt) (deg)
+#obliquity = 0
+
+# longitude of periastron (relative to equinox) (deg)
+# (stupid Earth convention applies)
+#longp = 0
+#########################################################
+# spin up and spin down start and stop. Set to negative value to ignore and run at full power
+# spins up/down with a multiplication factor following a sine from 0 to 1 from
+# start to stop of spin up time
+#dgrt_spinup_start = -1
+#dgrt_spinup_stop = -1
+#dgrt_spindown_start = -1
+#dgrt_spindown_stop = -1
+
+#########################################################
+#-- Boundary layer options (core_benchmark = NoBenchmark) --------------------#
+boundary_layer = false
+
+# type of boundary layer drag ('RayleighHS' or 'LocalMixL')
+bl_type = RayleighHS
+
+# strength of drag (bl_type = RayleighHS)
+#surf_drag = 1.157407e-5
+
+# boundary layer sigma (drag exist at this % of surface pressure and above)
+# (bl_type = RayleighHS)
+#bl_sigma = 0.7
+
+# surface roughness length (meters) (bl_type = LocalMixL)
+#z_rough = 3.21e-5
+
+# asymptotic scale length (ASL) for BL (meters) (bl_type = LocalMixL)
+#abl_asym_len = 150
+
+# asymptotic scale length (ASL) for free atmos (meters) (bl_type = LocalMixL)
+#free_asym_len = 30
+
+# height to transition from BL to free atmos for ASL (meters) (bl_type = LocalMixL)
+# set to -1 to use abl_asym_len in entire atmosphere
+#asl_transition_height = -1
+
+
+#-- Hot Jupiter chemistry ----------------------------------------------------#
+# chemical relaxation model from Mendonca+ 2018(b)
+# contains species CH4, CO, H2O, CO2, NH3
+chemistry = true
+
+# the files below contain the time scales and equilibrium values (from VULCAN)
+chem_time_file = src/physics/modules/ifile/solar_chem_time.txt
+chem_fEQ_file = src/physics/modules/ifile/solar_fEQ_THOR.txt
+
+
+#-- Device options -----------------------------------------------------------#
+# GPU ID number
+GPU_ID_N = 0

mjolnir/hamarr.py

@@ -896,6 +896,7 @@ def regrid(resultsf, simID, ntsi, nts, pgrid_ref='auto', overwrite=False, comp=4
                       'W': (output.Wh[:, :-1, 0] + (output.Wh[:, 1:, 0] - output.Wh[:, :-1, 0]) * interpz[None, :]) / output.Rho[:, :, 0],
                       'Rho': output.Rho[:, :, 0],
                       'Mh': output.Mh[:, :, :, 0],
+                      # 'diffmh': output.diffmh[:, :, :, 0],
                       'Pressure': output.Pressure[:, :, 0],
                       'Rd': output.Rd[:, :, 0],
                       'Cp': output.Cp[:, :, 0],