Merge pull request #8 from open-atmos/develop-debug-SIMPOL

jtgasparik · web-flow · commit b2e9b2af6f69 · 2024-10-15T17:46:43.000-04:00
Update unique names function to have is_at_surface flag
diff --git 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, num_is_at_surface)
+  function num_phase_instances(this, phase_name, is_at_surface)
     use camp_util,                                       only : i_kind
     import :: aero_rep_data_t
 
@@ -353,7 +356,7 @@ function num_phase_instances(this, phase_name, num_is_at_surface)
     !> Aerosol phase name
     character(len=*), intent(in) :: phase_name
     !> Indicates if aerosol phase is at the surface of particle
-    logical, intent(in), optional :: num_is_at_surface
+    logical, intent(in), optional :: is_at_surface
 
   end function num_phase_instances
 
@@ -511,10 +514,10 @@ function phase_ids(this, phase_name, is_at_surface)
     if (present(is_at_surface)) then
       if (is_at_surface) then
         num_instances = this%num_phase_instances(phase_name, &
-                 num_is_at_surface = .true.)
+                 is_at_surface = .true.)
       else
         num_instances = this%num_phase_instances(phase_name, &
-                 num_is_at_surface = .false.)
+                 is_at_surface = .false.)
       end if
     else
       num_instances = this%num_phase_instances(phase_name)
diff --git 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, num_is_at_surface)
+  function num_phase_instances(this, phase_name, is_at_surface)
 
     !> Number of instances of the aerosol phase
     integer(kind=i_kind) :: num_phase_instances
@@ -950,21 +968,20 @@ function num_phase_instances(this, phase_name, num_is_at_surface)
     !> Aerosol phase name
     character(len=*), intent(in) :: phase_name
     !> Indicates if aerosol phase is at the surface of particle
-    logical, intent(in), optional :: num_is_at_surface
+    logical, intent(in), optional :: is_at_surface
 
     integer(kind=i_kind) :: i_phase
 
-    if (present(num_is_at_surface)) then
-      if (num_is_at_surface) then
-        num_phase_instances = 0
+    num_phase_instances = 0
+    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
-        num_phase_instances = 0
         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
@@ -975,7 +992,6 @@ function num_phase_instances(this, phase_name, num_is_at_surface)
         end do 
       end if
     else
-      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
diff --git a/src/aero_reps/aero_rep_single_particle.F90 b/src/aero_reps/aero_rep_single_particle.F90
@@ -487,7 +487,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
@@ -500,6 +501,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
@@ -520,6 +523,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)
@@ -549,6 +556,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
@@ -641,27 +652,25 @@ end function spec_name
   !! 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, &
-                                                    num_is_at_surface)
+                                                    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 :: num_is_at_surface
+    logical, intent(in), optional :: is_at_surface
 
     integer(kind=i_kind) ::  i_phase, i_layer, phase_index
 
-
     num_phase_instances = 0
-    phase_index = 0
-    if (present(num_is_at_surface)) then
-      if (num_is_at_surface) then
+    if (present(is_at_surface)) then
+      if (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)) then
-              phase_index = phase_index + 1
+              num_phase_instances = num_phase_instances + 1
             end if
           end do
         end do
@@ -670,7 +679,7 @@ integer(kind=i_kind) function num_phase_instances(this, phase_name, &
           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. &
               .not. this%aero_phase_is_at_surface(i_phase)) then
-              phase_index = phase_index + 1
+              num_phase_instances = num_phase_instances + 1
             end if
           end do
         end do
@@ -679,12 +688,12 @@ integer(kind=i_kind) function num_phase_instances(this, phase_name, &
       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
+            num_phase_instances = num_phase_instances + 1
           end if
         end do
       end do
     end if
-    num_phase_instances = phase_index * MAX_PARTICLES_
+    num_phase_instances = num_phase_instances * MAX_PARTICLES_
 
   end function num_phase_instances
 
diff --git a/src/rxns/rxn_SIMPOL_phase_transfer.F90 b/src/rxns/rxn_SIMPOL_phase_transfer.F90
@@ -246,7 +246,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
@@ -308,12 +309,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_single_particle.F90 b/test/unit_aero_rep_data/test_aero_rep_single_particle.F90
@@ -215,19 +215,19 @@ subroutine test_config_read()
         max_part = aero_rep%maximum_computational_particles()
         jam_phase_instance = 0
         !check value
-        print *, "jam_instances ", aero_rep%num_phase_instances(phase_name_test, num_is_at_surface = .false.)
+        print *, "jam_instances ", aero_rep%num_phase_instances(phase_name_test, is_at_surface = .false.)
         call assert(417730478, 3 .eq. max_part)
         call assert(493602373, jam_phase_instance .eq. aero_rep%num_phase_instances(phase_name_test, &
-                                                         num_is_at_surface = .true.))
+                                                         is_at_surface = .true.))
 
         phase_name_test = "bread"
         num_bread = 1
         bread_phase_instance = num_bread * max_part
         !check value
         print *, bread_phase_instance
-        print *, "bread_instqnces ", aero_rep%num_phase_instances(phase_name_test, num_is_at_surface = .true.)
+        print *, "bread_instqnces ", aero_rep%num_phase_instances(phase_name_test, is_at_surface = .true.)
         call assert(734138496, bread_phase_instance .eq. aero_rep%num_phase_instances(phase_name_test, &
-                                                         num_is_at_surface = .true.))
+                                                         is_at_surface = .true.))
 
       class default
         call die_msg(519535557, rep_name)
