Create new module supporting ClimaOcean and ClimaLand use cases (WIP)

[kmdeck]

The Climate Modeling Alliance

Software Design Issue 📜

Purpose

1. Clima ocean has a parameterization that requires that the roughness lengths to be a function of u_star, q_star, \theta_star. They are not constant. The current SF.jl module assumes these are constant, and then solves for u_star, q_star, \theta_star.
2. For land, the "surface specific humidity" is not known/modeled. The same is true for the "surface temperature". In general, these are a combination of the states of the land model and must be parameterized. Currently, we pass in q_s = q_sat(T_component) for each component for the land model, and adjust the fluxes after the MOST solve to account for deviations in q_s from q_sat and in T_s from T_component. Note that this implies that the fluxes we use are not a solution of the MOST equations. Moreover, fluxes from each component must be computed in parallel. A more general solution is detailed here: https://www.overleaf.com/8479341197nzfdmtznnrgz#ed793a, where we (on the land side) parameterize q_s and T_s directly. These are typically of the form: q_s = r_land q_star u_star + q_sat. Therefore the evauation of q_s and T_s must be accounted for in the MOST solve (updated each iteration/new guess for u_star, q_star, theta_star). The nice thing about this formulation for land is that it allows us to do a lot of use cases with a single interface.

In summary, the purpose is to extend the MOST solve to let roughness lengths and surface state vary as a function of the solver state u_star, q_star, theta_star.

Cost/Benefits/Risks

costs: developer time and testing time benefits: allows land to compute parallel surface fluxes (atmos-land component) for each component with a single interface, allows land to compute atmos- composite land surface fluxes with the same interface. allows ocean to compute SF with the roughness parameterization they like, allows ocean and land to use the same SF interface risks: the solver is very slow to converge

People and Personnel

@akshaysridhar @kmdeck @glwagner @simone-silvestri

Components

New module in SF.jl which allows the user to pass in functions for roughness lengths and surface states. These functions could be of the form f(surface_properties, solver_state, atmos_state), where for land, surface_state would include resistances, q_sat(T_land component), while for ocean, surface_properties could contain whatever is needed to compute z_0.

Note that: the land bucket model uses q_a in estimated q_s, so I included the atmos_state as an input. If we only support the bucket model using the current SF module, this function can be instead: f(surface_properties, solver_state)

Inputs

Equations to solve from ocean and land

Results and Deliverables

Ocean and land using a single stable SF.jl package in offline (forced ERA5/reanalysis) simulations (OMIP/LMIP). Current bucket AMIP simulations unaffected by change

Scope of Work

### Tasks
- [ ] Complete design
- [ ] Write code (single PR in SF.jl adding a new module)
- [ ] Test convergence speed using existing SF.jl tests
- [ ] Ocean team tests using the new module from their repo
- [ ] Land team tests using the new module from their repo

SDI Revision Log

CC

@tapios @sriharshakandala @charleskawczynski @cmbengue

[glwagner]

I also want to point out that an evaluation of the surface temperature is also required for sea ice and is well-justified for the ocean as well. For example Fairall et al. 1996 showed that skin temperature errors can produce errors of ~ 10 W / m^2 in the latent heat flux (see also Zeng and Beljaars 2005). So I think it will be generally useful (and also is not very difficult) to be able to incorporate submodels for the surface state (as distinct from the prognostic interior state of the land/ocean model) into the flux solver.

[glwagner]

Here's are some links to ClimaOcean's implementation of a flux solver based on fixed point iteration:

https://github.com/CliMA/ClimaOcean.jl/blob/a4c405f6a731fadef33f49ac88e756c048e6f43b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl#L231

https://github.com/CliMA/ClimaOcean.jl/blob/a4c405f6a731fadef33f49ac88e756c048e6f43b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl#L267

https://github.com/CliMA/ClimaOcean.jl/blob/a4c405f6a731fadef33f49ac88e756c048e6f43b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl#L387

I think it makes sense to start with a solver based on fixed point iteration, which is very simple, and consider more complicated / efficient solvers once the design is matured a little bit and we have evaluated its relative performance.