;
pro plot_aitken_accum_large
;
;This procedure checks that size spectra used to initialize the parcel model are
;consistent with data reported to the Joint Research Centre of the European Commission by
;Erik Swietlicki of the University of Lund.
;
;The files pdh_*_excel_spectra.txt contain DMA
;size spectra data sampled at the PDH surface site and averaged over the Merlin flight interval.
;
;The files aps_*_excel_spectra.txt contain APS data sampled at the PDH surface site and averaged
;over the Merlin flight interval.
;
;The DMA and APS spectral densities are corrected to count mixing ratios at the
;standard P and T of the PDH site.
;
;The APS data was corrected (one step back) to account for the aerodynamic to sphere-equivalent
;sizing as discussed in Section 2.1.1 of part 2.
;
;The header of DMA and APS files are the fitting parameters corresponding to the bimode fit
;discussed in Section 2.1.1 of part 2.
;
;The third file (*_punto_nh42so4_co1__0_0.dry) contains the size spectrum used to initialize
;the parcel model.   This is based on the discretized DMA size spectra (sizes smaller than
;0.5 um) and on the bimode fit (sizes larger than 0.5 um).   Spectra are extended to sizes larger than
;4 um (see section 3.2.1 of part 2) so that comparisons can be made over the whole range of the data.
;
;APS concentrations at sizes larger than 10 um are zero (in the Swietlicki data).
;
;title4 = ' 9 July'
title4 = '26 June'
;
lu1      = 1
path     = ''
filename1 = 'pdh_fl21_excel_spectra.txt'
complete_path = path+filename1
print, complete_path
openr, lu1, complete_path
nchan_fine = 51
d_um            = fltarr(nchan_fine)
dndlog10d       = fltarr(nchan_fine)
dndlog10d_sigma = fltarr(nchan_fine)
dndlog10d_max   = fltarr(nchan_fine)
dndlog10d_min   = fltarr(nchan_fine)
;
readf, lu1, dbar1, sigmag1, ntot1, mass1
dbar1 = dbar1 / 1000.
readf, lu1, d_um
d_um = d_um / 1000.
readf, lu1, dndlog10d
dndlog10d = (8.314*273.15/(0.029*1.013e5))*dndlog10d
readf, lu1, dndlog10d_sigma
dndlog10d_sigma = (8.314*273.15/(0.029*1.013e5))*dndlog10d_sigma
readf, lu1, dndlog10d_max
dndlog10d_max = (8.314*273.15/(0.029*1.013e5))*dndlog10d_max
readf, lu1, dndlog10d_min
dndlog10d_min = (8.314*273.15/(0.029*1.013e5))*dndlog10d_min
close, lu1
;
lu2      = 2
path     = ''
filename2 = 'aps_fl21_excel_spectra.txt'
complete_path = path+filename2
print, complete_path
openr, lu2, complete_path
nchan_coarse = 46
d_um_aps            = fltarr(nchan_coarse)
dndlog10d_aps       = fltarr(nchan_coarse)
dndlog10d_sigma_aps = fltarr(nchan_coarse)
dndlog10d_max_aps   = fltarr(nchan_coarse)
dndlog10d_min_aps   = fltarr(nchan_coarse)
;
readf, lu2, dbar2, sigmag2, ntot2, mass2
readf, lu2, d_um_aps
readf, lu2, dndlog10d_aps
dndlog10d_aps = (8.314*273.15/(0.029*1.013e5))*dndlog10d_aps
readf, lu2, dndlog10d_sigma_aps
dndlog10d_sigma_aps = (8.314*273.15/(0.029*1.013e5))*dndlog10d_sigma_aps
readf, lu2, dndlog10d_max_aps
dndlog10d_max_aps = (8.314*273.15/(0.029*1.013e5))*dndlog10d_max_aps
readf, lu2, dndlog10d_min_aps
dndlog10d_min_aps = (8.314*273.15/(0.029*1.013e5))*dndlog10d_min_aps
close, lu2
;
lu3      = 3
path     = ''
filename3 = 'fl21_punto_nh42so4_co1__0_0.dry'
complete_path = path+filename3
print, complete_path
openr, lu3, complete_path
dummy = ''
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, dummy
readf, lu3, nchan_model_output
print, nchan_model_output
dndlogr_model_output = fltarr(nchan_model_output)
r_um_low_model_output = fltarr(nchan_model_output)
r_um_hig_model_output = fltarr(nchan_model_output)
readf, lu3, dndlogr_model_output
readf, lu3, r_um_low_model_output
readf, lu3, r_um_hig_model_output
d_um_ave_model_output = 2. * ((r_um_low_model_output + r_um_hig_model_output) / 2.)
print, 'output - max channel diameter fine (um)', d_um_ave_model_output(nchan_model_output-1)
;
readf, lu3, nchan_model_output_coarse
dndlogr_model_output_coarse = fltarr(nchan_model_output_coarse)
r_um_low_model_output_coarse = fltarr(nchan_model_output_coarse)
r_um_hig_model_output_coarse = fltarr(nchan_model_output_coarse)
readf, lu3, dndlogr_model_output_coarse
readf, lu3, r_um_low_model_output_coarse
readf, lu3, r_um_hig_model_output_coarse
d_um_ave_model_output_coarse = 2. * ((r_um_low_model_output_coarse + r_um_hig_model_output_coarse) / 2.)
print, 'output - max channel diameter coarse (um)', d_um_ave_model_output_coarse(nchan_model_output_coarse-1)
close, lu3
;
title1 = filename1+ $
         '  db= '+string(dbar1,format='(f4.2)')+ $
         '  sg= '+string(sigmag1,format='(f4.2)')+ $
         '  nt= '+string(ntot1,format='(f5.0)')+ $
         '  ma= '+string(mass1,format='(f3.0)')
title2 = filename2+ $
         '  db= '+string(dbar2,format='(f4.2)')+ $
         '  sg= '+string(sigmag2,format='(f4.2)')+ $
         '  nt= '+string(ntot2,format='(f5.0)')+ $
         '  ma= '+string(mass2,format='(f3.0)')
title3 = filename3
;
postscript_filename = 'initial_spectra_fl21_.ps'
set_plot, 'ps'
!p.font    =    0
my_xsize_inch   =  8.5
my_ysize_inch   = 11.0
my_yoffset_inch =  0.0
my_xoffset_inch =  0.5
journal_width_inch = 2.24
journal_width_normal = journal_width_inch / my_xsize_inch
journal_height_inch = 2.24
journal_height_normal = journal_height_inch / my_ysize_inch
x_coord_min_normal   = 0.10
y_coord_min_normal   = 0.10
x_coord_max_normal   = x_coord_min_normal + journal_width_normal
y_coord_max_normal   = y_coord_min_normal + journal_height_normal
device, /times, font_size = 9, filename = postscript_filename, inches = 1, yoffset = my_yoffset_inch, $
xoffset = my_xoffset_inch, xsize = my_xsize_inch, ysize = my_ysize_inch
!p.position   = [x_coord_min_normal,y_coord_min_normal,x_coord_max_normal,y_coord_max_normal]
;
ymin        = 0.01
index_count = 0
yplot = dndlog10d
index = where (yplot lt ymin, index_count)
if index_count gt 0 then yplot(index) = ymin / 10.
;
plot, d_um, yplot, /nodata, /xlog, /ylog, xrange = [0.01,10.], $
 ytitle = 'dN/dlog!I10!ND (mg!E-1!N) or dV/dlog!I10!ND (!9m!7m!E3!N mg!E-1!N)', $
 yrange = [ymin, 1000.], xstyle = 1, ystyle = 1, xtitle = 'Dry Diameter, !9m!7m'
;
x = fltarr(2)
y = fltarr(2)
;
xarray = fltarr(nchan_fine)
xarray[*] = 0.
yarray = fltarr(nchan_fine)
j_count = 0
for j = 0, nchan_fine-1, 2 do begin
 x(0) = d_um(j)
 y(0) = dndlog10d(j) + dndlog10d_sigma(j)
 x(1) = d_um(j)
 y(1) = dndlog10d(j) - dndlog10d_sigma(j)
 index = where (y lt ymin, index_count)
 if index_count gt 0 then  y(index) = ymin / 10.
 oplot, x, y
 xarray(j_count) = d_um(j)
 yarray(j_count) = yplot(j)
 j_count = j_count + 1
endfor
;
index = where (xarray gt 0., index_count)
if index_count gt 0 then oplot, xarray(index), yarray(index), psym=5, symsize=0.5
;
yplot = dndlog10d_aps
index = where (yplot lt ymin, index_count)
if index_count gt 0 then  yplot(index) = ymin / 10.
xarray = fltarr(nchan_coarse)
xarray[*] = 0.
yarray = fltarr(nchan_coarse)
j_count = 0
for j = 0, nchan_coarse-1, 2 do begin
 x(0) = d_um_aps(j)
 y(0) = dndlog10d_aps(j) + dndlog10d_sigma_aps(j)
 x(1) = d_um_aps(j)
 y(1) = dndlog10d_aps(j) - dndlog10d_sigma_aps(j)
 index = where (y lt ymin, index_count)
 if index_count gt 0 then  y(index) = ymin / 10.
 oplot, x, y
 xarray(j_count) = d_um_aps(j)
 yarray(j_count) = yplot(j)
 j_count = j_count + 1
endfor
;
index = where (xarray gt 0., index_count)
if index_count gt 0 then oplot, xarray(index), yarray(index), psym=6, symsize=0.5
;
dn_dlogd = fltarr(2*nchan_coarse+2)
dv_dlogd = fltarr(2*nchan_coarse+2)
dplot    = fltarr(2*nchan_coarse+2)
c1_1 = (8.314*273.15/(0.029*1.013e5)) * ntot1 / (sqrt(2.*!pi)*alog10(sigmag1))
c2_1 = 2. * alog(sigmag1) * alog(sigmag1)
c1_2 = (8.314*273.15/(0.029*1.013e5)) * ntot2 / (sqrt(2.*!pi)*alog10(sigmag2))
c2_2 = 2. * alog(sigmag2) * alog(sigmag2)
i = 0
dtest = 0.
while dtest lt dbar2*sigmag2^4 and i lt (2*nchan_coarse+2) - 1 do begin
 if i eq 0 then begin
  dplot(i) = dbar1 / sigmag1^4
 endif else begin
  dplot(i) = dplot(i-1) * 1.10
 endelse
 dn_dlogd(i) = c1_1 * exp ( -alog(dplot(i)/dbar1) * alog(dplot(i)/dbar1) / c2_1) + $
               c1_2 * exp ( -alog(dplot(i)/dbar2) * alog(dplot(i)/dbar2) / c2_2)
 dv_dlogd(i) = dn_dlogd(i) * dplot(i)^3 * !pi / 6.
 dtest = dplot(i)
 i = i + 1
endwhile
dma_upper_limit = 0.200
index_count = 0
index = where (dplot gt dma_upper_limit, index_count)
if index_count gt 0 then begin
 dplot = dplot(index)
 dn_dlogd = dn_dlogd(index)
 oplot, dplot, dn_dlogd
 dv_dlogd = dv_dlogd(index)
 oplot, dplot, dv_dlogd, linestyle = 1
endif
;
vi_upper_limit = 0.85
x(0) = vi_upper_limit
x(1) = vi_upper_limit
y(0) = ymin
y(1) = 10000.
oplot, x, y, linestyle = 2
;
xyouts, 0.12, 0.84, title1, size=0.6, /normal
xyouts, 0.12, 0.83, title2, size=0.6, /normal
xyouts, 0.12, 0.82, title3, size=0.6, /normal
xyouts, x_coord_max_normal - 0.01, y_coord_max_normal - 0.02, title4, size = 1.0, alignment = 1, /normal
;
oplot, d_um_ave_model_output, dndlogr_model_output, linestyle = 3
oplot, d_um_ave_model_output_coarse, dndlogr_model_output_coarse, linestyle = 3
;
device, /close
resetps2x
;
end