gravity_sboxspectral.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: gravity_sboxspectral.f90                                          #
!#                                                                           #
!# Copyright (C) 2015-2019                                                   #
!# Jannes Klee <jklee@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.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!
!! \todo Remove unnecessary allocations like ip_den, etc. At the moment there
!!       are certainly possibilities to save memory, because for every task
!!       a new array is allocated.
!
!! \attention The parallel version with this module needs a pencil decomposition
!!            of Mesh%INUM X Mesh%JNUM/nprocs for each process, in order to use
!!            the FFT-method with period boundaries in one direction.
!!
!! References:
!!
!! \cite gammie2001 , \cite gammiecode, \cite frigo2005, \cite mathkeisan2018
!!
!! \extends gravity_spectral
!! \ingroup gravity
!----------------------------------------------------------------------------!
MODULE gravity_sboxspectral_mod
  USE gravity_base_mod
  USE gravity_spectral_mod, ONLY : gravity_spectral
  USE boundary_base_mod
  USE fluxes_base_mod
  USE physics_base_mod
  USE mesh_base_mod
  USE logging_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE common_dict
  USE fftw
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  PRIVATE
  REAL, PARAMETER              :: sqrttwopi &
    = 2.50662827463100050241576528481104525300698674

  TYPE, EXTENDS(gravity_spectral) :: gravity_sboxspectral
#ifdef HAVE_FFTW

    REAL(C_DOUBLE), POINTER          :: mass2d(:,:)
    COMPLEX(C_DOUBLE_COMPLEX), POINTER &
                                     :: fmass2d(:,:)
    REAL, DIMENSION(:,:,:), POINTER  :: fmass2d_real
    INTEGER(C_INTPTR_T)              :: local_joff
    REAL,DIMENSION(:), POINTER       :: kx
    REAL,DIMENSION(:), POINTER       :: ky
    REAL                             :: shiftconst
    REAL, DIMENSION(:), POINTER      :: joff, jrem
    INTEGER                          :: order
    LOGICAL                          :: calc_fmass2d
#ifdef PARALLEL
    INTEGER(C_INTPTR_T)              :: c_inum, c_jnum
    INTEGER(C_INTPTR_T)              :: alloc_local, local_jnum
    TYPE(c_ptr)                      :: mass2d_pointer
    TYPE(c_ptr)                      :: fmass2d_pointer
#endif
#endif

  CONTAINS

    PROCEDURE :: initgravity_sboxspectral
    PROCEDURE :: updategravity_single
    PROCEDURE :: infogravity
    PROCEDURE :: calcdiskheight_single
    PROCEDURE :: finalize
#ifdef HAVE_FFTW
    PROCEDURE :: calcpotential
    PROCEDURE :: fft_forward
    PROCEDURE :: fft_backward
    PROCEDURE :: fieldshift
#endif
  END TYPE

  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       gravity_sboxspectral
  !--------------------------------------------------------------------------!

  CONTAINS

  SUBROUTINE initgravity_sboxspectral(this,Mesh,Physics,config,IO)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    TYPE(Dict_TYP),      POINTER    :: config,IO
    !------------------------------------------------------------------------!
    INTEGER             :: gravity_number
    INTEGER             :: err, valwrite, i
#if defined(HAVE_FFTW) && defined(PARALLEL)
    INTEGER(C_INTPTR_T) :: C_INUM,C_JNUM
    INTEGER(C_INTPTR_T) :: alloc_local, local_JNUM, local_joff
    INTEGER             :: nprocs
#endif
    !------------------------------------------------------------------------!
    CALL this%InitGravity(mesh,physics,"shearingbox spectral solver",config,io)

    !-------------- checks & warnings for initial conditions ----------------!
#if !defined(HAVE_FFTW)
    CALL this%Error("InitGravity_sboxspectral", &
         "No fftw package could be loaded.")
#else
    CALL getattr(config, "order", this%order, 2)

#ifdef PARALLEL
    CALL fftw_mpi_init()
    c_inum = mesh%INUM
    c_jnum = mesh%JNUM
#endif

    ! rotational symmetry is not allowed
    IF (mesh%ROTSYM.NE.0) &
      CALL this%Error("InitGravity_sboxspectral", &
         "Rotational symmetry not supported.")

    ! 1D simulations are not allowed
    IF (mesh%NDIMS.EQ.1) &
      CALL this%Error("InitGravity_sboxspectral", &
         "Only 2D (shearingsheet) and 3D (shearingbox) simulations allowed with this module.")

    ! check physics
    SELECT TYPE (phys => physics)
    CLASS IS(physics_eulerisotherm)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("InitGravity_sboxspectral", &
         "Physics modules with density necessary for this module.")
    END SELECT

    ! Check even number of cells
    IF(.NOT.(mod(mesh%JMAX-mesh%JMIN+1,2)==0)) THEN
      CALL this%Error("InitGravity_sboxspectral", &
                 "The spectral poisson solver needs an even number of "// &
                 "cells in the phi direction due to the discrete cosinus "// &
                 "transform.")
    END IF
    !------------------------------------------------------------------------!

    ! check the dimensions if fftw should be used in parallel in order to
    ! know how to allocate the local arrays
#ifdef PARALLEL
    CALL mpi_comm_size(mpi_comm_world, nprocs, err)
    IF (nprocs.GT.1 .AND. mesh%shear_dir.EQ.2) THEN
      CALL this%Error("InitGravity_sboxspectral", &
                 "Execution of the parallel Fourier solver is only possible with pencil "// &
                 "decomposition. The first dimension should be fully accessible by the "//&
                 "first core. The shear thereto needs to vary along the second " // &
                 "dimension. In order to achieve this set the boundaries accordingly. "// &
                 "Check InitBoundary where the shear direction is set.")
    END IF
    IF (mesh%dims(1).GT.1 .AND. mesh%shear_dir.EQ.1) THEN
      CALL this%Error("InitGravity_sboxspectral", &
                 "The first dimension needs to be fully accessible by each process. "// &
                 "Please change domain decomposition to pencil decompositon. " // &
                 "This needs to be done during initialization. In the Mesh dictionary " // &
                 "use the key 'decompositon' and the value (/ 1,-1, 1/) to force decomposition " // &
                 "along second direction")
    END IF
#endif
    ALLOCATE( &
             this%phi(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
#if !defined(PARALLEL)
             this%mass2D(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX), &
             this%Fmass2D(mesh%IMIN:mesh%IMAX/2+1,mesh%JMIN:mesh%JMAX), &
#else
             ! initialization handled by FFTW (see below)
#endif
             this%Fmass2D_real(mesh%IMIN:mesh%IMAX+2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
             this%kx(mesh%INUM), &
             this%ky(mesh%JNUM), &
             this%joff(mesh%IMIN:mesh%IMAX), &
             this%jrem(mesh%IMIN:mesh%IMAX), &
             this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
             stat=err)
    IF (err.NE.0) &
        CALL this%Error("InitGravity_sboxspectral","Memory allocation failed.")

    this%local_joff = 0

    ! use special allocation pattern from fftw when using MPI in order to
    ! assure good alignment
#if defined(PARALLEL)
    this%alloc_local     = fftw_mpi_local_size_2d(c_jnum, c_inum, &
                              mpi_comm_world, this%local_JNUM, this%local_joff)
    this%mass2D_pointer  = fftw_alloc_real(2*this%alloc_local)
    this%Fmass2D_pointer = fftw_alloc_complex(this%alloc_local)
    call c_f_pointer(this%mass2D_pointer,this%mass2D, &
                     [2*(c_inum/2+1),this%local_JNUM])
    call c_f_pointer(this%Fmass2D_pointer,this%Fmass2D, &
                     [c_inum/2+1,this%local_JNUM])
#endif

    !------------------- create plans for fftw ------------------------------!
    ! Pay attention to the argument order of the dimension (JNUM and INUM    !
    ! are switched because of C -> row-major, Fortran -> column-major),      !
    ! BUT ONLY in modern Fortran UNLIKE the legacy version                   !
    ! ------------ plans are allocated in dictionary ------------------------!
#if  !defined(PARALLEL)
    this%plan_r2c = fftw_plan_dft_r2c_2d(mesh%JMAX-mesh%JMIN+1, &
                                         mesh%IMAX-mesh%IMIN+1,this%mass2D, &
                                         this%Fmass2D,fftw_measure)
    this%plan_c2r = fftw_plan_dft_c2r_2d(mesh%JMAX-mesh%JMIN+1, &
                                         mesh%IMAX-mesh%IMIN+1, this%Fmass2D, &
                                         this%mass2D, fftw_measure)
    IF ((.NOT. c_associated(this%plan_r2c)) .OR. (.NOT. c_associated(this%plan_c2r))) THEN
       CALL this%Error("InitGravity_sboxspectral","FFT plan could not be created.")
    END IF
#elif defined(PARALLEL)
    this%plan_r2c = fftw_mpi_plan_dft_r2c_2d(c_jnum,c_inum, &
                                             this%mass2D,this%Fmass2D, &
                                             mpi_comm_world, fftw_measure)
    this%plan_c2r = fftw_mpi_plan_dft_c2r_2d(c_jnum,c_inum, &
                                             this%Fmass2D,this%mass2D, &
                                             mpi_comm_world, fftw_measure)
#endif
    !------------------------------------------------------------------------!
    ! set potential and acceleration to zero
    this%phi(:,:,:) = 0.
    this%mass2D(:,:) = 0.
    this%Fmass2D(:,:) = cmplx(0.,0)
    this%Fmass2D_real(:,:,:) = 0.
    this%joff(:) = 0.
    this%jrem(:) = 0.
    this%kx(:) = 0.
    this%ky(:) = 0.
    ! nullify not used potential explicitely
    NULLIFY(this%pot)

    ! precompute wavenumbers
    this%kx = cshift((/(i-(mesh%INUM+1)/2,i=0,mesh%INUM-1)/), &
                +(mesh%INUM+1)/2)*2.*pi/(mesh%XMAX-mesh%XMIN)
    this%ky = cshift((/(i-(mesh%JNUM+1)/2,i=0,mesh%JNUM-1)/), &
                +(mesh%JNUM+1)/2)*2.*pi/(mesh%YMAX-mesh%YMIN)

    ! precompute shift constant
    IF(mesh%shear_dir.EQ.2) THEN
      this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/mesh%dy
    ELSE IF (mesh%shear_dir.EQ.1) THEN
      this%shiftconst = mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)/mesh%dx
    ELSE
      CALL this%Error("InitGravity_sboxspectral", &
        "Either WE_shear or SN_shear need to be applied.")
    END IF

    !------------------------------- output ---------------------------------!
    valwrite = 0
    CALL getattr(config, "output/phi", valwrite, 0)
    IF (valwrite .EQ. 1) &
      CALL setattr(io, "phi", &
              this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    valwrite = 0
    this%calc_Fmass2D = .false.
    CALL getattr(config, "output/Fmass2D", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
      CALL setattr(io, "Fmass2D_real", &
              this%Fmass2D_real(mesh%IMIN:mesh%IMAX+2,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
      this%calc_Fmass2D = .true.
    END IF

#endif
  END SUBROUTINE initgravity_sboxspectral

  SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),            INTENT(IN) :: Mesh
    CLASS(physics_base),         INTENT(IN) :: Physics
    CLASS(fluxes_base),          INTENT(IN) :: Fluxes
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    REAL,                        INTENT(IN) :: time,dt
    !------------------------------------------------------------------------!
    INTEGER :: i,j,k
    REAL    :: w1,w2
    !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    ! calc potential first
    CALL this%CalcPotential(mesh,physics,time,pvar)

    IF (this%order.EQ.2) THEN
!NEC$ ivdep
      DO k = mesh%KMIN,mesh%KMAX
!NEC$ ivdep
        DO j = mesh%JMIN,mesh%JMAX
!NEC$ ivdep
          DO i = mesh%IMIN,mesh%IMAX
            ! second order approximation
            this%accel%data4d(i,j,k,1) = -1.0*(this%phi(i+1,j,k)-this%phi(i-1,j,k))/ &
                                 (2*mesh%dlx%data3d(i,j,k))
            this%accel%data4d(i,j,k,2) = -1.0*(this%phi(i,j+1,k)-this%phi(i,j-1,k))/ &
                                 (2*mesh%dly%data3d(i,j,k))
         END DO
        END DO
      END DO
    ELSE IF (this%order.EQ.4) THEN
      w1 = 3./48.
      w2 = 30./48.
!NEC$ ivdep
      DO k = mesh%KMIN,mesh%KMAX
!NEC$ ivdep
        DO j = mesh%JMIN,mesh%JMAX
!NEC$ ivdep
          DO i = mesh%IMIN,mesh%IMAX
            ! fourth order
            this%accel%data4d(i,j,k,1) = -1.0*(w1*this%phi(i-2,j,k)-w2*this%phi(i-1,j,k)+ &
                                        w2*this%phi(i+1,j,k)-w1*this%phi(i+2,j,k))/ &
                                       (mesh%dlx%data3d(i,j,k))
            this%accel%data4d(i,j,k,2) = -1.0*(w1*this%phi(i,j-2,k)-w2*this%phi(i,j-1,k)+ &
                                        w2*this%phi(i,j+1,k)-w1*this%phi(i,j+2,k)) / &
                                       (mesh%dly%data3d(i,j,k))
          END DO
        END DO
      END DO
    ELSE
      CALL this%Error("InitGravity_sboxspectral", &
        "The chosen order approximation does not exist (only 2 and 3).")
    END IF
#endif
  END SUBROUTINE updategravity_single

#ifdef HAVE_FFTW
  SUBROUTINE calcpotential(this,Mesh,Physics,time,pvar)
    USE physics_eulerisotherm_mod, ONLY : statevector_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    REAL,                INTENT(IN) :: time
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    !------------------------------------------------------------------------!
    INTEGER :: i,j,k
    REAL    :: K2max,K2,KO
    REAL    :: joff2,jrem2,delt,time0
    INTEGER :: nprocs
#ifdef PARALLEL
    INTEGER :: status(MPI_STATUS_SIZE),ierror,ier
    REAL    :: mpi_buf(Mesh%IMIN:Mesh%IMAX,Mesh%GJNUM,Mesh%KMIN:Mesh%KMAX) !TODO: ONLY 2D
#endif
    !------------------------------------------------------------------------!
    !---------------- fourier transformation of density ---------------------!
    ! calculate the shift of the indice at time t                            !
    ! shift density to pretend periodic behavior with interpolation          !
    nprocs=1
#ifdef PARALLEL
    CALL mpi_comm_size(mpi_comm_world, nprocs, ier)
#endif

    IF (mesh%OMEGA .EQ. 0) THEN
      time0 = 0.0
    ELSE IF (mesh%shear_dir.EQ.2) THEN
      time0 = nint(time*mesh%Q*mesh%OMEGA*(mesh%xmax-mesh%xmin)/ &
          (mesh%ymax-mesh%ymin))*(mesh%ymax-mesh%ymin)/(mesh%Q*mesh%OMEGA &
          *(mesh%xmax-mesh%xmin))
    ELSE IF (mesh%shear_dir.EQ.1) THEN
      time0 = nint(time*mesh%Q*mesh%OMEGA*(mesh%ymax-mesh%ymin)/ &
          (mesh%xmax-mesh%xmin))*(mesh%xmax-mesh%xmin)/(mesh%Q*mesh%OMEGA &
          *(mesh%ymax-mesh%ymin))
    END IF
    delt = time - time0

    !----------------- shift field to periodic point ------------------------!
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_eulerisotherm)
      CALL this%FieldShift(mesh,physics,delt, &
                p%density%data3d(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
                this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    CLASS DEFAULT
      CALL this%Error("gravity_sboxspectral::CalcPotential","unsupported state vector")
    END SELECT


#if defined(PARALLEL)
    this%mass2D(1:mesh%INUM,1:this%local_JNUM) = &
      reshape(this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), (/ mesh%INUM,mesh%JNUM/nprocs /))
#else
    this%mass2D(1:mesh%INUM,1:mesh%JNUM/nprocs) = &
      reshape(this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), (/ mesh%INUM,mesh%JNUM/nprocs /))
#endif

    !-------------- fourier transform (forward) of shifted density ----------!
! for performance checks
#ifdef _FTRACE
CALL ftrace_region_begin("forward_fft")
#endif

    CALL this%FFT_Forward(mesh,physics)

! performance checks
#ifdef _FTRACE
CALL ftrace_region_end("forward_fft")
#endif

    !------------- calculate potential in fourier domain --------------------!
    IF (mesh%dx .GT. mesh%dy) THEN
      k2max = 0.5*pi**2/(mesh%dx**2)
    ELSE
      k2max = 0.5*pi**2/(mesh%dy**2)
    END IF

    IF (mesh%shear_dir.EQ.2) THEN
!NEC$ IVDEP
      DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
        DO i = mesh%IMIN,mesh%IMAX/2+1
          k2 = (this%kx(i) + mesh%Q*mesh%OMEGA*this%ky(j)*delt)**2 + this%ky(j)**2
          ko = sqrt(k2)
          IF ((k2 .LT. k2max) .AND. (k2 .GT. 0)) THEN
            this%Fmass2D(i,j-this%local_joff) = -2.*pi*physics%Constants%GN* &
                                              this%Fmass2D(i,j-this%local_joff)/ko
          ELSE
            this%Fmass2D(i,j-this%local_joff) = 0.0
          END IF
        END DO
      END DO
    ELSE IF (mesh%shear_dir.EQ.1) THEN
!NEC$ IVDEP
      DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
        DO i = mesh%IMIN,mesh%IMAX/2+1
          k2 = this%kx(i)**2 + (this%ky(j)-mesh%Q*mesh%OMEGA*this%kx(i)*delt)**2
           ko = sqrt(k2)
          IF ((k2 .LT. k2max) .AND. (k2 .GT. 0)) THEN
            this%Fmass2D(i,j-this%local_joff) = -2.*pi*physics%Constants%GN* &
                                              this%Fmass2D(i,j-this%local_joff)/ko
          ELSE
            this%Fmass2D(i,j-this%local_joff) = 0.0
          END IF
        END DO
      END DO
    END IF

    !----------- fourier transform (backward) of shifted density -----------!
#ifdef _FTRACE
CALL ftrace_region_begin("backward_fft")
#endif

    CALL this%FFT_Backward(mesh,physics)

#ifdef _FTRACE
CALL ftrace_region_end("backward_fft")
#endif

    !------ calculate final potential with backshift and normalization ------!
!NEC$ IVDEP
    DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
      DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
        DO i = mesh%IMIN,mesh%IMAX
          this%phi(i,j,k) = this%mass2D(i,j-this%local_joff)/ &
               (mesh%JNUM*mesh%INUM)                        ! no norm. by FFTW !
        END DO
      END DO
    END DO
    !-------------------------- shift field back ----------------------------!
!!NEC$ IEXPAND
    CALL this%FieldShift(mesh,physics,-delt, &
            this%phi(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
            this%den_ip(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))

!NEC$ IVDEP
    DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
      DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
        DO i = mesh%IMIN,mesh%IMAX
          this%phi(i,j,k) = this%den_ip(i,j,k)
        END DO
      END DO
    END DO

    !----- copy values to boundaries in order to calculate acceleration -----!
    IF(mesh%shear_dir.EQ.2) THEN
!NEC$ SHORTLOOP
      DO j = 1,mesh%GJNUM
        ! southern northern (periodic)
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX)
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX)
      END DO
      joff2 = -this%shiftconst*delt
      jrem2 = joff2 - floor(joff2)
!NEC$ SHORTLOOP
      DO i = 1,mesh%GINUM
        ! western (periodic)
        this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMAX-i+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
        ! western integer shift
        this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)  = &
              cshift(this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX),floor(joff2))
        ! western residual shift
        this%phi(mesh%IMIN-i,mesh%JMAX+1,:) = this%phi(mesh%IMIN-i,mesh%JMIN,:)
        DO k=mesh%KMIN,mesh%KMAX
          DO j=mesh%JMIN,mesh%JMAX
            this%phi(mesh%IMIN-i,j,k) = &
                (1.0 - jrem2)*this%phi(mesh%IMIN-i,j,k) + jrem2*this%phi(mesh%IMIN-i,j+1,k)
          END DO
        END DO
      END DO
      joff2 = this%shiftconst*delt
      jrem2 = joff2 - floor(joff2)
!NEC$ SHORTLOOP
      DO i = 1,mesh%GINUM
        ! eastern (periodic)
        this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN+i-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
        ! eastern integer shift
        this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)  = &
              cshift(this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX),floor(joff2))
        ! eastern residual shift
        this%phi(mesh%IMAX+i,mesh%JMAX+1,:) = this%phi(mesh%IMAX+i,mesh%JMIN,:)
        DO k=mesh%KMIN,mesh%KMAX
          DO j=mesh%JMIN,mesh%JMAX
            this%phi(mesh%IMAX+i,j,k) = &
                (1.0 - jrem2)*this%phi(mesh%IMAX+i,j,k) + jrem2*this%phi(mesh%IMAX+i,j+1,k)
          END DO
        END DO
      END DO
! Only north-south direction has parallelization allowed
! Attention: The order of the copies plays a role. First the non-shifted direction needs to be
! copied, afterwards the shifted.
    ELSE IF(mesh%shear_dir.EQ.1) THEN
!NEC$ SHORTLOOP
      DO i = 1,mesh%GINUM
        ! western and eastern (always periodic)
        this%phi(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMAX-i+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
        this%phi(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN+i-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
      END DO
#ifdef PARALLEL
      IF(mesh%dims(2).GT.1) THEN
        mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX) = &
          this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-mesh%GJNUM+1:mesh%JMAX,mesh%KMIN:mesh%KMAX)
        CALL mpi_sendrecv_replace(&
          mpi_buf,&
          2*(mesh%IMAX-mesh%IMIN+1), &
          default_mpi_real, &
          mesh%neighbor(north), 53+north, &
          mesh%neighbor(south), mpi_any_tag, &
          mesh%comm_cart, status, ierror)
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JGMIN:mesh%JMIN-1,mesh%KMIN:mesh%KMAX) = &
          mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX)

        mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX) = &
          this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMIN+mesh%GJNUM-1,mesh%KMIN:mesh%KMAX)
        CALL mpi_sendrecv_replace(&
          mpi_buf,&
          2*(mesh%IMAX-mesh%IMIN+1), &
          default_mpi_real, &
          mesh%neighbor(south), 53+south, &
          mesh%neighbor(north), mpi_any_tag, &
          mesh%comm_cart, status, ierror)
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+1:mesh%JGMAX,mesh%KMIN:mesh%KMAX) = &
          mpi_buf(mesh%IMIN:mesh%IMAX,1:mesh%GJNUM,mesh%KMIN:mesh%KMAX)
      ELSE
!NEC$ SHORTLOOP
        DO j = 1,mesh%GJNUM
          ! southern northern (periodic in first step - further shift-treatment below)
          this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX)
          this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX)
        END DO
      END IF
#endif

#ifdef PARALLEL
      IF (mesh%mycoords(2).EQ.0) THEN
#endif
!NEC$ SHORTLOOP
      DO j = 1,mesh%GJNUM
        ! southern (shorn periodic) - residual and integer shift
#ifndef PARALLEL
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX)
#endif
        joff2 = this%shiftconst*delt
        jrem2 = joff2 - floor(joff2)
        !------- integral shift ---------------------------------------------!
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX)  = &
          cshift(this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX),floor(joff2))
        !------- residual shift ---------------------------------------------!
        this%phi(mesh%IMAX+1,mesh%JMIN-j,:) = this%phi(mesh%IMIN,mesh%JMIN-j,:)
        DO k = mesh%KMIN,mesh%KMAX
          DO i=mesh%IMIN,mesh%IMAX
            this%phi(i,mesh%JMIN-j,k) = &
              (1.0 - jrem2)*this%phi(i,mesh%JMIN-j,k) + jrem2*this%phi(i+1,mesh%JMIN-j,k)
          END DO
        END DO
      END DO
#ifdef PARALLEL
      END IF
#endif

#ifdef PARALLEL
      IF (mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
#endif
!NEC$ SHORTLOOP
      DO j = 1,mesh%GJNUM
#ifndef PARALLEL
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX) = this%phi(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX)
#endif
        ! northern (shorn periodic) - residual and integer shift
        joff2 = -this%shiftconst*delt
        jrem2 = joff2 - floor(joff2)
        !------- integral shift ---------------------------------------------!
        this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX)  = &
          cshift(this%phi(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX),floor(joff2))
        !------- residual shift ---------------------------------------------!
        this%phi(mesh%IMAX+1,mesh%JMAX+j,:) = this%phi(mesh%IMIN,mesh%JMAX+j,:)
        DO k = mesh%KMIN,mesh%KMAX
          DO i = mesh%IMIN,mesh%IMAX
            this%phi(i,mesh%JMAX+j,k) = &
              (1.0 - jrem2)*this%phi(i,mesh%JMAX+j,k) + jrem2*this%phi(i+1,mesh%JMAX+j,k)
          END DO
        END DO
      END DO
#ifdef PARALLEL
      END IF
#endif
    END IF
  END SUBROUTINE calcpotential


  SUBROUTINE fft_forward(this,Mesh,Physics)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT):: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    !------------------------------------------------------------------------!
    INTEGER :: i,j,k
#ifdef PARALLEL
    INTEGER :: nprocs,status(MPI_STATUS_SIZE)
    INTEGER :: ier
#endif
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    CALL mpi_comm_size(mpi_comm_world, nprocs, ier)
#endif

#ifdef _FTRACE
CALL ftrace_region_begin("forward FFT")
#endif

#if !defined(PARALLEL)
    CALL fftw_execute_dft_r2c(this%plan_r2c, this%mass2D, this%Fmass2D)
#else
    CALL fftw_mpi_execute_dft_r2c(this%plan_r2c,this%mass2D, this%Fmass2D)
#endif

    IF (this%calc_Fmass2D) THEN
!NEC$ IVDEP
      DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
        DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
          DO i = mesh%IMIN,mesh%IMAX/2+1
            this%Fmass2D_real(2*i-1,j,k) = REAL(REAL(this%Fmass2D(i,j-this%local_joff)))
            this%Fmass2D_real(2*i,j,k)   = REAL(AIMAG(this%Fmass2D(i,j-this%local_joff)))
          END DO
        END DO
      END DO
    END IF

#ifdef _FTRACE
CALL ftrace_region_end("foward FFT")
#endif

  END SUBROUTINE


  SUBROUTINE fft_backward(this,Mesh,Physics)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    INTEGER           :: nprocs,status(MPI_STATUS_SIZE)
    INTEGER           :: ier
#endif
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    CALL mpi_comm_size(mpi_comm_world, nprocs, ier)
#endif

#if defined(PARALLEL)
    CALL fftw_mpi_execute_dft_c2r(this%plan_c2r,this%Fmass2D, this%mass2D)
#else
    CALL fftw_execute_dft_c2r(this%plan_c2r, this%Fmass2D, this%mass2D)
#endif
  END SUBROUTINE
#endif


  SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),            INTENT(IN)    :: Mesh
    CLASS(physics_base),         INTENT(IN)    :: Physics
    CLASS(marray_compound),      INTENT(INOUT) :: pvar
    TYPE(marray_base),           INTENT(INOUT) :: bccsound,h_ext,height
    !------------------------------------------------------------------------!
    INTEGER           :: i,j,k
    REAL              :: cs2,p,q
    !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    ! pure self-gravitating shearing sheet with external point mass potential
!NEC$ COLLAPSE
    DO k=mesh%KGMIN,mesh%KGMAX
!NEC$ IVDEP
      DO j=mesh%JGMIN,mesh%JGMAX
!NEC$ IVDEP
        DO i=mesh%IGMIN,mesh%IGMAX
           cs2 = bccsound%data3d(i,j,k)*bccsound%data3d(i,j,k)
           p = -sqrttwopi*physics%Constants%GN*pvar%data4d(i,j,k,physics%DENSITY) &
                  /mesh%OMEGA**2.
           q = cs2/mesh%OMEGA**2.
           ! return the new disk height
           height%data3d(i,j,k) = p+sqrt(q+p*p)
        END DO
      END DO
    END DO
#endif
  END SUBROUTINE calcdiskheight_single

  SUBROUTINE infogravity(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(IN) :: this
    CLASS(mesh_base),            INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    CALL this%Info( "            FFT-Package:       FFTW " // &
#if defined(PARALLEL)
                                                    "- parallel mode")
#else
                                                    "- serial mode")
#endif
    CALL this%Info("            .. done initializing")
#endif
  END SUBROUTINE infogravity

#ifdef HAVE_FFTW
  SUBROUTINE fieldshift(this,Mesh,Physics,delt,field,mass2D)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)  :: Mesh
    CLASS(physics_base), INTENT(IN)  :: Physics
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
                         INTENT(IN)  :: field
#if defined(PARALLEL)
    REAL, DIMENSION(1:Mesh%INUM,1:this%local_JNUM), &
                         INTENT(OUT) :: mass2D
#else
    REAL, DIMENSION(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX), &
                         INTENT(OUT) :: mass2D
#endif
    !------------------------------------------------------------------------!
    INTEGER           :: i,j,k,local_joff
    REAL              :: delt
    !------------------------------------------------------------------------!
    local_joff = this%local_joff

    IF (mesh%shear_dir.EQ.1) THEN
!NEC$ IVDEP
      DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
        DO j = mesh%JMIN,mesh%JMAX
          this%joff(j)   = mesh%Q*mesh%OMEGA*mesh%bcenter(mesh%IMIN,j,k,2)* &
                      delt/mesh%dx
          this%jrem(j)  = this%joff(j) - floor(this%joff(j))
!NEC$ IVDEP
          DO i = mesh%IMIN,mesh%IMAX
            mass2d(i,j-local_joff) = &
             (1.0-this%jrem(j))*field(1+modulo(i-1+floor(this%joff(j)), &
              mesh%IMAX-mesh%IMIN+1),j,k) + this%jrem(j)*field(1+modulo(i+ &
              floor(this%joff(j)),mesh%IMAX-mesh%IMIN+1),j,k)
          END DO
        END DO
      END DO
    ELSE IF (mesh%shear_dir.EQ.2) THEN
!NEC$ IVDEP
      DO k = mesh%KMIN,mesh%KMAX
!NEC$ IVDEP
        DO i = mesh%IMIN,mesh%IMAX
          this%joff(i)   = -mesh%Q*mesh%OMEGA*mesh%bcenter(i,mesh%JMIN,k,1)* &
                      delt/mesh%dy
          this%jrem(i)  = this%joff(i) - floor(this%joff(i))
        END DO
        DO j = mesh%JMIN,mesh%JMAX
!NEC$ IVDEP
          DO i = mesh%IMIN,mesh%IMAX
            mass2d(i,j) = &
             (1.0-this%jrem(i))*field(i,1+modulo(j-1+floor(this%joff(i)), &
              mesh%JMAX-mesh%JMIN+1),k) + this%jrem(i)*field(i,1+modulo(j+ &
              floor(this%joff(i)),mesh%JMAX-mesh%JMIN+1),k)
          END DO
        END DO
      END DO
    END IF
  END SUBROUTINE fieldshift
#endif

  SUBROUTINE finalize(this)
   IMPLICIT NONE
   !------------------------------------------------------------------------!
   CLASS(gravity_sboxspectral), INTENT(INOUT) :: this
   !------------------------------------------------------------------------!
#ifdef HAVE_FFTW
    CALL fftw_destroy_plan(this%plan_r2c)
    CALL fftw_destroy_plan(this%plan_c2r)
#if defined(PARALLEL)
    CALL fftw_free(this%mass2D_pointer)
    CALL fftw_free(this%Fmass2D_pointer)
#else
    DEALLOCATE(this%mass2D,this%Fmass2D)
#endif
    ! Free memory
    DEALLOCATE(&
               this%phi, &
               this%Fmass2D_real, &
               this%kx, this%ky, &
               this%joff, &
               this%jrem, &
               this%den_ip &
               )
#endif
    END SUBROUTINE finalize

END MODULE gravity_sboxspectral_mod