allow coefficient for additive coag kernel to be set from spec file (#…

…187)

* allow coefficient for additive coag kernel to be set from spec file

* change beta_1 to additive_kernel_coefficient

---------

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

src/coag_kernel.F90

@@ -486,6 +486,9 @@ subroutine spec_file_read_coag_kernel_type(file, coag_kernel_type)
     !!     for the constant kernel; \c brown for the Brownian kernel,
     !!     or \c zero for no coagulation
     !!
+    !! If \c coag_kernel is \c additive, the kernel coefficient needs to be
+    !! provided using the \c additive_kernel_coeff parameter
+    !!
     !! See also:
     !!   - \ref spec_file_format --- the input file text format
 

src/coag_kernel_additive.F90

@@ -19,9 +19,6 @@ module pmc_coag_kernel_additive
   use pmc_aero_data
   use pmc_aero_particle
 
-  !> Scaling coefficient for constant kernel.
-  real(kind=dp), parameter :: beta_1 = 1000d0
-
 contains
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -41,7 +38,8 @@ subroutine kernel_additive(aero_particle_1, aero_particle_2, &
     !> Coagulation kernel.
     real(kind=dp), intent(out) :: k
 
-    k = beta_1 * (aero_particle_volume(aero_particle_1) &
+    k = env_state%additive_kernel_coefficient &
+         * (aero_particle_volume(aero_particle_1) &
          + aero_particle_volume(aero_particle_2))
 
   end subroutine kernel_additive
@@ -64,7 +62,7 @@ subroutine kernel_additive_minmax(v1, v2, aero_data, env_state, k_min, k_max)
     !> Coagulation kernel maximum value.
     real(kind=dp), intent(out) :: k_max
 
-    k_min = beta_1 * (v1 + v2)
+    k_min = env_state%additive_kernel_coefficient * (v1 + v2)
     k_max = k_min
 
   end subroutine kernel_additive_minmax
@@ -77,8 +75,8 @@ end subroutine kernel_additive_minmax
   !! Given input paramaters \f$R\f$ and \f$N_0\f$ we let the mean
   !! volume be \f$v_\mu = \frac{4\pi}{3} R^3\f$ and define the
   !! rescaled times \f$\tau = N_0 v_\mu \beta_1 t\f$ and \f$T = 1 -
-  !! e^{-\tau}\f$, where \f$\beta_1\f$ is the fixed kernel scaling
-  !! parameter. Then the solution is
+  !! e^{-\tau}\f$, where \f$\beta_1\f$ is the fixed additive kernel scaling
+  !! coefficient. Then the solution is
   !! \f[
   !!     n(D,t) \ {\rm d}\ln D
   !!     = \frac{\pi}{2} D^3
@@ -131,7 +129,7 @@ subroutine soln_additive_exp(bin_grid, aero_data, time, num_conc, &
                / mean_vol))
        end do
     else
-       tau = num_conc * mean_vol * beta_1 * time
+       tau = num_conc * mean_vol * env_state%additive_kernel_coefficient * time
        T = 1d0 - exp(-tau)
        do k = 1,bin_grid_size(bin_grid)
           rat_v = aero_data_rad2vol(aero_data, bin_grid%centers(k)) / mean_vol

src/env_state.F90

@@ -73,6 +73,8 @@ module pmc_env_state
      real(kind=dp) :: solar_zenith_angle
      !> Box height (m).
      real(kind=dp) :: height
+     !> Scaling coefficient for additive coagulation kernel.
+     real(kind=dp) :: additive_kernel_coefficient 
   end type env_state_t
 
 contains
@@ -385,7 +387,8 @@ integer function pmc_mpi_pack_size_env_state(val)
          + pmc_mpi_pack_size_integer(val%start_day) &
          + pmc_mpi_pack_size_real(val%elapsed_time) &
          + pmc_mpi_pack_size_real(val%solar_zenith_angle) &
-         + pmc_mpi_pack_size_real(val%height)
+         + pmc_mpi_pack_size_real(val%height) &
+         + pmc_mpi_pack_size_real(val%additive_kernel_coefficient)
 
   end function pmc_mpi_pack_size_env_state
 
@@ -429,6 +432,7 @@ subroutine pmc_mpi_pack_env_state(buffer, position, val)
     call pmc_mpi_pack_real(buffer, position, val%elapsed_time)
     call pmc_mpi_pack_real(buffer, position, val%solar_zenith_angle)
     call pmc_mpi_pack_real(buffer, position, val%height)
+    call pmc_mpi_pack_real(buffer, position, val%additive_kernel_coefficient)
     call assert(464101191, &
          position - prev_position <= pmc_mpi_pack_size_env_state(val))
 #endif
@@ -476,6 +480,7 @@ subroutine pmc_mpi_unpack_env_state(buffer, position, val)
     call pmc_mpi_unpack_real(buffer, position, val%elapsed_time)
     call pmc_mpi_unpack_real(buffer, position, val%solar_zenith_angle)
     call pmc_mpi_unpack_real(buffer, position, val%height)
+    call pmc_mpi_unpack_real(buffer, position, val%additive_kernel_coefficient)
     call assert(205696745, &
          position - prev_position <= pmc_mpi_pack_size_env_state(val))
 #endif

src/partmc.F90

@@ -647,6 +647,7 @@ subroutine partmc_exact(file)
     !!
     !! do_coagulation yes              # whether to do coagulation (yes/no)
     !! kernel additive                 # Additive coagulation kernel
+    !! additive_kernel_coeff 1000d0                   # Additive kernel constant
     !! </pre>
 
     ! only serial code here

src/run_exact.F90

@@ -142,6 +142,10 @@ subroutine spec_file_read_run_exact(file, run_exact_opt, aero_data, &
     if (run_exact_opt%do_coagulation) then
        call spec_file_read_coag_kernel_type(file, &
             run_exact_opt%coag_kernel_type)
+       if (run_exact_opt%coag_kernel_type == COAG_KERNEL_TYPE_ADDITIVE) then
+          call spec_file_read_real(file, 'additive_kernel_coeff', &
+               env_state%additive_kernel_coefficient)
+       end if
     else
        run_exact_opt%coag_kernel_type = COAG_KERNEL_TYPE_INVALID
     end if

src/run_part.F90

@@ -592,6 +592,10 @@ subroutine spec_file_read_run_part(file, run_part_opt, aero_data, &
          run_part_opt%do_coagulation)
     if (run_part_opt%do_coagulation) then
        call spec_file_read_coag_kernel_type(file, run_part_opt%coag_kernel_type)
+       if (run_part_opt%coag_kernel_type == COAG_KERNEL_TYPE_ADDITIVE) then
+          call spec_file_read_real(file, 'additive_kernel_coeff', &
+               env_state_init%additive_kernel_coefficient)
+       end if
     else
        run_part_opt%coag_kernel_type = COAG_KERNEL_TYPE_INVALID
     end if

src/run_sect.F90

@@ -257,6 +257,10 @@ subroutine spec_file_read_run_sect(file, run_sect_opt, aero_data, &
     if (run_sect_opt%do_coagulation) then
        call spec_file_read_coag_kernel_type(file, &
             run_sect_opt%coag_kernel_type)
+       if (run_sect_opt%coag_kernel_type == COAG_KERNEL_TYPE_ADDITIVE) then
+          call spec_file_read_real(file, 'additive_kernel_coeff', &
+               env_state%additive_kernel_coefficient)
+       end if
     else
        run_sect_opt%coag_kernel_type = COAG_KERNEL_TYPE_INVALID
     end if

test/additive/run_exact.spec

@@ -32,3 +32,4 @@ start_day 1                     # start day of year (UTC)
 
 do_coagulation yes              # whether to do coagulation (yes/no)
 coag_kernel additive            # coagulation kernel
+additive_kernel_coeff 1000      # additive kernel constant

test/additive/run_part.spec

@@ -38,6 +38,7 @@ start_day 1                     # start day of year (UTC)
 
 do_coagulation yes              # whether to do coagulation (yes/no)
 coag_kernel additive            # coagulation kernel
+additive_kernel_coeff 1000      # additive kernel constant
 do_condensation no              # whether to do condensation (yes/no)
 do_mosaic no                    # whether to do MOSAIC (yes/no)
 do_nucleation no                # whether to do nucleation (yes/no)