https://gcc.gnu.org/bugzilla/show_bug.cgi?id=125568
--- Comment #3 from Ivan Pribec <ivan.pribec at gmail dot com> ---
Thanks for the pointer. I can workaround with either of these variants that use
the product to fold the inner dimension:
! Variant A
integer, parameter :: c(0:nmax, 0:nmax) = &
product(reshape( &
[(((merge(j-ic+1, 1, ic<=k), ic=1, nmax), j=0, nmax), k=0, nmax)], &
[nmax, nmax+1, nmax+1]), dim=1)
! Variant B
integer, parameter :: c(0:nmax, 0:nmax) = &
reshape(product(reshape( &
(/( (/( [(ic, ic=j-k+1, j), (1, ic=1, nmax-k)], &
j=0, nmax) /), k=0, nmax) /), &
[nmax, (nmax+1)*(nmax+1)]), dim=1), &
shape=[nmax+1, nmax+1])