Changes in CLM elevation where SMB and Ts are computed (CESM_share portion)

[billsacks]

Moving this from the issue that @whlipscomb opened in CIME (ESMCI/cime#2809): The first part of the solution will be implemented in CMEPS (ESCOMP/CMEPS#251), and the second part here:

In JG/BG runs, Marcus Lofverstrom has found spatial and temporal discontinuities in the surface mass balance (SMB) and surface temperature (Ts) downscaled to CISM from CLM elevation classes (ECs). These discontinuities are related to the different treatment of real (nonzero-area) and virtual (zero-area) glacier columns in CLM. For real columns in a CLM grid cell, SMB and Ts are computed at the mean elevation of the ice-covered CISM cells in that EC within the domain of the CLM cell. For virtual columns, CLM computes SMB and Ts at the midpoint elevation of the EC.

To take a concrete example: Consider EC 9, which includes the range from 2500 to 3000 m. If a CLM cell has no ice in this range, it computes SMB and Ts at 2750 m. Now suppose the CISM topography is evolving in a BG simulation, and at some point there is a CISM cell at elevation 2999 m in the domain of the CLM cell. Then CLM abruptly introduces a real EC at 2999 m, replacing the virtual EC at 2750 m. This can result in abrupt changes in the values of SMB and Ts downscaled to EC 10 cells in CISM, because the linear interpolation in the coupler now uses a lower bound of 2999 m instead of 2750 m for EC 9. The changes are usually more pronounced for Ts than for SMB, since Ts has a stronger elevation dependence.

The proposed solution has two parts. First, in module map_glc2lnd_mod.F90, we smooth the transition from virtual to real elevations in CLM. Let h_k* = the elevation where CLM computes SMB and Ts for EC k. The new formula for this elevation is

h_k* = (a_0*H_k + a_k*h_k) / (a_0 + a_k),

where H_k is the fixed midpoint elevation of EC k; h_k is the mean elevation of CISM cells in EC k; a_k is the fractional ice area in EC k, a_0 is a prescribed constant. For a_0 = 0, this reduces to h_k* = h_k, which is the previous approach. When we set a_0 to a small nonzero value, ~ 0.05, the transition in h_k* takes place gradually as a_k increases from 0 to a finite fraction of the cell area.

Marcus has tested a_0 = 0.01, 0.05 and 0.10. He finds that 0.01 is too small to fix the problem, while results with 0.05 and 0.10 are similar. For now we will choose a default value a_0 = 0.05.

The second part of the solution, in module glc_elevclass_mod.F90, is to redefine the midpoint elevation H for EC 10. This value is arbitrary, since EC 10 has no upper bound. We have been using H_10 = 3500 m, but we will change to H_10 = 3001 m. Then when EC 10 changes from virtual to real, e.g. by thickening of CISM cells that were in EC 9 but now are slightly thicker than 3000 m, there will not be a sudden change in the elevation where CLM computes SMB and Ts.

In tests to date, these changes have helped smooth out temporal and spatial discontinuities in downscaled CISM fields. They are answer-changing in CLM, although the changes are expected to be modest.