[paulfield2024]

PR Summary

Sci/Tech Reviewer: Anthony Jones (@ajones244)
Code Reviewer: James Bruten (@james-bruten-mo)

• linked MetOffice/um#119

Code Quality Checklist

(Some checks are automatically carried out via the CI pipeline)

• I have performed a self-review of my own code
• My code follows the project's style guidelines
• Comments have been included that aid understanding and enhance the readability of the code
• My changes generate no new warnings

Testing

• If shared files have been modified, I have run the UM and LFRic Apps rose
  stem suites
• If any tests fail (rose-stem or CI) the reason is understood and
  acceptable (eg. kgo changes)
• I have added tests to cover new functionality as appropriate (eg. system
  tests, unit tests, etc.)

trac.log

Security Considerations

• I have reviewed my changes for potential security issues
• Sensitive data is properly handled (if applicable)
• Authentication and authorisation are properly implemented (if applicable)

Performance Impact

• Performance of the code has been considered and, if applicable, suitable performance measurements have been conducted

AI Assistance and Attribution

• Some of the content of this change has been produced with the assistance of Generative AI tool name (e.g., Met Office Github Copilot Enterprise, Github Copilot Personal, ChatGPT GPT-4, etc) and I have followed the Simulation Systems AI policy (including attribution labels)

Documentation

• Where appropriate I have updated documentation related to this change and confirmed that it builds correctly

Sci/Tech Review

• I understand this area of code and the changes being added
• The proposed changes correspond to the pull request description
• Documentation is sufficient (do documentation papers need updating)
• Sufficient testing has been completed

Please alert the code reviewer via a tag when you have approved the SR

Code Review

• All dependencies have been resolved
• Related Issues have been properly linked and addressed
• CLA compliance has been confirmed
• Code quality standards have been met
• Tests are adequate and have passed
• Documentation is complete and accurate
• Security considerations have been addressed
• Performance impact is acceptable

[paulfield2024]

This PR is linked to UM PR

Summary of changes on aerosol_interaction_vn14.2

Summary of differences: branch aerosol_interaction_vn14.2 vs main

main = a2066e6 ("Update GitHub workflows and templates (#18)")
HEAD = 5ef7601
Line numbers refer to the current (branch) version of each file.

racw -- src/accretion.F90

175-179 Aerosol transfer on rain-accretes-cloud rewritten.
Was: damass = dmass/cloud_mass*aerofields(k,i_am4)
Now: guarded by if (dmass /= 0.0_wp .and. qfields(k,i_ql) > 0.0_wp)
and divides by qfields(k,i_ql) instead of cloud_mass.
Reason: cloud_mass was the wrong/unset quantity here; the divide also
needed protecting against zero or negative cloud water.

activate -- src/activation.F90

206 smax_cloud = alphawtau -> MIN(0.10, alphawtau)
Reason: cap in-cloud supersaturation at 10% so large updraughts
cannot produce unphysical activation.

264-266 dactive = zero_real_wp / dmass_d = zero_real_wp moved out so they are
always initialised before the Smax_act test.
Reason: they could previously be used uninitialised.

269 Branch changed from if (iopt_inuc < 4) ... else to
if (iopt_inuc > 3) then
Reason: only the higher-order ice nucleation schemes now need the
dustphys-based values; initialisation is handled above.

examine_aerosol -- src/aerosol_routines.F90

18-19 Two blank lines added after the use statements (cosmetic).

214-215 aero_mact_mean_max: 2.3e-15 -> 1e30
dust_mact_mean_max: 3.1e-13 -> 1e30
Reason: the old limits were too aggressive for coarse soluble
aerosol; effectively disables the cap.

320-324 aeroact mean-mass clip replaced by:
if (aeroact(k)%mact_mean .gt. aero_mact_mean_max) then
aeroact(k)%mact_mean = min(...)
aeroact(k)%nact = aeroact(k)%mact/aero_mact_mean_max
ntot = aeroact(k)%nact
endif
Reason: previously only the mean mass was clipped, leaving mass and
number inconsistent; number is now recomputed to match.

463-467 Same change applied to dustact (nitot recomputed).
566-570 Same change applied to aeroice (nitot recomputed).
657-661 Same change applied to dustliq (ntot recomputed).

(deleted) The per-mode clipping lines
aeroact(k)%mact2_mean=min(...) aeroact(k)%mact1_mean=min(...)
dustact(k)%mact2_mean=min(...) dustact(k)%mact3_mean=min(...)
dustact(k)%mact1_mean=min(...) aeroice(k)%mact2_mean=min(...)
were removed (6 lines, previously around old lines 343, 351, 467,
476, 485, 569).
Reason: redundant now the cap is effectively unlimited and the
mass/number consistency fix is applied to the totals instead.

revp -- src/evaporation.F90

152 dmac=min(dmac,aeroact(k)%mact2/dt)
-> dmac=min(dmac,max(aeroact(k)%mact2,0.0_wp)/dt)
Reason: a negative mact2 would otherwise drive dmac negative; floor
the available mass at zero.

189-194 New local dnumber_a = dnumberaeroact(k)%nratio2, plus
if (dnumber_admac <= 0.0) then dnumber_a=0.0; dmac=0.0
before the which_mode call.
Reason: when mass and number changes have opposite signs the mode
transfer is ill-defined and small, so it is skipped.

218 Added dmacd=min(dmacd,max(dustliq(k)%mact2,0.0_wp)/dt)
Reason: same availability limiter for dust-in-liquid mass.

ihom_rain -- src/homogeneous_freezing.F90

160 dmac = dnumberaeroact(k)%mact2_mean
-> dnumberaeroact(k)%mact2_mean*aeroact(k)%nratio2
Reason: missing nratio2 factor; only the fraction of rain drops
containing this aerosol should be counted.
161 Added dmac=min(dmac,aeroact(k)%mact2/dt)
168 Added dmadl=min(dmadl,dustliq(k)%mact2/dt)
Reason: cannot remove more aerosol/dust than is available in dt.

ihom_droplets -- src/homogeneous_freezing.F90

300 dmac = dnumberaeroact(k)%mact1_mean
-> dnumberaeroact(k)%mact1_mean*aeroact(k)%nratio1
Reason: missing nratio1 factor (as above, for cloud droplets).
301 Added dmac=min(dmac,aeroact(k)%mact1/dt)
308 Added dmadl=min(dmadl,dustliq(k)%mact1/dt)
Reason: mass availability limiter.

iacc -- src/ice_accretion.F90

441 dmac=min(dmac,aeroact(k)%mact1/dt) (cloud branch, Y)
443 dmad=min(dmad,dustliq(k)%mact1/dt)
446 dmac=min(dmac,aeroact(k)%mact2/dt) (rain branch, Y)
448 dmad=min(dmad,dustliq(k)%mact2/dt)
451 dmac=min(dmac,aeroact(k)%mact1/dt) (cloud branch, X)
453 dmad=min(dmad,dustliq(k)%mact1/dt)
456 dmac=min(dmac,aeroact(k)%mact2/dt) (rain branch, X)
458 dmad=min(dmad,dustliq(k)%mact2/dt)
Reason: limit aerosol/dust mass moved by ice accretion to what is
actually present in the relevant mode during the timestep.

idep -- src/ice_deposition.F90

207 Added dmad=max(dmad,-dustact(k)%mact1/dt) (i_dsub)
211 Added dmad=max(dmad,-dustact(k)%mact2/dt) (i_dssub)
215 Added dmad=max(dmad,-dustact(k)%mact3/dt) (i_dgsub)
Reason: dust sublimation must not remove more than is available.

236 dmac=min(dmac,aeroice(k)%mact1/dt)
-> dmac=max(dmac,-aeroice(k)%mact1/dt)
240 dmac=min(dmac,aeroice(k)%mact2/dt)
-> dmac=max(dmac,-aeroice(k)%mact2/dt)
244 dmac=min(dmac,aeroice(k)%mact3/dt)
-> dmac=max(dmac,-aeroice(k)%mact3/dt)
Reason: bug fix - this block only runs when sublimating, so dmac is
negative; a min() against a positive limit was wrong and the limiter
must be a max() against the negative of the available mass.

inuc -- src/ice_nucleation.F90

390 Added dmac=min(dmac,aeroact(k)%mact1/dt)
394 Added dmadl=min(dmadl,dustliq(k)%mact1/dt)
Reason: mass availability limiter for frozen soluble aerosol and for
dust already in the liquid phase.

mphys_switches (module data) -- src/mphys_switches.F90

99 activate_in_cloud default 2 -> 1
Reason: revert to plain ARG in-cloud activation rather than the
"smaller of ARG smax and existing-droplet smax" option.

microphysics_common -- src/micro_main.F90

1522-1523 ensure_positive (cloud) iprocs_dependent extended from
(/i_aaut, i_aacw/) to
(/i_aaut, i_aacw, i_diacw, i_dsacw, i_dgacw, i_dhomc, i_dnuc/)

1528-1530 ensure_positive (ice) now passes
aeroprocs=aerosol_procs and
iprocs_dependent=(/i_draci, i_dsub, i_dimlt/)

1537 ensure_positive (rain) iprocs_dependent extended from
(/i_arevp/) to (/i_arevp, i_dsacr, i_dgacr, i_dhomr/)

1542-1544 ensure_positive (snow) now passes
aeroprocs=aerosol_procs and
iprocs_dependent=(/i_dsmlt, i_dsacr, i_dssub/)
Reason (all four): when a hydrometeor field is rescaled to keep it
positive, the dependent aerosol/dust processes must be rescaled by
the same factor, otherwise aerosol and hydrometeor budgets become
inconsistent.

1932-1936 ensure_positive_aerosol call for (/i_asedr, i_asedl/) commented out.
1944-1948 ensure_positive_aerosol call for
(/i_asedr, i_asedl, i_dsedi, i_dseds, i_dsedg/) commented out.
Reason: "No need to call ensure positive aerosol. Limiter added to
sedimentation" - positivity is now enforced inside sedr_1M_2M.

mphys_tidy (module use statements) -- src/mphys_tidy.F90

5-6 use thresholds now also imports aeromass_small, aeronumber_small.
18-19 use mphys_switches now also imports i_an1, i_am1, i_an2, i_am2,
i_an3, i_am3, i_am6, i_an6, i_am10, i_an10.
Reason: needed by the new small-value tidying block in tidy_ain.

qtidy -- src/mphys_tidy.F90

214-215 l_qsig(0)=.false. / l_qsneg(0)=.false.
-> l_qsig(:)=.false. / l_qsneg(:)=.false.
Reason: only element 0 was being initialised; the whole array must be
reset each call.

219 Added qr_reset=.false.
221 Added nl_reset=.false.
222 Added ql_reset=.false.
223 Added ni_reset=.false.
Reason: these reset flags were not being cleared at the top of the
k-loop and could carry stale values from the previous level.

tidy_ain -- src/mphys_tidy.F90

838 Blank line inserted (cosmetic).
841 When am4 is zeroed, its mass is added to am2 first.
846 When am9 is zeroed, its mass is added to am6 first.
853 When am5 is zeroed, its mass is added to am2 first.
860 When am7 is zeroed, its mass is added to am6 first.
867 When am8 is zeroed, its mass is added to am2 first.
Reason: previously the removed activated-aerosol/dust mass was simply
discarded; it is now returned to the corresponding unactivated
(accumulation / coarse-dust) mode so total aerosol mass is conserved.

880 Test condition aerofields(k,i_an12) > 0.0 -> aerofields(k,i_an11) > 0.0
Reason: copy/paste bug - the i_an11 branch was testing i_an12.

886-963 New block (78 lines) that zeroes each aerosol mode when its mass is
below aeromass_small or its number below aeronumber_small, for modes
am1/an1, am2/an2, am3/an3, am6/an6, am10/an10.
For the accumulation soluble mode (am2/an2) the mass is transferred
into am4 before zeroing (lines 907 and 939).
Reason: remove residual tiny/inconsistent mass-number pairs that
cause instability, while conserving the accumulation-mode mass.

ensure_positive_aerosol -- src/mphys_tidy.F90

1294-1295 Added declarations real(wp) :: drain and
real(wp) :: ratio_field(ntotala)
1330 Blank line added.
Reason: leftovers from development; both variables are currently
unused and could be removed. The routine body itself is unchanged
relative to main.

sedr_1M_2M -- src/sedimentation.F90

948-951, 957-960, 967-970, 976-979 (interior levels: cloud, rain)
988-995, 1005-1012, 1022-1029 (interior levels: ice, snow, graupel)
1050-1051, 1056-1057, 1063-1064,
1069-1070, 1078-1081, 1089-1092,
1101-1104 (top level, all species)

     Every dmac / dmad aerosol- and dust-mass flux-divergence expression
     was changed from a plain flux difference, e.g.
       dmac = (flux_n2(k+1)*nratio1*mact1_mean
              -flux_n2(k)  *nratio1*mact1_mean) * rdz_on_rho(k,...)
     to capping each level's contribution first:
       dmac = ( min(flux_n2(k+1)*nratio1*mact1_mean,
                    mact1(k+1)/(step_length*rdz_on_rho(k+1,...)))
               -min(flux_n2(k)  *nratio1*mact1_mean,
                    mact1(k)  /(step_length*rdz_on_rho(k,...))) )
              * rdz_on_rho(k,...)
     Reason: the sedimentation flux out of a level could exceed the
     aerosol/dust mass present there, driving fields negative. This is
     the limiter that replaces the ensure_positive_aerosol calls removed
     in micro_main.

casim_erfc -- src/special.F90

243 casim_erfc = erfg(x,c) -> casim_erfc = erfc(x)
Reason: bug in erfg for c=1 and x<1e-10; use the intrinsic erfc.

which_mode_to_use (module parameter) -- src/which_mode_to_use.F90

15 imethod = 2 -> 1
Reason: revert to the old method rather than the simple method.

which_mode -- src/which_mode_to_use.F90

83-87 r1 = min(max_accumulation_mean_radius, r1_in)
-> IF (r1_in > 0.0) THEN
r1 = min(max_accumulation_mean_radius, r1_in)
ELSE
r1 = max_accumulation_mean_radius
END IF
Reason: a non-positive r1_in would previously be passed straight
through, giving a zero/negative accumulation-mode radius downstream.

Other

CONTRIBUTORS.md 1 line added (contributor name).

[Anthony Jones (ajones244)]

Very small changes CASIM side. These are mostly bugfixes for a non-operational setup (aerosol-cloud coupling). I've followed the logic and it makes sense. There is a decision to revert to default ARG activation.

[Anthony Jones (ajones244)]

Will this change results in any operational configuration e.g. RAL3?

[Anthony Jones (ajones244)]

Just to comment, smax_cloud is the maximum supersaturation, we would not expect this to be greater than 10% (over marine regions, this is ~5%)

[Anthony Jones (ajones244)]

I think this can be removed, the only call to umMessage is commented out

[Anthony Jones (ajones244)]

suggest removing these commented out lines to improve readability?

[Anthony Jones (ajones244)]

This reformulates the mass conservation by calculating the ratio of total aerosol / aerosol lost or gained and then ensuring this lies between 0 and 1 (to avoid negative aerosol) and scaling accordingly

[Anthony Jones (ajones244)]

This logic is also on L353, should that be changed too?:
dmac=(flux_n2(k+1)aeroact(k+1)%nratio1aeroact(k+1)%mact1_mean - &
flux_n2(k)aeroact(k)%nratio1aeroact(k)%mact1_mean)* rdz_on_rho(k,ixy_inner)

[Anthony Jones (ajones244)]

NB: potentially this is another unused module and not sedr_1M_2M - consider if useful to copy the logic.

[Anthony Jones (ajones244)]

Test of units:
flux_n2 * nratio1 * mact1 = kg.m-2.s-1 * kg.kg-1 = kg.m-2.s-1
mact1 / (step*rdzontho) = kg.kg-1 * s-1 * kg.m-2= kg.m-2.s-1
So units wise all good. The test is that the flux of aerosol coming from above does not exceed the total mass above. I think this is fine. dmad is then the net incoming flux from above minus outgoing below.

[Anthony Jones (ajones244)]

Simple bug fix to avoid erroneous values of erf as commented

[Anthony Jones (ajones244)]

Simple check for reasonable r1, also reverts the default method

[Anthony Jones (ajones244)]

Hi Paul,
I have reviewed the code and I am satisfied that the changes are appropriate, mostly amounting to bug fixes, slight logic changes and conservation improvements.
To reiterate, the purpose of this ticket is to improve aerosol conservation when SOL/INSOL is coupled to CASIM, alongside associated bug fixes.
KGO for the SOL/INSOL-CASIM tests change because of the aerosol changes and the change of ndrop_pot (m-3 now rather than kg-1, which is what the rad_ctl wants)
Rose stem passes for the operational / higher-priority tests.
Passing to James Bruten (@james-bruten-mo) for CR. Happy to discuss this again next week if more is needed.
There are a larger number of scientific configuration improvements on the CASIM side which are meant to improve the aerosol conservation. You have demonstrated in offline tests (shown to me) that this conservation requires significant changes to the logic. I have hence reviewed the logic in detail to ensure the changes are consistent and am happy to approve as Sci tech reviewer)
Anthony

[James Bruten (james-bruten-mo)]

Unless there's a good reason to keep them (ie. imminently, they will be reinstated in a separate PR) it'd be better to keep the code tidy by removing the commented out sections.
Otherwise looks good

[James Bruten (james-bruten-mo)]

Is there a good reason to keep this code but commented out? Generally we prefer to not keep commented out code blocks

[paulfield2024]

Hello James - thanks - i've removed those bits in micro_main and mphys_tidy.

paulfield2024@0b85f9a

Please let me know if i missed anything.
thanks

[James Bruten (james-bruten-mo)]

Same comment as above

[paulfield2024]

Hello James - thanks - i've removed those bits in micro_main and mphys_tidy.

paulfield2024@0b85f9a

Please let me know if i missed anything.
thanks

[James Bruten (james-bruten-mo)]

Thanks Paul, looks good now