! Copyright (C) 2009-2013, 2016, 2017, 2021 Matthew West ! Licensed under the GNU General Public License version 2 or (at your ! option) any later version. See the file COPYING for details. !> \file !> The process program. !> Read NetCDF output files and process them. program process use pmc_output use pmc_stats character(len=PMC_MAX_FILENAME_LEN), parameter :: prefix & = "out/urban_plume" character(len=PMC_MAX_FILENAME_LEN) :: in_filename, out_filename type(bin_grid_t) :: diam_grid, bc_grid, sc_grid type(aero_data_t) :: aero_data type(aero_state_t) :: aero_state type(env_state_t) :: env_state integer :: ncid, index, repeat, i_index, i_repeat, n_index, n_repeat real(kind=dp) :: time, del_t, tot_num_conc, tot_mass_conc real(kind=dp) :: d_alpha, d_gamma, chi, d_alpha_a, d_gamma_a, chi_a character(len=PMC_UUID_LEN) :: uuid real(kind=dp), allocatable :: times(:), dry_diameters(:), num_concs(:), & dry_masses(:), masses(:), bc_masses(:), bc_fracs(:), & crit_rhs(:), scs(:), num_dist(:), & diam_bc_dist(:,:), diam_sc_dist(:,:) type(stats_1d_t) :: stats_num_dist, stats_d_alpha, stats_tot_num_conc, & stats_tot_mass_conc, stats_d_gamma, stats_chi, & stats_d_alpha_a, stats_d_gamma_a, stats_chi_a type(stats_2d_t) :: stats_diam_bc_dist, stats_diam_sc_dist character(len=AERO_NAME_LEN), allocatable :: mixing_state_groups(:,:) call pmc_mpi_init() call input_n_files(prefix, n_repeat, n_index) call bin_grid_make(diam_grid, BIN_GRID_TYPE_LOG, 180, 1d-9, 1d-3) call bin_grid_make(bc_grid, BIN_GRID_TYPE_LINEAR, 50, 0d0, 1d0) call bin_grid_make(sc_grid, BIN_GRID_TYPE_LOG, 50, 1d-4, 1d0) allocate(times(n_index)) allocate(mixing_state_groups(3, 4)) ! 3 groups, max 4 species per group mixing_state_groups(1,:) = ["OC ", "BC ", " ", " "] mixing_state_groups(2,:) = ["API1 ", "API2 ", "LIM1 ", "LIM2 "] mixing_state_groups(3,:) = ["SO4 ", "NO3 ", "NH4 ", " "] scs = [ real(kind=dp) :: ] ! silence compiler warnings bc_fracs = [ real(kind=dp) :: ] do i_index = 1,n_index do i_repeat = 1,n_repeat call make_filename(in_filename, prefix, ".nc", i_index, i_repeat) write(*,*) "Processing " // trim(in_filename) ! FIXME: add UUID check into input_state(), keyed off of index or ! time or something? Or init to "" and check if not this. call input_state(in_filename, index, time, del_t, repeat, & uuid, aero_data=aero_data, aero_state=aero_state, & env_state=env_state) times(i_index) = time dry_diameters = aero_state_dry_diameters(aero_state, aero_data) num_concs = aero_state_num_concs(aero_state, aero_data) num_dist = bin_grid_histogram_1d(diam_grid, dry_diameters, num_concs) call stats_1d_add(stats_num_dist, num_dist) tot_num_conc = sum(num_concs) call stats_1d_add_entry(stats_tot_num_conc, tot_num_conc, i_index) masses = aero_state_masses(aero_state, aero_data) tot_mass_conc = sum(masses * num_concs) call stats_1d_add_entry(stats_tot_mass_conc, tot_mass_conc, i_index) dry_masses = aero_state_masses(aero_state, aero_data, & exclude=(/"H2O"/)) bc_masses = aero_state_masses(aero_state, aero_data, & include=(/"BC"/)) bc_fracs = bc_masses / dry_masses diam_bc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & bc_grid, bc_fracs, num_concs) call stats_2d_add(stats_diam_bc_dist, diam_bc_dist) crit_rhs = aero_state_crit_rel_humids(aero_state, aero_data, & env_state) scs = crit_rhs - 1d0 diam_sc_dist = bin_grid_histogram_2d(diam_grid, dry_diameters, & sc_grid, scs, num_concs) call stats_2d_add(stats_diam_sc_dist, diam_sc_dist) call aero_state_mixing_state_metrics(aero_state, aero_data, & d_alpha, d_gamma, chi, groups=mixing_state_groups) call aero_state_mixing_state_metrics(aero_state, aero_data, & d_alpha_a, d_gamma_a, chi_a, exclude=(/"H2O"/)) call stats_1d_add_entry(stats_d_alpha, d_alpha, i_index) call stats_1d_add_entry(stats_d_gamma, d_gamma, i_index) call stats_1d_add_entry(stats_chi, chi, i_index) call stats_1d_add_entry(stats_d_alpha_a, d_alpha_a, i_index) call stats_1d_add_entry(stats_d_gamma_a, d_gamma_a, i_index) call stats_1d_add_entry(stats_chi_a, chi_a, i_index) end do 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, "1_urban_plume process") call bin_grid_output_netcdf(diam_grid, ncid, "diam", unit="m") call bin_grid_output_netcdf(bc_grid, ncid, "bc_frac", unit="1") call bin_grid_output_netcdf(sc_grid, ncid, "sc", unit="1") call stats_1d_output_netcdf(stats_num_dist, ncid, "num_dist", & dim_name="diam", unit="m^{-3}") call stats_1d_clear(stats_num_dist) call stats_2d_output_netcdf(stats_diam_bc_dist, ncid, "diam_bc_dist", & dim_name_1="diam", dim_name_2="bc_frac", unit="m^{-3}") call stats_2d_clear(stats_diam_bc_dist) call stats_2d_output_netcdf(stats_diam_sc_dist, ncid, "diam_sc_dist", & dim_name_1="diam", dim_name_2="sc", unit="m^{-3}") call stats_2d_clear(stats_diam_sc_dist) 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, "1_urban_plume 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="kg m^{-3}") call stats_1d_output_netcdf(stats_d_alpha, ncid, "d_alpha", & dim_name="time", unit="1") call stats_1d_output_netcdf(stats_d_gamma, ncid, & "d_gamma", dim_name="time", unit="1") call stats_1d_output_netcdf(stats_chi, ncid, "chi", & dim_name="time", unit="1") call stats_1d_output_netcdf(stats_d_alpha_a, ncid, "d_alpha_a", & dim_name="time", unit="1") call stats_1d_output_netcdf(stats_d_gamma_a, ncid, & "d_gamma_a", dim_name="time", unit="1") call stats_1d_output_netcdf(stats_chi_a, ncid, "chi_a", & dim_name="time", unit="1") call pmc_nc_close(ncid) call pmc_mpi_finalize() end program process