def pseudocoral(lat, lon, sst, sss=None, d18O=None, species='default', b1=
    0.3007062, b2=0.2619054, b3=0.436509, b4=0.1552032, b5=0.15):
    """
    DOCSTRING: Function 'pseudocoral' produces a d18O-coral record given SST, SSS, and global position.
       The function is based on the forward model published by [Thompson, 2011]:
       <Thompson, D. M., T. R. Ault, M. N. Evans, J. E. Cole, and J. Emile-Geay (2011),
       Comparison of observed and simulated tropical climate trends using a forward
       model of coral δ18O, Geophys.Res.Lett., 38, L14706, doi:10.1029/2011GL048224.>
       Returns a numpy array that is the same size and shape as the input vectors for SST, SSS.

    Input parameters:
        Latitude    [lat]       [-90, 90]
        Longitude   [lon]       [0, 360]
        SST         [SST]       SST ANOMALY Units = degrees Celsius
        SSS         [SSS]       SSS ANOMALY Units = PSU
        Please note that input files must be read in as 1-D vectors.

    Output:
        Returns a Numpy Array, ['coral'], which is saved in the main program.
        Optional Inputs: Please note that the function will use default params for
        d18O vs. SSS, a, and b unless specified and placed in the appropriate location as shown below.
        delta_18O [permil]: If the user does not have an input field for d18O, you must put -1 in its
        position in the call and use the equation for SSS below.
        **Note: [always use d18O (seawater) if available!]**
        The user may specify values for parameters a and b in the coral d18O forward model,
        which are currently set to the following default values:
        a = coral - SST regression, accounting for species-dependent variations in the SST-18O dependence.
        The code assigns this value universally at -.22, but published values vary by species/genus (i.e. Moses, 2006)
        Specify keyword input 'species=' to set slope:
        species =
                "Default":      a = -0.22
                "Porites_sp":   a = -.26178
                "Porites_lob":  a = -.19646
                "Porites_lut":  a = -.17391
                "Porites_aus":  a = -.21
                "Montast":      a = -.22124
                "Diploas":      a = -.14992
        b = d18O-SSS slope, as defined by location in: Legrande & Schmidt, (2006)
                <LeGrande, A. N., and G. A. Schmidt (2006), Global gridded data set of the oxygen isotopic          composition in seawater, Geophys. Res. Lett., 33, L12604, doi:10.1029/2006GL026011.>
                b1 = 0.31   [Red Sea]
                b2 = 0.27   [Tropical Pacific]
                b3 = 0.45   [South Pacific]
                b4 = 0.16   [Indian Ocean]
                b5 = 0.15   [Tropical Atlantic Ocean/Caribbean]
    Example of use: Call function 'pseudocoral' to calculate d18O_coral timeseries, lacking d18O sw
    but specifying own numeric parameters for a, b.
    pseudocoral(lat,lon,SST,SSS,-1,species="default",b1,b2,b3,b4,b5)
    """
    if species == 'Porites_sp':
        a = -0.26178
    elif species == 'Porites_lob':
        a = -0.19646
    elif species == 'Porites_lut':
        a = -0.17391
    elif species == 'Porites_aus':
        a = -0.21
    elif species == 'Montast':
        a = -0.22124
    elif species == 'Diploas':
        a = -0.14992
    elif species == 'default':
        a = -0.22
    else:
        a = -0.22
    V_corr = 0.97002
    slope1 = V_corr * b1
    slope2 = V_corr * b2
    slope3 = V_corr * b3
    slope4 = V_corr * b4
    slope5 = V_corr * b5
    if lon >= 32.83 and lon <= 43.5 and lat >= 12.38 and lat <= 28.5:
        b = slope1
    elif lon <= 120.0:
        b = slope4
    elif lon >= 270.0 and lat >= 10.0:
        b = slope5
    elif lat > -5.0 and lat <= 13.0:
        b = slope2
    elif lat <= -5.0:
        b = slope3
    else:
        b = slope2
    if d18O is None:
        coral = a * sst + b * sss
    else:
        coral = a * sst + 0.97002 * d18O
    return coral
