feldg Subroutine

private subroutine feldg(me, glat, glon, alt, bnorth, beast, bdown, Babs)

Calculates earth magnetic field from spherical harmonics model

Reference

ref: g. kluge, european space operations centre, internal note 61, 1970.

History

  • changes (d. bilitza, nov 87):
  • field coefficients in binary data files instead of block data
  • calculates dipol moment

Note

In the original code, [[feldc] and feldi were ENTRY points to this routine

Type Bound

shellig_type

Arguments

Type IntentOptional Attributes Name
class(shellig_type), intent(inout) :: me
real(kind=wp), intent(in) :: glat

geodetic latitude in degrees (north)
real(kind=wp), intent(in) :: glon

geodetic longitude in degrees (east)
real(kind=wp), intent(in) :: alt

altitude in km above sea level
real(kind=wp), intent(out) :: bnorth

components of the field with respect to the local geodetic coordinate system, with axis pointing in the tangential plane to the north, east and downward.
real(kind=wp), intent(out) :: beast

components of the field with respect to the local geodetic coordinate system, with axis pointing in the tangential plane to the north, east and downward.
real(kind=wp), intent(out) :: bdown

components of the field with respect to the local geodetic coordinate system, with axis pointing in the tangential plane to the north, east and downward.
real(kind=wp), intent(out) :: Babs

magnetic field strength in gauss

Called by

proc~~feldg~~CalledByGraph proc~feldg shellig_module::shellig_type%feldg proc~igrf shellig_module::shellig_type%igrf proc~igrf->proc~feldg proc~get_flux_g_ radbelt_module::radbelt_type%get_flux_g_ proc~get_flux_g_->proc~igrf none~get_flux radbelt_module::radbelt_type%get_flux none~get_flux->proc~get_flux_g_ proc~get_flux_c radbelt_module::get_flux_c proc~get_flux_c->none~get_flux proc~get_flux_g radbelt_module::get_flux_g proc~get_flux_g->none~get_flux proc~get_flux_g_c radbelt_c_module::get_flux_g_c proc~get_flux_g_c->none~get_flux interface~get_flux radbelt_module::get_flux interface~get_flux->proc~get_flux_c interface~get_flux->proc~get_flux_g

Source Code

    subroutine feldg(me, glat, glon, alt, bnorth, beast, bdown, babs)

        class(shellig_type), intent(inout) :: me
        real(wp), intent(in) :: glat !! geodetic latitude in degrees (north)
        real(wp), intent(in) :: glon !! geodetic longitude in degrees (east)
        real(wp), intent(in) :: alt !! altitude in km above sea level
        real(wp), intent(out) :: bnorth, beast, bdown !! components of the field with respect
                                                !! to the local geodetic coordinate system, with axis
                                                !! pointing in the tangential plane to the north, east
                                                !! and downward.
        real(wp), intent(out) :: Babs !! magnetic field strength in gauss

        real(wp) :: brho, bxxx, byyy, bzzz, cp, ct, d, f, rho, &
                    rlat, rlon, rq, s, sp, st, t, &
                    x, xxx, y, yyy, z, zzz
        integer :: i, ih, ihmax, il, imax, k, last, m

        ! same calculation as geo_to_cart, but not used here
        ! because the intermediate variables are also used below.
        rlat = glat * umr
        ct = sin(rlat)
        st = cos(rlat)
        d = sqrt(aquad - (aquad - bquad) * ct * ct)
        rlon = glon * umr
        cp = cos(rlon)
        sp = sin(rlon)
        zzz = (alt + bquad / d) * ct / era
        rho = (alt + aquad / d) * st / era
        xxx = rho * cp
        yyy = rho * sp

        rq = 1.0_wp / (xxx * xxx + yyy * yyy + zzz * zzz)
        me%xi = [xxx, yyy, zzz] * rq

        ihmax = me%nmax * me%nmax + 1
        last = ihmax + me%nmax + me%nmax
        imax = me%nmax + me%nmax - 1
        do i = ihmax, last
            me%h(i) = me%g(i)
        end do
        do k = 1, 3, 2
            i = imax
            ih = ihmax
            do
                il = ih - i
                f = 2.0_wp / real(i - k + 2, wp)
                x = me%xi(1) * f
                y = me%xi(2) * f
                z = me%xi(3) * (f + f)
                i = i - 2
                if ((i - 1) >= 0) then
                    if ((i - 1) > 0) then
                        do m = 3, i, 2
                            me%h(il + m + 1) = me%g(il + m + 1) + z * me%h(ih + m + 1) + x * (me%h(ih + m + 3) - &
                                                                       me%h(ih + m - 1)) - y * (me%h(ih + m + 2) + me%h(ih + m - 2))
                            me%h(il + m) = me%g(il + m) + z * me%h(ih + m) + x * (me%h(ih + m + 2) - &
                                                                       me%h(ih + m - 2)) + y * (me%h(ih + m + 3) + me%h(ih + m - 1))
                        end do
                    end if
                    me%h(il + 2) = me%g(il + 2) + z * me%h(ih + 2) + x * me%h(ih + 4) - y * (me%h(ih + 3) + me%h(ih))
                    me%h(il + 1) = me%g(il + 1) + z * me%h(ih + 1) + y * me%h(ih + 4) + x * (me%h(ih + 3) - me%h(ih))
                end if
                me%h(il) = me%g(il) + z * me%h(ih) + 2.0_wp * (x * me%h(ih + 1) + y * me%h(ih + 2))
                ih = il
                if (i < k) exit
            end do
        end do

        s = 0.5_wp * me%h(1) + 2.0_wp * (me%h(2) * me%xi(3) + me%h(3) * me%xi(1) + me%h(4) * me%xi(2))
        t = (rq + rq) * sqrt(rq)
        bxxx = t * (me%h(3) - s * xxx)
        byyy = t * (me%h(4) - s * yyy)
        bzzz = t * (me%h(2) - s * zzz)

        babs = sqrt(bxxx * bxxx + byyy * byyy + bzzz * bzzz)
        beast = byyy * cp - bxxx * sp
        brho = byyy * sp + bxxx * cp
        bnorth = bzzz * st - brho * ct
        bdown = -bzzz * ct - brho * st

    end subroutine feldg