Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions .gitmodules
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@

[submodule "carma"]
path = src/physics/carma/base
url = https://github.com/fvitt/CARMA_base.git
url = https://github.com/ESCOMP/CARMA_base.git
fxrequired = AlwaysRequired
fxtag = 4dbf023
fxtag = carma4_09
fxDONOTUSEurl = https://github.com/ESCOMP/CARMA_base.git

[submodule "pumas"]
Expand Down
902 changes: 902 additions & 0 deletions src/chemistry/aerosol/cloud_aqueous_chemistry.F90

Large diffs are not rendered by default.

152 changes: 152 additions & 0 deletions src/chemistry/aerosol/cloud_aqueous_chemistry_snapshot.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
! Captures inputs and output to cloud aqueous chemistry routines
module cloud_aqueous_chemistry_snapshot

use shr_kind_mod, only : r8 => shr_kind_r8

implicit none
private

public :: cloud_snapshot_init, cloud_snapshot_capture_input, cloud_snapshot_capture_output

contains

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine cloud_snapshot_init()

use cam_history, only : addfld, horiz_only

integer :: i_elem
character(len=10) :: index_string

call addfld( 'cloud_press_in', (/ 'lev' /), 'I', 'Pa', 'mid-point pressure' )
call addfld( 'cloud_pdel_in', (/ 'lev' /), 'I', 'Pa', 'pressure thickness of levels' )
call addfld( 'cloud_tfld_in', (/ 'lev' /), 'I', 'K', 'temperature' )
call addfld( 'cloud_mbar_in', (/ 'lev' /), 'I', 'AMU', 'mean wet atmopheric mass' )
call addfld( 'cloud_lwc_in', (/ 'lev' /), 'I', 'kg kg-1', 'cloud liquid water content' )
call addfld( 'cloud_cldfrc_in', (/ 'lev' /), 'I', 'unitless', 'cloud fraction' )
call addfld( 'cloud_cldnum_in', (/ 'lev' /), 'I', 'kg-1', 'droplet number concentration' )
call addfld( 'cloud_xhnm_in', (/ 'lev' /), 'I', 'cm-3', 'total atmospheric density' )
do i_elem = 1, 7
write(index_string, '(I10)') i_elem
call addfld( 'cloud_invariants_'//trim(adjustl(index_string))//'_in', (/ 'lev' /), 'I', 'unknown', 'invariant' )
end do
do i_elem = 1, 26
write(index_string, '(I10)') i_elem
call addfld( 'cloud_qcw_'//trim(adjustl(index_string))//'_in', (/ 'lev' /), 'I', 'vmr', 'cloud-borne aerosol' )
call addfld( 'cloud_qcw_'//trim(adjustl(index_string))//'_out', (/ 'lev' /), 'I', 'vmr', 'transported species' )
call addfld( 'cloud_qin_'//trim(adjustl(index_string))//'_in', (/ 'lev' /), 'I', 'vmr', 'transported species' )
call addfld( 'cloud_qin_'//trim(adjustl(index_string))//'_out', (/ 'lev' /), 'I', 'vmr', 'transported species' )
end do
call addfld( 'cloud_xphlwc_out', (/ 'lev' /), 'I', 'unitless', 'ph value multiplied by cloud liquid water content' )
do i_elem = 1, 4
write(index_string, '(I10)') i_elem
call addfld( 'cloud_aqso4_'//trim(adjustl(index_string))//'_out', horiz_only, 'I', 'unknown', 'aqueous phase chemistry' )
call addfld( 'cloud_aqh2so4_'//trim(adjustl(index_string))//'_out', horiz_only, 'I', 'unknown', 'aqueous phase chemistry' )
end do
call addfld( 'cloud_aqso4_h2o2_out', horiz_only, 'I', 'kg m-2', 'SO4 aqueous phase chemistry due to H2O2?' )
call addfld( 'cloud_aqso4_o3_out', horiz_only, 'I', 'kg m-2', 'SO4 aqueous phase chemistry due to O3?' )

end subroutine cloud_snapshot_init

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine cloud_snapshot_capture_input(ncol, lchnk, loffset, dtime, press, &
pdel, tfld, mbar, lwc, cldfrc, cldnum, xhnm, invariants, qcw, qin)

use cam_history, only : outfld
use cam_logfile, only : iulog
use spmd_utils, only : is_main_process => masterproc

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(in) :: qcw(:,:,:) ! cloud-borne aerosol (vmr)
real(r8), intent(in) :: qin(:,:,:) ! transported species ( vmr )

integer :: i_elem
character(len=10) :: index_string

#if 0
if (is_main_process) then
write(iulog,*) "*****************************"
write(iulog,*) "Cloud Chemistry scalar inputs"
write(iulog,*) "*****************************"
write(iulog,*) "ncol: ", ncol
write(iulog,*) "lchnk: ", lchnk
write(iulog,*) "loffset: ", loffset
write(iulog,*) "dtime: ", dtime
write(iulog,*) "press dims: ", size(press, dim=1), size(press, dim=2)
write(iulog,*) "invariants dims: ", size(invariants, dim=1), size(invariants, dim=2), size(invariants, dim=3)
write(iulog,*) "qcw dims: ", size(qcw, dim=1), size(qcw, dim=2), size(qcw, dim=3)
write(iulog,*) "qin dims: ", size(qin, dim=1), size(qin, dim=2), size(qin, dim=3)
end if
#endif
call outfld( 'cloud_press_in', press(:ncol,:), ncol, lchnk )
call outfld( 'cloud_pdel_in', pdel(:ncol,:), ncol, lchnk )
call outfld( 'cloud_tfld_in', tfld(:ncol,:), ncol, lchnk )
call outfld( 'cloud_mbar_in', mbar(:ncol,:), ncol, lchnk )
call outfld( 'cloud_lwc_in', lwc(:ncol,:), ncol, lchnk )
call outfld( 'cloud_cldfrc_in', cldfrc(:ncol,:), ncol, lchnk )
call outfld( 'cloud_cldnum_in', cldnum(:ncol,:), ncol, lchnk )
call outfld( 'cloud_xhnm_in', xhnm(:ncol,:), ncol, lchnk )
do i_elem = 1, 7
write(index_string, '(I10)') i_elem
call outfld( 'cloud_invariants_'//trim(adjustl(index_string))//'_in', invariants(:ncol,:,i_elem), ncol, lchnk )
end do
do i_elem = 1, 26
write(index_string, '(I10)') i_elem
call outfld( 'cloud_qcw_'//trim(adjustl(index_string))//'_in', qcw(:ncol,:,i_elem), ncol, lchnk )
call outfld( 'cloud_qin_'//trim(adjustl(index_string))//'_in', qin(:ncol,:,i_elem), ncol, lchnk )
end do

end subroutine cloud_snapshot_capture_input

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine cloud_snapshot_capture_output(ncol, lchnk, qcw, qin, xphlwc, &
aqso4, aqh2so4, aqso4_h2o2, aqso4_o3)

use cam_history, only : outfld

integer, intent(in) :: ncol ! num of columns in chunk
integer, intent(in) :: lchnk ! chunk id
real(r8), target, intent(in) :: qcw(:,:,:) ! cloud-borne aerosol (vmr)
real(r8), intent(in) :: qin(:,:,:) ! transported species ( vmr )
real(r8), intent(in) :: xphlwc(:,:) ! pH value multiplied by lwc

real(r8), intent(in) :: aqso4(:,:) ! aqueous phase chemistry
real(r8), intent(in) :: aqh2so4(:,:) ! aqueous phase chemistry
real(r8), intent(in) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2)
real(r8), intent(in) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2)

integer :: i_elem
character(len=10) :: index_string

call outfld( 'cloud_xphlwc_out', xphlwc(:ncol,:), ncol, lchnk )
call outfld( 'cloud_aqso4_h2o2_out', aqso4_h2o2(:ncol), ncol, lchnk )
call outfld( 'cloud_aqso4_o3_out', aqso4_o3(:ncol), ncol, lchnk )
do i_elem = 1, 4
write(index_string, '(I10)') i_elem
call outfld( 'cloud_aqso4_'//trim(adjustl(index_string))//'_out', aqso4(:ncol,i_elem), ncol, lchnk )
call outfld( 'cloud_aqh2so4_'//trim(adjustl(index_string))//'_out', aqh2so4(:ncol,i_elem), ncol, lchnk )
end do
do i_elem = 1, 26
write(index_string, '(I10)') i_elem
call outfld( 'cloud_qcw_'//trim(adjustl(index_string))//'_out', qcw(:ncol,:,i_elem), ncol, lchnk )
call outfld( 'cloud_qin_'//trim(adjustl(index_string))//'_out', qin(:ncol,:,i_elem), ncol, lchnk )
end do

end subroutine cloud_snapshot_capture_output

end module cloud_aqueous_chemistry_snapshot
149 changes: 149 additions & 0 deletions src/chemistry/aerosol/cloud_utilities.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
!----------------------------------------------------------------------------------
! low level utility module for cloud aerosols
!
! Created by Francis Vitt
!----------------------------------------------------------------------------------
module cloud_utilities

use shr_kind_mod, only : r8 => shr_kind_r8
use ppgrid, only : pcols, pver

implicit none
private

public :: cldaero_uptakerate
public :: cldaero_conc_t
public :: cldaero_allocate
public :: cldaero_deallocate

type cldaero_conc_t
real(r8), pointer :: so4c(:,:)
real(r8), pointer :: nh4c(:,:)
real(r8), pointer :: no3c(:,:)
real(r8), pointer :: xlwc(:,:)
real(r8) :: so4_fact
end type cldaero_conc_t

contains

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
function cldaero_allocate( ) result( cldconc )
type(cldaero_conc_t), pointer:: cldconc

allocate( cldconc )
allocate( cldconc%so4c(pcols,pver) )
allocate( cldconc%nh4c(pcols,pver) )
allocate( cldconc%no3c(pcols,pver) )
allocate( cldconc%xlwc(pcols,pver) )

cldconc%so4c(:,:) = 0._r8
cldconc%nh4c(:,:) = 0._r8
cldconc%no3c(:,:) = 0._r8
cldconc%xlwc(:,:) = 0._r8
cldconc%so4_fact = 2._r8

end function cldaero_allocate

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
subroutine cldaero_deallocate( cldconc )
type(cldaero_conc_t), pointer :: cldconc

if ( associated(cldconc%so4c) ) then
deallocate(cldconc%so4c)
nullify(cldconc%so4c)
endif

if ( associated(cldconc%nh4c) ) then
deallocate(cldconc%nh4c)
nullify(cldconc%nh4c)
endif

if ( associated(cldconc%no3c) ) then
deallocate(cldconc%no3c)
nullify(cldconc%no3c)
endif

if ( associated(cldconc%xlwc) ) then
deallocate(cldconc%xlwc)
nullify(cldconc%xlwc)
endif

deallocate( cldconc )
nullify( cldconc )

end subroutine cldaero_deallocate

!----------------------------------------------------------------------------------
! utility function for cloud-borne aerosols
!----------------------------------------------------------------------------------

function cldaero_uptakerate( xl, cldnum, cfact, cldfrc, tfld, press ) result( uptkrate )
use mo_constants, only : pi
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We use a constant from another module here, but in cloud_aqueous_chemistry.F90 we defined several constants within the code. We should probably push those constants into mo_constants if they don't exist


real(r8), intent(in) :: xl, cldnum, cfact, cldfrc, tfld, press

real(r8) :: uptkrate

real(r8) :: &
rad_cd, radxnum_cd, num_cd, &
gasdiffus, gasspeed, knudsen, &
fuchs_sutugin, volx34pi_cd

!-----------------------------------------------------------------------
! compute uptake of h2so4 and msa to cloud water
!
! first-order uptake rate is
! 4*pi*(drop radius)*(drop number conc)
! *(gas diffusivity)*(fuchs sutugin correction)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
! *(gas diffusivity)*(fuchs sutugin correction)
! *(gas diffusivity)*(fuchs sutugin correction)
! Fuchs Sutugin correction:
! Fuchs, N.A., Sutugin, A.G., 1971. HIGH-DISPERSED AEROSOLS, in: Hidy, G.M., Brock, J.R. (Eds.),
! Topics in Current Aerosol Research, International Reviews in Aerosol Physics and Chemistry. Pergamon, p. 1.
! https://doi.org/10.1016/B978-0-08-016674-2.50006-6

Is this a correct reference? Took me forever to track down

Copy link
Collaborator

@dmleung dmleung Dec 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes this is correct. A lot of papers on sulfur chemistry also used Fuchs' formulation. E.g., see
https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2005jd005870


! num_cd = (drop number conc in 1/cm^3)
num_cd = 1.0e-3_r8*cldnum*cfact/cldfrc
num_cd = max( num_cd, 0.0_r8 )

! rad_cd = (drop radius in cm), computed from liquid water and drop number,
! then bounded by 0.5 and 50.0 micrometers
Comment on lines +105 to +106
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but WHY is it bounded? Because a reference says we have to? Because numerically only values beteween 0.5 - 50 are valid?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to line 115 they applied the bounds at line 118 etc to avoid occasional unphysical values. I have no idea why that is the case though.

! radxnum_cd = (drop radius)*(drop number conc)
! volx34pi_cd = (3/4*pi) * (liquid water volume in cm^3/cm^3)

volx34pi_cd = xl*0.75_r8/pi
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the significance of 3/(4 pi)? WHY

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given xl (liquid water content) they want to estimate the droplet radius.
V = 4*pi/3 * R^3
So
V * 0.75/pi is related to the radius. See line 113 below.


! following holds because volx34pi_cd = num_cd*(rad_cd**3)
radxnum_cd = (volx34pi_cd*num_cd*num_cd)**0.3333333_r8

! apply bounds to rad_cd to avoid the occasional unphysical value
if (radxnum_cd .le. volx34pi_cd*4.0e4_r8) then
radxnum_cd = volx34pi_cd*4.0e4_r8
rad_cd = 50.0e-4_r8
else if (radxnum_cd .ge. volx34pi_cd*4.0e8_r8) then
radxnum_cd = volx34pi_cd*4.0e8_r8
rad_cd = 0.5e-4_r8
else
rad_cd = radxnum_cd/num_cd
end if
Comment on lines +116 to +124
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is 4.0e4_r8 arbitrary? Why this value? what makes radxnum_cd and why is 4.0e4_r8 related to that?

Copy link
Collaborator

@dmleung dmleung Dec 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think theses are just bounds. They didn't want radxnum_cd to be smaller than volx34pi_cd*4.0e4 and bigger than volx34pi_cd*4.0e8. Again, I don't know why. This is a question for Mary.


! gasdiffus = h2so4 gas diffusivity from mosaic code (cm^2/s)
! (pmid must be Pa)
gasdiffus = 0.557_r8 * (tfld**1.75_r8) / press

! gasspeed = h2so4 gas mean molecular speed from mosaic code (cm/s)
gasspeed = 1.455e4_r8 * sqrt(tfld/98.0_r8)

! knudsen number
knudsen = 3.0_r8*gasdiffus/(gasspeed*rad_cd)
Comment on lines +126 to +134
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is from mosaic, do we have a reference for where these values come from and why? Do we know if they are always valid or does mosaic assume something we don't?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can think of a reference for knudsen number (I am sure some textbooks will have this equation) but Mary probably knows better


! following assumes accomodation coefficient = 0.65
! (Adams & Seinfeld, 2002, JGR, and references therein)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it this paper? If so, let's add a full bibliographical citation that includes the DOI

! fuchs_sutugin = (0.75*accom*(1. + knudsen)) /
! (knudsen*(1.0 + knudsen + 0.283*accom) + 0.75*accom)
fuchs_sutugin = (0.4875_r8*(1._r8 + knudsen)) / &
(knudsen*(1.184_r8 + knudsen) + 0.4875_r8)

! instantaneous uptake rate
uptkrate = 12.56637_r8*radxnum_cd*gasdiffus*fuchs_sutugin

end function cldaero_uptakerate

end module cloud_utilities

Loading
Loading