Simpol multiple layers

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
 

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,23 +959,45 @@ 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
     !> Aerosol representation data
     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
 

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
 

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",

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 =

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.)

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"// &