37 #ifdef PMC_USE_SUNDIALS 79 real(kind=dp) :: v_comp
83 real(kind=dp) :: d_dry
85 real(kind=dp) :: kappa
94 real(kind=dp) :: hdot_i
96 real(kind=dp) :: hdot_env
98 real(kind=dp) :: dddot_dd
100 real(kind=dp) :: dddot_dh
102 real(kind=dp) :: dhdoti_dd
104 real(kind=dp) :: dhdoti_dh
106 real(kind=dp) :: dhdotenv_dd
108 real(kind=dp) :: dhdotenv_dh
148 env_state_final, del_t)
159 type(
env_state_t),
intent(inout) :: env_state_final
161 real(kind=dp),
intent(in) :: del_t
163 integer :: i_part, n_eqn, i_eqn
165 real(kind=dp) :: init_time, final_time
167 real(kind=dp) :: num_conc
169 real(kind=dp) :: water_vol_conc_initial, water_vol_conc_final
170 real(kind=dp) :: vapor_vol_conc_initial, vapor_vol_conc_final
171 real(kind=dp) :: d_water_vol_conc, d_vapor_vol_conc
172 real(kind=dp) :: V_comp_ratio, water_rel_error
173 #ifdef PMC_USE_SUNDIALS 176 type(c_ptr) :: state_f_p, abstol_f_p
177 integer(kind=c_int) :: n_eqn_f, solver_stat
178 real(kind=c_double) :: reltol_f, t_initial_f, t_final_f
181 #ifdef PMC_USE_SUNDIALS 182 #ifndef DOXYGEN_SKIP_DOC 185 reltol_f, t_initial_f, t_final_f) bind(c)
187 integer(kind=c_int),
value :: neq
188 type(c_ptr),
value :: x_f
189 type(c_ptr),
value :: abstol_f
190 real(kind=c_double),
value :: reltol_f
191 real(kind=c_double),
value :: t_initial_f
192 real(kind=c_double),
value :: t_final_f
199 water_vol_conc_initial = 0d0
201 num_conc = aero_weight_array_num_conc(aero_state%awa, &
202 aero_state%apa%particle(i_part), aero_data)
203 water_vol_conc_initial = water_vol_conc_initial &
204 + aero_state%apa%particle(i_part)%vol(aero_data%i_water) * num_conc
211 = (env_state_final%temp - env_state_initial%temp) / del_t
213 = (env_state_final%pressure - env_state_initial%pressure) / del_t
230 = aero_weight_array_num_conc(aero_state%awa, &
231 aero_state%apa%particle(i_part), aero_data)
233 aero_state%apa%particle(i_part), aero_data)
234 abs_tol_vector(i_part) = max(1d-30, &
240 #ifdef PMC_USE_SUNDIALS 243 n_eqn_f = int(n_eqn, kind=c_int)
244 reltol_f =
real(1d-8, kind=c_double)
245 t_initial_f =
real(0, kind=c_double)
246 t_final_f =
real(del_t, kind=c_double)
248 state_f(i_eqn) =
real(state(i_eqn), kind=c_double)
249 abstol_f(i_eqn) =
real(abs_tol_vector(i_eqn), kind=c_double)
251 state_f_p = c_loc(state_f)
252 abstol_f_p = c_loc(abstol_f)
255 solver_stat =
condense_solver(n_eqn_f, state_f_p, abstol_f_p, reltol_f, &
256 t_initial_f, t_final_f)
257 call condense_check_solve(solver_stat)
263 state(i_eqn) =
real(state_f(i_eqn), kind=
dp)
273 water_vol_conc_final = 0d0
277 num_conc = aero_weight_array_num_conc(aero_state%awa, &
278 aero_state%apa%particle(i_part), aero_data)
281 aero_state%apa%particle(i_part)%vol(aero_data%i_water) &
287 aero_state%apa%particle(i_part)%vol(aero_data%i_water) = max(0d0, &
288 aero_state%apa%particle(i_part)%vol(aero_data%i_water))
291 water_vol_conc_final = water_vol_conc_final &
292 + aero_state%apa%particle(i_part)%vol(aero_data%i_water) * num_conc
301 v_comp_ratio = env_state_final%temp * env_state_initial%pressure &
302 / (env_state_initial%temp * env_state_final%pressure)
303 vapor_vol_conc_initial = aero_data%molec_weight(aero_data%i_water) &
304 / (
const%univ_gas_const * env_state_initial%temp) &
306 * env_state_initial%rel_humid &
308 vapor_vol_conc_final = aero_data%molec_weight(aero_data%i_water) &
309 / (
const%univ_gas_const * env_state_final%temp) &
311 * env_state_final%rel_humid &
313 d_vapor_vol_conc = vapor_vol_conc_final - vapor_vol_conc_initial
314 d_water_vol_conc = water_vol_conc_final - water_vol_conc_initial
315 water_rel_error = (d_vapor_vol_conc + d_water_vol_conc) &
316 / (vapor_vol_conc_final + water_vol_conc_final)
318 "condensation water imbalance too high: " &
329 #ifdef PMC_USE_SUNDIALS 331 subroutine condense_check_solve(value)
334 integer(kind=c_int),
intent(in) :: value
339 call die_msg(123749472,
"condense_solver: " &
340 //
"failed to allocate y vector")
342 call die_msg(563665949,
"condense_solver: " &
343 //
"failed to allocate abstol vector")
345 call die_msg(700541443,
"condense_solver: " &
346 //
"failed to create the solver")
348 call die_msg(297559183,
"condense_solver: " &
349 //
"failure to initialize the solver")
351 call die_msg(848342417,
"condense_solver: " &
352 //
"failed to set tolerances")
354 call die_msg(275591501,
"condense_solver: " &
355 //
"failed to set maximum steps")
357 call die_msg(862254233,
"condense_solver: solver failed")
359 call die_msg(635697577,
"condense_solver: unknown return code: " &
363 end subroutine condense_check_solve
378 real(kind=dp) :: rho_w, M_w, P_0, dP0_dT_div_P0, rho_air, k_a, D_v, U
379 real(kind=dp) :: V, W, X, Y, Z, k_ap, dkap_dD, D_vp, dDvp_dD
380 real(kind=dp) :: a_w, daw_dD, delta_star, h, dh_ddelta, dh_dD
381 real(kind=dp) :: dh_dH, ddeltastar_dD, ddeltastar_dH
382 integer :: newton_step
384 rho_w =
const%water_density
385 m_w =
const%water_molec_weight
386 p_0 =
const%water_eq_vap_press &
387 * 10d0**(7.45d0 * (inputs%T -
const%water_freeze_temp) &
389 dp0_dt_div_p0 = 7.45d0 * log(10d0) * (
const%water_freeze_temp - 38d0) &
390 / (inputs%T - 38d0)**2
391 rho_air =
const%air_molec_weight * inputs%p &
392 / (
const%univ_gas_const * inputs%T)
394 k_a = 1d-3 * (4.39d0 + 0.071d0 * inputs%T)
395 d_v = 0.211d-4 / (inputs%p /
const%air_std_press) &
396 * (inputs%T / 273d0)**1.94d0
397 u =
const%water_latent_heat * rho_w / (4d0 * inputs%T)
398 v = 4d0 * m_w * p_0 / (rho_w *
const%univ_gas_const * inputs%T)
399 w =
const%water_latent_heat * m_w / (
const%univ_gas_const * inputs%T)
400 x = 4d0 * m_w *
const%water_surf_eng &
401 / (
const%univ_gas_const * inputs%T * rho_w)
402 y = 2d0 * k_a / (
const%accom_coeff * rho_air &
403 *
const%air_spec_heat) &
404 * sqrt(2d0 *
const%pi *
const%air_molec_weight &
405 / (
const%univ_gas_const * inputs%T))
406 z = 2d0 * d_v /
const%accom_coeff * sqrt(2d0 *
const%pi * m_w &
407 / (
const%univ_gas_const * inputs%T))
409 outputs%Hdot_env = - dp0_dt_div_p0 * inputs%Tdot * inputs%H &
410 + inputs%H * inputs%pdot / inputs%p
411 outputs%dHdotenv_dD = 0d0
412 outputs%dHdotenv_dH = - dp0_dt_div_p0 * inputs%Tdot &
413 + inputs%pdot / inputs%p
415 if (inputs%D <= inputs%D_dry)
then 416 k_ap = k_a / (1d0 + y / inputs%D_dry)
418 d_vp = d_v / (1d0 + z / inputs%D_dry)
423 delta_star = u * v * d_vp * inputs%H / k_ap
425 outputs%Ddot = k_ap * delta_star / (u * inputs%D_dry)
426 outputs%Hdot_i = - 2d0 *
const%pi / (v * inputs%V_comp) &
427 * inputs%D_dry**2 * outputs%Ddot
431 dh_dh = - u * v * d_vp
433 ddeltastar_dd = - dh_dd / dh_ddelta
434 ddeltastar_dh = - dh_dh / dh_ddelta
436 outputs%dDdot_dD = 0d0
437 outputs%dDdot_dH = k_ap / (u * inputs%D_dry) * ddeltastar_dh
438 outputs%dHdoti_dD = - 2d0 *
const%pi / (v * inputs%V_comp) &
439 * inputs%D_dry**2 * outputs%dDdot_dD
440 outputs%dHdoti_dH = - 2d0 *
const%pi / (v * inputs%V_comp) &
441 * inputs%D_dry**2 * outputs%dDdot_dH
446 k_ap = k_a / (1d0 + y / inputs%D)
447 dkap_dd = k_a * y / (inputs%D + y)**2
448 d_vp = d_v / (1d0 + z / inputs%D)
449 ddvp_dd = d_v * z / (inputs%D + z)**2
450 a_w = (inputs%D**3 - inputs%D_dry**3) &
451 / (inputs%D**3 + (inputs%kappa - 1d0) * inputs%D_dry**3)
452 daw_dd = 3d0 * inputs%D**2 * inputs%kappa * inputs%D_dry**3 &
453 / (inputs%D**3 + (inputs%kappa - 1d0) * inputs%D_dry**3)**2
461 delta_star = delta_star - h / dh_ddelta
462 h = k_ap * delta_star - u * v * d_vp &
463 * (inputs%H - a_w / (1d0 + delta_star) &
464 * exp(w * delta_star / (1d0 + delta_star) &
465 + (x / inputs%D) / (1d0 + delta_star)))
467 k_ap - u * v * d_vp * a_w / (1d0 + delta_star)**2 &
468 * (1d0 - w / (1d0 + delta_star) &
469 + (x / inputs%D) / (1d0 + delta_star)) &
470 * exp(w * delta_star / (1d0 + delta_star) &
471 + (x / inputs%D) / (1d0 + delta_star))
474 abs(h) < 1d3 * epsilon(1d0) * abs(u * v * d_vp * inputs%H), &
475 "condensation newton loop did not satisfy convergence tolerance")
477 outputs%Ddot = k_ap * delta_star / (u * inputs%D)
478 outputs%Hdot_i = - 2d0 *
const%pi / (v * inputs%V_comp) &
479 * inputs%D**2 * outputs%Ddot
481 dh_dd = dkap_dd * delta_star &
482 - u * v * ddvp_dd * inputs%H + u * v &
483 * (a_w * ddvp_dd + d_vp * daw_dd &
484 - d_vp * a_w * (x / inputs%D**2) / (1d0 + delta_star)) &
485 * (1d0 / (1d0 + delta_star)) &
486 * exp((w * delta_star) / (1d0 + delta_star) &
487 + (x / inputs%D) / (1d0 + delta_star))
488 dh_dh = - u * v * d_vp
490 ddeltastar_dd = - dh_dd / dh_ddelta
491 ddeltastar_dh = - dh_dh / dh_ddelta
493 outputs%dDdot_dD = dkap_dd * delta_star / (u * inputs%D) &
494 + k_ap * ddeltastar_dd / (u * inputs%D) &
495 - k_ap * delta_star / (u * inputs%D**2)
496 outputs%dDdot_dH = k_ap / (u * inputs%D) * ddeltastar_dh
497 outputs%dHdoti_dD = - 2d0 *
const%pi / (v * inputs%V_comp) &
498 * (2d0 * inputs%D * outputs%Ddot + inputs%D**2 * outputs%dDdot_dD)
499 outputs%dHdoti_dH = - 2d0 *
const%pi / (v * inputs%V_comp) &
500 * inputs%D**2 * outputs%dDdot_dH
506 #ifdef PMC_USE_SUNDIALS 509 subroutine condense_vf_f(n_eqn, time, state_p, state_dot_p) bind(c)
512 integer(kind=c_int),
value,
intent(in) :: n_eqn
514 real(kind=c_double),
value,
intent(in) :: time
516 type(c_ptr),
value,
intent(in) :: state_p
518 type(c_ptr),
value,
intent(in) :: state_dot_p
520 real(kind=c_double),
pointer :: state(:)
521 real(kind=c_double),
pointer :: state_dot(:)
522 real(kind=dp) :: Hdot
529 call c_f_pointer(state_p, state, (/ n_eqn /))
530 call c_f_pointer(state_dot_p, state_dot, (/ n_eqn /))
538 inputs%H = state(n_eqn)
541 do i_part = 1,(n_eqn - 1)
542 inputs%D = state(i_part)
544 inputs%V_comp = (inputs%T &
550 state_dot(i_part) = outputs%Ddot
551 hdot = hdot + outputs%Hdot_i
553 hdot = hdot + outputs%Hdot_env
555 state_dot(n_eqn) = hdot
557 end subroutine condense_vf_f
562 #ifdef PMC_USE_SUNDIALS 566 subroutine condense_jac(n_eqn, time, state_p, dDdot_dD, dDdot_dH, &
570 integer(kind=c_int),
intent(in) :: n_eqn
572 real(kind=c_double),
intent(in) :: time
574 type(c_ptr),
intent(in) :: state_p
576 real(kind=dp),
intent(out) :: dDdot_dD(n_eqn - 1)
578 real(kind=dp),
intent(out) :: dDdot_dH(n_eqn - 1)
580 real(kind=dp),
intent(out) :: dHdot_dD(n_eqn - 1)
582 real(kind=dp),
intent(out) :: dHdot_dH
584 real(kind=c_double),
pointer :: state(:)
589 call c_f_pointer(state_p, state, (/ n_eqn /))
597 inputs%H = state(n_eqn)
600 do i_part = 1,(n_eqn - 1)
601 inputs%D = state(i_part)
603 inputs%V_comp = (inputs%T &
609 dddot_dd(i_part) = outputs%dDdot_dD
610 dddot_dh(i_part) = outputs%dDdot_dH
611 dhdot_dd(i_part) = outputs%dHdoti_dD + outputs%dHdotenv_dD
612 dhdot_dh = dhdot_dh + outputs%dHdoti_dH
614 dhdot_dh = dhdot_dh + outputs%dHdotenv_dH
616 end subroutine condense_jac
621 #ifdef PMC_USE_SUNDIALS 625 subroutine condense_jac_solve_f(n_eqn, time, state_p, state_dot_p, &
626 rhs_p, gamma) bind(c)
629 integer(kind=c_int),
value,
intent(in) :: n_eqn
631 real(kind=c_double),
value,
intent(in) :: time
633 type(c_ptr),
value,
intent(in) :: state_p
635 type(c_ptr),
value,
intent(in) :: state_dot_p
637 type(c_ptr),
value,
intent(in) :: rhs_p
639 real(kind=c_double),
value,
intent(in) :: gamma
641 real(kind=c_double),
pointer :: state(:), state_dot(:), rhs(:)
642 real(kind=c_double) :: soln(n_eqn)
643 real(kind=dp) :: dDdot_dD(n_eqn - 1), dDdot_dH(n_eqn - 1)
644 real(kind=dp) :: dHdot_dD(n_eqn - 1), dHdot_dH
645 real(kind=dp) :: lhs_n, rhs_n
646 real(kind=c_double) :: residual(n_eqn)
647 real(kind=dp) :: rhs_norm, soln_norm, residual_norm
652 call condense_jac(n_eqn, time, state_p, dddot_dd, dddot_dh, &
655 call c_f_pointer(state_p, state, (/ n_eqn /))
656 call c_f_pointer(state_dot_p, state_dot, (/ n_eqn /))
657 call c_f_pointer(rhs_p, rhs, (/ n_eqn /))
660 lhs_n = 1d0 - gamma * dhdot_dh
662 do i_part = 1,(n_eqn - 1)
663 lhs_n = lhs_n - (- gamma * dddot_dh(i_part)) &
664 * (- gamma * dhdot_dd(i_part)) / (1d0 - gamma * dddot_dd(i_part))
665 rhs_n = rhs_n - (- gamma * dhdot_dd(i_part)) * rhs(i_part) &
666 / (1d0 - gamma * dddot_dd(i_part))
668 soln(n_eqn) = rhs_n / lhs_n
670 do i_part = 1,(n_eqn - 1)
671 soln(i_part) = (rhs(i_part) &
672 - (- gamma * dddot_dh(i_part)) * soln(n_eqn)) &
673 / (1d0 - gamma * dddot_dd(i_part))
680 residual(n_eqn) = sum(dhdot_dd * soln(1:(n_eqn-1))) &
681 + dhdot_dh * soln(n_eqn)
682 residual(1:(n_eqn-1)) = dddot_dd * soln(1:(n_eqn-1)) &
683 + dddot_dh * soln(n_eqn)
685 residual = rhs - (soln - gamma * residual)
686 rhs_norm = sqrt(sum(rhs**2))
687 soln_norm = sqrt(sum(soln**2))
688 residual_norm = sqrt(sum(residual**2))
689 write(0,*)
'rhs, soln, residual, residual/rhs = ', &
690 rhs_norm, soln_norm, residual_norm, residual_norm / rhs_norm
695 end subroutine condense_jac_solve_f
711 real(kind=dp) :: X, kappa, D_dry, D, g, dg_dD, a_w, daw_dD
712 integer :: newton_step
714 x = 4d0 *
const%water_molec_weight *
const%water_surf_eng &
715 / (
const%univ_gas_const * env_state%temp &
716 *
const%water_density)
725 do newton_step = 1,20
727 a_w = (d**3 - d_dry**3) / (d**3 + (kappa - 1d0) * d_dry**3)
728 daw_dd = 3d0 * d**2 * kappa * d_dry**3 &
729 / (d**3 + (kappa - 1d0) * d_dry**3)**2
730 g = env_state%rel_humid - a_w * exp(x / d)
731 dg_dd = - daw_dd * exp(x / d) + a_w * exp(x / d) * (x / d**2)
734 "convergence problem in equilibriation")
759 aero_state%valid_sort = .false.
765 aero_state%apa%particle(i_part))
integer, parameter pmc_condense_solver_init_y
Result code indicating failure to allocate y vector.
type(env_state_t) condense_saved_env_state_initial
Internal-use variable for storing the initial environment state during calls to the ODE solver...
logical, parameter condense_do_test_counts
Whether to print call-counts for helper routines during execution (for debugging only).
real(kind=dp) function aero_particle_solute_volume(aero_particle, aero_data)
Returns the total solute volume (m^3).
integer, parameter pmc_condense_solver_svtol
Result code indicating failure to set tolerances.
real(kind=dp) elemental function aero_data_diam2vol(aero_data, d)
Convert geometric diameter (m) to mass-equivalent volume (m^3).
integer, parameter pmc_condense_solver_init_cvode
Result code indicating failure to initialize the solver.
The env_state_t structure and associated subroutines.
integer, parameter pmc_condense_solver_init_cvode_mem
Result code indicating failure to create the solver.
Water condensation onto aerosol particles.
subroutine condense_equilib_particle(env_state, aero_data, aero_particle)
Determine the water equilibrium state of a single particle.
The aero_particle_t structure and associated subroutines.
integer, parameter pmc_condense_solver_success
Result code indicating successful completion.
elemental integer function aero_state_n_part(aero_state)
Return the current number of particles.
subroutine condense_equilib_particles(env_state, aero_data, aero_state)
Call condense_equilib_particle() on each particle in the aerosol to ensure that every particle has it...
real(kind=dp) function aero_particle_solute_kappa(aero_particle, aero_data)
Returns the average of the solute kappas (1).
real(kind=dp), dimension(:), allocatable condense_saved_kappa
Internal-use variable for storing the per-particle kappa values during calls to the ODE solver...
The aero_state_t structure and assocated subroutines.
real(kind=dp) function aero_particle_water_density(aero_data)
Returns the water density (kg/m^3).
logical, parameter condense_do_test_jac_solve
Whether to numerically test the Jacobian-solve function during execution (for debugging only)...
subroutine condense_particles(aero_state, aero_data, env_state_initial, env_state_final, del_t)
Do condensation to all the particles for a given time interval, including updating the environment to...
integer, save condense_count_vf
Internal-use variable for counting calls to the vector field subroutine.
subroutine aero_state_reweight(aero_state, aero_data, reweight_num_conc)
Reweight all particles after their constituent volumes have been altered.
Current environment state.
type(aero_data_t) condense_saved_aero_data
Internal-use variable for storing the aerosol data during calls to the ODE solver.
subroutine warn_assert_msg(code, condition_ok, warning_msg)
Prints a warning message if condition_ok is false.
subroutine die_msg(code, error_msg)
Error immediately.
int condense_solver(int neq, double *x_f, double *abstol_f, double reltol_f, double t_initial_f, double t_final_f)
Call the ODE solver.
elemental real(kind=dp) function aero_particle_diameter(aero_particle, aero_data)
Total diameter of the particle (m).
integer, parameter pmc_condense_solver_set_max_steps
Result code indicating failure to set maximum steps.
The current collection of aerosol particles.
character(len=pmc_util_convert_string_len) function real_to_string(val)
Convert a real to a string format.
The aero_data_t structure and associated subroutines.
Single aerosol particle data structure.
integer, parameter pmc_condense_solver_init_abstol
Result code indicating failure to allocate abstol vector.
subroutine aero_state_num_conc_for_reweight(aero_state, aero_data, reweight_num_conc)
Save the correct number concentrations for later use by aero_state_reweight().
real(kind=dp), dimension(:), allocatable condense_saved_d_dry
Internal-use variable for storing the per-particle dry diameters during calls to the ODE solver...
real(kind=dp) condense_saved_tdot
Internal-use variable for storing the rate of change of the temperature during calls to the ODE solve...
integer, parameter dp
Kind of a double precision real number.
Internal-use structure for storing the outputs from the rate-calculation function.
real(kind=dp), dimension(:), allocatable condense_saved_num_conc
Internal-use variable for storing the per-particle number concentrations during calls to the ODE solv...
integer, save condense_count_solve
Internal-use variable for counting calls to the Jacobian-solving subroutine.
real(kind=dp) elemental function aero_data_vol2diam(aero_data, v)
Convert mass-equivalent volume (m^3) to geometric diameter (m).
type(const_t), save const
Fixed variable for accessing the constant's values.
real(kind=dp) condense_saved_pdot
Internal-use variable for storing the rate of change of the pressure during calls to the ODE solver...
character(len=pmc_util_convert_string_len) function integer_to_string(val)
Convert an integer to a string format.
Aerosol material properties and associated data.
integer, parameter pmc_condense_solver_fail
Result code indicating failure of the solver.
Common utility subroutines.
subroutine condense_rates(inputs, outputs)
Compute the rate of change of particle diameter and relative humidity for a single particle...
real(kind=dp) function env_state_sat_vapor_pressure(env_state)
Computes the current saturation vapor pressure (Pa).