Additional Modifications for negative snow

This issue continues work on issues related to negative amounts of snow discussed under FCM tickets [1396](https://code.metoffice.gov.uk/trac/jules/ticket/1396) and [1516](https://code.metoffice.gov.uk/trac/jules/ticket/1516).

The additional purely technical changes required to run with this code in the UM and LFRic are discussed under issues [um#140](https://github.com/MetOffice/um/issues/140) and [lfric_apps#797](https://github.com/MetOffice/lfric_apps/issues/797) respectively.

Under [1396](https://code.metoffice.gov.uk/trac/jules/ticket/1396), modifications were made to the code to address problems in the code that could lead to the generation of negative amounts of snow. These modifications were enabled using a logical, l_fix_neg_snow. In the original version of the code, changes to the amount of snow were calculated as a rate of change of melting or sublimation multiplied by the timestep. To avoid issues of numerical precision and guarantee that amounts of snow would not be negative, the code was reformulated to use increments to snow. It was technically infeasible to do this only when l_fix_neg_snow = .T., so bit-reproducibility was lost, even if l_fix_neg_now=.F.

The loss of bit-reproducibility obscured that fact that a bug had been introduced into sf_melt by this change. Additionally, it was subsequently discovered that on some compilers, not used in standard tests, numerical issues with amounts of snow close to underflowing could lead to crashes in UM runs. This issue also affected LFRic runs more generally. [1516](https://code.metoffice.gov.uk/trac/jules/ticket/1516) was created to address these issues.

An attempt to include this was made at version 7.6. Unfortunately, one test of bit comparison failed in a zonal mean diagnostic, even though the underlying 2D field did bit-compare. The change was not accepted at 7.6. (I believe this indicates that the code to calculate zonal means in the UM is not absolutely guaranteed to give identical results over different processor decompositions.)

At 7.9 the algorithm was very slightly reformulated to avoid the peculiar circumstances leading to this failure. This exposed a previously unidentified deficiency in the code whereby even the tiniest amounts of snow on a needle-leaved tree canopy can suppress evaporation of water. Because we could not be sure that fixing this would not have a systematic impact, it was decided not to include the branch at 7.9. Subsequently, other more urgent priorities prevented further work on this issue. Meanwhile, LFRic was being developed with a branch containing equivalent to those not included at 7.6 or 7.9. For the convenience of developing LFRic, it is desired to consolidate this branch on to the trunk.

We are therefore in the position of needing to support several options of the fix to consolidate the code being used by different systems. The purpose of this modification is to support these different options. l_fix_neg_snow will therefore be replaced by an integer option, i_fix_neg_snow that will take the following values:
* 0: This will maintain the behaviour of the trunk with l_fix_neg_snow=.F. (e.g. as used in standard UM jobs and UM AMIP runs)
* 1: This will be equivalent to l_fix_neg_snow=.F., but with the bug in sf_melt corrected. This is intended only as a historical convenience for people upgrading from older versions of JULES
* 2: This will maintain the behaviour of the trunk with l_fix_neg_snow=.T. (e.g. as used in standard standalone JULES jobs).
* 3: This will include the used currently used on a branch in LFRic and is what has been used in developing GC6.
* 4: This is an enhanced version of the previous option designed to remove the sensitivity of canopy evaporation to small amounts of snow. This is still developmental and may be modified at a future release.