Source code for pycif.plugins.domains.chimere.utils.make_vcoord

import numpy as np
from .....utils.path import init_dir




def make_vcoord(domain):

    # Pressures in hPa
    nlev = domain.nlev
    pmax = domain.pmax
    pref = 1000
    p1 = domain.p1

    dp1 = pref - p1
    a = (pref - pmax) / (pref - p1)

    # Add or not a buffer layer
    sigma_a = []
    sigma_b = []
    if p1 / pref < domain.buffer_layer < 1:
        sigma_b = [domain.buffer_layer] + sigma_b
        sigma_a = [0] + sigma_a
        nlev -= 1

    # Scaling of thickness
    x = 1.01
    for i in range(30):
        f = 1 - a + np.sum(x ** np.arange(1, nlev))
        fp = 1 + np.sum((1 + np.arange(2, nlev)) * x ** np.arange(1, nlev - 1))
        x -= f / fp

    # Generate sigma_a and sigma_b
    p = pref
    dp = dp1
    for i in range(nlev):
        p -= dp
        dp *= x
        ca = max(0, (pmax / pref) * (p1 - p) / (p1 - pmax))
        cb = max(0, (p1 / pref) * (p - pmax) / (p1 - pmax))
        sigma_a.append(ca)
        sigma_b.append(cb)

    # Changing sigma_a to traditional CHIMERE unit
    sigma_a = np.array(sigma_a)  # / 1000

    # Dump to VCOORD file
    init_dir(f"{domain.workdir}/domain/VCOORD/")
    np.savetxt(
        f"{domain.workdir}/domain/VCOORD/VCOORD_{domain.nlev}_{p1}_{pmax}",
        np.array([sigma_a, sigma_b]).T,
        fmt="%10.5f")

    # Add sigma_a and sigma_b to domain
    domain.sigma_a = np.array(sigma_a)
    domain.sigma_b = np.array(sigma_b)