diff --git a/src/aero_rep_data.F90 b/src/aero_rep_data.F90 index 160d758b4..c68d4081f 100644 --- a/src/aero_rep_data.F90 +++ b/src/aero_rep_data.F90 @@ -276,7 +276,8 @@ end function get_size !> Get a list of unique names for each element on the !! \c camp_camp_state::camp_state_t::state_var array for this aerosol !! representation. - function unique_names(this, phase_name, tracer_type, spec_name) + function unique_names(this, phase_name, tracer_type, spec_name, & + phase_is_at_surface) use camp_util, only : string_t, i_kind import :: aero_rep_data_t @@ -290,6 +291,8 @@ function unique_names(this, phase_name, tracer_type, spec_name) integer(kind=i_kind), optional, intent(in) :: tracer_type !> Aerosol-phase species name character(len=*), optional, intent(in) :: spec_name + !> Indicates if aerosol phase is at the surface of particle + logical, optional, intent(in) :: phase_is_at_surface end function unique_names @@ -342,7 +345,7 @@ end function spec_name !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Get the number of instances of a specified aerosol phase - function num_phase_instances(this, phase_name) + function num_phase_instances(this, phase_name, is_at_surface) use camp_util, only : i_kind import :: aero_rep_data_t @@ -352,6 +355,8 @@ function num_phase_instances(this, phase_name) class(aero_rep_data_t), intent(in) :: this !> Aerosol phase name character(len=*), intent(in) :: phase_name + !> Indicates if aerosol phase is at the surface of particle + logical, intent(in), optional :: is_at_surface end function num_phase_instances @@ -506,7 +511,8 @@ function phase_ids(this, phase_name, is_at_surface) integer(kind=i_kind) :: num_instances, i_instance, i_phase - num_instances = this%num_phase_instances(phase_name) + num_instances = this%num_phase_instances(phase_name, is_at_surface) + allocate(phase_ids(num_instances)) if (present(is_at_surface)) then if (is_at_surface) then @@ -537,6 +543,7 @@ function phase_ids(this, phase_name, is_at_surface) end if end do end if + call assert(642387392, num_instances == i_instance-1) end function phase_ids diff --git a/src/aero_reps/aero_rep_modal_binned_mass.F90 b/src/aero_reps/aero_rep_modal_binned_mass.F90 index 7ae58d7c3..afb405bc8 100644 --- a/src/aero_reps/aero_rep_modal_binned_mass.F90 +++ b/src/aero_reps/aero_rep_modal_binned_mass.F90 @@ -734,7 +734,8 @@ end function get_size !! !! ... and for modes are: !! - "mode name.phase name.species name" - function unique_names(this, phase_name, tracer_type, spec_name) + function unique_names(this, phase_name, tracer_type, spec_name, & + phase_is_at_surface) !> List of unique names type(string_t), allocatable :: unique_names(:) @@ -746,6 +747,8 @@ function unique_names(this, phase_name, tracer_type, spec_name) integer(kind=i_kind), optional, intent(in) :: tracer_type !> Aerosol-phase species name character(len=*), optional, intent(in) :: spec_name + !> Aerosol-phase species is at surface + logical, optional, intent(in) :: phase_is_at_surface integer(kind=i_kind) :: num_spec, i_spec, j_spec, i_phase, j_phase, & i_section, i_bin @@ -764,6 +767,12 @@ function unique_names(this, phase_name, tracer_type, spec_name) if (phase_name.ne.curr_phase_name) cycle end if + ! Filter by phase is at surface + if (present(phase_is_at_surface)) then + if (phase_is_at_surface .neqv. & + this%aero_phase_is_at_surface(i_phase)) cycle + end if + ! Filter by spec name and/or tracer type if (present(spec_name).or.present(tracer_type)) then spec_names = this%aero_phase(i_phase)%val%get_species_names() @@ -810,6 +819,15 @@ function unique_names(this, phase_name, tracer_type, spec_name) end if end if + ! Filter by phase is at surface + if (present(phase_is_at_surface)) then + if (phase_is_at_surface .neqv. & + this%aero_phase_is_at_surface(i_phase)) then + i_phase = i_phase + NUM_BINS_(i_section) + cycle + end if + end if + ! Get the species names in this phase spec_names = this%aero_phase(i_phase)%val%get_species_names() @@ -941,7 +959,7 @@ end function spec_name !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Get the number of instances of a specified aerosol phase. - function num_phase_instances(this, phase_name) + function num_phase_instances(this, phase_name, is_at_surface) !> Number of instances of the aerosol phase integer(kind=i_kind) :: num_phase_instances @@ -949,15 +967,37 @@ function num_phase_instances(this, phase_name) class(aero_rep_modal_binned_mass_t), intent(in) :: this !> Aerosol phase name character(len=*), intent(in) :: phase_name + !> Indicates if aerosol phase is at the surface of particle + logical, intent(in), optional :: is_at_surface integer(kind=i_kind) :: i_phase num_phase_instances = 0 - do i_phase = 1, size(this%aero_phase) - if (this%aero_phase(i_phase)%val%name().eq.phase_name) then - num_phase_instances = num_phase_instances + 1 + if (present(is_at_surface)) then + if (is_at_surface) then + do i_phase = 1, size(this%aero_phase) + if (this%aero_phase(i_phase)%val%name().eq.phase_name .and. & + this%aero_phase_is_at_surface(i_phase)) then + num_phase_instances = num_phase_instances + 1 + end if + end do + else + do i_phase = 1, size(this%aero_phase) + if (this%aero_phase(i_phase)%val%name().eq.phase_name .and. & + .not. this%aero_phase_is_at_surface(i_phase)) then + call die_msg(507144607, "Species must exist at surface "// & + "in modal/binned mass aerosol representation '"// & + this%rep_name//"'") + end if + end do end if - end do + else + do i_phase = 1, size(this%aero_phase) + if (this%aero_phase(i_phase)%val%name().eq.phase_name) then + num_phase_instances = num_phase_instances + 1 + end if + end do + end if end function num_phase_instances diff --git a/src/aero_reps/aero_rep_single_particle.F90 b/src/aero_reps/aero_rep_single_particle.F90 index 274920ece..7d364a66c 100644 --- a/src/aero_reps/aero_rep_single_particle.F90 +++ b/src/aero_reps/aero_rep_single_particle.F90 @@ -375,7 +375,6 @@ subroutine initialize(this, aero_phase_set, spec_state_id) curr_phase = 1 do i_layer = 1, size(ordered_layer_id) j_layer = ordered_layer_id(i_layer) - ! Loop through the phases and make sure they exist call phases(j_layer)%val_%iter_reset() do i_phase = 1, phases(j_layer)%val_%size() @@ -387,7 +386,6 @@ subroutine initialize(this, aero_phase_set, spec_state_id) layer_names_unordered(j_layer)%string// & "' in single-particle layer aerosol representation '"// & this%rep_name//"'") - ! find phase and set pointer and indices found = .false. do j_phase = 1, size(aero_phase_set) @@ -411,6 +409,7 @@ subroutine initialize(this, aero_phase_set, spec_state_id) exit end if end do + call assert(373124707, found) call phases(j_layer)%val_%iter_next() end do @@ -485,7 +484,8 @@ end function per_particle_size !! For a single particle representation, the unique names will be a 'P' !! followed by the computational particle number, a '.', the phase name, !! another '.', and the species name. - function unique_names(this, phase_name, tracer_type, spec_name) + function unique_names(this, phase_name, tracer_type, spec_name, & + phase_is_at_surface) use camp_util, only : integer_to_string !> List of unique names @@ -498,6 +498,8 @@ function unique_names(this, phase_name, tracer_type, spec_name) integer(kind=i_kind), optional, intent(in) :: tracer_type !> Aerosol-phase species name character(len=*), optional, intent(in) :: spec_name + !> Indicates if aerosol phase is at the surface of particle + logical, optional, intent(in) :: phase_is_at_surface integer :: i_particle, i_layer, i_phase, i_spec, j_spec integer :: num_spec, curr_tracer_type @@ -518,6 +520,10 @@ function unique_names(this, phase_name, tracer_type, spec_name) if (present(phase_name)) then if(phase_name .ne. this%aero_phase(i_phase)%val%name()) cycle end if + if (present(phase_is_at_surface)) then + if (phase_is_at_surface .neqv. & + this%aero_phase_is_at_surface(i_phase)) cycle + end if if (present(spec_name) .or. present(tracer_type)) then spec_names = this%aero_phase(i_phase)%val%get_species_names() do j_spec = 1, size(spec_names) @@ -547,6 +553,10 @@ function unique_names(this, phase_name, tracer_type, spec_name) if (present(phase_name)) then if(phase_name .ne. this%aero_phase(i_phase)%val%name()) cycle end if + if (present(phase_is_at_surface)) then + if (phase_is_at_surface .neqv. & + this%aero_phase_is_at_surface(i_phase)) cycle + end if spec_names = this%aero_phase(i_phase)%val%get_species_names() do j_spec = 1, this%aero_phase(i_phase)%val%size() if (present(spec_name)) then @@ -638,25 +648,39 @@ end function spec_name !> Get the number of instances of a specified aerosol phase. In the single !! particle representation with layers, a phase can exist in multiple layers !! in one particle. - integer(kind=i_kind) function num_phase_instances(this, phase_name) + integer(kind=i_kind) function num_phase_instances(this, phase_name, & + is_at_surface) !> Aerosol representation data class(aero_rep_single_particle_t), intent(in) :: this !> Aerosol phase name character(len=*), intent(in) :: phase_name + !> Indicates if aerosol phase is at the surface of particle + logical, intent(in), optional :: is_at_surface integer(kind=i_kind) :: i_phase, i_layer, phase_index num_phase_instances = 0 - phase_index = 0 - do i_layer = 1, NUM_LAYERS_ - do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer) - if (this%aero_phase(i_phase)%val%name() .eq. phase_name) then - phase_index = phase_index + 1 - end if + if (present(is_at_surface)) then + do i_layer = 1, NUM_LAYERS_ + do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer) + if (this%aero_phase(i_phase)%val%name() .eq. phase_name .and. & + this%aero_phase_is_at_surface(i_phase) .eqv. is_at_surface) then + num_phase_instances = num_phase_instances + 1 + end if + end do end do - end do - num_phase_instances = phase_index * MAX_PARTICLES_ + + else + do i_layer = 1, NUM_LAYERS_ + do i_phase = LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer) + if (this%aero_phase(i_phase)%val%name() .eq. phase_name) then + num_phase_instances = num_phase_instances + 1 + end if + end do + end do + end if + num_phase_instances = num_phase_instances * MAX_PARTICLES_ end function num_phase_instances diff --git a/src/aero_reps/aero_rep_single_particle.c b/src/aero_reps/aero_rep_single_particle.c index 047d1badc..5e7902de6 100644 --- a/src/aero_reps/aero_rep_single_particle.c +++ b/src/aero_reps/aero_rep_single_particle.c @@ -296,8 +296,8 @@ void aero_rep_single_particle_get_aero_phase_mass__kg_m3( double *state = (double *)(model_data->grid_cell_state); state += i_part * PARTICLE_STATE_SIZE_ + PHASE_STATE_ID_(i_layer,i_phase); double mw; - aero_phase_get_mass__kg_m3(model_data, aero_phase_idx, state, - aero_phase_mass, &mw, partial_deriv, NULL); + aero_phase_get_mass__kg_m3(model_data, PHASE_MODEL_DATA_ID_(i_layer,i_phase), + state, aero_phase_mass, &mw, partial_deriv, NULL); if (partial_deriv) partial_deriv += PHASE_NUM_JAC_ELEM_(i_layer,i_phase); } else if (partial_deriv) { for (int i_spec = 0; i_spec < PHASE_NUM_JAC_ELEM_(i_layer,i_phase); ++i_spec) @@ -345,8 +345,8 @@ void aero_rep_single_particle_get_aero_phase_avg_MW__kg_mol( double *state = (double *)(model_data->grid_cell_state); state += i_part * PARTICLE_STATE_SIZE_ + PHASE_STATE_ID_(i_layer,i_phase); double mass; - aero_phase_get_mass__kg_m3(model_data, aero_phase_idx, state, &mass, - aero_phase_avg_MW, NULL, partial_deriv); + aero_phase_get_mass__kg_m3(model_data, PHASE_MODEL_DATA_ID_(i_layer,i_phase), + state, &mass, aero_phase_avg_MW, NULL, partial_deriv); if (partial_deriv) partial_deriv += PHASE_NUM_JAC_ELEM_(i_layer,i_phase); } else if (partial_deriv) { for (int i_spec = 0; i_spec < PHASE_NUM_JAC_ELEM_(i_layer,i_phase); ++i_spec) @@ -423,6 +423,8 @@ void aero_rep_single_particle_print(int *aero_rep_int_data, printf("\nParticle state size: %d", PARTICLE_STATE_SIZE_); for(int i_layer = 0; i_layer < NUM_LAYERS_; ++i_layer){ printf("\nLayer: %d", i_layer); + printf("\n Start phase: %d End phase: %d", LAYER_PHASE_START_(i_layer), LAYER_PHASE_END_(i_layer)); + printf("\n Number of phases: %d", NUM_PHASES_(i_layer)); printf("\n\n - Phases -"); for (int i_phase = 0; i_phase < NUM_PHASES_(i_layer); ++i_phase) { printf("\n state id: %d model data id: %d num Jac elements: %d", diff --git a/src/camp_solver.c b/src/camp_solver.c index 9f5ff6da2..bd3f33ff5 100644 --- a/src/camp_solver.c +++ b/src/camp_solver.c @@ -199,7 +199,7 @@ void *solver_new(int n_state_var, int n_cells, int *var_type, int n_rxn, // If there are no reactions, flag the solver not to run sd->no_solve = (n_rxn == 0); - // Allocate space for the aerosol phase data and st the number + // Allocate space for the aerosol phase data and set the number // of aerosol phases (including one int for the number of // phases) sd->model_data.aero_phase_int_data = diff --git a/src/rxns/rxn_HL_phase_transfer.F90 b/src/rxns/rxn_HL_phase_transfer.F90 index 0557653b4..a954f86c8 100644 --- a/src/rxns/rxn_HL_phase_transfer.F90 +++ b/src/rxns/rxn_HL_phase_transfer.F90 @@ -238,9 +238,11 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells) ! Get the unique names in this aerosol representation for the ! partitioning species and aerosol-phase water unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = aero_spec_name) + phase_name = phase_name, spec_name = aero_spec_name, & + phase_is_at_surface = .true.) unique_water_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = water_name) + phase_name = phase_name, spec_name = water_name, & + phase_is_at_surface = .true.) ! Skip aerosol representations that do not contain this phase if (.not.allocated(unique_spec_names)) cycle @@ -311,9 +313,11 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells) ! Get the unique names in this aerosol representation for the ! partitioning species and aerosol-phase water unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = aero_spec_name) + phase_name = phase_name, spec_name = aero_spec_name, & + phase_is_at_surface = .true.) unique_water_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = water_name) + phase_name = phase_name, spec_name = water_name, & + phase_is_at_surface = .true.) ! Get the phase ids for this aerosol phase phase_ids = aero_rep(i_aero_rep)%val%phase_ids(phase_name, is_at_surface=.true.) diff --git a/src/rxns/rxn_SIMPOL_phase_transfer.F90 b/src/rxns/rxn_SIMPOL_phase_transfer.F90 index 8e0d8725b..14f134086 100644 --- a/src/rxns/rxn_SIMPOL_phase_transfer.F90 +++ b/src/rxns/rxn_SIMPOL_phase_transfer.F90 @@ -245,7 +245,8 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells) ! Get the unique names in this aerosol representation for the ! partitioning species unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = aero_spec_name) + phase_name = phase_name, spec_name = aero_spec_name, & + phase_is_at_surface = .true.) ! Skip aerosol representations that do not contain this phase if (.not.allocated(unique_spec_names)) cycle @@ -306,12 +307,14 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells) ! Get the unique names in this aerosol representation for the ! partitioning species unique_spec_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = aero_spec_name) + phase_name = phase_name, spec_name = aero_spec_name, & + phase_is_at_surface = .true.) ! Find the corresponding activity coefficients, if specified if (has_act_coeff) then unique_act_names = aero_rep(i_aero_rep)%val%unique_names( & - phase_name = phase_name, spec_name = act_name) + phase_name = phase_name, spec_name = act_name, & + phase_is_at_surface = .true.) call assert_msg(236251734, size(unique_act_names).eq. & size(unique_spec_names), & "Mismatch of SIMPOL species and activity coeffs"// & diff --git a/test/unit_aero_rep_data/test_aero_rep_modal_binned_mass.F90 b/test/unit_aero_rep_data/test_aero_rep_modal_binned_mass.F90 index 49243871f..81851ca2f 100644 --- a/test/unit_aero_rep_data/test_aero_rep_modal_binned_mass.F90 +++ b/test/unit_aero_rep_data/test_aero_rep_modal_binned_mass.F90 @@ -110,6 +110,9 @@ logical function build_aero_rep_data_set_test() integer(kind=i_kind) :: i_spec, j_spec, i_phase, spec_id character(len=:), allocatable :: rep_name, spec_name, phase_name type(string_t), allocatable :: file_list(:), unique_names(:) + character(len=:), allocatable :: phase_name_test + integer, allocatable :: last_phase_id_correct(:) + integer, dimension(9) :: last_phase_id #ifdef CAMP_USE_MPI type(string_t), allocatable :: rep_names(:) type(aero_rep_factory_t) :: aero_rep_factory @@ -159,6 +162,28 @@ logical function build_aero_rep_data_set_test() call assert(208973841, unique_names(25)%string.eq."binned aerosol.6.my test phase one.species b") call assert(386300586, unique_names(47)%string.eq."binned aerosol.8.my last test phase.species b") + phase_name_test = "my last test phase" + last_phase_id(1) = 3 + last_phase_id(2) = 12 + last_phase_id(3) = 13 + last_phase_id(4) = 14 + last_phase_id(5) = 15 + last_phase_id(6) = 16 + last_phase_id(7) = 17 + last_phase_id(8) = 18 + last_phase_id(9) = 19 + !check values + last_phase_id_correct = aero_rep%phase_ids(phase_name_test, is_at_surface = .true.) + call assert(087831392, last_phase_id(1) .eq. last_phase_id_correct(1)) + call assert(597703093, last_phase_id(2) .eq. last_phase_id_correct(2)) + call assert(803195649, last_phase_id(3) .eq. last_phase_id_correct(3)) + call assert(516324200, last_phase_id(4) .eq. last_phase_id_correct(4)) + call assert(160854212, last_phase_id(5) .eq. last_phase_id_correct(5)) + call assert(626595159, last_phase_id(6) .eq. last_phase_id_correct(6)) + call assert(886729286, last_phase_id(7) .eq. last_phase_id_correct(7)) + call assert(787382607, last_phase_id(8) .eq. last_phase_id_correct(8)) + call assert(132691589, last_phase_id(9) .eq. last_phase_id_correct(9)) + ! Set the species concentrations phase_name = "my test phase one" spec_name = "species a" diff --git a/test/unit_aero_rep_data/test_aero_rep_single_particle.F90 b/test/unit_aero_rep_data/test_aero_rep_single_particle.F90 index f9d6d8b5d..31c8bd77e 100644 --- a/test/unit_aero_rep_data/test_aero_rep_single_particle.F90 +++ b/test/unit_aero_rep_data/test_aero_rep_single_particle.F90 @@ -28,7 +28,7 @@ program camp_test_aero_rep_data implicit none ! Test computational particle - integer(kind=i_kind), parameter :: TEST_PARTICLE = 2 + integer(kind=i_kind), parameter :: TEST_PARTICLE_1 = 2 ! Total computational particles integer(kind=i_kind), parameter :: NUM_COMP_PARTICLES = 3 @@ -36,8 +36,11 @@ program camp_test_aero_rep_data ! Number of aerosol phases per particle integer(kind=i_kind), parameter :: NUM_AERO_PHASE = 4 - ! Index for the test phase (test-particle phase 2) - integer(kind=i_kind), parameter :: AERO_PHASE_IDX = ((TEST_PARTICLE-1)*NUM_AERO_PHASE+2) + ! Index for the test phase + ! (test-particle 1 phase 2 : middle layer, jam) + integer(kind=i_kind), parameter :: AERO_PHASE_IDX_1 = ((TEST_PARTICLE_1-1)*NUM_AERO_PHASE+2) + ! (test-particle 1 phase 4 : top layer, top bread) + integer(kind=i_kind), parameter :: AERO_PHASE_IDX_2 = ((TEST_PARTICLE_1-1)*NUM_AERO_PHASE+4) ! Number of expected Jacobian elements for the test phase integer(kind=i_kind), parameter :: NUM_JAC_ELEM = 12 @@ -118,8 +121,6 @@ subroutine test_ordered_layer_ids() type(string_t), dimension(4) :: layer_names, correct_layer_names type(string_t), dimension(4) :: cover_names - !type(string_t), allocatable :: layer_names(:), correct_layer_names(:) - !type(string_t), allocatable :: cover_names(:) integer, allocatable :: ordered_ids(:) ! Assemble input arguments @@ -136,8 +137,6 @@ subroutine test_ordered_layer_ids() correct_layer_names(3)%string = 'jam' correct_layer_names(4)%string = 'top bread' - !layer_names(1)%string = 'one layer' - !correct_layer_names(1)%string = 'one layer' ! Call the function and enter inputs ordered_ids = ordered_layer_ids(layer_names, cover_names) ! check individual values: @@ -157,9 +156,12 @@ subroutine test_config_read() type(camp_core_t), pointer :: camp_core class(aero_rep_data_t), pointer :: aero_rep - type(string_t), allocatable :: file_list(:), unique_names(:) + type(string_t), allocatable :: file_list(:), unique_names(:), unique_names_surface(:) character(len=:), allocatable :: rep_name, phase_name_test integer :: i_name, num_bread, max_part, bread_phase_instance + integer :: num_jam, jam_phase_instance + integer, dimension(3) :: bread_phase_id, jam_phase_id + integer, allocatable :: jam_phase_id_correct(:), bread_phase_id_correct(:) ! create the camp core from test data allocate(file_list(1)) @@ -207,15 +209,39 @@ subroutine test_config_read() call assert(008760891, aero_rep%aero_phase_is_at_surface(12) .eqv. .true.) call assert(768528274, aero_rep%phase_state_size(layer=3,phase=1) .eq. 3) - ! test num_phase_instances funtion - phase_name_test = "bread" - num_bread = 2 + + ! test phase_ids and num_phase_instances function + ! test for jam (one instance in particle) + phase_name_test = "jam" + num_jam = 0 + jam_phase_id(1) = 2 + jam_phase_id(2) = 6 + jam_phase_id(3) = 10 max_part = aero_rep%maximum_computational_particles() - bread_phase_instance = num_bread * max_part + jam_phase_instance = 0 !check value call assert(417730478, 3 .eq. max_part) - call assert(734138496, bread_phase_instance .eq. aero_rep%num_phase_instances(phase_name_test)) - + jam_phase_id_correct = aero_rep%phase_ids(phase_name_test, is_at_surface = .false.) + call assert(757273007, jam_phase_id(1) .eq. jam_phase_id_correct(1)) + call assert(396819026, jam_phase_id(2) .eq. jam_phase_id_correct(2)) + call assert(777863671, jam_phase_id(3) .eq. jam_phase_id_correct(3)) + call assert(493602373, jam_phase_instance .eq. aero_rep%num_phase_instances(phase_name_test, & + is_at_surface = .true.)) + ! check for bread (two instances but only one at surface) + phase_name_test = "bread" + num_bread = 1 + bread_phase_id(1) = 4 + bread_phase_id(2) = 8 + bread_phase_id(3) = 12 + bread_phase_instance = num_bread * max_part + !check values + bread_phase_id_correct = aero_rep%phase_ids(phase_name_test, is_at_surface = .true.) + call assert(942495687, bread_phase_id(1) .eq. bread_phase_id_correct(1)) + call assert(283157050, bread_phase_id(2) .eq. bread_phase_id_correct(2)) + call assert(799763568, bread_phase_id(3) .eq. bread_phase_id_correct(3)) + call assert(734138496, bread_phase_instance .eq. aero_rep%num_phase_instances(phase_name_test, & + is_at_surface = .true.)) + class default call die_msg(519535557, rep_name) end select @@ -259,6 +285,19 @@ subroutine test_config_read() call assert(779745112, unique_names(35)%string .eq. "P3.top bread.bread.water") call assert(109646110, unique_names(36)%string .eq. "P3.top bread.bread.salt") + ! Test the unique names function with phase_is_at_surface flag + phase_name_test = "bread" + unique_names_surface = aero_rep%unique_names(phase_name=phase_name_test, phase_is_at_surface=.true.) + call assert(516157019, size( unique_names_surface ) .eq. 9) + call assert(071532611, unique_names_surface(1)%string .eq. "P1.top bread.bread.wheat") + call assert(911371605, unique_names_surface(2)%string .eq. "P1.top bread.bread.water") + call assert(537662837, unique_names_surface(3)%string .eq. "P1.top bread.bread.salt") + call assert(952202112, unique_names_surface(4)%string .eq. "P2.top bread.bread.wheat") + call assert(019187427, unique_names_surface(5)%string .eq. "P2.top bread.bread.water") + call assert(471502051, unique_names_surface(6)%string .eq. "P2.top bread.bread.salt") + call assert(623071633, unique_names_surface(7)%string .eq. "P3.top bread.bread.wheat") + call assert(862917237, unique_names_surface(8)%string .eq. "P3.top bread.bread.water") + call assert(521426951, unique_names_surface(9)%string .eq. "P3.top bread.bread.salt") #endif end subroutine test_config_read @@ -329,7 +368,7 @@ logical function build_aero_rep_data_set_test() end do end do - ! Set the species concentrations + ! Set the species concentrations (P1) phase_name = "top bread" spec_name = "wheat" unique_names = aero_rep%unique_names() @@ -384,6 +423,117 @@ logical function build_aero_rep_data_set_test() call assert_msg(149863598, i_spec.gt.0, rep_name) camp_state%state_var(i_spec) = 12.5 + ! Set the species concentrations (P2) + phase_name = "top bread" + spec_name = "wheat" + unique_names = aero_rep%unique_names() + i_spec = aero_rep%spec_state_id(unique_names(13)%string) + call assert_msg(752656273, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 1.5 + spec_name = "water" + i_spec = aero_rep%spec_state_id(unique_names(14)%string) + call assert_msg(491172859, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 2.5 + spec_name = "salt" + i_spec = aero_rep%spec_state_id(unique_names(15)%string) + call assert_msg(903979796, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 3.5 + phase_name = "jam" + spec_name = "rasberry" + i_spec = aero_rep%spec_state_id(unique_names(16)%string) + call assert_msg(094588550, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 4.5 + spec_name = "honey" + i_spec = aero_rep%spec_state_id(unique_names(17)%string) + call assert_msg(203341280, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 5.5 + spec_name = "sugar" + i_spec = aero_rep%spec_state_id(unique_names(18)%string) + call assert_msg(010697815, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 6.5 + spec_name = "lemon" + i_spec = aero_rep%spec_state_id(unique_names(19)%string) + call assert_msg(437481598, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 7.5 + phase_name = "almond butter" + spec_name = "almonds" + i_spec = aero_rep%spec_state_id(unique_names(20)%string) + call assert_msg(471657965, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 8.5 + spec_name = "sugar" + i_spec = aero_rep%spec_state_id(unique_names(21)%string) + call assert_msg(494073104, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 9.5 + phase_name = "bottom bread" + spec_name = "wheat" + i_spec = aero_rep%spec_state_id(unique_names(22)%string) + call assert_msg(390101312, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 10.5 + spec_name = "water" + i_spec = aero_rep%spec_state_id(unique_names(23)%string) + call assert_msg(320239819, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 11.5 + spec_name = "salt" + i_spec = aero_rep%spec_state_id(unique_names(24)%string) + call assert_msg(501204067, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 12.5 + + ! Set the species concentrations (P3) + phase_name = "top bread" + spec_name = "wheat" + unique_names = aero_rep%unique_names() + i_spec = aero_rep%spec_state_id(unique_names(25)%string) + call assert_msg(189910318, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 1.5 + spec_name = "water" + i_spec = aero_rep%spec_state_id(unique_names(26)%string) + call assert_msg(704399585, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 2.5 + spec_name = "salt" + i_spec = aero_rep%spec_state_id(unique_names(27)%string) + call assert_msg(418977803, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 3.5 + phase_name = "jam" + spec_name = "rasberry" + i_spec = aero_rep%spec_state_id(unique_names(28)%string) + call assert_msg(888754404, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 4.5 + spec_name = "honey" + i_spec = aero_rep%spec_state_id(unique_names(29)%string) + call assert_msg(006633311, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 5.5 + spec_name = "sugar" + i_spec = aero_rep%spec_state_id(unique_names(30)%string) + call assert_msg(701211482, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 6.5 + spec_name = "lemon" + i_spec = aero_rep%spec_state_id(unique_names(31)%string) + call assert_msg(818017683, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 7.5 + phase_name = "almond butter" + spec_name = "almonds" + i_spec = aero_rep%spec_state_id(unique_names(32)%string) + call assert_msg(735591255, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 8.5 + spec_name = "sugar" + i_spec = aero_rep%spec_state_id(unique_names(33)%string) + call assert_msg(467138849, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 9.5 + phase_name = "bottom bread" + spec_name = "wheat" + i_spec = aero_rep%spec_state_id(unique_names(34)%string) + call assert_msg(373107037, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 10.5 + spec_name = "water" + i_spec = aero_rep%spec_state_id(unique_names(35)%string) + call assert_msg(498732740, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 11.5 + spec_name = "salt" + i_spec = aero_rep%spec_state_id(unique_names(36)%string) + call assert_msg(165886615, i_spec.gt.0, rep_name) + camp_state%state_var(i_spec) = 12.5 + + end do rep_name = "AERO_REP_BAD_NAME" @@ -477,15 +627,18 @@ logical function eval_c_func(camp_core) result(passed) ! Check the number of Jacobian elements for the test phase call assert_msg(153137613, & - aero_rep%num_jac_elem(AERO_PHASE_IDX) .eq. NUM_JAC_ELEM, & + aero_rep%num_jac_elem(AERO_PHASE_IDX_1) .eq. NUM_JAC_ELEM, & + rep_name) + call assert_msg(994566346, & + aero_rep%num_jac_elem(AERO_PHASE_IDX_2) .eq. NUM_JAC_ELEM, & rep_name) ! Update external properties - call update_number%set_number__n_m3( TEST_PARTICLE, 12.3d0 ) + call update_number%set_number__n_m3( TEST_PARTICLE_1, 12.3d0 ) call camp_core%update_data( update_number ) ! Test re-setting number concentration - call update_number%set_number__n_m3( TEST_PARTICLE, PART_NUM_CONC ) + call update_number%set_number__n_m3( TEST_PARTICLE_1, PART_NUM_CONC ) call camp_core%update_data( update_number ) passed = run_aero_rep_single_particle_c_tests( & diff --git a/test/unit_aero_rep_data/test_aero_rep_single_particle.c b/test/unit_aero_rep_data/test_aero_rep_single_particle.c index 307567c90..6a1916ef1 100644 --- a/test/unit_aero_rep_data/test_aero_rep_single_particle.c +++ b/test/unit_aero_rep_data/test_aero_rep_single_particle.c @@ -17,7 +17,7 @@ #define AERO_REP_IDX 0 // test computational particle -#define TEST_PARTICLE 2 +#define TEST_PARTICLE_1 2 // number of computational particles in the test #define N_COMP_PARTICLES 3 @@ -25,8 +25,11 @@ // number of aerosol phases per particle #define NUM_AERO_PHASE 4 -// index for the test phase (test-particle phase 2) -#define AERO_PHASE_IDX ((TEST_PARTICLE-1)*NUM_AERO_PHASE+1) +// index for the test phase +// (test-particle phase 2 : middle layer, jam) +#define AERO_PHASE_IDX_1 ((TEST_PARTICLE_1-1)*NUM_AERO_PHASE+1) +// (test-particle phase 4 : top layer, bread) +#define AERO_PHASE_IDX_2 ((TEST_PARTICLE_1-1)*NUM_AERO_PHASE+3) // number of Jacobian elements used for the test phase #define N_JAC_ELEM 12 @@ -40,10 +43,6 @@ #define CONC_sugar 6.0 #define CONC_lemon 7.0 #define CONC_almonds 8.0 -#define CONC_sugarB 9.0 -#define CONC_wheatB 10.0 -#define CONC_waterB 11.0 -#define CONC_saltB 12.0 // Molecular weight (kg/mol) of test species (must match json file) #define MW_wheat 1.0 @@ -79,13 +78,17 @@ int test_effective_radius(ModelData * model_data, N_Vector state) { int ret_val = 0; - double partial_deriv[N_JAC_ELEM+2]; - double eff_rad = -999.9; + double partial_deriv_1[N_JAC_ELEM+2]; + double partial_deriv_2[N_JAC_ELEM+2]; + double eff_rad_1 = -999.9; + double eff_rad_2 = -999.9; - for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv[i] = 999.9; + for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv_1[i] = 999.9; aero_rep_get_effective_radius__m(model_data, AERO_REP_IDX, - AERO_PHASE_IDX, &eff_rad, &(partial_deriv[1])); + AERO_PHASE_IDX_1, &eff_rad_1, &(partial_deriv_1[1])); + aero_rep_get_effective_radius__m(model_data, AERO_REP_IDX, + AERO_PHASE_IDX_2, &eff_rad_2, &(partial_deriv_2[1])); double volume_density = ( CONC_wheat / DENSITY_wheat + CONC_water / DENSITY_water + @@ -95,45 +98,47 @@ int test_effective_radius(ModelData * model_data, N_Vector state) { CONC_sugar / DENSITY_sugar + CONC_lemon / DENSITY_lemon + CONC_almonds / DENSITY_almonds + - CONC_sugarB / DENSITY_sugar + - CONC_wheatB / DENSITY_wheat + - CONC_waterB / DENSITY_water + - CONC_saltB / DENSITY_salt ); // volume density (m3/m3) + CONC_sugar / DENSITY_sugar + + CONC_wheat / DENSITY_wheat + + CONC_water / DENSITY_water + + CONC_salt / DENSITY_salt ); // volume density (m3/m3) double eff_rad_expected = pow( ( 3.0 / 4.0 / 3.14159265359 * volume_density ), 1.0/3.0 ); - ret_val += ASSERT_MSG(fabs(eff_rad-eff_rad_expected) < 1.0e-6*eff_rad_expected, + ret_val += ASSERT_MSG(fabs(eff_rad_1-eff_rad_expected) < 1.0e-6*eff_rad_expected, + "Bad effective radius"); + ret_val += ASSERT_MSG(fabs(eff_rad_2-eff_rad_expected) < 1.0e-6*eff_rad_expected, "Bad effective radius"); - ret_val += ASSERT_MSG(partial_deriv[0] = 999.9, + ret_val += ASSERT_MSG(partial_deriv_1[0] = 999.9, "Bad Jacobian (-1)"); double d_eff_rad_dx = 1.0 / 4.0 / 3.14159265359 * pow( 3.0 / 4.0 / 3.14159265359 * volume_density, -2.0/3.0 ); - ret_val += ASSERT_MSG(fabs(partial_deriv[1] - d_eff_rad_dx / DENSITY_wheat) < - 1.0e-10 * partial_deriv[1], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[2] - d_eff_rad_dx / DENSITY_water) < - 1.0e-10 * partial_deriv[2], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[3] - d_eff_rad_dx / DENSITY_salt) < - 1.0e-10 * partial_deriv[3], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[4] - d_eff_rad_dx / DENSITY_rasberry) < - 1.0e-10 * partial_deriv[4], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[5] - d_eff_rad_dx / DENSITY_honey) < - 1.0e-10 * partial_deriv[5], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[6] - d_eff_rad_dx / DENSITY_sugar) < - 1.0e-10 * partial_deriv[6], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[7] - d_eff_rad_dx / DENSITY_lemon) < - 1.0e-10 * partial_deriv[7], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[8] - d_eff_rad_dx / DENSITY_almonds) < - 1.0e-10 * partial_deriv[8], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[9] - d_eff_rad_dx / DENSITY_sugar) < - 1.0e-10 * partial_deriv[9], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[10] - d_eff_rad_dx / DENSITY_wheat) < - 1.0e-10 * partial_deriv[10], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[11] - d_eff_rad_dx / DENSITY_water) < - 1.0e-10 * partial_deriv[11], "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[12] - d_eff_rad_dx / DENSITY_salt) < - 1.0e-10 * partial_deriv[12], "Bad Jacobian element"); - ret_val += ASSERT_MSG(partial_deriv[N_JAC_ELEM+1] = 999.9, + ret_val += ASSERT_MSG(fabs(partial_deriv_1[1] - d_eff_rad_dx / DENSITY_wheat) < + 1.0e-10 * partial_deriv_1[1], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[2] - d_eff_rad_dx / DENSITY_water) < + 1.0e-10 * partial_deriv_1[2], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[3] - d_eff_rad_dx / DENSITY_salt) < + 1.0e-10 * partial_deriv_1[3], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[4] - d_eff_rad_dx / DENSITY_rasberry) < + 1.0e-10 * partial_deriv_1[4], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[5] - d_eff_rad_dx / DENSITY_honey) < + 1.0e-10 * partial_deriv_1[5], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[6] - d_eff_rad_dx / DENSITY_sugar) < + 1.0e-10 * partial_deriv_1[6], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[7] - d_eff_rad_dx / DENSITY_lemon) < + 1.0e-10 * partial_deriv_1[7], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[8] - d_eff_rad_dx / DENSITY_almonds) < + 1.0e-10 * partial_deriv_1[8], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[9] - d_eff_rad_dx / DENSITY_sugar) < + 1.0e-10 * partial_deriv_1[9], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[10] - d_eff_rad_dx / DENSITY_wheat) < + 1.0e-10 * partial_deriv_1[10], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[11] - d_eff_rad_dx / DENSITY_water) < + 1.0e-10 * partial_deriv_1[11], "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_1[12] - d_eff_rad_dx / DENSITY_salt) < + 1.0e-10 * partial_deriv_1[12], "Bad Jacobian element"); + ret_val += ASSERT_MSG(partial_deriv_1[N_JAC_ELEM+1] = 999.9, "Bad Jacobian (end+1)"); return ret_val; @@ -154,8 +159,9 @@ int test_number_concentration(ModelData * model_data, N_Vector state) { for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv[i] = 999.9; aero_rep_get_number_conc__n_m3(model_data, AERO_REP_IDX, - AERO_PHASE_IDX, &num_conc, &(partial_deriv[1])); + AERO_PHASE_IDX_1, &num_conc, &(partial_deriv[1])); + printf("\nnumber conc : %f", num_conc); ret_val += ASSERT_MSG(fabs(num_conc-PART_NUM_CONC) < 1.0e-10*PART_NUM_CONC, "Bad number concentration"); @@ -179,38 +185,65 @@ int test_number_concentration(ModelData * model_data, N_Vector state) { int test_aero_phase_mass(ModelData * model_data, N_Vector state) { int ret_val = 0; - double partial_deriv[N_JAC_ELEM+2]; - double phase_mass = -999.9; + double partial_deriv_1[N_JAC_ELEM+2]; + double partial_deriv_2[N_JAC_ELEM+2]; + double phase_mass_1 = -999.9; + double phase_mass_2 = -999.9; - for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv[i] = 999.9; + for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv_1[i] = 999.9; - aero_rep_get_aero_phase_mass__kg_m3(model_data, AERO_REP_IDX, AERO_PHASE_IDX, - &phase_mass, &(partial_deriv[1])); + for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv_2[i] = 999.9; - double mass = CONC_rasberry + CONC_honey + CONC_sugar + CONC_lemon; + aero_rep_get_aero_phase_mass__kg_m3(model_data, AERO_REP_IDX, AERO_PHASE_IDX_1, + &phase_mass_1, &(partial_deriv_1[1])); - ret_val += ASSERT_MSG(fabs(phase_mass-mass) < 1.0e-10*mass, + aero_rep_get_aero_phase_mass__kg_m3(model_data, AERO_REP_IDX, AERO_PHASE_IDX_2, + &phase_mass_2, &(partial_deriv_2[1])); + + double mass_1 = CONC_rasberry + CONC_honey + CONC_sugar + CONC_lemon; + double mass_2 = CONC_wheat + CONC_water + CONC_salt; + + ret_val += ASSERT_MSG(fabs(phase_mass_1-mass_1) < 1.0e-10*mass_1, "Bad aerosol phase mass"); - ret_val += ASSERT_MSG(partial_deriv[0] = 999.9, + ret_val += ASSERT_MSG(partial_deriv_1[0] = 999.9, "Bad Jacobian (-1)"); for( int i = 1; i < 4; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ZERO, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ZERO, "Bad Jacobian element"); for( int i = 4; i < 8; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ONE, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ONE, "Bad Jacobian element"); for( int i = 8; i < 10; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ZERO, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ZERO, "Bad Jacobian element"); for( int i = 10; i < N_JAC_ELEM+1; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ZERO, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ZERO, "Bad Jacobian element"); - ret_val += ASSERT_MSG(partial_deriv[N_JAC_ELEM+1] = 999.9, + ret_val += ASSERT_MSG(partial_deriv_1[N_JAC_ELEM+1] = 999.9, "Bad Jacobian (end+1)"); - return ret_val; + ret_val += ASSERT_MSG(fabs(phase_mass_2-mass_2) < 1.0e-10*mass_2, + "Bad aerosol phase mass"); + printf("\nphase_mass_2 : %f", phase_mass_2); + ret_val += ASSERT_MSG(partial_deriv_2[0] = 999.9, + "Bad Jacobian (-1)"); + for( int i = 1; i < 4; ++i ) + ret_val += ASSERT_MSG(partial_deriv_2[i] == ZERO, + "Bad Jacobian element"); + for( int i = 4; i < 8; ++i ) + ret_val += ASSERT_MSG(partial_deriv_2[i] == ZERO, + "Bad Jacobian element"); + for( int i = 8; i < 10; ++i ) + ret_val += ASSERT_MSG(partial_deriv_2[i] == ZERO, + "Bad Jacobian element"); + for( int i = 10; i < N_JAC_ELEM+1; ++i ) + ret_val += ASSERT_MSG(partial_deriv_2[i] == ONE, + "Bad Jacobian element"); + ret_val += ASSERT_MSG(partial_deriv_2[N_JAC_ELEM+1] = 999.9, + "Bad Jacobian (end+1)"); + return ret_val; } /** \brief Test the aerosol phase average molecular weight function @@ -222,44 +255,73 @@ int test_aero_phase_mass(ModelData * model_data, N_Vector state) { int test_aero_phase_avg_MW(ModelData * model_data, N_Vector state) { int ret_val = 0; - double partial_deriv[N_JAC_ELEM+2]; - double avg_mw = -999.9; - - for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv[i] = 999.9; - - aero_rep_get_aero_phase_avg_MW__kg_mol(model_data, AERO_REP_IDX, AERO_PHASE_IDX, - &avg_mw, &(partial_deriv[1])); - - double mass = CONC_rasberry + CONC_honey + CONC_sugar + CONC_lemon; - double moles = CONC_rasberry / MW_rasberry + CONC_honey / MW_honey + CONC_sugar / MW_sugar + CONC_lemon / MW_lemon; - double avg_mw_real = mass / moles; - double dMW_drasberry = ONE / moles - mass / (moles * moles * MW_rasberry); - double dMW_dhoney = ONE / moles - mass / (moles * moles * MW_honey); - double dMW_dsugar = ONE / moles - mass / (moles * moles * MW_sugar); - double dMW_dlemon = ONE / moles - mass / (moles * moles * MW_lemon); - - ret_val += ASSERT_MSG(fabs(avg_mw-avg_mw_real) < 1.0e-10*avg_mw_real, + double partial_deriv_1[N_JAC_ELEM+2]; + double partial_deriv_2[N_JAC_ELEM+2]; + double avg_mw_1 = -999.9; + double avg_mw_2 = -999.9; + + for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv_1[i] = 999.9; + for( int i = 0; i < N_JAC_ELEM+2; ++i ) partial_deriv_2[i] = 999.9; + + aero_rep_get_aero_phase_avg_MW__kg_mol(model_data, AERO_REP_IDX, AERO_PHASE_IDX_1, + &avg_mw_1, &(partial_deriv_1[1])); + aero_rep_get_aero_phase_avg_MW__kg_mol(model_data, AERO_REP_IDX, AERO_PHASE_IDX_2, + &avg_mw_2, &(partial_deriv_2[1])); + + double mass_1 = CONC_rasberry + CONC_honey + CONC_sugar + CONC_lemon; + double moles_1 = CONC_rasberry / MW_rasberry + CONC_honey / MW_honey + CONC_sugar / MW_sugar + CONC_lemon / MW_lemon; + double avg_mw_real_1 = mass_1 / moles_1; + double dMW_drasberry = ONE / moles_1 - mass_1 / (moles_1 * moles_1 * MW_rasberry); + double dMW_dhoney = ONE / moles_1 - mass_1 / (moles_1 * moles_1 * MW_honey); + double dMW_dsugar = ONE / moles_1 - mass_1 / (moles_1 * moles_1 * MW_sugar); + double dMW_dlemon = ONE / moles_1 - mass_1 / (moles_1 * moles_1 * MW_lemon); + + double mass_2 = CONC_wheat + CONC_water + CONC_salt; + double moles_2 = CONC_wheat / MW_wheat + CONC_water / MW_water + CONC_salt / MW_salt; + double avg_mw_real_2 = mass_2 / moles_2; + double dMW_dwheat = ONE / moles_2 - mass_2 / (moles_2 * moles_2 * MW_wheat); + double dMW_dwater = ONE / moles_2 - mass_2 / (moles_2 * moles_2 * MW_water); + double dMW_dsalt = ONE / moles_2 - mass_2 / (moles_2 * moles_2 * MW_salt); + + ret_val += ASSERT_MSG(fabs(avg_mw_1-avg_mw_real_1) < 1.0e-10*avg_mw_real_1, "Bad average MW"); - ret_val += ASSERT_MSG(partial_deriv[0] = 999.9, + ret_val += ASSERT_MSG(partial_deriv_1[0] = 999.9, "Bad Jacobian (-1)"); for( int i = 1; i < 4; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ZERO, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ZERO, "Bad Jacobian element"); - ret_val += ASSERT_MSG(fabs(partial_deriv[4]-dMW_drasberry) < 1.0e-10*fabs(dMW_drasberry), + ret_val += ASSERT_MSG(fabs(partial_deriv_1[4]-dMW_drasberry) < 1.0e-10*fabs(dMW_drasberry), "Bad Jacobian (-1)"); - ret_val += ASSERT_MSG(fabs(partial_deriv[5]-dMW_dhoney) < 1.0e-10*fabs(dMW_dhoney), + ret_val += ASSERT_MSG(fabs(partial_deriv_1[5]-dMW_dhoney) < 1.0e-10*fabs(dMW_dhoney), "Bad Jacobian (-1)"); - ret_val += ASSERT_MSG(fabs(partial_deriv[6]-dMW_dsugar) < 1.0e-10*fabs(dMW_dsugar), + ret_val += ASSERT_MSG(fabs(partial_deriv_1[6]-dMW_dsugar) < 1.0e-10*fabs(dMW_dsugar), "Bad Jacobian (-1)"); - ret_val += ASSERT_MSG(fabs(partial_deriv[7]-dMW_dlemon) < 1.0e-10*fabs(dMW_dlemon), + ret_val += ASSERT_MSG(fabs(partial_deriv_1[7]-dMW_dlemon) < 1.0e-10*fabs(dMW_dlemon), "Bad Jacobian (-1)"); for( int i = 8; i < N_JAC_ELEM+1; ++i ) - ret_val += ASSERT_MSG(partial_deriv[i] == ZERO, + ret_val += ASSERT_MSG(partial_deriv_1[i] == ZERO, "Bad Jacobian element"); - ret_val += ASSERT_MSG(partial_deriv[N_JAC_ELEM+1] = 999.9, + ret_val += ASSERT_MSG(partial_deriv_1[N_JAC_ELEM+1] = 999.9, "Bad Jacobian (end+1)"); + + ret_val += ASSERT_MSG(fabs(avg_mw_2-avg_mw_real_2) < 1.0e-10*avg_mw_real_2, + "Bad average MW"); + + ret_val += ASSERT_MSG(partial_deriv_2[0] = 999.9, + "Bad Jacobian (-1)"); + for( int i = 1; i < 9; ++i ) + ret_val += ASSERT_MSG(partial_deriv_2[i] == ZERO, + "Bad Jacobian element"); + ret_val += ASSERT_MSG(fabs(partial_deriv_2[10]-dMW_dwheat) < 1.0e-10*fabs(dMW_dwheat), + "Bad Jacobian (-1)"); + ret_val += ASSERT_MSG(fabs(partial_deriv_2[11]-dMW_dwater) < 1.0e-10*fabs(dMW_dwater), + "Bad Jacobian (-1)"); + ret_val += ASSERT_MSG(fabs(partial_deriv_2[12]-dMW_dsalt) < 1.0e-10*fabs(dMW_dsalt), + "Bad Jacobian (-1)"); + ret_val += ASSERT_MSG(partial_deriv_2[N_JAC_ELEM+1] = 999.9, + "Bad Jacobian (end+1)"); return ret_val; } @@ -295,7 +357,7 @@ int run_aero_rep_single_particle_c_tests(void *solver_data, double *state, doubl for (int i_var=0; i_vargrid_cell_aero_rep_env_data = model_data->aero_rep_env_data; @@ -327,6 +389,8 @@ int run_aero_rep_single_particle_c_tests(void *solver_data, double *state, doubl aero_rep_update_env_state(model_data); aero_rep_update_state(model_data); + aero_rep_print_data(sd); + // Run the property tests ret_val += test_effective_radius(model_data, solver_state); ret_val += test_aero_phase_mass(model_data, solver_state); diff --git a/test/unit_rxn_data/test_HL_phase_transfer.json b/test/unit_rxn_data/test_HL_phase_transfer.json index 155efe6e6..168367058 100644 --- a/test/unit_rxn_data/test_HL_phase_transfer.json +++ b/test/unit_rxn_data/test_HL_phase_transfer.json @@ -79,6 +79,13 @@ "phases": [ "aqueous aerosol" ], + "covers": "two layer" + }, + { + "name": "two layer", + "phases": [ + "aqueous aerosol" + ], "covers": "none" } ] diff --git a/test/unit_rxn_data/test_SIMPOL_phase_transfer.json b/test/unit_rxn_data/test_SIMPOL_phase_transfer.json index 2586fab40..429a6e9d3 100644 --- a/test/unit_rxn_data/test_SIMPOL_phase_transfer.json +++ b/test/unit_rxn_data/test_SIMPOL_phase_transfer.json @@ -40,7 +40,7 @@ "maximum computational particles" : 1, "layers": [ { - "name": "one layer", + "name": "one_layer", "phases": [ "aqueous aerosol" ], @@ -58,6 +58,13 @@ "phases": [ "aqueous aerosol" ], + "covers": "two layer" + }, + { + "name": "two layer", + "phases": [ + "aqueous aerosol" + ], "covers": "none" } ] @@ -68,7 +75,7 @@ "maximum computational particles" : 1, "layers": [ { - "name": "one layer", + "name": "one_layer", "phases": [ "aqueous aerosol" ], diff --git a/test/unit_rxn_data/test_rxn_HL_phase_transfer.F90 b/test/unit_rxn_data/test_rxn_HL_phase_transfer.F90 index 4b99de9f2..1d828f527 100644 --- a/test/unit_rxn_data/test_rxn_HL_phase_transfer.F90 +++ b/test/unit_rxn_data/test_rxn_HL_phase_transfer.F90 @@ -97,10 +97,13 @@ logical function run_HL_phase_transfer_test(scenario) type(chem_spec_data_t), pointer :: chem_spec_data class(aero_rep_data_t), pointer :: aero_rep_ptr real(kind=dp), allocatable, dimension(:,:) :: model_conc, true_conc - integer(kind=i_kind) :: idx_phase, idx_aero_rep, idx_O3, idx_O3_aq, & - idx_H2O2, idx_H2O2_aq, idx_H2O_aq, idx_HNO3, idx_HNO3_aq, & + integer(kind=i_kind) :: idx_phase, idx_aero_rep, idx_HNO3, idx_HNO3_aq, & idx_NH3, idx_NH3_aq, idx_H2O, idx_Na_p, idx_Cl_m, idx_Ca_pp, & i_time, i_spec + integer(kind=i_kind) :: idx_O3, idx_O3_aq, & + idx_O3_aq_layer1, idx_O3_aq_layer2, idx_H2O2, idx_H2O2_aq, & + idx_H2O2_aq_layer1, idx_H2O2_aq_layer2, idx_H2O_aq, & + idx_H2O_aq_layer1, idx_H2O_aq_layer2 character(len=:), allocatable :: key, idx_prefix real(kind=dp) :: time_step, time, n_star, del_H, del_S, del_G, alpha, & mfp, Kn, corr, M_to_ppm, kgm3_to_ppm, K_eq_O3, K_eq_H2O2, & @@ -127,6 +130,7 @@ logical function run_HL_phase_transfer_test(scenario) type(aero_rep_update_data_modal_binned_mass_GMD_t) :: update_data_GMD type(aero_rep_update_data_modal_binned_mass_GSD_t) :: update_data_GSD + print *, "did we make it here?" call assert_msg(144071521, scenario.ge.1 .and. scenario.le.2, & "Invalid scenario specified: "//to_string( scenario ) ) @@ -134,8 +138,8 @@ logical function run_HL_phase_transfer_test(scenario) ! Allocate space for the results if (scenario.eq.1) then - allocate(model_conc(0:NUM_TIME_STEP, 11)) - allocate(true_conc(0:NUM_TIME_STEP, 11)) + allocate(model_conc(0:NUM_TIME_STEP, 14)) + allocate(true_conc(0:NUM_TIME_STEP, 14)) else if (scenario.eq.2) then allocate(model_conc(0:NUM_TIME_STEP, 24)) allocate(true_conc(0:NUM_TIME_STEP, 24)) @@ -210,20 +214,44 @@ logical function run_HL_phase_transfer_test(scenario) ! Get species indices if (scenario.eq.1) then - idx_prefix = "P1.one layer." + idx_prefix = "P1." + key = "O3" + idx_O3 = chem_spec_data%gas_state_id(key); + key = idx_prefix//"one layer.aqueous aerosol.O3_aq" + print *, "key", key + idx_O3_aq_layer1 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"two layer.aqueous aerosol.O3_aq" + idx_O3_aq_layer2 = aero_rep_ptr%spec_state_id(key); + key = "H2O2" + idx_H2O2 = chem_spec_data%gas_state_id(key); + key = idx_prefix//"one layer.aqueous aerosol.H2O2_aq" + idx_H2O2_aq_layer1 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"two layer.aqueous aerosol.H2O2_aq" + idx_H2O2_aq_layer2 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"one layer.aqueous aerosol.H2O_aq" + idx_H2O_aq_layer1 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"two layer.aqueous aerosol.H2O_aq" + idx_H2O_aq_layer2 = aero_rep_ptr%spec_state_id(key); + ! Make sure the expected species are in the model + call assert(167389209, idx_O3.gt.0) + call assert(156210838, idx_O3_aq_layer1.gt.0) + call assert(802672616, idx_O3_aq_layer2.gt.0) + call assert(844092221, idx_H2O2.gt.0) + call assert(097661255, idx_H2O2_aq_layer1.gt.0) + call assert(429463219, idx_H2O2_aq_layer2.gt.0) + call assert(312315501, idx_H2O_aq_layer1.gt.0) + call assert(448976280, idx_H2O_aq_layer2.gt.0) else if (scenario.eq.2) then idx_prefix = "the mode." - end if + key = "O3" + idx_O3 = chem_spec_data%gas_state_id(key); + key = idx_prefix//"aqueous aerosol.O3_aq" + idx_O3_aq = aero_rep_ptr%spec_state_id(key); + key = "H2O2" + idx_H2O2 = chem_spec_data%gas_state_id(key); + key = idx_prefix//"aqueous aerosol.H2O2_aq" + idx_H2O2_aq = aero_rep_ptr%spec_state_id(key); - key = "O3" - idx_O3 = chem_spec_data%gas_state_id(key); - key = idx_prefix//"aqueous aerosol.O3_aq" - idx_O3_aq = aero_rep_ptr%spec_state_id(key); - key = "H2O2" - idx_H2O2 = chem_spec_data%gas_state_id(key); - key = idx_prefix//"aqueous aerosol.H2O2_aq" - idx_H2O2_aq = aero_rep_ptr%spec_state_id(key); - if (scenario.eq.2) then key = "HNO3" idx_HNO3 = chem_spec_data%gas_state_id(key); key = idx_prefix//"aqueous aerosol.HNO3_aq" @@ -240,23 +268,20 @@ logical function run_HL_phase_transfer_test(scenario) idx_Cl_m = aero_rep_ptr%spec_state_id(key); key = idx_prefix//"aqueous aerosol.Ca_pp" idx_Ca_pp = aero_rep_ptr%spec_state_id(key); - end if - key = idx_prefix//"aqueous aerosol.H2O_aq" - idx_H2O_aq = aero_rep_ptr%spec_state_id(key); - - ! Make sure the expected species are in the model - call assert(202593939, idx_O3.gt.0) - call assert(262338032, idx_O3_aq.gt.0) - call assert(374656377, idx_H2O2.gt.0) - call assert(704441571, idx_H2O2_aq.gt.0) - if (scenario.eq.2) then + key = idx_prefix//"aqueous aerosol.H2O_aq" + idx_H2O_aq = aero_rep_ptr%spec_state_id(key); + ! Make sure the expected species are in the model + call assert(202593939, idx_O3.gt.0) + call assert(262338032, idx_O3_aq.gt.0) + call assert(374656377, idx_H2O2.gt.0) + call assert(704441571, idx_H2O2_aq.gt.0) call assert(274044966, idx_HNO3.gt.0) call assert(386363311, idx_HNO3_aq.gt.0) call assert(833731157, idx_NH3.gt.0) call assert(946049502, idx_NH3_aq.gt.0) call assert(688163399, idx_H2O.gt.0) + call assert(758921357, idx_H2O_aq.gt.0) end if - call assert(758921357, idx_H2O_aq.gt.0) #ifdef CAMP_USE_MPI ! pack the camp core @@ -283,10 +308,18 @@ logical function run_HL_phase_transfer_test(scenario) ! broadcast the species ids call camp_mpi_bcast_integer(idx_O3) - call camp_mpi_bcast_integer(idx_O3_aq) call camp_mpi_bcast_integer(idx_H2O2) - call camp_mpi_bcast_integer(idx_H2O2_aq) - if (scenario.eq.2) then + if (scenario.eq.1) then + call camp_mpi_bcast_integer(idx_O3_aq_layer1) + call camp_mpi_bcast_integer(idx_O3_aq_layer2) + call camp_mpi_bcast_integer(idx_H2O2_aq_layer1) + call camp_mpi_bcast_integer(idx_H2O2_aq_layer2) + call camp_mpi_bcast_integer(idx_H2O_aq_layer1) + call camp_mpi_bcast_integer(idx_H2O_aq_layer2) + else if (scenario.eq.2) then + call camp_mpi_bcast_integer(idx_O3_aq) + call camp_mpi_bcast_integer(idx_H2O2_aq) + call camp_mpi_bcast_integer(idx_H2O_aq) call camp_mpi_bcast_integer(idx_HNO3) call camp_mpi_bcast_integer(idx_HNO3_aq) call camp_mpi_bcast_integer(idx_NH3) @@ -296,7 +329,6 @@ logical function run_HL_phase_transfer_test(scenario) call camp_mpi_bcast_integer(idx_Cl_m) call camp_mpi_bcast_integer(idx_Ca_pp) end if - call camp_mpi_bcast_integer(idx_H2O_aq) call camp_mpi_bcast_integer(i_sect_unused) call camp_mpi_bcast_integer(i_sect_the_mode) @@ -355,10 +387,13 @@ logical function run_HL_phase_transfer_test(scenario) true_conc(:,:) = 0.0 if (scenario.eq.1) then true_conc(0,idx_O3) = 0.0 - true_conc(0,idx_O3_aq) = 1.0e-3 + true_conc(0,idx_O3_aq_layer1) = 1.0e-3 + true_conc(0,idx_O3_aq_layer2) = 1.0e-3 true_conc(0,idx_H2O2) = 1.0 - true_conc(0,idx_H2O2_aq) = 0.0 - true_conc(0,idx_H2O_aq) = 1.4e-2 + true_conc(0,idx_H2O2_aq_layer1) = 0.0 + true_conc(0,idx_H2O2_aq_layer2) = 0.0 + true_conc(0,idx_H2O_aq_layer1) = 1.4e-2 + true_conc(0,idx_H2O_aq_layer2) = 1.4e-2 else if (scenario.eq.2) then true_conc(0,idx_O3) = 0.0 true_conc(0,idx_O3_aq) = 1.0 @@ -380,11 +415,15 @@ logical function run_HL_phase_transfer_test(scenario) if (scenario.eq.1) then number_conc = 1.3e6 ! particle number concentration (#/cc) ! single particle aerosol mass concentrations are per particle - true_conc(0,idx_O3_aq) = true_conc(0,idx_O3_aq) / number_conc - true_conc(0,idx_H2O_aq) = true_conc(0,idx_H2O_aq) / number_conc + true_conc(0,idx_O3_aq_layer1) = true_conc(0,idx_O3_aq_layer1) / number_conc + true_conc(0,idx_O3_aq_layer2) = true_conc(0,idx_O3_aq_layer2) / number_conc + true_conc(0,idx_H2O_aq_layer1) = true_conc(0,idx_H2O_aq_layer1) / number_conc + true_conc(0,idx_H2O_aq_layer2) = true_conc(0,idx_H2O_aq_layer2) / number_conc ! radius (m) calculated based on particle mass - radius = ( ( true_conc(0,idx_O3_aq) / 1000.0 + & - true_conc(0,idx_H2O_aq) / 1000.0 ) & + radius = ( ( true_conc(0,idx_O3_aq_layer1) / 1000.0 + & + true_conc(0,idx_O3_aq_layer2) / 1000.0 + & + true_conc(0,idx_H2O_aq_layer1) / 1000.0 + & + true_conc(0,idx_H2O_aq_layer2) / 1000.0 ) & * 3.0 / 4.0 / 3.14159265359 )**(1.0/3.0) else if (scenario.eq.2) then ! radius (m) @@ -418,8 +457,13 @@ logical function run_HL_phase_transfer_test(scenario) end if ! Determine the M -> ppm conversion using the total aerosol water - M_to_ppm = number_conc * 1.0d6 * true_conc(0,idx_H2O_aq) * & - const%univ_gas_const * temp / pressure + if (scenario.eq.1) then + M_to_ppm = number_conc * 1.0d6 * true_conc(0,idx_H2O_aq_layer1) * & + const%univ_gas_const * temp / pressure + else if (scenario.eq.2) then + M_to_ppm = number_conc * 1.0d6 * true_conc(0,idx_H2O_aq) * & + const%univ_gas_const * temp / pressure + end if ! O3 rate constants n_star = 1.89d0 @@ -464,21 +508,39 @@ logical function run_HL_phase_transfer_test(scenario) ! [A_aero] = [A_total] / (K_HL + 1) kgm3_to_ppm = const%univ_gas_const * 1.0e6 * temp & / (MW_O3 * pressure) - equil_O3 = (true_conc(0,idx_O3) + & - true_conc(0,idx_O3_aq)*number_conc*kgm3_to_ppm) / & - (K_eq_O3*M_to_ppm + 1.0d0) - equil_O3_aq = (true_conc(0,idx_O3)/kgm3_to_ppm/number_conc + & - true_conc(0,idx_O3_aq)) / & - (1.0d0 + 1.0d0/(K_eq_O3*M_to_ppm)) + if (scenario.eq.1) then + equil_O3 = (true_conc(0,idx_O3) + & + (true_conc(0,idx_O3_aq_layer1)+true_conc(0,idx_O3_aq_layer1)) * & + number_conc*kgm3_to_ppm) / (K_eq_O3*M_to_ppm + 1.0d0) + equil_O3_aq = (true_conc(0,idx_O3)/kgm3_to_ppm/number_conc + & + true_conc(0,idx_O3_aq_layer1) + true_conc(0,idx_O3_aq_layer2)) / & + (1.0d0 + 1.0d0/(K_eq_O3*M_to_ppm)) + else if (scenario.eq.2) then + equil_O3 = (true_conc(0,idx_O3) + & + true_conc(0,idx_O3_aq)*number_conc*kgm3_to_ppm) / & + (K_eq_O3*M_to_ppm + 1.0d0) + equil_O3_aq = (true_conc(0,idx_O3)/kgm3_to_ppm/number_conc + & + true_conc(0,idx_O3_aq)) / & + (1.0d0 + 1.0d0/(K_eq_O3*M_to_ppm)) + end if kgm3_to_ppm = const%univ_gas_const * 1.0e6 * temp & / (MW_H2O2 * pressure) - equil_H2O2 = (true_conc(0,idx_H2O2) + & - true_conc(0,idx_H2O2_aq)*number_conc*kgm3_to_ppm) / & - (K_eq_H2O2*M_to_ppm + 1.0d0) - equil_H2O2_aq = (true_conc(0,idx_H2O2)/kgm3_to_ppm/number_conc + & - true_conc(0,idx_H2O2_aq)) / & - (1.0d0 + 1.0d0/(K_eq_H2O2*M_to_ppm)) + if (scenario.eq.1) then + equil_H2O2 = (true_conc(0,idx_H2O2) + & + (true_conc(0,idx_H2O2_aq_layer1) + true_conc(0,idx_H2O2_aq_layer2)) * & + number_conc*kgm3_to_ppm) / (K_eq_H2O2*M_to_ppm + 1.0d0) + equil_H2O2_aq = (true_conc(0,idx_H2O2)/kgm3_to_ppm/number_conc + & + true_conc(0,idx_H2O2_aq_layer1) + true_conc(0,idx_H2O2_aq_layer2)) / & + (1.0d0 + 1.0d0/(K_eq_H2O2*M_to_ppm)) + else if (scenario.eq.2) then + equil_H2O2 = (true_conc(0,idx_H2O2) + & + true_conc(0,idx_H2O2_aq)*number_conc*kgm3_to_ppm) / & + (K_eq_H2O2*M_to_ppm + 1.0d0) + equil_H2O2_aq = (true_conc(0,idx_H2O2)/kgm3_to_ppm/number_conc + & + true_conc(0,idx_H2O2_aq)) / & + (1.0d0 + 1.0d0/(K_eq_H2O2*M_to_ppm)) + end if ! Set the initial state in the model camp_state%state_var(:) = model_conc(0,:) @@ -514,16 +576,30 @@ logical function run_HL_phase_transfer_test(scenario) ! [A_aero] = ([A_init_aero] - [A_eq_aero]) ! * exp(-t * (k_f + k_b)) + [A_eq_aero] time = i_time * time_step - true_conc(i_time,idx_O3) = (true_conc(0,idx_O3) - equil_O3) * & - exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3 - true_conc(i_time,idx_O3_aq) = (true_conc(0,idx_O3_aq) - equil_O3_aq) * & - exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3_aq - true_conc(i_time,idx_H2O2) = (true_conc(0,idx_H2O2) - equil_H2O2) * & - exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2 - true_conc(i_time,idx_H2O2_aq) = & - (true_conc(0,idx_H2O2_aq) - equil_H2O2_aq) * & - exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2_aq - true_conc(i_time,idx_H2O_aq) = true_conc(0,idx_H2O_aq) + if (scenario.eq.1) then + true_conc(i_time,idx_O3) = (true_conc(0,idx_O3) - equil_O3) * & + exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3 + true_conc(i_time,idx_O3_aq_layer1) = (true_conc(0,idx_O3_aq_layer1) - equil_O3_aq) * & + exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3_aq + true_conc(i_time,idx_H2O2) = (true_conc(0,idx_H2O2) - equil_H2O2) * & + exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2 + true_conc(i_time,idx_H2O2_aq_layer1) = & + (true_conc(0,idx_H2O2_aq_layer1) - equil_H2O2_aq) * & + exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2_aq + true_conc(i_time,idx_H2O_aq_layer1) = true_conc(0,idx_H2O_aq_layer1) + true_conc(i_time,idx_H2O_aq_layer2) = true_conc(0,idx_H2O_aq_layer2) + else if (scenario.eq.2) then + true_conc(i_time,idx_O3) = (true_conc(0,idx_O3) - equil_O3) * & + exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3 + true_conc(i_time,idx_O3_aq) = (true_conc(0,idx_O3_aq) - equil_O3_aq) * & + exp(-time * (k_O3_forward + k_O3_backward)) + equil_O3_aq + true_conc(i_time,idx_H2O2) = (true_conc(0,idx_H2O2) - equil_H2O2) * & + exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2 + true_conc(i_time,idx_H2O2_aq) = & + (true_conc(0,idx_H2O2_aq) - equil_H2O2_aq) * & + exp(-time * (k_H2O2_forward + k_H2O2_backward)) + equil_H2O2_aq + true_conc(i_time,idx_H2O_aq) = true_conc(0,idx_H2O_aq) + end if end do ! Save the results @@ -539,14 +615,20 @@ logical function run_HL_phase_transfer_test(scenario) write(7,*) i_time*time_step, & ' ', true_conc(i_time, idx_O3), & ' ', model_conc(i_time, idx_O3), & - ' ', true_conc(i_time, idx_O3_aq),& - ' ', model_conc(i_time, idx_O3_aq), & + ' ', true_conc(i_time, idx_O3_aq_layer1),& + ' ', model_conc(i_time, idx_O3_aq_layer1), & + ' ', true_conc(i_time, idx_O3_aq_layer2),& + ' ', model_conc(i_time, idx_O3_aq_layer2), & ' ', true_conc(i_time, idx_H2O2), & ' ', model_conc(i_time, idx_H2O2), & - ' ', true_conc(i_time, idx_H2O2_aq), & - ' ', model_conc(i_time, idx_H2O2_aq), & - ' ', true_conc(i_time, idx_H2O_aq), & - ' ', model_conc(i_time, idx_H2O_aq) + ' ', true_conc(i_time, idx_H2O2_aq_layer1), & + ' ', model_conc(i_time, idx_H2O2_aq_layer1), & + ' ', true_conc(i_time, idx_H2O2_aq_layer2), & + ' ', model_conc(i_time, idx_H2O2_aq_layer2), & + ' ', true_conc(i_time, idx_H2O_aq_layer1), & + ' ', model_conc(i_time, idx_H2O_aq_layer1), & + ' ', true_conc(i_time, idx_H2O_aq_layer2), & + ' ', model_conc(i_time, idx_H2O_aq_layer2) end do else if (scenario.eq.2) then do i_time = 0, NUM_TIME_STEP diff --git a/test/unit_rxn_data/test_rxn_SIMPOL_phase_transfer.F90 b/test/unit_rxn_data/test_rxn_SIMPOL_phase_transfer.F90 index b914b0450..1b36807d0 100644 --- a/test/unit_rxn_data/test_rxn_SIMPOL_phase_transfer.F90 +++ b/test/unit_rxn_data/test_rxn_SIMPOL_phase_transfer.F90 @@ -109,7 +109,9 @@ logical function run_SIMPOL_phase_transfer_test(scenario) class(aero_rep_data_t), pointer :: aero_rep_ptr real(kind=dp), allocatable, dimension(:,:) :: model_conc, true_conc integer(kind=i_kind) :: idx_phase, idx_aero_rep - integer(kind=i_kind) :: idx_ethanol, idx_ethanol_aq, idx_H2O_aq + integer(kind=i_kind) :: idx_ethanol, idx_ethanol_aq_layer1, idx_H2O_aq_layer1 + integer(kind=i_kind) :: idx_ethanol_aq_layer2, idx_H2O_aq_layer2 + integer(kind=i_kind) :: idx_ethanol_aq, idx_H2O_aq integer(kind=i_kind) :: i_time, i_spec real(kind=dp) :: time_step, time real(kind=dp) :: n_star, del_H, del_S, del_G, alpha, kgm3_to_ppm @@ -147,8 +149,8 @@ logical function run_SIMPOL_phase_transfer_test(scenario) ! Allocate space for the results if (scenario.eq.1) then - allocate(model_conc(0:NUM_TIME_STEP, 7)) - allocate(true_conc(0:NUM_TIME_STEP, 7)) + allocate(model_conc(0:NUM_TIME_STEP, 9)) + allocate(true_conc(0:NUM_TIME_STEP, 9)) else if (scenario.eq.2) then allocate(model_conc(0:NUM_TIME_STEP, 6)) allocate(true_conc(0:NUM_TIME_STEP, 6)) @@ -221,21 +223,36 @@ logical function run_SIMPOL_phase_transfer_test(scenario) ! Get species indices if (scenario.eq.1) then - idx_prefix = "P1.one layer." + idx_prefix = "P1." + key = "ethanol" + idx_ethanol = chem_spec_data%gas_state_id(key); + key = idx_prefix//"one layer.aqueous aerosol.ethanol_aq" + idx_ethanol_aq_layer1 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"one layer.aqueous aerosol.H2O_aq" + idx_H2O_aq_layer1 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"two layer.aqueous aerosol.ethanol_aq" + idx_ethanol_aq_layer2 = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"two layer.aqueous aerosol.H2O_aq" + idx_H2O_aq_layer2 = aero_rep_ptr%spec_state_id(key); + ! Make sure the expected species are in the model + call assert(884352514, idx_ethanol.gt.0) + call assert(556235354, idx_ethanol_aq_layer1.gt.0) + call assert(709481366, idx_H2O_aq_layer1.gt.0) + call assert(638424219, idx_ethanol_aq_layer2.gt.0) + call assert(954729633, idx_H2O_aq_layer2.gt.0) else if (scenario.eq.2) then idx_prefix = "the mode." + key = "ethanol" + idx_ethanol = chem_spec_data%gas_state_id(key); + key = idx_prefix//"aqueous aerosol.ethanol_aq" + idx_ethanol_aq = aero_rep_ptr%spec_state_id(key); + key = idx_prefix//"aqueous aerosol.H2O_aq" + idx_H2O_aq = aero_rep_ptr%spec_state_id(key); + ! Make sure the expected species are in the model + call assert(640274151, idx_ethanol.gt.0) + call assert(554593868, idx_ethanol_aq.gt.0) + call assert(041265027, idx_H2O_aq.gt.0) endif - key = "ethanol" - idx_ethanol = chem_spec_data%gas_state_id(key); - key = idx_prefix//"aqueous aerosol.ethanol_aq" - idx_ethanol_aq = aero_rep_ptr%spec_state_id(key); - key = idx_prefix//"aqueous aerosol.H2O_aq" - idx_H2O_aq = aero_rep_ptr%spec_state_id(key); - - ! Make sure the expected species are in the model - call assert(884352514, idx_ethanol.gt.0) - call assert(379146109, idx_ethanol_aq.gt.0) - call assert(208989205, idx_H2O_aq.gt.0) #ifdef CAMP_USE_MPI ! pack the camp core @@ -262,6 +279,10 @@ logical function run_SIMPOL_phase_transfer_test(scenario) ! broadcast the species ids call camp_mpi_bcast_integer(idx_ethanol) + call camp_mpi_bcast_integer(idx_ethanol_aq_layer1) + call camp_mpi_bcast_integer(idx_H2O_aq_layer1) + call camp_mpi_bcast_integer(idx_ethanol_aq_layer2) + call camp_mpi_bcast_integer(idx_H2O_aq_layer2) call camp_mpi_bcast_integer(idx_ethanol_aq) call camp_mpi_bcast_integer(idx_H2O_aq) call camp_mpi_bcast_integer(i_sect_unused) @@ -320,19 +341,31 @@ logical function run_SIMPOL_phase_transfer_test(scenario) ! Save the initial concentrations true_conc(:,:) = 0.0 - true_conc(0,idx_ethanol) = 1.0e-1 - true_conc(0,idx_ethanol_aq) = 1.0e-8 - true_conc(0,idx_H2O_aq) = 1.4e-2 + if (scenario.eq.1) then + true_conc(0,idx_ethanol) = 1.0e-1 + true_conc(0,idx_ethanol_aq_layer1) = 1.0e-8 + true_conc(0,idx_H2O_aq_layer1) = 1.4e-2 + true_conc(0,idx_ethanol_aq_layer2) = 0.1e-8 + true_conc(0,idx_H2O_aq_layer2) = 0.1e-2 + else if (scenario.eq.2) then + true_conc(0,idx_ethanol) = 1.0e-1 + true_conc(0,idx_ethanol_aq) = 1.0e-8 + true_conc(0,idx_H2O_aq) = 1.4e-2 + end if ! Calculate the radius and number concentration to use if (scenario.eq.1) then number_conc = 1.3e6 ! particle number concentration (#/cc) ! single particle aerosol mass concentrations are per particle - true_conc(0,idx_ethanol_aq) = true_conc(0,idx_ethanol_aq) / number_conc - true_conc(0,idx_H2O_aq) = true_conc(0,idx_H2O_aq) / number_conc + true_conc(0,idx_ethanol_aq_layer1) = true_conc(0,idx_ethanol_aq_layer1) / number_conc + true_conc(0,idx_H2O_aq_layer1) = true_conc(0,idx_H2O_aq_layer1) / number_conc + true_conc(0,idx_ethanol_aq_layer2) = true_conc(0,idx_ethanol_aq_layer2) / number_conc + true_conc(0,idx_H2O_aq_layer2) = true_conc(0,idx_H2O_aq_layer2) / number_conc ! radius (m) calculated based on particle mass - radius = ( ( true_conc(0,idx_ethanol_aq) / 1000.0 + & - true_conc(0,idx_H2O_aq) / 1000.0 ) & + radius = ( ( true_conc(0,idx_ethanol_aq_layer1) / 1000.0 + & + true_conc(0,idx_H2O_aq_layer1) / 1000.0 + & + true_conc(0,idx_ethanol_aq_layer2) / 1000.0 + & + true_conc(0,idx_H2O_aq_layer2) / 1000.0 ) & * 3.0 / 4.0 / 3.14159265359 )**(1.0/3.0) else if (scenario.eq.2) then ! radius (m) @@ -346,7 +379,7 @@ logical function run_SIMPOL_phase_transfer_test(scenario) model_conc(0,:) = true_conc(0,:) - ! Update the aerosol representation (single partile only) + ! Update the aerosol representation (single particle only) if (scenario.eq.1) then call number_update%set_number__n_m3(1, number_conc) call camp_core%update_data(number_update) @@ -392,7 +425,7 @@ logical function run_SIMPOL_phase_transfer_test(scenario) if (scenario.eq.1) then ! aerosol mass concentrations are per particle total_mass = true_conc(0,idx_ethanol)/kgm3_to_ppm + & - true_conc(0,idx_ethanol_aq)*number_conc ! (kg/m3) + (true_conc(0,idx_ethanol_aq_layer1))*number_conc! (kg/m3) else if (scenario.eq.2) then ! aerosol mass concentrations for the total mode total_mass = true_conc(0,idx_ethanol)/kgm3_to_ppm + & @@ -446,37 +479,68 @@ logical function run_SIMPOL_phase_transfer_test(scenario) ! [A_aero] = ([A_init_aero] - [A_eq_aero]) * exp(-t * (k_f + k_b)) + ! [A_eq_aero] time = i_time * time_step - true_conc(i_time,idx_ethanol) = & - (true_conc(0,idx_ethanol) - equil_ethanol) * & - exp(-time * (k_forward + k_backward)) + equil_ethanol - true_conc(i_time,idx_ethanol_aq) = & - (true_conc(0,idx_ethanol_aq) - equil_ethanol_aq) * & - exp(-time * (k_forward + k_backward)) + equil_ethanol_aq - true_conc(i_time,idx_H2O_aq) = true_conc(0,idx_H2O_aq) + if (scenario.eq.1) then + true_conc(i_time,idx_ethanol) = & + (true_conc(0,idx_ethanol) - equil_ethanol) * & + exp(-time * (k_forward + k_backward)) + equil_ethanol + true_conc(i_time,idx_ethanol_aq_layer1) = & + (true_conc(0,idx_ethanol_aq_layer1) - equil_ethanol_aq) * & + exp(-time * (k_forward + k_backward)) + equil_ethanol_aq + true_conc(i_time,idx_ethanol_aq_layer2) = true_conc(0,idx_ethanol_aq_layer2) + true_conc(i_time,idx_H2O_aq_layer1) = true_conc(0,idx_H2O_aq_layer1) + true_conc(i_time,idx_H2O_aq_layer2) = true_conc(0,idx_H2O_aq_layer2) + else if (scenario.eq.2) then + true_conc(i_time,idx_ethanol) = & + (true_conc(0,idx_ethanol) - equil_ethanol) * & + exp(-time * (k_forward + k_backward)) + equil_ethanol + true_conc(i_time,idx_ethanol_aq) = & + (true_conc(0,idx_ethanol_aq) - equil_ethanol_aq) * & + exp(-time * (k_forward + k_backward)) + equil_ethanol_aq + true_conc(i_time,idx_H2O_aq) = true_conc(0,idx_H2O_aq) + end if end do + ! Save the results if (scenario.eq.1) then - open(unit=7, file="out/SIMPOL_phase_transfer_results.txt", & + open(unit=9, file="out/SIMPOL_phase_transfer_results.txt", & status="replace", action="write") else if (scenario.eq.2) then open(unit=7, file="out/SIMPOL_phase_transfer_results_2.txt", & status="replace", action="write") endif do i_time = 0, NUM_TIME_STEP - write(7,*) i_time*time_step, & - ' ', true_conc(i_time, idx_ethanol), & - ' ', model_conc(i_time, idx_ethanol), & - ' ', true_conc(i_time, idx_ethanol_aq), & - ' ', model_conc(i_time, idx_ethanol_aq), & - ' ', true_conc(i_time, idx_H2O_aq), & - ' ', model_conc(i_time, idx_H2O_aq) + if (scenario.eq.1) then + write(9,*) i_time*time_step, & + ' ', true_conc(i_time, idx_ethanol), & + ' ', model_conc(i_time, idx_ethanol), & + ' ', true_conc(i_time, idx_ethanol_aq_layer1), & + ' ', model_conc(i_time, idx_ethanol_aq_layer1), & + ' ', true_conc(i_time, idx_ethanol_aq_layer2), & + ' ', model_conc(i_time, idx_ethanol_aq_layer2), & + ' ', true_conc(i_time, idx_H2O_aq_layer1), & + ' ', model_conc(i_time, idx_H2O_aq_layer1), & + ' ', true_conc(i_time, idx_H2O_aq_layer2), & + ' ', model_conc(i_time, idx_H2O_aq_layer2) + else if (scenario.eq.2) then + write(7,*) i_time*time_step, & + ' ', true_conc(i_time, idx_ethanol), & + ' ', model_conc(i_time, idx_ethanol), & + ' ', true_conc(i_time, idx_ethanol_aq), & + ' ', model_conc(i_time, idx_ethanol_aq), & + ' ', true_conc(i_time, idx_H2O_aq), & + ' ', model_conc(i_time, idx_H2O_aq) + end if end do - close(7) + if (scenario.eq.1) then + close(9) + else if (scenario.eq.2) then + close(7) + endif ! Analyze the results ! - ! The particle radius changes as ethanol condenses/evaporates, so an + ! The particle radius changes as ethanol condenses/evaporates, so ! an exact solution is not calculated. The tolerances on the comparison ! with "true" values are higher to account for this. do i_time = 1, NUM_TIME_STEP @@ -491,7 +555,8 @@ logical function run_SIMPOL_phase_transfer_test(scenario) (model_conc(i_time, i_spec).lt.1.2*model_conc(NUM_TIME_STEP, i_spec).and. & true_conc(i_time, i_spec).lt.1.2*true_conc(NUM_TIME_STEP, i_spec)).or. & (model_conc(i_time, i_spec).lt.1e-2*model_conc(1, i_spec).and. & - true_conc(i_time, i_spec).lt.1e-2*true_conc(1, i_spec)), & + true_conc(i_time, i_spec).lt.1e-2*true_conc(1, i_spec)) .and. & + model_conc(i_time,idx_ethanol_aq_layer2).eq.model_conc(i_time,idx_ethanol_aq_layer2), & "time: "//trim(to_string(i_time))//"; species: "// & trim(to_string(i_spec))//"; mod: "// & trim(to_string(model_conc(i_time, i_spec)))//"; true: "// &