From f9e0d8dcb84537c1b5c337a8177124c600855fc4 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 20 Feb 2025 17:50:58 -0800 Subject: [PATCH 01/12] rename, remove boilerplate --- .../aerosol/cloud_aqueous_chemistry.F90 | 253 +++++++++--------- .../test_cloud_aqueous_chemistry.F90 | 14 +- 2 files changed, 138 insertions(+), 129 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 48cf96251a..2ddeda9961 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -1,5 +1,13 @@ !---------------------------------------------------------------------------------- ! Cloud aqueous chemistry +! +! The purpose of this module is to calculate the sulfate formation due to various +! oxidation/condensation pathways of cloud chemistry. These pathways include: +! - SO2 oxidation by O3 +! - SO2 oxidation by H2O2 +! - H2SO4 condensation +! Updated gas and aqueous species concentrations along with cloud pH changes are +! calculated. !---------------------------------------------------------------------------------- module cloud_aqueous_chemistry @@ -14,25 +22,37 @@ module cloud_aqueous_chemistry implicit none private - public :: sox_inti, setsox - public :: has_sox + public :: initialize, calculate + public :: do_cloud_aqueous_chemistry - logical :: inv_o3 - integer :: id_msa + logical :: do_cloud_aqueous_chemistry = .false. - integer :: id_so2, id_nh3, id_hno3, id_h2o2, id_o3, id_ho2 - integer :: id_so4, id_h2so4 + integer, parameter :: CLOUD_INDEX_UNDEFINED = -1 - logical :: has_sox = .true. - logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2 + ! Cloud chemistry species information + type :: cloud_species_t + logical :: is_constant_ = .false. + integer :: state_index_ = CLOUD_INDEX_UNDEFINED + character(len=:), allocatable :: name_ + contains + procedure :: exists => cloud_species_exists + procedure :: mixing_ratio => cloud_species_get_mixing_ratio + end type cloud_species_t + interface cloud_species_t + module procedure :: cloud_species_constructor + end interface cloud_species_t + + type(cloud_species_t) :: so2, nh3, hno3, h2o2, o3, ho2, msa, so4, h2so4 + + ! TODO: Figure out what this flag is for logical :: cloud_borne = .false. contains !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine sox_inti + subroutine initialize() !----------------------------------------------------------------------- ! ... initialize the hetero sox routine !----------------------------------------------------------------------- @@ -57,82 +77,42 @@ subroutine sox_inti ! ... get species indicies !----------------------------------------------------------------- + so2 = cloud_species_t( 'SO2' ) + nh3 = cloud_species_t( 'NH3' ) + hno3 = cloud_species_t( 'HNO3' ) + h2o2 = cloud_species_t( 'H2O2' ) + o3 = cloud_species_t( 'O3' ) + ho2 = cloud_species_t( 'HO2' ) + msa = cloud_species_t( 'MSA' ) + so4 = cloud_species_t( 'SO4' ) + h2so4 = cloud_species_t( 'H2SO4' ) + + do_cloud_aqueous_chemistry = so2%exists() .and. h2o2%exists() .and. & + o3%exists() .and. ho2%exists() if (cloud_borne) then - id_h2so4 = get_spc_ndx( 'H2SO4' ) - else - id_so4 = get_spc_ndx( 'SO4' ) - endif - id_msa = get_spc_ndx( 'MSA' ) - - inv_so2 = .false. - id_so2 = get_inv_ndx( 'SO2' ) - inv_so2 = id_so2 > 0 - if ( .not. inv_so2 ) then - id_so2 = get_spc_ndx( 'SO2' ) - endif - - inv_NH3 = .false. - id_NH3 = get_inv_ndx( 'NH3' ) - inv_NH3 = id_NH3 > 0 - if ( .not. inv_NH3 ) then - id_NH3 = get_spc_ndx( 'NH3' ) - endif - - inv_HNO3 = .false. - id_HNO3 = get_inv_ndx( 'HNO3' ) - inv_HNO3 = id_hno3 > 0 - if ( .not. inv_HNO3 ) then - id_HNO3 = get_spc_ndx( 'HNO3' ) - endif - - inv_H2O2 = .false. - id_H2O2 = get_inv_ndx( 'H2O2' ) - inv_H2O2 = id_H2O2 > 0 - if ( .not. inv_H2O2 ) then - id_H2O2 = get_spc_ndx( 'H2O2' ) - endif - - inv_HO2 = .false. - id_HO2 = get_inv_ndx( 'HO2' ) - inv_HO2 = id_HO2 > 0 - if ( .not. inv_HO2 ) then - id_HO2 = get_spc_ndx( 'HO2' ) - endif - - inv_o3 = get_inv_ndx( 'O3' ) > 0 - if (inv_o3) then - id_o3 = get_inv_ndx( 'O3' ) + do_cloud_aqueous_chemistry = do_cloud_aqueous_chemistry .and. & + h2so4%exists() else - id_o3 = get_spc_ndx( 'O3' ) - endif - inv_ho2 = get_inv_ndx( 'HO2' ) > 0 - if (inv_ho2) then - id_ho2 = get_inv_ndx( 'HO2' ) - else - id_ho2 = get_spc_ndx( 'HO2' ) - endif - - has_sox = (id_so2>0) .and. (id_h2o2>0) .and. (id_o3>0) .and. (id_ho2>0) - if (cloud_borne) then - has_sox = has_sox .and. (id_h2so4>0) - else - has_sox = has_sox .and. (id_so4>0) .and. (id_nh3>0) + do_cloud_aqueous_chemistry = do_cloud_aqueous_chemistry .and. & + so4%exists() .and. nh3%exists() endif if (masterproc) then - write(iulog,*) 'sox_inti: has_sox = ',has_sox + write(iulog,*) 'cloud_aqueous_chemistry_init: '//& + 'do_cloud_aqueous_chemistry = ',& + do_cloud_aqueous_chemistry endif - if( has_sox ) then + if( do_cloud_aqueous_chemistry ) then if (masterproc) then write(iulog,*) '-----------------------------------------' - write(iulog,*) ' mo_setsox will do sox aerosols' + write(iulog,*) ' cloud aqueous chemistry is active' write(iulog,*) '-----------------------------------------' endif else if (masterproc) then write(iulog,*) '-----------------------------------------' - write(iulog,*) ' mo_setsox will not do sox aerosols' + write(iulog,*) ' cloud aqueous chemistry is inactive' write(iulog,*) '-----------------------------------------' endif return @@ -140,11 +120,11 @@ subroutine sox_inti call sox_cldaero_init() - end subroutine sox_inti + end subroutine initialize !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine setsox( state, & + subroutine calculate( state, & pbuf, & ncol, & lchnk, & @@ -335,57 +315,20 @@ subroutine setsox( state, & xso4c => cldconc%so4c xnh4c => cldconc%nh4c xno3c => cldconc%no3c - xso4(:,:) = 0._r8 xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 - - do k = 1,pver - xph(:,k) = xph0 ! initial PH value - - if ( inv_so2 ) then - xso2 (:,k) = invariants(:,k,id_so2)/xhnm(:,k) ! mixing ratio - else - xso2 (:,k) = qin(:,k,id_so2) ! mixing ratio - endif - - if (id_hno3 > 0) then - xhno3(:,k) = qin(:,k,id_hno3) - else - xhno3(:,k) = 0.0_r8 - endif - - if ( inv_h2o2 ) then - xh2o2 (:,k) = invariants(:,k,id_h2o2)/xhnm(:,k) ! mixing ratio - else - xh2o2 (:,k) = qin(:,k,id_h2o2) ! mixing ratio - endif - - if (id_nh3 > 0) then - xnh3 (:,k) = qin(:,k,id_nh3) - else - xnh3 (:,k) = 0.0_r8 - endif - - if ( inv_o3 ) then - xo3 (:,k) = invariants(:,k,id_o3)/xhnm(:,k) ! mixing ratio - else - xo3 (:,k) = qin(:,k,id_o3) ! mixing ratio - endif - if ( inv_ho2 ) then - xho2 (:,k) = invariants(:,k,id_ho2)/xhnm(:,k)! mixing ratio - else - xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio - endif - - if (cloud_borne) then - xh2so4(:,k) = qin(:,k,id_h2so4) - else - xso4 (:,k) = qin(:,k,id_so4) ! mixing ratio - endif - if (id_msa > 0) xmsa (:,k) = qin(:,k,id_msa) - - end do + call so2%mixing_ratio( qin, invariants, xhnm, xso2 ) + call nh3%mixing_ratio( qin, invariants, xhnm, xnh3 ) + call hno3%mixing_ratio( qin, invariants, xhnm, xhno3 ) + call h2o2%mixing_ratio( qin, invariants, xhnm, xh2o2 ) + call o3%mixing_ratio( qin, invariants, xhnm, xo3 ) + call ho2%mixing_ratio( qin, invariants, xhnm, xho2 ) + call msa%mixing_ratio( qin, invariants, xhnm, xmsa ) + call so4%mixing_ratio( qin, invariants, xhnm, xso4 ) + call h2so4%mixing_ratio( qin, invariants, xhnm, xh2so4 ) + + xph(:,:) = xph0 ! initial PH value !----------------------------------------------------------------- ! ... Temperature dependent Henry constants @@ -750,7 +693,7 @@ subroutine setsox( state, & ! ... nh3 !------------------------------------------------------------------------ px = henh3(i,k) * Ra * tz * xl - if (id_nh3>0) then + if (nh3%exists()) then nh3g(i,k) = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) else nh3g(i,k) = 0._r8 @@ -896,7 +839,71 @@ subroutine setsox( state, & call sox_cldaero_destroy_obj(cldconc) - end subroutine setsox + end subroutine calculate + +!------------------------------------------------------------------------------- +! +! Support routines to be moved to separate modules +! +!------------------------------------------------------------------------------- + + !------------------------------------------------------------------------------- + ! Creates a cloud species object with state indices + function cloud_species_constructor( species_name ) result( this ) + + use mo_chem_utls, only : get_spc_ndx, get_inv_ndx + + type(cloud_species_t) :: this + character(len=*), intent(in) :: species_name + + this%name_ = species_name + this%state_index_ = get_inv_ndx( species_name ) + this%is_constant_ = this%state_index_ > 0 + if ( .not. this%is_constant_ ) & + this%state_index_ = get_spc_ndx( species_name ) + if ( this%state_index_ <= 0 ) this%state_index_ = CLOUD_INDEX_UNDEFINED + + end function cloud_species_constructor + + !------------------------------------------------------------------------------- + ! Returns whether a cloud species is defined + logical function cloud_species_exists( this ) + + class(cloud_species_t), intent(in) :: this + + cloud_species_exists = this%state_index_ .ne. CLOUD_INDEX_UNDEFINED + + end function cloud_species_exists + + !------------------------------------------------------------------------------- + ! Returns the mixing ratio for a cloud species + ! + ! Constant species use the fixed number concentration and the air density + ! to calculate the mixing ratio. Time-varying species simply return the + ! mixing ratio from the state array. + ! + ! For species that are not defined, the mixing ratio is set to zero. + subroutine cloud_species_get_mixing_ratio( this, mixing_ratios, & + fixed_concentrations, air_density, mixing_ratio ) + + class(cloud_species_t), intent(in) :: this + real(r8), intent(in) :: mixing_ratios(:,:,:) ! all mixing ratios (kg kg-1) [column, layer, species] + real(r8), intent(in) :: fixed_concentrations(:,:,:) ! all fixed concentrations (kg m-3) [column, layer, species] + real(r8), intent(in) :: air_density(:,:) ! air density (kg m-3) [column, layer] + real(r8), intent(out) :: mixing_ratio(:,:) ! species mixing ratio (kg kg-1) [column, layer] + + if ( this%state_index_ == CLOUD_INDEX_UNDEFINED ) then + mixing_ratio(:,:) = 0._r8 + return + end if + if ( this%is_constant_ ) then + mixing_ratio(:,:) = fixed_concentrations(:,:,this%state_index_) & + / air_density(:,:) + else + mixing_ratio(:,:) = mixing_ratios(:,:,this%state_index_) + end if + + end subroutine cloud_species_get_mixing_ratio end module cloud_aqueous_chemistry diff --git a/test/chemistry/cloud_aqueous_chemistry/test_cloud_aqueous_chemistry.F90 b/test/chemistry/cloud_aqueous_chemistry/test_cloud_aqueous_chemistry.F90 index c56e3f3f2d..66033bf166 100644 --- a/test/chemistry/cloud_aqueous_chemistry/test_cloud_aqueous_chemistry.F90 +++ b/test/chemistry/cloud_aqueous_chemistry/test_cloud_aqueous_chemistry.F90 @@ -57,7 +57,8 @@ subroutine test_as_primary_rank_against_snapshot() use spmd_utils, only: masterproc use physics_buffer, only: physics_buffer_desc use physics_types, only: physics_state - use cloud_aqueous_chemistry, only: sox_inti, setsox + use cloud_aqueous_chemistry, only: sox_inti => initialize, & + setsox => calculate type(chemistry_args), allocatable :: chem_args(:), expected_outputs(:) type(physics_buffer_desc), pointer :: pbuf(:) @@ -111,7 +112,8 @@ subroutine test_as_other_rank_against_snapshot() use spmd_utils, only: masterproc use physics_buffer, only: physics_buffer_desc use physics_types, only: physics_state - use cloud_aqueous_chemistry, only: sox_inti, setsox + use cloud_aqueous_chemistry, only: sox_inti => initialize, & + setsox => calculate type(chemistry_args), allocatable :: chem_args(:), expected_outputs(:) type(physics_buffer_desc), pointer :: pbuf(:) @@ -165,8 +167,8 @@ subroutine test_as_primary_against_original_module() use spmd_utils, only: masterproc use physics_buffer, only: physics_buffer_desc use physics_types, only: physics_state - use cloud_aqueous_chemistry, only: new_sox_inti => sox_inti, & - new_setsox => setsox + use cloud_aqueous_chemistry, only: new_sox_inti => initialize, & + new_setsox => calculate use mo_setsox, only: old_sox_inti => sox_inti, old_setsox => setsox type(chemistry_args), allocatable :: new_args(:), old_args(:) @@ -248,8 +250,8 @@ subroutine test_as_other_rank_against_original_module() use spmd_utils, only: masterproc use physics_buffer, only: physics_buffer_desc use physics_types, only: physics_state - use cloud_aqueous_chemistry, only: new_sox_inti => sox_inti, & - new_setsox => setsox + use cloud_aqueous_chemistry, only: new_sox_inti => initialize, & + new_setsox => calculate use mo_setsox, only: old_sox_inti => sox_inti, old_setsox => setsox type(chemistry_args), allocatable :: new_args(:), old_args(:) From 90e3f1736b5cddb8bebc7433eedc721a38b665e4 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 10:53:20 -0800 Subject: [PATCH 02/12] rename function args --- .../aerosol/cloud_aqueous_chemistry.F90 | 187 +++++++++--------- 1 file changed, 95 insertions(+), 92 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 2ddeda9961..039d49ca30 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -124,31 +124,31 @@ end subroutine initialize !----------------------------------------------------------------------- !----------------------------------------------------------------------- - subroutine calculate( state, & - pbuf, & - ncol, & - lchnk, & - loffset,& - dtime, & - press, & - pdel, & - tfld, & - mbar, & - lwc, & - cldfrc, & - cldnum, & - xhnm, & - invariants, & - qcw, & - qin, & - xphlwc, & - aqso4, & - aqh2so4,& - aqso4_h2o2, & - aqso4_o3, & - yph_in, & - aqso4_h2o2_3d, & - aqso4_o3_3d & + subroutine calculate( state, & + pbuf, & + ncol, & + lchnk, & + loffset, & + time_step, & + midpoint_pressure, & + pressure_thickness, & + temperature, & + mean_mass, & + cloud_water, & + cloud_fraction, & + cloud_droplet_number, & + atmospheric_number_density, & + fixed_concentrations, & + cloud_borne_aerosol_vmr, & + species_vmr, & + ph_times_cloud_water, & + aq_so4_production, & + aq_h2so4_production, & + aq_so4_production_from_h2o2, & + aq_so4_production_from_o3, & + specified_ph, & + aq_so4_production_from_h2o2_3d, & + aq_so4_production_from_o3_3d & ) !----------------------------------------------------------------------- @@ -172,10 +172,12 @@ subroutine calculate( state, & use physconst, only : mwdry, gravit use mo_constants, only : pi #ifdef USE_MAM - use mam_clouds, only : sox_cldaero_update, sox_cldaero_create_obj, sox_cldaero_destroy_obj + use mam_clouds, only : sox_cldaero_update, sox_cldaero_create_obj, & + sox_cldaero_destroy_obj #endif #ifdef USE_CARMA - use carma_clouds, only : sox_cldaero_update, sox_cldaero_create_obj, sox_cldaero_destroy_obj + use carma_clouds, only : sox_cldaero_update, sox_cldaero_create_obj, & + sox_cldaero_destroy_obj #endif use cloud_utilities, only : cldaero_conc_t @@ -183,30 +185,30 @@ subroutine calculate( state, & !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- - integer, intent(in) :: ncol ! num of columns in chunk - integer, intent(in) :: lchnk ! chunk id - integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array - real(r8), intent(in) :: dtime ! time step (sec) - real(r8), intent(in) :: press(:,:) ! midpoint pressure ( Pa ) - real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa) - real(r8), intent(in) :: tfld(:,:) ! temperature - real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu ) - real(r8), target, intent(in) :: lwc(:,:) ! cloud liquid water content (kg/kg) - real(r8), target, intent(in) :: cldfrc(:,:) ! cloud fraction - real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg) - real(r8), intent(in) :: xhnm(:,:) ! total atms density ( /cm**3) - real(r8), intent(in) :: invariants(:,:,:) - real(r8), target, intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) - real(r8), intent(inout) :: qin(:,:,:) ! transported species ( vmr ) - real(r8), intent(out) :: xphlwc(:,:) ! pH value multiplied by lwc - - real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry - real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) - real(r8), intent(in), optional :: yph_in ! ph value - real(r8), intent(out), optional :: aqso4_h2o2_3d(:, :) ! 3D SO4 aqueous phase chemistry due to H2O2 (kg/m2) - real(r8), intent(out), optional :: aqso4_o3_3d(:, :) ! 3D SO4 aqueous phase chemistry due to O3 (kg/m2) + integer, intent(in) :: ncol ! num of columns in chunk + integer, intent(in) :: lchnk ! chunk id + integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array + real(r8), intent(in) :: time_step ! time step (sec) + real(r8), intent(in) :: midpoint_pressure(:,:) ! midpoint pressure ( Pa ) + real(r8), intent(in) :: pressure_thickness(:,:) ! pressure thickness of levels (Pa) + real(r8), intent(in) :: temperature(:,:) ! temperature (K) [tfld elsewhere in CAM] + real(r8), intent(in) :: mean_mass(:,:) ! mean wet atmospheric mass ( amu ) + real(r8), target, intent(in) :: cloud_water(:,:) ! cloud liquid water content (kg/kg) + real(r8), target, intent(in) :: cloud_fraction(:,:) ! cloud fraction (unitless) + real(r8), intent(in) :: cloud_droplet_number(:,:) ! droplet number concentration (#/kg) + real(r8), intent(in) :: atmospheric_number_density(:,:) ! total atmospheric number density (/cm**3) + real(r8), intent(in) :: fixed_concentrations(:,:,:) ! fixed concentrations (/cm**3) + real(r8), target, intent(inout) :: cloud_borne_aerosol_vmr(:,:,:) ! cloud-borne aerosol (vmr) mol/mol = m3/m3 + real(r8), intent(inout) :: species_vmr(:,:,:) ! transported species (vmr) mol/mol = m3/m3 + real(r8), intent(out) :: ph_times_cloud_water(:,:) ! pH value multiplied by lwc + + real(r8), intent(out) :: aq_so4_production(:,:) ! aqueous phase production of SO4 (kg/m2/s) + real(r8), intent(out) :: aq_h2so4_production(:,:) ! aqueous phase production of H2SO4 (kg/m2/s) + real(r8), intent(out) :: aq_so4_production_from_h2o2(:) ! SO4 aqueous phase production due to H2O2 (kg/m2/s) + real(r8), intent(out) :: aq_so4_production_from_o3(:) ! SO4 aqueous phase production due to O3 (kg/m2/s) + real(r8), intent(in), optional :: specified_ph ! specified pH value. If present, this value will be used instead of calculated pH + real(r8), intent(out), optional :: aq_so4_production_from_h2o2_3d(:, :) ! 3D SO4 aqueous phase production due to H2O2 (kg/m2/s) + real(r8), intent(out), optional :: aq_so4_production_from_o3_3d(:, :) ! 3D SO4 aqueous phase production due to O3 (kg/m2/s) type(physics_state), intent(in) :: state ! Physics state variables @@ -305,28 +307,29 @@ subroutine calculate( state, & xph0 = 10._r8**(-ph0) ! initial PH value do k = 1,pver - cfact(:,k) = xhnm(:,k) & ! /cm3(a) + cfact(:,k) = atmospheric_number_density(:,k) & ! /cm3(a) * 1.e6_r8 & ! /m3(a) * 1.38e-23_r8/287._r8 & ! Kg(a)/m3(a) * 1.e-3_r8 ! Kg(a)/L(a) end do - cldconc => sox_cldaero_create_obj( cldfrc,qcw,lwc, cfact, ncol, loffset ) + cldconc => sox_cldaero_create_obj( cloud_fraction, cloud_borne_aerosol_vmr, & + cloud_water, cfact, ncol, loffset ) xso4c => cldconc%so4c xnh4c => cldconc%nh4c xno3c => cldconc%no3c xso4(:,:) = 0._r8 xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 - call so2%mixing_ratio( qin, invariants, xhnm, xso2 ) - call nh3%mixing_ratio( qin, invariants, xhnm, xnh3 ) - call hno3%mixing_ratio( qin, invariants, xhnm, xhno3 ) - call h2o2%mixing_ratio( qin, invariants, xhnm, xh2o2 ) - call o3%mixing_ratio( qin, invariants, xhnm, xo3 ) - call ho2%mixing_ratio( qin, invariants, xhnm, xho2 ) - call msa%mixing_ratio( qin, invariants, xhnm, xmsa ) - call so4%mixing_ratio( qin, invariants, xhnm, xso4 ) - call h2so4%mixing_ratio( qin, invariants, xhnm, xh2so4 ) + call so2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xso2 ) + call nh3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xnh3 ) + call hno3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xhno3 ) + call h2o2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xh2o2 ) + call o3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xo3 ) + call ho2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xho2 ) + call msa%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xmsa ) + call so4%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xso4 ) + call h2so4%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xh2so4 ) xph(:,:) = xph0 ! initial PH value @@ -336,15 +339,15 @@ subroutine calculate( state, & ver_loop0: do k = 1,pver !! pver loop for STEP 0 col_loop0: do i = 1,ncol - if (cloud_borne .and. cldfrc(i,k)>0._r8) then - xso4(i,k) = xso4c(i,k) / cldfrc(i,k) - xnh4(i,k) = xnh4c(i,k) / cldfrc(i,k) - xno3(i,k) = xno3c(i,k) / cldfrc(i,k) + if (cloud_borne .and. cloud_fraction(i,k)>0._r8) then + xso4(i,k) = xso4c(i,k) / cloud_fraction(i,k) + xnh4(i,k) = xnh4c(i,k) / cloud_fraction(i,k) + xno3(i,k) = xno3c(i,k) / cloud_fraction(i,k) endif xl = cldconc%xlwc(i,k) if( xl >= 1.e-8_r8 ) then - work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8 + work1(i) = 1._r8 / temperature(i,k) - 1._r8 / 298._r8 !----------------------------------------------------------------- ! 21-mar-2011 changes by rce @@ -360,10 +363,10 @@ subroutine calculate( state, & !----------------------------------------------------------------- !----------------------------------------------------------------- - pz = .01_r8*press(i,k) !! pressure in mb - tz = tfld(i,k) + pz = .01_r8*midpoint_pressure(i,k) !! pressure in mb + tz = temperature(i,k) patm = pz/1013._r8 - xam = press(i,k)/(1.38e-23_r8*tz) !air density /M3 + xam = midpoint_pressure(i,k)/(1.38e-23_r8*tz) !air density /M3 !----------------------------------------------------------------- ! ... hno3 @@ -453,7 +456,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! ... so4 effect !----------------------------------------------------------------- - Eso4 = xso4(i,k)*xhnm(i,k) & ! /cm3(a) + Eso4 = xso4(i,k)*atmospheric_number_density(i,k) & ! /cm3(a) *const0/xl @@ -469,7 +472,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- do iter = 1,itermax - if (.not. present(yph_in)) then + if (.not. present(specified_ph)) then if (iter == 1) then ! 1st iteration ph = lower bound value yph_lo = 2.0_r8 @@ -484,7 +487,7 @@ subroutine calculate( state, & yph = 0.5_r8*(yph_lo + yph_hi) end if else - yph = yph_in + yph = specified_ph end if ! calc current [H+] from ph @@ -594,13 +597,13 @@ subroutine calculate( state, & !============================================================== ver_loop1: do k = 1,pver col_loop1: do i = 1,ncol - work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8 - tz = tfld(i,k) + work1(i) = 1._r8 / temperature(i,k) - 1._r8 / 298._r8 + tz = temperature(i,k) xl = cldconc%xlwc(i,k) - patm = press(i,k)/101300._r8 ! press is in pascal - xam = press(i,k)/(1.38e-23_r8*tz) ! air density /M3 + patm = midpoint_pressure(i,k)/101300._r8 ! press is in pascal + xam = midpoint_pressure(i,k)/(1.38e-23_r8*tz) ! air density /M3 !----------------------------------------------------------------------- ! ... hno3 @@ -658,7 +661,7 @@ subroutine calculate( state, & endif if ( .not. cloud_borne) then ! this seems to be specific to aerosols that are not cloud borne - xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production + xh2o2(i,k) = xh2o2(i,k) + r2h2o2*time_step ! updated h2o2 by het production endif !----------------------------------------------- @@ -755,10 +758,10 @@ subroutine calculate( state, & pso4 = pso4 & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / const0 & ! [/L(a)/s] - / xhnm(i,k) + / atmospheric_number_density(i,k) - ccc = pso4*dtime + ccc = pso4*time_step ccc = max(ccc, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) @@ -803,9 +806,9 @@ subroutine calculate( state, & pso4 = pso4 & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / const0 & ! [/L(a)/s] - / xhnm(i,k) ! [mixing ratio/s] + / atmospheric_number_density(i,k) ! [mixing ratio/s] - ccc = pso4*dtime + ccc = pso4*time_step ccc = max(ccc, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) @@ -824,15 +827,15 @@ subroutine calculate( state, & end do ver_loop1 call sox_cldaero_update( state, & - pbuf, ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & - xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & - aqso4, aqh2so4, aqso4_h2o2, aqso4_o3, aqso4_h2o2_3d=aqso4_h2o2_3d, aqso4_o3_3d=aqso4_o3_3d ) + pbuf, ncol, lchnk, loffset, time_step, mean_mass, pressure_thickness, midpoint_pressure, temperature, cloud_droplet_number, cloud_fraction, cfact, cldconc%xlwc, & + xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, cloud_borne_aerosol_vmr, species_vmr, & + aq_so4_production, aq_h2so4_production, aq_so4_production_from_h2o2, aq_so4_production_from_o3, aqso4_h2o2_3d=aq_so4_production_from_h2o2_3d, aqso4_o3_3d=aq_so4_production_from_o3_3d ) - xphlwc(:,:) = 0._r8 + ph_times_cloud_water(:,:) = 0._r8 do k = 1, pver do i = 1, ncol - if (cldfrc(i,k)>=1.e-5_r8 .and. lwc(i,k)>=1.e-8_r8) then - xphlwc(i,k) = -1._r8*log10(xph(i,k)) * lwc(i,k) + if (cloud_fraction(i,k)>=1.e-5_r8 .and. cloud_water(i,k)>=1.e-8_r8) then + ph_times_cloud_water(i,k) = -1._r8*log10(xph(i,k)) * cloud_water(i,k) endif end do end do @@ -887,10 +890,10 @@ subroutine cloud_species_get_mixing_ratio( this, mixing_ratios, & fixed_concentrations, air_density, mixing_ratio ) class(cloud_species_t), intent(in) :: this - real(r8), intent(in) :: mixing_ratios(:,:,:) ! all mixing ratios (kg kg-1) [column, layer, species] - real(r8), intent(in) :: fixed_concentrations(:,:,:) ! all fixed concentrations (kg m-3) [column, layer, species] - real(r8), intent(in) :: air_density(:,:) ! air density (kg m-3) [column, layer] - real(r8), intent(out) :: mixing_ratio(:,:) ! species mixing ratio (kg kg-1) [column, layer] + real(r8), intent(in) :: mixing_ratios(:,:,:) ! all mixing ratios (mol mol-1) [column, layer, species] + real(r8), intent(in) :: fixed_concentrations(:,:,:) ! all fixed concentrations (# cm-3) [column, layer, species] + real(r8), intent(in) :: air_density(:,:) ! air density (# cm-3) [column, layer] + real(r8), intent(out) :: mixing_ratio(:,:) ! species mixing ratio (mol mol-1) [column, layer] if ( this%state_index_ == CLOUD_INDEX_UNDEFINED ) then mixing_ratio(:,:) = 0._r8 From 2254f72474a1e75b7d2b627593ea393b58c5f43b Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 15:11:54 -0800 Subject: [PATCH 03/12] update constants --- .../aerosol/cloud_aqueous_chemistry.F90 | 83 +++++++++---------- 1 file changed, 38 insertions(+), 45 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 039d49ca30..67778246e4 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -48,6 +48,11 @@ module cloud_aqueous_chemistry ! TODO: Figure out what this flag is for logical :: cloud_borne = .false. + ! Constants that should be moved to a common module + real(r8), parameter :: AVOGADRO = 6.02214129e23_r8 + real(r8), parameter :: PASCAL_TO_ATM = 1.0_r8/101325.0_r8 + real(r8), parameter :: GAS_CONSTANT_L_ATM_MOL_K = 8314.46261815324_r8*PASCAL_TO_ATM + contains !----------------------------------------------------------------------- @@ -219,24 +224,15 @@ subroutine calculate( state, & ! ! xhno3 ... in mixing ratio !----------------------------------------------------------------------- - integer, parameter :: itermax = 20 + integer, parameter :: max_iterations = 20 real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES - real(r8), parameter :: const0 = 1.e3_r8/6.023e23_r8 - real(r8), parameter :: xa0 = 11._r8 - real(r8), parameter :: xb0 = -.1_r8 - real(r8), parameter :: xa1 = 1.053_r8 - real(r8), parameter :: xb1 = -4.368_r8 - real(r8), parameter :: xa2 = 1.016_r8 - real(r8), parameter :: xb2 = -2.54_r8 - real(r8), parameter :: xa3 = .816e-32_r8 - real(r8), parameter :: xb3 = .259_r8 + real(r8), parameter :: const0 = 1.e3_r8/AVOGADRO real(r8), parameter :: kh0 = 9.e3_r8 ! HO2(g) -> Ho2(a) real(r8), parameter :: kh1 = 2.05e-5_r8 ! HO2(a) -> H+ + O2- real(r8), parameter :: kh2 = 8.6e5_r8 ! HO2(a) + ho2(a) -> h2o2(a) + o2 real(r8), parameter :: kh3 = 1.e8_r8 ! HO2(a) + o2- -> h2o2(a) + o2 - real(r8), parameter :: Ra = 8314._r8/101325._r8 ! universal constant (atm)/(M-K) - real(r8), parameter :: xkw = 1.e-14_r8 ! water acidity + real(r8), parameter :: hplus_scaling_factor = 1.e-14_r8 ! [H+] (mol/L) for pH=14 ! real(r8) :: xdelso4hp(ncol,pver) @@ -263,17 +259,14 @@ subroutine calculate( state, & real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s - real(r8), dimension(ncol,pver) :: & - xhno3, xh2o2, xso2, xso4, xno3, & - xnh3, xnh4, xo3, & - cfact, & - xph, xho2, & - xh2so4, xmsa, xso4_init, & - hehno3, & ! henry law const for hno3 - heh2o2, & ! henry law const for h2o2 - heso2, & ! henry law const for so2 - henh3, & ! henry law const for nh3 - heo3 !!, & ! henry law const for o3 + ! volume mixing ratios for cloud chemistry species + real(r8), dimension(ncol,pver) :: xhno3, xh2o2, xso2, xso4, xno3, & + xnh3, xnh4, xo3, xph, xho2, xh2so4, xmsa, xso4_init + + real(r8), dimension(ncol,pver) :: cfact ! TODO: what is this? + + ! Effective Henry's Law constants + real(r8), dimension(ncol,pver) :: hehno3, heh2o2, heso2, henh3, heo3 real(r8) :: patm_x @@ -373,7 +366,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! previous code ! hehno3(i,k) = xk*(1._r8 + xe/xph(i,k)) - ! px = hehno3(i,k) * Ra * tz * xl + ! px = hehno3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl ! hno3g = xhno3(i,k)/(1._r8 + px) ! Ehno3 = xk*xe*hno3g *patm ! equivalent new code @@ -388,7 +381,7 @@ subroutine calculate( state, & xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) ) xe = 15.4_r8 fact1_hno3 = xk*xe*patm*xhno3(i,k) - fact2_hno3 = xk*ra*tz*xl + fact2_hno3 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl fact3_hno3 = xe !----------------------------------------------------------------- @@ -396,7 +389,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! previous code ! heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k))) - ! px = heso2(i,k) * Ra * tz * xl + ! px = heso2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl ! so2g = xso2(i,k)/(1._r8+ px) ! Eso2 = xk*xe*so2g *patm ! equivalent new code @@ -412,7 +405,7 @@ subroutine calculate( state, & xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) ) x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) ) fact1_so2 = xk*xe*patm*xso2(i,k) - fact2_so2 = xk*ra*tz*xl + fact2_so2 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl fact3_so2 = xe fact4_so2 = x2 @@ -420,30 +413,30 @@ subroutine calculate( state, & ! ... nh3 !----------------------------------------------------------------- ! previous code - ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw) - ! px = henh3(i,k) * Ra * tz * xl + ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/hplus_scaling_factor) + ! px = henh3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl ! nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) - ! Enh3 = xk*xe*nh3g/xkw *patm + ! Enh3 = xk*xe*nh3g/hplus_scaling_factor *patm ! equivalent new code - ! henh3 = xk + xk*xe*hplus/xkw + ! henh3 = xk + xk*xe*hplus/hplus_scaling_factor ! nh3g = xnh34/(1 + px) ! = xnh34/(1 + henh3*ra*tz*xl) - ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw) - ! enh3 = nh3g*xk*xe*patm/xkw - ! = ((xk*xe*patm/xkw)*xnh34)/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw) + ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/hplus_scaling_factor) + ! enh3 = nh3g*xk*xe*patm/hplus_scaling_factor + ! = ((xk*xe*patm/hplus_scaling_factor)*xnh34)/(1 + xk*ra*tz*xl*(1 + xe*hplus/hplus_scaling_factor) ! = ( fact1_nh3 )/(1 + fact2_nh3 *(1 + fact3_nh3*hplus) ! [nh4+] = enh3*hplus xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) ) - fact1_nh3 = (xk*xe*patm/xkw)*(xnh3(i,k)+xnh4(i,k)) - fact2_nh3 = xk*ra*tz*xl - fact3_nh3 = xe/xkw + fact1_nh3 = (xk*xe*patm/hplus_scaling_factor)*(xnh3(i,k)+xnh4(i,k)) + fact2_nh3 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl + fact3_nh3 = xe/hplus_scaling_factor !----------------------------------------------------------------- ! ... h2o effects !----------------------------------------------------------------- - Eh2o = xkw + Eh2o = hplus_scaling_factor !----------------------------------------------------------------- ! ... co2 effects @@ -470,7 +463,7 @@ subroutine calculate( state, & ! yposnet_lo and yposnet_hi = net positive ions for ! yph_lo and yph_hi !----------------------------------------------------------------- - do iter = 1,itermax + do iter = 1, max_iterations if (.not. present(specified_ph)) then if (iter == 1) then @@ -634,7 +627,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) ) - henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw) + henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/hplus_scaling_factor) !----------------------------------------------------------------- ! ... o3 @@ -671,31 +664,31 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------- - px = hehno3(i,k) * Ra * tz * xl + px = hehno3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl hno3g(i,k) = (xhno3(i,k)+xno3(i,k))/(1._r8 + px) !------------------------------------------------------------------------ ! ... h2o2 !------------------------------------------------------------------------ - px = heh2o2(i,k) * Ra * tz * xl + px = heh2o2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl h2o2g = xh2o2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... so2 !------------------------------------------------------------------------ - px = heso2(i,k) * Ra * tz * xl + px = heso2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl so2g = xso2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... o3 !------------------------------------------------------------------------ - px = heo3(i,k) * Ra * tz * xl + px = heo3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl o3g = xo3(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... nh3 !------------------------------------------------------------------------ - px = henh3(i,k) * Ra * tz * xl + px = henh3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl if (nh3%exists()) then nh3g(i,k) = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) else From 9d1f90f4cbde219256db963ba4970e481190f278 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 15:57:35 -0800 Subject: [PATCH 04/12] renaming --- .../aerosol/cloud_aqueous_chemistry.F90 | 57 +++++++++---------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 67778246e4..0832934e6b 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -49,9 +49,10 @@ module cloud_aqueous_chemistry logical :: cloud_borne = .false. ! Constants that should be moved to a common module - real(r8), parameter :: AVOGADRO = 6.02214129e23_r8 + real(r8), parameter :: AVOGADRO = 6.02214129e23_r8 ! mol-1 real(r8), parameter :: PASCAL_TO_ATM = 1.0_r8/101325.0_r8 real(r8), parameter :: GAS_CONSTANT_L_ATM_MOL_K = 8314.46261815324_r8*PASCAL_TO_ATM + real(r8), parameter :: BOLTZMANN = 1.380649e-23_r8 ! J K-1 contains @@ -202,10 +203,10 @@ subroutine calculate( state, & real(r8), target, intent(in) :: cloud_fraction(:,:) ! cloud fraction (unitless) real(r8), intent(in) :: cloud_droplet_number(:,:) ! droplet number concentration (#/kg) real(r8), intent(in) :: atmospheric_number_density(:,:) ! total atmospheric number density (/cm**3) - real(r8), intent(in) :: fixed_concentrations(:,:,:) ! fixed concentrations (/cm**3) - real(r8), target, intent(inout) :: cloud_borne_aerosol_vmr(:,:,:) ! cloud-borne aerosol (vmr) mol/mol = m3/m3 - real(r8), intent(inout) :: species_vmr(:,:,:) ! transported species (vmr) mol/mol = m3/m3 - real(r8), intent(out) :: ph_times_cloud_water(:,:) ! pH value multiplied by lwc + real(r8), intent(in) :: fixed_concentrations(:,:,:) ! fixed concentrations (/cm**3) [invariants elsewhere in CAM] + real(r8), target, intent(inout) :: cloud_borne_aerosol_vmr(:,:,:) ! cloud-borne aerosol (vmr) mol/mol = m3/m3 [qcw elsewhere in CAM] + real(r8), intent(inout) :: species_vmr(:,:,:) ! transported species (vmr) mol/mol = m3/m3 [qin elsewhere in CAM] + real(r8), intent(out) :: ph_times_cloud_water(:,:) ! pH value multiplied by cloud liquid water content real(r8), intent(out) :: aq_so4_production(:,:) ! aqueous phase production of SO4 (kg/m2/s) real(r8), intent(out) :: aq_h2so4_production(:,:) ! aqueous phase production of H2SO4 (kg/m2/s) @@ -234,13 +235,13 @@ subroutine calculate( state, & real(r8), parameter :: kh3 = 1.e8_r8 ! HO2(a) + o2- -> h2o2(a) + o2 real(r8), parameter :: hplus_scaling_factor = 1.e-14_r8 ! [H+] (mol/L) for pH=14 - ! - real(r8) :: xdelso4hp(ncol,pver) + ! Change in aqueous sulfate volume mixing ratio over current time step (mol mol-1) + real(r8) :: change_in_aq_so4_mixing_ratio(ncol,pver) integer :: k, i, iter, file real(r8) :: wrk, delta real(r8) :: xph0, aden, xk, xe, x2 - real(r8) :: tz, xl, px, qz, pz, es, qs, patm + real(r8) :: xl, px, qz, pz, es, qs, patm real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3 real(r8) :: so2g, h2o2g, co2g, o3g real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a @@ -254,7 +255,6 @@ subroutine calculate( state, & ! schwartz JGR, 1984, 11589 !----------------------------------------------------------------------- real(r8) :: kh4 ! kh2+kh3 - real(r8) :: xam ! air density /cm3 real(r8) :: ho2s ! ho2s = ho2(a)+o2- real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s @@ -357,9 +357,8 @@ subroutine calculate( state, & !----------------------------------------------------------------- pz = .01_r8*midpoint_pressure(i,k) !! pressure in mb - tz = temperature(i,k) - patm = pz/1013._r8 - xam = midpoint_pressure(i,k)/(1.38e-23_r8*tz) !air density /M3 + ! This should be divided by 101325, not 101300, but fixing this breaks the tests + patm = midpoint_pressure(i,k)/101300._r8 !----------------------------------------------------------------- ! ... hno3 @@ -381,7 +380,7 @@ subroutine calculate( state, & xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) ) xe = 15.4_r8 fact1_hno3 = xk*xe*patm*xhno3(i,k) - fact2_hno3 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl + fact2_hno3 = xk*GAS_CONSTANT_L_ATM_MOL_K*temperature(i,k)*xl fact3_hno3 = xe !----------------------------------------------------------------- @@ -405,7 +404,7 @@ subroutine calculate( state, & xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) ) x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) ) fact1_so2 = xk*xe*patm*xso2(i,k) - fact2_so2 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl + fact2_so2 = xk*GAS_CONSTANT_L_ATM_MOL_K*temperature(i,k)*xl fact3_so2 = xe fact4_so2 = x2 @@ -430,7 +429,7 @@ subroutine calculate( state, & xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) ) fact1_nh3 = (xk*xe*patm/hplus_scaling_factor)*(xnh3(i,k)+xnh4(i,k)) - fact2_nh3 = xk*GAS_CONSTANT_L_ATM_MOL_K*tz*xl + fact2_nh3 = xk*GAS_CONSTANT_L_ATM_MOL_K*temperature(i,k)*xl fact3_nh3 = xe/hplus_scaling_factor !----------------------------------------------------------------- @@ -591,12 +590,10 @@ subroutine calculate( state, & ver_loop1: do k = 1,pver col_loop1: do i = 1,ncol work1(i) = 1._r8 / temperature(i,k) - 1._r8 / 298._r8 - tz = temperature(i,k) - xl = cldconc%xlwc(i,k) - patm = midpoint_pressure(i,k)/101300._r8 ! press is in pascal - xam = midpoint_pressure(i,k)/(1.38e-23_r8*tz) ! air density /M3 + ! This should be dividing by 101325, not 101300, but changing it breaks the tests + patm = midpoint_pressure(i,k) / 101300._r8 ! press is in pascal !----------------------------------------------------------------------- ! ... hno3 @@ -646,11 +643,11 @@ subroutine calculate( state, & if ( cloud_borne ) then r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s / const0*1.e+6_r8 & ! correct a bug here ???? - / xam + / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) else r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s * const0 & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s - / xam ! /cm3(a)/s / air-den = mix-ratio/s + / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) ! /cm3(a)/s / air-den = mix-ratio/s endif if ( .not. cloud_borne) then ! this seems to be specific to aerosols that are not cloud borne @@ -664,31 +661,31 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------- - px = hehno3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl + px = hehno3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * temperature(i,k) * xl hno3g(i,k) = (xhno3(i,k)+xno3(i,k))/(1._r8 + px) !------------------------------------------------------------------------ ! ... h2o2 !------------------------------------------------------------------------ - px = heh2o2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl + px = heh2o2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * temperature(i,k) * xl h2o2g = xh2o2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... so2 !------------------------------------------------------------------------ - px = heso2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl + px = heso2(i,k) * GAS_CONSTANT_L_ATM_MOL_K * temperature(i,k) * xl so2g = xso2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... o3 !------------------------------------------------------------------------ - px = heo3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl + px = heo3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * temperature(i,k) * xl o3g = xo3(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... nh3 !------------------------------------------------------------------------ - px = henh3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl + px = henh3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * temperature(i,k) * xl if (nh3%exists()) then nh3g(i,k) = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) else @@ -710,8 +707,8 @@ subroutine calculate( state, & !------------------------------------------------------------------------ ! ... S(IV)+ O3 !------------------------------------------------------------------------ - rao3 = 4.39e11_r8 * EXP(-4131._r8/tz) & - + 2.56e3_r8 * EXP(-996._r8 /tz) /xph(i,k) + rao3 = 4.39e11_r8 * EXP(-4131._r8/temperature(i,k)) & + + 2.56e3_r8 * EXP(-996._r8 /temperature(i,k)) /xph(i,k) !----------------------------------------------------------------- ! ... Prediction after aqueous phase @@ -788,7 +785,7 @@ subroutine calculate( state, & END IF if (cloud_borne) then - xdelso4hp(i,k) = xso4(i,k) - xso4_init(i,k) + change_in_aq_so4_mixing_ratio(i,k) = xso4(i,k) - xso4_init(i,k) endif !........................... ! S(IV) + O3 = S(VI) @@ -821,7 +818,7 @@ subroutine calculate( state, & call sox_cldaero_update( state, & pbuf, ncol, lchnk, loffset, time_step, mean_mass, pressure_thickness, midpoint_pressure, temperature, cloud_droplet_number, cloud_fraction, cfact, cldconc%xlwc, & - xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, cloud_borne_aerosol_vmr, species_vmr, & + change_in_aq_so4_mixing_ratio, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, cloud_borne_aerosol_vmr, species_vmr, & aq_so4_production, aq_h2so4_production, aq_so4_production_from_h2o2, aq_so4_production_from_o3, aqso4_h2o2_3d=aq_so4_production_from_h2o2_3d, aqso4_o3_3d=aq_so4_production_from_o3_3d ) ph_times_cloud_water(:,:) = 0._r8 From 755e4063d948b3c49cb62a8887b954376ac6e7b3 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 16:29:11 -0800 Subject: [PATCH 05/12] remove unused variables --- .../aerosol/cloud_aqueous_chemistry.F90 | 90 +++++++++---------- 1 file changed, 41 insertions(+), 49 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 0832934e6b..159f0aa1e3 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -53,6 +53,7 @@ module cloud_aqueous_chemistry real(r8), parameter :: PASCAL_TO_ATM = 1.0_r8/101325.0_r8 real(r8), parameter :: GAS_CONSTANT_L_ATM_MOL_K = 8314.46261815324_r8*PASCAL_TO_ATM real(r8), parameter :: BOLTZMANN = 1.380649e-23_r8 ! J K-1 + real(r8), parameter :: GAS_CONSTANT_DRY_AIR_J_KG_K = 287.0_r8 ! J kg-1 K-1 contains @@ -143,7 +144,7 @@ subroutine calculate( state, & cloud_water, & cloud_fraction, & cloud_droplet_number, & - atmospheric_number_density, & + air_number_density, & fixed_concentrations, & cloud_borne_aerosol_vmr, & species_vmr, & @@ -202,7 +203,7 @@ subroutine calculate( state, & real(r8), target, intent(in) :: cloud_water(:,:) ! cloud liquid water content (kg/kg) real(r8), target, intent(in) :: cloud_fraction(:,:) ! cloud fraction (unitless) real(r8), intent(in) :: cloud_droplet_number(:,:) ! droplet number concentration (#/kg) - real(r8), intent(in) :: atmospheric_number_density(:,:) ! total atmospheric number density (/cm**3) + real(r8), intent(in) :: air_number_density(:,:) ! total atmospheric number density (/cm**3) real(r8), intent(in) :: fixed_concentrations(:,:,:) ! fixed concentrations (/cm**3) [invariants elsewhere in CAM] real(r8), target, intent(inout) :: cloud_borne_aerosol_vmr(:,:,:) ! cloud-borne aerosol (vmr) mol/mol = m3/m3 [qcw elsewhere in CAM] real(r8), intent(inout) :: species_vmr(:,:,:) ! transported species (vmr) mol/mol = m3/m3 [qin elsewhere in CAM] @@ -227,7 +228,6 @@ subroutine calculate( state, & !----------------------------------------------------------------------- integer, parameter :: max_iterations = 20 real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES - real(r8), parameter :: const0 = 1.e3_r8/AVOGADRO real(r8), parameter :: kh0 = 9.e3_r8 ! HO2(g) -> Ho2(a) real(r8), parameter :: kh1 = 2.05e-5_r8 ! HO2(a) -> H+ + O2- @@ -238,15 +238,12 @@ subroutine calculate( state, & ! Change in aqueous sulfate volume mixing ratio over current time step (mol mol-1) real(r8) :: change_in_aq_so4_mixing_ratio(ncol,pver) - integer :: k, i, iter, file - real(r8) :: wrk, delta - real(r8) :: xph0, aden, xk, xe, x2 - real(r8) :: xl, px, qz, pz, es, qs, patm + integer :: k, i, iter + real(r8) :: xk, xe, x2 + real(r8) :: xl, px, patm real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3 real(r8) :: so2g, h2o2g, co2g, o3g - real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a real(r8) :: rah2o2, rao3, pso4, ccc - real(r8) :: cnh3, chno3, com, com1, com2, xra real(r8) :: hno3g(ncol,pver), nh3g(ncol,pver) ! @@ -263,7 +260,7 @@ subroutine calculate( state, & real(r8), dimension(ncol,pver) :: xhno3, xh2o2, xso2, xso4, xno3, & xnh3, xnh4, xo3, xph, xho2, xh2so4, xmsa, xso4_init - real(r8), dimension(ncol,pver) :: cfact ! TODO: what is this? + real(r8), dimension(ncol,pver) :: air_mass_density_kg_l ! kg L-1 ! Effective Henry's Law constants real(r8), dimension(ncol,pver) :: hehno3, heh2o2, heso2, henh3, heo3 @@ -297,34 +294,32 @@ subroutine calculate( state, & ! ... Initial values ! The values of so2, so4 are after (1) SLT, and CHEM !----------------------------------------------------------------- - xph0 = 10._r8**(-ph0) ! initial PH value - do k = 1,pver - cfact(:,k) = atmospheric_number_density(:,k) & ! /cm3(a) - * 1.e6_r8 & ! /m3(a) - * 1.38e-23_r8/287._r8 & ! Kg(a)/m3(a) - * 1.e-3_r8 ! Kg(a)/L(a) + air_mass_density_kg_l(:,k) = & + air_number_density(:,k) & ! /cm3(a) + * 1.e3_r8 & ! /L(a) + * BOLTZMANN/GAS_CONSTANT_DRY_AIR_J_KG_K ! Kg(a)/L(a) end do cldconc => sox_cldaero_create_obj( cloud_fraction, cloud_borne_aerosol_vmr, & - cloud_water, cfact, ncol, loffset ) + cloud_water, air_mass_density_kg_l, ncol, loffset ) xso4c => cldconc%so4c xnh4c => cldconc%nh4c xno3c => cldconc%no3c xso4(:,:) = 0._r8 xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 - call so2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xso2 ) - call nh3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xnh3 ) - call hno3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xhno3 ) - call h2o2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xh2o2 ) - call o3%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xo3 ) - call ho2%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xho2 ) - call msa%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xmsa ) - call so4%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xso4 ) - call h2so4%mixing_ratio( species_vmr, fixed_concentrations, atmospheric_number_density, xh2so4 ) - - xph(:,:) = xph0 ! initial PH value + call so2%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xso2 ) + call nh3%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xnh3 ) + call hno3%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xhno3 ) + call h2o2%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xh2o2 ) + call o3%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xo3 ) + call ho2%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xho2 ) + call msa%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xmsa ) + call so4%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xso4 ) + call h2so4%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xh2so4 ) + + xph(:,:) = 10._r8**(-ph0) ! initial PH value !----------------------------------------------------------------- ! ... Temperature dependent Henry constants @@ -356,7 +351,6 @@ subroutine calculate( state, & !----------------------------------------------------------------- !----------------------------------------------------------------- - pz = .01_r8*midpoint_pressure(i,k) !! pressure in mb ! This should be divided by 101325, not 101300, but fixing this breaks the tests patm = midpoint_pressure(i,k)/101300._r8 @@ -448,8 +442,8 @@ subroutine calculate( state, & !----------------------------------------------------------------- ! ... so4 effect !----------------------------------------------------------------- - Eso4 = xso4(i,k)*atmospheric_number_density(i,k) & ! /cm3(a) - *const0/xl + Eso4 = xso4(i,k)*air_number_density(i,k) & ! /cm3(a) + * (1.e3_r8/AVOGADRO) / xl !----------------------------------------------------------------- @@ -575,8 +569,7 @@ subroutine calculate( state, & end do ! iter if( .not. converged ) then - write(iulog,*) 'setsox: pH failed to converge @ (',i,',',k,'), % change=', & - 100._r8*delta + write(iulog,*) 'setsox: pH failed to converge @ (',i,',',k,')' end if else xph(i,k) = 1.e-7_r8 @@ -616,8 +609,7 @@ subroutine calculate( state, & xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) ) x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) ) - wrk = xe/xph(i,k) - heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k))) + heso2(i,k) = xk*(1._r8 + xe/xph(i,k)*(1._r8 + x2/xph(i,k))) !----------------------------------------------------------------- ! ... nh3 @@ -641,12 +633,12 @@ subroutine calculate( state, & r1h2o2 = kh4*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s if ( cloud_borne ) then - r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s - / const0*1.e+6_r8 & ! correct a bug here ???? + r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s + / (1.e3_r8/AVOGADRO)*1.e+6_r8 & ! correct a bug here ???? / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) else - r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s - * const0 & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s + r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s + * (1.e3_r8/AVOGADRO) & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) ! /cm3(a)/s / air-den = mix-ratio/s endif @@ -747,8 +739,8 @@ subroutine calculate( state, & pso4 = pso4 & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] - / const0 & ! [/L(a)/s] - / atmospheric_number_density(i,k) + / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] + / air_number_density(i,k) ccc = pso4*time_step @@ -793,10 +785,10 @@ subroutine calculate( state, & pso4 = rao3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s] - pso4 = pso4 & ! [M/s] = [mole/L(w)/s] - * xl & ! [mole/L(a)/s] - / const0 & ! [/L(a)/s] - / atmospheric_number_density(i,k) ! [mixing ratio/s] + pso4 = pso4 & ! [M/s] = [mole/L(w)/s] + * xl & ! [mole/L(a)/s] + / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] + / air_number_density(i,k) ! [mixing ratio/s] ccc = pso4*time_step ccc = max(ccc, 1.e-30_r8) @@ -817,7 +809,7 @@ subroutine calculate( state, & end do ver_loop1 call sox_cldaero_update( state, & - pbuf, ncol, lchnk, loffset, time_step, mean_mass, pressure_thickness, midpoint_pressure, temperature, cloud_droplet_number, cloud_fraction, cfact, cldconc%xlwc, & + pbuf, ncol, lchnk, loffset, time_step, mean_mass, pressure_thickness, midpoint_pressure, temperature, cloud_droplet_number, cloud_fraction, air_mass_density_kg_l, cldconc%xlwc, & change_in_aq_so4_mixing_ratio, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, cloud_borne_aerosol_vmr, species_vmr, & aq_so4_production, aq_h2so4_production, aq_so4_production_from_h2o2, aq_so4_production_from_o3, aqso4_h2o2_3d=aq_so4_production_from_h2o2_3d, aqso4_o3_3d=aq_so4_production_from_o3_3d ) @@ -877,12 +869,12 @@ end function cloud_species_exists ! ! For species that are not defined, the mixing ratio is set to zero. subroutine cloud_species_get_mixing_ratio( this, mixing_ratios, & - fixed_concentrations, air_density, mixing_ratio ) + fixed_concentrations, air_number_density, mixing_ratio ) class(cloud_species_t), intent(in) :: this real(r8), intent(in) :: mixing_ratios(:,:,:) ! all mixing ratios (mol mol-1) [column, layer, species] real(r8), intent(in) :: fixed_concentrations(:,:,:) ! all fixed concentrations (# cm-3) [column, layer, species] - real(r8), intent(in) :: air_density(:,:) ! air density (# cm-3) [column, layer] + real(r8), intent(in) :: air_number_density(:,:) ! air density (# cm-3) [column, layer] real(r8), intent(out) :: mixing_ratio(:,:) ! species mixing ratio (mol mol-1) [column, layer] if ( this%state_index_ == CLOUD_INDEX_UNDEFINED ) then @@ -891,7 +883,7 @@ subroutine cloud_species_get_mixing_ratio( this, mixing_ratios, & end if if ( this%is_constant_ ) then mixing_ratio(:,:) = fixed_concentrations(:,:,this%state_index_) & - / air_density(:,:) + / air_number_density(:,:) else mixing_ratio(:,:) = mixing_ratios(:,:,this%state_index_) end if From a84989f89ca58d3d1a0ff0f4500f8eca89fb84e7 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 16:45:17 -0800 Subject: [PATCH 06/12] remove unused dependencies --- .../aerosol/cloud_aqueous_chemistry.F90 | 21 +++++++------------ 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 159f0aa1e3..133e358430 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -14,9 +14,9 @@ module cloud_aqueous_chemistry #define USE_MAM #undef USE_CARMA - use shr_kind_mod, only : r8 => shr_kind_r8 - use cam_logfile, only : iulog - use physics_buffer,only: physics_buffer_desc, pbuf_get_index, pbuf_add_field, dtype_r8 + use shr_kind_mod, only : r8 => shr_kind_r8 + use cam_logfile, only : iulog + use physics_buffer, only: physics_buffer_desc use physics_types, only: physics_state implicit none @@ -64,15 +64,14 @@ subroutine initialize() ! ... initialize the hetero sox routine !----------------------------------------------------------------------- - use mo_chem_utls, only : get_spc_ndx, get_inv_ndx - use spmd_utils, only : masterproc - use phys_control, only : phys_getopts + use spmd_utils, only : masterproc + use phys_control, only : phys_getopts use carma_flags_mod, only : carma_do_cloudborne #ifdef USE_MAM - use mam_clouds, only : sox_cldaero_init + use mam_clouds, only : sox_cldaero_init #endif #ifdef USE_CARMA - use carma_clouds, only : sox_cldaero_init + use carma_clouds, only : sox_cldaero_init #endif logical :: modal_aerosols @@ -173,11 +172,7 @@ subroutine calculate( state, & ! (d) PREDICTION !----------------------------------------------------------------------- ! - use ppgrid, only : pcols, pver - use chem_mods, only : gas_pcnst, nfs - use chem_mods, only : adv_mass - use physconst, only : mwdry, gravit - use mo_constants, only : pi + use ppgrid, only : pver #ifdef USE_MAM use mam_clouds, only : sox_cldaero_update, sox_cldaero_create_obj, & sox_cldaero_destroy_obj From 388ec25e74731179663dbf529d2171dfd0ec55f0 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Fri, 21 Feb 2025 17:05:35 -0800 Subject: [PATCH 07/12] renaming --- .../aerosol/cloud_aqueous_chemistry.F90 | 51 ++++++++++--------- 1 file changed, 26 insertions(+), 25 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index 133e358430..adae0a76e5 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -238,7 +238,10 @@ subroutine calculate( state, & real(r8) :: xl, px, patm real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3 real(r8) :: so2g, h2o2g, co2g, o3g - real(r8) :: rah2o2, rao3, pso4, ccc + real(r8) :: k_siv_h2o2 ! rate constant for reaction of S(IV) with H2O2 + real(r8) :: k_siv_o3 ! rate constant for reaction of S(IV) with O3 + real(r8) :: dso4_dt ! rate of change of SO4 + real(r8) :: delta_concentration real(r8) :: hno3g(ncol,pver), nh3g(ncol,pver) ! @@ -688,13 +691,13 @@ subroutine calculate( state, & !------------------------------------------------------------------------ ! ... S(IV) (HSO3) + H2O2 !------------------------------------------------------------------------ - rah2o2 = 8.e4_r8 * EXP( -3650._r8*work1(i) ) & + k_siv_h2o2 = 8.e4_r8 * EXP( -3650._r8*work1(i) ) & / (.1_r8 + xph(i,k)) !------------------------------------------------------------------------ ! ... S(IV)+ O3 !------------------------------------------------------------------------ - rao3 = 4.39e11_r8 * EXP(-4131._r8/temperature(i,k)) & + k_siv_o3 = 4.39e11_r8 * EXP(-4131._r8/temperature(i,k)) & + 2.56e3_r8 * EXP(-996._r8 /temperature(i,k)) /xph(i,k) !----------------------------------------------------------------- @@ -723,28 +726,26 @@ subroutine calculate( state, & endif if (cloud_borne) then - - pso4 = rah2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm_x & + dso4_dt = k_siv_h2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm_x & * 1.23_r8 *EXP(3120._r8*work1(i)) * so2g * patm_x else - pso4 = rah2o2 * heh2o2(i,k) * h2o2g * patm_x & + dso4_dt = k_siv_h2o2 * heh2o2(i,k) * h2o2g * patm_x & * heso2(i,k) * so2g * patm_x ! [M/s] - endif - pso4 = pso4 & ! [M/s] = [mole/L(w)/s] + dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] / air_number_density(i,k) - ccc = pso4*time_step - ccc = max(ccc, 1.e-30_r8) + delta_concentration = dso4_dt*time_step + delta_concentration = max(delta_concentration, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) IF (xh2o2(i,k) .gt. xso2(i,k)) THEN - if (ccc .gt. xso2(i,k)) then + if (delta_concentration .gt. xso2(i,k)) then xso4(i,k)=xso4(i,k)+xso2(i,k) if (cloud_borne) then xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) @@ -754,20 +755,20 @@ subroutine calculate( state, & xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) endif else - xso4(i,k) = xso4(i,k) + ccc - xh2o2(i,k) = xh2o2(i,k) - ccc - xso2(i,k) = xso2(i,k) - ccc + xso4(i,k) = xso4(i,k) + delta_concentration + xh2o2(i,k) = xh2o2(i,k) - delta_concentration + xso2(i,k) = xso2(i,k) - delta_concentration end if ELSE - if (ccc .gt. xh2o2(i,k)) then + if (delta_concentration .gt. xh2o2(i,k)) then xso4(i,k)=xso4(i,k)+xh2o2(i,k) xso2(i,k)=xso2(i,k)-xh2o2(i,k) xh2o2(i,k)=1.e-20_r8 else - xso4(i,k) = xso4(i,k) + ccc - xh2o2(i,k) = xh2o2(i,k) - ccc - xso2(i,k) = xso2(i,k) - ccc + xso4(i,k) = xso4(i,k) + delta_concentration + xh2o2(i,k) = xh2o2(i,k) - delta_concentration + xso2(i,k) = xso2(i,k) - delta_concentration end if END IF @@ -778,24 +779,24 @@ subroutine calculate( state, & ! S(IV) + O3 = S(VI) !........................... - pso4 = rao3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s] + dso4_dt = k_siv_o3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s] - pso4 = pso4 & ! [M/s] = [mole/L(w)/s] + dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] / air_number_density(i,k) ! [mixing ratio/s] - ccc = pso4*time_step - ccc = max(ccc, 1.e-30_r8) + delta_concentration = dso4_dt*time_step + delta_concentration = max(delta_concentration, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) - if (ccc .gt. xso2(i,k)) then + if (delta_concentration .gt. xso2(i,k)) then xso4(i,k) = xso4(i,k) + xso2(i,k) xso2(i,k) = 1.e-20_r8 else - xso4(i,k) = xso4(i,k) + ccc - xso2(i,k) = xso2(i,k) - ccc + xso4(i,k) = xso4(i,k) + delta_concentration + xso2(i,k) = xso2(i,k) - delta_concentration end if END IF !! WHEN CLOUD IS PRESENTED From 0635d424f0fed951a3296e4d546656e44c513a57 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Mon, 24 Feb 2025 09:42:34 -0800 Subject: [PATCH 08/12] more clean-up --- .../aerosol/cloud_aqueous_chemistry.F90 | 70 +++++++++++-------- 1 file changed, 39 insertions(+), 31 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index adae0a76e5..d1eb34e31b 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -32,7 +32,7 @@ module cloud_aqueous_chemistry ! Cloud chemistry species information type :: cloud_species_t logical :: is_constant_ = .false. - integer :: state_index_ = CLOUD_INDEX_UNDEFINED + integer :: state_index_ = CLOUD_INDEX_UNDEFINED ! index in the state vector or the fixed concentrations array character(len=:), allocatable :: name_ contains procedure :: exists => cloud_species_exists @@ -61,7 +61,9 @@ module cloud_aqueous_chemistry !----------------------------------------------------------------------- subroutine initialize() !----------------------------------------------------------------------- - ! ... initialize the hetero sox routine + ! ... prepare for cloud aqueous chemistry + ! - Look up cloud chemistry species + ! - Determine if enough species are present to perform cloud chemistry !----------------------------------------------------------------------- use spmd_utils, only : masterproc @@ -129,6 +131,10 @@ subroutine initialize() end subroutine initialize !----------------------------------------------------------------------- +! Calculates the formation of sulfate and updates sulfate concentrations +! due to cloud aqueous chemistry. Also outputs the production rates +! (kg m-2 s-1) of various reactions of sulfur species with various +! oxidants (e.g., H2O2(aq), H2SO4(aq)). !----------------------------------------------------------------------- subroutine calculate( state, & pbuf, & @@ -191,10 +197,10 @@ subroutine calculate( state, & integer, intent(in) :: lchnk ! chunk id integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array real(r8), intent(in) :: time_step ! time step (sec) - real(r8), intent(in) :: midpoint_pressure(:,:) ! midpoint pressure ( Pa ) - real(r8), intent(in) :: pressure_thickness(:,:) ! pressure thickness of levels (Pa) + real(r8), intent(in) :: midpoint_pressure(:,:) ! midpoint pressure (Pa) + real(r8), intent(in) :: pressure_thickness(:,:) ! pressure thickness of levels (Pa) [pdel elsewhere in CAM] real(r8), intent(in) :: temperature(:,:) ! temperature (K) [tfld elsewhere in CAM] - real(r8), intent(in) :: mean_mass(:,:) ! mean wet atmospheric mass ( amu ) + real(r8), intent(in) :: mean_mass(:,:) ! mean wet atmospheric mass (amu) real(r8), target, intent(in) :: cloud_water(:,:) ! cloud liquid water content (kg/kg) real(r8), target, intent(in) :: cloud_fraction(:,:) ! cloud fraction (unitless) real(r8), intent(in) :: cloud_droplet_number(:,:) ! droplet number concentration (#/kg) @@ -218,21 +224,25 @@ subroutine calculate( state, & !----------------------------------------------------------------------- ! ... Local variables - ! - ! xhno3 ... in mixing ratio !----------------------------------------------------------------------- - integer, parameter :: max_iterations = 20 - real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES + integer, parameter :: MAX_ITERATIONS = 20 + real(r8), parameter :: INITIAL_PH = 5.0_r8 ! Initial pH value + real(r8), parameter :: MINIMUM_CLOUD_LIQUID_WATER = 1.e-8_r8 ! Minimum cloud liquid water content (kg/kg) + ! Effective Henry's Law constants for HO2 partitioning + ! TODO: skipping remnaming of these in anticipation of a partitioning struct real(r8), parameter :: kh0 = 9.e3_r8 ! HO2(g) -> Ho2(a) real(r8), parameter :: kh1 = 2.05e-5_r8 ! HO2(a) -> H+ + O2- real(r8), parameter :: kh2 = 8.6e5_r8 ! HO2(a) + ho2(a) -> h2o2(a) + o2 real(r8), parameter :: kh3 = 1.e8_r8 ! HO2(a) + o2- -> h2o2(a) + o2 - real(r8), parameter :: hplus_scaling_factor = 1.e-14_r8 ! [H+] (mol/L) for pH=14 + + ! Water dissociation constant (mol2/L2) [H+][OH-] + real(r8), parameter :: water_dissociation_constant = 1.e-14_r8 ! Change in aqueous sulfate volume mixing ratio over current time step (mol mol-1) real(r8) :: change_in_aq_so4_mixing_ratio(ncol,pver) + ! TODO: Skipping renaming of these in anticipation of partitioning/pH estimation structs integer :: k, i, iter real(r8) :: xk, xe, x2 real(r8) :: xl, px, patm @@ -249,7 +259,6 @@ subroutine calculate( state, & ! for Ho2(g) -> H2o2(a) formation ! schwartz JGR, 1984, 11589 !----------------------------------------------------------------------- - real(r8) :: kh4 ! kh2+kh3 real(r8) :: ho2s ! ho2s = ho2(a)+o2- real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s @@ -263,16 +272,19 @@ subroutine calculate( state, & ! Effective Henry's Law constants real(r8), dimension(ncol,pver) :: hehno3, heh2o2, heso2, henh3, heo3 + ! TODO: Figure out what this actually represents (it differs based on the value of cloud_borne) real(r8) :: patm_x real(r8), dimension(ncol) :: work1 logical :: converged + ! Cloud-borne species volume mixing ratios real(r8), pointer :: xso4c(:,:) real(r8), pointer :: xnh4c(:,:) real(r8), pointer :: xno3c(:,:) type(cldaero_conc_t), pointer :: cldconc + ! TODO: Skipping renaming of these in anticipation of partitioning/pH estimation structs real(r8) :: fact1_hno3, fact2_hno3, fact3_hno3 real(r8) :: fact1_so2, fact2_so2, fact3_so2, fact4_so2 real(r8) :: fact1_nh3, fact2_nh3, fact3_nh3 @@ -282,10 +294,6 @@ subroutine calculate( state, & real(r8) :: yph, yph_lo, yph_hi real(r8) :: ynetpos, ynetpos_lo, ynetpos_hi - !----------------------------------------------------------------- - ! ... NOTE: The press array is in pascals and must be - ! mutiplied by 10 to yield dynes/cm**2. - !----------------------------------------------------------------- !================================================================== ! ... First set the PH !================================================================== @@ -317,7 +325,7 @@ subroutine calculate( state, & call so4%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xso4 ) call h2so4%mixing_ratio( species_vmr, fixed_concentrations, air_number_density, xh2so4 ) - xph(:,:) = 10._r8**(-ph0) ! initial PH value + xph(:,:) = 10._r8**(-INITIAL_PH) ! initial PH value !----------------------------------------------------------------- ! ... Temperature dependent Henry constants @@ -332,7 +340,7 @@ subroutine calculate( state, & endif xl = cldconc%xlwc(i,k) - if( xl >= 1.e-8_r8 ) then + if( xl >= MINIMUM_CLOUD_LIQUID_WATER ) then work1(i) = 1._r8 / temperature(i,k) - 1._r8 / 298._r8 !----------------------------------------------------------------- @@ -404,30 +412,31 @@ subroutine calculate( state, & ! ... nh3 !----------------------------------------------------------------- ! previous code - ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/hplus_scaling_factor) + ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/water_dissociation_constant) ! px = henh3(i,k) * GAS_CONSTANT_L_ATM_MOL_K * tz * xl ! nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) - ! Enh3 = xk*xe*nh3g/hplus_scaling_factor *patm + ! Enh3 = xk*xe*nh3g/water_dissociation_constant *patm ! equivalent new code - ! henh3 = xk + xk*xe*hplus/hplus_scaling_factor + ! henh3 = xk + xk*xe*hplus/water_dissociation_constant ! nh3g = xnh34/(1 + px) ! = xnh34/(1 + henh3*ra*tz*xl) - ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/hplus_scaling_factor) - ! enh3 = nh3g*xk*xe*patm/hplus_scaling_factor - ! = ((xk*xe*patm/hplus_scaling_factor)*xnh34)/(1 + xk*ra*tz*xl*(1 + xe*hplus/hplus_scaling_factor) + ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/water_dissociation_constant) + ! enh3 = nh3g*xk*xe*patm/water_dissociation_constant + ! = ((xk*xe*patm/water_dissociation_constant)*xnh34)/ + ! (1 + xk*ra*tz*xl*(1 + xe*hplus/water_dissociation_constant) ! = ( fact1_nh3 )/(1 + fact2_nh3 *(1 + fact3_nh3*hplus) ! [nh4+] = enh3*hplus xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) ) - fact1_nh3 = (xk*xe*patm/hplus_scaling_factor)*(xnh3(i,k)+xnh4(i,k)) + fact1_nh3 = (xk*xe*patm/water_dissociation_constant)*(xnh3(i,k)+xnh4(i,k)) fact2_nh3 = xk*GAS_CONSTANT_L_ATM_MOL_K*temperature(i,k)*xl - fact3_nh3 = xe/hplus_scaling_factor + fact3_nh3 = xe/water_dissociation_constant !----------------------------------------------------------------- ! ... h2o effects !----------------------------------------------------------------- - Eh2o = hplus_scaling_factor + Eh2o = water_dissociation_constant !----------------------------------------------------------------- ! ... co2 effects @@ -454,7 +463,7 @@ subroutine calculate( state, & ! yposnet_lo and yposnet_hi = net positive ions for ! yph_lo and yph_hi !----------------------------------------------------------------- - do iter = 1, max_iterations + do iter = 1, MAX_ITERATIONS if (.not. present(specified_ph)) then if (iter == 1) then @@ -614,7 +623,7 @@ subroutine calculate( state, & !----------------------------------------------------------------- xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) ) - henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/hplus_scaling_factor) + henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/water_dissociation_constant) !----------------------------------------------------------------- ! ... o3 @@ -626,9 +635,8 @@ subroutine calculate( state, & ! ... for Ho2(g) -> H2o2(a) formation ! schwartz JGR, 1984, 11589 !------------------------------------------------------------------------ - kh4 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2) ho2s = kh0*xho2(i,k)*patm*(1._r8 + kh1/xph(i,k)) ! ho2s = ho2(a)+o2- - r1h2o2 = kh4*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s + r1h2o2 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2)*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s if ( cloud_borne ) then r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s @@ -717,7 +725,7 @@ subroutine calculate( state, & ! S(IV) + H2O2 = S(VI) !............................ - IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED + IF (XL .ge. MINIMUM_CLOUD_LIQUID_WATER) THEN !! WHEN CLOUD IS PRESENTED if (cloud_borne) then patm_x = patm From 9dd2e86a8067e79a8da4624174dd91c1ce4728bc Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 25 Feb 2025 09:18:32 -0800 Subject: [PATCH 09/12] address review comments --- .../aerosol/cloud_aqueous_chemistry.F90 | 76 +++++++++---------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index d1eb34e31b..cd62b81b69 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -31,12 +31,12 @@ module cloud_aqueous_chemistry ! Cloud chemistry species information type :: cloud_species_t - logical :: is_constant_ = .false. - integer :: state_index_ = CLOUD_INDEX_UNDEFINED ! index in the state vector or the fixed concentrations array character(len=:), allocatable :: name_ + integer :: state_index_ = CLOUD_INDEX_UNDEFINED ! index in the state vector or the fixed concentrations array + logical :: is_constant_ = .false. contains - procedure :: exists => cloud_species_exists - procedure :: mixing_ratio => cloud_species_get_mixing_ratio + procedure :: exists => cloud_species_exists + procedure :: mixing_ratio => cloud_species_get_mixing_ratio end type cloud_species_t interface cloud_species_t @@ -76,14 +76,10 @@ subroutine initialize() use carma_clouds, only : sox_cldaero_init #endif - logical :: modal_aerosols + logical :: is_modal_aerosols - call phys_getopts( prog_modal_aero_out=modal_aerosols ) - cloud_borne = modal_aerosols .or. carma_do_cloudborne - - !----------------------------------------------------------------- - ! ... get species indicies - !----------------------------------------------------------------- + call phys_getopts( prog_modal_aero_out=is_modal_aerosols ) + cloud_borne = is_modal_aerosols .or. carma_do_cloudborne so2 = cloud_species_t( 'SO2' ) nh3 = cloud_species_t( 'NH3' ) @@ -97,12 +93,16 @@ subroutine initialize() do_cloud_aqueous_chemistry = so2%exists() .and. h2o2%exists() .and. & o3%exists() .and. ho2%exists() - if (cloud_borne) then - do_cloud_aqueous_chemistry = do_cloud_aqueous_chemistry .and. & - h2so4%exists() - else - do_cloud_aqueous_chemistry = do_cloud_aqueous_chemistry .and. & - so4%exists() .and. nh3%exists() + if (do_cloud_aqueous_chemistry) then + if (cloud_borne) then + if (.not. h2so4%exists()) then + do_cloud_aqueous_chemistry = .false. + endif + else + if (.not. (so4%exists() .and. nh3%exists())) then + do_cloud_aqueous_chemistry = .false. + endif + endif endif if (masterproc) then @@ -111,20 +111,18 @@ subroutine initialize() do_cloud_aqueous_chemistry endif - if( do_cloud_aqueous_chemistry ) then - if (masterproc) then - write(iulog,*) '-----------------------------------------' - write(iulog,*) ' cloud aqueous chemistry is active' - write(iulog,*) '-----------------------------------------' - endif - else - if (masterproc) then - write(iulog,*) '-----------------------------------------' - write(iulog,*) ' cloud aqueous chemistry is inactive' - write(iulog,*) '-----------------------------------------' - endif - return + if (masterproc) then + if( do_cloud_aqueous_chemistry ) then + write(iulog,*) '-----------------------------------------' + write(iulog,*) ' cloud aqueous chemistry is active' + write(iulog,*) '-----------------------------------------' + else + write(iulog,*) '-----------------------------------------' + write(iulog,*) ' cloud aqueous chemistry is inactive' + write(iulog,*) '-----------------------------------------' + end if end if + if (.not. do_cloud_aqueous_chemistry) return call sox_cldaero_init() @@ -149,7 +147,7 @@ subroutine calculate( state, & cloud_water, & cloud_fraction, & cloud_droplet_number, & - air_number_density, & + air_number_density, & fixed_concentrations, & cloud_borne_aerosol_vmr, & species_vmr, & @@ -193,6 +191,8 @@ subroutine calculate( state, & !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- + type(physics_state), intent(in) :: state ! Physics state variables + type(physics_buffer_desc), pointer, intent(inout) :: pbuf(:) ! Physics buffer integer, intent(in) :: ncol ! num of columns in chunk integer, intent(in) :: lchnk ! chunk id integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array @@ -218,9 +218,6 @@ subroutine calculate( state, & real(r8), intent(out), optional :: aq_so4_production_from_h2o2_3d(:, :) ! 3D SO4 aqueous phase production due to H2O2 (kg/m2/s) real(r8), intent(out), optional :: aq_so4_production_from_o3_3d(:, :) ! 3D SO4 aqueous phase production due to O3 (kg/m2/s) - type(physics_state), intent(in) :: state ! Physics state variables - - type(physics_buffer_desc), pointer :: pbuf(:) !----------------------------------------------------------------------- ! ... Local variables @@ -302,9 +299,9 @@ subroutine calculate( state, & !----------------------------------------------------------------- do k = 1,pver air_mass_density_kg_l(:,k) = & - air_number_density(:,k) & ! /cm3(a) - * 1.e3_r8 & ! /L(a) - * BOLTZMANN/GAS_CONSTANT_DRY_AIR_J_KG_K ! Kg(a)/L(a) + air_number_density(:,k) & ! molecules(air) cm-3 + * 1.e3_r8 & ! L-1 + * BOLTZMANN/GAS_CONSTANT_DRY_AIR_J_KG_K ! kg(air) L-1 end do cldconc => sox_cldaero_create_obj( cloud_fraction, cloud_borne_aerosol_vmr, & @@ -837,7 +834,10 @@ end subroutine calculate !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- - ! Creates a cloud species object with state indices + ! Creates a cloud species object with the given name + ! + ! The name is used to determine the species index in the state arrays. + ! If the species is not found, the index is set to CLOUD_INDEX_UNDEFINED. function cloud_species_constructor( species_name ) result( this ) use mo_chem_utls, only : get_spc_ndx, get_inv_ndx From 62e287ff16afeca93bd70903ffc6001dd522a126 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 25 Feb 2025 09:23:16 -0800 Subject: [PATCH 10/12] remove duplicate log info --- src/chemistry/aerosol/cloud_aqueous_chemistry.F90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index cd62b81b69..e26e256822 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -105,12 +105,6 @@ subroutine initialize() endif endif - if (masterproc) then - write(iulog,*) 'cloud_aqueous_chemistry_init: '//& - 'do_cloud_aqueous_chemistry = ',& - do_cloud_aqueous_chemistry - endif - if (masterproc) then if( do_cloud_aqueous_chemistry ) then write(iulog,*) '-----------------------------------------' From 31efc1f8d5f6c732ccac55c726f1c6cbdd0a22b4 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Tue, 25 Feb 2025 11:58:28 -0800 Subject: [PATCH 11/12] fix units in description --- src/chemistry/aerosol/cloud_aqueous_chemistry.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index e26e256822..a68b8220d2 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -294,7 +294,7 @@ subroutine calculate( state, & do k = 1,pver air_mass_density_kg_l(:,k) = & air_number_density(:,k) & ! molecules(air) cm-3 - * 1.e3_r8 & ! L-1 + * 1.e3_r8 & ! molecules(air) L-1 * BOLTZMANN/GAS_CONSTANT_DRY_AIR_J_KG_K ! kg(air) L-1 end do From 32068a09af4a7801436b276347e79cabb6bb3a14 Mon Sep 17 00:00:00 2001 From: Matt Dawson Date: Thu, 27 Feb 2025 09:31:05 -0800 Subject: [PATCH 12/12] address review comments --- .../aerosol/cloud_aqueous_chemistry.F90 | 33 ++++++++++--------- .../cloud_chemistry_dependencies.F90 | 4 +-- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 index a68b8220d2..813ec82eb9 100644 --- a/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 +++ b/src/chemistry/aerosol/cloud_aqueous_chemistry.F90 @@ -54,6 +54,7 @@ module cloud_aqueous_chemistry real(r8), parameter :: GAS_CONSTANT_L_ATM_MOL_K = 8314.46261815324_r8*PASCAL_TO_ATM real(r8), parameter :: BOLTZMANN = 1.380649e-23_r8 ! J K-1 real(r8), parameter :: GAS_CONSTANT_DRY_AIR_J_KG_K = 287.0_r8 ! J kg-1 K-1 + real(r8), parameter :: SMALL_NUMBER = 1.e-30_r8 ! unitless contains @@ -251,8 +252,8 @@ subroutine calculate( state, & ! schwartz JGR, 1984, 11589 !----------------------------------------------------------------------- real(r8) :: ho2s ! ho2s = ho2(a)+o2- - real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s - real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s + real(r8) :: dh2o2_dt_mol_L_s ! prod(h2o2) by ho2 in mole/L(w)/s + real(r8) :: dh2o2_dt_vmr_s ! prod(h2o2) by ho2 in mix/s ! volume mixing ratios for cloud chemistry species real(r8), dimension(ncol,pver) :: xhno3, xh2o2, xso2, xso4, xno3, & @@ -627,20 +628,20 @@ subroutine calculate( state, & ! schwartz JGR, 1984, 11589 !------------------------------------------------------------------------ ho2s = kh0*xho2(i,k)*patm*(1._r8 + kh1/xph(i,k)) ! ho2s = ho2(a)+o2- - r1h2o2 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2)*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s + dh2o2_dt_mol_L_s = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2)*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s if ( cloud_borne ) then - r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s - / (1.e3_r8/AVOGADRO)*1.e+6_r8 & ! correct a bug here ???? + dh2o2_dt_vmr_s = dh2o2_dt_mol_L_s*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s + / (1.e3_r8/AVOGADRO)*1.e+6_r8 & ! correct a bug here ???? / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) else - r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s - * (1.e3_r8/AVOGADRO) & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s + dh2o2_dt_vmr_s = dh2o2_dt_mol_L_s*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s + * (1.e3_r8/AVOGADRO) & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s / (midpoint_pressure(i,k)/(BOLTZMANN*temperature(i,k))) ! /cm3(a)/s / air-den = mix-ratio/s endif if ( .not. cloud_borne) then ! this seems to be specific to aerosols that are not cloud borne - xh2o2(i,k) = xh2o2(i,k) + r2h2o2*time_step ! updated h2o2 by het production + xh2o2(i,k) = xh2o2(i,k) + dh2o2_dt_vmr_s*time_step ! updated h2o2 based on heterogeneous production endif !----------------------------------------------- @@ -732,14 +733,14 @@ subroutine calculate( state, & * heso2(i,k) * so2g * patm_x ! [M/s] endif - dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] - * xl & ! [mole/L(a)/s] + dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] + * xl & ! [mole/L(a)/s] / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] / air_number_density(i,k) delta_concentration = dso4_dt*time_step - delta_concentration = max(delta_concentration, 1.e-30_r8) + delta_concentration = max(delta_concentration, SMALL_NUMBER) xso4_init(i,k)=xso4(i,k) @@ -748,9 +749,9 @@ subroutine calculate( state, & xso4(i,k)=xso4(i,k)+xso2(i,k) if (cloud_borne) then xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) - xso2(i,k)=1.e-20_r8 + xso2(i,k)=1.e-20_r8 ! TODO: See if SMALL_NUMBER is more appropriate else ! ???? bug ???? - xso2(i,k)=1.e-20_r8 + xso2(i,k)=1.e-20_r8 ! TODO: See if SMALL_NUMBER is more appropriate xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) endif else @@ -780,19 +781,19 @@ subroutine calculate( state, & dso4_dt = k_siv_o3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s] - dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] + dso4_dt = dso4_dt & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / (1.e3_r8/AVOGADRO) & ! [/L(a)/s] / air_number_density(i,k) ! [mixing ratio/s] delta_concentration = dso4_dt*time_step - delta_concentration = max(delta_concentration, 1.e-30_r8) + delta_concentration = max(delta_concentration, SMALL_NUMBER) xso4_init(i,k)=xso4(i,k) if (delta_concentration .gt. xso2(i,k)) then xso4(i,k) = xso4(i,k) + xso2(i,k) - xso2(i,k) = 1.e-20_r8 + xso2(i,k) = 1.e-20_r8 ! TODO: See if SMALL_NUMBER is more appropriate else xso4(i,k) = xso4(i,k) + delta_concentration xso2(i,k) = xso2(i,k) - delta_concentration diff --git a/test/chemistry/cloud_aqueous_chemistry/cloud_chemistry_dependencies.F90 b/test/chemistry/cloud_aqueous_chemistry/cloud_chemistry_dependencies.F90 index 63b3d3f5a1..dbabae341a 100644 --- a/test/chemistry/cloud_aqueous_chemistry/cloud_chemistry_dependencies.F90 +++ b/test/chemistry/cloud_aqueous_chemistry/cloud_chemistry_dependencies.F90 @@ -23,14 +23,14 @@ module physics_buffer integer, parameter :: dtype_r8 = 1 integer, parameter :: pbuf_get_index = 1 integer, parameter :: pbuf_add_field = 1 - type, public :: physics_buffer_desc + type :: physics_buffer_desc end type physics_buffer_desc end module physics_buffer module physics_types implicit none public :: physics_state - type, public :: physics_state + type :: physics_state end type physics_state end module physics_types