;+ ; NAME: ; vt3d_temp_terms_local_iter ; PURPOSE: (one line) ; Return temperature in balance with mean absorbed solar flux ; DESCRIPTION: ; Return temperature in balance with mean absorbed solar flux ; CATEGORY: ; Volatile Transport ; CALLING SEQUENCE: ; temp_terms = vt3d_temp_terms_local_iter(sol_terms,flux_int,emis,freq,therminertia, thermcond, $ ; is_volatile, mass_volatile, specheat_volatile, gravacc, name_species) ; INPUTS: ; sol_terms : absorbed solar flux, complex or real, erg/cm^2/s: ; complex[n_term+1] or complex[n_loc,n_term+1] ; sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) ; + sol_terms[2]*exp(2 *_i freq t) + .. ; ; flux_int : internal hear flux, erg/cm^2/s: float or float[n_loc] ; emis : thermal emissivity, unitless : float or float[n_loc] ; freq : frequency of forcing ; therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc] ; thermcond : thermal conductivity, erg K-1 : float or float[n_loc] ; is_volatile ; 1 if has a volatile, 0 if bare. flag. scalar ; mass_slab ; mass of the volatile slab, g cm^-2, scalar ; specheat_slab ; 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 ; OPTIONAL INPUTS: ; mflux_esc_terms ; escape rate in g cm^-2 s^-1 ; complex[n_term+1] or complex[n_loc,n_term+1] ; defined so that (if _i = sqrt(-1) ) ; mflux_esc = mflux_esc_terms[0] + mflux_esc_terms[1]*exp( _i * freq * t) ; + mflux_esc_terms[2]*exp(2 *_i * freq * t) + .. ; OUTPUTS: ; temp_t: terms of approximate temperture, K, s.t. ; temp = temp_t[0] + temp_t[1]*exp( _i * freq * t + sqrt( _i * zeta) ) ; + temp_t[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + .. ; temp_terms is complex[n_term+1] or complex[n_loc,n_term+1] ; RESTRICTIONS: ; If defined, mflux_esc_terms must be the same dimensions as sol_terms ; PROCEDURE: ; Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D) ; Resubmitted to Icarus. ; Eq. 3.2-10 ; MODIFICATION HISTORY: ; Written 2011 Dec 31, by Leslie Young, SwRI ; 2016 Mar 24 LAY. Modified for inclusion in vt3d library. ;- function vt3d_temp_terms_local_iter, sol_terms,flux_int,emis,freq,therminertia, thermcond, $ is_volatile, mass_volatile, specheat_volatile, gravacc, name_species, mflux_esc_terms _i = complex(0,1) if n_params() lt 12 then mflux_esc_terms = 0. * sol_terms sol_0 = float(extract_last(sol_terms,0,s_dim,n_term)) mflux_esc_0 = float(extract_last(mflux_esc_terms,0,s_dim,n_term)) n_term -= 1 ; temp_0, first without escape, then with temp_0 = vt3d_temp_term0_bare(sol_0, flux_int, emis) vp = bz1980_vp(name_species, temp_0, latheat) temp_0 = vt3d_temp_term0_local(sol_0, flux_int, emis, latheat, mflux_esc_0) ; the Phi's phi_e = vt3d_dfluxdtemp_emit(emis, temp_0) phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia) phi_v = vt3d_dfluxdtemp_slab(freq, mass_volatile, specheat_volatile) frac_varea = 1. ; if not interacting, fraction covered by voaltiles for each location is 100$ phi_a = vt3d_dfluxdtemp_atm(freq, temp_0, frac_varea, gravacc, name_species, latheat) ; latheat is an output ; get dimensionality foo = temp_0 * phi_s t_2d = (size(foo, /n_dim) gt 0) or (s_dim eq 2) ; 1 if temp_terms is a 2-D array n_loc = n_elements(foo) if t_2d then temp_terms = complexarr(n_loc,n_term+1) else temp_terms = complexarr(n_term+1) if t_2d then temp_terms[*,0] = temp_0 else temp_terms[0] = temp_0 for m = 1, n_term do begin sol_m = extract_last(sol_terms,m) mflux_esc_m = extract_last(mflux_esc_terms,m) if t_2d then begin temp_terms[*,m] = (sol_m - latheat*mflux_esc_m) / (phi_e + sqrt(_i * m) * phi_s + (_i*m) * (phi_v + phi_a) * is_volatile ) endif else begin temp_terms[m] = (sol_m - latheat*mflux_esc_m) / (phi_e + sqrt(_i * m) * phi_s + (_i*m) * (phi_v + phi_a) * is_volatile ) endelse endfor dtemp_dzeta = flux_int/thermcond if t_2d then begin if size(dtemp_dzeta,/n_dim) eq 0 then dtemp_dzeta = replicate(dtemp_dzeta, n_loc) if size(emis,/n_dim) eq 0 then em = replicate(emis, n_loc) else em=emis endif in = re(sol_0 + flux_int) n_t = 12 * (n_term + 1) temp_zeta0 = vt3d_thermwave(temp_terms, findgen(n_t)*2*!pi/n_t, dtemp_dzeta, 1., 0.) if not(t_2d) then begin out = mean(emis * !phys.sigma * temp_zeta0^4) i_iter = 0 max_iter = 10 while abs(out-in)/in gt 0.0001 and (i_iter lt max_iter) do begin if temp_0 eq 0. then stop temp_0 = temp_0 + (in-out)/(4*emis * !phys.sigma * temp_0^3) temp_terms[0] = temp_0 phi_e = vt3d_dfluxdtemp_emit(emis, temp_0) phi_a = vt3d_dfluxdtemp_atm(freq, temp_0, frac_varea, gravacc, name_species, latheat) ; latheat is an output for m = 1, n_term do begin ; sol_m is the the m-th term of sol_terms sol_m = extract_last(sol_terms,m) mflux_esc_m = extract_last(mflux_esc_terms,m) temp_terms[m] = (sol_m - latheat*mflux_esc_m) / (phi_e + sqrt(_i * m) * phi_s + (_i*m) * (phi_v + phi_a) * is_volatile ) endfor temp_zeta0 = vt3d_thermwave(temp_terms, findgen(n_t)*2*!pi/n_t, dtemp_dzeta, 1., 0.) out = mean(emis * !phys.sigma * temp_zeta0^4) i_iter += 1 endwhile endif else begin out = em * !phys.sigma * total(temp_zeta0^4,1)/n_t if size(out,/n_dim) eq 0 then stop in = re(sol_0) indx = where(re(sol_0) gt 0., n_indx) for ii = 0, n_indx-1 do begin i = indx[ii] i_iter = 0 max_iter = 10 while abs(out[i]-in[i])/in[i] gt 0.0001 and (i_iter lt max_iter) do begin if temp_0[i] eq 0. then stop temp_0[i] = temp_0[i] + (in[i]-out[i])/(4*em[i] * !phys.sigma * temp_0[i]^3) phi_e[i] = vt3d_dfluxdtemp_emit(em[i], temp_0[i]) phi_a[i] = vt3d_dfluxdtemp_atm(freq, temp_0[i], frac_varea, gravacc, name_species, latheat) ; latheat is an output temp_terms[*,0] = temp_0 for m = 1, n_term do begin sol_m = extract_last(sol_terms,m) mflux_esc_m = extract_last(mflux_esc_terms,m) temp_terms[*,m] = (sol_m - latheat*mflux_esc_m) / (phi_e + sqrt(_i * m) * phi_s + (_i*m) * (phi_v + phi_a) * is_volatile ) endfor temp_zeta0_i = vt3d_thermwave(temp_terms[i,*], findgen(n_t)*2*!pi/n_t, dtemp_dzeta[i], 1., 0.) out[i] = mean(em[i] * !phys.sigma * temp_zeta0_i^4) i_iter += 1 endwhile endfor endelse return, temp_terms end