Add minimal modal run to model dry deposition

.dockerignore

@@ -25,3 +25,5 @@ scenarios/4_chamber/out
 scenarios/5_coag_brownian/out
 !scenarios/6_camp
 scenarios/6_camp/out
+!scenarios/7_drydep/
+scenarios/7_drydep/out

CMakeLists.txt

@@ -31,6 +31,7 @@ option(ENABLE_TCHEM "Enable TChem chemistry support" OFF)
 option(ENABLE_MPI "Enable MPI parallel support" OFF)
 option(ENABLE_SUNDIALS "Enable SUNDIALS solver for condensation support" OFF)
 option(ENABLE_C_SORT "Enable C sorting routines" OFF)
+option(ENABLE_QUADPACK "Enable QUADPACK for numerical integration (dry dep)" OFF)
 
 ######################################################################
 # CPack
@@ -179,6 +180,25 @@ if(ENABLE_SUNDIALS)
   add_definitions(-DPMC_USE_SUNDIALS)
 endif()
 
+######################################################################
+# QUADPACK
+
+if (ENABLE_QUADPACK)
+  find_path(QUADPACK_INCLUDE_DIR quadpack_double.mod
+    DOC "QUADPACK include directory"
+    PATHS $ENV{QUADPACK_HOME}/include)
+  find_library(QUADPACK_LIB quadpack
+    DOC "QUADPACK library"
+    PATHS $ENV{QUADPACK_HOME}/lib)
+  if(EXISTS "${QUADPACK_INCLUDE_DIR}/quadpack_double.mod" AND EXISTS "${QUADPACK_LIB}")
+      include_directories(${QUADPACK_INCLUDE_DIR})
+      add_definitions(-DPMC_USE_QUADPACK)
+  else()
+      message(WARNING "QUADPACK library not found - disabling QUADPACK support")
+      set(ENABLE_QUADPACK OFF)
+  endif()
+endif()
+
 ######################################################################
 # tests
 
@@ -236,14 +256,15 @@ add_library(partmclib src/aero_state.F90 src/integer_varray.F90
   src/nucleate.F90 src/condense.F90 src/fractal.F90 src/chamber.F90
   src/camp_interface.F90 src/photolysis.F90 src/sys.F90
   src/ice_nucleation.F90
+  src/run_modal.F90
   src/tchem_interface.F90
   src/aero_component.F90
   ${SUNDIALS_SRC} ${GSL_SRC} ${TCHEM_SRC}
   ${C_SORT_SRC})
 
 target_link_libraries(partmclib ${NETCDF_LIBS} ${SUNDIALS_LIBS}
   ${MOSAIC_LIB} ${GSL_LIBS} ${CAMP_LIB} ${TCHEM_LIB} ${YAML_LIB}
-  ${KOKKOS_LIB} ${KOKKOSKERNEL_LIB} ${TINES_LIB} ${CPP_LIB} ${LAPACK_LIB})
+  ${KOKKOS_LIB} ${KOKKOSKERNEL_LIB} ${TINES_LIB} ${CPP_LIB} ${LAPACK_LIB} ${QUADPACK_LIB})
 
 if (ENABLE_TCHEM)
    find_package(OpenMP REQUIRED)
@@ -431,3 +452,10 @@ add_executable(chamber_process
 target_link_libraries(chamber_process partmclib)
 
 ######################################################################
+# scenarios/7_drydep/drydep_modal_process
+
+add_executable(drydep_modal_process
+  scenarios/7_drydep/drydep_modal_process.F90)
+target_link_libraries(drydep_modal_process partmclib)
+
+######################################################################

scenarios/7_drydep/1_run_modal.sh

@@ -0,0 +1,12 @@
+#!/bin/sh
+
+# exit on error
+set -e
+# turn on command echoing
+set -v
+
+mkdir -p data
+
+../../build/partmc drydep_modal.spec
+
+# Now run ./2_process_drydep_modal.sh to process the data

scenarios/7_drydep/2_process_modal.sh

@@ -0,0 +1,10 @@
+#!/bin/sh
+
+# exit on error
+set -e
+# turn on command echoing
+set -v
+
+# The data should have already been generated by ./1_run_dry_dep_modal.sh
+
+../../build/drydep_modal_process

scenarios/7_drydep/aero_back.dat

@@ -0,0 +1,6 @@
+# time (s)
+# rate (s^{-1})
+# aerosol distribution filename
+time    0
+rate    0
+dist    aero_back_dist.dat

scenarios/7_drydep/aero_back_dist.dat

@@ -0,0 +1 @@
+# no background distribution

scenarios/7_drydep/aero_data.dat

@@ -0,0 +1,2 @@
+#     dens (kg/m^3)   ions in soln (1)    molec wght (kg/mole)   kappa (1)
+SO4            1800                0                   96d-3      0.65

scenarios/7_drydep/aero_emit.dat

@@ -0,0 +1,6 @@
+# time (s)
+# rate (s^{-1})
+# aerosol distribution filename   
+time    0
+rate    0 
+dist    aero_emit_dist.dat

scenarios/7_drydep/aero_emit_dist.dat

@@ -0,0 +1 @@
+# no emissions

scenarios/7_drydep/aero_init_comp.dat

@@ -0,0 +1,2 @@
+#       mass fractions
+SO4              1

scenarios/7_drydep/aero_init_dist.dat

@@ -0,0 +1,7 @@
+mode_name init_mode
+mass_frac aero_init_comp.dat      # composition proportions of species
+diam_type geometric               # type of diameter specified
+mode_type log_normal              # type of distribution
+num_conc 3.2e9                    # particle number concentration (#/m^3)
+geom_mean_diam 1e-8               # geometric mean diameter (m)
+log10_geom_std_dev 0.041          # log_10 of geometric std dev of diameter

scenarios/7_drydep/drydep_modal.spec

@@ -0,0 +1,33 @@
+run_type modal                       # modal run
+output_prefix data/modal_dpg_00001000000000000000_sig_2_5_emerson_grass
+
+t_max 28800                          # total simulation time (s)
+del_t 60                             # timestep (s)
+t_output 3600                        # output interval (0 disables) (s)
+t_progress 0                         # progress printing interval (0 disables) (s)
+
+n_bin 1000                           # number of bins (for processing purposes)
+d_min 4e-8                           # minimum diameter (m)
+d_max 2.5e-3                         # maximum diameter (m)
+
+gas_data gas_data.dat                # file containing gas data
+aerosol_data aero_data.dat           # file containing aerosol data
+do_fractal no                        # whether to do fractal treatment 
+aerosol_init aero_init_dist.dat      # aerosol initial condition file
+
+temp_profile temp.dat                # temperature profile file
+pressure_profile pres.dat            # pressure profile file
+height_profile height.dat            # height profile file
+gas_emissions gas_emit.dat           # gas emissions file
+gas_background gas_back.dat          # background gas concentrations file
+aero_emissions aero_emit.dat         # aerosol emissions file
+aero_background aero_back.dat        # aerosol background file
+loss_function drydep                 # loss function specification
+drydep_params drydep_params.dat      # dry deposition parameters
+
+rel_humidity 0.95                    # initial relative humidity (1)
+latitude 0                           # latitude (degrees, -90 to 90)
+longitude 0                          # longitude (degrees, -180 to 180)
+altitude 0                           # altitude (m)
+start_time 21600                     # start time (s since 00:00 UTC)
+start_day 200                        # start day of year (UTC)

scenarios/7_drydep/drydep_modal_process.F90

@@ -0,0 +1,110 @@
+!> Read NetCDF output files from a modal runn and process them.
+
+program process
+
+  use pmc_output
+  use pmc_stats
+  use pmc_aero_dist
+  use pmc_scenario
+
+  implicit none
+
+  character(len=PMC_MAX_FILENAME_LEN) :: prefix
+  character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename
+  type(bin_grid_t) :: bin_grid
+  type(aero_dist_t) :: aero_dist
+  type(aero_binned_t) :: aero_binned
+  type(aero_data_t) :: aero_data
+  type(env_state_t) :: env_state
+  type(gas_data_t) :: gas_data
+  type(gas_state_t) :: gas_state
+  type(scenario_t) :: scenario
+  integer :: ncid, index, i_mode, i_index, dum, n_index
+  real(kind=dp) :: time, del_t, tot_num_conc, density, tot_mass_conc
+  character(len=PMC_UUID_LEN) :: uuid
+  real(kind=dp), allocatable :: times(:), num_conc(:)
+  type(stats_1d_t) :: stats_tot_num_conc, stats_num_conc, stats_velocities_0, &
+                      stats_velocities_3, stats_tot_mass_conc
+  real(kind=dp), allocatable :: velocities_0(:), velocities_3(:)
+  character(len=PMC_MAX_FILENAME_LEN), allocatable :: file_list(:)
+
+  call pmc_mpi_init()
+
+  call get_command_argument(1, prefix)
+if (len_trim(prefix) == 0) &
+  call die_msg(192837465, "usage: drydep_modal_process <output_prefix>")
+
+  call input_filename_list(prefix, file_list)
+  n_index = size(file_list)
+
+  allocate(times(n_index))
+
+  do i_index = 1,n_index
+    call make_filename(in_filename, prefix, ".nc", i_index)
+    write(*,*) "Processing " // trim(in_filename)
+    call input_modal(in_filename, index, time, del_t, uuid, aero_dist=aero_dist, &
+        aero_binned=aero_binned, aero_data=aero_data, env_state=env_state, &
+        gas_data=gas_data, gas_state=gas_state, bin_grid=bin_grid, scenario=scenario)
+    times(i_index) = time
+
+    density = aero_data%density(1)
+
+    num_conc = aero_binned%num_conc * bin_grid%widths
+    tot_num_conc = sum(num_conc)
+
+    tot_mass_conc = sum(aero_binned%vol_conc(:,1) * bin_grid%widths) * aero_data%density(1)
+
+    if (.not. allocated(velocities_0)) allocate(velocities_0(aero_dist_n_mode(aero_dist)))
+    if (.not. allocated(velocities_3)) allocate(velocities_3(aero_dist_n_mode(aero_dist)))
+
+    call scenario_modal_drydep_velocities(scenario, aero_dist, 0.0d0, density, &
+      env_state, velocities_0)
+    call stats_1d_add(stats_velocities_0, velocities_0)
+    call scenario_modal_drydep_velocities(scenario, aero_dist, 3.0d0, density, &
+      env_state, velocities_3)
+    call stats_1d_add(stats_velocities_3, velocities_3)
+
+    call stats_1d_add(stats_num_conc, num_conc)
+    call stats_1d_add_entry(stats_tot_num_conc, tot_num_conc, i_index)
+
+    call stats_1d_add_entry(stats_tot_mass_conc, tot_mass_conc, i_index)
+
+    call make_filename(out_filename, prefix, "_process.nc", index)
+    write(*,*) "Writing " // trim(out_filename)
+    call pmc_nc_open_write(out_filename, ncid)
+    call pmc_nc_write_info(ncid, uuid, "7_drydep modal process")
+    call env_state_output_netcdf(env_state, ncid)
+    call aero_data_output_netcdf(aero_data, ncid)
+    call aero_dist_output_netcdf(aero_dist, ncid)
+    call aero_binned_output_netcdf(aero_binned, ncid, bin_grid, aero_data)
+
+    call stats_1d_output_netcdf(stats_num_conc, ncid, "num concs. per bin", &
+      dim_name="diam", unit="m^{-3}")
+    call stats_1d_output_netcdf(stats_velocities_0, ncid, "loss_velocities_0", &
+      dim_name="modes", unit="m s^{-1}")
+    call stats_1d_output_netcdf(stats_velocities_3, ncid, "loss_velocities_3", &
+      dim_name="modes", unit="m s^{-1}")
+
+    call aero_binned_zero(aero_binned)
+
+    call stats_1d_clear(stats_num_conc)
+    call stats_1d_clear(stats_velocities_0)
+    call stats_1d_clear(stats_velocities_3)
+
+    call pmc_nc_close(ncid) 
+  end do
+
+  call make_filename(out_filename, prefix, "_process.nc")
+  write(*,*) "Writing " // trim(out_filename)
+  call pmc_nc_open_write(out_filename, ncid)
+  call pmc_nc_write_info(ncid, uuid, "7_drydep modal process")
+  call pmc_nc_write_real_1d(ncid, times, "time", dim_name="time", unit="s")
+  call stats_1d_output_netcdf(stats_tot_num_conc, ncid, "tot_num_conc", &
+      dim_name="time", unit="m^{-3}")
+  call stats_1d_output_netcdf(stats_tot_mass_conc, ncid, "tot_mass_conc", &
+      dim_name="time", unit="ug m^{-3}")
+  call pmc_nc_close(ncid)
+
+  call pmc_mpi_finalize()
+
+end program process

scenarios/7_drydep/drydep_params.dat

@@ -0,0 +1,12 @@
+z_ref 10.0   # reference height (m)
+u_mean 2.0   # mean wind speed at the reference height (m)
+z_rough 3.8  # surface roughness length (m)
+A 0.0024     # characteristic radius of collectors (m)
+alpha 5.0    # parameter used in impaction efficiency
+eps_0 6.0    # emprical constant used in surface resistance
+gamma .756   # exponent for Brownian diffusion collection efficiency
+C_B 0.82     # coefficient for Brownian diffusion
+C_IN 92.5    # coefficient for interception
+C_IM 10.4    # coefficient for impaction
+nu 11.8      # exponent for interception collection efficiency
+beta 12.7    # exponent for impaction collection efficiency

scenarios/7_drydep/gas_back.dat

@@ -0,0 +1,79 @@
+# time (s)
+# rate (s^{-1})
+# concentrations (ppb)
+time            0
+rate         1.5e-5
+NO           0.1E+00
+NO2          1.0E+00
+NO3          0.0E+00
+N2O5         0.0E+00
+HONO         0.0E+00
+HNO3         1.0E+00
+HNO4         0.0E+00
+O3           5.0E+01
+O1D          0.0E+00
+O3P          0.0E+00
+OH           0.0E+00
+HO2          0.0E+00
+H2O2         1.1E+00
+CO           2.1E+02
+SO2          0.8E+00
+H2SO4        0.0E+00
+NH3          0.5E+00
+HCl          0.7E+00
+CH4          2.2E+03
+C2H6         1.0E+00
+CH3O2        0.0E+00
+ETHP         0.0E+00
+HCHO         1.2E+00
+CH3OH        1.2E-01
+CH3OOH       0.5E+00
+ETHOOH       0.0E+00
+ALD2         1.0E+00
+HCOOH        0.0E+00
+PAR          2.0E+00
+AONE         1.0E+00
+MGLY         0.0E+00
+ETH          0.2E+00
+OLET         2.3E-02
+OLEI         3.1E-04
+TOL          0.1E+00
+XYL          0.1E+00
+CRES         0.0E+00
+TO2          0.0E+00
+CRO          0.0E+00
+OPEN         0.0E+00
+ONIT         0.1E+00
+PAN          0.8E+00
+RCOOH        0.2E+00
+ROOH         2.5E-02
+C2O3         0.0E+00
+RO2          0.0E+00
+ANO2         0.0E+00
+NAP          0.0E+00
+ARO1         0.0E+00
+ARO2         0.0E+00
+ALK1         0.0E+00
+OLE1         0.0E+00
+XO2          0.0E+00
+XPAR         0.0E+00
+ISOP         0.5E+00
+API          0.0E+00
+LIM          0.0E+00
+API1         0.0E+00
+API2         0.0E+00
+LIM1         0.0E+00
+LIM2         0.0E+00
+ISOPRD       0.0E+00
+ISOPP        0.0E+00
+ISOPN        0.0E+00
+ISOPO2       0.0E+00
+DMS          0.0E+00
+MSA          0.0E+00
+DMSO         0.0E+00
+DMSO2        0.0E+00
+CH3SO2H      0.0E+00
+CH3SCH2OO    0.0E+00
+CH3SO2       0.0E+00
+CH3SO3       0.0E+00
+CH3SO2OO     0.0E+00