Unify vertical diffusion calculations in CAM #534

[JulioTBacmeister]

Vertical diffusion is calculated in several places within CAM physics: CLUBB, several times within the GW codes, and possibly other places when molecular diffusion is on. We should accumulate the various diffusivities and do a single diffusion calculation for each species. This would make UBCs for example easier to track as well as being more economical.

[adamrher]

@tilmes @fvitt @dan800 we are trying to list all the places where there is subgrid vertical transport of species. Could you double check my understanding below? This is initially to understand what is currently going on for the default CAM6 physics (not CARMA).
(pls consult NCAR/amwg_dev#76 on different tracer categories: water tracers, chemical tracers, aerosol tracers)

Pure Vertical Diffusion
CLUBB - diffuses all tracers except aerosols and NUMLIQ. Does not diffuse T/thl or Q/qt (prognostic!)
MAM (subroutine explmix in ndrop.F90) - diffuses aerosols
GW - @JulioTBacmeister?
vertical_diffusion.F90 - molecular diffusivity (?)

Transport from ZM convection
ZM convtran phase 1 - water tracers
ZM convtran phase 2 - chemical tracers (+ aerosols if convproc_do_aer=F?)
MAM - aerosols if convproc_do_aer=T

Surface Emissions
vertical_diffusion.F90 - all tracer emissions (?) dumped into lowest model level (except Q->LHFLX, BC in CLUBB)
[edit - Hui Wan's showed that this should occur upstream of subgrid vertical transport, and this mixing should be upstream of dry deposition. We currently don't do this and the result is emissions dumped into lowest level get deposited by dry dep before being given the opportunity to be mixed vertically.]

[tilmes]

Here is what I understand from the code: Convective transport for both MAM and CARMA (even though not of interest currently) are only performed in src/chemistry/modal_aero/modal_aero_convproc.F90 if convproc_do_aer is set to true (in the namelist).
The setttings are done in chemistry/modal_aero/modal_aero_data.F90: call cnst_set_convtran2(numptr_amode(m), cam_do_aero_conv). Else, aerosols are convected in the zm routine.

I agree that ndrop is doing diffusion of aerosols.

[tilmes]

@adamrher regarding the diffusion in ndrop, we need to check whether all aerosols or only cloud-borne aerosols experience diffusion, I talked to Xiaohong and he mentioned that the later would be true.. need to check!

[adamrher]

This could not be any more confusing. CLUBB is taking special care to set the MAM aerosols to false, so as to not double-count mixing already done in ndrop.F90. The exact aerosols it sets to false could be changing on the fly, as they are activated in cloud droplets ... I'd need to run w/ some print statements to be sure. See how CLUBB omits the MAM aerosols here:

CAM/src/physics/cam/clubb_intr.F90

Lines 952 to 974 in 0404b72

	if (prog_modal_aero) then
	! Turn off modal aerosols and decrement edsclr_dim accordingly
	call rad_cnst_get_info(0, nmodes=nmodes)
	do m = 1, nmodes
	call rad_cnst_get_mode_num_idx(m, lptr)
	lq(lptr)=.false.
	edsclr_dim = edsclr_dim-1
	call rad_cnst_get_info(0, m, nspec=nspec)
	do l = 1, nspec
	call rad_cnst_get_mam_mmr_idx(m, l, lptr)
	lq(lptr)=.false.
	edsclr_dim = edsclr_dim-1
	end do
	end do
	! In addition, if running with MAM, droplet number is transported
	! in dropmixnuc, therefore we do NOT want CLUBB to apply transport
	! tendencies to avoid double counted. Else, we apply tendencies.
	lq(ixnumliq) = .false.
	edsclr_dim = edsclr_dim-1
	endif

@XiZhao1013 can you help us understand this better, or provide Xiaohong's git handle?

[tilmes]

In ndrop, explmix is called and applies vertical diffusion to both cloud borne and interstitial aerosols, just checked.

[XiZhao1013]

@adamrher Thanks, I have forwarded this to Xiaohong and our group, they may add more comments later.

[XiZhao1013]

Yes, as @adamrher said. The vertical mixing of aerosol is handled in MAM (ndrop), not in CLUBB.
Thanks, Zheng Lu (at Xiaohong's Group), and @mingxuanwupnnl for talking about this with me.

[adamrher]

In the spirit of unifying diffusion operations in CAM, can anyone see a problem with computing the diffusion coefficients in ndrop.F90, but applying them in CLUBB's diffusion solver? The only issue I could see is that CLUBB is located upstream of where ndrop is called in the time-loop (microp_aero_run I think).

@huiwanpnnl bringing you in to this convo in case you're interested.

[tilmes]

Aerosols are changed from interstitial to cloud borne during the diffusion process, so I guess, you could bring this part to CLUBB but have to add those processes there.

[XiZhao1013]

Xiaohong instructed me to post his comments here:
"I know some years ago, only mixing for cloud droplet number is treated in ndrop (dropmixnuc). It is done to have implicit solver of cloud droplet number change due to mixing with environment (related to resuspension of cloud borne aerosols). At that time Steve Ghan and Dick Easter were specifically interested in cloud droplet number, and mixing of cloud water mixing ratio and aerosols are handled by the turbulence scheme.
It is worth to check if ndrop treats all MAM aerosols and cloud droplet number now, and CLUBB does all else."

[huiwanpnnl]

Thank you @adamrher for involving me. Indeed interested. Two caveats about my comments below: (1) I don't have access to CAM6 or newer CAM codes, so I can only speak about EAM. (2) My understanding of dropmixnuc is limited: from past conversations with Dick Easter and @kaizhangpnl, I got a rough idea about what dropmixnuc does but I have not spent time to understand all details of the equations and numerics.

What I see in dropmixnuc in EAM and heard from Dick Easter and @kaizhangpnl seems different from @XiZhao1013's notes about the old treatment, possibly because of version differences (evolution of code). In EAM,

• dropmixnuc calculates the changes of interstitial aerosol mixing ratios, cloud-borne aerosol mixing ratios, and stratiform cloud liquid number mixing ratio caused by turbulent mixing and aerosol activation/resuspension.
• The numerical method for solving the mixing-activation-resuspension equations is explicit, as the subroutine name explmix suggests.
• The turbulent diffusivity coefficient Kvh is retrieved from pbuf. The value in pbuf comes from CLUBB.

[huiwanpnnl]

In the spirit of unifying diffusion operations in CAM, can anyone see a problem with computing the diffusion coefficients in ndrop.F90, but applying them in CLUBB's diffusion solver? The only issue I could see is that CLUBB is located upstream of where ndrop is called in the time-loop (microp_aero_run I think).

@adamrher, I think dropmixnuc uses diffusion coefficient provided by CLUBB, so it is good that it is located downstream of CLUBB.

Note that dropmixnuc uses an explicit discretization method which is computationally inefficient. The solver in CLUBB is implicit hence allows for longer time steps. To unify turbulent mixing in CLUBB and dropmixnuc, one could

• separate tracer mixing from CLUBB and move it into dropmixnuc (i.e., merge with explmix but replace the explicit algorithm in explmix), or
• as @tilmes pointed out, bring the calculations of aerosol activation/resuspension into CLUBB and let CLUBB mix aerosols and cloud droplet number mixing ratio, too.

Also note that the aerosol activation/resuspension calculation needs information about cloud fraction change, which needs the liquid cloud fraction calculated by CLUBB. This results in the requirement that dropmixnuc be calculated downstream of CLUBB's PDF closure.

[JulioTBacmeister]

I don't think dropmixnuc should be in charge of diffusing all trace constituents. I would prefer this to be done somewhere else, maybe not even in CLUBB but in a separate vertical diffusion solver that accumulates diffusivities from all relevant processes.

[adamrher]

@huiwanpnnl thanks for providing your understanding. I also think that @XiZhao1013's notes are probably for an older version; I believe dropmixnuc is mixing all the interstitial and cloud-borne aerosols, as well as cloud liquid number. I'm trying to double check this with print statements in that module, but this is a work in progress. I can confirm that clubb intentionally does not mix all 19 aerosols + NUMLIQ for the specific compset I just ran ... So this is consistent with it being mixed in dropmixnuc.

@JulioTBacmeister says

I don't think dropmixnuc should be in charge of diffusing all trace constituents. I would prefer this to be done somewhere else, maybe not even in CLUBB but in a separate vertical diffusion solver that accumulates diffusivities from all relevant processes.

I'm trying to figure out if this is even possible. @tilmes statement:

Aerosols are changed from interstitial to cloud borne during the diffusion process

Is true. I'm wondering whether we can separate out the diffusion process, so it can be done elsewhere, from the activation/resuspension. For example, here is a typical call to explmix in dropmixnuc

CAM/src/physics/cam/ndrop.F90

Lines 1051 to 1059 in 0404b72

	call explmix( &
	raercol_cw(:,mm,nnew), source, ekkp, ekkm, overlapp, &
	overlapm, raercol_cw(:,mm,nsav), zero, zero, pver, &
	dtmix, .false.)
	call explmix( &
	raercol(:,mm,nnew), source, ekkp, ekkm, overlapp, &
	overlapm, raercol(:,mm,nsav), zero, flxconv, pver, &
	dtmix, .true., raercol_cw(:,mm,nsav))

I believe the source input argument is the tendency due to activation. Not sure about "resuspension" I'm out of my wheelhouse here. But nonetheless, this mixing routine seems to just diffuse raercol using CLUBB's KVH, while also adding the increment due to the srouce term. Could we instead just set the new state raercol = raercol + dt*source, and omit the call to explmix entirely, so that we can do this mixing somewhere else? Is this too crazy of an idea? I am thinking the diffusion tendencies would be rather small, and that activation is the more important process that probably should be updated within the CLUBB/MG subcycles, as the clouds evolve.