;+ ; NAME: ; vt3d_1loc_diurnal_local ; PURPOSE: (one line) ; Calculate the diurnal variation at a single non-interacting location ; DESCRIPTION: ; Calculate the diurnal variation at a single non-interacting location ; CATEGORY: ; Volatile Transport ; CALLING SEQUENCE: ; vt3d_1loc_diurnal_bare, constants, input, grid, program, output ; INPUTS: ; constants: constants structure, with ; i ; sqrt(-1) ; deg ; radian per deg ; hour ; s per hour ; year ; s per year ; km ; cm per km ; ti_SI ; [cgs thermal inertia] per [SI thermal inertia] ; m_u ; atomic mass constant, g ; N_a ; Avogadro constant, mol^-1 ; k ; Boltzmann constant, erg K^-1 ; sigma ; Stefan-Boltzmann constant erg cm^-2 K^-4 s^-1 ; G ; gravitational constant, cm^3 g^-1 s^-2 ; sol_norm_1au ; solar flux at 1 AU, erg/cm^2/s ; input : inputs defining the problem ; dist_sol_au ; heliocentric distance in AU ; albedo ; bond albedo, assumed equal to hemispheric ; emis ; emissivity, scalar ; flux_int ; internal heat flux, erg cm^-2 s^-1 ; therminertia ; thermal inertia in erg cm^-2 s^1/2 K^-1, scalar ; dens ; density, g cm^-3, scalar ; specheat ; specific heat in erg g^-1 K^-1, scalar ; thermcond ; thermal conductivity in erg cm^-1 s^-1 K^-1, scalar ; lon_sol_ref ; reference sub-solar longitude, radian ; t_lon_sol_ref ; time of lon_sol_ref, s ; dlon_sol_dt ; rate of change of lon_sol, radian/s ; ; (negative for rotation N pole, E longitude) ; lat_sol ; sub-solar latitude, radian ; is_volatile ; 1 if has a volatile, 0 if bare. flag. scalar ; mass_volatile ; mass of the volatile slab, g cm^-2, scalar ; specheat_volatile ; specific heat of the volatile slab, erg g^-1 K^-1, scalar ; gravacc ; gravitational acceleration, cm s^-2 ; name_species ; name of the species. String. "N2", "CO", or "CH4" ; mflux_esc ; escape rate in g cm^-2 s^-1, constant with time ; grid : grids ; lat ; latitude of spot, radian, scalar ; lon ; longitude of spot, radian, scalar ; t_start ; time of calculation start, s, scalar ; t_end ; time of calculation end, s, scalar ; tau ; unitless timestep, (delta time) * (frequency) ; zeta ; unitless depths of layers, cm, [nz] ; program : program control ; n_terms ; number of terms in the sinusoidal expansion ; iterate_temp_terms ; 1 to iterate the explicit temp terms ; iterate_dfluxdtemp_e ; iterate the flux-per-temperature for emission ; explicit ; 1 if explicit timestep, 0 if implicit (Crank-Nicholson) ; debug ; 1 to print debug statements, 0 otherwise ; OUTPUT: ; output: output ; sol_terms : terms in solar expansion, erg cm^-2 s^-1, complex[nterm+1] ; temp_terms : terms in temperature expansion, K, complex[nterm+1] ; s : array of states ; t : time ; temp_0 : temperature of layer 0, K, scalar ; temp_1J : temperature of layers 1..J, K, [nz] ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D) ; Resubmitted to Icarus. ; Section 3 (is_volatile=0) & 4 (is_volatile=1) ;- pro vt3d_1loc_diurnal_local, c, i, g, p, o freq = abs(i.dlon_sol_dt) ; frequency, Hz period = 2.d*!dpi/freq ; period, s ; ====================================== ; Step 1: get the analytic form ; ; sub-solar longitude at time 0 lon_sol_0 = i.lon_sol_ref + i.dlon_sol_dt * (i.t_lon_sol_ref - g.t_start) ; hour angle at time 0 h_phase0 = (lon_sol_0 - g.lon) * sgn(i.dlon_sol_dt) ; subsolar solar flux, erg cm^-2 s^-1 sol_ss = c.sol_norm_1au * (1-i.albedo) / i.dist_sol_au^2 ; decomposition of insolation (_terms for terms) flux_sol_terms = vt3d_sol_terms_diurnal(i.dist_sol_au, i.albedo, g.lat, h_phase0, i.lat_sol, p.n_terms, $ sol_norm_1au=c.sol_norm_1au) ; decomposition of escape (_terms for terms) ; Only the constant term is non-zero mflux_esc_terms = complexarr(p.n_terms + 1) mflux_esc_terms[0] = i.mflux_esc ; decomposition of temperature case p.iterate_temp_terms of 0: temp_terms = vt3d_temp_terms_local(flux_sol_terms,i.flux_int,i.emis,freq,i.therminertia, $ i.is_volatile, i.mass_volatile, i.specheat_volatile, i.gravacc, i.name_species, mflux_esc_terms) 1: temp_terms = vt3d_temp_terms_local_iter(flux_sol_terms,i.flux_int,i.emis,freq,i.therminertia, i.thermcond, $ i.is_volatile, i.mass_volatile, i.specheat_volatile, i.gravacc, i.name_species, mflux_esc_terms) 2: flux_emis_t = vt3d_eflux_terms_local(flux_sol_terms,i.flux_int,i.emis,freq,i.therminertia, $ i.is_volatile, i.mass_volatile, i.specheat_volatile, i.gravacc, i.name_species, mflux_esc_terms, $ temp_terms=temp_terms) ; flux_emis_t is not used, only the output temp_terms endcase if p.debug then begin print_traceback,'flux_int, emis,freq, therminertia, thermcond & temp_terms', /routine,/cr ; %%% help, i.flux_int,i.emis,freq,i.therminertia, i.thermcond print, temp_terms ; %%% endif temp_terms_0 = temp_terms[0] dtemp_dzeta = -i.flux_int/(sqrt(freq)*i.therminertia) z_skin = vt3d_skindepth(i.dens, i.specheat, i.thermcond, freq) ; ====================================== ; Step 2: get the initial temperature, assuming phase = 0 t_delta = g.tau / freq lon_delta = i.dlon_sol_dt * t_delta n_per_period = round(2.d*!dpi/g.tau) i_t = 0L t = g.t_start + i_t * t_delta lon_sol = lon_sol_0 + i.dlon_sol_dt * t_delta/2.d + (i_t mod n_per_period) * lon_delta phase = 0. n_z = n_elements(g.zeta) temp = vt3d_thermwave_1loc_p0_nz(temp_terms, dtemp_dzeta, g.zeta) temp_0 = temp[0] temp_1J = temp[1:n_z-1]; ; solar longitude lon_sol_beg = lon_sol_0 + (t-g.t_start)*i.dlon_sol_dt lon_sol_end = lon_sol_beg + i.dlon_sol_dt * t_delta lon_sol = lon_sol_beg + i.dlon_sol_dt * t_delta/2 ; absorbed solar flux sol_beg = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol_beg) sol_end = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol_end) sol = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol) s_curr = { $ t:t, $ i_t_mod:0, $ lon_sol_beg:lon_sol_beg, $ lon_sol_end:lon_sol_end, $ lon_sol:lon_sol, $ sol_beg:sol_beg, $ sol_end:sol_end, $ sol : sol, $ temp_0 : temp_0, $ temp_1J : temp_1J $ } ; ====================================== ; Step 3: calculate those things that don't change with time ; Call the "analytic" variables "dfluxdemp" ; and the "numerical" variables "phi" vt3d_zdelta, g.zeta, zeta_delta, zeta_delta_above, zeta_delta_below n_z = n_elements(g.zeta) dfluxdtemp_s = vt3d_dfluxdtemp_substrate(freq, i.therminertia) ; erg cm^-2 s^-1 K^-1 : float dfluxdtemp_v= vt3d_dfluxdtemp_slab(freq, i.mass_volatile, i.specheat_volatile) phi_H_J = zeta_delta[n_z - 1] * dfluxdtemp_s / g.tau gamma_J = i.flux_int / phi_H_J ; lower boundary %%% alpha_1J = vt3d_alpha_interior(g.tau, zeta_delta, zeta_delta_above) ; float[n_z-1] beta_1J = vt3d_beta_interior(g.tau, zeta_delta, zeta_delta_below) ; float[n_z-1] eta_1J = 1 - alpha_1J - beta_1J ; ====================================== ; Step 4: make the room for the output n_t = ceil((g.t_end - t) / t_delta) s = replicate(s_curr, n_t) ; ====================================== ; Step 5: step ; Each loop is from time time_n to time_n + delta_t ; The temperature is defined at the start of the loop at time_n ; The step finds the new temperature at time_n + delta_t while i_t le n_t - 1 do begin ; t = t + t_delta t_beg = s_curr.t ; end of last time step, start of this one t = t_beg + t_delta ; end of this time step t_mid = t - t_delta/2. ; beta_0 and gamma_0 are the only two elements that ; change inside the time loop ; set up the matrix elements dfluxdtemp_e = vt3d_dfluxdtemp_emit(i.emis, temp_0) ; dfluxdtemp_s and dfluxdtemp_v do not change and are computed outside the loop dfluxdtemp_a = vt3d_dfluxdtemp_atm(freq, temp_0, 1.0, i.gravacc, i.name_species, latheat) ; latheat is an output phi_KB_0 = dfluxdtemp_s/zeta_delta_below[0]; Young 2016 3.3-10 phi_H_0 = zeta_delta[0]*dfluxdtemp_s/g.tau ; Young 2016 3.3-5 phi_E = dfluxdtemp_e/2. ; Young 2016 3.3-8 phi_V = (dfluxdtemp_v / g.tau) * i.is_volatile phi_A = (dfluxdtemp_a / g.tau) * i.is_volatile; Young 2016 3.3-5 phi_tot = phi_H_0 + phi_E + phi_V + phi_A; erg cm^-1 s^-1 K^-1 : float[n_loc] beta_0 = phi_KB_0 /phi_tot ; vector elements lon_sol_beg = lon_sol_0 + (t_beg - g.t_start)*i.dlon_sol_dt lon_sol_end = lon_sol_0 + (t - g.t_start)*i.dlon_sol_dt lon_sol = (lon_sol_beg + lon_sol_end)/2. sol_beg = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol_beg) sol_end = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol_end) sol = sol_ss * vt3d_solar_mu(g.lat, g.lon, i.lat_sol, lon_sol) gamma_0 = (sol - c.sigma * i.emis * (temp_0^4.)) / phi_tot ; calculate the new temperatures if p.iterate_dfluxdtemp_e eq 1 then begin temp_0_n = temp_0 temp_1J_n = temp_1J endif vt3d_step_cn_1loc, alpha_1J, beta_0, beta_1J, gamma_0, gamma_J, $ temp_0, temp_1J if p.iterate_dfluxdtemp_e eq 1 then begin temp_0_mid = (temp_0 + temp_0_n)/2. dfluxdtemp_e = dfluxdtemp_emit(i.emis, temp_0_mid) phi_tot = zeta_delta[0]*dfluxdtemp_s/g.tau + dfluxdtemp_e/2. ; erg cm^-1 s^-1 K^-1 : float[n_loc] beta_0 = (dfluxdtemp_s/zeta_delta_below[0]) /phi_tot emission = c.sigma * i.emis * ( (temp_0_mid^4.) - 2 * (temp_0_mid^3) * (temp_0 - temp_0_n)) gamma_0 = (sol - emission) / phi_tot temp_0 = temp_0_n ; reset to temp_1J = temp_1J_n ; previous timestep vt3d_step_cn_1loc, alpha_1J, beta_0, beta_1J, gamma_0, gamma_J, $ temp_0, temp_1J endif ; assign it to the state array s_curr.t = t s_curr.i_t_mod = (i_t mod n_per_period) s_curr.lon_sol_beg = lon_sol_beg s_curr.lon_sol_end = lon_sol_end s_curr.lon_sol = lon_sol s_curr.sol_beg = sol_beg s_curr.sol_end = sol_end s_curr.sol = sol s_curr.temp_0 = temp_0 s_curr.temp_1J = temp_1J s[i_t] = s_curr i_t = i_t + 1 endwhile o = {$ flux_sol_terms:flux_sol_terms, $ temp_terms:temp_terms, $ dtemp_dzeta:dtemp_dzeta, $ s : s $ ; array of states } end