;+ ; NAME: ; vty16_fig5_7_func ; PURPOSE: (one line) ; Make the data for a Pluto movie for Young 2013 ApJL resubmission ; DESCRIPTION: ; Make the data for a Pluto movie for Young 2013 ApJL resubmission ; CATEGORY: ; VT ; CALLING SEQUENCE: ; res = vty16_fig5_7_func(run, av, ev, as, es, ti, mvtot, n_off, res_all) ; INPUTS: ; run - name of this run. If "", then construct from parameters ; av - hemispheric albedo of the volatile ; ev - emissivity of the volatile ; as - hemispheric albedo of the substrate ; es - emissivity of the substrate ; ti - thermal inertia of the substrate, cgs units ; mvtot - total inventory of volatiles ; n_off - number of time periods to offset start by. ; If n_off = 0, start at aphelion. ; If n_off > 0, start n_off timesteps before aphelion ; OPTIONAL INPUT PARAMETERS: ; none ; KEYWORD INPUT PARAMETERS: ; none ; KEYWORD OUTPUT PARAMETERS: ; none ; OUTPUTS: ; res: structure with the results of the run, suitable for plotting. ; ALBEDO_COMP FLOAT Array[3] albedo of composition types, unitless ; EMIS_COMP FLOAT Array[3] albedo of composition types, unitless ; ; AV FLOAT albedo of volatile, unitless ; EV FLOAT emissivity of volatile, unitless ; AS FLOAT albedo of substrate, unitless ; ES FLOAT emissivity of substrate, unitless ; THERMINERTIA DOUBLE Array[60, 19] ; thermal inertia, [loc, depth], cgs ; MVTOT INT 32 total volatile, g/cm^2 ; LAT DOUBLE Array[60] latitude of locations, radian ; N_LAT_EDGE LONG 59 number of latitude edges ; N_TIME_PER_PERIOD ; LONG 240 timesteps per period ; N_TIME LONG 720 number of timesteps ; A FLOAT Array[240] geometric albedo, [time] ; TV FLOAT Array[240] temperature of volatile, K, [time] ; TS FLOAT Array[60, 240] temperature of layer 0, K, [loc,time] ; MS FLOAT Array[60, 240] ; mass of slab, g/cm^2, [loc, time] ; V FLOAT Array[59, 240] ; velocity across lat boundary, cm/s, [lat_edge] ; VS FLOAT Array[59, 240] ; sound speed at lat boundary, cm/s, [lat_edge] ; P FLOAT Array[240] pressure, microbar, [time] ; YR DOUBLE Array[240] time, year, [time] ; LAT0 FLOAT Array[240] subsolar lat, radian, [time] ; R FLOAT Array[240] helioc. dist, AU, [time] ; ANG FLOAT Array[240] true anomoly, radian, [time] ; ISV INT Array[60, 240] ; is volatile-covered, boolean, [loc, time] ; res_all: structure with the results of the run, suitable for full output. ; All the tags in res, PLUS ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; RESTRICTIONS: ; None ; PROCEDURE: ; MODIFICATION HISTORY: ; Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro) ; Oct 21 2012. Make function (pluto_mssearch_func.pro) ; Dec 13 2012. Add start phase as options ; 2016 Apr 06 LAY. Modified for inclusion in vty16 library. ;- function vty16_fig5_7_func, run, av, ev, as, es, ti, mvtot, n_off, res_all ;------------------------- ; constants and flags ;------------------------- physconstants deg = !pi/180. day = 24. * 3600. year = 365.25 * day km = 1e5 ; cm per km au = 149597870.691e5 ; cm per AU sol_norm_1au = 1367.6e3 ; solar constant at 1 AU, erg/cm^2/s time_day = 6.3865509d * day ; %%% time_year = 247.92065d * year ; %%% gravacc = 80. ; cm/s^2 radius = 1163.e5 ; Sicardy et al 2011 molweight = 28. ; molecular weight of N2 specheat_ratio = 1.4 ; usually called gamma. cp/cv for gas crosssection = 0.43e-14 ; cm^2. Tables of physical and chemical constants, Longman, London, 1973 _i = complex(0,1) ; _i = sqrt(-1) flag_volatileslab = 1 ; Should be 1 to include specific heat of volatile slab flag_atm = 1 ; Should be 1 to include atmospheric "breathting" flag_stepscheme = 0 ; 0 for forward, 1 for C-N, 2 for backward name_stepscheme = ['Forward', 'Crank-Nicholson', 'Backward'] flag_volatile_init = 'global' ; %%% ; set the run name if strlen(run) eq 0 then begin run = 'pluto'+ $ '_' + strtrim(string(av,fo='(F4.2)'),2) + $ '_' + strtrim(string(ev,fo='(F4.2)'),2) + $ '_' + strtrim(string(as,fo='(F4.2)'),2) + $ '_' + strtrim(string(es,fo='(F4.2)'),2) + $ '_' + strtrim(string(ti,fo='(E8.2)'),2) + $ '_' + strtrim(string(mvtot,fo='(F5.1)'),2) + $ '_' + strtrim(string(n_off,fo='(I03)'),2) endif jd_peri = 2447823.7 ; Julian date of perihelion name_species= 'N2' ;------------------------- ; surface grid ;------------------------- ; define n_loc - number of locations ; lon double[n_loc] - longitude, radian ; lat double[n_loc] - longitude, radian ; angarea_delta double[n_loc] - angular area of each location, ster name_loc = 'latband_60' name_loc2map = name_loc+'->map_60_30' locgrid, name_loc, n_loc, lon, lat, angarea_delta ;------- ; Set up latitude bands for winds n_lat_edge = n_loc - 1 ; number of latitude-band boundaries ; lat_edge = latitude of latitude-band boundaries lat_edge = (lat[0:n_lat_edge-1] + lat[1:n_lat_edge])/2. ; indx_lat_edge ; lat[0..indx_lat_edge[m]] are all less than lat_edge[m] indx_lat_edge = indgen(n_lat_edge) ;------------------------- ; surface quantities ;------------------------- ; volatile and composition ; is_volatile byte[n_loc] - 1 if the location is volatile covered case flag_volatile_init of ; 2011 conditions -- Southward (IAU) of -20 and northward of 60 '2011': is_volatile = ( (lat le -20*deg) or (lat ge 60*deg) ) ; aphelion conditions - northward (old IAU) of +20 'aphelion' : is_volatile = (lat ge 20*deg) ; global 'global' : is_volatile = (lat ge -100*deg) endcase indx_volatile = where(is_volatile, nindx_volatile, compl=indx_not_volatile, ncompl=nindx_not_volatile) name_comp = ['','volatile','substrate'] n_comp = n_elements(name_comp) indxcomp = 2-is_volatile ; is_volatile=1 -> indx=1, is_volatile=0 -> indx=2 ; albedo. volatile av=albedo of volatile, substrate as = albedo of substrate albedo_comp = [0,av,as] ; %%% albedo_geom_comp = albedo_comp * [1,1/.9, 1/.4] albedo = albedo_comp[indxcomp] albedo_geom = albedo_geom_comp[indxcomp] ; emissivity emis_comp = [1., ev, es] ; %%%% emis = emis_comp[indxcomp] ; mass of volatile slab mass_slab_inventory = mvtot ; %%%%% mass_slab_avg = mass_slab_inventory * 4 * !pi / total(angarea_delta * is_volatile) ;print, run, ': mass_slab_avg ', mass_slab_avg, 'g', fo=formf mass_slab = mass_slab_avg * is_volatile ; volatile specific heat specheat_volatile = 1.3e7 ; %%%% ; beta fraction frac_beta = 1. ; internal heat flux eflux_int = 6. ; %%%% ;------------------------- ; period ;------------------------- time_period = time_year freq = 2.*!pi/(time_period) ;------------------------- ; thermo-physical params ;------------------------- dens=0.93 ; %%%%% ; Assume the skin depth is 1500 cm (15 m) and scale from there, using ; Z = ti / (dens specheat sqrt(freq)) specheat = (ti)/(1500.*sqrt(freq))/dens ; %%%%% thermcond = (ti^2)/(dens * specheat) ; %%%%% therminertia = sqrt(dens * specheat * thermcond) skindepth = vt3d_skindepth(dens, specheat, thermcond, freq) dens_comp = replicate(dens,n_comp) specheat_comp = replicate(specheat,n_comp) thermcond_comp = replicate(thermcond,n_comp) therminertia_comp = replicate(therminertia,n_comp) skindepth_comp = replicate(skindepth,n_comp) ;------------------------- ; a word about naming convension ;------------------------- ; in general, start with the quantity ; eg ; eflux_sol for energy flux, solar ; temp for temperature ; mu0 for cos(illumination angle) ; modify with units, if not cgs ; eg year in phase_year, au in dist_sol_au ; Other modifiers include: ; app : approximate (calculations for the sinusoidal expansion) ; ftc : Fourier Transform, complex ; ftr : Fourier Transform, real ; (fti, fta, ftp = Fourier Transform, imaginary, amplitude, phase, ; but theres aren't actually used) ; time : calculated at each n_time ; *app*volatile - average over the volatiles ; *app*tavg - average over time ;------------------------- ; approximate solar forcing ;------------------------- time_delta = time_period/240. n_time = 240 ; number of timesteps phase_delta = 2.*!pi/n_time ; change in phase per timestep (radian) phase_app = (indgen(n_time)) * phase_delta ; list of phases jd0 = jd_peri ; first julian date ; et = ephemeris time (aka Barycentric Dynamical Time in seconds since ; J2000. See SPICE documentation for details. double[n_time], units seconds et = (jd0 -2451545.0000d)*day + phase_app * time_period / (2.*!pi) - time_period /2. ; phase_year double[n_time], units year. phase_year = (et / year) + 2000. ; phys is an array of structures, n_time long ; The relevant fields are: ; soldist - heliocentric distance in km ; subsollat - subsolar latitude in radian phys = pluto_phys(9,et) dist_sol_au = phys.soldist*km/au ; double[n_time], AU lat0 = phys.subsollat ; double[n_time], radian ; since we're doing diurnal average, the lon0 doesn't matter ; as long as it's the right length lon0 = replicate(180., n_time) ; float[n_time], radian ; Get the diurnal average of the solar incidence angle ; FIRST for each latitude band ; eflux_sol_app_time : double[n_loc, n_time], erg/(cm^2 s) ; eflux_sol_app_ftc : complex[n_loc, n_freq+1],erg/(cm^2 s) ; eflux_sol_app_fta : float[n_loc, n_freq+1], erg/(cm^2 s) ; eflux_sol_app_ftr : float[n_loc, n_freq+1], erg/(cm^2 s) ; eflux_sol_app_tavg : float[n_loc], erg/(cm^2 s) mu0_time = vt3d_solar_mu(lat,lon,lat0, lon0, /lon) ; double[n_loc, n_time], unitless ; diurnal average of the absorbed solar flux ; for each latitude band and each time eflux_sol_app_time = sol_norm_1au * ((1-albedo) # (1/dist_sol_au^2)) * mu0_time ; Decompose the solar forcing n_freq = 2 eflux_sol_app_ftc = complexarr(n_loc, n_freq+1) ; complex[n_loc, n_freq+1], erg/(cm^2 s) ; "sft" = "Slow Fourier Transform" for i_loc = 0, n_loc-1 do eflux_sol_app_ftc[i_loc,*] = sft(phase_app,eflux_sol_app_time[i_loc,*], n_freq) ; c2r outputs the phase (/ophase is set), but only the real is used ; eflux_sol_app_ftr : float[n_loc, n_freq+1], erg/(cm^2 s) eflux_sol_app_ftp = c2r(eflux_sol_app_ftc, a=eflux_sol_app_fta,r=eflux_sol_app_ftr,i=eflux_sol_app_fti,/ophase) eflux_sol_app_tavg = eflux_sol_app_fta[*,0] ; NEXT for averages over the volatiles ; This might appear needlessly complicated, but it is here for the ; legecy of different values of flag_volatile_init, ; for which there are areas that are initially bare ; eflux_sol_app_time_volatile : double[n_time], erg/(cm^2 s) ; eflux_sol_app_ftc_volatile : complex[n_freq+1],erg/(cm^2 s) ; eflux_sol_app_fta_volatile : float[n_freq+1], erg/(cm^2 s) ; eflux_sol_app_tavg_volatile : float, erg/(cm^2 s) eflux_sol_app_time_volatile = vt3d_locavg(eflux_sol_app_time, angarea_delta, indx=indx_volatile) eflux_sol_app_ftc_volatile = sft(phase_app,eflux_sol_app_time_volatile, n_freq) eflux_sol_app_ftp_volatile = c2r(eflux_sol_app_ftc_volatile, a=eflux_sol_app_fta_volatile,r=eflux_sol_app_ftr_volatile,$ i=eflux_sol_app_fti_volatile,/op) eflux_sol_app_tavg_volatile = eflux_sol_app_fta_volatile[0] ;------------------------- ; approximate thermal response ;------------------------- ; FIRST for each latitude band ; temp_surf_app_tavg : float[n_loc], K ; temp_surf_app_ftc : complex[n_loc, n_freq+1], K ; temp_surf_app_ftr : float[n_loc, n_freq+1], K temp_surf_app_tavg = ( (eflux_sol_app_tavg + eflux_int)/(!phys.sigma * emis) ) ^ 0.25 ; Calculate the thermal parameters for the substrate and volatile slab ; for those locations where approximate average temperature is positive indx_tpos = where(temp_surf_app_tavg gt 0) thermparam_substrate_tavg = fltarr(n_loc) thermparam_substrate_tavg[indx_tpos] = therminertia[indx_tpos] * sqrt(freq) / $ (!phys.sigma * emis[indx_tpos] * temp_surf_app_tavg[indx_tpos]^3) thermparam_volatileslab_tavg = fltarr(n_loc) thermparam_volatileslab_tavg = freq * mass_slab[indx_tpos] * specheat_volatile[indx_tpos] / $ (!phys.sigma * emis[indx_tpos] * temp_surf_app_tavg[indx_tpos]^3) temp_surf_app_ftc = complexarr(n_loc, n_freq+1) temp_surf_app_ftc[0] = temp_surf_app_tavg for i_freq = 1, n_freq do begin temp_surf_app_ftc[indx_tpos,i_freq] = (eflux_sol_app_ftc[indx_tpos,i_freq]/$ (!phys.sigma*emis[indx_tpos]*temp_surf_app_tavg[indx_tpos]^3)) / $ (4. + sqrt(_i*i_freq)*thermparam_substrate_tavg[indx_tpos] + $ (_i*i_freq)*(thermparam_volatileslab_tavg[indx_tpos] * flag_volatileslab) ) endfor temp_surf_app_ftp = c2r(temp_surf_app_ftc, a=temp_surf_app_fta,r=temp_surf_app_ftr,i=temp_surf_app_fti,/op) ; NEXT for averages over the volatiles ; temp_volatile_app_tavg : float, K ; temp_volatile_app_ftc : complex[n_freq+1], K ; temp_volatile_app_ftr : float[n_freq+1], K emis_volatile = vt3d_locavg(emis, angarea_delta, indx=indx_volatile) therminertia_volatile = vt3d_locavg(therminertia, angarea_delta, indx=indx_volatile) mass_slab_volatile = vt3d_locavg(mass_slab, angarea_delta, indx=indx_volatile) eflux_int_volatile = vt3d_locavg(eflux_int, angarea_delta, indx=indx_volatile) temp_volatile_app_tavg = ( (eflux_sol_app_tavg_volatile + eflux_int_volatile)/(!phys.sigma * emis_volatile) ) ^ 0.25 thermparam_substrate_tavg_volatile = therminertia * sqrt(freq) / (!phys.sigma * emis_volatile * temp_volatile_app_tavg^3) thermparam_volatileslab_tavg_volatile = freq * mass_slab_volatile * specheat_volatile / (!phys.sigma * emis_volatile * temp_volatile_app_tavg^3) pres_app_tavg = n2vp(temp_volatile_app_tavg) latheat_app_tavg = n2latentheat(temp_volatile_app_tavg) thermparam_atm_tavg_volatile = vt3d_dfluxdtemp_atm(freq, temp_volatile_app_tavg, 4*!pi, gravacc, "N2", latheat_app_tavg)/$ (!phys.sigma * emis_volatile * temp_volatile_app_tavg^3) temp_volatile_app_ftc = complexarr(n_freq+1) temp_volatile_app_ftc[0] = temp_volatile_app_tavg for i_freq = 1, n_freq do begin temp_volatile_app_ftc[i_freq] = (eflux_sol_app_ftc_volatile[i_freq]/(!phys.sigma*emis_volatile*temp_volatile_app_tavg^3)) / $ (4. + sqrt(_i*i_freq)*thermparam_substrate_tavg_volatile + $ (_i*i_freq)*(thermparam_volatileslab_tavg_volatile * flag_volatileslab + $ thermparam_atm_tavg_volatile * flag_atm) ) endfor temp_volatile_app_ftp = c2r(temp_volatile_app_ftc, a=temp_volatile_app_fta,r=temp_volatile_app_ftr,i=temp_volatile_app_fti,/op) ;------------------------- ; zgrid ;------------------------- z = vt3d_zgrid(replicate(skindepth,n_loc),z_delta,n_z, /med) ;------------------------- ; promote the dimension of the thermophysical params ;------------------------- therminertia = loczmat(therminertia, n_loc, n_z) thermcond = loczmat(thermcond, n_loc, n_z) dens = loczmat(dens, n_loc, n_z) specheat = loczmat(specheat, n_loc, n_z) eflux_int = locarr(eflux_int, n_loc) ;------------------------- ; Set up for timestepping - the times, phases and years ;------------------------- ; n_time_per_period : number of timesteps per period (long) n_time_per_period = 240L ; %%% ; n_period : number of complete periods to run (long) n_period = 3L ; %%% ; n_time : number of timesteps, including the optional offset from aphelion n_time = n_time_per_period * n_period + n_off ; ### time_delta = time_period/n_time_per_period ; float, s phase_delta = 2.*!pi*time_delta/time_period ; float, radian ; phase_time : double[n_time], radian ; phase_year : double[n_time], year ; phase_time is used to initilize the temperature ; phase_year is used for plotting, and is not used again in this calculation phase_time = ( (findgen(n_time) + n_time_per_period - n_off) * phase_delta) mod (2. * !pi) phase_year = (1. / year)*((jd0 -2451545.0000d)*day + phase_time * time_period / (2.*!pi) - time_period /2.) + 2000. ;------------------------- ; analytic temp vs. time, and initial temp ;------------------------- ; temp_anal_time[n_loc, n_z, n_time] is the analytic approximation to ; the temperature over the previous period. It is used to calculate ; the initial temperature, temp[n_loc, n_z]. ; We could have jumped straight to temp. Going through temp_anal_time ; is a vestige of the diagnostics in the development of this code. temp_anal_time = fltarr(n_loc, n_z, n_time) temp_volatile_anal_time = fltarr(n_time) ; bare areas use temp_surf_app_ftr and temp_surf_app_ftv indx = indx_not_volatile if indexcount(indx) gt 0 then begin for i_time = 0, n_time-1 do begin w = exp(_i*phase_time[i_time]+sqrt(_i*freq)*(therminertia[indx,*])*z[indx,*]/thermcond[indx,*]) temp_anal_time[indx,*,i_time] = -((eflux_int[indx]#onev(n_z))/thermcond[indx,*])*z[indx,*] $ + temp_surf_app_ftr[indx,0]#onev(n_z) + $ float((temp_surf_app_ftc[indx,1]#onev(n_z) )*w) endfor endif ; volatile areas use temp_volatile_app_ftr and temp_volatile_app_ftv indx = indx_volatile if indexcount(indx) gt 0 then begin for i_time = 0, n_time-1 do begin w = exp(_i*phase_time[i_time]+sqrt(_i*freq)*(therminertia[indx,*])*z[indx,*]/thermcond[indx,*]) temp_anal_time[indx,*,i_time] = -((eflux_int[indx]#onev(n_z))/thermcond[indx,*])*z[indx,*] $ + temp_volatile_app_ftr[0] + $ float(temp_volatile_app_ftc[1]*w) temp_volatile_anal_time[i_time] = temp_anal_time[indx[0],0,i_time] endfor endif ; temp is the temperature at the start of the time step temp = temp_anal_time[*,*,0] ; float[n_loc, n_z] ;------------------------- ; things we want to record and output ;------------------------- time = fltarr(n_time) temp_surf_time = fltarr(n_loc, n_time) temp_time = fltarr(n_loc, n_z, n_time) temp_volatile_time = fltarr(n_time) eflux_sol_time = fltarr(n_loc, n_time) is_volatile_time = intarr(n_loc, n_time) is_xport_time = intarr(n_loc, n_time) mflux_time = fltarr(n_loc, n_time) mass_slab_time = fltarr(n_loc, n_time) trueanom_time = fltarr(n_time) dist_sol_au_time = fltarr(n_time) lat0_time = fltarr(n_time) albedo_geom_time = fltarr(n_time) vel_time = fltarr(n_lat_edge, n_time) vel_sound_time = fltarr(n_lat_edge, n_time) z_atm_time = fltarr(n_lat_edge, n_time) z_mfp_time = fltarr(n_lat_edge, n_time) ;------------------------- ; tri-diag matrices ;------------------------- ; The different layer depths ; z_delta_top is Delta_A in Young 2016 ; z_delta_bot is Delta_B in Young 2016 z_delta_top = fltarr(n_loc,n_z) z_delta_top[*,1] = z_delta[*,1]/2 + (thermcond[*,1:n_z-1]/thermcond[*,0:n_z-2])*z_delta[*,0] z_delta_top[*,2:n_z-1] = z_delta[*,2:n_z-1]/2 + (thermcond[*,2:n_z-1]/thermcond[*,1:n_z-2])*z_delta[*,1:n_z-2]/2 z_delta_bot = fltarr(n_loc,n_z) z_delta_bot[*,0] = z_delta[*,0] + (thermcond[*,0]/thermcond[*,1])*z_delta[*,1]/2 z_delta_bot[*,1:n_z-2] = z_delta[*,1:n_z-2]/2 + (thermcond[*,1:n_z-2]/thermcond[*,2:n_z-1])*z_delta[*,2:n_z-1]/2 ; The elements of the tridiagonal ; alpha_top is alpha in Young 2016 ; alpha_bot is beta in Young 2016 ; alpha_mid equals 1 - eta in Young 2016 diffusion = thermcond/(dens * specheat) alpha_top = fltarr(n_loc,n_z) alpha_bot = fltarr(n_loc,n_z) alpha_top[*,1:n_z-1] = diffusion[*,1:n_z-1] * time_delta / (z_delta_top[*,1:n_z-1]*z_delta[*,1:n_z-1]) alpha_bot[*,0:n_z-2] = diffusion[*,0:n_z-2] * time_delta / (z_delta_bot[*,0:n_z-2]*z_delta[*,0:n_z-2]) alpha_mid = (alpha_top + alpha_bot)/2 ; beta is gamma in Young 2016 beta = fltarr(n_loc, n_z) ; temp_next is the temperature at the end of the time step temp_next = fltarr(n_loc, n_z) for i_time = 0, n_time-1 do begin ; vvvvvvvvvvvvv ; the phase (0 to pi) at start of time step phase = phase_time[i_time] ; calculate solar fluxes at each location %%%% ; for start of time step, eflux_sol_this double[n_loc] ; and end of time step eflux_sol_next double[n_loc] et = (jd0 -2451545.0000d)*day + phase * time_period / (2.*!pi) - time_period /2. time[i_time] = et phys = pluto_phys(9,et,trueanom=trueanom) dist_sol_au = phys.soldist*km/au trueanom_time[i_time] = trueanom dist_sol_au_time[i_time] = dist_sol_au lat0 = phys.subsollat lat0_time[i_time] = lat0 lon0 = 180*deg solar_mu = vt3d_solar_mu(lat,lon,lat0, phase, /lon) albedo_geom_time[i_time] = total(angarea_delta*solar_mu*albedo_geom) / !pi eflux_sol_this = sol_norm_1au * ((1-albedo)/dist_sol_au^2) * solar_mu eflux_sol_next = sol_norm_1au * ((1-albedo)/dist_sol_au^2) * $ vt3d_solar_mu(lat,lon,lat0, phase+phase_delta, /lon) ; store values at start of time step into arrays for output eflux_sol_time[*,i_time] = eflux_sol_this mass_slab_time[*,i_time] = mass_slab is_volatile_time[*,i_time] = is_volatile ; thermal emission temp_surf = temp[*,0] ; surface temperature at start of time step, float temp_surf_time[*,i_time] = temp_surf ; save it for output temp_time[*,*,i_time] = temp ; save it for output is_volatile_time[*,i_time] = is_volatile vty16_fig5_7_mat, flag_volatileslab, time_delta, n_loc, n_z, $ emis, temp_surf, (eflux_sol_this+eflux_sol_next)/2, mass_slab, specheat_volatile, $ z_delta, z_delta_bot, dens, specheat, thermcond, eflux_int, $ beta, alpha_bot, alpha_mid, denom, eflux_net iter = 0 is_xport = is_volatile repeat begin vty16_fig5_7_timestep, $ flag_atm, freq, time_delta, gravacc, name_species, flag_stepscheme, n_loc, lat, n_z, z_delta,z_delta_bot, $ is_xport, angarea_delta,$ emis, temp_surf, (eflux_sol_this+eflux_sol_next)/2, mass_slab, specheat_volatile, $ dens,specheat,thermcond, $ alpha_top, alpha_mid,alpha_bot,beta, denom, eflux_net,$ temp, temp_volatile, $ temp_next, temp_volatile_next, angarea_atm, mflux temp_volatile_time[i_time] = temp_volatile ; velocities for i_edge = 0,n_lat_edge-1 do begin ; %%%%%%% KLUDGE FOR LAT BAND - ; adjsut mflux so that it's zero everywhere indx_mf = where(mflux ne 0) mfluxnet = mflux mfluxnet[indx_mf] = mflux[indx_mf] - total(mflux[indx_mf] * angarea_delta[indx_mf])/total(angarea_delta[indx_mf]) temp_atm = (temp_next[i_edge,0] + temp_next[i_edge+1,0]) / 2 if temp_atm lt 3. then stop z_atm = !phys.k * temp_atm / (molweight * !phys.m_u * gravacc) ; scale height numdens_atm = n2vp(temp_atm)/(!phys.k * temp_atm) dens_atm = molweight * !phys.m_u * numdens_atm vel_time[i_edge, i_time] = -total(mfluxnet[0:i_edge]*angarea_delta[0:i_edge]) * radius/$ (2*!pi*z_atm*dens_atm*cos(lat_edge[i_edge])) vel_sound = sqrt(specheat_ratio * !phys.k * temp_atm / (molweight * !phys.m_u)) vel_sound_time[i_edge, i_time] = vel_sound z_atm_time[i_edge, i_time] = z_atm z_mfp_time[i_edge, i_time] = 1/(numdens_atm * crosssection) endfor ; no velocity beyond where there's no xtransport indx_xport = where(is_xport, nindx_xport) lat_minmax = minmax(lat[indx_xport]) indx_noatm = where(lat_edge lt lat_minmax[0] or lat_edge gt lat_minmax[1], nindx_noatm) if nindx_noatm gt 0 then vel_time[indx_noatm,i_time] = 0 mach = abs(vel_time[*,i_time] / vel_sound_time[*,i_time]) if max(mach) gt 1 then begin indx_xport = where(is_xport, nindx_xport, compl=indx_not_volatile, ncompl=nindx_not_volatile) mflux_max = max(mfluxnet[indx_xport], ii) is_xport[indx_xport[ii]] = 0 indx_xport = where(is_xport, nindx_xport, compl=indx_not_volatile, ncompl=nindx_not_volatile) endif iter += 1 endrep until max(mach) lt 1 or iter gt 50 or nindx_xport le 2 ;stop is_xport_time[*,i_time] = is_xport mflux_time[*,i_time] = mflux mass_slab_next = mass_slab + mflux * time_delta ; update is_volatile is_volatile_next = is_volatile indx_newbare = where(is_volatile and (mass_slab_next lt 0), nindx_newbare) indx_stillvolatile = where(is_volatile and (mass_slab_next ge 0), nindx_stillvolatile) if nindx_newbare gt 0 then begin ; mass_gain is always negative mass_gain = total(mass_slab_next[indx_newbare] * angarea_delta[indx_newbare]) mass_stillvolatile = total(mass_slab_next[indx_stillvolatile] * angarea_delta[indx_stillvolatile]) mass_slab_next[indx_stillvolatile] = mass_slab_next[indx_stillvolatile] * (1 + mass_gain / mass_stillvolatile) mass_slab_next[indx_newbare] = 0. is_volatile_next[indx_newbare] = 0 endif indx_newvolatile = where( (is_volatile eq 0) and (temp_next[*,0] lt temp_volatile_next), nindx_newvolatile) if nindx_newvolatile gt 0 then begin temp_next[indx_newvolatile,*] = temp_volatile_next is_volatile_next[indx_newvolatile] = 1 endif mass_atm_next = 4*!pi*n2vp(temp_volatile_next) / gravacc temp = temp_next mass_slab = mass_slab_next is_volatile = is_volatile_next indx_volatile = where(is_volatile, nindx_volatile, compl=indx_not_volatile, ncompl=nindx_not_volatile) ; update composition %%% could be at end of time-step indxcomp = 2-is_volatile ; is_volatile=1 -> indx=1, is_volatile=0 -> indx=2 albedo = albedo_comp[indxcomp] albedo_geom = albedo_geom_comp[indxcomp] emis = emis_comp[indxcomp] endfor ; ^^^^^^^^^^^ ;print, run, ' COMPLETED NORMALLY' ;----------------------- ; mass slab ;----------------------- mass_atm_time = 4*!pi*n2vp(temp_volatile_time) / gravacc mass_slab_tot_time = total(mass_slab_time * (angarea_delta#onev(n_time)),1) ;window,1 tv = temp_volatile_time[n_time - n_time_per_period:n_time-1] ts = temp_surf_time[*,n_time - n_time_per_period:n_time-1] res = { $ albedo_comp:albedo_comp, emis_comp:emis_comp, $ av:av, ev:ev, as:as, es:es, therminertia:therminertia, mvtot:mvtot, n_off:n_off, $ lat:lat, $ n_lat_edge:n_lat_edge, $ n_time_per_period:n_time_per_period, $ n_time:n_time, $ a : albedo_geom_time[n_time - n_time_per_period:n_time-1], $ tv : tv, ts:ts, $ ms : mass_slab_time[*,n_time - n_time_per_period:n_time-1], $ v: vel_time[*,n_time - n_time_per_period:n_time-1], $ vs : vel_sound_time[*,n_time - n_time_per_period:n_time-1], $ p : n2vp(tv), $ yr : phase_year[n_time - n_time_per_period:n_time-1], $ lat0:lat0_time[n_time - n_time_per_period:n_time-1], $ r:dist_sol_au_time[n_time - n_time_per_period:n_time-1], $ ang:trueanom_time[n_time - n_time_per_period:n_time-1], $ isv:is_volatile_time[*,n_time - n_time_per_period:n_time-1] $ } res_all = { $ albedo_comp:albedo_comp, emis_comp:emis_comp, $ av:av, ev:ev, as:as, es:es, therminertia:therminertia, mvtot:mvtot, n_off:n_off, $ lat:lat, $ n_lat_edge:n_lat_edge, $ n_time_per_period:n_time_per_period, $ n_time:n_time, $ a : albedo_geom_time, $ tv : tv, ts:ts, $ ms : mass_slab_time, $ v: vel_time, $ vs : vel_sound_time, $ p : n2vp(tv), $ yr : phase_year, $ lat0:lat0_time, $ r:dist_sol_au_time, $ ang:trueanom_time, $ isv:is_volatile_time, $ time:time, $ z:z, $ temp_surf_time:temp_surf_time, $ temp_time:temp_time, $ temp_volatile_time:temp_volatile_time, $ eflux_sol_time:eflux_sol_time, $ is_volatile_time:is_volatile_time, $ is_xport_time:is_xport_time, $ mflux_time:mflux_time, $ mass_slab_time:mass_slab_time, $ trueanom_time: trueanom_time, $ dist_sol_au_time:dist_sol_au_time, $ lat0_time:lat0_time, $ albedo_geom_time:albedo_geom_time, $ vel_time:vel_time, $ vel_sound_time:vel_sound_time , $ z_atm_time:z_atm_time, $ z_mfp_time:z_mfp_time $ } save, res, file=run+'.sav' save, res_all, file=run+'_all.sav' return, res end