;
;jeff snider, univeristy of wyoming, march 2006
;
pro init
;
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  aerosol, r_m, ncum_mixrat
;
my_path = 'c:\jeff\parcel_model_2006_1\'
;
for j_updraft = 0, 0 do begin
;
 set_common, j_updraft, w_passed
;
 filename   = my_path+aerosol_type+'_'+salt_type
 filename2  = filename+'_'+string(j_updraft,format='(i1)')+'.out'
 filename4  = filename+'_w.all'
 filename5  = filename+'_'+string(j_updraft,format='(i1)')+'.ccn'
 filename7  = filename+'_mass.all'
 filenameps = filename+'_'+string(j_updraft,format='(i1)')+'.ps'
;
 r_m         = dblarr(1)
 ncum_mixrat = dblarr(1)
;
;..aerosol initialization has largest particles at index = 0, smallest at index = nchan
;..the sizing set up here is used for both "twomey" and "lognorm"
;
 r_m(0)      = rmax
 test        = 1
;
 while r_m(n_elements(r_m)-1) ge rmin do begin
  if not test then r_m = [r_m, r_m(n_elements(r_m)-1) / coef]
  x_pass_1             = (alog(r_m(n_elements(r_m)-1))-alog(r_1)) / alog(geo_sigma_1)
  gauss_integral_1     = 1.0d - gaussint(x_pass_1)
  x_pass_2             = (alog(r_m(n_elements(r_m)-1))-alog(r_2)) / alog(geo_sigma_2)
  gauss_integral_2     = 1.0d - gaussint(x_pass_2)
  x_pass_3             = (alog(r_m(n_elements(r_m)-1))-alog(r_3)) / alog(geo_sigma_3)
  gauss_integral_3     = 1.0d - gaussint(x_pass_3)
  if not test then begin
   ncum_mixrat         = [ncum_mixrat, gauss_integral_1*mixrat_tot_1 + gauss_integral_2*mixrat_tot_2 + gauss_integral_3*mixrat_tot_3]
  endif else begin
   ncum_mixrat         = gauss_integral_1*mixrat_tot_1 + gauss_integral_2*mixrat_tot_2 + gauss_integral_3*mixrat_tot_3
   test                = 0
  endelse
 endwhile
;
;..number of channels is one less than the number of particle sizes
;
 nchan   = n_elements(r_m) - 1
;
 title            = ''
 r_n_twomey       = -999.
 r_n_model        = -999.
 r_n_twomey_local = -999.
 r_n_model_local  = -999.
;
 set_plot, 'ps'
 device, filename = filenameps, xoffset=1.5, yoffset=1.0, xsize=18.0, ysize=27.0,$
  /Portrait, /color, /encapsul
 tvlct, [   0, 255,   0,   0, 255], $
        [   0,   0, 255,   0,   0], $
        [   0,   0,   0, 250, 255]
 !p.multi = [0,2,2,0,0]
 !x.style = 1
 !y.style = 1
;
 openw, 2, filename2
 openw, 4, filename4, append = 1
 openw, 5, filename5
 openw, 7, filename7
 title = 'p_start= '+string(p_start,format='(e10.3)')+$
        ' tk_start= '+string(tk_start,format='(f6.2)')+$
        ' sratio_start= '+string(sratio_start,format='(f5.3)')+$
        ' p_base= '+string(p_base,format='(e10.3)')+$
        ' tk_base= '+string(tk_base,format='(f6.2)')
 printf, 2, format='(a)', title
 printf, 4, format='(a)', title
 printf, 5, format='(a)', title
 printf, 7, format='(a)', title
 for i = 2, 7 do begin
  if            i eq 2 then begin
   printf, 2, format='(a)', filename2
  endif else if i eq 4 then begin
   printf, 4, format='(a)', filename4
  endif else if i eq 5 then begin
   printf, 5, format='(a)', filename5
  endif else if i eq 7 then begin
   printf, 7, format='(a)', filename7
  endif
 endfor
;
 main, w_passed, title, r_n_twomey, r_n_model, j_updraft
 r_n_twomey_local = r_n_twomey
 r_n_model_local  = r_n_model
;
 device, /close
 close, 2
 close, 4
 close, 5
 close, 7
;
endfor
;
end