Add immersion freezing #205

tangwhiap · 2024-11-27T05:58:05Z

Add the immersion freezing simulation codes, freezing test case, and freezing scenario.

remove old directories

jcurtis2

Just a few more things I've noted while taking a pass through for PyPartMC.

src/ice_nucleation.F90

src/env_state.F90

jcurtis2 · 2025-08-29T01:35:13Z

src/env_state.F90

+
+  !> compute saturated vapor pressure (units : Pa) with respective to water
+  !> Formula (10) from [Murphy & Koop, 2004]  (https://doi.org/10.1256/qj.04.94)
+  real(kind=dp) function env_state_saturated_vapor_pressure_water(T)


I also wonder why env_state isn't passed to this function given it is an env_state function?

Because this function only requires the temperature, T is the only necessary input. I thought you had previously suggested that I use T instead of env_state as the input—should I change it back to env_state?

Yes, you are right - I see where I suggested that. This was probably the simplest suggestion at the time as you had two functions that did the same thing (which isn't a great idea), one taking temperature and the other taking env_state.

I think that made this pretty confusing and in hindsight, we probably should apply a different fix. Not 100% sure what it should be at the moment.

The fact is, that all env_state functions take env_state as an input so this function is an outlier. So I think we need it to take env_state. However we also shouldn't have two functions that do exactly the same thing.

Two paths that I can think of:

Make an env_state in the theoretical_freezing.F90 code and then call those functions by passing an env_state that has it the temperature set. Then switch to have env_state as the input rather than temperature.

A function that takes env_state that is just a wrapper to call another function that just takes temperature (and this function name shouldn't start with env_state). However, I don't know that we have any cases in the code that does something quite like that. And I'm not quite sure what file this function would actually exist in (it doesn't quite feel like it belongs in util.F90).

Thanks for your reply! I prefer the second approach since it wouldn’t complicate theoretical_freezing.F90. However, it would still be somewhat of an outlier. We can discuss this further.

src/ice_nucleation.F90

slayoo · 2025-08-29T16:06:17Z

src/run_part.F90

+    if (run_part_opt%do_immersion_freezing .and. &
+         (run_part_opt%immersion_freezing_scheme_type .eq. &
+         IMMERSION_FREEZING_SCHEME_SINGULAR)) then
+       call ice_nucleation_singular_initialize(aero_state, aero_data)


this is now not yet called if using timestep or timeblock variants -> should be moved to the timestep logic

src/ice_nucleation.F90

Co-authored-by: Jeffrey Curtis <jcurtis2@illinois.edu>

src/run_part.F90

…tate_saturated_vapor_pressure_wrt_water

slayoo · 2025-09-15T15:36:23Z

closing and reopening to trigger codecov analysis in further commits

codecov · 2025-09-15T15:36:43Z

Codecov Report

❌ Patch coverage is 92.99363% with 22 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.40%. Comparing base (7e7ac4c) to head (6d962a3).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/aero_particle.F90	66.66%	10 Missing ⚠️
src/ice_nucleation.F90	94.64%	9 Missing ⚠️
src/aero_data.F90	87.50%	3 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           master     #205       +/-   ##
===========================================
+ Coverage        0   77.40%   +77.40%     
===========================================
  Files           0       55       +55     
  Lines           0     9427     +9427     
===========================================
+ Hits            0     7297     +7297     
- Misses          0     2130     +2130

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

…c files.

slayoo · 2025-09-25T01:18:24Z

CMakeLists.txt

  add_test(test_tchem ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh tchem)
 endif()
 add_test(test_weighting ${CMAKE_BINARY_DIR}/test_run/run_test_directory.sh weighting)
+#>>>>>>> 2341ef410d6f49f3169b8461b5fa8c89dbd3c7a2


this seems to be a merge-conflict leftover to be removed?

slayoo · 2025-09-25T01:19:53Z

doc/freezing/freezing.tex

would it be better to make it a part of the web-publishable doxygen html?

slayoo · 2025-09-25T01:21:27Z

scenarios/6_camp/monarch_mod37/custom_species.json

IIUC, the lack of these additions was not detected by CI, it would be great to add a mechanism at least parsing all the scenario input on CI

Yeah, this is a fairly common occurrence. The urban_plume scenarios are often fixed but the scenarios that often require more libraries (SUNDIALs for condense, CAMP) are often overlooked. This is often as simple as missing a spec file line change/addition. In this case, it would be a lack of adding things in aero_data.dat for most scenarios or a more special case in terms of the CAMP scenario. The CAMP test caught this, it was fixed in the test but not in the scenario til later.

I think that something that simply attempts to run a 0 second simulation will work in terms of reading in all inputs.

jcurtis2 · 2026-01-20T17:08:01Z

src/ice_nucleation.F90

+    aerosol_diameter =  aero_particle_dry_diameter(aero_particle, aero_data)
+    immersed_surface_area = const%pi * aerosol_diameter **2
+
+    total_vol = 0d0


Suggested change

total_vol = 0d0

jcurtis2 · 2026-01-20T17:09:14Z

src/ice_nucleation.F90

+       if (i_spec == aero_data%i_water) cycle
+       abifm_m = aero_data%abifm_m(i_spec)
+       abifm_c = aero_data%abifm_c(i_spec)
+       j_het = 10 ** (abifm_m  * (1 - a_w_ice) + abifm_c) * 10000


Suggested change

j_het = 10 ** (abifm_m * (1 - a_w_ice) + abifm_c) * 10000

j_het = 10d0 ** (abifm_m * (1d0 - a_w_ice) + abifm_c) * 10000

Also define what the 10000 is.

jcurtis2 · 2026-01-20T17:10:46Z

src/ice_nucleation.F90

+  !> Initialization for the sigular scheme, sampling the freezing temperature
+  !> for each particles.


Suggested change

!> Initialization for the sigular scheme, sampling the freezing temperature

!> for each particles.

!> Initialization for the singular scheme, sampling the freezing

!> temperature for each particle.

jcurtis2 · 2026-01-20T17:10:58Z

src/ice_nucleation.F90

+  !> for each particles.
+  subroutine ice_nucleation_singular_initialize(aero_state, aero_data, &
+       INAS_a, INAS_b)
+    implicit none


Suggested change

implicit none

jcurtis2 · 2026-01-20T17:15:18Z

src/ice_nucleation.F90

+  subroutine ice_nucleation_immersion_freezing_time_dependent(aero_state, &
+       aero_data, env_state, del_t, immersion_freezing_scheme_type, &
+       freezing_rate)
+
+    !> Aerosol state.
+    type(aero_state_t), intent(inout) :: aero_state
+    !> Aerosol data.
+    type(aero_data_t), intent(in) :: aero_data
+    !> Environment state.
+    type(env_state_t), intent(inout) :: env_state
+    !> Total time to integrate.
+    real(kind=dp), intent(in) :: del_t
+    !> Freezing rate (only used for the constant rate scheme).
+    real(kind=dp), intent(in) :: freezing_rate
+    integer, intent(in) :: immersion_freezing_scheme_type
+
+    integer :: i_part, i_bin, i_class, n_bins, n_class
+    real(kind=dp) :: a_w_ice, pis, pvs
+    real(kind=dp) :: p_freeze
+
+    real(kind=dp) :: p_freeze_max, radius_max, diameter_max
+
+    integer :: k_th, n_parts_in_bin
+    real(kind=dp) :: rand
+    real(kind=dp), allocatable :: H2O_masses(:), total_masses(:), &
+         H2O_frac(:)
+    integer :: i_spec_max
+    real(kind=dp) :: j_het_max
+    integer :: rand_geo
+
+
+    allocate(total_masses(aero_state_n_part(aero_state)))
+    allocate(H2O_masses(aero_state_n_part(aero_state)))
+    allocate(H2O_frac(aero_state_n_part(aero_state)))
+
+    call aero_state_sort(aero_state, aero_data)
+
+    total_masses = aero_state_masses(aero_state, aero_data)
+    H2O_masses = aero_state_masses(aero_state, aero_data, include=["H2O"])
+    H2O_frac = H2O_masses / total_masses
+    pvs = env_state_saturated_vapor_pressure_wrt_water(env_state%temp)
+    pis = env_state_saturated_vapor_pressure_wrt_ice(env_state%temp)
+    a_w_ice = pis / pvs
+    if (immersion_freezing_scheme_type == IMMERSION_FREEZING_SCHEME_ABIFM) then
+       call ABIFM_max_spec(aero_data, a_w_ice, i_spec_max, j_het_max)
+    endif
+
+    n_bins = aero_sorted_n_bin(aero_state%aero_sorted)
+    n_class = aero_sorted_n_class(aero_state%aero_sorted)
+
+    loop_bins: do i_bin = 1, n_bins
+       loop_classes: do i_class = 1, n_class
+          n_parts_in_bin = integer_varray_n_entry(&
+               aero_state%aero_sorted%size_class%inverse(i_bin, i_class))
+          radius_max = aero_state%aero_sorted%bin_grid%edges(i_bin + 1)
+          diameter_max = radius_max * 2
+          if (immersion_freezing_scheme_type == &
+               IMMERSION_FREEZING_SCHEME_ABIFM) then
+             p_freeze_max = ABIFM_Pfrz_max(diameter_max, aero_data, j_het_max, &
+                   del_t)
+          else if (immersion_freezing_scheme_type == &
+               IMMERSION_FREEZING_SCHEME_CONST) then
+             p_freeze_max = 1 - exp(freezing_rate * del_t)
+          endif
+
+          k_th = n_parts_in_bin + 1
+          loop_choosed_particles: do while(.TRUE.)
+             rand_geo = rand_geometric(p_freeze_max)
+             k_th = k_th - rand_geo
+             if (k_th <= 0) then
+                EXIT loop_choosed_particles
+             endif
+             i_part = aero_state%aero_sorted%size_class &
+                  %inverse(i_bin, i_class)%entry(k_th)
+             if (aero_state%apa%particle(i_part)%frozen) then
+                cycle
+             end if
+             if (H2O_frac(i_part) < const%imf_water_threshold) then
+                cycle
+             end if
+             if (immersion_freezing_scheme_type == &
+                  IMMERSION_FREEZING_SCHEME_ABIFM) then
+                p_freeze = ABIFM_Pfrz_particle( & 
+                     aero_state%apa%particle(i_part), aero_data, a_w_ice, del_t)
+                call warn_assert_msg(301184565, p_freeze <= p_freeze_max,&
+                     "p_freeze > p_freeze_max.")
+                rand = pmc_random()
+                if (rand < p_freeze / p_freeze_max) then
+                   aero_state%apa%particle(i_part)%frozen = .TRUE.
+                   aero_state%apa%particle(i_part)%den_ice = &
+                        const%reference_ice_density
+                   aero_state%apa%particle(i_part)%ice_shape_phi = 1d0
+
+                endif
+             else
+                aero_state%apa%particle(i_part)%frozen = .TRUE.
+                aero_state%apa%particle(i_part)%den_ice = &
+                     const%reference_ice_density
+                aero_state%apa%particle(i_part)%ice_shape_phi = 1d0
+             endif
+
+          enddo loop_choosed_particles
+       enddo loop_classes
+    enddo loop_bins
+
+    deallocate(total_masses)
+    deallocate(H2O_masses)
+    deallocate(H2O_frac)
+
+  end subroutine ice_nucleation_immersion_freezing_time_dependent


Fix all enddo and endif to be end do and end if

jcurtis2 · 2026-01-20T17:24:26Z

src/ice_nucleation.F90

+                   del_t)
+          else if (immersion_freezing_scheme_type == &
+               IMMERSION_FREEZING_SCHEME_CONST) then
+             p_freeze_max = 1 - exp(freezing_rate * del_t)


Suggested change

p_freeze_max = 1 - exp(freezing_rate * del_t)

p_freeze_max = 1d0 - exp(freezing_rate * del_t)

jcurtis2 · 2026-01-20T17:26:50Z

src/ice_nucleation.F90

+         H2O_frac(:)
+    real(kind=dp) :: rand
+
+


Suggested change

jcurtis2 · 2026-01-20T20:07:05Z

src/ice_nucleation.F90

+       else if (immersion_freezing_scheme_type == &
+            IMMERSION_FREEZING_SCHEME_CONST) then
+          p_freeze = 1 - exp(freezing_rate * del_t)
+       end if


p_freeze for IMMERSION_FREEZING_SCHEME_CONST doesn't rely on the particle at all. I think it would be worth considering moving outside the loop.

jcurtis2 · 2026-01-20T20:11:03Z

src/ice_nucleation.F90

+          else if (immersion_freezing_scheme_type == &
+               IMMERSION_FREEZING_SCHEME_CONST) then
+             p_freeze_max = 1 - exp(freezing_rate * del_t)
+          endif


If I understand correctly, consider refactoring p_freeze_max outside all these loops for the constant case.

jcurtis2 · 2026-01-20T20:16:58Z

src/aero_state.F90

+    aero_state_frozen_fraction = sum(pack(particle_num_concs, particle_frozen)) &
+         / sum(particle_num_concs)


It appears that I incorrectly suggested pack here and that we should avoid it when we aren't doing anything extra with the subset of particles. Sum with a mask is much more clear and efficient.

Suggested change

aero_state_frozen_fraction = sum(pack(particle_num_concs, particle_frozen)) &

/ sum(particle_num_concs)

aero_state_frozen_fraction = sum(particle_num_concs, mask=particle_frozen) &

/ sum(particle_num_concs)

tangwhiap and others added 30 commits August 29, 2023 15:24

change aero_data

a692f71

add freezing folder

44a075a

development of freezing process

70e4710

ice nucleation (immersion freezing) development

a31d570

freezing scenario

8f16bf0

modify freezing scenario

7ec13bc

Delete scenarios/7_freezing/run_scripts directory

2cf3381

remove old directories

Delete scenarios/7_freezing/aero_init_tools directory

ff65ded

remove old directories

Delete scenarios/7_freezing/visual_tools directory

f78b80a

remove old directories

Delete scenarios/7_freezing/out_sese_default1 directory

5c95091

remove old directories

Delete scenarios/7_freezing/data_temp directory

6143948

remove old directories

Delete scenarios/7_freezing/__pycache__ directory

efb3cee

remove old directories

remove old directories

ec418fe

Resolve merge conflicts by accepting deletions

8851a31

remove files

c647593

Delete scenarios/7_freezing/freezing_run.sh

8fa3e1c

Delete scenarios/7_freezing/hist.py

ef0fa57

Delete scenarios/7_freezing/make_pres.py

f713e41

Delete scenarios/7_freezing/make_temp.py

5bc0996

Delete scenarios/7_freezing/plot_temp1.py

ceaf677

Update .gitignore to exclude some freezing scripts

4b8a632

Merge branch 'imf' of github.com:tangwhiap/partmc into imf

bcad5c3

update and modify test cases for freezing

8688c83

add fortran codes for extracting data and theoretical calculation

c01b77c

delete old freezing codes

43a8ceb

remove the freezing timing output

d4b99bb

remove the forward ice growth parameters for condense subroutine

a3aa800

modify 1_run.sh for freezing scenario

40075be

modify visualization code

3ab971d

add freezing document (tex file)

89d1740

fortitude changes. break long lines

085ad81

jcurtis2 requested changes Aug 29, 2025

View reviewed changes

slayoo reviewed Aug 29, 2025

View reviewed changes

slayoo and others added 2 commits September 6, 2025 05:22

Merge branch 'master' into HEAD

dc8f24e

add abifm parameters to camp aerosol species

d8c7cf0

jcurtis2 reviewed Sep 12, 2025

View reviewed changes

src/ice_nucleation.F90 Show resolved Hide resolved

tangwhiap and others added 4 commits September 14, 2025 18:41

Simplify the freezing scenario, keep only exp8

00234bb

Resolve conflicts

73fe08f

Update src/ice_nucleation.F90

e79f406

Co-authored-by: Jeffrey Curtis <jcurtis2@illinois.edu>

Update src/ice_nucleation.F90

2d003f7

Co-authored-by: Jeffrey Curtis <jcurtis2@illinois.edu>

jcurtis2 reviewed Sep 15, 2025

View reviewed changes

src/run_part.F90 Outdated Show resolved Hide resolved

tangwhiap added 5 commits September 14, 2025 19:21

Change function name env_state_saturated_vapor_pressure_water-> env_s…

fbc9ca7

…tate_saturated_vapor_pressure_wrt_water

Remove duplicated function: env_state_saturated_vapor_pressure_water_2

6dd5b62

Modifications based on Jeff's comments

e7d0543

Modify the document for spec_file_read_immersion_freezing_scheme_type

6493ea7

Remove unreachable branch

6629ee5

slayoo closed this Sep 15, 2025

slayoo reopened this Sep 15, 2025

tangwhiap and others added 6 commits September 18, 2025 09:43

Add freezing test cases for singular and const

0fefef5

Move do_freezing_naive to spec file and add tests for naive algorithrm

5efbf4f

Add test_freezing_4.sh and test_freezing_5.sh, including relative spe…

03ecfb1

…c files.

reduce number of particles

b6032c1

reduce output frequency. reduce number of runs

ce1b709

fix the fixed rand_init issue

6d962a3

slayoo reviewed Sep 25, 2025

View reviewed changes

jcurtis2 mentioned this pull request Sep 25, 2025

Test/CI to parse all scenario inputs to ensure inputs remain valid #216

Open

Stop tracking 7_freezing scenario

6b13604

jcurtis2 requested changes Jan 20, 2026

View reviewed changes

	j_het = 10 ** (abifm_m * (1 - a_w_ice) + abifm_c) * 10000
	j_het = 10d0 ** (abifm_m * (1d0 - a_w_ice) + abifm_c) * 10000

		!> Initialization for the sigular scheme, sampling the freezing temperature
		!> for each particles.

	p_freeze_max = 1 - exp(freezing_rate * del_t)
	p_freeze_max = 1d0 - exp(freezing_rate * del_t)

		aero_state_frozen_fraction = sum(pack(particle_num_concs, particle_frozen)) &
		/ sum(particle_num_concs)

Add immersion freezing #205

Are you sure you want to change the base?

Add immersion freezing #205

Uh oh!

Conversation

tangwhiap commented Nov 27, 2024

Uh oh!

jcurtis2 left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

slayoo commented Sep 15, 2025

Uh oh!

codecov bot commented Sep 15, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Sep 15, 2025 •

edited

Loading