trara1 Subroutine

    subroutine trara1(me, descr, map, fl, bb0, e, f, n)

        class(trm_type), intent(inout) :: me
        integer, intent(in) :: n !! number of energies
        integer, intent(in) :: descr(8) !! header of specified trapped radition model
        real(wp), intent(in) :: e(n) !! array of energies in mev
        real(wp), intent(in) :: fl !! l-value
        real(wp), intent(in) :: bb0 !! =b/b0  magnetic field strength normalized
                                 !! to field strength at magnetic equator
        integer, intent(in) :: map(*) !! map of trapped radition model
                                !! (descr and map are explained at the begin
                                !! of the main program model)
        real(wp), intent(out) :: f(n) !! decadic logarithm of integral fluxes in
                                !! particles/(cm*cm*sec)

        real(wp) :: e0, e1, e2, escale, f0, fscale, xnl
        real(wp) :: bb0_ !! local copy of `bb0`. in the original code
                    !! this was modified by this routine.
                    !! added this so `bb0` could be `intent(in)`
        integer :: i0, i1, i2, i3, ie, l3, nb, nl
        logical :: s0, s1, s2

        e0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings
        f0 = 0.0_wp ! to avoid -Wmaybe-uninitialized warnings
        i0 = 0 ! to avoid -Wmaybe-uninitialized warnings
        s0 = .false. ! to avoid -Wmaybe-uninitialized warnings  -- but not sure what default value here should be !  -JW

        bb0_ = bb0
        me%fistep = descr(7) / descr(2)
        escale = descr(4)
        fscale = descr(7)
        xnl = min(15.6_wp, abs(fl))
        nl = int(xnl * descr(5))
        if (bb0_ < 1.0_wp) bb0_ = 1.0_wp
        nb = int((bb0_ - 1.0_wp) * descr(6))

        ! i2 is the number of elements in the flux map for the first energy.
        ! i3 is the index of the last element of the second energy map.
        ! l3 is the length of the map for the third energy.
        ! e1 is the energy of the first energy map (unscaled)
        ! e2 is the energy of the second energy map (unscaled)
        i1 = 0
        i2 = map(1)
        i3 = i2 + map(i2 + 1)
        l3 = map(i3 + 1)
        e1 = map(i1 + 2) / escale
        e2 = map(i2 + 2) / escale

        ! s0, s1, s2 are logical variables which indicate whether the flux for
        ! a particular e, b, l point has already been found in a previous call
        ! to function trara2. if not, s.. =.true.
        s1 = .true.
        s2 = .true.

        ! energy loop

        do ie = 1, n

            ! for each energy e(i) find the successive energies e0,e1,e2 in
            ! model map, which obey  e0 < e1 < e(i) < e2 .

            do while ((e(ie) > e2) .and. (l3 /= 0))
                i0 = i1
                i1 = i2
                i2 = i3
                i3 = i3 + l3
                l3 = map(i3 + 1)
                e0 = e1
                e1 = e2
                e2 = map(i2 + 2) / escale
                s0 = s1
                s1 = s2
                s2 = .true.
                f0 = me%f1
                me%f1 = me%f2
            end do

            ! call trara2 to interpolate the flux-maps for e1,e2 in l-b/b0-
            ! space to find fluxes f1,f2 [if they have not already been
            ! calculated for a previous e(i)].

            if (s1) me%f1 = me%trara2(map(i1 + 3), nl, nb) / fscale
            if (s2) me%f2 = me%trara2(map(i2 + 3), nl, nb) / fscale
            s1 = .false.
            s2 = .false.

            ! finally, interpolate in energy.

            f(ie) = me%f1 + (me%f2 - me%f1) * (e(ie) - e1) / (e2 - e1)
            if (me%f2 <= 0.0_wp) then
                if (i1 /= 0) then
                    ! --------- special interpolation ---------------------------------
                    ! if the flux for the second energy cannot be found (i.e. f2=0.0),
                    ! and the zeroth energy map has been defined (i.e. i1 not equal 0),
                    ! then interpolate using the flux maps for the zeroth and first
                    ! energy and choose the minimum of this interpolations and the
                    ! interpolation that was done with f2=0.
                    if (s0) f0 = me%trara2(map(i0 + 3), nl, nb) / fscale
                    s0 = .false.
                    f(ie) = min(f(ie), f0 + (me%f1 - f0) * (e(ie) - e0) / (e1 - e0))
                end if
            end if

            ! the logarithmic flux is always kept greater or equal zero.

            f(ie) = max(f(ie), 0.0_wp)
        end do
    end subroutine trara1