Degenerate Conic

Command Line Arguments

A standard way to access a program's command line arguments wasn't added to Fortran until 2003. Command line arguments were one of those things (like environment variables, the file system, and the operating system shell) that the Fortran standards committee pretended didn't exist or wasn't important. Of course, that was and is ridiculous. So for decades, all Fortran compilers provided a non-standard way to get the command line arguments (for example, GETARG in gfortran). Currently, the standards committee still pretends that the file system doesn't exist (for example, there's no standard way to change an application's working directory, or list the files in a directory, or even to check for the existence of directories). But I digress.

The two routines you need to know are COMMAND_ARGUMENT_COUNT and GET_COMMAND_ARGUMENT. Frequently in examples on the internet (even the gfortran example) you will see people just declare a really long string and use GET_COMMAND_ARGUMENT to get the command line argument and hope for the best. That's terrible of course, since this routine does have the ability to return the actual length of the argument, so you should use that to allocate a string with the correct size to hold it.

Consider the following example:

module command_line_module

    implicit none

    private

    type,public :: string
        !! a variable-length string type
        character(len=:),allocatable :: str
    end type

    public :: get_command_arguments

contains

subroutine get_command_arguments(nargs, args)

    !! return all the command line arguments.
    !! (each one is the correct size)

    implicit none

    integer,intent(out) :: nargs !! number of arguments
    type(string),dimension(:),allocatable,intent(out) :: args

    integer :: i, ilen

    nargs = command_argument_count()
    allocate(args(0:nargs))

    do i = 0, nargs
        call get_command_argument(i,length=ilen)
        allocate(character(len=ilen) :: args(i)%str)
        call get_command_argument(i,value=args(i)%str)
    end do

end subroutine get_command_arguments

end module command_line_module

Here we have to define our own string type (string), so we can create an array of deferred-length strings (since a useful varying-length string type is another thing Fortran pretends is not important). See StringiFor for a full-featured variable-string type, but for our purposes here the simple one is fine. The procedure in the get_command_arguments() routine is to check for the size of each argument, allocate a string large enough to hold it, and then retrieve the value. It returns an array of string variables containing the command line arguments. With a wrapper module like this, getting command line arguments is effortless, and there is no reason anymore to use the old non-standard routines.

Another One Bites the Dust: GRAM Atmosphere Model

It looks like the Fortran community has lost another venerable NASA Fortran library to C++. I recently noticed that the latest release of NASA's Global Reference Atmosphere Model (GRAM) is now C++. From the release link:

Earth-GRAM 2016 is now available as an open-source C++ computer code that can run on a variety of platforms including PCs and UNIX stations. The software provides a model that offers values for atmospheric parameters such as density, temperature, winds, and constituents for any month and at any altitude and location within the Earth's atmosphere. Earth-GRAM 2010 is available in FORTRAN. Similar computer models of Mars, Venus, Titan, and Neptune are available for diverse mission applications through NASA’s Space Environments and Effects (SEE) Program, a partnership with industry, academia, and other government agencies that seeks to develop more reliable, more effective spacecraft.

The Earth-GRAM model, originally written in Fortran 77, is produced by MSFC and has been around for a while, with many revisions (i.e., GRAM-86, GRAM-88, GRAM-90, GRAM-95, GRAM-99, GRAM 2007, and GRAM 2010). In the 2010 revision, they converted the old Fortran 77 code to Fortran 90, which I thought was a step in the right direction (although maybe 20 years too late). It seems like they could have used Fortran 2003 and provided a C interface using the standard Fortran/C Interoperability feature, which would have enabled calling it from a variety of languages (including C++).

I remain utterly unconvinced that C++ is the programming language that engineers and scientists should be using. Are CS people taking over programming jobs once done by engineers? Is this a failure of university engineering departments? Although it seems like even CS people don't want to be using C++ anymore. Every time I turn around a new programming language (Java, Rust, Go, D, C#, Swift, ...) is created by CS people fed up with the C++ dumpster fire. So I don't know who is telling engineers that C++ is the future. And the perception that Fortran is an obsolete programming language is very strong, even within the engineering/technical disciplines that it was designed for. Sometimes this is due to sheer ignorance, but sometimes it is due to ancient Fortran 77 legacy code that people have to deal with. I sympathize with them, Fortran 77 code is terrible and needs to go away, but it should be modernized, not thrown out and replaced with C++. Think of the children! (and the memory leaks!)

Parting words of wisdom

Security tips when programming in C (2017 edition):

1) Stop typing

2) Delete what you've already typed

— ryan huber (@ryanhuber) June 21, 2017

See also

• Earth-GRAM 2016 is now available!
• Earth Global Reference Atmospheric Model (Earth-GRAM 2010) -- The last Fortran 90 release. [request here]
• F.W. Leslie, C.G. Justus, "The NASA Marshall Space Flight Center Earth Global Reference Atmospheric Model—2010 Version", NASA/TM—2011–216467, 2011
• ModelWeb Catalogue and Archive -- other atmosphere models listed on GSFC ModelWeb
• D. L. Johnson, R. Blocker, C. G. Justus, Global Reference Atmosphere Model (GRAM), NASA Tech Brief, MFS-28397, Jan 01, 1993

JSON + SPICE

I have mentioned various kinds of configuration file formats used in Fortran here before. One that I haven't mentioned is the text PCK file format used by the NAIF SPICE Toolkit. This is a format that is similar in some ways to Fortran namelists, but with a better API for reading the file and querying the contents. Variables read from a PCK file are inserted into the SPICE variable pool. An example PCK file (pck00010.tpc) can be found here.

The PCK file format has some limitations. It doesn't allow inline comments (only block comments separate from the variable declarations). Integers are stored as rounded doubles, and logical variables are not supported at all (you have to use 0.0 or 1.0 for this). One particular annoyance is that the files are not cross platform (the line breaks must match the platform they are being used on). However, SPICELIB also provides routines to programmatically enter variables into the pool. This allows us to create files in other formats that can be used in a SPICE application. For example, we can use JSON-Fortran to read variables in JSON files and insert them into the SPICE pool.

First, let's declare the interfaces to the SPICE routines we need (since SPICELIB is straight up Fortran 77, they are not in a module):

interface
    subroutine pcpool ( name, n, cvals )
    implicit none
    character(len=*) :: name
    integer :: n
    character(len=*) :: cvals ( * )
    end subroutine pcpool

    subroutine pdpool ( name, n, values )
    import :: wp
    implicit none
    character(len=*) :: name
    integer :: n
    real(wp) :: values ( * )
    end subroutine pdpool

    subroutine pipool ( name, n, ivals )
    implicit none
    character(len=*) :: name
    integer :: n
    integer :: ivals ( * )
    end subroutine pipool

    subroutine setmsg ( msg )
    implicit none
    character(len=*) :: msg
    end subroutine setmsg

    subroutine sigerr ( msg )
    implicit none
    character(len=*) :: msg
    end subroutine sigerr
end interface

Then, we can use the following routine:

subroutine json_to_spice(jsonfile)

implicit none

character(len=*),intent(in) :: jsonfile
!! the JSON file to load

integer :: i,j
type(json_core) :: json
type(json_value),pointer :: p,p_var,p_element
real(wp) :: real_val
integer :: var_type,int_val,n_vars,n_children,element_var_type
logical :: found
integer,dimension(:),allocatable :: int_vec
real(wp),dimension(:),allocatable :: real_vec
character(len=:),dimension(:),allocatable :: char_vec
character(len=:),allocatable :: char_val,name
integer,dimension(:),allocatable :: ilen

nullify(p)

! allow for // style comments:
call json%initialize(comment_char='/')

! read the file:
call json%parse(file=jsonfile, p=p)

if (json%failed()) then
    ! we will use the SPICE error handler:
    call setmsg ( 'json_to_spice: Could not load file: '&
                  trim(jsonfile) )
    call sigerr ( 'SPICE(INVALIDJSONFILE)' )
else

    ! how many variables are in the file:
    call json%info(p,n_children=n_vars)
    main : do i = 1, n_vars
        call json%get_child(p, i, p_var, found)

        ! what kind of variable is it?:
        call json%info(p_var,var_type=var_type,name=name)
        if (var_type == json_array) then

            ! how many elements in this array?:
            call json%info(p_var,n_children=n_children)

            ! first make sure all the variables are the same type
            ! [must be integer, real, or character]
            do j = 1, n_children
                call json%get_child(p_var, j, p_element, found)
                call json%info(p_element,var_type=element_var_type)
                if (j==1) then
                    var_type = element_var_type
                else
                    if (var_type /= element_var_type) then
                        call setmsg ( 'json_to_spice: Invalid array ('&
                                      trim(name)//') in file: '//trim(jsonfile) )
                                      call sigerr ( 'SPICE(INVALIDJSONVAR)' )
                        exit main
                    end if
                end if
            end do

            ! now we know the var type, so get as a vector:
            select case (var_type)
            case(json_integer); call json%get(p_var,int_vec)
            case(json_double ); call json%get(p_var,real_vec)
            case(json_string ); call json%get(p_var,char_vec,ilen)
            case default
                call setmsg ( 'json_to_spice: Invalid array ('&
                              trim(name)//') in file: '//trim(jsonfile) )
                call sigerr ( 'SPICE(INVALIDJSONVAR)' )
                exit main
            end select

        else

            ! scalar:
            n_children = 1
            select case (var_type)
            case(json_integer)
                call json%get(p_var,int_val); int_vec = [int_val]
            case(json_double )
                call json%get(p_var,real_val); real_vec = [real_val]
            case(json_string )
                call json%get(p_var,char_val); char_vec = [char_val]
            case default
                call setmsg ( 'json_to_spice: Invalid variable ('&
                              trim(name)//') in file: '//trim(jsonfile) )
                call sigerr ( 'SPICE(INVALIDJSONVAR)' )
                exit main
            end select

        end if

        ! now, add them to the pool:
        select case (var_type)
        case(json_integer)
            call pipool ( name, n_children, int_vec )
        case(json_double )
            call pdpool ( name, n_children, real_vec )
        case(json_string )
            call pcpool ( name, n_children, char_vec )
        end select

    end do main

end if

call json%destroy(p)

end subroutine json_to_spice

Here we are using the SPICE routines PIPOOL, PDPOOL, and PCPOOL to insert variables into the pool. We are also using the built-in SPICE error handling routines (SETMSG and SIGERR) to manage errors. The code works for scalar and array variables. It can be used to read the following example JSON file:

{
    "str": "qwerty", // a scalar string
    "strings": ["a", "ab", "abc", "d"], // a vector of strings
    "i": 8, // a scalar integer
    "ints": [255,127,0], // a vector of integers
    "r": 123.345, // a scalar real
    "reals": [999.345, 5671.8] // a vector of reals
}

After reading it, we can verify that the variables were added to the pool by printing the pool contents using the routine WRPOOL. This produces:

\begindata

str     = 'qwerty'
i       = 0.80000000000000000D+01
strings = ( 'a',
            'ab',
            'abc',
            'd' )
ints    = ( 0.25500000000000000D+03,
            0.12700000000000000D+03,
            0.00000000000000000D+00 )
r       = 0.12334500000000000D+03
reals   = ( 0.99934500000000003D+03,
            0.56718000000000002D+04 )

\begintext

See also

• NASA's Navigation and Ancillary Information Facility
• SPICE Toolkit, N0066, released April 10, 2017

IAU SOFA Version 14

The IAU Standards of Fundamental Astronomy (SOFA) library implements standard models used in fundamental astronomy. Version 14 ( 2018-01-30) has just been released. According to the release notes, this update includes the following:

• Change in the copyright status of the iau_DAT routine. This is to provide for a user-supplied mechanism for updating the number of leap seconds. SOFA doesn't provide any such mechanism, but they are now allowing you to replace or modify this routine to provide your own without having to rename the routine (all other SOFA routines cannot be modified without renaming them).
• Implementation of two new categories of routines:
  • Three new routines for the horizon/equatorial plane coordinates.
    • iau_AE2HD — (azimuth, altitude) to (hour angle, declination)
    • iau_HD2AE — (hour angle, declination) to (azimuth, altitude)
    • iau_HD2PA — parallactic angle
  • Six new routines dealing with gnomonic (tangent plane) projections.
    • iau_TPORS — solve for tangent point, spherical
    • iau_TPORV — solve for tangent point, vector
    • iau_TPSTS — project tangent plane to celestial, spherical
    • iau_TPSTV — project tangent plane to celestial, vector
    • iau_TPXES — project celestial to tangent plane, spherical
    • iau_TPXEV — project celestial to tangent plane, vector
• The Astrometry Tools Cookbook, the test program and other supporting files have also been updated.
• Other minor documentation/typographical corrections to various files.

IAU SOFA has a Fortran 77 and a C version. While this is a great resource, I do wonder why the IAU insists on using a source format for the Fortran version that was rendered obsolete 30 years ago? There's no reason for this. This is new code that is regularly kept up to date, so why not modernize it for those of us living in the 21st century?

Consider this SOFA routine:

      SUBROUTINE iau_TR ( R, RT )

      IMPLICIT NONE

      DOUBLE PRECISION R(3,3), RT(3,3)

      DOUBLE PRECISION WM(3,3)
      INTEGER I, J

      DO 2 I=1,3
        DO 1 J=1,3
          WM(I,J) = R(J,I)
1       CONTINUE
2     CONTINUE
      CALL iau_CR ( WM, RT )

      END SUBROUTINE iau_TR

All this is doing is transposing a 3x3 matrix. Well, since the late 1990s, to do that in Fortran all you need is:

rt = transpose(r)

The original routine is just awful: all caps, line numbers, CONTINUE statements, and 100% unnecessary. Unfortunately, it's code like this that gives Fortran a bad name (reinforcing the perception that Fortran is an arcane and obsolete programming language). Anyone who is compiling this code is probably using a compiler that is at least Fortran 95 compatible. Who on earth is using straight up Fortran 77 compilers at this point? So why restrict such a useful library to a programming style that is so totally out of date? I think it is even probably affecting the efficiency of the code, since the example above required using a temporary array and a function call to do what can now be done with a (presumably much more efficient) intrinsic routine. It wouldn't even be that hard to convert this code to modern Fortran style. It's just low-level math routines, it doesn't have to be object-oriented or anything fancy like that. A while back, I wrote a little Python script to merge all the SOFA files into a single Fortran module (another innovation from the 1990s) so they can take advantage of the automatic interface checking that modules provide. There are also any number of tools out that could be used to automate the fixed to free-form conversion.

See also

• IAU Standards of Fundamental Astronomy
• SOFA Library Issue 2018-01-30 for Fortran 77

Fortran 2018

The upcoming Fortran standard formerly known as Fortran 2015 has a new name: Fortran 2018. It was decided to change it in order to match the expected year of publication. This makes sense. The previous standard (Fortran 2008) was published in 2010.

Waiting for an updated Fortran standard is an exercise in Zen-like patience. Almost a decade after Fortran 2008, we'll get a fairly minor update to the core language. And it will be years after that before it's fully supported by any compiler that most users will have available (gfortran still doesn't have a bug-free implementation of all of Fortran 2008 or even 2003). Fortran was essentially reborn in the Fortran 2003 standard, which was an amazing update that brought Fortran into the modern world. It's a terrific programming language for scientific and technical computing. However, the limitations are all too clear and have been for a long time:

• We need better facilities for generic programming. It's impossible to do some things without having to duplicate code or use "tricks" like preprocessing or include files (see JSON-Fortran for examples).
• We need some kind of exception handling. Fortran does have a floating-point exception handling feature, but honestly, it's somewhat half-baked.
• We need a better implementation of strings. Allocatable strings (introduced in Fortran 2003) are great, but not enough, since they can't be used in all instances where strings are needed.
• We need the language to be generally less verbose. I'm tired to having to type multiple nested SELECT TYPE statements to do something that is a one liner in Python (sure I know Fortran will never be as succinct as Python, but some of the verbosity required for object-oriented Fortran is just perverse).
• We need any number of new features to make it easier to extend the language with third-party libraries (so we don't have to wait two decades for a feature we want).
• We also need the language development process to embrace a more open collaborative model using modern tools (Usenet is not the future). I guess the recent survey was unprecedented, but it's not enough.

Fortran is a programming language that needs a better PR department. Legacy Fortran codes are being rewritten in C++, Python, or even Julia. NASA frequently throws massive Fortran 77 libraries with decades of heritage (e.g., DPTRAJ/ODP, GTDS, SPICELIB) into the trash in order to rewrite it all from the ground up in C++, without ever considering Fortran 2003+ (or maybe not realizing it exists?). The information about modern Fortran on the internet is spotty at best, outdated, or downright wrong (what is the deal with REAL*8?). In popular consciousness Fortran is mostly a punchline (usually something to do with punchcards and your granddad). A language like Python (which was never designed for technical computing) is now seen by many as a superior solution for technical computing. Matlab refers to itself as "the only top programming language dedicated to mathematical and technical computing"! The Julia website lists somewhat misleading benchmarks than implies that C, Julia, and even Lua are faster than Fortran.

Now, get off my lawn, you kids!

See also

• N2144 Results of Straw Poll on Changes to the Fortran Draft Standard [JTC1/SC22/WG5]
• S. Lionel, "Eighteen is the new Fifteen", November 20, 2017 [Intel]
• Fortran standards committee Fortran 202x Feature Suggestions Survey [Surveymonkey]
• When will see Fortran Compilers conformant to Fortran 2015? [comp.lang.fortran]

Intel Fortran Compiler 18.0

Intel has just released version 18 of the Intel Fortran Compiler (part of Intel Parallel Studio XE 2018). At long last, this release includes full support for the Fortran 2008 standard. The updates since the previous compiler release include:

• COMPILER_OPTIONS and COMPILER_VERSION in ISO_FORTRAN_ENV
• Complex arguments to trigonometric and hyperbolic intrinsic functions
• FINDLOC intrinsic function
• Optional argument BACK in MAXLOC and MINLOC intrinsic functions
• Multiple type-bound procedures in a PROCEDURE list
• Passing a non-pointer data item to a pointer dummy argument
• Polymorphic assignment with allocatable Left Hand Side (LHS)
• Allocatable components of recursive type and forward reference
• Data statement restrictions removed

In addition, the new release also includes support for all the features from "Technical Specification 29113 Further Interoperability with C", planned for inclusion in Fortran 2015. These include:

• Assumed type (TYPE(*))
• Assumed rank (DIMENSION(..))
• Relaxed restrictions on interoperable dummy arguments
• ISO_Fortran_binding.h C include file for use by C code manipulating "C descriptors" used by Fortran

Hopefully, it won't take so long to get the compiler up to full Fortran 2015 compliance (see a previous post for a list of new Fortran 2015 features).

See also

• Igor V. (Intel), Intel® Visual Fortran Compiler 18.0 for Windows Release Notes for Intel Parallel Studio XE 2018, Updated April 5, 2017.

Time Conversions with SPICE

JPL's SPICE Toolkit (SPICELIB) is the premier software library for computations related to solar system geometry. It is freely distributed, and is also one of the best-documented libraries I have ever come across. SPICELIB also includes a comprehensive set of routines for date and time conversions. An example is shown here:

program spice_test

use iso_fortran_env, only: wp => real64

implicit none

interface
    ! SPICELIB routines
    subroutine timout ( et, pictur, output )
        import :: wp
        implicit none
        real(wp),intent(in) :: et
        character(len=*),intent(in) :: pictur
        character(len=*),intent(out) :: output
    end subroutine timout
    subroutine str2et ( string, et )
        import :: wp
        implicit none
        character(len=*),intent(in) :: string
        real(wp),intent(out) :: et
    end subroutine str2et
    subroutine furnsh ( file )
        implicit none
        character(len=*),intent(in) :: file
    end subroutine furnsh
end interface

character(len=*),parameter :: time_in = &
    '2017 Aug 12 00:00:00 TDB'
character(len=*),parameter :: pictur = &
    'Mon DD,YYYY HR:MN:SC.#### UTC ::UTC'
real(wp) :: et
character(len=100) :: time_out

! load the leap second kernel:
call furnsh('naif0012.tls')

! example conversion:
call str2et(time_in, et)
call timout(et, pictur, time_out)

write(*,*) 'time_in: ', time_in
write(*,*) 'et: ', et
write(*,*) 'time_out: ', time_out

end program spice_test

A few things to note:

• Here we are using the SPICE routines str2et and timout to convert a string from a TDB calendar date to ephemeris time and then to a UTC calendar date. These routines are very flexible and can convert a wide range of date formats. Other routines are available to do other transformations.
• The base time system of SPICE is Barycentric Dynamical Time (TDB). "Ephemeris time" is a count of TDB seconds since the J2000 epoch (Jan 1, 2000 12:00:00).
• We have to load the latest leap second kernel (naif0012.tls in this case), which is necessary to define UTC.
• The SPICE routines are not in a module (the code is Fortran 77), and so have no explicit interfaces. Thus it is good practice to specify them as I do here.

The output of this example is:

time_in:  2017 Aug 12 00:00:00 TDB
et:       555768000.00000000
time_out: Aug 11,2017 23:58:50.8169 UTC

See also

• Time Routines in SPICE [NAIF]
• NAIF SPICE Toolkit Hypertext Documentation [NAIF]

Another One Bites the Dust

JPL recently released an update to their awesome SPICE Toolkit (it is now at version N66). The major new feature in this release is the Digital Shape Kernel (DSK) capability to define the shapes of bodies (such as asteroids) via tessellated plate models.

Unfortunately for Fortran users, they also announced that they have decided to reimplement the entire library in C++. SPICELIB is currently written in Fortran 77, which they f2c to provide a C version (which is also callable from IDL, Matlab, and Python, among others). Their reason for this "upgrade" is to provide thread safety and object oriented features. Of course, modern Fortran can be thread safe and object oriented, and upgrading the code to modern standards could be done in a fraction of the time it will take to rewrite everything from scratch in C++. SPICELIB is extremely well-written Fortran 77 code, and is not infested with COMMON blocks, EQUIVALENCE statements, etc. I actually don't think it would take much effort to modernize it. In addition, Fortran/C interoperability could be employed to easily provide an interface that is callable from C without source transformation.

However, I guess it isn't meant to be, and the science/engineering community will lose another Fortran code to C++ like many times before, in spite of C++ being a terrible language for scientists and engineers.

See also

• The SPICE Toolkit [JPL]
• SPICE Newsletter, April 10, 2017

The New Features of Fortran 2015

The glacially slow pace of Fortran language development continues! The next standard, Fortran 2015, mainly consists of updates for Fortran/C interoperability and new coarray features such as teams. In addition, there are a bunch of minor changes and discrepancy fixes. A few of the new features are:

• The venerable implicit none statement has been updated to allow for some additional use related to external procedures.
• The stop code in error stop can now be any integer or character expression.
• An out_of_range intrinsic was added to allow for testing whether a real or integer value can be safely converted to a different real or integer type and kind.
• You can now declare the kind of the loop variable inside an implied do loop. For example: iarray = [(2*i, integer :: i=1,n)].
• All procedures are now recursive by default. This is an interesting change. Ever since recursion was added to the language in Fortran 90, a procedure has had to be explicitly declared as recursive. Now, you have to use non_recursive if you don't want to allow a procedure to be used recursively.
• Some new syntax to indicate the locality status of variables within a do concurrent loop.
• There are a lot of new IEEE intrinsic routines for some reason.

Fortran 2015 is expected to be published in 2018.

See also

• J. Reid, The new features of Fortran 2015, ISO/IEC JTC1/SC22/WG5 N2127, June 23, 2017 [Note: this link is weird, it downloads to my Mac as a '.pdf.flv' file. Just remove the '.flv' extension, it's just a pdf]
• S. Lionel, Doctor Fortran in "Are We There Yet?", August 6, 2017
• JTC1/SC22/WG5 (New website for the home of Fortran standards)

Fortran in the Cloud

You can now try OpenCoarrays and Gfortran in the cloud, courtesy of Zaak Beekman and the Sourcery Institute. Just navigate to http://bit.ly/TryCoarrays and then click "Launch". This awesome project is enabled by various other awesome tools like Binder, Jupyter, and GitHub. Truly, we are living in the future.

Coarrays are the parallel processing component built into the Fortran language (standardized in Fortran 2008). It uses the Partitioned Global Address Space (PGAS) and Single-Program-Multiple-Data (SPMD) programming models. OpenCoarrays is an open source library to enable coarray usage in Gfortran.

See also

• Coarray Fortran Jupyter notebook kernel [GitHub]
• J. Reid, Coarrays in the next Fortran Standard, June 15, 2009