Cloud chemistry module review

[mattldawson]

Reimplements four cloud chemistry modules with new file and module names so that the code can be reviewed.

closes #1

From issue #1, things to consider in review:

• missing descriptions for units for variables
• missing descriptions of functions and function arguments
• unused variables
• cryptically named variables
• dead code blocks
• potential to break up large code blocks, and introduce reusable/testable functions
• compare MAM and CARMA code for consistency
• identify host-model vs aerosol-model portions of the code
• separate configuration data from algorithms
• magic numbers

[dwfncar]

This loop can be condensed and re-ordered.

[tilmes]

not sure how this would work?

[dwfncar]

Could split subroutine into sox_cldaero_tend to compute tendencies and sox_cldaero_update to apply tendencies.

[dwfncar]

Index gymnastics is confusing here.

[tilmes]

This indexing can be removed, similar to what is done in carma_clouds.F90 Line 72 to 102.

[K20shores]

Reviewing cloud_aqueous_chemistry.F90

[K20shores]

What does cloud_borne logically mean?

[mattldawson]

I agree, this needs to be described clearly in the source code. Recent conversations with people who have worked on CAM for years revealed that there is a lot of confusion about what this actually means.

[K20shores]

1 / avogadro?

[K20shores]

Do we maybe have a description of these, where they are from (reference?), units?

[dmleung]

I believe none of these variables are used (grep them and there is only one occurrence -- here). I have no idea why they are here. I think we can remove them. Other people can confirm this.

[K20shores]

Do we have a reference for these or is this just known in the ether of science? Is it applicable forever and always or only for a certain range of environmental conditions?

[dmleung]

I am not sure about this. This module is Mary's expertise and she can confirm and add references.

[K20shores]

Suggested change

	real(r8), parameter :: xkw = 1.e-14_r8 ! water acidity
	real(r8), parameter :: water_acidity = 1.e-14_r8 ! water acidity

So let's call it that

[K20shores]

well did the bug get corrected? Is the correction const0 or the multiplication of const0

[dmleung]

I have no idea about this. Someone needs to figure this out... Maybe in our regular Tue meetings.

[K20shores]

are we sure it's specific? If so, why is it?

[dmleung]

I remember Simone and Mary talked about this a while ago and nobody knew what happened here. I am not sure if they now have an updated view on this...

[K20shores]

do we have a citation for these values?

[K20shores]

also these values

[K20shores]

presented what? why is 1e-8 the value we care about? Is it empirical? Numerical? If numerical, was it derived when we only had single precision and we could make it smaller now?

[dmleung]

Yea it's a small positive definite. You can use greater than zero instead of great or equal to a positive definite, I suppose. There are a lot of these small numbers in different parts of CAM-Chem. There are 1e-20, 1e-30, and so on.

[mattldawson]

@tilmes @marybarthbrock - Could we get a high-level description of the purpose of this module in science terms?

[dmleung]

I would say the purpose of this module is to calculate the sulfate formation due to various oxidation/condensation pathways (e.g., SO2 oxidized by ozone, SO2 oxidized by H2O2, H2SO4 condensation) of cloud chemistry. In the cloud chemistry subroutine setsox(), gas-phase species will be consumed and cloud-borne aerosol species be produced, and the main outputs are these two arrays: qin() and qcw(). Sulfate production rates (e.g., aqso4) and cloud pH changes (xphlwc) are also calculated and outputted by setsox().

[marybarthbrock]

Looks good to me

[mattldawson]

We should state what array(s) these are indices for

[dmleung]

We can write something like
"these indices are mainly for two arrays: qin and qcw -- qin is the volume mixing ratio (mol/mol or m3/m3) for chemical tracers (gas and aerosol), and qcw is the volume mixing ratio of cloud-borne aerosols. e.g., qin(i,k,id_so2) where i is grid, k is vertical level, and id_so2 is the index for SO2."
Feel free to rewrite this in a succinct way.

[mattldawson]

this should have a description of what it means and how it is used

[dmleung]

Mary or Francis will have a better idea on this

[mattldawson]

Suggested change

	subroutine sox_inti
	subroutine initialize_cloud_aqueous_chemistry

[mattldawson]

we should include a brief description of what is being initialized

[mattldawson]

there is repetitive but inconsistent logic used here (e.g., different names for sulfate, some species checked for get_inv_ndx() and others not). I would suggest a data structure for species-specific data like this, with some common functionality for identifying indices, etc., and standardization of logic across all species for consistency.

[mattldawson]

@tilmes @marybarthbrock - could we get a brief (2-3 sentence) high-level description of what this function is doing in science terms? (there is a comment inside the function that might be a starting point, but I wasn't sure if this describes the entirety of what this function does)

We should also rename the function to be more descriptive (e.g., solve_cloud_aqueous_chemistry() or similar)

[mattldawson]

Maybe we can consider using Doxygen style comments

[dmleung]

For me, a very brief description of setsox will be
"This subroutine calculates the formation of sulfate and updates sulfate concentration due to cloud aqueous chemistry (i.e., updating qin / qcw). This subroutine also outputs the production rates (in kg/m2/s) of various reactions of sulfur species with various oxidants (i.e., outputting things like aqso4_h2o2 and aqh2so4)."
Something along the line... Mary and Simone et al. can add onto this...

[mattldawson]

How are the CPP conditional blocks being used in conjunction with the abstract aerosol interface? It seems like we might want to make sure we're using a consistent approach to permitting the use of multiple aerosol modules.

[mattldawson]

after some discussion, we should look into how USE_CLOUDS and USE_MAM are being set and used

[dmleung]

Units of these production rates should be kg/m2/s instead of kg/m2.

[K20shores]

Suggested change

	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(:,:) ! aqueous phase chemistry (kg/m2/s)
	real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry (kg/m2/s)
	real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2/s)

[K20shores]

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

[K20shores]

cloud_utilities.F90 reivew

[K20shores]

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

[dmleung]

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.

[K20shores]

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

[K20shores]

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

[dmleung]

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.

[K20shores]

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

[dmleung]

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.

[K20shores]

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?

[dmleung]

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

[K20shores]

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

[dmleung]

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

[dmleung]

The unit of this is also kg/m2/s

[dmleung]

These two have a unit of kg/m2/s too.

[K20shores]

carma_clouds.F90

[K20shores]

Can we just delete this since it's commented out?

[dmleung]

Yes I think so

[K20shores]

In other places (the backward euler solver) small is set to 1e-40. I imagine there are other versions of small throughout CAM. Do we happen to know why this particular value of small is more useful for clouds?

[dmleung]

I think it depends on the order of magnitudes of the variables concerned. in the code here the order of magnitudes of chemical tracers (e.g., SO2, H2O2) are typically anywhere between 1e-14 (0.01 ppt) to 1e-5 or (10 ppm) in volume mixing ratios. People then want to disregard the reactions once the concentrations are low enough.
The small_value can definitely change depending on the situation and variable concerned. E.g., the small_value for relative humidity should be 0.1 instead of 1e-20. For solving differential equations, I think it makes sense that people try to make things as accurate as possible by using smaller small_values even for chemical tracers. But I think there's always wiggle room, e.g., setting 1e-30 in the solver probably doesn't hurt at all.

[K20shores]

Perhaps we can delete commented out code?

Suggested change

!integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr

[dmleung]

I guess Mary is the best to answer this. Mary needs to review the code changes...

[K20shores]

Suggested change

	!st if (id_nh3>0) then
	!st delnh3 = nh3g(i,k) - xnh3(i,k)
	!st delnh4 = - delnh3
	!st endif

[K20shores]

Suggested change

	!st dmsadt_gasuptk_tomsa = dmsadt_gasuptk
	!st if (ntot_msa_c == 0) then
	!st dmsadt_gasuptk_tomsa = 0.0_r8
	!st dmsadt_gasuptk_toso4 = dmsadt_gasuptk
	!st end if

[K20shores]

Suggested change

!qin(i,k,id_h2so4) = MAX( qin(i,k,id_h2so4), small_value )

[K20shores]

Suggested change

! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) )

[K20shores]

Suggested change

! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )

[K20shores]

Suggested change

	!st if (id_nh3>0) then
	!st dqdt_aq = delnh3/dtime*cldfrc(i,k)
	!st dqdt = dqdt_aq
	!st qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime
	!st endif

[boulderdaze]

What does inv_* stand for?

[dmleung]

According to Francis and Simone, inv_* stands for invariant. So, it sounds like in certain compsets people holds these tracers (SO2, NH3, H2O2, etc.) constant. In most cases nowadays I don't think we use invariant tracers, so the code mostly uses id_so2 instead of inv_so2, for instance. But, it's used for some compsets.

Francis is the one to answer this question.

[boulderdaze]

What does masterproc mean? (There is a typo: sox_inti -> sox_init.)
It doesn't appear to provide useful information, and it also duplicates the following logs.

[dmleung]

masterproc comes from utils/spmd_utils.F90 and is a flag called when scientists want to write() something in a log file.
Yes, that's a typo and can be changed to sox_init.

[boulderdaze]

Some variables (line435, line 492 and so on) that are constant should be defined outside the loop.

[dmleung]

I am not sure about this. Line 435 and line 458 are different, so the existing approach is to declare xe as real(r8) and update xe for different reactions. xe is a constant in line 435 but get updated in line 458 as a function of work1 (which is dependent on temperature tfld).

[boulderdaze]

Should we verify if this code is valid?

[dmleung]

I think these lines mean that if you remove lines 821-824, lines 818-820 sometimes will not work. This is probably because xso2 was too small (smaller than 1.e-20_r8), but I'm not sure why that is the case. This is something that Mary should comment on.

BTW, I think line 822 can be
xso2(i,k)=small_value

[boulderdaze]

This function is large and difficult to test. I think it would be beneficial to break it down into smaller functions for each step, or a series of steps to make it more testable.

[dmleung]

Sure, I think since the code is being refactored it is up to the software engineers to design how the code looks like, as long as it looks good to Mary (this is her science after all).

[marybarthbrock]

I am a-okay with breaking the code down into smaller functions.

[boulderdaze]

My understanding of this function is that it sets variables to update cldaero by calling sox_cldaero_update(). There are also variables with intent(out) (eg. xphlwc in line 880) that aren't used for updating cldaero. What is the purpose of those variables, and how are they used?

[dmleung]

I think there are two questions here:

1. sox_cldaero_update() is a function that calculates and outputs the production rates of sulfate (so4) with cloud chemistry through several pathways (e.g., so4 formed due to o3 oxidation, due to h2o2 oxidation, etc). So, looking at sox_cldaero_update() under carmo_aero/carma_cloud.F90, the only intent(out) variables are aqso4, aqh2so4, aqso4_h2o2, ... which are sulfate production rates (kg/m2/s). sox_cldaero_update() is called here in aerosol/cloud_aqueous_chemistry.F90, so setsox() gets these arrays as intent(out) too (see line 215). And, when setsox() is called in the aerosol interface, we will see that they are outputted by setsox() just to be set as history variables, as follows:

     call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk)
     call outfld( 'AQSO4_O3',   aqso4_o3(:ncol),   ncol, lchnk)
     call outfld( 'XPH_LWC',    xphlwc(:ncol,:),   ncol, lchnk )
   

2. I think lines 875-878 and lines 880-887 are two different things. xphlwc is a separate intent(out) variable that setsox() wants to output when setsox() is called. It is also a history variable in the aerosol interface, as follow:

     call outfld( 'XPH_LWC',    xphlwc(:ncol,:),   ncol, lchnk )
   

=====

I copied the above code for outfld() from the CESM2.2 version. In the older CESM2.2, modal_aero/aero_model.F90 calls setsox(), setsox() calculates xphlwc and calls sox_cldaero_update() that calculates aqso4 etc, and aero_model.F90 gets them and do outfld(), so scientists get the history fields for plotting and analyses. I'm not sure where setsox() is called in this new CAM tag, and Francis and Simone will know better where the aerosol interface is. You can also grep it in this new CAM code.

But basically in one sentence, the purpose of sox_cldaero_update() in CESM2.2 is to output sulfate formation rates to aero_model.F90 so scientists get the history fields.

[boulderdaze]

If fwetrem is zero, then dqdt_wr should also be zero, which seems to be unnecessary logic. Could you confirm?

[dmleung]

Agreed. I'm not sure why h2o2 and so2 wet removals were decided to set as zero. I think Mary should make a decision on whether fwetrem should be non-zero (i.e., using line 424 instead of line 425).

But yes, if we just keep the code as is, I would say this can be deleted or commented out (preferred).

[marybarthbrock]

I believe fwetrem for h2o2 and so2 is set to zero because wet removal (or wet scavenging) is computed elsewhere in CAM-chem (the Neu & Prather scheme, mo_neu_wetdep.F90). Let's confirm this with Francis and Simone. If it's correct, then I suggest removing these lines of code.