From 8bb82d0313b0646170a0927028f4f8e2729a4507 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 11 Jan 2024 16:39:47 +0100 Subject: [PATCH] Fortran API: computing stencils for StructuredGrid --- doc/CMakeLists.txt | 4 + doc/example-fortran/CMakeLists.txt | 10 + .../structured-grid-stencils/CMakeLists.txt | 12 + .../structured-grid-stencils/program.F90 | 208 +++++++++++ src/atlas/CMakeLists.txt | 2 + .../detail/grid/StencilComputerInterface.cc | 48 +++ .../detail/grid/StencilComputerInterface.h | 27 ++ src/atlas_f/CMakeLists.txt | 7 +- src/atlas_f/atlas_module.F90 | 5 + .../grid/atlas_StencilComputer_module.F90 | 339 ++++++++++++++++++ 10 files changed, 661 insertions(+), 1 deletion(-) create mode 100644 doc/example-fortran/CMakeLists.txt create mode 100644 doc/example-fortran/structured-grid-stencils/CMakeLists.txt create mode 100644 doc/example-fortran/structured-grid-stencils/program.F90 create mode 100644 src/atlas/grid/detail/grid/StencilComputerInterface.cc create mode 100644 src/atlas/grid/detail/grid/StencilComputerInterface.h create mode 100644 src/atlas_f/grid/atlas_StencilComputer_module.F90 diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 6c7fa7d3a..05c754c4b 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -6,6 +6,10 @@ # granted to it by virtue of its status as an intergovernmental organisation nor # does it submit to any jurisdiction. +if(atlas_HAVE_FORTRAN) + add_subdirectory(example-fortran) +endif() + if(atlas_HAVE_DOCS) add_subdirectory(user-guide) diff --git a/doc/example-fortran/CMakeLists.txt b/doc/example-fortran/CMakeLists.txt new file mode 100644 index 000000000..564a289f0 --- /dev/null +++ b/doc/example-fortran/CMakeLists.txt @@ -0,0 +1,10 @@ +# (C) Copyright 2024- ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + +add_subdirectory(structured-grid-stencils) + diff --git a/doc/example-fortran/structured-grid-stencils/CMakeLists.txt b/doc/example-fortran/structured-grid-stencils/CMakeLists.txt new file mode 100644 index 000000000..c3a221662 --- /dev/null +++ b/doc/example-fortran/structured-grid-stencils/CMakeLists.txt @@ -0,0 +1,12 @@ +# (C) Copyright 2024- ECMWF. +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation nor +# does it submit to any jurisdiction. + +ecbuild_add_executable( TARGET example-fortran-structured-grid-stencils + OUTPUT_NAME program NOINSTALL + SOURCES program.F90 + LIBS atlas_f ) diff --git a/doc/example-fortran/structured-grid-stencils/program.F90 b/doc/example-fortran/structured-grid-stencils/program.F90 new file mode 100644 index 000000000..179fcd51a --- /dev/null +++ b/doc/example-fortran/structured-grid-stencils/program.F90 @@ -0,0 +1,208 @@ +! (C) Copyright 2023 ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation nor +! does it submit to any jurisdiction. + +! ------------------------------------------------------------------------------------------------ +! Example program on how to access and use StructuredColumns functionspace +! Also includes computation of stencils around a given coordinate +! ------------------------------------------------------------------------------------------------ + +subroutine run + +use, intrinsic :: iso_fortran_env, only : real64, real32 +use atlas_module, only : & + & atlas_StructuredGrid ,& + & atlas_functionspace_StructuredColumns ,& + & atlas_Field, & + & ATLAS_KIND_IDX, & + & atlas_StructuredGrid_ComputeNorth ,& + & atlas_StructuredGrid_ComputeWest ,& + & atlas_StructuredGrid_ComputeStencil ,& + & atlas_StructuredGrid_Stencil, & + & atlas_real + +use fckit_mpi_module, only : fckit_mpi + +! ------------------------------------------------ +implicit none + +character(len=:), allocatable :: gridname +type(atlas_StructuredGrid) :: grid +type(atlas_functionspace_StructuredColumns) :: functionspace +real(real64) :: zlon, zlat +integer(ATLAS_KIND_IDX) :: jlon,jlat, jgp +real(real64), pointer :: xy(:,:) +integer :: log_rank + +! ------------------------------------------------ + +log_rank = 1 ! First partition has rank=0 !!! +if (fckit_mpi%size() == 1) then + log_rank = 0 +endif + +gridname = "O320" + +! Create grid and write some properties +grid = atlas_StructuredGrid(gridname) + +if (fckit_mpi%rank() == log_rank) then + write(0,*) "grid%name() : ", grid%name() + write(0,*) "grid%size() : ", grid%size() + write(0,*) "grid%ny() : ", grid%ny() + write(0,*) "grid%nxmax() : ", grid%nxmax() +endif + +! Create functionspace of grid, distributed across MPI tasks with default partitioner, and with halo of 2 wide +! The default partitioner is "equal_regions" +functionspace = atlas_functionspace_StructuredColumns(grid, halo=2) + +! Latitude bounds without and with halos of this partition +if (fckit_mpi%rank() == log_rank) then + write(0,*) "functionspace%size_owned() : ", functionspace%size_owned() + write(0,*) "functionspace%size() : ", functionspace%size() + write(0,'(A,I0,A,I0)') " halo points between index ", functionspace%size_owned()+1 , " and " , functionspace%size() + write(0,*) "functionspace%j_begin(), functionspace%j_end() : ", & + & functionspace%j_begin(), functionspace%j_end() + write(0,*) "functionspace%j_begin_halo(), functionspace%j_end_halo() : ", & + & functionspace%j_begin_halo(), functionspace%j_end_halo() +endif + +! Longitude bounds without and with halos of this partition for the first latitude +jlat = functionspace%j_begin() +if (fckit_mpi%rank() == log_rank) then + write(0,*) "functionspace%i_begin(jlat), functionspace%i_end(jlat) : ", functionspace%i_begin(jlat), & + & functionspace%i_end(jlat) + write(0,*) "functionspace%i_begin_halo(jlat), functionspace%i_end_halo(jlat) : ", functionspace%i_begin_halo(jlat),& + & functionspace%i_end_halo(jlat) +endif + +! Access i,j indices of grid points +block + type(atlas_Field) :: field_index_i, field_index_j + integer(ATLAS_KIND_IDX), pointer :: index_i(:), index_j(:) + field_index_i = functionspace%index_i() + field_index_j = functionspace%index_j() + call field_index_i%data(index_i) + call field_index_j%data(index_j) + if (fckit_mpi%rank() == log_rank) then + write(0,*) "i,j of first partition point :", index_i(1), index_j(1) + endif + call field_index_i%final() + call field_index_j%final() +end block + +! Access to xy coordinates +block + type(atlas_Field) :: field_xy + field_xy = functionspace%xy() + call field_xy%data(xy) + call field_xy%final() +end block + +! Creating a horizontal field and perform halo exchange +block + type(atlas_Field) :: field + real(real32), pointer :: view(:) + field = functionspace%create_field(name="myfield", kind=atlas_real(real32)) + call field%data(view) + if (fckit_mpi%rank() == log_rank) then + write(0,*) "shape( horizontal field )", shape(view) + endif + view(1:functionspace%size_owned()) = 1. + view(functionspace%size_owned()+1:functionspace%size()) = 0 + call field%set_dirty() ! This marks that halos are in "dirty" state + call field%halo_exchange() ! Only exchanges halos when halos are marked as "dirty" + if (fckit_mpi%rank() == log_rank) then + write(0,*) "halo exhange success : ", all(view == 1.) + endif + call field%final() +end block + +! Creating a horizontal/vertical field +block + type(atlas_Field) :: field + real(real32), pointer :: view(:,:) + field = functionspace%create_field(name="myfield", kind=atlas_real(real32), levels=10) + call field%data(view) + if (fckit_mpi%rank() == log_rank) then + write(0,*) "shape( horizontal/vertical field )", shape(view) + endif + call field%final() +end block + +! Set a coordinate somewhere south-east of the first grid point, but still within this partition +jgp = 1 +zlon = xy(1,jgp) + 0.1 +zlat = xy(2,jgp) - 0.1 + +! Compute nearest points to the north-west of a coordinate +block + type(atlas_StructuredGrid_ComputeNorth) :: compute_north + type(atlas_StructuredGrid_ComputeWest) :: compute_west + compute_north = atlas_StructuredGrid_ComputeNorth(grid, halo=2) + compute_west = atlas_StructuredGrid_ComputeWest(grid, halo=2) + + jlat = compute_north%execute(zlat) + jlon = compute_west%execute(zlon, jlat) + jgp = functionspace%index(jlon,jlat) ! gridpoint index in 1-dimensional array + + if (fckit_mpi%rank() == log_rank) then + write(0,'(A,F6.2,A,I0)') "compute_north%execute(y=",zlat,") : ", jlat + write(0,'(A,F6.2,A,I0,A,I0)') "compute_west%execute(x=",zlon,", j=", jlat,") : ", jlon + write(0,'(A,I0,A,F6.2,A,F6.2,A)') "gridpoint north-west: jpg=",jgp," (lon,lat)=(", xy(1,jgp), ",", xy(2,jgp), ")" + endif + call compute_west%final() + call compute_north%final() +end block + +! Stencil computer of 4x4 stencil around a coordinate +block + type(atlas_StructuredGrid_ComputeStencil) :: compute_stencil + type(atlas_StructuredGrid_Stencil) :: stencil + call compute_stencil%setup(grid, stencil_width=4) + call compute_stencil%execute(zlon,zlat,stencil) + if (fckit_mpi%rank() == log_rank) then + write(0,'(A,A,dt)') 'stencil:', new_line('a'), stencil + endif + if (fckit_mpi%rank() == log_rank) then + write(0,'(A,A,dt)') 'stencil gridpoints:' + do jlat=1,stencil%width + do jlon=1,stencil%width + write(0,'(I8)',advance='no') functionspace%index(stencil%i(jlon,jlat),stencil%j(jlat)) + enddo + write(0,'(A)') "" + enddo + + write(0,'(A,A,dt)') 'stencil coordinates:' + do jlat=1,stencil%width + jgp = functionspace%index(stencil%i(1,jlat),stencil%j(jlat)) + write(0,'(A, F6.2,A)', advance='no') " lat : ", xy(2,jgp), " --- lon : " + do jlon=1,stencil%width + jgp = functionspace%index(stencil%i(jlon,jlat),stencil%j(jlat)) + write(0,'(F8.2)',advance='no') xy(1,jgp) + enddo + write(0,'(A)') "" + enddo + + endif + call compute_stencil%final() +end block + +call functionspace%final() +call grid%final() +end subroutine + +! ------------------------------------------------------------------------------------------------ + +program main +use atlas_module, only : atlas_library +implicit none +call atlas_library%initialize() +call run() +call atlas_library%finalize() +end program \ No newline at end of file diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 3636d74d6..1f8306cf9 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -173,6 +173,8 @@ grid/detail/grid/RegionalVariableResolution.h grid/detail/grid/RegionalVariableResolution.cc grid/detail/grid/Healpix.h grid/detail/grid/Healpix.cc +grid/detail/grid/StencilComputerInterface.h +grid/detail/grid/StencilComputerInterface.cc grid/detail/tiles/Tiles.cc grid/detail/tiles/Tiles.h diff --git a/src/atlas/grid/detail/grid/StencilComputerInterface.cc b/src/atlas/grid/detail/grid/StencilComputerInterface.cc new file mode 100644 index 000000000..4e1df4765 --- /dev/null +++ b/src/atlas/grid/detail/grid/StencilComputerInterface.cc @@ -0,0 +1,48 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "StencilComputerInterface.h" + +extern "C" { + atlas::grid::ComputeNorth* atlas__grid__ComputeNorth__new(const atlas::StructuredGrid::Implementation* grid, int halo) { + return new atlas::grid::ComputeNorth{atlas::StructuredGrid{grid}, halo}; + } + + void atlas__grid__ComputeNorth__delete(atlas::grid::ComputeNorth* This) { + delete This; + } + + atlas::idx_t atlas__grid__ComputeNorth__execute_real32(const atlas::grid::ComputeNorth* This, float y) { + return This->operator()(y); + } + + atlas::idx_t atlas__grid__ComputeNorth__execute_real64(const atlas::grid::ComputeNorth* This, double y) { + return This->operator()(y); + } + + + atlas::grid::ComputeWest* atlas__grid__ComputeWest__new(const atlas::StructuredGrid::Implementation* grid, int halo) { + return new atlas::grid::ComputeWest{atlas::StructuredGrid{grid}, halo}; + } + + void atlas__grid__ComputeWest__delete(atlas::grid::ComputeWest* This) { + delete This; + } + + atlas::idx_t atlas__grid__ComputeWest__execute_real32(const atlas::grid::ComputeWest* This, float x, atlas::idx_t j) { + return This->operator()(x, j); + } + + atlas::idx_t atlas__grid__ComputeWest__execute_real64(const atlas::grid::ComputeWest* This, double x, atlas::idx_t j) { + return This->operator()(x, j); + } + + +} diff --git a/src/atlas/grid/detail/grid/StencilComputerInterface.h b/src/atlas/grid/detail/grid/StencilComputerInterface.h new file mode 100644 index 000000000..87e4172d2 --- /dev/null +++ b/src/atlas/grid/detail/grid/StencilComputerInterface.h @@ -0,0 +1,27 @@ +/* + * (C) Copyright 2023 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include "atlas/grid/StencilComputer.h" +#include "atlas/grid/StructuredGrid.h" + +extern "C" { + atlas::grid::ComputeNorth* atlas__grid__ComputeNorth__new(const atlas::StructuredGrid::Implementation* grid, int halo); + void atlas__grid__ComputeNorth__delete(atlas::grid::ComputeNorth* This); + atlas::idx_t atlas__grid__ComputeNorth__execute_real32(const atlas::grid::ComputeNorth* This, float y); + atlas::idx_t atlas__grid__ComputeNorth__execute_real64(const atlas::grid::ComputeNorth* This, double y); + + atlas::grid::ComputeWest* atlas__grid__ComputeWest__new(const atlas::StructuredGrid::Implementation* grid, int halo); + void atlas__grid__ComputeWest__delete(atlas::grid::ComputeWest* This); + atlas::idx_t atlas__grid__ComputeWest__execute_real32(const atlas::grid::ComputeWest* This, float x, atlas::idx_t j); + atlas::idx_t atlas__grid__ComputeWest__execute_real64(const atlas::grid::ComputeWest* This, double x, atlas::idx_t j); + +} diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index 3c3d396e0..bcfd373e8 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/src/atlas_f/CMakeLists.txt @@ -49,7 +49,7 @@ function(generate_fortran_bindings output filename) OUTPUT ${outfile} COMMAND ${PYTHON_EXECUTABLE} ${PROJECT_SOURCE_DIR}/tools/c2f.py ${CMAKE_CURRENT_SOURCE_DIR}/${filename} -o ${outfile} -m ${_PAR_MODULE} - -t '{"idx_t":"${F_IDX}","gidx_t":"${F_GIDX}"}' + -t '{"idx_t":"${F_IDX}","gidx_t":"${F_GIDX}","atlas::idx_t":"${F_IDX}","atlas::gidx_t":"${F_GIDX}"}' DEPENDS ${filename} ) set_source_files_properties(${outfile} PROPERTIES GENERATED TRUE) endfunction() @@ -137,6 +137,10 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/redistribution/detail/Redist MODULE atlas_redistribution_c_binding OUTPUT redistribution_c_binding.f90) +generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/grid/detail/grid/StencilComputerInterface.h + MODULE atlas_grid_StencilComputer_c_binding + OUTPUT grid_StencilComputer_c_binding.f90 ) + if( atlas_HAVE_ATLAS_NUMERICS ) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/numerics/Nabla.h) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/numerics/Nabla.h) @@ -230,6 +234,7 @@ ecbuild_add_library( TARGET atlas_f grid/atlas_GridDistribution_module.F90 grid/atlas_Vertical_module.F90 grid/atlas_Partitioner_module.F90 + grid/atlas_StencilComputer_module.F90 mesh/atlas_MeshBuilder_module.F90 mesh/atlas_MeshGenerator_module.F90 mesh/atlas_Mesh_module.F90 diff --git a/src/atlas_f/atlas_module.F90 b/src/atlas_f/atlas_module.F90 index 97ddb45ae..24de358e7 100644 --- a/src/atlas_f/atlas_module.F90 +++ b/src/atlas_f/atlas_module.F90 @@ -89,6 +89,11 @@ module atlas_module & atlas_RegionalGrid use atlas_Vertical_module, only :& & atlas_Vertical +use atlas_StencilComputer_module, only: & + & atlas_StructuredGrid_ComputeNorth, & + & atlas_StructuredGrid_ComputeWest, & + & atlas_StructuredGrid_ComputeStencil, & + & atlas_StructuredGrid_Stencil use atlas_functionspace_EdgeColumns_module, only: & & atlas_functionspace_EdgeColumns use atlas_functionspace_CellColumns_module, only: & diff --git a/src/atlas_f/grid/atlas_StencilComputer_module.F90 b/src/atlas_f/grid/atlas_StencilComputer_module.F90 new file mode 100644 index 000000000..8ec0a1b0f --- /dev/null +++ b/src/atlas_f/grid/atlas_StencilComputer_module.F90 @@ -0,0 +1,339 @@ +! (C) Copyright 2013 ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation nor +! does it submit to any jurisdiction. + +#include "atlas/atlas_f.h" + +module atlas_StencilComputer_module + +use fckit_shared_object_module, only: fckit_shared_object, fckit_c_deleter +use atlas_kinds_module, only : ATLAS_KIND_IDX +use, intrinsic :: iso_c_binding, only : c_ptr, c_int + +implicit none + +private :: fckit_shared_object, fckit_c_deleter +private :: c_ptr, c_int +private :: ATLAS_KIND_IDX + +public :: atlas_StructuredGrid_ComputeNorth +public :: atlas_StructuredGrid_ComputeWest +public :: atlas_StructuredGrid_ComputeStencil +public :: atlas_StructuredGrid_Stencil + +private + +!------------------------------------------------------------------------------ +TYPE, extends(fckit_shared_object) :: atlas_StructuredGrid_ComputeNorth + +! Purpose : +! ------- +! *atlas_StructuredGrid_ComputeNorth* : To compute latitude index north of latitude + +! Methods : +! ------- + +! Author : +! ------ +! 9-Jan-2024 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ +contains + + procedure, private :: atlas_StructuredGrid_ComputeNorth__execute_real32 + procedure, private :: atlas_StructuredGrid_ComputeNorth__execute_real64 + generic :: execute => atlas_StructuredGrid_ComputeNorth__execute_real32, & + & atlas_StructuredGrid_ComputeNorth__execute_real64 + +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_StructuredGrid_ComputeNorth__final_auto +#endif + +END TYPE atlas_StructuredGrid_ComputeNorth + +interface atlas_StructuredGrid_ComputeNorth + module procedure atlas_StructuredGrid_ComputeNorth__ctor +end interface + +!------------------------------------------------------------------------------ + +!------------------------------------------------------------------------------ +TYPE, extends(fckit_shared_object) :: atlas_StructuredGrid_ComputeWest + +! Purpose : +! ------- +! *atlas_StructuredGrid_ComputeWest* : To compute longitude index west of longitude at given latitude index + +! Methods : +! ------- + +! Author : +! ------ +! 9-Jan-2024 Willem Deconinck *ECMWF* + +!------------------------------------------------------------------------------ +contains + + procedure, private :: atlas_StructuredGrid_ComputeWest__execute_real32 + procedure, private :: atlas_StructuredGrid_ComputeWest__execute_real64 + generic :: execute => atlas_StructuredGrid_ComputeWest__execute_real32, & + & atlas_StructuredGrid_ComputeWest__execute_real64 + +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_StructuredGrid_ComputeWest__final_auto +#endif + +END TYPE atlas_StructuredGrid_ComputeWest + +interface atlas_StructuredGrid_ComputeWest + module procedure atlas_StructuredGrid_ComputeWest__ctor +end interface + +!------------------------------------------------------------------------------ + +integer(c_int), parameter, private :: STENCIL_MAX_WIDTH = 6 !-> maximum 6x6 stencil + +TYPE atlas_StructuredGrid_Stencil + integer(ATLAS_KIND_IDX) :: width + integer(ATLAS_KIND_IDX) :: j_begin + integer(ATLAS_KIND_IDX) :: i_begin(STENCIL_MAX_WIDTH) ! on stack +contains + procedure, pass :: write => atlas_StructuredGrid_Stencil__write + generic, public :: write(FORMATTED) => write + procedure, public :: i => atlas_StructuredGrid_Stencil__i + procedure, public :: j => atlas_StructuredGrid_Stencil__j +END TYPE atlas_StructuredGrid_Stencil + +TYPE :: atlas_StructuredGrid_ComputeStencil + integer(c_int) :: halo + integer(ATLAS_KIND_IDX) :: stencil_width + integer(ATLAS_KIND_IDX) :: stencil_offset + type(atlas_StructuredGrid_ComputeNorth) :: compute_north + type(atlas_StructuredGrid_ComputeWest) :: compute_west +contains + procedure :: setup => atlas_StructuredGrid_ComputeStencil__setup + procedure :: execute => atlas_StructuredGrid_ComputeStencil__execute_real64 + + procedure :: assignment_operator => atlas_StructuredGrid_ComputeStencil__assignment + generic, public :: assignment(=) => assignment_operator + + procedure :: final => atlas_StructuredGrid_ComputeStencil__final +END TYPE + +! Better not use, use setup member function instead ! +!interface atlas_StructuredGrid_ComputeStencil +! module procedure atlas_StructuredGrid_ComputeStencil__ctor +!end interface + + +!======================================================== +contains +!======================================================== + + +! ----------------------------------------------------------------------------- +! Destructor + +#if FCKIT_FINAL_NOT_INHERITING +ATLAS_FINAL subroutine atlas_StructuredGrid_ComputeNorth__final_auto(this) + type(atlas_StructuredGrid_ComputeNorth), intent(inout) :: this +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine +#endif + +#if FCKIT_FINAL_NOT_INHERITING +ATLAS_FINAL subroutine atlas_StructuredGrid_ComputeWest__final_auto(this) + type(atlas_StructuredGrid_ComputeWest), intent(inout) :: this +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine +#endif + + +! ----------------------------------------------------------------------------- +! Constructors + +function atlas_StructuredGrid_ComputeNorth__ctor(grid, halo) result(this) + use, intrinsic :: iso_c_binding, only : c_int + use fckit_c_interop_module, only: c_str + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__new, atlas__grid__ComputeNorth__delete + use atlas_grid_module, only : atlas_StructuredGrid + implicit none + type(atlas_StructuredGrid_ComputeNorth) :: this + type(atlas_StructuredGrid), intent(in) :: grid + integer(c_int), intent(in) :: halo + call this%reset_c_ptr( atlas__grid__ComputeNorth__new(grid%CPTR_PGIBUG_B, halo), & + & fckit_c_deleter(atlas__grid__ComputeNorth__delete) ) + call this%return() +end function + +function atlas_StructuredGrid_ComputeWest__ctor(grid, halo) result(this) + use, intrinsic :: iso_c_binding, only : c_int + use fckit_c_interop_module, only: c_str + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__new, atlas__grid__ComputeWest__delete + use atlas_grid_module, only : atlas_StructuredGrid + implicit none + type(atlas_StructuredGrid_ComputeWest) :: this + type(atlas_StructuredGrid), intent(in) :: grid + integer(c_int), intent(in) :: halo + call this%reset_c_ptr( atlas__grid__ComputeWest__new(grid%CPTR_PGIBUG_B, halo), & + & fckit_c_deleter(atlas__grid__ComputeWest__delete) ) + call this%return() +end function + +subroutine atlas_StructuredGrid_ComputeStencil__setup(this, grid, stencil_width) + use, intrinsic :: iso_c_binding, only : c_double + use atlas_grid_module, only : atlas_StructuredGrid + implicit none + class(atlas_StructuredGrid_ComputeStencil) :: this + type(atlas_StructuredGrid), intent(in) :: grid + integer(ATLAS_KIND_IDX), intent(in) :: stencil_width + this%stencil_width = stencil_width + this%halo = (stencil_width + 1) / 2 + this%stencil_offset = stencil_width - floor(real(stencil_width,c_double) / 2._c_double + 1._c_double, ATLAS_KIND_IDX) + this%compute_north = atlas_StructuredGrid_ComputeNorth(grid, this%halo) + this%compute_west = atlas_StructuredGrid_ComputeWest(grid, this%halo) +end subroutine + + + +! ---------------------------------------------------------------------------------------- + +function atlas_StructuredGrid_ComputeNorth__execute_real32(this, y) result(index) + use, intrinsic :: iso_c_binding, only : c_float + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__execute_real32 + implicit none + integer(ATLAS_KIND_IDX) :: index + class(atlas_StructuredGrid_ComputeNorth), intent(in) :: this + real(c_float), intent(in) :: y + index = atlas__grid__ComputeNorth__execute_real32(this%CPTR_PGIBUG_B, y) + 1 +end function + +function atlas_StructuredGrid_ComputeNorth__execute_real64(this, y) result(index) + use, intrinsic :: iso_c_binding, only : c_double + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeNorth__execute_real64 + implicit none + integer(ATLAS_KIND_IDX) :: index + class(atlas_StructuredGrid_ComputeNorth), intent(in) :: this + real(c_double), intent(in) :: y + index = atlas__grid__ComputeNorth__execute_real64(this%CPTR_PGIBUG_B, y) + 1 +end function +! ---------------------------------------------------------------------------------------- + +function atlas_StructuredGrid_ComputeWest__execute_real32(this, x, j) result(index) + use, intrinsic :: iso_c_binding, only : c_float + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__execute_real32 + implicit none + integer(ATLAS_KIND_IDX) :: index + class(atlas_StructuredGrid_ComputeWest), intent(in) :: this + real(c_float), intent(in) :: x + integer(ATLAS_KIND_IDX), intent(in) :: j + index = atlas__grid__ComputeWest__execute_real32(this%CPTR_PGIBUG_B, x, j-int(1,ATLAS_KIND_IDX)) + 1 +end function + +function atlas_StructuredGrid_ComputeWest__execute_real64(this, x, j) result(index) + use, intrinsic :: iso_c_binding, only : c_double + use atlas_grid_StencilComputer_c_binding, only : atlas__grid__ComputeWest__execute_real64 + implicit none + integer(ATLAS_KIND_IDX) :: index + class(atlas_StructuredGrid_ComputeWest), intent(in) :: this + real(c_double), intent(in) :: x + integer(ATLAS_KIND_IDX), intent(in) :: j + index = atlas__grid__ComputeWest__execute_real64(this%CPTR_PGIBUG_B, x, j-int(1,ATLAS_KIND_IDX)) + 1 +end function +! ---------------------------------------------------------------------------------------- + + +subroutine atlas_StructuredGrid_ComputeStencil__execute_real64(this, x, y, stencil) + use, intrinsic :: iso_c_binding, only : c_double + implicit none + class(atlas_StructuredGrid_ComputeStencil), intent(in) :: this + real(c_double), intent(in) :: x, y + type(atlas_StructuredGrid_Stencil), intent(inout) :: stencil + + integer(ATLAS_KIND_IDX) :: jj + stencil%width = this%stencil_width + + stencil%j_begin = this%compute_north%execute(y) - this%stencil_offset + do jj = 1_ATLAS_KIND_IDX, this%stencil_width + stencil%i_begin(jj) = this%compute_west%execute(x, stencil%j_begin + jj - 1) - this%stencil_offset + enddo +end subroutine +! ---------------------------------------------------------------------------------------- + +subroutine atlas_StructuredGrid_Stencil__write (stencil, unit, iotype, v_list, iostat, iomsg) + implicit none + class(atlas_StructuredGrid_Stencil), intent(in) :: stencil + INTEGER, INTENT(IN) :: unit + CHARACTER(*), INTENT(IN) :: iotype + INTEGER, INTENT(IN) :: v_list(:) + INTEGER, INTENT(OUT) :: iostat + CHARACTER(*), INTENT(INOUT) :: iomsg + integer(ATLAS_KIND_IDX) :: jlat, jlon + do jlat = 1, stencil%width + write(unit,'(a,I0,a)',advance='no',IOSTAT=iostat) " j: ", stencil%j(jlat), " i = " + do jlon = 1, stencil%width + write(unit,'(I3)',advance='no',IOSTAT=iostat) stencil%i(jlon,jlat) + enddo + write(0,'(a)',IOSTAT=iostat) new_line('a') + enddo +end subroutine + +function atlas_StructuredGrid_Stencil__j(this, j_index) result(j) + integer(ATLAS_KIND_IDX) :: j + class(atlas_StructuredGrid_Stencil), intent(in) :: this + integer(ATLAS_KIND_IDX) :: j_index + j = this%j_begin + (j_index-1) +end function + +function atlas_StructuredGrid_Stencil__i(this, i_index, j_index) result(i) + integer(ATLAS_KIND_IDX) :: i + class(atlas_StructuredGrid_Stencil), intent(in) :: this + integer(ATLAS_KIND_IDX) :: i_index + integer(ATLAS_KIND_IDX) :: j_index + i = this%i_begin(j_index) + (i_index-1) +end function + +! ---------------------------------------------------------------------------------------- + +function atlas_StructuredGrid_ComputeStencil__ctor(grid, stencil_width) result(this) + use atlas_grid_module, only : atlas_StructuredGrid + implicit none + type(atlas_StructuredGrid_ComputeStencil) :: this + type(atlas_StructuredGrid), intent(in) :: grid + integer(ATLAS_KIND_IDX), intent(in) :: stencil_width + call this%setup(grid, stencil_width) +end function + + +subroutine atlas_StructuredGrid_ComputeStencil__assignment(this, other) +implicit none +class(atlas_StructuredGrid_ComputeStencil), intent(inout) :: this +class(atlas_StructuredGrid_ComputeStencil), intent(in) :: other +call this%final() +write(0,*) "owners = ", other%compute_north%owners() +this%compute_north = other%compute_north +this%compute_west = other%compute_west +this%stencil_width = other%stencil_width +this%stencil_offset = other%stencil_offset +this%halo = other%halo +end subroutine + +subroutine atlas_StructuredGrid_ComputeStencil__final(this) + class(atlas_StructuredGrid_ComputeStencil), intent(inout) :: this + call this%compute_north%final() + call this%compute_west%final() +end subroutine + +! ---------------------------------------------------------------------------------------- + +end module atlas_StencilComputer_module