Climate model code is so outdated, MIT starts from scratch

Here is a Fortran equivalent, which I showed earlier. After defining some general-purpose operators as in APL, it is a two-liner.

module apl_mod
implicit none
interface operator(.i.)
   module procedure irange
end interface
!
interface operator(.d.)
   module procedure drop
end interface
!
interface operator(.o.)
   module procedure outer_product
end interface
!
interface operator (.ni.)
   module procedure not_in
end interface
!
contains
pure function irange(i) result(v)
! return integers from 1 to i
integer, intent(in) :: i
integer             :: v(i)
integer             :: j
do j=1,i
   v(j) = j
end do
end function irange
!
pure function drop(v) result(dv)
! return v(:) without the first element
integer, intent(in) :: v(:)
integer             :: dv(size(v)-1)
if (size(v) > 1) dv = v(2:)
end function drop
!
pure function outer_product(v) result(m)
! multiplication outer product
integer, intent(in) :: v(:)
integer             :: m(size(v),size(v))
integer             :: i,j
do i=1,size(v)
   do j=1,size(v)
      m(i,j) = v(i)*v(j)
   end do
end do
end function outer_product
!
pure function not_in(a,b) result(c)
! return in c the elements of a(:) not in b(:)
integer, intent(in)  :: a(:),b(:)
integer, allocatable :: c(:)
logical              :: amask(size(a))
integer              :: i
do i=1,size(a)
   amask(i) = .not. any(b == a(i)) 
end do
c = pack(a,amask)
end function not_in
end module apl_mod
!
program main
use apl_mod
implicit none
associate (t => .d. (.i. 100) )
   print "(20(1x,i0))",t .ni. [.o. t]
end associate
end program main

output:

 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
 73 79 83 89 97

Fortran lets you define your own operators, unlike C++ or Python I believe. In C++ you can overload existing operators. For interactive use as in LFortran, abbreviations for transpose and matmul would be convenient.

There is a free 818-page book Mastering Dyalog APL. I wrote a

Fortran program -- compiles with ifort
module m
implicit none
interface operator (.plus.)
   module procedure plus
end interface operator (.plus.)
interface operator (.times.)
   module procedure times
end interface operator (.times.)
interface operator (.in.)
   module procedure in_int
end interface operator (.in.)
interface
   elemental function f_int_int(i,j) result(k)
   integer, intent(in) :: i,j
   integer             :: k
   end function f_int_int
end interface
contains
!
function in_int(x,y) result(x_in_y)
integer, intent(in) :: x(:),y(:)
logical             :: x_in_y(size(x))
integer             :: i
do i=1,size(x)
   x_in_y(i) = any(y == x(i))
end do
end function in_int
!
pure function pos(text,char) result(ipos)
! find the postions of char in text
character (len=*), intent(in) :: text
character (len=*), intent(in) :: char
integer, allocatable          :: ipos(:)
integer                       :: i,nfound,nlen
nlen = len(text)
if (nlen == 0 .or. len(char) /= 1) then
   allocate (ipos(0))
   return
end if
allocate (ipos(nlen),source=0)
nfound = 0
do i=1,nlen
   if (text(i:i) == char) then
      nfound = nfound + 1
      ipos(nfound) = i
   end if
end do
ipos = ipos(:nfound)
end function pos
!
subroutine disp(mat)
integer, intent(in) :: mat(:,:)
integer             :: i
do i=1,size(mat,1)
   print "(*(1x,i6))",mat(i,:)
end do
end subroutine disp  
!
elemental function mult(i,j) result(k)
integer, intent(in) :: i,j
integer             :: k
k = i*j
end function mult
!
function outerprod(i,j,func) result(k)
integer, intent(in)  :: i(:) ! (ni)
integer, intent(in)  :: j(:) ! (nj)
procedure(f_int_int) :: func
integer              :: k(size(i),size(j))
integer              :: ii,jj
do ii=1,size(i)
   do jj=1,size(j)
      k(ii,jj) = func(i(ii),j(jj))
   end do
end do
end function outerprod
!
function tpos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
   if (tf(i)) then
      j = j + 1
      ipos(j) = i
   end if
end do
end function tpos
!
pure integer function minus(v) result(z)
! subtract from right to left
integer, intent(in) :: v(:)
integer             :: i,n,w
n = size(v)
if (n == 0) then
   z = 0
   return
end if
z = v(n)
do i=n-1,1,-1
   z = v(i) - z
end do
end function minus
!
pure integer function plus(i,j)
integer, intent(in) :: i,j
plus = i+j
end function plus
!
pure integer function times(i,j)
integer, intent(in) :: i,j
times = i*j
end function times
!
pure function iota(n) result(v) ! p12 built-in APL function
integer, intent(in) :: n
integer             :: v(n)
integer             :: i
v = [(i,i=1,n)]
end function iota
!
elemental function int(tf) result(i)
logical, intent(in) :: tf
integer :: i
i = merge(1,0,tf)
end function int
!
end module m

program main
use m
implicit none
integer, parameter :: n1 = 4, n2 = 2
integer, allocatable :: val(:)
character (len=*), parameter :: fmt_ci = "(a,*(1x,i0))", fmt_cr = "(a,*(1x,f0.2))"
! p5
associate (price => [5.2, 11.5, 3.6, 4.0, 8.45], &
           qty   => [2, 1, 3, 6, 2], vat => 19.6)
associate (costs => price * qty)
print fmt_cr,"costs",costs
print fmt_cr,"sum(costs)",sum(costs)
associate (forecast => reshape( &
   [150,300,100,50,200,330,250,120],[n1,n2]), &
   actual => reshape([141,321,118,43,188,306,283,91],[n1,n2])) ! p8
print fmt_ci,"error",actual-forecast
print "(a,*(1x,f0.4))","VAT tax",price*vat/100 ! p7
print*,"max(75.6,87.3) =",max(75.6,87.3) ! p7
print fmt_ci,"max([11,28,52,14],[30,10,50,20]) =",max([11,28,52,14],[30,10,50,20])
print fmt_ci,"max([11,28,52,14],20) =",max([11,28,52,14],20)
associate (vec => [21,45,18,27,11]) ! p9-10
print fmt_ci,"vec",vec
print fmt_ci,"product(vec)",product(vec)
print fmt_ci,"maxval(vec)",maxval(vec)
val = [22,37,41,19,54,11,34]
print fmt_ci,"val",val
print fmt_cr,"mean(val)",sum(val)/real(size(val))
print fmt_cr,"mean(val)",mean(val) ! using a function
print fmt_ci,"(3 .plus. 6) .times. (5 .plus. 2)",(3 .plus. 6) .times. (5 .plus. 2) ! p11
print fmt_ci,"reduce(val,plus)",reduce(val,plus) ! p11
print fmt_ci,"val(4)=",val(4) ! p11
print fmt_ci,"val([2,4,7,1,4])=",val([2,4,7,1,4]) ! p11
val([3,5,1]) = 0 ! p12
print fmt_ci,"val",val ! p12
print fmt_ci,"val(1:3)",val(1:3) ! p12
print fmt_ci,"val(iota(3))",val(iota(3)) ! p12
associate (salaries => [4225,1619,3706,2240,2076,1389,3916,3918,4939,2735], &
           categories => [3,1,3,2,2,1,3,3,3,2], rates => 0.01*[8,5,2])
associate (salary_inc => salaries*rates(categories)) 
print fmt_ci,"salaries:",salaries
print fmt_cr,"salary increases:",salary_inc
print fmt_cr,"sum of salary increases:",sum(salary_inc) ! p13
print fmt_ci,"salaries > 3000:",int(salaries > 3000) ! p14
print fmt_ci,"# salaries < 2500:",count(salaries < 2500) ! p15
print fmt_ci,"categories == 3 .and. salaries < 4000:",int(categories == 3 .and. salaries < 4000) ! p14
print fmt_ci,"[0,0,1,1] /= [0,1,0,1]:",int([0,0,1,1] /= [0,1,0,1]) ! p14
print fmt_ci,"pack([23,55,17,46,81,82,83],[1,1,0,1,0,0,1]==1):",pack([23,55,17,46,81,82,83],[1,1,0,1,0,0,1]==1) ! p15
print fmt_ci,"pack(salaries, categories == 2):",pack(salaries, categories == 2) ! p15
associate (val => [22,37,41,19,54,11,34])
print fmt_ci,"val:",val
print fmt_ci,"positions with val > 35:",tpos(val > 35) ! p16
print "(a)","outer mult of [5,4,10,3] and [8,5,15,9,11,40]:"
print fmt_ci,"pos('Panama is a canal between Atlantic and Pacific','a'):",pos("Panama is a canal between Atlantic and Pacific","a") ! p16
print fmt_ci,"int([5,7,2,8,4,9] .in. [3,4,5,6]):",int([5,7,2,8,4,9] .in. [3,4,5,6]) ! p16
call disp(outerprod([5,4,10,3],[8,5,15,9,11,40],mult)) ! p23
associate (contents => [12,56,78,74,85,96,30,22,44,66,82,27]) 
print fmt_ci,"contents:",contents
print*,"all(contents > 20), any(contents < 30), count(contents < 30) =", &
        all(contents > 20), any(contents < 30), count(contents < 30) ! p106
print fmt_ci,"minus([45,9,11,2,5]):",minus([45,9,11,2,5]) ! p106
end associate ; end associate ; end associate ; end associate ; end associate 
end associate ; end associate ; end associate
!
contains
!
pure real function mean(v)
integer, intent(in) :: v(:)
mean = sum(v)/real(size(v))
end function mean
!
end program main

to do the calculations found in the first 106 pages with similar brevity and intend to continue.
Q is an APL descendant which can be used with a regular keyboard, and the 32-bit version is free. I want to see how much of Q for Mortals (also videos) can be translated.

2 Likes