Visibility

def UPPVIS(forecast):

    # Get maximum values between 1 and 2nd model levels
    value = {
        "QCLOUD" : np.max([forecast["qcloud"], forecast["qcloud", 1]], axis=0),
        "QVAPOR" : forecast["qvapor"],
        "QRAIN" : np.max([forecast["qrain"], forecast["qrain", 1]], axis=0),
        "QICE" : np.max([forecast["qice"], forecast["qice", 1]], axis=0),
        "QSNOW" : np.max([forecast["qsnow"], forecast["qsnow", 1]], axis=0),
        "rh" : np.max([forecast["rh"], forecast["rh", 1]], axis=0),
        "RHO" : forecast["rho"]
    }

    # Calculate the RUC visibility
    min = np.full((forecast.size.lat, forecast.size.lon), 0.8)
    zeros = np.full((forecast.size.lat, forecast.size.lon), 0.8)
    qrh = np.max( [ zeros, np.min( [ min, ( value["rh"] / 100. - 0.15 ) ], axis=0 ) ], axis=0 )
    visrh = 60. * np.exp( -2.5 * qrh )

    # Calculate SW99 visibility
    vis = SW99VIS( value["QVAPOR"], value["QCLOUD"], value["QRAIN"], value["QICE"], value["QSNOW"], value["RHO"] ) / 1000.

    # Return the minimum value
    return np.min([vis, visrh], axis=0)


# Calculates visibility by SW99 method
def SW99VIS(qv, qc, qr, qi, qs, rho, max = 24.135):
    """
    Calculates visibility based on the UPP algorithm.

    See documentation in UPPV2.2/src/unipost/CALVIS.f for the description of
    input arguments and references.
    """
    Rd = 287.
    COEFLC = 144.7
    COEFLP = 2.24
    COEFFC = 327.8
    COEFFP = 10.36
    EXPLC = 0.88
    EXPLP = 0.75
    EXPFC = 1.
    EXPFP = 0.7776

    #Tv = T * (1 + 0.61 * qv)  # Virtual temperature

    #rhoa = p / (Rd * Tv)  # Air density [kg m^-3]
    rhoa = rho
    rhow = 1e3  # Water density [kg m^-3]
    rhoi = 0.917e3  # Ice density [kg m^-3]

    vovmd = (1 + qv) / rhoa + (qc + qr) / rhow + (qi + qs) / rhoi

    conc_lc = 1e3 * qc / vovmd
    conc_lp = 1e3 * qr / vovmd
    conc_fc = 1e3 * qi / vovmd
    conc_fp = 1e3 * qs / vovmd

    # Make sure all concentrations are positive
    conc_lc[conc_lc < 0] = 0
    conc_lp[conc_lp < 0] = 0
    conc_fc[conc_fc < 0] = 0
    conc_fp[conc_fp < 0] = 0

    betav = COEFFC * conc_fc ** EXPFC \
            + COEFFP * conc_fp ** EXPFP \
            + COEFLC * conc_lc ** EXPLC \
            + COEFLP * conc_lp ** EXPLP + 1E-10

    vis = -numpy.log(0.02) / betav  # Visibility [km]
    vis[vis > max] = max

    return 1000 * vis