;+ ; NAME: ; vt3d_thermpro ; PURPOSE: (one line) ; Calculate temperatures on a surface in a style ; similar to thermprojrs ; DESCRIPTION: ; vt3d_thermpro ; CATEGORY: ; VT3D ; CALLING SEQUENCE: ; vt3d_thermpro,tsurf,tod,tdeep,teq,balance,tdarr,zarr0 ; INPUTS: ; rhel heliocentric distance, AU ; alb bolometric albedo ; ti thermal inertia, erg-cgs. Can be an array, in which case it gives the ti of each slab ; rot Rotation period, days ; KEYWORDS: ; emvty emissivity (default 1.0) ; rho density, g cm-3. Can be an array, giving the density of each slab. (default 1.0) ; cp specific heat, erg g-1 K-1 (default H20) ; plat latitude (radians) (default 0.0) ; sunlat Subsolar latitude (radians) (default 0.0) ; heatflow Endogenic heat flow, erg cm-2 s-1 (default 0.0) ; zarr Array of depths to the BASE of each slab (NOT slab thicknesses), in cm. ; First element must be non-zero. ; zeta array to CENTER of each slab ; nslab Number of slabs (ignored if zarr set) (default 30) ; ntinc Time increments per day (default 5000) ; nday Number of days per run (default 4) ; skdepths Model depth, in skindepths (ignored if zarr set) (default 6) ; sol_norm_1au solar constant, erg cm-2 s-1 ; silent If set, sends no output to the terminal (except stability warnings) ; OUTPUT ; tsurf: diurnal surface temperature array ; tod: time of day array (in radians) ; tdeep: final deep temperature ; teq: Instantaneous-equilibrium temperature array ; balance: Thermal balance: (power out - heatflow)/(power in) ; tdarr: temperature/depth array, with 200 diurnal increments ; (set this to a defined variable before entering thermprojrs). ; zarr0: for convenience, returns the array of slab-bottom depths ; (same as zarr) ; MODIFICATION HISTORY: ; Written 2012 Oct 12, by Leslie Young, SwRI ; 2016 Mar 24 LAY. Modified for inclusion in layoung library. ;- pro vt3d_thermpro,tsurf,tod,tdeep,teq,balance,temp_mat,zarr0, $ soltod = soltod, $ ; Required physical keyword parameters rhel=rhel, $ ; heliocentric distance, AU alb=alb, $ ; bolometric albedo ti=ti, $ ; thermal inertia, erg-cgs. Can be an array, in ; which case it gives the ti of each slab rot=rot, $ ; Rotation period, days ; Optional physical keyword parameters: emvty=emvty, $ ; emissivity (default 1.0) rho=rho, $ ; density, g cm-3. Can be an array, giving the density of each slab. (default 1.0) cp=cp, $ ; specific heat, erg g-1 K-1 (default H20) plat=plat, $ ; latitude (radians) (default 0.0) sunlat=sunlat, $ ; Subsolar latitude (radians) (default 0.0) heatflow=heatflow, $ ; Endogenic heat flow, erg cm-2 s-1 (default 0.0) ; Optional simulation keyword parameters: zarr=zarr, $ ; Array of depths to the BASE of each slab (NOT slab thicknesses), in cm. First element must be non-zero. zeta=zeta, $ ; array to CENTER of each slab nslab=nslab, $ ; Number of slabs (ignored if zarr set) (default 30) ntinc=ntinc, $ ; Time increments per day (default 5000) nday=nday, $ ; Number of days per run (default 4) skdepths=skdepths, $ ; Model depth, in skindepths (ignored if zarr set) (default 6) sol_norm_1au = sol_norm_1au, $ ; solar constant, erg cm-2 s-1 ; Optional debugging keywords temp_terms = temp_terms, $ ; Complex ; Silent - last keyword silent=silent ; If set, sends no output to the terminal (except stability warnings) ; 1-dimensional surface thermal model for a rotating body (arbitrary ; insolation dependence with time is also catered for) ; Output variables (Position-dependent: later items in the list can be omitted): ; tsurf: diurnal surface temperature array ; tod: time of day array (in radians) ; tdeep: final deep temperature ; teq: Instantaneous-equilibrium temperature array ; balance: Thermal balance: (power out - heatflow)/(power in) ; tdarr: temperature/depth array, with 200 diurnal increments ; (set this to a defined variable before entering thermprojrs). ; zarr0: for convenience, returns the array of slab-bottom depths ; (same as zarr) ; Required physical keyword parameters ;------------------------- ; this function name ;------------------------- funcname = 'vt3d_thermpro' ;------------------------- ; constants & defaults ;------------------------- physconstants deg = !pi/180. if not keyword_set(sol_norm_1au) then $ sol_norm_1au = 1374.e3 ; solar constant at 1 AU, erg/cm^2/s ; sol_norm_1au = 1367.6e3 ; solar constant at 1 AU, erg/cm^2/s s_per_day = 24.d * 3600.d s_per_year = 365.25d * s_per_day cm_per_au = 149597870.691e5 ; cm per AU km_per_au = 149597870.691 ; km per AU sigma=5.670e-5 ; Stefan/Boltzmann, erg cm-2 s-1 K-4 _i = complex(0,1) if not keyword_set(heatflow) then heatflow=0.0 if not keyword_set(plat) then plat=0.0 if n_elements(sunlat) eq 0 then sunlat = 0. if not keyword_set(emvty) then emvty=1.0 if not keyword_set(rho) then rho=1.0 if not keyword_set(cp) then cp=1.5e7 ; Basalt flows, from Davies 1996 if not keyword_set(nslab) then nslab=45 if not keyword_set(ntinc) then ntinc=5000 if not keyword_set(nday) then nday=4 if not keyword_set(skdepths) then skdepths=9 if not keyword_set(silent) then silent=0 ;------------------------- ; Initialize insolation on an bare location ;------------------------- period = double(rot) * s_per_day lat = plat ; latitude in radians lon = 0. ; east longitude in radians lat_sol = sunlat albedo = alb dist_sol_au = rhel ; absorbed solar at sub-solar (ss) point sol_ss = sol_norm_1au * (1-albedo) / dist_sol_au^2 ; erg/cm^2/s ; Calculate the sinusoidal expansion h_phase0 = - !pi ; hour angle at phase 0 lon_sol_phase0 = lon - h_phase0 flux_sol_t = vty16_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 7, $ sol_norm_1au=sol_norm_1au) ; Calculate the temperature terms freq = 2.d * !dpi / period flux_int = heatflow emis = emvty therminertia = ti dens = rho specheat = cp thermcond = therminertia^2/(dens * specheat) z_skin = vty16_skindepth(dens, specheat, thermcond, freq) ;temp_t = temp_terms_bare(flux_sol_t,flux_int,emis,freq,therminertia) temp_t = vty16_temp_terms_bare_iter(flux_sol_t,flux_int,emis,freq,therminertia, thermcond) ;temp_t = temp_terms_bare_flux(flux_sol_t,flux_int,emis,freq,therminertia, thermcond) ;;print_traceback,/routine,/cr ; %%% ;;help, flux_int,emis,freq,therminertia, thermcond ; %%% ;;print, temp_t ; %%% ;dtemp_dzeta = flux_int / thermcond dtemp_dzeta = -flux_int / (sqrt(freq) * therminertia) ; phase, defined to 0 at start n_per_period = ntinc n_t = nday * n_per_period delta_t = period/n_per_period tau = delta_t * freq ; change in phase per time step, radian phase = dindgen(n_t) * tau phase_tavg = phase + 0.5 * tau temp_zeta0 = vty16_thermwave(temp_t, findgen(n_per_period)*2*!pi/n_per_period, dtemp_dzeta, z_skin, 0.) balance = mean(emis * !phys.sigma * temp_zeta0^4) / $ (sol_ss * vty16_solar_mu(lat, lon, lat_sol, 0., /lonavg) ) del_typ = float(skdepths)/float(nslab) zeta = -findgen(skdepths/del_typ+1) * del_typ temp_wave = vty16_thermwave(temp_t, phase, dtemp_dzeta, z_skin, zeta) temp_phase0 = vty16_thermwave(temp_t, 0., dtemp_dzeta, z_skin, zeta) ;------------------------- ; things that don't change with time ; the interior alpha and beta matrix elements ;------------------------- n_z = n_elements(zeta) del_a = replicate(del_typ, n_z) ;normalized layer width Above, float[n_z] del_b = replicate(del_typ, n_z) ;normalized layer width Below, float[n_z] del = replicate(del_typ, n_z) ; normalized layer width, float[n_z] del[0] = del[0]/2. phi_s = vty16_dfluxdtemp_substrate(freq, therminertia) ; erg cm^-2 s^-1 K^-1 : float[n_loc] gamma_J = flux_int / (del[n_z-1] * phi_s / tau) ; lower boundary alpha_i = vty16_alpha_interior(tau, del, del_a) beta_i = vty16_beta_interior(tau, del, del_b) eta_i = 1 - alpha_i - beta_i ;initialize surface and interior temperatures temp_0 = temp_phase0[0] temp_i = temp_phase0[1:n_z-1] ;------------------------- ; The counter is i_phase ;------------------------- ; assign arrays to save temperatures temp_mat = fltarr(n_t,n_z) ; temp at time_n temp_mat[0,0] = temp_0 temp_mat[0,1:n_z-1] = temp_i sol_arr = dblarr(n_t) flux_sol_tavg_arr = fltarr(n_t) ; insolation from time_n to time_n+delta_t for i_t = 0L, n_t-2 do begin ; 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 phi_e = vty16_dfluxdtemp_emit(emis, temp_0) ; erg cm^-1 s^-1 K^-1 : float[n_loc] ; total derivative of flux w/r.to temperature at the surface ; in this timestep phi_tot = del[0]*phi_s/tau + phi_e/2. ; erg cm^-1 s^-1 K^-1 : float[n_loc] ; the surface beta matrix elements beta_0 = (phi_s/del_b[0])/phi_tot ; float[n_loc] ; the surface gamma matrix element phase_mid = phase_tavg[i_t] ; (i_t + 0.5)*tau sol_tavg = sol_ss * vty16_solar_mu(lat, lon, lat_sol, lon_sol_phase0 - phase_mid) sol_arr[i_t] = sol_tavg flux_sol_tavg_arr[i_t] = sol_tavg ; Forcing array element, K gamma_0 = (sol_tavg - emis*!phys.sigma*( temp_0)^4.)/ phi_tot ; vty16_step_expl_single, alpha_i, beta_0, beta_i, gamma_0, gamma_J, $ temp_0, temp_i ; vtstep_cn_single, alpha_i, beta_0, beta_i, gamma_0, gamma_J, $ ; temp_0, temp_i temp_mat[i_t+1,0] = temp_0 temp_mat[i_t+1,1:n_z-1] = temp_i if ((i_t+1) mod n_per_period) eq 0 then begin balance = mean(emis * !phys.sigma * temp_mat[i_t+1-(n_per_period-1):i_t+1,0]^4) / $ ( sol_ss * vty16_solar_mu(lat, lon, lat_sol, 0., /lonavg) ) endif endfor tsurf = temp_mat[n_t-ntinc:n_t-1,0] ;tdeep = temp_mat[n_t-ntinc:n_t-1, n_z-2] ;tod = lon_sol_phase0 - phase[n_t-ntinc-1:n_t-1] - !pi/2 tod = phase[n_t-ntinc:n_t-1] - (nday - 1) * 2.d*!dpi zarr0 = (zeta - del_typ/2.) * z_skin ;soltod = sol_ss * solar_mu(lat, lon, lat_sol+0*tod, tod+!pi) soltod = sol_arr[n_t-ntinc:n_t-1] teq=((soltod+heatflow)/emvty/sigma)^0.25 avein=total(soltod)/ntinc aveout=(emvty*sigma*total(tsurf^4))/ntinc tmean=total(tsurf)/ntinc tbase = temp_mat[n_t-ntinc:n_t-1, n_z-2] tdeepfinal=total(tbase)/ntinc tdeep = tdeepfinal netpower=aveout-avein balance=(aveout-heatflow)/avein ; make this assignment last temp_terms = temp_t temp_mat = temp_mat[n_t-ntinc:n_t-1,*] ; This next line is a debug, to see the "initial" phase ;temp_mat = thermwave(temp_t, phase[n_t-ntinc:n_t-1], dtemp_dzeta, z_skin, zeta) if not silent then print,'X',avein/1e3,aveout/1e3,balance,netpower/1e3,tmean,tdeepfinal,format="(a4,e11.3,e11.3,f8.4,e11.3,f7.2,f7.2)" end