SIMPOL with two layers and some mass in second layer

src/rxns/rxn_SIMPOL_phase_transfer.F90

@@ -1,4 +1,4 @@
-! Copyright (C) 2021 Barcelona :/repSupercomputing Center and University of
+! Copyright (C) 2021 Barcelona Supercomputing Center and University of
 ! Illinois at Urbana-Champaign
 ! SPDX-License-Identifier: MIT
 
@@ -221,7 +221,6 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
             this%property_set%get_string(key_name, aero_spec_name), &
             "Missing aerosol-phase species in SIMPOL.1 phase-transfer "// &
             "reaction")
-    print *, "species name ", aero_spec_name
 
     ! Set up a general error message
     error_msg = " for SIMPOL.1 phase transfer of gas species '"// &
@@ -258,7 +257,6 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
       ! Get the number of Jacobian elements for calculations of mass, volume,
       ! number, etc. for this partitioning into this phase
       phase_ids = aero_rep(i_aero_rep)%val%phase_ids(phase_name, is_at_surface=.true.)
-      print *, "phase_ids", phase_ids
       do i_phase = 1, size(phase_ids)
         n_aero_jac_elem = n_aero_jac_elem + &
                 aero_rep(i_aero_rep)%val%num_jac_elem(phase_ids(i_phase))
@@ -425,18 +423,12 @@ subroutine initialize(this, chem_spec_data, aero_rep, n_cells)
     ! Check the sizes of the data arrays
     tmp_size = PHASE_INT_LOC_(i_aero_id - 1) + 1 + &
                2*NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1) - 1
-    print *, "PHASE_INT_LOC_", PHASE_INT_LOC_(i_aero_id - 1)
-    print *, "NUM_AERO_PHASE_JAC_ELEM_", NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1)
-    print *, "tmp_size_1", tmp_size
-    print *, "size condensed data int",size(this%condensed_data_int)
     call assert_msg(625802519, size(this%condensed_data_int) .eq. tmp_size, &
                     "int array size mismatch"//error_msg)
     tmp_size = PHASE_REAL_LOC_(i_aero_id - 1) + &
                4*NUM_AERO_PHASE_JAC_ELEM_(i_aero_id - 1) - 1
     call assert_msg(391089510, size(this%condensed_data_real) .eq. tmp_size, &
                     "real array size mismatch"//error_msg)
-    print *, "tmp_size_2", tmp_size
-    print *, "size condensed data real",size(this%condensed_data_real)
 
   end subroutine initialize
 

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"
                     ],
@@ -54,14 +54,14 @@
     "maximum computational particles" : 1,
     "layers": [
             {
-                    "name": "one layer",
+                    "name": "one_layer",
                     "phases": [
                                "aqueous aerosol"
                     ],
-                    "covers": "two layer"
+                    "covers": "two_layer"
             },
             {
-                    "name": "two layer",
+                    "name": "two_layer",
                     "phases": [
                                "aqueous aerosol"
                     ],
@@ -75,7 +75,7 @@
     "maximum computational particles" : 1,
     "layers": [
             {
-                    "name": "one layer",
+                    "name": "one_layer",
                     "phases": [
                                "aqueous aerosol"
                     ],

test/unit_rxn_data/test_rxn_SIMPOL_phase_transfer.F90

@@ -226,13 +226,13 @@ logical function run_SIMPOL_phase_transfer_test(scenario)
         idx_prefix = "P1."
         key = "ethanol"
         idx_ethanol = chem_spec_data%gas_state_id(key);
-        key = idx_prefix//"one layer.aqueous aerosol.ethanol_aq"
+        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"
+        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"
+        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"
+        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)
@@ -345,8 +345,8 @@ logical function run_SIMPOL_phase_transfer_test(scenario)
         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.0e0
-        true_conc(0,idx_H2O_aq_layer2) = 0.0e0
+        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
@@ -367,7 +367,6 @@ logical function run_SIMPOL_phase_transfer_test(scenario)
                      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)
-        print *, 'radius', radius
       else if (scenario.eq.2) then
         ! radius (m)
         radius = 9.37e-7 / 2.0 * exp(5.0/2.0 * log(0.9) * log(0.9))
@@ -426,9 +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_layer1)+ &
-                      true_conc(0,idx_ethanol_aq_layer2))*number_conc! (kg/m3)
-        print *, 'total mass', total_mass
+                     (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 + &
@@ -503,27 +500,26 @@ logical function run_SIMPOL_phase_transfer_test(scenario)
         end if
       end do
 
-      print *, "layer 1 ethanol conc", model_conc(:,idx_ethanol_aq_layer1)
 
       ! Save the results
       if (scenario.eq.1) then
-        open(unit=11, 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
         if (scenario.eq.1) then
-          write(11,*) i_time*time_step, &
+          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),       &
+                ' ', 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
@@ -536,15 +532,19 @@ logical function run_SIMPOL_phase_transfer_test(scenario)
                 ' ', model_conc(i_time, idx_H2O_aq)
         end if    
       end do
-      close(11)
+      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 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
-        do i_spec = 1, 7
+        do i_spec = 1, 5
           ! scenario 1 - Only check the second phase
           if (scenario.eq.1.and.i_spec.ge.2.and.i_spec.le.3) cycle
           ! scenario 2 - exclude activity coefficient