Outputting individual cloud microphysics process rates (source/sink terms) for Thompson and WSM6 in MPAS-Atmosphere

[X1a0wu]

Dear MPAS developers and users,

I am running MPAS-Atmosphere [v7.1] to simulate a [heavy snowfall event], and I would like to diagnose the cloud microphysics budgets of cloud ice, snow, and graupel in detail — for example, the separate contributions from deposition/sublimation, riming, freezing, and melting, similar to the microphysical budget analysis in Ko et al. (2022, GMD). I have been using config_microp_scheme = 'mp_thompson' (and I would also like to do the same with 'mp_wsm6').

So far, I have found the following fields in the Registry, which give the net tendency for each hydrometeor species:

• rqimpten, rqsmpten, rqgmpten (tendency of cloud ice/snow/graupel mixing ratio due to cloud microphysics)
• rthmpten (potential temperature tendency)
• rnimpten, rnrmpten (number concentration tendencies, Thompson only)

These are very useful, but they combine all source and sink processes into a single value. What I need is the individual process terms (e.g., deposition growth of snow, riming growth of graupel, melting of snow/graupel), such as those listed in Table 4 of Ko et al. (2022).

From reading the source code, I understand that these process rates are already computed internally in module_mp_thompson.F (e.g., prs_sde, prg_gde, prs_scw, prg_gcw, ...) and in module_mp_wsm6.F, but they are local variables and are not written to the model output.

Could you please advise:

1. Is there any existing option in MPAS-Atmosphere (namelist, stream field, or compile flag) to output these individual microphysical process rates, or has anyone already implemented such a diagnostic?
2. If not, what would be the recommended (minimal) way to add them? My current plan is to modify the Registry, the microphysics scheme module, the microphysics driver, and the stream file — is that the right approach, and are there any pitfalls to be aware of (e.g., the sub-stepping of the microphysics call)?

Thank you very much for your time.

Best regards

 

[Ming Chen]

Hi X1a0wu,

Please follow the steps below to add those process rates to your MPAS output files:

(1) in src/core_atmosphere/Registry.xml, add prs_sde, prg_gde, prs_scw, prg_gcw, etc. under the variable structure "tend_physics", i.e.,

<var_struct name="tend_physics" time_levs="1">

<var name="tend_prs_sde" type="real" dimensions="nVertLevels nEdges Time" units="kg m^{-2} s^{-1}"
description="Total horizontal momentum tendency from physics"
persistence="scratch" />

(2) in physics/physics_wrf/module_mp_thompson.F, get tend_prs_sde and add it as an output argument.

(3) in physics/mpas_atmphys_driver_microphysics.F, obtain tend_prs_sde from "call mp_gt_driver(...)"

Note that we need to pay attention to other details related to the newly added microphysics rates, e.g., initialization, array allocation, deallocation, etc.

It may be better to start with adding just one rate such as prs_sde. If it works, then we can add more rates.

Please try and let me know if you have any issues during the process. Thanks.

 

[X1a0wu]

Hi X1a0wu,

Please follow the steps below to add those process rates to your MPAS output files:

(1) in src/core_atmosphere/Registry.xml, add prs_sde, prg_gde, prs_scw, prg_gcw, etc. under the variable structure "tend_physics", i.e.,

<var_struct name="tend_physics" time_levs="1">

<var name="tend_prs_sde" type="real" dimensions="nVertLevels nEdges Time" units="kg m^{-2} s^{-1}"
description="Total horizontal momentum tendency from physics"
persistence="scratch" />

(2) in physics/physics_wrf/module_mp_thompson.F, get tend_prs_sde and add it as an output argument.

(3) in physics/mpas_atmphys_driver_microphysics.F, obtain tend_prs_sde from "call mp_gt_driver(...)"

Note that we need to pay attention to other details related to the newly added microphysics rates, e.g., initialization, array allocation, deallocation, etc.

It may be better to start with adding just one rate such as prs_sde. If it works, then we can add more rates.

Please try and let me know if you have any issues during the process. Thanks.

I have successfully resolved the issue. Thank you for your help. However, I encountered the following error while running WSM6:
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
atmosphere_model 0000000000ACB864 for__signal_handl Unknown Unknown
libpthread-2.17.s 00002AE86B9D65D0 Unknown Unknown Unknown
atmosphere_model 000000000067FFD0 Unknown Unknown Unknown
atmosphere_model 0000000000537093 Unknown Unknown Unknown
atmosphere_model 00000000004C9FDA Unknown Unknown Unknown
atmosphere_model 00000000004C0D6F Unknown Unknown Unknown
atmosphere_model 00000000004A3F59 Unknown Unknown Unknown
atmosphere_model 00000000004127E7 Unknown Unknown Unknown
atmosphere_model 0000000000412757 Unknown Unknown Unknown
atmosphere_model 00000000004126DE Unknown Unknown Unknown
libc-2.17.so 00002AE86BF073D5 __libc_start_main Unknown Unknown
atmosphere_model 00000000004125E9 Unknown Unknown Unknown Could it be that my time step is too large / not small enough? I set it to 15 s for the 60–3 km configuration. Or is there a problem with the namelist settings（
config_physics_suite = 'convection_permitting'
config_microp_scheme = 'mp_wsm6'
）? Or is it because the memory on the nodes I requested is insufficient?

 

[Ming Chen]

@ X1a0wu,

How long did the case run before it crashed? If you run the same case but using the un-modified code, can it run successfully?

Please send me your namelist.atmosphere and streams.atmosphere to take a look. Thanks.

 

[X1a0wu]

@ X1a0wu,

How long did the case run before it crashed? If you run the same case but using the un-modified code, can it run successfully?

Please send me your namelist.atmosphere and streams.atmosphere to take a look. Thanks.

Hello,

Thank you for your previous help. I have successfully run the simulation. My current question is: does WSM6 require different modifications compared to Thompson for outputting microphysical rate variables?

I am attempting to add 9 3D microphysics process rate diagnostic variables (e.g., snow deposition/sublimation, snow riming, graupel melting, etc.) that are output for both the Thompson and WSM6 schemes using unified Registry variable names.

What I did:

1. Registry.xml: Registered 9 diagnostic variables (e.g., tend_prs_sde, tend_prs_scw, ...) in the tend var_struct with dimensions nVertLevels nCells Time.
2. mp_wsm6.F: Added 9 output parameters (2D arrays, intent(out), dimension its:ite × kts:kte) to the mp_wsm6_run subroutine signature. At the end of the subroutine (after all process rate calculations), I added a copy loop:

do k = kts, kte<br> do i = its, ite<br> prs_sde_out_2d(i,k) = psdep(i,k) - psevp(i,k)<br> prs_scw_out_2d(i,k) = psacw(i,k)<br> prs_iau_out_2d(i,k) = psaut(i,k)<br> pri_ide_out_2d(i,k) = pidep(i,k)<br> prs_sci_out_2d(i,k) = psaci(i,k)<br> prr_sml_out_2d(i,k) = psmlt(i,k)<br> prg_gcw_out_2d(i,k) = pgacw(i,k)<br> prg_rcs_out_2d(i,k) = pgacs(i,k)<br> prr_gml_out_2d(i,k) = pgmlt(i,k)<br> enddo<br>enddo

1. module_mp_wsm6.F (MPAS wrapper): Added 9 3D output parameters to the wsm6 subroutine, declared 9 local 2D work arrays, passed them to mp_wsm6_run via keyword arguments, then copied the 2D results back to the 3D arrays ((ims:ime, kms:kme, jms:jme)).
2. MPAS driver/interface (mpas_atmphys_driver_microphysics.F, mpas_atmphys_interface.F, mpas_atmphys_vars.F): Allocated the 9 local 3D arrays, passed them to wsm6 as keyword arguments, and copied the results back to the MPAS pool fields via mpas_pool_get_array(tend, 'tend_prs_sde', ...).

Results:

• Thompson experiment: All 9 variables output correctly with non-zero values. ✓
• WSM6 experiment: The model runs to completion without crashes or NaN, but all 9 process rate variables are identically zero everywhere in the domain and for all output times. ✗

Verification I have done:

• Confirmed the modified mp_wsm6.F and module_mp_wsm6.F are on the server (checked with grep).
• Confirmed the atmosphere_model executable was compiled AFTER the source changes (checked file timestamps).
• The Thompson experiment using the same Registry, interface, and pool mechanism works correctly, so the Registry/pool/interface setup is functional.

Questions:

1. Is there something structurally different about how WSM6 computes its process rates that would cause them to remain zero even in a snowstorm case?
2. Are the internal process rate variables (psdep, psacw, psmlt, etc.) in mp_wsm6_run populated during normal execution, or are they computed only under certain conditional branches that might not be triggered?
3. Is there a known issue or recommended approach for exposing WSM6 process rate diagnostics in MPAS?

Any advice would be greatly appreciated.

Thank you."