Memory ordering solver is not optimal for any form of 3D radiation

[Chiil]

The solver is executed in blocks of 16 n_col_block (hard coded value resulting from performance tuning in CPU version, see for instance in src_test/Radiation_solve.cxx:378) columns and then for all layers and gpoints. This setup is chosen, because in the RRTMGP part of the code, ngpt is the inner dimension, in order to have the most efficient table lookups. In the end of the RRTMGP solver the data is transposed, in order to make ngpt the outermost dimension (for instance src/Gas_optics_rrtmgp.cpp:1280). This means that the code constantly works with fields that are of the order of ngpt*nlay*n_col_block in size.

For a solver that could be extended to 3D (slanted column, n-stream, ray tracing) it is essential to have the full 3D field in memory, so ngpt should be the outer dimension under all conditions. In this case, the full solver could be run per gpoint.