marray_base.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: marray_base.f90                                                   #
!#                                                                           #
!# Copyright (C) 2018                                                        #
!# 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  :: mask1d(:)
     LOGICAL, POINTER  :: mask2d(:,:)
     LOGICAL, POINTER, CONTIGUOUS :: mask3d(:,:,:)
  CONTAINS
     PROCEDURE :: cuboid
     PROCEDURE :: everything
     PROCEDURE :: destroy_selection
     generic :: destroy => destroy_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 :: 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 :: 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 :: destroy
  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 (idx_init) THEN
      IF (PRESENT(m)) THEN
        new_ma%DIMS(1) = m
        IF (PRESENT(n)) THEN
          new_ma%DIMS(2) = n
          new_ma%RANK = 2
        ELSE
          new_ma%DIMS(2) = 1
          new_ma%RANK = 1
        END IF
      ELSE
        new_ma%DIMS(:) = 1
        new_ma%RANK = 0
      END IF
    END IF
  END FUNCTION createmarray
  
  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
  
  SUBROUTINE assignpointers(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    IF (ASSOCIATED(this%data1d).AND.SIZE(this%data1d).GT.0) THEN
      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
      END SELECT
    END IF
  END SUBROUTINE 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
    !------------------------------------------------------------------------!
    INTEGER :: err
    !------------------------------------------------------------------------!
    IF (ASSOCIATED(ma%data1d)) THEN
      ! copy meta data
      this%RANK    = ma%RANK
      this%DIMS(:) = ma%DIMS(:)
      IF (ASSOCIATED(this%data1d)) THEN
        IF (SIZE(this%data1d).EQ.0) THEN
          ! this%data1d allocated with size zero
          DEALLOCATE(this%data1d)
          ! allocate and assign data
          ALLOCATE(this%data1d(SIZE(ma%data1d,1)),source=ma%data1d,stat=err)
          IF (err.NE.0) THEN
            print *,"ERROR in marray_base::AssignMArray_0: memory allocation failed"
            stop 1
          END IF
        ELSE IF (SIZE(this%data1d).EQ.SIZE(ma%data1d)) THEN
          ! just copy the data
          this%data1d(:) = ma%data1d(:)
        ELSE
          print *,"ERROR in marray_base::AssignMArray_0: size of input and output do not match"
          stop 1
        END IF
      ELSE
        ! allocate and copy data
        ALLOCATE(this%data1d(SIZE(ma%data1d,1)),source=ma%data1d,stat=err)
        IF (err.NE.0) THEN
          print *,"ERROR in marray_base::AssignMArray_0: memory allocation failed"
          stop 1
        END IF
      END IF
      CALL this%AssignPointers()
    ELSE IF ((.NOT.ASSOCIATED(ma%data1d)).AND.(ma%RANK.GE.0)) THEN
      this%RANK    = ma%RANK
      this%DIMS(:) = ma%DIMS(:)
      ! just allocate data of newly created marray without copying
      ALLOCATE(this%data1d(inum*jnum*knum*this%DIMS(1)*this%DIMS(2)),stat=err)
      IF (err.NE.0) THEN
        print *,"ERROR in marray_base::AssignMArray_0: memory allocation failed"
        stop 1
      END IF
      CALL this%AssignPointers()
    ELSE IF (ASSOCIATED(this%data1d)) THEN
      ! both of same type ma%data1d not associated, but this%data1d associated
      print *,"ERROR in marray_base::AssignMArray_0: ma%data1d not associated but this%data1d is"
      stop 1
    ELSE
      ! both of same type, but none of the data1d arrays allocated
      print *,"ERROR in marray_base::AssignMArray_0: ma%data1d and this%data1d not associated"
      stop 1
    END IF
  END SUBROUTINE assignmarray_0
  
  PURE SUBROUTINE assignmarray_1(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(INUM*JNUM*KNUM*this%DIMS(1)*this%DIMS(2)), INTENT(IN) :: a
    !------------------------------------------------------------------------!
    this%data1d(:) = a(:)
  END SUBROUTINE assignmarray_1

  PURE SUBROUTINE assignmarray_2(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(SIZE(this%data2d,1),SIZE(this%data2d,2)), INTENT(IN) :: a
    !------------------------------------------------------------------------!
    this%data2d(:,:) = a(:,:)
  END SUBROUTINE assignmarray_2

  PURE SUBROUTINE assignmarray_3(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(SIZE(this%data3d,1),SIZE(this%data3d,2),SIZE(this%data3d,3)), &
                          INTENT(IN) :: a
    !------------------------------------------------------------------------!
    this%data3d(:,:,:) = a(:,:,:)
  END SUBROUTINE assignmarray_3

  PURE SUBROUTINE assignmarray_4(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(SIZE(this%data4d,1),SIZE(this%data4d,2),SIZE(this%data4d,3), &
                SIZE(this%data4d,4)), INTENT(IN) :: a
    !------------------------------------------------------------------------!
    this%data4d(:,:,:,:) = a(:,:,:,:)
  END SUBROUTINE assignmarray_4

  PURE SUBROUTINE assignmarray_5(this,a)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(INOUT) :: this
    REAL, DIMENSION(SIZE(this%data5d,1),SIZE(this%data5d,2),SIZE(this%data5d,3), &
                SIZE(this%data5d,4),SIZE(this%data5d,5)), INTENT(IN) :: a
    !------------------------------------------------------------------------!
    this%data5d(:,:,:,:,:) = a(:,:,:,:,:)
  END SUBROUTINE assignmarray_5

  PURE FUNCTION addmarray_0(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a,b
    REAL, DIMENSION(SIZE(a%data1d)) :: c
    !------------------------------------------------------------------------!
    IF (SIZE(a%data1d).EQ.SIZE(b%data1d)) &
        c(:) = a%data1d(:) + b%data1d(:)
  END FUNCTION addmarray_0

  PURE FUNCTION addmarray_1(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data1d)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data1d)) :: c
    !------------------------------------------------------------------------!
    c(:) = a%data1d(:) + b(:)
  END FUNCTION addmarray_1

  PURE FUNCTION addmarray_2(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data2d,1),SIZE(a%data2d,2)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data2d,1),SIZE(a%data2d,2)) :: c
    !------------------------------------------------------------------------!
    c(:,:) = a%data2d(:,:) + b(:,:)
  END FUNCTION addmarray_2

  PURE FUNCTION addmarray_3(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data3d,1),SIZE(a%data3d,2),SIZE(a%data3d,3)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data3d,1),SIZE(a%data3d,2),SIZE(a%data3d,3)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:) = a%data3d(:,:,:) + b(:,:,:)
  END FUNCTION addmarray_3

  PURE FUNCTION addmarray_4(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data4d,1),SIZE(a%data4d,2),SIZE(a%data4d,3), &
          SIZE(a%data4d,4)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data4d,1),SIZE(a%data4d,2),SIZE(a%data4d,3), &
          SIZE(a%data4d,4)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:,:) = a%data4d(:,:,:,:) + b(:,:,:,:)
  END FUNCTION addmarray_4

  PURE FUNCTION addmarray_5(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data5d,1),SIZE(a%data5d,2),SIZE(a%data5d,3), &
          SIZE(a%data5d,4),SIZE(a%data5d,5)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data5d,1),SIZE(a%data5d,2),SIZE(a%data5d,3), &
          SIZE(a%data5d,4),SIZE(a%data5d,5)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:,:,:) = a%data5d(:,:,:,:,:) + b(:,:,:,:,:)
  END FUNCTION addmarray_5

  PURE FUNCTION multmarray_0(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a,b
    REAL, DIMENSION(SIZE(a%data1d)) :: c
    !------------------------------------------------------------------------!
    IF (SIZE(a%data1d).EQ.SIZE(b%data1d)) &
        c(:) = a%data1d(:) * b%data1d(:)
  END FUNCTION multmarray_0

  PURE FUNCTION multmarray_1(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data1d)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data1d)) :: c
    !------------------------------------------------------------------------!
    c(:) = a%data1d(:) * b(:)
  END FUNCTION multmarray_1

  PURE FUNCTION multmarray_2(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data2d,1),SIZE(a%data2d,2)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data2d,1),SIZE(a%data2d,2)) :: c
    !------------------------------------------------------------------------!
    c(:,:) = a%data2d(:,:) * b(:,:)
  END FUNCTION multmarray_2

  PURE FUNCTION multmarray_3(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data3d,1),SIZE(a%data3d,2),SIZE(a%data3d,3)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data3d,1),SIZE(a%data3d,2),SIZE(a%data3d,3)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:) = a%data3d(:,:,:) * b(:,:,:)
  END FUNCTION multmarray_3

  PURE FUNCTION multmarray_4(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data4d,1),SIZE(a%data4d,2),SIZE(a%data4d,3), &
          SIZE(a%data4d,4)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data4d,1),SIZE(a%data4d,2),SIZE(a%data4d,3), &
          SIZE(a%data4d,4)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:,:) = a%data4d(:,:,:,:) * b(:,:,:,:)
  END FUNCTION multmarray_4

  PURE FUNCTION multmarray_5(a,b) RESULT(c)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(marray_base),INTENT(IN) :: a
    REAL, DIMENSION(SIZE(a%data5d,1),SIZE(a%data5d,2),SIZE(a%data5d,3), &
          SIZE(a%data5d,4),SIZE(a%data5d,5)),INTENT(IN) :: b
    REAL, DIMENSION(SIZE(a%data5d,1),SIZE(a%data5d,2),SIZE(a%data5d,3), &
          SIZE(a%data5d,4),SIZE(a%data5d,5)) :: c
    !------------------------------------------------------------------------!
    c(:,:,:,:,:) = a%data5d(:,:,:,:,:) * b(:,:,:,:,:)
  END FUNCTION multmarray_5

  SUBROUTINE destroy(this)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    CLASS(marray_base) :: this
    !-------------------------------------------------------------------!
    IF (ASSOCIATED(this%data1d)) DEALLOCATE(this%data1d)
    NULLIFY(this%data1d,this%data2d,this%data3d,this%data4d,this%data5d)
    this%RANK = 0
    this%DIMS = 0
  END SUBROUTINE destroy

  FUNCTION createselection(idx) RESULT(new_sel)
    IMPLICIT NONE
    !-------------------------------------------------------------------!
    TYPE(selection_base) :: new_sel
    INTEGER, OPTIONAL    :: idx(6)
    !-------------------------------------------------------------------!
    INTEGER :: err
    !-------------------------------------------------------------------!
    IF (idx_init) THEN
      ! allocate 1D mask array
      ALLOCATE(new_sel%mask3d(igmin:igmax,jgmin:jgmax,kgmin:kgmax),stat=err)
      IF (err.NE.0) THEN
        print *,"ERROR in marray_base::CreateSelection: memory allocation failed"
        stop 1
      END IF
      ! create 2D & 3D pointers into the 1D mask array
      new_sel%mask1d(1:inum*jnum*knum) => new_sel%mask3d
      new_sel%mask2d(1:inum*jnum,kgmin:kgmax) => new_sel%mask3d
      IF (PRESENT(idx)) THEN
        CALL new_sel%Cuboid(idx(1),idx(2),idx(3),idx(4),idx(5),idx(6))
      ELSE
        CALL new_sel%Everything()
      END IF
    END IF
  END FUNCTION createselection

  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%mask3d)) DEALLOCATE(this%mask3d)
    NULLIFY(this%mask1d,this%mask2d)
  END SUBROUTINE destroy_selection

END MODULE marray_base_mod