LILAC: coupling CTSM to atmospheric boundary layers schemes

In the current WRF-CTSM via LILAC coupling, we chose to use a popular, relatively-simple WRF boundary layer scheme (bl_pbl_physics == 1 or the so called YSU scheme). This scheme has a flux lower boundary condition similar to what is used in CAM. In addition to sensible and latent heat fluxes, this scheme also requires inputs that normally come from WRF's surface layer scheme, e.g., psim, psih, bulk Richardson number, roughness length, friction velocity. Most LSMs in WRF don't have multiple tiles and don't calculate their own surface exchange.

Sub-Issue 1: should we average these fields from the different tiles in CTSM and use those in the boundary layer scheme? If so, how to average since many are quite non-linear?

Sub-Issue 2: alternatively, we could call the surface layer scheme after CTSM using the bulk fluxes to create bulk estimates of these fields. To be consistent with CTSM approach(es) this may require accessing these CTSM modules.

Sub-Issue 3: WRF currently has at least 13 boundary layer scheme options, other atmospheric models will have additional options with different requirements. Coupling with other boundary layer schemes may require additional fields be passed out of CTSM to the atmosphere. For example, another popular WRF boundary layer scheme requires a flux-relevant surface temperature and exchange coefficient, so how to calculate average values of these?