Thanks to your suggestions on the block matrices, my fully coupled solver is proceeding very well! 

I am now about to take advantage of the block structure of the matrix using PCFIELDSPLIT. I have learned a bit from the user manual and followed with interest this discussion in the mailing list: https://lists.mcs.anl.gov/pipermail/petsc-users/2015-February/024154.html

which is actually the exact same situation I am in, so I guess most of the command line options will be copy and paste.

I would like however to code them in fortran, as I usually provide some default implementation alongside the command line options. 

While coding some of the options I got an error here in PCFieldSplitSetFields() which looks to be undefined. I am importing petscksp, do I need to import something else maybe? 

Thank you!