Compiler Explorer

Source code

module stdlib_kinds
  !use iso_fortran_env, only: int8, int16, int32, int64
  !use iso_c_binding, only: c_bool, c_char
  implicit none
  private
  public :: sp, dp, lk!, xdp, qp, int8, int16, int32, int64, lk, c_bool, c_char

!> Single precision real numbers
  integer, parameter :: sp = selected_real_kind(6)

!> Double precision real numbers
  integer, parameter :: dp = selected_real_kind(15)

!> Default logical kind parameter
  integer, parameter :: lk = kind(.true.)

end module stdlib_kinds

module differentials
  use stdlib_kinds
  
  contains

subroutine dx1(f,e1,dfdx1)
      real(dp), intent(in   ), dimension(:,:) :: f, e1
      real(dp), intent(  out), dimension(:,:) :: dfdx1
      integer :: i,j
      dfdx1 = 0._dp

do i=2,size(dfdx1,1)
        do j=1,size(dfdx1,2)
          dfdx1(i,j) = ( ( f(i,j) - f(i-1,j) ) / e1(i,j) ) 
        end do
      end do
    end subroutine dx1

subroutine dx2(f,e2,dfdx2)
      real(dp), intent(in   ), dimension(:,:) :: f, e2
      real(dp), intent(  out), dimension(:,:) :: dfdx2
      integer :: i,j
      dfdx2 = 0._dp

do i=1,size(dfdx2,1)
        do j=2,size(dfdx2,2)
          dfdx2(i,j) = ( f(i,j) - f(i,j-1) ) / e2(i,j)  
        end do
      end do
    end subroutine dx2

subroutine divergence(u,v,e1,e2,div)
      real(dp), intent(in   ), dimension(:,:) :: u,v,e1,e2
      real(dp), intent(  out), dimension(:,:) :: div
      integer :: i,j
      div = 0._dp

do i=2,size(div,1)
        do j=2,size(div,2)
          div(i,j) = ( ( u(i,j) - u(i-1,j) ) / e1(i,j) ) + &
                     ( ( v(i,j) - v(i,j-1) ) / e2(i,j) ) 
        end do
      end do
    end subroutine divergence

end module differentials

!
! Scratches for a grid and topology model.
!
module grid2d
  use stdlib_kinds
  type :: grid
    !
    real(dp), pointer, dimension(:,:) :: e1 => null()
    real(dp), pointer, dimension(:,:) :: e2 => null()
    !
    type(grid), pointer :: upper => null()
    type(grid), pointer :: lower => null()
    type(grid), pointer :: right => null()
    type(grid), pointer :: left => null()
    !
  contains
    procedure :: div2  ! procedure declaration
  end type
  !
  ! Here I should have some association routine
  !
contains
  !
  subroutine div2(grid_in,u,v,diver)
    use differentials, only: divergence
    class(grid), intent(in   )                 :: grid_in
    real(dp),    intent(in   ), dimension(:,:) :: u,v
    real(dp),    intent(  out), dimension(:,:) :: diver

if (.not.associated(grid_in%right)) error stop "grid_in%right not associated"
    if (.not.associated(grid_in%upper)) error stop "grid_in%upper not associated"

if (.not.associated(grid_in%right%e1)) error stop "grid_in%right%e1 not associated"
    if (.not.associated(grid_in%upper%e2)) error stop "grid_in%upper%e2 not associated"

call divergence(u,v,grid_in%right%e1,grid_in%upper%e2,diver)
    
  end subroutine
  !
end module grid2d

module topology
  use stdlib_kinds
  use grid2d

type :: ArakawaC
    type(grid) :: Ugrid
    type(grid) :: Vgrid
    type(grid) :: Tgrid
  contains
    procedure :: div  ! procedure declaration
  end type

contains
  subroutine div(topol,u,v,diver)
    use differentials, only: divergence
    class(ArakawaC), intent(in   )                 :: topol
    real(dp),         intent(in   ), dimension(:,:) :: u,v
    real(dp),         intent(  out), dimension(:,:) :: diver
    diver = 0._dp
    call divergence(u,v,topol%Ugrid%e1,topol%Vgrid%e2,diver)
  end subroutine

end module topology

module fields
  use grid2d
  implicit none
  type :: field2D
    real(dp), allocatable, dimension(:,:) :: data
    type(grid), pointer :: gridtype => null()
  end type
contains
  
end module fields

module dims
  integer, parameter :: nx = 1000
  integer, parameter :: ny = 1000
end module dims

program main
  use dims
  use grid2d
  use topology
  use differentials

! Here I define metric terms
  real(dp), target, dimension(nx,ny) :: dx
  real(dp), target, dimension(nx,ny) :: dy
  ! Here I define the velocity field
  real(dp), dimension(nx,ny) :: u
  real(dp), dimension(nx,ny) :: v
  ! Here I define the divergence
  real(dp), dimension(nx,ny) :: divU

real(dp) :: tic, toc, meantime
  real(dp) :: rand
  integer :: i,j,t
  integer :: Ntimes

! Here I define the topology
  type(ArakawaC) :: topo
  type(grid), target :: Ugrid
  type(grid), target :: Vgrid
  type(grid), target :: Tgrid
  character(len=:), allocatable :: type_operation

meantime=0._dp

dx = 1._dp
  dy = 1._dp

topo%Ugrid%e1 => dx
  topo%Vgrid%e2 => dy
  Ugrid%e1 => dx
  Vgrid%e2 => dy
  Tgrid%upper => Vgrid
  Tgrid%right => Ugrid
  !
  Ntimes = 100

!type_operation="classical"
  !type_operation="Arakawa"    
  type_operation="derived"

do t=1,Ntimes
  call random_number(rand)
  do i = 1, nx
    v(i,:) = [(rand*t*j, j=1,ny)]
  end do
  do j = 1, ny
    u(:,j) = [(rand*t*i, i=1,nx)]
  end do  
  select case ( type_operation )
    case("classical")
      call cpu_time(tic)
      call divergence(u,v,dx,dy,divU)
      call cpu_time(toc)
      !print *, "Passing array:", toc-tic
      meantime = meantime + toc-tic
    case("Arakawa")
      call cpu_time(tic)
      call divergence(u,v,topo%Ugrid%e1,topo%Vgrid%e2,divU)
      !!call topo%div(u,v,divU)
      call cpu_time(toc)
      !print *, "Arakawa:", toc-tic
      meantime = meantime + toc-tic
    case("derived")
      call cpu_time(tic)
      call divergence(u,v,Ugrid%e1,Vgrid%e2,divU)
      !!call Tgrid%div2(u,v,divU)
      call cpu_time(toc)
      !print *, "Derived:", toc-tic
      meantime = meantime + toc-tic
  end select
end do

select case ( type_operation )
    case("classical")
      print *, "Passing array average time:", meantime/Ntimes
    case("derived")
      print *, "Arakawa average time:", meantime/Ntimes
    case("staggered")
      print *, "Derived average time:", meantime/Ntimes
  end select

stop
end program main