marray_base.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: marray_base.f90                                                   #
!#                                                                           #
!# Copyright (C) 2018-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE marray_base_mod
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  INTEGER, SAVE          :: igmin,igmax
  INTEGER, SAVE          :: jgmin,jgmax
  INTEGER, SAVE          :: kgmin,kgmax
  INTEGER, SAVE          :: inum,jnum,knum
  LOGICAL, SAVE          :: idx_init = .false.            
  TYPE selection_base
     INTEGER           :: imin,imax
     INTEGER           :: jmin,jmax
     INTEGER           :: kmin,kmax
     LOGICAL, POINTER, CONTIGUOUS :: &
                           mask1d(:) => null(), &    
                           mask2d(:,:) => null(), &  
                           mask3d(:,:,:) => null()   
  CONTAINS
     PROCEDURE :: init => init_selection
     PROCEDURE :: assignselection
     PROCEDURE :: cuboid
     PROCEDURE :: everything
     PROCEDURE :: destroy_selection
     generic :: ASSIGNMENT (=) => assignselection
     generic :: destroy => destroy_selection
     final   :: destructor_selection
  END TYPE selection_base
  TYPE :: marray_base
    INTEGER       :: rank = -1,dims(2) = 0
    REAL, POINTER, CONTIGUOUS :: data1d(:) => null()
    REAL, POINTER, CONTIGUOUS :: data2d(:,:) => null()
    REAL, POINTER, CONTIGUOUS :: data3d(:,:,:) => null()
    REAL, POINTER, CONTIGUOUS :: data4d(:,:,:,:) => null()
    REAL, POINTER, CONTIGUOUS :: data5d(:,:,:,:,:) => null()
  CONTAINS
    PROCEDURE :: init
    PROCEDURE :: assignpointers
    PROCEDURE :: remapbounds_0
    PROCEDURE :: remapbounds_1
    PROCEDURE :: remapbounds_2
    generic   :: remapbounds => remapbounds_0, remapbounds_1, remapbounds_2
    PROCEDURE :: assignmarray_0
    PROCEDURE :: assignmarray_1
    PROCEDURE :: assignmarray_2
    PROCEDURE :: assignmarray_3
    PROCEDURE :: assignmarray_4
    PROCEDURE :: assignmarray_5
    generic   :: ASSIGNMENT (=) => assignmarray_0 , assignmarray_1, assignmarray_2, &
                                   assignmarray_3 , assignmarray_4, assignmarray_5
    PROCEDURE :: shapesmatch
    generic   :: OPERATOR(.match.) => shapesmatch
    PROCEDURE :: addmarray_0
    PROCEDURE :: addmarray_1
    PROCEDURE :: addmarray_2
    PROCEDURE :: addmarray_3
    PROCEDURE :: addmarray_4
    PROCEDURE :: addmarray_5
    generic   :: OPERATOR (+) => addmarray_0, addmarray_1, addmarray_2, &
                                 addmarray_3, addmarray_4, addmarray_5
    PROCEDURE :: multmarray_0
    PROCEDURE :: multmarray_1
    PROCEDURE :: multmarray_2
    PROCEDURE :: multmarray_3
    PROCEDURE :: multmarray_4
    PROCEDURE :: multmarray_5
    generic   :: OPERATOR (*) => multmarray_0, multmarray_1, multmarray_2, &
                                 multmarray_3, multmarray_4, multmarray_5
    PROCEDURE :: crossproduct_0
    generic   :: OPERATOR (.x.) => crossproduct_0
    PROCEDURE :: destroy
    final     :: finalize
  END TYPE
  INTERFACE marray_base
    MODULE PROCEDURE createmarray
  END INTERFACE
  INTERFACE selection_base
    MODULE PROCEDURE createselection
  END INTERFACE
  !--------------------------------------------------------------------------!
  PUBLIC :: marray_base, &
    selection_base, &
    initmeshproperties, &
    closemeshproperties
 
  CONTAINS
 
  FUNCTION createmarray(m,n) RESULT(new_ma)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    TYPE(marray_base) :: new_ma
    INTEGER, OPTIONAL, INTENT(IN) :: m,n
    !-------------------------------------------------------------------!
    INTEGER :: err
    !-------------------------------------------------------------------!
    IF (new_ma%Init(m,n)) return ! immediately return if successful
    ! something went wrong
#ifdef DEBUG
    print *,"ERROR in marray_base::CreateMArray: marray initialization failed"
    stop 1
#endif
  END FUNCTION createmarray
 
  FUNCTION init(this,m,n) RESULT(success)
    !-------------------------------------------------------------------!
    CLASS(marray_base), INTENT(INOUT) :: this
    INTEGER, OPTIONAL, INTENT(IN) :: m,n
    !-------------------------------------------------------------------!
    LOGICAL :: success
    INTEGER :: err
    !-------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::Init: marray initialization"
#endif
    success = .false.
    IF (.NOT.idx_init) return ! with success == .false.
    ! check input parameter
    ! ATTENTION: m=0 and n=0 is permitted to allow for creation of empty
    ! mesh_arrays which is required, e.g., during initialization of compounds
    IF (PRESENT(m)) THEN
      IF (m.LT.0) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::Init: 1st dimension of mesh array should be >= 0"
#endif
        return ! with success == .false.
      END IF
      this%DIMS(1) = m
      IF (PRESENT(n)) THEN
        IF (n.LT.0) THEN
#ifdef DEBUG
          print *,"ERROR in marray_base::Init: 2nd dimension of mesh array should be >= 0"
#endif
          return ! with success == .false.
        END IF
        this%DIMS(2) = n
        this%RANK = 2
      ELSE
        this%DIMS(2) = 1
        this%RANK = 1
      END IF
    ELSE
      this%DIMS(:) = 1
      this%RANK = 0
    END IF
 
    IF (.NOT.ASSOCIATED(this%data1d)) THEN
#if DEBUG > 2
      print '(A,I2,A,2(I4))',"  creating marray with rank ",this%RANK," and dimensions ",this%DIMS(1:2)
#endif
      ! allocate memory for data1d array
      ALLOCATE(this%data1d(inum*jnum*knum*this%DIMS(1)*this%DIMS(2)),stat=err)
      IF (err.NE.0) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::Init: memory allocation failed for data1d array"
#endif
        return ! with success == .false.
      ELSE
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::Init: memory allocated for data1d, size=",SIZE(this%data1d)
#endif
 
      END IF
    ELSE
      IF (SIZE(this%data1d).NE.inum*jnum*knum*this%DIMS(1)*this%DIMS(2)) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::Init: data1d array size mismatch"
#endif
        return ! with success == .false.
      END IF
    END IF
    ! assign the 2d,3d,... pointers
    IF (SIZE(this%data1d).GT.0) THEN ! do not try pointer assignment with zero-size data1d array
      IF (.NOT.this%AssignPointers()) return ! pointer assignment failed -> return immediately
    END IF
    ! report success
    success=.true.
  END FUNCTION init
 
  SUBROUTINE initmeshproperties(igmin_,igmax_,jgmin_,jgmax_,kgmin_,kgmax_)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    INTEGER, INTENT(IN) :: igmin_,igmax_,jgmin_,jgmax_,kgmin_,kgmax_
    !-------------------------------------------------------------------!
    IF (.NOT.idx_init) THEN
      igmin = igmin_
      igmax = igmax_
      jgmin = jgmin_
      jgmax = jgmax_
      kgmin = kgmin_
      kgmax = kgmax_
      inum = igmax-igmin+1
      jnum = jgmax-jgmin+1
      knum = kgmax-kgmin+1
      idx_init = .true.
    END IF
  END SUBROUTINE initmeshproperties
  
  SUBROUTINE closemeshproperties
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    IF (idx_init) THEN
      idx_init = .false.
    END IF
  END SUBROUTINE closemeshproperties
  
  FUNCTION assignpointers(this) RESULT(success)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    LOGICAL :: success
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignPointers: assign 2d,3d,... pointers"
#endif
    success=.false.
    IF (.NOT.ASSOCIATED(this%data1d)) return
    ! assign pointers depending on rank
    IF (SIZE(this%data1d).GT.0) THEN ! exclude initialized but empty marrays (see compounds)
      SELECT CASE(this%RANK)
      CASE(0)
        this%data2d(1:inum*jnum,kgmin:kgmax) => this%data1d
        this%data3d(igmin:igmax,jgmin:jgmax,kgmin:kgmax) => this%data1d
        this%data4d(igmin:igmax,jgmin:jgmax,kgmin:kgmax,1:1) => this%data1d
        this%data5d(igmin:igmax,jgmin:jgmax,kgmin:kgmax,1:1,1:1) => this%data1d
      CASE(1)
        this%data2d(1:inum*jnum*knum,1:this%DIMS(1)) => this%data1d
        this%data3d(1:inum*jnum,kgmin:kgmax,1:this%DIMS(1)) => this%data1d
        this%data4d(igmin:igmax,jgmin:jgmax,kgmin:kgmax,1:this%DIMS(1)) => this%data1d
        this%data5d(igmin:igmax,jgmin:jgmax,kgmin:kgmax,1:this%DIMS(1),1:1) => this%data1d
      CASE(2)
        this%data2d(1:inum*jnum*knum*this%DIMS(1),1:this%DIMS(2)) => this%data1d
        this%data3d(1:inum*jnum*knum,1:this%DIMS(1),1:this%DIMS(2)) => this%data1d
        this%data4d(1:inum*jnum,kgmin:kgmax,1:this%DIMS(1),1:this%DIMS(2)) => this%data1d
        this%data5d(igmin:igmax,jgmin:jgmax,kgmin:kgmax,1:this%DIMS(1),1:this%DIMS(2)) => this%data1d
      CASE DEFAULT
#ifdef DEBUG
        print *,"ERROR in marray_base::AssignPointers: rank must be in {0,1,2}"
#endif
        return ! wrong rank
      END SELECT
    END IF
    ! report success
    success=.true.
  END FUNCTION assignpointers
 
  FUNCTION remapbounds_0(this,array) RESULT(ptr)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base)                :: this
    REAL, DIMENSION(IGMIN:,JGMIN:,KGMIN:), TARGET :: array
    !------------------------------------------------------------------------!
    REAL, DIMENSION(:,:,:), POINTER :: ptr
    !------------------------------------------------------------------------!
    INTENT(IN)                      :: array
    !------------------------------------------------------------------------!
    ptr => array
  END FUNCTION remapbounds_0
 
  FUNCTION remapbounds_1(this,array) RESULT(ptr)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base)                  :: this
    REAL, DIMENSION(IGMIN:,JGMIN:,KGMIN:,:), TARGET &
                    :: array
    !------------------------------------------------------------------------!
    REAL, DIMENSION(:,:,:,:), POINTER :: ptr
    !------------------------------------------------------------------------!
    INTENT(IN)                        :: array
    !------------------------------------------------------------------------!
    ptr => array
  END FUNCTION remapbounds_1
 
  FUNCTION remapbounds_2(this,array) RESULT(ptr)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base) :: this
    REAL, DIMENSION(IGMIN:,JGMIN:,KGMIN:,:,:), TARGET &
                       :: array
    !------------------------------------------------------------------------!
    REAL, DIMENSION(:,:,:,:,:), POINTER :: ptr
    !------------------------------------------------------------------------!
    INTENT(IN)                          :: array
    !------------------------------------------------------------------------!
    ptr => array
  END FUNCTION remapbounds_2
 
  SUBROUTINE assignmarray_0(this,ma)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    CLASS(marray_base),INTENT(IN)    :: ma
    !------------------------------------------------------------------------!
#if !defined(__GFORTRAN__) || (defined(__GFORTRAN__) && __GNUC__ >= 13)
    INTEGER :: err
#endif
    LOGICAL :: LHS_initialized
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_0: marray assignment"
#endif
    IF (.NOT.ASSOCIATED(ma%data1d).OR.ma%rank.LT.0) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::AssignMArray_0: rhs of assignment not initialized"
#endif
      return
    END IF
 
    ! allways copy rank & shape
    this%RANK = ma%RANK
    this%DIMS(:) = ma%DIMS(:)
    ! check whether lhs is already initialized
    lhs_initialized = ASSOCIATED(this%data1d)
    ! lhs of assignment uninitialized -> initialize new mesh array
    ! ATTENTION: finalization of derived types works different for
    !   GNU Fortran version < 13 , hence to prevent memory leaks, one has to point
    !   the data1d array of the lhs (this%data1d) to the already associated
    !   data1d array of the rhs (ma%data1d).
    !   Other compilers, e.g., ifort (intel) & nfort (NEC) require allocation
    !   of memory for new data1d array on the lhs, because ma on the rhs
    !   is finalized on exit and the data1d array is deallocated
#if !defined(__GFORTRAN__) || (defined(__GFORTRAN__) && __GNUC__ > 12)
    ! all compilers except for gfortran version < 13
    IF (.NOT.lhs_initialized) THEN
      ALLOCATE(this%data1d,source=ma%data1d,stat=err)
      IF (err.NE.0) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::AssignMArray_0: marray initialization failed"
#endif
        RETURN
#if DEBUG > 2
      ELSE
        print *,"DEBUG INFO in marray_base::AssignMArray_0: memory allocated for data1d, size=",SIZE(this%data1d)
#endif
      END IF
    END IF
#else
    ! only gfortran < 13
    IF (lhs_initialized) THEN
#endif
      IF (.NOT.(this.match.ma)) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::AssignMArray_0: shape mismatch"
#endif
        RETURN
      END IF
      IF (SIZE(this%data1d).NE.SIZE(ma%data1d)) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::AssignMArray_0: size mismatch"
#endif
        RETURN
      END IF
      ! copy data
      this%data1d(:) = ma%data1d(:)
#ifdef __GFORTRAN__
#if __GNUC__ < 13
    ELSE
      ! pointer assignment: only gfortran < 13
      this%data1d => ma%data1d
    END IF
#else
    ! only gfortran >= 13
    IF (.NOT.lhs_initialized) THEN
      ! destroy LHS explicitely if LHS was not initialized and ma is of derived class
      SELECT TYPE(ma)
      TYPE IS(marray_base) ! do nothing
      CLASS DEFAULT
        CALL ma%Destroy()
      END SELECT
    END IF
#endif
#endif __GFORTRAN__
    IF (.NOT.this%AssignPointers()) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::AssignMArray_0: pointer reassignment failed"
#endif
      RETURN
    END IF
  END SUBROUTINE assignmarray_0
 
#ifndef DEBUG
  PURE &
#endif
  SUBROUTINE assignmarray_1(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(:), INTENT(IN) :: a
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_1: assigning 1D Fortran array"
#endif
#ifdef DEBUG
    IF (SIZE(this%data1d).NE.SIZE(a)) THEN
      print *,"ERROR in marray_base::AssignMArray_1: size mismatch ",SIZE(this%data1d)," != ",SIZE(a)
      stop 1
    END IF
#endif
    this%data1d(:) = a(:)
  END SUBROUTINE assignmarray_1
 
#ifndef DEBUG
  PURE &
#endif
  SUBROUTINE assignmarray_2(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(:,:), INTENT(IN) :: a
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_2: assigning 2D Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data2d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AssignMArray_2: shape mismatch"
      stop 1
    END IF
#endif
    this%data2d(:,:) = a(:,:)
  END SUBROUTINE assignmarray_2
 
#ifndef DEBUG
  PURE &
#endif
  SUBROUTINE assignmarray_3(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(:,:,:), INTENT(IN) :: a
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_3: assigning 3D Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data3d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AssignMArray_3: shape mismatch"
      stop 1
    END IF
#endif
    this%data3d(:,:,:) = a(:,:,:)
  END SUBROUTINE assignmarray_3
 
#ifndef DEBUG
  PURE &
#endif
  SUBROUTINE assignmarray_4(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(:,:,:,:), INTENT(IN) :: a
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_4: assigning 4D Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data4d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AssignMArray_4: shape mismatch"
      stop 1
    END IF
#endif
    this%data4d(:,:,:,:) = a(:,:,:,:)
  END SUBROUTINE assignmarray_4
 
#ifndef DEBUG
  PURE &
#endif
  SUBROUTINE assignmarray_5(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(:,:,:,:,:), INTENT(IN) :: a
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignMArray_5: assigning 5D Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data5d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AssignMArray_5: shape mismatch"
      stop 1
    END IF
#endif
    this%data5d(:,:,:,:,:) = a(:,:,:,:,:)
  END SUBROUTINE assignmarray_5
 
  PURE FUNCTION shapesmatch(this,that) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this, that
    !------------------------------------------------------------------------!
    LOGICAL                       :: res
    !------------------------------------------------------------------------!
    res = (this%rank.EQ.that%rank).AND.all(this%dims(:).EQ.that%dims(:))
  END FUNCTION shapesmatch
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_0(this,that) RESULT(data1d)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this,that
    REAL, DIMENSION(SIZE(this%data1d)) :: data1d
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_0: adding 2 marrays"
#endif
#ifdef DEBUG
    IF (.NOT.ASSOCIATED(this%data1d)) THEN
      print *,"ERROR in marray_base::AddMArray_0: 1nd argument not initialized"
      stop 1
    END IF
    IF (.NOT.ASSOCIATED(that%data1d)) THEN
      print *,"ERROR in marray_base::AddMArray_0: 2nd argument not initialized"
      stop 1
    END IF
    IF (.NOT.(this.match.that)) THEN
      print *,"ERROR in marray_base::AddMArray_0: shape mismatch"
      stop 1
    END IF
#endif
    data1d(:) = this%data1d(:) + that%data1d(:)
  END FUNCTION addmarray_0
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_1(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data1d)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_1: adding marray to 1d Fortran array"
#endif
#ifdef DEBUG
    IF (SIZE(this%data1d).NE.SIZE(a)) THEN
      print *,"ERROR in marray_base::AddMArray_1: size mismatch"
      stop 1
    END IF
#endif
    b(:) = this%data1d(:) + a(:)
  END FUNCTION addmarray_1
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_2(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(SIZE(this%data2d,1),SIZE(this%data2d,2)),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data2d,1),SIZE(this%data2d,2)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_2: adding marray to 2d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data2d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AddMArray_2: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:) = this%data2d(:,:) + a(:,:)
  END FUNCTION addmarray_2
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_3(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:), INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data3d,1),SIZE(this%data3d,2),SIZE(this%data3d,3)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_3: adding marray to 3d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data3d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AddMArray_3: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:) = this%data3d(:,:,:) + a(:,:,:)
  END FUNCTION addmarray_3
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_4(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data4d,1),SIZE(this%data4d,2),SIZE(this%data4d,3), &
          SIZE(this%data4d,4)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_4: adding marray to 4d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data4d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AddMArray_4: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:,:) = this%data4d(:,:,:,:) + a(:,:,:,:)
  END FUNCTION addmarray_4
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION addmarray_5(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data5d,1),SIZE(this%data5d,2),SIZE(this%data5d,3), &
          SIZE(this%data5d,4),SIZE(this%data5d,5)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AddMArray_5: adding marray to 5d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data5d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::AddMArray_5: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:,:,:) = this%data5d(:,:,:,:,:) + a(:,:,:,:,:)
  END FUNCTION addmarray_5
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_0(this,that) RESULT(data1d)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    CLASS(marray_base), INTENT(IN) :: that
    REAL, DIMENSION(SIZE(this%data1d)) :: data1d
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_0: multiply 2 marrays"
#endif
#ifdef DEBUG
    IF (.NOT.ASSOCIATED(this%data1d)) THEN
      print *,"ERROR in marray_base::MultMArray_0: 1nd argument not initialized"
    END IF
    IF (.NOT.ASSOCIATED(that%data1d)) THEN
      print *,"ERROR in marray_base::MultMArray_0: 2nd argument not initialized"
    END IF
    IF (.NOT.(this.match.that)) THEN
      print *,"ERROR in marray_base::MultMArray_0: shape mismatch"
    END IF
#endif
    data1d(:) = this%data1d(:) * that%data1d(:)
  END FUNCTION multmarray_0
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_1(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data1d)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_1: multiply marray with 1d Fortran array"
#endif
#ifdef DEBUG
    IF (SIZE(this%data1d).NE.SIZE(a)) THEN
      print *,"ERROR in marray_base::MultMArray_1: size mismatch"
      stop 1
    END IF
#endif
    b(:) = this%data1d(:) * a(:)
  END FUNCTION multmarray_1
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_2(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data2d,1),SIZE(this%data2d,2)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_2: multiply marray with 2d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data2d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::MultMArray_2: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:) = this%data2d(:,:) * a(:,:)
  END FUNCTION multmarray_2
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_3(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data3d,1),SIZE(this%data3d,2),SIZE(this%data3d,3)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_3: multiply marray with 3d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data3d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::MultMArray_3: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:) =  this%data3d(:,:,:) * a(:,:,:)
  END FUNCTION multmarray_3
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_4(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data4d,1),SIZE(this%data4d,2),SIZE(this%data4d,3), &
          SIZE(this%data4d,4)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_4: multiply marray with 4d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data4d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::MultMArray_4: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:,:) =  this%data4d(:,:,:,:) * a(:,:,:,:)
  END FUNCTION multmarray_4
 
#ifndef DEBUG
  PURE &
#endif
  FUNCTION multmarray_5(this,a) RESULT(b)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    REAL, DIMENSION(:,:,:,:,:),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(this%data5d,1),SIZE(this%data5d,2),SIZE(this%data5d,3), &
          SIZE(this%data5d,4),SIZE(this%data5d,5)) :: b
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::MultMArray_5: multiply marray with 5d Fortran array"
#endif
#ifdef DEBUG
    IF (.NOT.all(shape(this%data5d).EQ.shape(a))) THEN
      print *,"ERROR in marray_base::MultMArray_5: shape mismatch"
      stop 1
    END IF
#endif
    b(:,:,:,:,:) =  this%data5d(:,:,:,:,:) * a(:,:,:,:,:)
  END FUNCTION multmarray_5
 
#ifndef DEBUG
  PURE &
#endif
    FUNCTION crossproduct_0(this,that) RESULT(data2d)
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: this
    CLASS(marray_base), INTENT(IN) :: that
    REAL, DIMENSION(SIZE(this%data2d,DIM=1),3) :: data2d
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::CrossProduct_0: compute cross product of two 3D vector like marrays"
#endif
#ifdef DEBUG
    IF (.NOT.ASSOCIATED(this%data1d)) THEN
      print *,"ERROR in marray_base::CrossProduct_0: 1nd argument not initialized"
    END IF
    IF (.NOT.ASSOCIATED(that%data1d)) THEN
      print *,"ERROR in marray_base::CrossProduct_0: 2nd argument not initialized"
    END IF
#endif
    IF (.NOT.(this.match.that)) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::CrossProduct_0: shape mismatch"
#endif
      RETURN
    END IF
    IF (this%DIMS(this%RANK).NE.3) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::CrossProduct_0: not a 3D vector"
#endif
      RETURN
    END IF
    data2d(:,1) = axb1(this%data2d(:,2),this%data2d(:,3),that%data2d(:,2),that%data2d(:,3))
    data2d(:,2) = axb1(this%data2d(:,3),this%data2d(:,1),that%data2d(:,3),that%data2d(:,1))
    data2d(:,3) = axb1(this%data2d(:,1),this%data2d(:,2),that%data2d(:,1),that%data2d(:,2))
 
    CONTAINS
      ELEMENTAL FUNCTION axb1(a2,a3,b2,b3)
        IMPLICIT NONE
        !------------------------------------------------------------------------!
        REAL, INTENT(IN)  :: a2,a3,b2,b3
        REAL :: axb1
        !------------------------------------------------------------------------!
        axb1 = a2*b3 - a3*b2
      END FUNCTION axb1
  END FUNCTION crossproduct_0
 
  SUBROUTINE destroy(this,called_from_finalize)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    CLASS(marray_base) :: this
    LOGICAL, OPTIONAL :: called_from_finalize
    !-------------------------------------------------------------------!
#if DEBUG > 2
    WRITE(*,'(A)',advance='NO') " DEBUG INFO in marray_base::Destroy called"
    IF (PRESENT(called_from_finalize)) THEN
      IF (called_from_finalize) WRITE(*,'(A)') " from Finalize"
    END IF
#endif
    IF (ASSOCIATED(this%data1d)) THEN
#if DEBUG > 2
      print *,"DEBUG INFO in marray_base::Destroy: deallocating data1d, size=",SIZE(this%data1d)
#endif
      DEALLOCATE(this%data1d)
    END IF
    NULLIFY(this%data1d,this%data2d,this%data3d,this%data4d,this%data5d)
    this%rank    =-1
    this%dims(:) = 0
  END SUBROUTINE destroy
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    TYPE(marray_base) :: this
    !-------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::Finalize called"
#endif
    CALL this%Destroy(.true.)
  END SUBROUTINE finalize
 
  FUNCTION createselection(idx) RESULT(new_sel)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    TYPE(selection_base) :: new_sel
    INTEGER, OPTIONAL    :: idx(6)
    !-------------------------------------------------------------------!
    INTEGER :: err
    !-------------------------------------------------------------------!
    IF (new_sel%Init(idx)) return ! immediately return if successful
#ifdef DEBUG
    print *,"ERROR in selection_base::CreateSelection: initialization failed"
    stop 1
#endif
  END FUNCTION createselection
 
  FUNCTION init_selection(this,idx) RESULT(success)
    !-------------------------------------------------------------------!
    CLASS(selection_base), INTENT(INOUT) :: this
    INTEGER, OPTIONAL    :: idx(6)
    !-------------------------------------------------------------------!
    LOGICAL :: success
    INTEGER :: err
    !-------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::Init_selection: selection initialization"
#endif
    success = .false.
    IF (.NOT.idx_init) return ! with success == .false.
    ! allocate 1D mask array
    ALLOCATE(this%mask1d(1:inum*jnum*knum),stat=err)
    IF (err.NE.0) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::Init_selection: memory allocation failed"
#endif
      return ! with success == .false.
    END IF
    ! set 2D & 3D pointers
    this%mask2d(1:inum*jnum,kgmin:kgmax) => this%mask1d
    this%mask3d(igmin:igmax,jgmin:jgmax,kgmin:kgmax) => this%mask1d
    IF (PRESENT(idx)) THEN
      CALL this%Cuboid(idx(1),idx(2),idx(3),idx(4),idx(5),idx(6))
    ELSE
      CALL this%Everything()
    END IF
    ! report success
    success=.true.
  END FUNCTION init_selection
 
  SUBROUTINE assignselection(this,sel)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(selection_base),INTENT(INOUT) :: this
    CLASS(selection_base),INTENT(IN)    :: sel
    !------------------------------------------------------------------------!
#if DEBUG > 2
    print *,"DEBUG INFO in marray_base::AssignSelection: selection assignment"
#endif
    IF (.NOT.ASSOCIATED(sel%mask1d)) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::AssignSelection: rhs of assignment not initialized"
#endif
      return
    END IF
 
    IF (.NOT.ASSOCIATED(this%mask1d)) THEN
      ! lhs of assignment uninitialized -> initialize new selection
      ! ATTENTION: finalization of derived types works different for
      !   GNU Fortran, hence to prevent memory leaks, one has to point
      !   the mask1d array of the lhs (this%mask1d) to the already associated
      !   mask1d array of the rhs (ma%mask1d).
      !   Other compilers, e.g., ifort (intel) & nfort (NEC) require generation
      !   of a new selection with mask1d array which is destroyed on exit.
#if defined(__GFORTRAN__) && __GNUC__ < 13
      this%imin = sel%imin
      this%imax = sel%imax
      this%jmin = sel%jmin
      this%jmax = sel%jmax
      this%kmin = sel%kmin
      this%kmax = sel%kmax
      this%mask1d => sel%mask1d
      ! set 2D & 3D pointers
      this%mask2d(1:inum*jnum,kgmin:kgmax) => this%mask1d
      this%mask3d(igmin:igmax,jgmin:jgmax,kgmin:kgmax) => this%mask1d
      ! immediately return
      return
#else
      IF (.NOT.this%Init((/sel%imin,sel%imax,sel%jmin,sel%jmax,sel%kmin,sel%kmax/))) THEN
#ifdef DEBUG
        print *,"ERROR in marray_base::AssignSelection: initialization failed"
#endif
        return
      END IF
#endif
    ELSE
      this%imin = sel%imin
      this%imax = sel%imax
      this%jmin = sel%jmin
      this%jmax = sel%jmax
      this%kmin = sel%kmin
      this%kmax = sel%kmax
    END IF
    IF (.NOT.SIZE(this%mask1d).EQ.SIZE(sel%mask1d)) THEN
#ifdef DEBUG
      print *,"ERROR in marray_base::AssignSelection: size mismatch"
#endif
      return
    END IF
    ! copy data
    this%mask1d(:) = sel%mask1d(:)
  END SUBROUTINE assignselection
 
  SUBROUTINE cuboid(this,imin,imax,jmin,jmax,kmin,kmax)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    CLASS(selection_base) :: this
    INTEGER, INTENT(IN)   :: imin,imax,jmin,jmax,kmin,kmax
    !-------------------------------------------------------------------!
    INTEGER :: i,j,k
    !-------------------------------------------------------------------!
    this%imin = min(max(imin,igmin),igmax)
    this%imax = max(min(imax,igmax),igmin)
    this%jmin = min(max(jmin,jgmin),jgmax)
    this%jmax = max(min(jmax,jgmax),jgmin)
    this%kmin = min(max(kmin,kgmin),kgmax)
    this%kmax = max(min(kmax,kgmax),kgmin)
    DO k=kgmin,kgmax
      DO j=jgmin,jgmax
        DO i=igmin,igmax
          IF ((i.GE.this%imin.AND.i.LE.this%imax).AND. &
              (j.GE.this%jmin.AND.j.LE.this%jmax).AND. &
              (k.GE.this%kmin.AND.k.LE.this%kmax)) THEN
            this%mask3d(i,j,k) = .true.
          ELSE
            this%mask3d(i,j,k) = .false.
          END IF
        END DO
      END DO
    END DO
  END SUBROUTINE cuboid
 
  SUBROUTINE everything(this)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    CLASS(selection_base) :: this
    !-------------------------------------------------------------------!
    this%imin = igmin
    this%imax = igmax
    this%jmin = jgmin
    this%jmax = jgmax
    this%kmin = kgmin
    this%imax = kgmax
    this%mask1d(:) = .true.
  END SUBROUTINE everything
 
  SUBROUTINE destroy_selection(this)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    CLASS(selection_base) :: this
    !-------------------------------------------------------------------!
    IF (ASSOCIATED(this%mask1d)) DEALLOCATE(this%mask1d)
    NULLIFY(this%mask2d,this%mask3d)
  END SUBROUTINE destroy_selection
 
  SUBROUTINE destructor_selection(this)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    TYPE(selection_base) :: this
    !-------------------------------------------------------------------!
    CALL this%Destroy_selection()
  END SUBROUTINE destructor_selection
 
END MODULE marray_base_mod