; ;jeff snider, univeristy of wyoming, march 2006 ; pro main, w, title, r_n_twomey, r_n_model, j_updraft ; common constants, gascon, rho_h2o, mw_nh42so4, rho_nh42so4, sigma_0, $ sigma_t, dsigmadmolality_nh42so4, dsigmadmolality_nacl, mw_h2o, $ mw_nacl, rho_nacl, tkmelt, c_nh42so4, a_nh42so4, c_nacl, a_nacl, $ cp_air_o, cp_h2o_o, mw_air, gravity, p0, alpha, beta_local, $ gascon_h2o, cw_h2o_o, tk_o, lv_o, ew_o, tk_aerosol, salt_type, $ epsilon, c_twomey, k_twomey, aerosol_type, right_tail, $ gascon_air, tc_base, tk_base, p_base, h_start, hmax, $ r0, dt, rmin, rmax, coef, epsilon_aerosol, specific_volume_meas, $ nchan, mw_salt, rho_salt, a_salt, c_salt, specific_volume_base, $ dsigmadmolality_salt, rho_insoluble, volume_insoluble_to_soluble, $ mixrat_tot_1,r_1,geo_sigma_1,mixrat_tot_2,r_2,geo_sigma_2, $ tk_start,p_start,sratio_start,mixrat_tot_3,r_3,geo_sigma_3 ; common graphics, write_flag ; common aerosol, r_m, ncum_mixrat ; print, format = '(2(a,x,e15.4))', 'p_start, Pa = ',p_start,' tk_start = ',tk_start ; icolor = 1 p_parcel_previous = 0.0d ;pressure at previous time step p_parcel = p_start ;initial pressure tk_parcel = tk_start ;initial temperature tc_parcel = tk_start - tkmelt height = 0.0 ;msl height h_start = 0.0 ;initial msl height (initialized below) ssmax = 0.0 ;maximum environmental supersaturation sratio_parcel = sratio_start ;initial saturation ratio sratio_previous = 0. nstep_real = (hmax/w) / dt ;number of steps nstep = long (nstep_real) height_above_cldbase = 200. ;number of meters above thermodynamic cloud base cut_equivalent_radius = 0.425d-6 z_past_cldbase = 0.0 ;height above the cloud base flag_past_cldbase = 0 flag_base = 0 flag_ssmax = 0 write_flag = 0 ;used for positioning of text on graphics r_l_isentrope = 0.d ;isentropic value of liquid water mixing ratio clock = 0. ;time from start of simulation plot_text1 = '' plot_text2 = '' ; ;Matrices de sortie ; r_n = fltarr(nstep) ;number of droplets larger than r0 per kg of dry air r_n_activated = fltarr(nstep) ;number of droplets formally activated per kg of dry air ss = fltarr(nstep) ;Sursaturation h = fltarr(nstep) ;Hauteur time = fltarr(nstep) ;Temps lwc = fltarr(nstep) ;liquid water content rmode = fltarr(nstep) ;mode diameter of droplets dbz = fltarr(nstep) ;reflectivity in dBZ p = fltarr(nstep) ;pressure v = fltarr(nstep) ;specific volume alwc = fltarr(nstep) ;isentropic liquid water content tau_growth = dblarr(nchan+1) ;characteristic time for droplet growth r_eq_now = dblarr(nchan+1) ;equilibrium wet size at current time step r_eq_previous = dblarr(nchan+1) ;equilibrium wet size at previous time step tau_equilibrium = fltarr(nchan+1) ;characteristic time for following the kohler curve growth_flag = intarr(nchan+1) ;flag set when tau_quilibrium is less than tau_growth wgtf_h2o = dblarr(nchan) ;weight fraction of water in droplet rho = dblarr(nchan) ;density of water/salt solution r_bin_ave = dblarr(nchan) ;average size of bin r_crit_ave = dblarr(nchan) ;average critical radius rdry_bin_ave = dblarr(nchan) r1 = fltarr(nchan) r2 = fltarr(nchan) dmixingratiodlogr = fltarr(nchan) growth_fact = fltarr(nchan) ; ;this flag delineates if growth is equilibium or kinetic ; growth_flag(0:nchan) = 0 ; ;use the twomey equation as described in Johnson (JAS, 215-217, 1981) ; w_passed = float(w) tk_passed = float(tk_base) p_passed = float(p_base) twomey, w_passed, tk_passed, p_passed, r_n_twomey ; ;initialization of thermodynamics ; ew = ew_jeff(tk_parcel) e = ew*sratio_parcel r_v = epsilon * e / (p_parcel - e) q_v = r_v / (1.d + r_v) gascon_mix = (1.d - q_v)*gascon_air + q_v*gascon_h2o cp_mix = (1.d - q_v)*cp_air_o + q_v*cp_h2o_o specific_volume = (tk_parcel/p_parcel)*(gascon_air + r_v*gascon_h2o) ;volume occupied by parcel containing 1 kg of dry air specific_volume_initial = specific_volume lv_t = lv_o + (cw_h2o_o - cp_h2o_o)*(tk_o - tk_parcel) ; ;evaluate the initial msl height assuming hydrostatics and adiabatic layer with uniform himidity ; chi_mix = gascon_mix / cp_mix tk0 = tk_parcel * (p0 / p_parcel)^chi_mix dt_dz = gravity / cp_mix height = float ((tk0 - tk_parcel) / dt_dz) h_start = height ; ;aerosol initialization: the procedure ccn_spectrum assumes that aerosol/ccn ;is expressed as a mixing ratio (count per kg of dry air). The ;procedure ccn_spectrum returns the aerosol as a concentration at the initial ;state of the parcel model. Conversion from mixing ratio to concentration is ;performed using the parameter specific_volume_initial ; ccn_spectrum, specific_volume_initial, nconc_initial, rsalt, rinsoluble, rdry, r, $ sr_crit, s_sfc, rho_solution_vector, r_crit ; ;the array nconc_initial is the concentration initially. This is diluted ;because of expansion. Parameters used in the dilution are the current value ;of specific volume (specific_volume) and the ;initial value of the specific volume (specific_volume_initial) ; nconc_parcel = nconc_initial ; ;store the equilibrium size (a.k.a Kohler size) that initilizes the model ; r_eq_now = r ; ;using the Zou and Fukuta growth law (Atmospheric Research, 52, 115-141, 1999) ; alpha_prime = 2.d*alpha / (2.d - alpha) beta_prime = 2.d*beta_local / (2.d - beta_local) ; ;initialize moments of the wet size distribution...bin-averaged sizes ; r_l_insoluble = 0.d r_l = 0.d sixth_moment_radar = 0.d lwc_gm_per_cubic_meter = 0.d mass_dry_ug_per_cubic_meter_meas = 0.d average_density = (1.d / epsilon_aerosol) / (1.d/rho_salt + (1./epsilon_aerosol - 1.d)/rho_insoluble) for i = 0, nchan - 1 do begin r_crit_ave(i) = (r_crit(i) + r_crit(i+1)) / 2.d r_bin_ave(i) = (r(i) + r(i+1)) / 2.d rdry_bin_ave(i) = (rdry(i) + rdry(i+1)) / 2.d wet_size = r_bin_ave(i) salt_size = (rsalt(i) + rsalt(i+1)) / 2.d insoluble_size = (rinsoluble(i) + rinsoluble(i+1)) / 2.d crit_flag = 0 sratio_bdry, crit_flag, wet_size, salt_size, insoluble_size, sratio, loop_count, mw_on_mtot, rho_solution, $ water_activity, molality, wgtf, sigma, kelvin wgtf_h2o(i) = mw_on_mtot rho(i) = rho_solution r_l_insoluble = r_l_insoluble + (4.d*!pi/3.d)*rho_insoluble*insoluble_size^3.d*nconc_parcel(i) * specific_volume r_l_increment = (4.d*!pi/3.d)*specific_volume*nconc_parcel(i)*wgtf_h2o(i)*rho(i)*(r_bin_ave(i)^3.d - insoluble_size^3.d) r_l = r_l_increment + r_l if r_l_increment lt 0.d then begin print,'*parcel* condensate mass problem ',i,r_bin_ave(i),insoluble_size stop endif sixth_moment_radar = nconc_parcel(i)*(r_bin_ave(i)*2.d3)^6.d + sixth_moment_radar lwc_gm_per_cubic_meter = 1.d3*(4.d*!pi/3.d)*nconc_parcel(i)*wgtf_h2o(i)*rho(i)*(r_bin_ave(i)^3.d - insoluble_size^3.d) + lwc_gm_per_cubic_meter if rdry_bin_ave(i) lt cut_equivalent_radius then begin mass_dry_ug_per_cubic_meter_meas = 1.d9*specific_volume_initial*nconc_initial(i)*average_density* $ (4.d*!pi/3.d)*rdry_bin_ave(i)^3.d / (tkmelt*gascon_air / p0) + $ mass_dry_ug_per_cubic_meter_meas endif endfor ; printf, 7, specific_volume_initial printf, 7, average_density printf, 7, mass_dry_ug_per_cubic_meter_meas ; ;water budget initilization ; r_l_previous = r_l ; ;increment in liquid water mixing ratio is (assumed) zero (no latent heat effects initially) ; dr_l = 0.d r_tot = r_v + r_l ;sum of vapor and condensate cp_tot = cp_air_o + r_v*cp_h2o_o + r_l*cw_h2o_o ;heat capacity per unit mass dry air ; ;write and plot the initialization ; initial_output, icolor, rdry, rsalt, rinsoluble, r, nconc_parcel, sixth_moment_radar, $ lwc_gm_per_cubic_meter, sratio_parcel, tk_parcel, p_parcel, $ tc_parcel, ssmax, w, clock, sr_crit, specific_volume, plot_text1, $ plot_text2, title ; ;begin the calculation loop ; for ii = long (0), nstep - 1 do Begin ; if ii gt 0 then begin ; print, format = '(a,i10,f10.1,f10.6)', '*parcel*', ii, h(ii-1), sratio_parcel ; print, format = '(a,i10,9(e11.3))','*parcel*',ii,tk_parcel,p_parcel,ew,r_v,r_l,dr_l,cp_tot,gravity,dt_dz if (sratio_parcel lt sratio_start or sratio_parcel gt 1.10d) then begin print, format = '(10(x,e10.3))', s_sfc(0:nchan) print, format = '(10(x,e10.3))', r(0:nchan) print, format = '(10(x,i10))', growth_flag(0:nchan) print,format='(a,e10.3))','Saturation ratio unrealistic ',sratio_parcel stop endif endif ; ;..evaluate mixing ratio of activated droplets, liquid water volume mixing ratio of activated droplets ;..(ignoring insoluble volume), and the volume medium size. Note: vbar is actually vbar/(4.d*!pi/3.d) ; r1 = r(1:nchan) r2 = r(0:nchan-1) r_bin_ave_test = (r1 + r2) / 2.d index_count = 0 index = where (r_bin_ave_test(*) gt r0, index_count) if index_count gt 0 then begin r_n(ii) = specific_volume * total(nconc_parcel(index)) r_n_model = r_n(ii) vbar = total(nconc_parcel(index)*r_bin_ave_test(index)^3.d) / total(nconc_parcel(index)) rmode(ii) = (vbar)^(1.d/3.d) if sixth_moment_radar gt 0. then begin dbz(ii) = 10.*alog10(sixth_moment_radar) endif else begin dbz(ii) = 0. endelse endif else if index_count le 0 then begin r_n(ii) = 0. r_n_model = 0. rmode(ii) = 0. endif lwc(ii) = r_l / specific_volume ; ;if index_count gt 0 then begin ; print,format = '(3(x,i10),6(x,e10.3),(x,f10.4))', $ ; nchan,index(0),growth_flag(index(0)), $ ; tau_growth(index(0)),tau_equilibrium(index(0)),rsalt(index(0)), $ ; r_bin_ave(index(0)),r0,r_n(ii),(sratio_parcel-1.)*100. ;endif ; ;calculate the concentration of droplets that have exceeded their critical size ; index_count = 0 index = where (r_bin_ave(*) gt r_crit_ave(*), index_count) if index_count gt 0 then begin r_n_activated(ii) = specific_volume*total(nconc_parcel(index)) r_n_activated_now = r_n_activated(ii) endif else if index_count le 0 then begin r_n_activated(ii) = 0. r_n_activated_now = 0. endif ; ;calculate the isentropic water content ; start_satd_isentrope, p_parcel, r_l_isentrope, tk_isentrope, r_v_isentrope specific_volume_isentrope = (tk_isentrope/p_parcel)*(gascon_air + r_v_isentrope*gascon_h2o) alwc(ii) = r_l_isentrope / specific_volume_isentrope ; ;store P,h, ss in arrays ; ss(ii) = (sratio_parcel-1.)*100. p(ii) = p_parcel time(ii) = clock h(ii) = height v(ii) = specific_volume ; ;height increment, hydrostatic and thermodynamic response ; dz = w*dt clock = clock + dt height = height + dz effective_density = p_parcel / (gascon_mix*tk_parcel) + r_l / specific_volume p_parcel = p_parcel - effective_density*gravity*dz tk_parcel = tk_parcel + (-(1.d + r_v + r_l)*gravity*dz + lv_t*dr_l) / cp_tot tc_parcel = tk_parcel - tkmelt diffusivity = dif_h2o(p_parcel,tk_parcel) conductivity = cond_air(tc_parcel) ; ;vertical rate of temperature change used to check the physics, but not in parcel integration ; dt_dz = ((-(1.d + r_v + r_l)*gravity*dz + lv_t*dr_l) / cp_tot) / dz ; ;account for vertical advection of the aerosol ; specific_volume = (tk_parcel/p_parcel)*(gascon_air + r_v*gascon_h2o) nconc_parcel = nconc_initial * specific_volume_initial / specific_volume ; ;account for vertical advection of the water vapor ; e = r_v*p_parcel / (epsilon + r_v) ; ;thermodynamics ; ew = ew_jeff(tk_parcel) sratio_parcel = e / ew lv_t = lv_o + (cw_h2o_o - cp_h2o_o)*(tk_o - tk_parcel) ; for i = 0, nchan do begin ; ;if i eq 19 then begin ; print,format='(i5,4(i5,x,f7.3))', i,growth_flag(i-1),r(i-1)*1.e6,growth_flag(i),r(i)*1.e6,growth_flag(i+1),r(i+1)*1.e6,growth_flag(i+2),r(i+2)*1.e6 ; print,format='(8(e10.2))',tau_equilibrium(i-1),tau_growth(i-1),tau_equilibrium(i),tau_growth(i),tau_equilibrium(i+1),tau_growth(i+1),tau_equilibrium(i+2),tau_growth(i+2) ;endif ; ;evaluate the time constant for kinetically-controlled droplet growth ; if growth_flag(i) eq 0 then begin ; ;calculate growth time for kinetic growth ; growth_coef_cond = rho_solution_vector(i) * (lv_t /(gascon_h2o*tk_parcel) - 1.d) * (lv_t / (conductivity*tk_parcel)) length_cond = conductivity*(2.d*!pi*gascon_air*tk_parcel)^0.5d / (3.d*alpha_prime*p_parcel*gascon_air) f_alpha = r(i) / (r(i) + length_cond) growth_coef_dif = rho_solution_vector(i)*gascon_h2o*tk_parcel / (ew*diffusivity) length_dif = (diffusivity / beta_prime) * (2.d*!pi / (gascon_h2o*tk_parcel))^0.5d f_beta = r(i) / (r(i) + length_dif) dr_dt = abs(sratio_parcel - s_sfc(i)) / (r(i)*(growth_coef_cond/f_alpha + growth_coef_dif/f_beta)) tau_growth(i) = r(i) / dr_dt ; ;calculate growth time for particles following the kohler curve ; if sr_crit(i) lt sratio_parcel then begin growth_flag(i) = 2 ; print, nchan, i, growth_flag(i), tau_growth(i), tau_equilibrium(i), rsalt(i) ; wait, 5 endif else begin r_eq_previous(i) = r_eq_now(i) salt_size = rsalt(i) insoluble_size = rinsoluble(i) kohler_size, sratio_parcel, salt_size, insoluble_size, loop_count, wet_size, mw_on_mtot, rho_solution r_eq_now(i) = wet_size ; if r_eq_now(i) le r_eq_previous(i) then begin ; print,'*parcel* equilibrium evaporation',rsalt(i),r_eq_now(i),r_eq_previous(i) ; endif tau_equilibrium(i) = 0.5 * (r_eq_previous(i) + r_eq_now(i)) * dt / abs ((r_eq_now(i)-r_eq_previous(i))) endelse endif ; ;large droplets and kinetically-limited droplets require a growth law see cooper et al. (JAS, 1449-1469, 1997) ;note that the values of tau_equilibrium(i) are unrealistic near Smax ; if (tau_growth(i) ge tau_equilibrium(i) and tau_equilibrium(i) lt 1.e5 and growth_flag(i) eq 0) $ or (growth_flag(i) gt 0) then begin if growth_flag(i) eq 0 then begin growth_flag(i) = 1 ;print, nchan, i, growth_flag(i), tau_growth(i), tau_equilibrium(i), rsalt(i) endif ddt = (dt / 10.d) integration_time = 0.d while integration_time lt dt do begin integration_time = integration_time + ddt salt_size = rsalt(i) wet_size = r(i) insoluble_size = rinsoluble(i) crit_flag = 0 sratio_bdry, crit_flag, wet_size, salt_size, insoluble_size, sratio, loop_count, mw_on_mtot, rho_solution, $ water_activity, molality, wgtf, sigma, kelvin s_sfc(i) = sratio growth_coef_cond = rho_solution * (lv_t /(gascon_h2o*tk_parcel) - 1.d) * (lv_t / (conductivity*tk_parcel)) length_cond = conductivity * (2.d*!pi*gascon_air*tk_parcel)^0.5d/(3.d*alpha_prime*p_parcel*gascon_air) f_alpha = r(i) / (r(i) + length_cond) growth_coef_dif = rho_solution*gascon_h2o*tk_parcel / (ew*diffusivity) length_dif = (diffusivity / beta_prime) * (2.d*!pi / (gascon_h2o*tk_parcel))^0.5d f_beta = r(i) / (r(i) + length_dif) dr = ddt*(sratio_parcel-s_sfc(i)) / (r(i)*(growth_coef_cond/f_alpha + growth_coef_dif/f_beta)) r(i) = dr + r(i) ;print, sratio, r(i), rdry(i) ; if (i eq nchan-5 and integration_time eq ddt) then begin ; print, format = '(a,i10,11(e11.2))','*parcel*',ii,tk_parcel,p_parcel,r(i),rho_solution,lv_t,conductivity,growth_coef_cond,alpha_prime,length_cond,f_alpha,growth_coef_cond/f_alpha ; print, format = '(a,i10,11(e11.2))','*parcel*',ii,tk_parcel,p_parcel,r(i),rho_solution,lv_t,diffusivity,growth_coef_dif,beta_prime,length_dif,f_beta,growth_coef_dif/f_beta ; endif endwhile if growth_flag(i) ge 1 and i gt 0 then begin if r(i-1) le r(i) then begin print,format='(a,i5,2(i5,x,f7.3))','*****parcel***** growth transition problem ', $ i,growth_flag(i-1),r(i-1)*1.e6,growth_flag(i),r(i)*1.e6 print,format='(a,4(e10.2))','*****parcel***** growth transition problem ', $ tau_equilibrium(i-1),tau_growth(i-1),tau_equilibrium(i),tau_growth(i) stop endif endif endif else begin r(i) = r_eq_now(i) endelse endfor ; ;end of droplet growth section ; ;calculate the liquid water mixing ratio ; r_l = 0.d sixth_moment_radar = 0.d lwc_gm_per_cubic_meter = 0.d for i = 0, nchan - 1 do begin salt_size = (rsalt(i) + rsalt(i+1)) / 2.d r_bin_ave(i) = (r(i) + r(i+1)) / 2.d wet_size = r_bin_ave(i) insoluble_size = (rinsoluble(i) + rinsoluble(i+1)) / 2.d crit_flag = 0 sratio_bdry, crit_flag, wet_size, salt_size, insoluble_size, sratio, loop_count, mw_on_mtot, rho_solution, $ water_activity, molality, wgtf, sigma, kelvin wgtf_h2o(i) = mw_on_mtot rho(i) = rho_solution r_l_increment = (4.d*!pi/3.d)*specific_volume*nconc_parcel(i)*wgtf_h2o(i)*rho(i)*(r_bin_ave(i)^3.d - insoluble_size^3.d) r_l = r_l_increment + r_l if r_l_increment lt 0.d then begin print,format='(a,i5,2e17.5))','*parcel* condensate mass problem ',i,r_bin_ave(i),insoluble_size stop endif sixth_moment_radar = nconc_parcel(i)*(r_bin_ave(i)*2.d3)^6.d + sixth_moment_radar lwc_gm_per_cubic_meter = 1.d3*(4.d*!pi/3.d)*nconc_parcel(i)*wgtf_h2o(i)*rho(i)*(r_bin_ave(i)^3.d - insoluble_size^3.d) + lwc_gm_per_cubic_meter endfor ; ;water mass budget ; dr_l = r_l - r_l_previous r_v = r_tot - r_l r_l_previous = r_l ; ;thermodynamic properties used in next turn of the loop ; q_v = r_v / (1.d + r_v) gascon_mix = (1.d - q_v)*gascon_air + q_v*gascon_h2o cp_tot = cp_air_o + r_v*cp_h2o_o + r_l*cw_h2o_o ; ;capture the spectra at the thermodynamic cloud base and at prescribed height above cloud base ; if ii ge 1 then begin ; if (p_parcel_previous gt p_base and p_parcel lt p_base) or (flag_base gt 0) then begin flag_base = 1 + flag_base z_past_cldbase = z_past_cldbase + dz if z_past_cldbase gt height_above_cldbase then begin ;if lwc_gm_per_cubic_meter gt 0.05 then begin flag_past_cldbase = 1 + flag_past_cldbase endif endif ; if sratio_parcel lt sratio_previous then begin flag_ssmax = flag_ssmax + 1 if flag_ssmax eq 1 then ssmax = (sratio_parcel-1.)*100. endif ; ;write and plot ; if flag_base eq 1 or flag_past_cldbase eq 1 or flag_ssmax eq 1 then begin plot_text1 = +$ 'p(pa)= '+string(p_parcel,format='(e10.3)')+$ ' T(C)= '+string(tc_parcel,format='(f5.1)')+$ ' ss(%)= '+string((sratio_parcel-1.)*100.,format='(f8.4)')+$ ' ssmax(%)= '+string(ssmax,format='(f7.4)')+$ ' w(cm/s)= '+string(w*100.,format='(f6.1)')+$ ' time(s)= '+string(time(ii),format='(f6.1)')+$ ' CDNC(cm-3)= '+string(1.e-6*r_n_model/specific_volume,format='(e10.3)')+$ ' CDNC_act(cm-3)= '+string(1.e-6*r_n_activated_now/specific_volume,format='(e10.3)')+$ ' nconc_twomey(cm-3)= '+string(1.e-6*r_n_twomey/specific_volume,format='(e10.3)')+$ ' alpha= '+string(alpha,format='(f5.3)')+$ ' beta_local= '+string(beta_local,format='(f5.3)')+$ ' LWC(gm/m^3) = '+string(lwc_gm_per_cubic_meter,format='(f10.5)')+$ ' dBZ= '+string(10.*alog10(sixth_moment_radar),format='(f10.1)') plot_text2 = +$ 'epsilon_aerosol= '+string(epsilon_aerosol,format='(f5.3)')+$ ' '+salt_type+$ ' '+aerosol_type+$ ' right_tail = '+string(right_tail,format='(i2)')+$ ' c_twomey(cm-3)= '+string(1.e-6*c_twomey,format='(f5.0)')+$ ' k_twomey= '+string(k_twomey,format='(f5.3)')+$ ' conc_tot_1(cm-3)= '+string(1.e-6*mixrat_tot_1/specific_volume_meas,format='(f6.0)')+$ ' r_1(nm)= '+string(r_1*1.e9,format='(f5.0)')+$ ' geo_sigma_1= '+string(geo_sigma_1,format='(f4.2)')+$ ' conc_tot_2(cm-3)= '+string(1.e-6*mixrat_tot_2/specific_volume_meas,format='(f6.0)')+$ ' r_2(nm)= '+string(r_2*1.e9,format='(f5.0)')+$ ' geo_sigma_2= '+string(geo_sigma_2,format='(f4.2)')+$ ' conc_tot_3(cm-3)= '+string(1.e-6*mixrat_tot_3/specific_volume_meas,format='(f7.2)')+$ ' r_3(nm)= '+string(r_3*1.e9,format='(f5.0)')+$ ' geo_sigma_3= '+string(geo_sigma_3,format='(f4.2)')+$ ' rmin= '+string(rmin,format='(e9.2)')+$ ' rmax= '+string(rmax,format='(e9.2)') title = plot_text1+' '+plot_text2 printf, 4, format='(a)', title r1 = r(1:nchan) r2 = r(0:nchan-1) dmixingratiodlogr = nconc_parcel * specific_volume / (alog10 (r2/r1)) printf, 2, format = '(a)', title printf, 2, format = '(i10)', nchan printf, 2, format = '(10(x,f10.3))', 1.e-6*dmixingratiodlogr(0:nchan-1) printf, 2, format = '(10(x,f10.3))', 1.e6*r(0:nchan-1) printf, 2, format = '(10(x,f10.3))', 1.e6*r(1:nchan) icolor = icolor + 1 plot_droplets, plot_text1, plot_text2, nconc_parcel, r, specific_volume, icolor endif endif ; ;recompute droplet surface conditions, do this for both the kinetically-limited and equilibrium droplets ; for i = 0, nchan do begin wet_size = r(i) salt_size = rsalt(i) insoluble_size = rinsoluble(i) crit_flag = 0 sratio_bdry, crit_flag, wet_size, salt_size, insoluble_size, sratio, loop_count, mw_on_mtot, rho_solution, $ water_activity, molality, wgtf, sigma, kelvin s_sfc(i) = sratio rho_solution_vector(i) = rho_solution endfor ; ;update the value of p_parcel_previous ; p_parcel_previous = p_parcel sratio_previous = sratio_parcel ; ;end of parcel calculation ; endfor ; ;supersaturation versus altitude ; !p.position = [0.15,0.20,0.45,0.50] yrange1 = h_start - 5. yrange2 = h_start + 995. plot, /noerase, ss, h, psym = 0, xticks = 4, xrange = [-6.0, 2.0], yticks = 4, yminor = 5, $ yrange = [yrange1,yrange2], xtitle='Supersaturation, %', ytitle='Height, m', color = 0, $ xminor = 2 ; ;liquid water content, both adiabatic and actual, versus altitude ; !p.position = [0.60,0.20,0.90,0.50] plot, /noerase, lwc*1.e3, h, psym = 0, line = 0, xticks = 5, xrange = [-0.2,2.3], yticks = 4, yminor = 5, $ ytickname = [' ',' ',' ',' ',' '], yrange = [yrange1,yrange2], $ xtitle='Liquid Water Content, g m!E-3!N', color = 0, xminor = 5 oplot, alwc*1.e3, h, line = 1, color = 0 ; ;write P, z, and S to the *.out file ; for ii = long (0), nstep - 1 do begin if ii eq 0 then begin printf, 2, format = '(a)', 'pres(mb) time(s) z(m) ssat(%) lwc(g/m^3) CDNC(cm-3) actCDNC(cm-3) rmode(m) dBz' endif printf, 2, format = '(f10.2,f10.2,f10.2,f10.3,f10.3,f10.2,f10.2,e12.2,e12.2,e12.2)', $ p(ii)/100., time(ii), h(ii), ss(ii), 1.e3*lwc(ii), 1.e-6*r_n(ii)/v(ii), 1.e-6*r_n_activated(ii)/v(ii), rmode(ii), dbz(ii) endfor ; ;clear the window ; erase ; return ; end