OpenMP run failure: subroutine atm_rk_integration_setup in mpas_atm_time_integration.F

[MicroTed]

I was experimenting using OpenMP with an ideal periodic domain and was getting random instabilities in W right at the start of the second OMP thread. (e.g., at column 451 of 900 grid points). This suggested something being used before getting initialized, and it appears to stem from a thread getting past atm_rk_integration_setup and into atm_compute_moist_coefficients before the other OMP thread. (atm_compute_moist_coefficients uses values that are initialized in atm_rk_integration_setup.)

Two possible solutions (mpas_atm_time_integration.F):

1. Don't thread atm_rk_integration_setup at all because the only work is copying arrays, which like does not benefit from threading at all.
2. Put an OMP BARRIER at the end of atm_rk_integration_setup to prevent any threads from getting ahead into the next routine. (Or alternatively add a barrier at the beginning of atm_compute_moist_coefficients.)

It doesn't seem to work to put a barrier in the calling routine (atm_srk3) between the two calls because it isn't within a parallel section.

Details: MacOS with gfortran 12.4, openmpi/4.1.2, hdf5/1.10.6, netcdf/4.7.4, pio/2.5.7
Openmp-c libomp from homebrew:
"CFLAGS_OMP = -Xpreprocessor -fopenmp -I/opt/homebrew/Cellar/libomp/21.1.8/include -L/opt/homebrew/opt/llvm/lib -lomp"

Compile-time options:
Build target: gnu
OpenMP support: yes
OpenACC support: no
Default real precision: single
Compiler flags: optimize
I/O layer: PIO 2.x

Run-time settings:
MPI task count: 1
OpenMP max threads: 2

Note: I have not checked bit-for-bit results with and without omp threading. As I recall, the bug did not show up with normal debug options until I add -O2. I thought that maybe the mpas_timer_stop/start timer calls would sync the threads (act as a barrier), but apparently not?

[MicroTed]

It appears than a single omp_barrier at the end of atm_rk_integration_setup is not sufficient. I had to add a barrier at the beginning as well. (Or perhaps only at the beginning might work, but I haven't tested that yet.)

Additionally, the results with different numbers of OMP threads do not match. (In constrast to matching results with 1 or more MPI tasks). The attached archive has an input condition (idealized, periodic) with some land surface settings for running full physics. It has a weak thermal pert to initiate small vertical motions. The namelist turns on radiation (rrtmg), pbl (ysu), land (noahmp), microphysics (thompson).

In addition to not matching, the boundary layer evolves faster as the number of OMP threads increases. Cloud water appears earlier and earlier. The attached image shows cloud water at the same time for 1, 2, and 4 threads (left to right). It seems to act like some rate might be getting applied globally by each thread, but that is just speculation.

One more oddity is that most timings are reduced by adding threads, except for the monotonic scalar advection, which actually increases:
atm_advance_scalars_mono 19.17885 (omp=1) 27.36747 (omp=2) 33.08643 (omp=3)

I don't see anything obvious in the code to cause this, but I haven't really scrutinized it yet. (A test of 1 vs. 4 threads with mono turned off had the same kind of differences.)

This is using the current develop branch version (commit 447de92, Dec. 31 2025) with just the one change to add barriers in atm_rk_integration_setup.

mpas_omp_test.tgz