;
;jeff snider, univeristy of wyoming, march 2006
;
;calculate the wet size from passed values of the satruation ratio
;temperature, salt type, dry size, and insoluble_size
;use the thermo data in Tang and Munkelwitz (1994) and in Tang (1996)
;
pro kohler_size, sratio, salt_size, insoluble_size, loop_count, wet_size, mw_on_mtot, rho_solution
;
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
;
tk = tk_aerosol
;
;check that passed sratio larger than crystalization relative humidity
;this is the limit of the Tang data
;
if salt_type eq 'nacl' and sratio lt 0.48d then begin
 print, '*kohler_size*  ',salt_type,' requested sratio too small', sratio
 return
endif else if salt_type eq 'nh42so4' and sratio lt 0.40d then begin
 print, '*kohler_size*  ',salt_type,' requested sratio too small', sratio
 return
endif
;
;  there are three unknowns: weight percent composition (wgtp),
;  solution density (rho), and wet size (wet_size).  These
;  variables are related by the activity(wgtp) and rho(wgtp) parameterizations
;  given in Tang (1996) and Tang and Munkelwitz (1994), and by
;  the geometric relationship wgtp(rho,solution_size,salt_size,dry_density)
;  (assuming spheres)
;
;  match ambient rh (sratio) to the product of the solute term
;  and the kelvin term.  See Chylek and Wong (1998) for
;  formulation of the kohler function
;
;  Include the temperature correction for surface tension of pure water
;  note: there is no temperature-dependent correction for activity
;  and the Tang data is only good for 298.15, but, the composition
;  variables are independent of temperature.  Accounting for
;  temperature requires that the temp dependence of rho and
;  activity be parameterized
;
sigma_at_tk = sigma_0 + sigma_t * (tkmelt - tk)
;
;initilization
;
if salt_type eq 'nacl' then begin
 wet_size      = (5.d*salt_size^3.d + insoluble_size^3.d)^(1.d/3.d)
endif else if salt_type eq 'nh42so4' then begin
 wet_size      = (3.d*salt_size^3.d + insoluble_size^3.d)^(1.d/3.d)
endif
;
solution_size = (wet_size^3.d - insoluble_size^3.d)^(1.d/3.d)
rho = rho_h2o
wgtp = (100.d*rho_salt*salt_size^3.d/(rho*solution_size^3.d))
funct = 1.d
loop_count = 0
;
while (abs(funct) gt 0.0000001d*sratio) and (loop_count lt 80) do begin
 rho = 1000.d*(0.9971d + a_salt(0)*wgtp + a_salt(1)*wgtp^2.d + a_salt(2)*wgtp^3.d + a_salt(3)*wgtp^4.d)
 solution_size = (100.d*rho_salt*salt_size^3.d/(rho*wgtp))^(1.d/3.d)
 wet_size      = (solution_size^3.d + insoluble_size^3.d)^(1.d/3.d)
 loop_count = loop_count + 1
 wgtf = wgtp / 100.d
 dwgtfdwgtp = 1.d / 100.d
 molality = wgtf / (mw_salt*(1.d - wgtf))
 sigma = sigma_at_tk + dsigmadmolality_salt * molality
 kelvin = exp (2.d*sigma*mw_h2o/(rho*gascon*tk*wet_size))
 dkelvindsigma = kelvin*(2.d*mw_h2o/(rho*gascon*tk*wet_size))
 dkelvindrho = -kelvin*(2.d*sigma*mw_h2o/(gascon*tk*wet_size))*(1.d/rho^2.d)
 dkelvindwet_size = -kelvin*(2.d*sigma*mw_h2o/(rho*gascon*tk))*(1.d/wet_size^2.d)
 dmolalitydwgtf = ((1.d/mw_salt))/(1.d - wgtf)^2.d
 dwet_sizedwgtp = -(wet_size^3.d - insoluble_size^3.d) / (3.d*wet_size^2.d*wgtp)
 drhodwgtp = 1000.d*(a_salt(0) + 2.d*a_salt(1)*wgtp + 3.d*a_salt(2)*wgtp^2.d + 4.*a_salt(3)*wgtp^3.d)
 water_activity = 1.d + c_salt(0)*wgtp + c_salt(1)*wgtp^2.d + c_salt(2)*wgtp^3.d + c_salt(3)*wgtp^4.d
 funct = water_activity * kelvin - sratio
 term1 = (c_salt(0) + 2.d*c_salt(1)*wgtp + 3.d*c_salt(2)*wgtp^2.d + 4.d*c_salt(3)*wgtp^3.d) * kelvin
 term2 = water_activity * dkelvindsigma * dsigmadmolality_salt * dmolalitydwgtf * dwgtfdwgtp
 term3 = water_activity * dkelvindwet_size * dwet_sizedwgtp
 term4 = water_activity * dkelvindrho * drhodwgtp
 dfunctdwgtp = term1 + term2 + term3 + term4
 wgtp_new = wgtp - funct / dfunctdwgtp
;
 if (wgtp_new gt 100.d) then begin
  print, '*kohler_size* wgtp>100%',wgtp,wgtp_new,loop_count,salt_size,wet_size
  wgtp_new = wgtp * 1.2d
 endif else if (wgtp_new lt 0.d) then begin
  print, '*kohler_size* wgtp<  0% (may not be a solution!)',wgtp,wgtp_new,loop_count,salt_size,wet_size
  wgtp_new = wgtp * 0.8d
 endif
;
 wgtp  = wgtp_new
;
endwhile
;
rho = 1000.d*(0.9971d + a_salt(0)*wgtp + a_salt(1)*wgtp^2.d + a_salt(2)*wgtp^3.d + a_salt(3)*wgtp^4.d)
solution_size = (100.d*rho_salt*salt_size^3.d/(rho*wgtp))^(1.d/3.d)
wet_size      = (solution_size^3.d + insoluble_size^3.d)^(1.d/3.d)
wgtf = wgtp / 100.d
molality = wgtf / (mw_salt*(1.d - wgtf))
sigma = sigma_at_tk + dsigmadmolality_salt * molality
kelvin = exp (2.d*sigma*mw_h2o/(rho*gascon*tk*wet_size))
water_activity = 1.0d + c_salt(0)*wgtp + c_salt(1)*wgtp^2.d + c_salt(2)*wgtp^3.d + c_salt(3)*wgtp^4.d
sratio_test = water_activity * kelvin
mw_on_mtot = 1.d - wgtf
rho_solution = rho
;
if sratio_test gt 1.1d then begin
 print, '*kohler_size* ',sratio_test,' requested sratio too large', salt_size, loop_count, wgtp, rho_solution
 stop
endif
;
if wet_size lt 0.d then begin
 print, '*kohler_size* ',wet_size,' wet_size less than zero ',salt_size, loop_count, wgtp, rho_solution
 stop
endif
;
if wet_size gt 60.d-6 then begin
 print, '*kohler_size* ',wet_size,' wet_size greater than 60 um ',salt_size, loop_count, wgtp, rho_solution, sratio, sratio_test
 stop
endif
;
return
;
end
;