def doigrf(lon, lat, alt, date, **kwargs):
    """
    Calculates the interpolated (<2015) or extrapolated (>2015) main field and
    secular variation coefficients and passes them to the Malin and Barraclough
    routine (function pmag.magsyn) to calculate the field from the coefficients.

    Parameters:
    -----------
    lon  : east longitude in degrees (0 to 360 or -180 to 180)
    lat   : latitude in degrees (-90 to 90)
    alt   : height above mean sea level in km (itype = 1 assumed)
    date  : Required date in years and decimals of a year (A.D.)

    Optional Parameters:
    -----------
    coeffs : if True, then return the gh coefficients
    mod  : model to use ('arch3k','cals3k','pfm9k','hfm10k','cals10k.2','cals10k.1b','shadif14k')
        arch3k (Korte et al., 2009)
        cals3k (Korte and Constable, 2011)
        cals10k.1b (Korte et al., 2011)
        pfm9k  (Nilsson et al., 2014)
        hfm.OL1.A1 (Constable et al., 2016)
        cals10k.2 (Constable et al., 2016)
        shadif14k (Pavon-Carrasco et al. (2014)
          NB : the first four of these models, are constrained to agree
               with gufm1 (Jackson et al., 2000) for the past four centuries
    Return
    -----------
    x : north component of the magnetic field in nT
    y : east component of the magnetic field in nT
    z : downward component of the magnetic field in nT
    f : total magnetic field in nT

    By default, igrf12 coefficients are used between 1900 and 2020
    from http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html.


    To check the results you can run the interactive program at the NGDC
    www.ngdc.noaa.gov/geomag-web
    """
    from . import coefficients as cf
    gh, sv = [], []
    colat = 90. - lat
#! convert to colatitude for MB routine
    if lon < 0:
        lon = lon + 360.
# ensure all positive east longitudes
    itype = 1
    models, igrf12coeffs = cf.get_igrf12()
    if 'mod' in list(kwargs.keys()):
        if kwargs['mod'] == 'arch3k':
            psvmodels, psvcoeffs = cf.get_arch3k()  # use ARCH3k coefficients
        elif kwargs['mod'] == 'cals3k':
            # use CALS3K_4b coefficients between -1000,1940
            psvmodels, psvcoeffs = cf.get_cals3k()
        elif kwargs['mod'] == 'pfm9k':
            # use PFM9k (Nilsson et al., 2014), coefficients from -7000 to 1900
            psvmodels, psvcoeffs = cf.get_pfm9k()
        elif kwargs['mod'] == 'hfm10k':
            # use HFM.OL1.A1 (Constable et al., 2016), coefficients from -8000
            # to 1900
            psvmodels, psvcoeffs = cf.get_hfm10k()
        elif kwargs['mod'] == 'cals10k.2':
            # use CALS10k.2 (Constable et al., 2016), coefficients from -8000
            # to 1900
            psvmodels, psvcoeffs = cf.get_cals10k_2()
        elif kwargs['mod'] == 'shadif14k':
            # use CALS10k.2 (Constable et al., 2016), coefficients from -8000
            # to 1900
            psvmodels, psvcoeffs = cf.get_shadif14k()
        else:
            # Korte and Constable, 2011;  use prior to -1000, back to -8000
            psvmodels, psvcoeffs = cf.get_cals10k()
# use geodetic coordinates
    if 'models' in kwargs:
        if 'mod' in list(kwargs.keys()):
            return psvmodels, psvcoeffs
        else:
            return models, igrf12coeffs
    if date < -12000:
        print('too old')
        return
    if 'mod' in list(kwargs.keys()) and kwargs['mod'] == 'shadif14k':
        if date < -10000:
            incr = 100
        else:
            incr = 50
        model = date - date % incr
        gh = psvcoeffs[psvmodels.index(int(model))]
        sv = old_div(
            (psvcoeffs[psvmodels.index(int(model + incr))] - gh), float(incr))
        x, y, z, f = magsyn(gh, sv, model, date, itype, alt, colat, lon)
    elif date < -1000:
        incr = 10
        model = date - date % incr
        gh = psvcoeffs[psvmodels.index(int(model))]
        sv = old_div(
            (psvcoeffs[psvmodels.index(int(model + incr))] - gh), float(incr))
        x, y, z, f = magsyn(gh, sv, model, date, itype, alt, colat, lon)
    elif date < 1900:
        if kwargs['mod'] == 'cals10k':
            incr = 50
        else:
            incr = 10
        model = date - date % incr
        gh = psvcoeffs[psvmodels.index(model)]
        if model + incr < 1900:
            sv = old_div(
                (psvcoeffs[psvmodels.index(model + incr)] - gh), float(incr))
        else:
            field2 = igrf12coeffs[models.index(1940)][0:120]
            sv = old_div((field2 - gh), float(1940 - model))
        x, y, z, f = magsyn(gh, sv, model, date, itype, alt, colat, lon)
    else:
        model = date - date % 5
        if date < 2015:
            gh = igrf12coeffs[models.index(model)]
            sv = old_div((igrf12coeffs[models.index(model + 5)] - gh), 5.)
            x, y, z, f = magsyn(gh, sv, model, date, itype, alt, colat, lon)
        else:
            gh = igrf12coeffs[models.index(2015)]
            sv = igrf12coeffs[models.index(2015.20)]
            x, y, z, f = magsyn(gh, sv, model, date, itype, alt, colat, lon)
    if 'coeffs' in list(kwargs.keys()):
        return gh
    else:
        return x, y, z, f