Fortran from punch card era

I was working on one geophysical software, when I found one of the libraries to compute Love and Rayleigh waves is still from the times of punch cards and tapes. Its source fits the standard of fortran pads (here you may take a look on the original release of this utility in the year 1978, including "a translation of a Russian article" to "help understand the theory used").

And yes, this piece of code is still used. I tried to look for another library, but haven't found any (at least, in public domain). All the other source codes I've seen use the same tool. That's actually why do I appreciate Fortran so much: please, show me any other programming language, so that both widely used compilers (GCC and ifort) can assemble a working program from a file consisting a mixture of standards and codes written with 40 years difference! And also having a native datatype of complex variable and functions of is, which is very needed by any wave theory math.

However, it's difficult to manage this code and improve it (for example, to add seismic wave attenuation corrections etc). Therefore, I decided to rewrite it to modern day standard. 80 kilobytes is not so lot.

This code was written when jumps between code pieces was very cheap relatively to the actual computations of some value. Eigenfunctions in this code have up to 200 variables to store sines and cosines. However, there no conditional cases or cycles in our modern-day meaning. There are jumps to certain labels of positions in code. For example:

if (number) 401,402,403

if (number1.lt.number2) goto 505

If number is less than zero, code goes to label 401, if it is zero, it goes to 402, and to 403 in the rest of values. The second expression jumps to label 505 is number1 is less than number2. And also:

10       CONTINUE

Means to continue execution (go to the next line). In the case of cycle it return to the first line of it:

DO 10 I=1,20

Do it for all the i-s in between 1 and 20, and the last statement is on the line with label 10.

Now I'd like to all the people, as well as linguistics/literature necrophiliacs, who want to feel the spirit of the good old days to translate this code to any modern language:

C********************************************************************
C---------------start of loop for Period k--------------------------S
        do 9998 k=1,kmax
        call flat1(d,rho,a,b,mmax,kind)
C--------------------loop for modes---------------------------------S
        do 9997 iq=1,kmode
        if(k-1) 605,605,599
599     if(iq-2) 600,601,601
600     c1=0.90*c(k-1,1)
        goto 605
601     if(dc*(c(k-1,iq)-c(k,iq-1))) 602,603,604
602     c1=c(k,iq-1)+0.01*dc
        goto 605
603     c1=c(k,iq-1)+0.01*dc
        goto 605
604     c1=c(k-1,iq)
605     continue
        idrop=0
999     del1=dltar(c1,t1,ifunc)
80      c2=c1+dc
        idrop=0
        del2=dltar(c2,t1,ifunc)
        if(sign(1.,del1).ne.sign(1.,del2)) goto 54
81      c1=c2
        del1=del2
        qbq=0.8*b(1)
        if(c1-0.8*b(1)) 250,251,251 
251     if(c1-(b(mmax)+0.3)) 252,250,250
252     goto 80
54      continue
        call nevill (t1,c1,c2,del1,del2,ifunc,cn) 
        if(lstop.eq.0)                  go to 3006
        lstop=0
        if(k.eq.1.and.iq.eq.1)          go to 3005
        if(k.eq.1)                      go to 3004
        do 3003 i=1,k-1
        do 3000 j=1,mode
        if(i.gt.imax(j))                go to 3001
3000    continue
        jmax=mode
        go to 3002
3001    jmax=j-1
3002    continue
3003    continue
3004    if(iq.eq.1)                     go to 3005
        jmax=iq-1
3005    continue 
        go to 9999    ! to the end of the subroutine
3006    c1=cn
        if(c1-b(mmax)) 121,121,250
121     c(k,iq)=c1
C----------------ifunc=kind : Rayleigh if =2 or Love if =1-----------
        goto (6001,6002),ifunc
6002    continue
        ratio(k,iq)=dltar(c1,t1,3)
        goto 6003
6001    continue
6003    continue
        c1=c1+0.01*dc
        imax(iq)=k
        goto 9997
250     if(k*iq-1) 256,256,255
256     continue 
        write(*,*) ncount
        if(ncount.eq.1)then
        PRINT 258
        do i=1,mmax 
        print *, b(i), a(i) ,rho(i), d(i), qs_ref(i)
        enddo
258     format(1x,/22x,'Improper initial value. No zero found.')
        endif
        goto 9999  ! this goes to the end of the subroutine
9997    continue 
C-----------------end of loop for mode------------------------------E
C--------------------search of root is finished----------------------
        goto 9998
255     kmode=iq-1
        if(kmode) 9996,9996,9998
9998    continue
C-----------------end of loop for period k--------------------------E
9996    continue
C********************************************************************

That's a pretty interesting exercise which gives a very good feeling of the long way made by computers from Turing's machine to the modern day gadgets. Persons who have coded in assembly language may also see how close is this code to it. You may also notice the concept of re-use: they tried to reuse any single piece of code. That makes sense: while one can make typos in one of the lines, leading to a hard-to-find bugs (like -0.001 instead of +0.001), usage of a wrong block will be clearly seen because of weird results.
For the aforementioned piece of code I go this listing (a bit optimised for the routine is uses):

! start of loop for Period k
   do k=1,kmax
      ! kmode is = mode = 1 in our case (fundamental modes)
      do iq=1,kmode
         if (k-1 > 0) then
            if (iq-2 < 0) then
               c1=0.90*cc(k-1,1)
            else
               if(dc*(cc(k-1,iq)-cc(k,iq-1)) <= 0.) then
                  c1=cc(k,iq-1)+0.01*dc
               else
                  c1=cc(k-1,iq)
               endif
            endif
         endif
         idrop=0 ; del1=dltar(c1,t1,ifunc)
         goin = .false.
         do
            c2=c1+dc
            idrop=0 ; del2=dltar(c2,t1,ifunc)
            if (sign(1., real(del1))*sign(1., real(del2)) < 0.) then
               call nevill (t1,c1,c2,del1,del2,ifunc,cn)
               goin = .true.
               exit
            endif
            c1=c2
            del1=del2
            qbq=0.8*b(1)
            if (c1-0.8*b(1) < 0. .and. c1-(b(mmax)+0.3) >= 0.) exit
         enddo
         if (goin) then
            if (lstop /= 0) then
               lstop=0
               if(k /= 1) then
                  do i=1,k-1
                     jmax=mode
                     do j=1,mode
                        if(i > imax(j)) then
                           jmax=j-1
                           exit
                        endif
                     enddo
                  enddo
               endif
               if(iq.eq.1) return
               jmax=iq-1
            endif
            c1 = cn
            if(c1-b(mmax) <= 0) then
               cc(k,iq)=c1
               if (ifunc == 2) ratio(k,iq)=dltar(c1,t1,3) ! 2 = Ray, 1 = Love
               c1=c1+0.01*dc
               imax(iq)=k
               cycle
            endif
         endif
         if(k*iq-1>0) exit
         write(*,*) ncount
         if(ncount == 1)then
            PRINT *, 'Improper initial value. No zero found in calcul.f'
            do i=1,mmax
               print *, b(i), a(i) ,rho(i), d(i), qs_ref(i)
            enddo
         endif
         return 
      enddo
   enddo ! end of loop for k

I think humankind needs a new wave of epic poetry. A programmer will bring to the Machine a list of code captured in a fight with with hordes of wild punch cards - similar to a knight bringing a head of defeated dragon. The hero will be also a necromancer and practice sorcery.

And yeah, have a nice Saturday evening!