feldcof Subroutine

    subroutine feldcof(me, year, dimo)

        class(shellig_type), intent(inout) :: me
        real(wp), intent(in) :: year !! decimal year for which geomagnetic field is to
                               !! be calculated (e.g.:1995.5 for day 185 of 1995)
        real(wp), intent(out) :: dimo !! geomagnetic dipol moment in gauss (normalized
                                !! to earth's radius) at the time (year)

        real(wp) :: dte1, dte2, erad, gha(144), sqrt2
        integer :: i, ier, j, l, m, n, iyea
        character(len=:), allocatable :: fil2
        real(wp) :: x, f0, f !! these were double precision in original
                          !! code while everything else was single precision

        ! changed to conform with IGRF 45-95, also FILMOD, DTEMOD arrays +1
        character(len=filename_len), dimension(17), parameter :: filmod = [ &
                                                               'dgrf1945.dat ', 'dgrf1950.dat ', 'dgrf1955.dat ', 'dgrf1960.dat ', &
                                                               'dgrf1965.dat ', 'dgrf1970.dat ', 'dgrf1975.dat ', 'dgrf1980.dat ', &
                                                               'dgrf1985.dat ', 'dgrf1990.dat ', 'dgrf1995.dat ', 'dgrf2000.dat ', &
                                                               'dgrf2005.dat ', 'dgrf2010.dat ', 'dgrf2015.dat ', 'igrf2020.dat ', &
                                                                 'igrf2020s.dat']
        real(wp), dimension(17), parameter :: dtemod = [1945.0_wp, 1950.0_wp, 1955.0_wp, &
                                                        1960.0_wp, 1965.0_wp, 1970.0_wp, &
                                                        1975.0_wp, 1980.0_wp, 1985.0_wp, &
                                                        1990.0_wp, 1995.0_wp, 2000.0_wp, &
                                                        2005.0_wp, 2010.0_wp, 2015.0_wp, &
                                                        2020.0_wp, 2025.0_wp]
        integer, parameter :: numye = size(dtemod) - 1 ! number of 5-year priods represented by IGRF
        integer, parameter :: is = 0 !! * is=0 for schmidt normalization
                               !! * is=1 gauss normalization

        logical :: read_file

        !-- determine igrf-years for input-year
        me%time = year
        iyea = int(year / 5.0_wp) * 5
        read_file = iyea /= me%iyea ! if we have to read the file
        me%iyea = iyea
        l = (me%iyea - 1945) / 5 + 1
        if (l < 1) l = 1
        if (l > numye) l = numye
        dte1 = dtemod(l)
        me%name = me%get_data_file_dir()//trim(filmod(l))
        dte2 = dtemod(l + 1)
        fil2 = me%get_data_file_dir()//trim(filmod(l + 1))
        if (read_file) then
            ! get igrf coefficients for the boundary years
            ! [if they have not ready been loaded]
            call getshc(me%name, me%nmax1, erad, me%g, ier)
            if (ier /= 0) error stop 'error reading file: '//trim(me%name)
            me%g_cache = me%g ! because it is modified below, we have to cache the original values from the file
            call getshc(fil2, me%nmax2, erad, me%gh2, ier)
            if (ier /= 0) error stop 'error reading file: '//trim(fil2)
        else
            me%g = me%g_cache
        end if
        !-- determine igrf coefficients for year
        if (l <= numye - 1) then
            call intershc(year, dte1, me%nmax1, me%g, dte2, me%nmax2, me%gh2, me%nmax, gha)
        else
            call extrashc(year, dte1, me%nmax1, me%g, me%nmax2, me%gh2, me%nmax, gha)
        end if
        !-- determine magnetic dipol moment and coeffiecients g
        f0 = 0.0_wp
        do j = 1, 3
            f = gha(j) * 1.0e-5_wp
            f0 = f0 + f * f
        end do
        dimo = sqrt(f0)

        me%g(1) = 0.0_wp
        i = 2
        f0 = 1.0e-5_wp
        if (is == 0) f0 = -f0
        sqrt2 = sqrt(2.0_wp)

        do n = 1, me%nmax
            x = n
            f0 = f0 * x * x / (4.0_wp * x - 2.0_wp)
            if (is == 0) f0 = f0 * (2.0_wp * x - 1.0_wp) / x
            f = f0 * 0.5_wp
            if (is == 0) f = f * sqrt2
            me%g(i) = gha(i - 1) * f0
            i = i + 1
            do m = 1, n
                f = f * (x + m) / (x - m + 1.0_wp)
                if (is == 0) f = f * sqrt((x - m + 1.0_wp) / (x + m))
                me%g(i) = gha(i - 1) * f
                me%g(i + 1) = gha(i) * f
                i = i + 2
            end do
        end do

    end subroutine feldcof