how to update the model for earth magnetic field

An international body with acrony IAGA produces formulae for the earth's magnetic field (https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) and oce tracks these changes for use in its magdec function.

The formulae are provided in a standalone Fortran program (https://www.ngdc.noaa.gov/IAGA/vmod/igrf13.f) and this must manipulated to get it into a form that can be used by oce. This wiki page outlines the procedure for that manipulation.

In March 2020, I updated from IAGA version 12 to version 13. This wiki page is an approximate record of my work.

The steps are as follows. But, before proceeding, the reader is cautioned that these instructions are for the developers, or for others who are comfortable programming in Fortran and are aware of its evolution over the past 60 years, for the code has remnants that go back to 1960s style of coding.

1. In a unix shell:

   1. git checkout develop # switch to develop branch
   2. git pull # get the latest updates
   3. git checkout magneticField # switch branch
   4. `git pull # get updates
   5. git merge --no-ff develop # now we have the "develop" updates and are ready to code
   6. cd src # got to right directory

2. In R:

   1. download.file("https://www.ngdc.noaa.gov/IAGA/vmod/igrf13.f", "igrf13.f")

3. In a unix shell:

   1. git add igrf13.f
   2. git commit -m "download igrf13.f"

4. Change igrf13.f in stages. To make it easier to follow changes, it makes sense to store explanations in git comments, and appropriate git-commit comments are given.

   1. delete the main program, the DMDDEC subroutine, and the DDECDM subroutine. Commit note: "delete MAIN, DMDDEC and DDECDM".
   2. notice that in line 58 and many lines thereafter, there is a comment at the end. On line 58, it is 1900 (space, then the number). This is following an old Fortran convention that any characters appearing in columns 73 to 80 is ignored. When R tries to compile the code, a warning will be issued for every one of these lines. Since the CRAN team does not approve of warnings, these lines should all be edited. Do not remove material in column 72, because that will entirely break the fortran code. Be careful, since you are editing about 500 lines, and doing this without a decent programming editor is a recipe for errors. Commit note: "remove column 73-80 comments".

5. At this stage, the code should compile without errors, but it is likely to have lots of warnings, e.g. the following. You ought to check that you get only these warnings, because if more crop up, that will mean some more thought will be required.

igrf13.f:670:72:

  670 |    10 m     = m + 1
      |                                                                        1
Warning: Fortran 2018 deleted feature: DO termination statement which is not END DO or CONTINUE with label 10 at (1)
igrf13.f:546:14:

  546 |       ll    = t
      |              1
Warning: Possible change of value in conversion from REAL(8) to INTEGER(4) at (1) [-Wconversion]
igrf13.f:560:15:

  560 |        ll    = 0.2*(date - 1995.0)
      |               1
Warning: Possible change of value in conversion from REAL(8) to INTEGER(4) at (1) [-Wconversion]
igrf13.f:646:0:

  646 |        three = (fn + gn)/one
      | 
Warning: 'gn' may be used uninitialized in this function [-Wmaybe-uninitialized]
igrf13.f:631:0:

  631 |        fn    = n
      | 
Warning: 'fn' may be used uninitialized in this function [-Wmaybe-uninitialized]

6. Since we want oce to build without warnings, these must be addressed. It's a bit pointless to explain the details of how to fix these (it's a Fortran lesson) so I'll just give the diff, which will tell what I did:

diff --git a/src/igrf13.f b/src/igrf13.f
index fd31699c..d227a008 100644
--- a/src/igrf13.f
+++ b/src/igrf13.f
@@ -536,6 +536,15 @@ c
       x     = 0.0
       y     = 0.0
       z     = 0.0
+c *** oce: give several fn, three and gn initial values, to prevent
+c *** oce: warnings issued during compilation by R. (It seems likely
+c *** oce: that the original code was either using knowledge that
+c *** oce: things would always be assigned, or was using the
+c *** oce: zero-initialization property of fortran.)
+      fn = 0.0
+      three = 0.0
+      gn = 0.0
+
       if (date.lt.1900.0.or.date.gt.2030.0) go to 11
       if (date.gt.2025.0) write (6,960) date
   960 format (/' This version of the IGRF is intended for use up',
@@ -543,7 +552,7 @@ c
      2        ' but may be of reduced accuracy'/)
       if (date.ge.2020.0) go to 1
       t     = 0.2*(date - 1900.0)                                             
-      ll    = t
+      ll    = int(t)
       one   = ll
       t     = t - one
 c
@@ -557,7 +566,7 @@ c
       else
        nmx   = 13
        nc    = nmx*(nmx+2)
-       ll    = 0.2*(date - 1995.0)
+       ll    = int(0.2*(date - 1995.0))
 c
 c     19 is the number of SH models that extend to degree 10
 c
@@ -623,7 +632,7 @@ c
       p(3)  = st
       q(1)  = 0.0
       q(3)  =  ct
-      do 10 k=2,kmx                                                       
+      do k=2,kmx                                                       
        if (n.ge.m) go to 4
        m     = 0
        n     = n + 1
@@ -668,6 +677,7 @@ c
        z     = z - (fn + 1.0)*one*p(k)
        l     = l + 1
    10 m     = m + 1
+       end do
 c
 c     conversion to coordinate system specified by itype
 c

We must remove all code that writes messages to the terminal output. There is a way to do that in C, so I suspect there is a way also in Fortran, but I cannot be bothered to look that up, and instead I will just comment all those lines out. (I am commenting them out instead of deleting them because if bugs arise, we might want to reinvigorate these lines.) Commit note: " remove fortran warnings, comment-out write()" 7. Now, it should build without warnings. If not, you're on your own, maties, but please let me know what you've tried 8. If things are compiling without warnings, edit src/magdec.f to handle the new function. Doing that should be obvious to someone familiar with fortran. Commit note: "let fortran see new igrf13.f code")

This ought to be the end of the coding. Now, clean and rebuild the whole thing, and run the test suite, which ought to give no errors.