# Climate model code is so outdated, MIT starts from scratch

**URL:** <https://fortran-lang.discourse.group/t/climate-model-code-is-so-outdated-mit-starts-from-scratch/3256>\
**Category:** Uncategorized\
**Created:** [April 18, 2022, 2:19am UTC](https://fortran-lang.discourse.group/t/climate-model-code-is-so-outdated-mit-starts-from-scratch/3256 "2022-04-18T02:19:49Z")\
**Posts on this page:** 1\
**Showing post:** 25

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [April 22, 2022, 2:32pm UTC](https://fortran-lang.discourse.group/t/climate-model-code-is-so-outdated-mit-starts-from-scratch/3256/25 "2022-04-22T14:32:51Z")

</div>

> [@mecej4](#):
>
> For a rather extreme example of _expressiveness_ versus _readability_, consider a program to sieve and print primes up to 100. In APL, the program is
> 
> ```auto
> (2=+⌿0=(⍳X)∘.|⍳X)/⍳X←100
> 
> ```
> 
> The corresponding C program is about forty lines, and you can see an explanation of the APL code at [SieveOfEratosthenes - APL Wiki](https://old.aplwiki.com/SieveOfEratosthenes) .

Here is a Fortran equivalent, which I [showed earlier](https://fortran-lang.discourse.group/t/79-languages-speed-competition-can-we-make-fortran-win/3038/24). After defining some general-purpose operators as in APL, it is a two-liner.

```auto
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:

```auto
 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](https://www.dyalog.com/mastering-dyalog-apl.htm). I wrote a

> **Fortran program -- compiles with ifort**
>
> ```auto
> 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](https://kx.com/download/). I want to see how much of [Q for Mortals](https://code.kx.com/q4m3/) (also [videos](https://www.youtube.com/watch?v=8eoysfqO3UY)) can be translated.

---

_[View the full topic](https://fortran-lang.discourse.group/t/climate-model-code-is-so-outdated-mit-starts-from-scratch/3256)._