introducing P-CSI solver

[lj-mu]

What changes does this PR introduce?

New feature. This P-CSI solver (implemented in CESM, Huang et al., 2016, doi:10.5194/gmd-9-4209-2016) is a potentially faster alternative to the CG2D solver.

What is the current behaviour?

At very high resolution S/R CG2D can be a bottleneck (see #530)

What is the new behaviour

• the P-CSI solver can be faster by a factor of 2 in specific configurations and architectures
• the P-CSI solver needs global communications only for checking convergence
• The code is compiled with PCSI by defining ALLOW_PCSI in the CPP_OPTIONS.h.
• To use, set usePcsiSolver=.TRUE. and select values for PcsiMaxIters and PcsiTargetResidual in data
• new flag useCG2DPrecond=.TRUE. by default to use the full cg2d-preconditioner (also for P-CSI), set to .FALSE. to only use the main diagonal (for tuning and testing)
• global_ocean.90x40x15.idemix tests the new feature

Does this PR introduce a breaking change?

No

Other information:

• to fool testreport, solve_for_pressure prints pcsi_init_res as cg2d_init_res
• Yuhu Chen (@xiaohumengdie) contributed most of the code

Suggested addition to tag-index

• new P-CSI solver for solve_for_pressure, to compile define ALLOW_PCSI, to use set usePcsiSolver=.TRUE. and select values for PcsiMaxIters and PcsiTargetResidual in namelist data&PARM02