docs: 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