text:
!===============================================================================
! Copyright 2015-2017 Intel Corporation All Rights Reserved.
!
! The source code,  information  and material  ("Material") contained  herein is
! owned by Intel Corporation or its  suppliers or licensors,  and  title to such
! Material remains with Intel  Corporation or its  suppliers or  licensors.  The
! Material  contains  proprietary  information  of  Intel or  its suppliers  and
! licensors.  The Material is protected by  worldwide copyright  laws and treaty
! provisions.  No part  of  the  Material   may  be  used,  copied,  reproduced,
! modified, published,  uploaded, posted, transmitted,  distributed or disclosed
! in any way without Intel's prior express written permission.  No license under
! any patent,  copyright or other  intellectual property rights  in the Material
! is granted to  or  conferred  upon  you,  either   expressly,  by implication,
! inducement,  estoppel  or  otherwise.  Any  license   under such  intellectual
! property rights must be express and approved by Intel in writing.
!
! Unless otherwise agreed by Intel in writing,  you may not remove or alter this
! notice or  any  other  notice   embedded  in  Materials  by  Intel  or Intel's
! suppliers or licensors in any way.
!===============================================================================

! Generated automatically.  DO NOT EDIT!

  include 'fftw3.f03'

  integer(C_INTPTR_T), parameter :: FFTW_MPI_DEFAULT_BLOCK = 0
  integer(C_INT), parameter :: FFTW_MPI_SCRAMBLED_IN = 134217728
  integer(C_INT), parameter :: FFTW_MPI_SCRAMBLED_OUT = 268435456
  integer(C_INT), parameter :: FFTW_MPI_TRANSPOSED_IN = 536870912
  integer(C_INT), parameter :: FFTW_MPI_TRANSPOSED_OUT = 1073741824

  type, bind(C) :: fftw_mpi_ddim
     integer(C_INTPTR_T) n, ib, ob
  end type fftw_mpi_ddim

  interface
    subroutine fftw_mpi_init() bind(C, name='fftw_mpi_init')
      import
    end subroutine fftw_mpi_init

    subroutine fftw_mpi_cleanup() bind(C, name='fftw_mpi_cleanup')
      import
    end subroutine fftw_mpi_cleanup

    integer(C_INTPTR_T) function fftw_mpi_local_size_many_transposed(rnk,n,howmany,block0,block1,comm,local_n0,local_0_start, &
                                                                     local_n1,local_1_start) &
                                 bind(C, name='fftw_mpi_local_size_many_transposed_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INTPTR_T), value :: block1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftw_mpi_local_size_many_transposed

    integer(C_INTPTR_T) function fftw_mpi_local_size_many(rnk,n,howmany,block0,comm,local_n0,local_0_start) &
                                 bind(C, name='fftw_mpi_local_size_many_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftw_mpi_local_size_many

    integer(C_INTPTR_T) function fftw_mpi_local_size_transposed(rnk,n,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftw_mpi_local_size_transposed_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftw_mpi_local_size_transposed

    integer(C_INTPTR_T) function fftw_mpi_local_size(rnk,n,comm,local_n0,local_0_start) bind(C, name='fftw_mpi_local_size_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftw_mpi_local_size

    integer(C_INTPTR_T) function fftw_mpi_local_size_many_1d(n0,howmany,comm,sign,flags,local_ni,local_i_start,local_no, &
                                                             local_o_start) bind(C, name='fftw_mpi_local_size_many_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: howmany
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
      integer(C_INTPTR_T), intent(out) :: local_ni
      integer(C_INTPTR_T), intent(out) :: local_i_start
      integer(C_INTPTR_T), intent(out) :: local_no
      integer(C_INTPTR_T), intent(out) :: local_o_start
    end function fftw_mpi_local_size_many_1d

    integer(C_INTPTR_T) function fftw_mpi_local_size_1d(n0,comm,sign,flags,local_ni,local_i_start,local_no,local_o_start) &
                                 bind(C, name='fftw_mpi_local_size_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
      integer(C_INTPTR_T), intent(out) :: local_ni
      integer(C_INTPTR_T), intent(out) :: local_i_start
      integer(C_INTPTR_T), intent(out) :: local_no
      integer(C_INTPTR_T), intent(out) :: local_o_start
    end function fftw_mpi_local_size_1d

    integer(C_INTPTR_T) function fftw_mpi_local_size_2d(n0,n1,comm,local_n0,local_0_start) &
                                 bind(C, name='fftw_mpi_local_size_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftw_mpi_local_size_2d

    integer(C_INTPTR_T) function fftw_mpi_local_size_2d_transposed(n0,n1,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftw_mpi_local_size_2d_transposed_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftw_mpi_local_size_2d_transposed

    integer(C_INTPTR_T) function fftw_mpi_local_size_3d(n0,n1,n2,comm,local_n0,local_0_start) &
                                 bind(C, name='fftw_mpi_local_size_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftw_mpi_local_size_3d

    integer(C_INTPTR_T) function fftw_mpi_local_size_3d_transposed(n0,n1,n2,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftw_mpi_local_size_3d_transposed_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftw_mpi_local_size_3d_transposed

    type(C_PTR) function fftw_mpi_plan_many_transpose(n0,n1,howmany,block0,block1,in,out,comm,flags) &
                         bind(C, name='fftw_mpi_plan_many_transpose_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INTPTR_T), value :: block1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_many_transpose

    type(C_PTR) function fftw_mpi_plan_transpose(n0,n1,in,out,comm,flags) bind(C, name='fftw_mpi_plan_transpose_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_transpose

    type(C_PTR) function fftw_mpi_plan_many_dft(rnk,n,howmany,block,tblock,in,out,comm,sign,flags) &
                         bind(C, name='fftw_mpi_plan_many_dft_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block
      integer(C_INTPTR_T), value :: tblock
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_many_dft

    type(C_PTR) function fftw_mpi_plan_dft(rnk,n,in,out,comm,sign,flags) bind(C, name='fftw_mpi_plan_dft_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft

    type(C_PTR) function fftw_mpi_plan_dft_1d(n0,in,out,comm,sign,flags) bind(C, name='fftw_mpi_plan_dft_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_1d

    type(C_PTR) function fftw_mpi_plan_dft_2d(n0,n1,in,out,comm,sign,flags) bind(C, name='fftw_mpi_plan_dft_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_2d

    type(C_PTR) function fftw_mpi_plan_dft_3d(n0,n1,n2,in,out,comm,sign,flags) bind(C, name='fftw_mpi_plan_dft_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_3d

    type(C_PTR) function fftw_mpi_plan_many_r2r(rnk,n,howmany,iblock,oblock,in,out,comm,kind,flags) &
                         bind(C, name='fftw_mpi_plan_many_r2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_many_r2r

    type(C_PTR) function fftw_mpi_plan_r2r(rnk,n,in,out,comm,kind,flags) bind(C, name='fftw_mpi_plan_r2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_r2r

    type(C_PTR) function fftw_mpi_plan_r2r_2d(n0,n1,in,out,comm,kind0,kind1,flags) bind(C, name='fftw_mpi_plan_r2r_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_r2r_2d

    type(C_PTR) function fftw_mpi_plan_r2r_3d(n0,n1,n2,in,out,comm,kind0,kind1,kind2,flags) bind(C, name='fftw_mpi_plan_r2r_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_FFTW_R2R_KIND), value :: kind2
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_r2r_3d

    type(C_PTR) function fftw_mpi_plan_many_dft_r2c(rnk,n,howmany,iblock,oblock,in,out,comm,flags) &
                         bind(C, name='fftw_mpi_plan_many_dft_r2c_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_many_dft_r2c

    type(C_PTR) function fftw_mpi_plan_dft_r2c(rnk,n,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_r2c_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_r2c

    type(C_PTR) function fftw_mpi_plan_dft_r2c_2d(n0,n1,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_r2c_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_r2c_2d

    type(C_PTR) function fftw_mpi_plan_dft_r2c_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_r2c_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_r2c_3d

    type(C_PTR) function fftw_mpi_plan_many_dft_c2r(rnk,n,howmany,iblock,oblock,in,out,comm,flags) &
                         bind(C, name='fftw_mpi_plan_many_dft_c2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_many_dft_c2r

    type(C_PTR) function fftw_mpi_plan_dft_c2r(rnk,n,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_c2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_c2r

    type(C_PTR) function fftw_mpi_plan_dft_c2r_2d(n0,n1,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_c2r_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_c2r_2d

    type(C_PTR) function fftw_mpi_plan_dft_c2r_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_c2r_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_c2r_3d

    subroutine fftw_mpi_gather_wisdom(comm_) bind(C, name='fftw_mpi_gather_wisdom_f03')
      import
      integer(C_INT32_T), value :: comm_
    end subroutine fftw_mpi_gather_wisdom

    subroutine fftw_mpi_broadcast_wisdom(comm_) bind(C, name='fftw_mpi_broadcast_wisdom_f03')
      import
      integer(C_INT32_T), value :: comm_
    end subroutine fftw_mpi_broadcast_wisdom

    subroutine fftw_mpi_execute_dft(p,in,out) bind(C, name='fftw_mpi_execute_dft')
      import
      type(C_PTR), value :: p
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftw_mpi_execute_dft

    subroutine fftw_mpi_execute_dft_r2c(p,in,out) bind(C, name='fftw_mpi_execute_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftw_mpi_execute_dft_r2c

    subroutine fftw_mpi_execute_dft_c2r(p,in,out) bind(C, name='fftw_mpi_execute_dft_c2r')
      import
      type(C_PTR), value :: p
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
    end subroutine fftw_mpi_execute_dft_c2r

    subroutine fftw_mpi_execute_r2r(p,in,out) bind(C, name='fftw_mpi_execute_r2r')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
    end subroutine fftw_mpi_execute_r2r

  end interface

  type, bind(C) :: fftwf_mpi_ddim
     integer(C_INTPTR_T) n, ib, ob
  end type fftwf_mpi_ddim

  interface
    subroutine fftwf_mpi_init() bind(C, name='fftwf_mpi_init')
      import
    end subroutine fftwf_mpi_init

    subroutine fftwf_mpi_cleanup() bind(C, name='fftwf_mpi_cleanup')
      import
    end subroutine fftwf_mpi_cleanup

    integer(C_INTPTR_T) function fftwf_mpi_local_size_many_transposed(rnk,n,howmany,block0,block1,comm,local_n0,local_0_start, &
                                                                      local_n1,local_1_start) &
                                 bind(C, name='fftwf_mpi_local_size_many_transposed_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INTPTR_T), value :: block1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftwf_mpi_local_size_many_transposed

    integer(C_INTPTR_T) function fftwf_mpi_local_size_many(rnk,n,howmany,block0,comm,local_n0,local_0_start) &
                                 bind(C, name='fftwf_mpi_local_size_many_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftwf_mpi_local_size_many

    integer(C_INTPTR_T) function fftwf_mpi_local_size_transposed(rnk,n,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftwf_mpi_local_size_transposed_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftwf_mpi_local_size_transposed

    integer(C_INTPTR_T) function fftwf_mpi_local_size(rnk,n,comm,local_n0,local_0_start) bind(C, name='fftwf_mpi_local_size_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftwf_mpi_local_size

    integer(C_INTPTR_T) function fftwf_mpi_local_size_many_1d(n0,howmany,comm,sign,flags,local_ni,local_i_start,local_no, &
                                                              local_o_start) bind(C, name='fftwf_mpi_local_size_many_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: howmany
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
      integer(C_INTPTR_T), intent(out) :: local_ni
      integer(C_INTPTR_T), intent(out) :: local_i_start
      integer(C_INTPTR_T), intent(out) :: local_no
      integer(C_INTPTR_T), intent(out) :: local_o_start
    end function fftwf_mpi_local_size_many_1d

    integer(C_INTPTR_T) function fftwf_mpi_local_size_1d(n0,comm,sign,flags,local_ni,local_i_start,local_no,local_o_start) &
                                 bind(C, name='fftwf_mpi_local_size_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
      integer(C_INTPTR_T), intent(out) :: local_ni
      integer(C_INTPTR_T), intent(out) :: local_i_start
      integer(C_INTPTR_T), intent(out) :: local_no
      integer(C_INTPTR_T), intent(out) :: local_o_start
    end function fftwf_mpi_local_size_1d

    integer(C_INTPTR_T) function fftwf_mpi_local_size_2d(n0,n1,comm,local_n0,local_0_start) &
                                 bind(C, name='fftwf_mpi_local_size_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftwf_mpi_local_size_2d

    integer(C_INTPTR_T) function fftwf_mpi_local_size_2d_transposed(n0,n1,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftwf_mpi_local_size_2d_transposed_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftwf_mpi_local_size_2d_transposed

    integer(C_INTPTR_T) function fftwf_mpi_local_size_3d(n0,n1,n2,comm,local_n0,local_0_start) &
                                 bind(C, name='fftwf_mpi_local_size_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
    end function fftwf_mpi_local_size_3d

    integer(C_INTPTR_T) function fftwf_mpi_local_size_3d_transposed(n0,n1,n2,comm,local_n0,local_0_start,local_n1,local_1_start) &
                                 bind(C, name='fftwf_mpi_local_size_3d_transposed_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      integer(C_INT32_T), value :: comm
      integer(C_INTPTR_T), intent(out) :: local_n0
      integer(C_INTPTR_T), intent(out) :: local_0_start
      integer(C_INTPTR_T), intent(out) :: local_n1
      integer(C_INTPTR_T), intent(out) :: local_1_start
    end function fftwf_mpi_local_size_3d_transposed

    type(C_PTR) function fftwf_mpi_plan_many_transpose(n0,n1,howmany,block0,block1,in,out,comm,flags) &
                         bind(C, name='fftwf_mpi_plan_many_transpose_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block0
      integer(C_INTPTR_T), value :: block1
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_many_transpose

    type(C_PTR) function fftwf_mpi_plan_transpose(n0,n1,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_transpose_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_transpose

    type(C_PTR) function fftwf_mpi_plan_many_dft(rnk,n,howmany,block,tblock,in,out,comm,sign,flags) &
                         bind(C, name='fftwf_mpi_plan_many_dft_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: block
      integer(C_INTPTR_T), value :: tblock
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_many_dft

    type(C_PTR) function fftwf_mpi_plan_dft(rnk,n,in,out,comm,sign,flags) bind(C, name='fftwf_mpi_plan_dft_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft

    type(C_PTR) function fftwf_mpi_plan_dft_1d(n0,in,out,comm,sign,flags) bind(C, name='fftwf_mpi_plan_dft_1d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_1d

    type(C_PTR) function fftwf_mpi_plan_dft_2d(n0,n1,in,out,comm,sign,flags) bind(C, name='fftwf_mpi_plan_dft_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_2d

    type(C_PTR) function fftwf_mpi_plan_dft_3d(n0,n1,n2,in,out,comm,sign,flags) bind(C, name='fftwf_mpi_plan_dft_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_3d

    type(C_PTR) function fftwf_mpi_plan_many_r2r(rnk,n,howmany,iblock,oblock,in,out,comm,kind,flags) &
                         bind(C, name='fftwf_mpi_plan_many_r2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_many_r2r

    type(C_PTR) function fftwf_mpi_plan_r2r(rnk,n,in,out,comm,kind,flags) bind(C, name='fftwf_mpi_plan_r2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_r2r

    type(C_PTR) function fftwf_mpi_plan_r2r_2d(n0,n1,in,out,comm,kind0,kind1,flags) bind(C, name='fftwf_mpi_plan_r2r_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_r2r_2d

    type(C_PTR) function fftwf_mpi_plan_r2r_3d(n0,n1,n2,in,out,comm,kind0,kind1,kind2,flags) &
                         bind(C, name='fftwf_mpi_plan_r2r_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_FFTW_R2R_KIND), value :: kind2
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_r2r_3d

    type(C_PTR) function fftwf_mpi_plan_many_dft_r2c(rnk,n,howmany,iblock,oblock,in,out,comm,flags) &
                         bind(C, name='fftwf_mpi_plan_many_dft_r2c_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_many_dft_r2c

    type(C_PTR) function fftwf_mpi_plan_dft_r2c(rnk,n,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_r2c_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_r2c

    type(C_PTR) function fftwf_mpi_plan_dft_r2c_2d(n0,n1,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_r2c_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_r2c_2d

    type(C_PTR) function fftwf_mpi_plan_dft_r2c_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_r2c_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_r2c_3d

    type(C_PTR) function fftwf_mpi_plan_many_dft_c2r(rnk,n,howmany,iblock,oblock,in,out,comm,flags) &
                         bind(C, name='fftwf_mpi_plan_many_dft_c2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      integer(C_INTPTR_T), value :: howmany
      integer(C_INTPTR_T), value :: iblock
      integer(C_INTPTR_T), value :: oblock
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_many_dft_c2r

    type(C_PTR) function fftwf_mpi_plan_dft_c2r(rnk,n,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_c2r_f03')
      import
      integer(C_INT), value :: rnk
      integer(C_INTPTR_T), dimension(*), intent(in) :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_c2r

    type(C_PTR) function fftwf_mpi_plan_dft_c2r_2d(n0,n1,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_c2r_2d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_c2r_2d

    type(C_PTR) function fftwf_mpi_plan_dft_c2r_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftwf_mpi_plan_dft_c2r_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT32_T), value :: comm
      integer(C_INT), value :: flags
    end function fftwf_mpi_plan_dft_c2r_3d

    subroutine fftwf_mpi_gather_wisdom(comm_) bind(C, name='fftwf_mpi_gather_wisdom_f03')
      import
      integer(C_INT32_T), value :: comm_
    end subroutine fftwf_mpi_gather_wisdom

    subroutine fftwf_mpi_broadcast_wisdom(comm_) bind(C, name='fftwf_mpi_broadcast_wisdom_f03')
      import
      integer(C_INT32_T), value :: comm_
    end subroutine fftwf_mpi_broadcast_wisdom

    subroutine fftwf_mpi_execute_dft(p,in,out) bind(C, name='fftwf_mpi_execute_dft')
      import
      type(C_PTR), value :: p
      complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftwf_mpi_execute_dft

    subroutine fftwf_mpi_execute_dft_r2c(p,in,out) bind(C, name='fftwf_mpi_execute_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftwf_mpi_execute_dft_r2c

    subroutine fftwf_mpi_execute_dft_c2r(p,in,out) bind(C, name='fftwf_mpi_execute_dft_c2r')
      import
      type(C_PTR), value :: p
      complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
    end subroutine fftwf_mpi_execute_dft_c2r

    subroutine fftwf_mpi_execute_r2r(p,in,out) bind(C, name='fftwf_mpi_execute_r2r')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
    end subroutine fftwf_mpi_execute_r2r

  end interface

max_stars_repo_licenses: Intel MIT

max_issues_repo_licenses: Intel MIT

max_forks_repo_licenses: Intel MIT


text:
MODULE constrain_SALSA
  USE mo_progn_state, ONLY : a_naerop, a_naerot, a_ncloudp, a_ncloudt, a_nprecpp, a_nprecpt,   &
                             a_maerop, a_maerot, a_mcloudp, a_mcloudt, a_mprecpp, a_mprecpt,   &
                             a_nicep,  a_nicet,  a_micep,   a_micet, a_gaerop, a_indefp
  USE mo_diag_state, ONLY : a_rc, a_srp, a_snrp, a_rh, a_temp, a_ri, a_riri, a_rhi
  USE mo_aux_state, ONLY : aetot
  USE mo_submctl, ONLY : spec, nlim, prlim, ice_theta_dist
  USE grid, ONLY : level
  USE mo_structured_datatypes
  IMPLICIT NONE

  CONTAINS

   SUBROUTINE tend_constrain2
     USE grid, ONLY : dtlt
     USE mo_field_state, ONLY : SALSA_tracers_4d
     IMPLICIT NONE

      TYPE(FloatArray4d), POINTER :: varp => NULL(),vart => NULL()
      INTEGER :: nv

      DO nv = 1,SALSA_tracers_4d%count
         CALL SALSA_tracers_4d%getData(1,varp,nv)
         CALL SALSA_tracers_4d%getData(2,vart,nv)
         vart%d(:,:,:,:) = MAX( (-1.+1.e-8)*varp%d(:,:,:,:)/dtlt, vart%d(:,:,:,:)  )
      END DO
      varp => NULL(); vart => NULL()

      IF (ice_theta_dist .AND. level > 4)   &
           a_indefp%d(:,:,:,:) = MIN( MAX( a_indefp%d(:,:,:,:),0. ),1. )

   END SUBROUTINE tend_constrain2

   !
   ! ---------------------------------------------------------------------
   ! SALSA_diagnostics: Update properties for the current timestep:
   !                    E.g. if enough water has evaporated from droplets,
   !                    deplete the cloud droplet bins and move CCN material
   !                    back to the aerosol regime.
   !                    In addition, update the diagnostic scalars for total grid-cell
   !                    liquid water contents.
   !
   ! Juha Tonttila, FMI, 2014
   ! Tomi Raatikainen, FMI, 2016

   SUBROUTINE SALSA_diagnostics(onlyDiag)
     USE grid, ONLY : nxp,nyp,nzp,    &
                      level
     USE mo_derived_procedures, ONLY : binMixrat ! Maybe at some point one could use the actual variables instead of this subroutine directly?
     USE mo_submctl, ONLY : nbins,ncld,nprc,nice,ica,fca,icb,fcb,ira,fra,              &
                            in1a,fn2b,                        &
                            nice,iia,fia,                    &
                            surfw0, rg, pi, &
                            lscndgas, pi6, avog
     USE util, ONLY : getMassIndex
     IMPLICIT NONE

     LOGICAL, INTENT(in), OPTIONAL :: onlyDiag ! If true, only update diagnostic concentrations and don't do anything else. Default value is FALSE.

     REAL, PARAMETER :: massTH = 1.e-25! Minimum mass threshold; corresponds to about a 1 nm particle == does not make sense

     INTEGER :: i,j,k,bc,ba,s,sc,sa,str,end,nc,nspec
     INTEGER :: mi,mi2 !"mass index" variables

     REAL :: zvol
     REAL :: zdh2o
     REAL :: ns, zbb, zaa ! Number of moles, Raoult effect, Kelvin effect; For calculating the critical radius
     REAL :: cdcld,cdprc, cdice! Critical diameter for cloud droplets and precipitation
     REAL :: zdrms, zwams
     REAL :: mice_tot   ! Variable for total ice mass
     REAL :: vice_tot   ! Variable for total ice vol

     LOGICAL :: l_onlyDiag

     nspec = spec%getNSpec(type="wet") ! Aerosol species + water. For rime add +1

     IF (PRESENT(onlyDiag)) THEN
        l_onlyDiag = onlyDiag
     ELSE
        l_onlyDiag = .FALSE.
     END IF

     ! Remove negative values
     a_naerop%d = MAX(0.,a_naerop%d)
     a_ncloudp%d = MAX(0.,a_ncloudp%d)
     a_nprecpp%d = MAX(0.,a_nprecpp%d)
     a_maerop%d = MAX(0.,a_maerop%d)
     a_mcloudp%d = MAX(0.,a_mcloudp%d)
     a_mprecpp%d = MAX(0.,a_mprecpp%d)

     IF (level == 5) THEN
        a_nicep%d = MAX(0.,a_nicep%d)
        a_micep%d = MAX(0.,a_micep%d)
     END IF

     IF ( .NOT. l_onlyDiag) THEN

        ! Remove particles that have number but not mass
        DO j = 3,nyp-2
           DO i = 3,nxp-2
              DO k = 1,nzp
                 ! Aerosols
                 DO bc = 1,nbins
                    mi = getMassIndex(nbins,bc,nspec)
                    IF (a_naerop%d(k,i,j,bc) > 0. .AND. SUM(a_maerop%d(k,i,j,bc:mi:nbins)) <= 0.) THEN
                       a_naerop%d(k,i,j,bc) = 0.
                       a_maerop%d(k,i,j,bc:mi:nbins) = 0.
                    END IF
                 END DO

                 ! Clouds
                 DO bc = 1,ncld
                    mi = getMassIndex(ncld,bc,nspec)
                    IF (a_ncloudp%d(k,i,j,bc) > 0. .AND. SUM(a_mcloudp%d(k,i,j,bc:mi:ncld)) <= 0.) THEN
                       a_ncloudp%d(k,i,j,bc) = 0.
                       a_mcloudp%d(k,i,j,bc:mi:ncld) = 0.
                    END IF
                 END DO ! ncld

                 ! Precipitation
                 DO bc = 1,nprc
                    mi = getMassIndex(nprc,bc,nspec)
                    IF (a_nprecpp%d(k,i,j,bc) > 0. .AND. a_mprecpp%d(k,i,j,mi) <= 0.) THEN
                       a_nprecpp%d(k,i,j,bc) = 0.
                       a_mprecpp%d(k,i,j,bc:mi:nprc) = 0.
                    END IF
                 END DO ! nprc

                 ! Ice
                 IF (level<5) CYCLE
                 DO bc = 1,nice
                    mi = getMassIndex(nice,bc,nspec)     ! Pristine
                    mi2 = getMassIndex(nice,bc,nspec+1)  ! Rime
                    mice_tot = a_micep%d(k,i,j,mi) + a_micep%d(k,i,j,mi2)
                    IF (a_nicep%d(k,i,j,bc) > 0. .AND. mice_tot < massTH) THEN
                       a_nicep%d(k,i,j,bc) = 0.
                       a_micep%d(k,i,j,bc:mi:nice) = 0.
                       a_micep%d(k,i,j,bc:mi2:nice) = 0.
                    END IF
                 END DO ! nice

              END DO !k
           END DO !i
        END DO !j

        ! Check evaporation of particles
        ! --------------------------------------
        ! Loop over cloud droplet bins
        DO bc = ica%cur,fcb%cur
           DO j = 3,nyp-2
              DO i = 3,nxp-2
                 DO k = 1,nzp

                    mi = getMassIndex(ncld,bc,nspec)    ! all compounds
                    mi2 = getMassIndex(ncld,bc,nspec-1) ! Dry aerosol part
                    IF ( a_ncloudp%d(k,i,j,bc) > nlim .AND. a_rh%d(k,i,j)<0.999 .AND.   &
                         a_mcloudp%d(k,i,j,mi) < 1.e-5  ) THEN

                       ! Critical diameter
                       ns = SUM( spec%diss(1:nspec-1)*a_mcloudp%d(k,i,j,bc:mi2:ncld)/spec%MM(1:nspec-1) ) / &
                            a_ncloudp%d(k,i,j,bc)
                       zbb = 6.*spec%mwa*ns/(pi*spec%rhowa)
                       zaa = 4.*spec%mwa*surfw0/(rg*spec%rhowa*a_temp%d(k,i,j))
                       cdcld = SQRT(3.*zbb/zaa)

                       ! Wet diameter
                       zvol = SUM( a_mcloudp%d(k,i,j,bc:mi:ncld)/spec%rholiq(1:nspec) )/a_ncloudp%d(k,i,j,bc)
                       zdh2o = (zvol/pi6)**(1./3.)

                       ! Lose the droplets if small
                       IF ( zdh2o < MIN(0.7*cdcld,1.3e-5) .OR.     &
                            a_mcloudp%d(k,i,j,mi) < massTH*a_ncloudp%d(k,i,j,bc) )  THEN

                          IF (bc<=fca%cur) THEN
                             ba = ica%par + (bc-ica%cur) ! Index for parallel aerosol bin
                          ELSE
                             ba = icb%par + (bc-icb%cur) ! Index for parallel aerosol bin
                          END IF
                          ! Move the number of particles from cloud to aerosol bins
                          a_naerop%d(k,i,j,ba) = a_naerop%d(k,i,j,ba) + a_ncloudp%d(k,i,j,bc)
                          a_ncloudp%d(k,i,j,bc) = 0.

                          ! Move ccn material back to aerosol regime (including water)
                          DO s = 1,nspec
                             sc = getMassIndex(ncld,bc,s)
                             sa = getMassIndex(nbins,ba,s)
                             a_maerop%d(k,i,j,sa) = a_maerop%d(k,i,j,sa) + a_mcloudp%d(k,i,j,sc)
                             a_mcloudp%d(k,i,j,sc) = 0.
                          END DO

                       END IF ! critical diameter
                    END IF  ! nlim

                 END DO
              END DO
           END DO
        END DO ! bc

        ! Loop over precipitation bins
        DO bc = 1,nprc
           DO j = 3,nyp-2
              DO i = 3,nxp-2
                 DO k = 1,nzp

                    mi = getMassIndex(nprc,bc,nspec)    ! all
                    mi2 = getMassIndex(nprc,bc,nspec-1) ! Dry aerosol part
                    IF ( a_nprecpp%d(k,i,j,bc) > prlim .AND. a_rh%d(k,i,j)<0.999 .AND.  &
                         a_mprecpp%d(k,i,j,mi) < 1.e-6 ) THEN

                       ! Critical diameter
                       ns = SUM( spec%diss(1:nspec-1)*a_mprecpp%d(k,i,j,bc:mi2:nprc)/spec%MM(1:nspec-1) ) / &
                            a_nprecpp%d(k,i,j,bc)
                       zbb = 6.*spec%mwa*ns/(pi*spec%rhowa)
                       zaa = 4.*spec%mwa*surfw0/(rg*spec%rhowa*a_temp%d(k,i,j))
                       cdprc = SQRT(3.*zbb/zaa)

                       ! Wet diameter
                       zvol = SUM( a_mprecpp%d(k,i,j,bc:mi:nprc)/spec%rholiq(1:nspec) )/a_nprecpp%d(k,i,j,bc)
                       zdh2o = (zvol/pi6)**(1./3.)

                       ! Lose the droplets if small
                       IF ( zdh2o < MIN(0.1*cdprc,1.3e-5) .OR.   &
                            a_mprecpp%d(k,i,j,mi) < massTH*a_nprecpp%d(k,i,j,bc) ) THEN

                          ba = findDry4Wet(a_nprecpp,a_mprecpp,nprc,bc,k,i,j)

                          ! Move the number of particles from cloud to aerosol bins
                          a_naerop%d(k,i,j,ba) = a_naerop%d(k,i,j,ba) + a_nprecpp%d(k,i,j,bc)
                          a_nprecpp%d(k,i,j,bc) = 0.

                          ! Move ccn material back to aerosol regime (including water)
                          DO s = 1,nspec
                             sc = getMassIndex(nprc,bc,s)
                             sa = getMassIndex(nbins,ba,s)
                             a_maerop%d(k,i,j,sa) = a_maerop%d(k,i,j,sa) + a_mprecpp%d(k,i,j,sc)
                             a_mprecpp%d(k,i,j,sc) = 0.
                          END DO

                       END IF ! Critical diameter

                    END IF ! prlim

                 END DO
              END DO
           END DO
        END DO ! bc

        ! Loop over ice bins
        IF (level == 5) THEN
           DO bc = 1,nice
              DO j = 3,nyp-2
                 DO i = 3,nxp-2
                    DO k = 1,nzp

                       mi = getMassIndex(nice,bc,nspec)     ! Pristine ice
                       mi2 = getMassIndex(nice,bc,nspec+1)  ! rimed ice
                       mice_tot = a_micep%d(k,i,j,mi) + a_micep%d(k,i,j,mi2)
                       vice_tot = SUM(a_micep%d(k,i,j,bc:mi2:nice)/spec%rhoice(1:nspec+1))

                       IF ( a_nicep%d(k,i,j,bc) > prlim .AND. a_rhi%d(k,i,j)<0.999 .AND.   &
                            mice_tot < 1.e-15 ) THEN

                          ! Ice and snow don't have a critical size, but lose particles when water content becomes low enough or size is < 2e-6m

                          ! Diameter
                          cdice = vice_tot / a_nicep%d(k,i,j,bc)
                          cdice = (cdice/pi6)**(1./3.)

                          ! Lose ice when dry to total mass ratio is more than 0.5
                          CALL binMixrat("ice","dry",bc,i,j,k,zdrms)
                          CALL binMixrat("ice","wet",bc,i,j,k,zwams)
                          zvol = zdrms/zwams

                          IF ( zvol>0.5 .OR. cdice < 2.e-6 ) THEN

                             ! Find the aerosol bin corresponding to the composition of the IN the current bin
                             ba = findDry4Wet(a_nicep,a_micep,nice,bc,k,i,j)

                             ! Move the number of particles from ice to aerosol bins
                             a_naerop%d(k,i,j,ba) = a_naerop%d(k,i,j,ba) + a_nicep%d(k,i,j,bc)
                             a_nicep%d(k,i,j,bc) = 0.

                             ! Move mass material back to aerosol regime (including water)
                             DO s = 1,nspec
                                sc = getMassIndex(nice,bc,s)
                                sa = getMassIndex(nbins,ba,s)
                                a_maerop%d(k,i,j,sa) = a_maerop%d(k,i,j,sa) + a_micep%d(k,i,j,sc)
                                a_micep%d(k,i,j,sc) = 0.
                             END DO
                             ! Rimed ice
                             sa = getMassIndex(nbins,ba,nspec)
                             sc = getMassIndex(nice,bc,nspec+1)
                             a_maerop%d(k,i,j,sa) = a_maerop%d(k,i,j,sa) + a_micep%d(k,i,j,sc)
                             a_micep%d(k,i,j,sc) = 0.
                          END IF

                       END IF  ! prlim

                    END DO
                 END DO
              END DO
           END DO ! bc
        END IF ! level = 5

        ! Loop over aerosol bins
        DO ba = 1,nbins
           DO j = 3,nyp-2
              DO i = 3,nxp-2
                 DO k = 1,nzp

                    mi = getMassIndex(nbins,ba,nspec)    ! all
                    mi2 = getMassIndex(nbins,ba,nspec-1) ! dry part
                    IF (a_naerop%d(k,i,j,ba) > nlim) THEN
                       zvol = SUM( a_maerop%d(k,i,j,ba:mi2:nbins)/spec%rholiq(1:nspec-1) )/ &
                            a_naerop%d(k,i,j,ba) ! Dry volume

                       ! Particles smaller then 0.1 nm diameter are set to zero
                       IF ( zvol < pi6*1.e-10**3 ) THEN
                          ! Volatile species to the gas phase
                          IF (spec%isUsed('SO4') .AND. lscndgas) THEN
                             nc = spec%getIndex('SO4')
                             s = getMassIndex(nbins,ba,nc)
                             a_gaerop%d(k,i,j,1) = a_gaerop%d(k,i,j,1) + a_maerop%d(k,i,j,s) / spec%msu * avog
                          END IF
                          IF (spec%isUsed('OC') .AND. lscndgas) THEN
                             nc = spec%getIndex('OC')
                             s = getMassIndex(nbins,ba,nc)
                             a_gaerop%d(k,i,j,5) = a_gaerop%d(k,i,j,5) + a_maerop%d(k,i,j,s) / spec%moc * avog
                          END IF
                          IF (spec%isUsed('NO') .AND. lscndgas) THEN
                             nc = spec%getIndex('NO')
                             s = getMassIndex(nbins,ba,nc)
                             a_gaerop%d(k,i,j,2) = a_gaerop%d(k,i,j,2) + a_maerop%d(k,i,j,s) / spec%mno * avog
                          END IF
                          IF (spec%isUsed('NH') .AND. lscndgas) THEN
                             nc = spec%getIndex('NH')
                             s = getMassIndex(nbins,ba,nc)
                             a_gaerop%d(k,i,j,3) = a_gaerop%d(k,i,j,3) + a_maerop%d(k,i,j,s) / spec%mnh * avog
                          END IF

                          ! Mass and number to zero (insolube species and water are lost)
                          a_maerop%d(k,i,j,ba:mi:nbins) = 0.
                          a_naerop%d(k,i,j,ba) = 0.
                       END IF
                    END IF

                 END DO
              END DO
           END DO
        END DO

     END IF ! onlyDiag

    !!!!!!!!!!!!!!!!!!!!!!!
    ! Update diagnostic tracers
    !!!!!!!!!!!!!!!!!!!!!!!

    ! Liquid water content
    nc = spec%getIndex('H2O')
    ! Aerosols, regimes a and b
    str = getMassIndex(nbins,in1a,nc)
    end = getMassIndex(nbins,fn2b,nc)
    a_rc%d(:,:,:) = SUM(a_maerop%d(:,:,:,str:end),DIM=4)
    ! Clouds, regime a and b
    str = getMassIndex(ncld,ica%cur,nc)
    end = getMassIndex(ncld,fcb%cur,nc)
    a_rc%d(:,:,:) = a_rc%d(:,:,:) + SUM(a_mcloudp%d(:,:,:,str:end),DIM=4)
    ! Precipitation
    str = getMassIndex(nprc,ira,nc)
    end = getMassIndex(nprc,fra,nc)
    a_srp%d(:,:,:) = SUM(a_mprecpp%d(:,:,:,str:end),DIM=4)
    a_snrp%d(:,:,:) = SUM(a_nprecpp%d(:,:,:,ira:fra),DIM=4)
    ! ice
    IF (level == 5) THEN
       str = getMassIndex(nice,iia,nc)
       end = getMassIndex(nice,fia,nc)
       a_ri%d(:,:,:) = SUM(a_micep%d(:,:,:,str:end),DIM=4)
       ! Rimed ice
       nc = spec%getIndex("rime")
       str = getMassIndex(nice,iia,nc)
       end = getMassIndex(nice,fia,nc)
       a_riri%d(:,:,:) = SUM(a_micep%d(:,:,:,str:end),DIM=4)
    END IF

  END SUBROUTINE SALSA_diagnostics

  !
  ! ----------------------------------------------------------------------
  !

  INTEGER FUNCTION findDry4Wet(nevap,mevap,nb,ib,iz,ix,iy)
    USE mo_submctl, ONLY : in2a,fn2a,nbins,pi6
    USE util, ONLY :  calc_correlation, getMassIndex
    IMPLICIT NONE

    INTEGER, INTENT(in)            :: nb  ! Number of bins and number of species in the evaporating particle class
    INTEGER, INTENT(in)            :: ib,iz,ix,iy  ! Current bin and grid indices
    TYPE(FloatArray4d), INTENT(in) :: mevap
    TYPE(FloatArray4d), INTENT(in) :: nevap

    REAL :: cd        ! dry particle diameter
    INTEGER :: mi,mi2 ! Mass indices
    INTEGER :: ba,bb  ! Corresponding aerosol indices for regimes a and b
    REAL :: ra, rb ! Correlation coefficients for a and b aerosol bins

    INTEGER :: ndry

    ! This function finds a suitable aerosol bin for evaporating
    ! precipitation or ice bins

    ndry = spec%getNSpec(type="dry")

    mi = getMassIndex(nb,ib,ndry)  ! Final dry mass index

    ! 1) Find the closest matching bin based on dry particle diameter (a and b bins)
    !    -- note that spec%rholiq for all evaporating categories is ok here,
    !    -- since the aerosol densities are always the same
    cd = SUM( mevap%d(iz,ix,iy,ib:mi:nb)/spec%rholiq(1:ndry) ) / &
              nevap%d(iz,ix,iy,ib)
    cd = (cd/pi6)**(1./3.) ! Dry diameter

    ba=in2a !Ignore 1a and note that "aerobins" contains the lower limit of bin dry diameter
    ba = MIN( MAX(in2a + COUNT( cd > aetot%d(ba:fn2a) )-1, in2a), fn2a )

    ! Corresponding b bin is ba+(fn2a-fn1a)
    bb = fn2a + (ba - in2a)

    ! 2) Select a or b bin
    IF (a_naerop%d(iz,ix,iy,bb) <= nlim) THEN
       ! Empty b bin so select a, i.e. do nothing here
       findDry4Wet = ba
    ELSE IF (a_naerop%d(iz,ix,iy,ba) <= nlim) THEN
       ! Empty a bin so select b
       findDry4Wet = bb
    ELSE
       mi = getMassIndex(nbins,ba,ndry) ! Index of the last dry species for current bin in aerosol regime a
       mi2 = getMassIndex(nb,ib,ndry)   ! The same for the evaporating particle

       ! Both are present - find bin based on compositional similarity
       ra = calc_correlation(a_maerop%d(iz,ix,iy,ba:mi:nbins),  &
                             mevap%d(iz,ix,iy,ib:mi2:nb),   &
                             ndry                       )

       mi = getMassIndex(nbins,bb,ndry) ! Index of the last dry species for current bin in aerosol regime b
       rb = calc_correlation(a_maerop%d(iz,ix,iy,bb:mi:nbins),  &
                             mevap%d(iz,ix,iy,ib:mi2:nb),   &
                             ndry                       )

       findDry4Wet = ba
       IF (ra<rb) findDry4Wet = bb
    END IF

  END FUNCTION findDry4Wet

END MODULE constrain_SALSA

max_stars_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_issues_repo_name: UCLALES-SALSA/UCLALES-SALSA

max_forks_repo_name: UCLALES-SALSA/UCLALES-SALSA


text: Many people in Davis play with flying discs, most commonly known by the brand specific frisbee.  From casual games of toss at local parks, possibly with dogs a canine companion, to serious sports like ultimate Frisbee and Frisbee golf, there are opportunities for all levels and interests of players.  If youre a students student at UCD, you might consider the UC Davis Womens Frisbee Team or UC Davis Mens Frisbee Team, both part of the UCD Intramural Sports & Sport Clubs.  Campus also has the UCD Frisbee Golf Course.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program outputwincorr
use myomp
use head
implicit none

integer(4)::LL(3),LLk(3)
real(4),allocatable::vc3(:,:,:)

real(4)::boxinfo(2,3)
real(4)::fsvec(3)

integer(4)::i,j,k,dir

real(8),allocatable::pk8(:,:)
real(8),allocatable::pk2d(:,:,:)
real(4)::fs
real(4)::npeff
integer(4)::pp
logical(4)::debug=.false.

real(4),parameter::pi=4.*atan(1.)
real(4)::kbasic

LL=[L,L,L]
LLk=[L+2,L,L]
boxinfo=reshape([0,LL(1),0,LL(2),0,LL(3)],[2,3])

call binscheme
call memo(0)

if (sncorr) then
  pp=0
  if (second) then
    if (wincorr) pp=pp+painter
    if (wincorr2) pp=pp+painter2
  else
    pp=pp+painter*2
  endif
  if (interlace) then
    call getshotnoise_interlace_new(vc3,LL,1.,pp)
  else
    call getshotnoise(vc3,LL,1.,pp)
  endif
  call calcpkell(vc3,LL,pk8,kbin,outsnname)
else
  write(*,*) 'shot noise is not corrected'
endif

if (wincorr.or.(wincorr2.and.second)) then
  pp=0
  if (second) then
    if (wincorr) pp=pp+painter
    if (wincorr2) pp=pp+painter2
  else
    if (wincorr) pp=pp+painter*2
  endif

  call getnormalization(vc3,LL,pp)
  call calcpkell(vc3,LL,pk8,kbin,outwinname)
endif

call memo(1)

contains

subroutine binscheme
  implicit none
  kbasic=2*pi/box
  write(*,'(a,f8.3,a,f8.3)') 'kbasic=',kbasic,'  knyquist=',kbasic*L/2
  if (kbin.eq.0.and.dk.eq.0.) then
    write(*,*) 'ERROR: kbin and dk cannot to 0 in the same time'
    stop
  endif

  if (usedefaulkminkmax) then
    if (logbin) then
      kmax=alog10(kbasic*L/2)
      kmin=alog10(kbasic*1.28)
    else
      kmax=kbasic*(L+1)/2
      kmin=kbasic/2
    endif
  else
    if (logbin) then
      kmax=alog10(kmax)
      if (kmin.eq.0) then
        write(*,*) 'ERROR: kmin=0, set kmin=1.28*kbasic'
        kmin=alog10(kbasic*1.28)
      else
        kmin=alog10(kmin)
      endif
    endif
  endif

  if (kbin.ne.0) then
    dk=(kmax-kmin)/kbin
  else
    kbin=floor((kmax-kmin)/dk)
    kmax=kmin+dk*kbin
  endif
  if (logbin) then
    write(*,'(a,f8.3,a,f8.3)') 'kmin=',10.**kmin,'  kmax=',10.**kmax
    write(*,'(a,f8.3,a,i8)') 'dlogk=',dk,'  kbin=',kbin
  else
    write(*,'(a,f8.3,a,f8.3)') 'kmin=',kmin,'  kmax=',kmax
    write(*,'(a,f8.3,a,i8)') 'dk=',dk,'  kbin=',kbin
  endif
endsubroutine binscheme

subroutine memo(command)
  implicit none
  integer(4)::command
  if (command.eq.0) then
    write(*,*) 'memo for vc3:',(float(LL(1)+2)*LL(2)*LL(3))*4/1024.**3,'G'
    allocate(vc3(L+2,L,L))
    allocate(pk8(9,kbin))
  elseif(command.eq.1) then
    deallocate(vc3)
    deallocate(pk8)
  endif
endsubroutine memo

subroutine getnormalization(vc,LL,pp)
  implicit none
  integer(4)::LL(3)
  real(4)::vc(LL(1)+2,LL(2),LL(3))
  real(4)::hpi(3),a1,a2,a3,d2,ri,rj,rk
  real(4)::wx(LL(1)/2+1),wy(LL(2)),wz(LL(3))
  integer(4)::pp
  if (pp.eq.0) then
    write(*,*) 'no window function correction'
    return
  endif
  write(*,*) 'window function is deconvolved, pp=',pp
  hpi=pi/float(LL)

  wz(0)=1
  !$omp parallel do default(private) shared(hpi,LL,wz) schedule(static)
  do k=2,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    a3=hpi(3)*rk
    wz(k)=sin(a3)/a3
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(hpi,LL,wy) schedule(static)
  do j=1,LL(2)
    if (j.le.LL(2)/2+1) then
      rj=j-1
    else
      rj=j-1-LL(2)
    endif
    a2=hpi(2)*rj
    wy(j)=sin(a2)/a2
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(hpi,LL,wx) schedule(static)
  do i=1,LL(1)+2,2
    ri=i/2
    a1=hpi(1)*ri
    wx(i/2+1)=sin(a1)/a1
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(vc,wx,wy,wz,pp) schedule(static)
  do k=1,LL(3)
    do j=1,LL(2)
      do i=1,LL(1)+2,2
        d2=(wx(i/2+1)*wy(j)*wz(k))**pp
        vc(i:i+1,j,k)=d2
      enddo
    enddo
  enddo
  !$omp end parallel do
endsubroutine getnormalization

subroutine getshotnoise(vc,LL,npeff,pp)
  implicit none
  integer(4)::LL(3)
  real(4)::vc(LL(1)+2,LL(2),LL(3))
  real(4)::npeff
  integer(4)::pp
  real(4)::hpi(3),a1,a2,a3,d2,ri,rj,rk
  real(4)::wx(LL(1)/2+1),wy(LL(2)),wz(LL(3))
  if (npeff.eq.0) then
    write(*,*) 'can not correct shot noise with unknown particle number'
    return
  endif
  write(*,*) 'shot noise is corrected'
  hpi=pi/float(LL)

  !$omp parallel do default(private) shared(hpi,LL,pp,wz) schedule(static)
  do k=1,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    a3=sin(hpi(3)*rk)
    wz(k)=wcorr(a3,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(hpi,LL,pp,wy) schedule(static)
  do j=1,LL(2)
    if (j.le.LL(2)/2+1) then
      rj=j-1
    else
      rj=j-1-LL(2)
    endif
    a2=sin(hpi(2)*rj)
    wy(j)=wcorr(a2,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(hpi,LL,pp,wx) schedule(static)
  do i=1,LL(1)+2,2
    ri=i/2
    a1=sin(hpi(1)*ri)
    wx(i/2+1)=wcorr(a1,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(vc,wx,wy,wz,LL,npeff) schedule(static)
  do k=1,LL(3)
    do j=1,LL(2)
      do i=1,LL(1)+2,2
        d2 = wx(i/2+1)*wy(j)*wz(k)
        vc(i:i+1,j,k)=d2/npeff
      enddo
    enddo
  enddo
  !$omp end parallel do
endsubroutine getshotnoise

function wcorr(aa,pp)
  implicit none
  real(4)::wcorr
  real(4)::aa  ! sinx
  integer(4)::pp
  real(4)::bb  ! cosx
  if (mod(pp,2).ne.0) bb=sqrt(1.-aa**2)
  if (pp.eq.2) then
    wcorr=1
  elseif (pp.eq.3) then
    wcorr=bb
  elseif (pp.eq.4) then
    wcorr=1-2./3*aa**2
  elseif (pp.eq.5) then
    wcorr=bb*(1-1./3*aa**2)
  elseif (pp.eq.6) then
    wcorr=1-aa**2+2./15*aa**4
  elseif (pp.eq.7) then
    wcorr=bb*(1-2./3*aa**2+2./45*aa**4)
  elseif(pp.eq.8) then
    wcorr=1-4./3*aa**2+2./5*aa**4-4./315*aa**6
  else
    wcorr=1
  endif
endfunction wcorr

subroutine getshotnoise_interlace(vc,LL,npeff,pp)
  implicit none
  integer(4)::LL(3)
  real(4)::vc(LL(1)+2,LL(2),LL(3))
  real(4)::npeff
  integer(4)::pp
  real(4)::hpi(3),a1,a2,a3,d2,ri,rj,rk
  real(4)::wx(LL(1)/2+1),wy(LL(2)),wz(LL(3))
  if (npeff.eq.0) then
    write(*,*) 'can not correct shot noise with unknown particle number'
    return
  endif
  write(*,*) 'shot noise is corrected (interlaced)'
  write(*,*) 'pp=',pp
  write(*,*) 'npeff=',npeff
  hpi=pi/float(LL)/2
  !$omp parallel do default(private) shared(vc,hpi,LL,npeff,painter,pp) schedule(static)
  do k=1,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    a3=sin(hpi(3)*rk)
    a3=wcorr_interlace(a3,pp)

    do j=1,LL(2)
      if (j.le.LL(2)/2+1) then
        rj=j-1
      else
        rj=j-1-LL(2)
      endif
      a2=sin(hpi(2)*rj)
      a2=wcorr_interlace(a2,pp)

      do i=1,LL(1)+2,2
        ri=i/2
        a1=sin(hpi(1)*ri)
        a1=wcorr_interlace(a1,pp)

        vc(i:i+1,j,k)=a1*a2*a3/npeff
      enddo
    enddo
  enddo
  !$omp end parallel do
endsubroutine getshotnoise_interlace

function wcorr_interlace(aa,pp)
  implicit none
  real(4)::wcorr_interlace
  real(4)::aa  ! sinx
  integer(4)::pp
  real(4)::bb  ! cosx
  bb=sqrt(1.-aa**2)
  if (pp.eq.2) then
    wcorr_interlace=bb**2
  elseif (pp.eq.3) then
    wcorr_interlace=bb**4
  elseif (pp.eq.4) then
    wcorr_interlace=(1-2./3*aa**2)*bb**4
  elseif (pp.eq.5) then
    wcorr_interlace=(1-1./3*aa**2)*bb**6
  elseif (pp.eq.6) then
    wcorr_interlace=(1-aa**2+2./15*aa**4)*bb**6
  elseif (pp.eq.7) then
    wcorr_interlace=(1-2./3*aa**2+2./45*aa**4)*bb**8
  elseif(pp.eq.8) then
    wcorr_interlace=(1-4./3*aa**2+2./5*aa**4-4./315*aa**6)*bb**8
  else
    wcorr_interlace=1
  endif
endfunction wcorr_interlace

subroutine getshotnoise_interlace_new(vc,LL,npeff,pp)
  implicit none
  integer(4)::LL(3)
  real(4)::vc(LL(1)+2,LL(2),LL(3))
  real(4)::npeff
  integer(4)::pp
  real(4)::hpi(3),a1,a2,a3,d2,ri,rj,rk
  real(4)::wx(LL(1)/2+1,2),wy(LL(2),2),wz(LL(3),2)
  if (npeff.eq.0) then
    write(*,*) 'can not correct shot noise with unknown particle number'
    return
  endif
  write(*,*) 'shot noise is corrected (interlaced)'
  write(*,*) 'pp=',pp
  write(*,*) 'npeff=',npeff
  hpi=pi/float(LL)

  !$omp parallel do default(private) shared(LL,hpi,wz,pp) schedule(static)
  do k=1,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    a3=hpi(3)*rk/2
    wz(k,1:2)=wcorr_interlace_new(a3,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(LL,hpi,wy,pp) schedule(static)
  do j=1,LL(2)
    if (j.le.LL(2)/2+1) then
      rj=j-1
    else
      rj=j-1-LL(2)
    endif
    a2=hpi(2)*rj/2
    wy(j,1:2)=wcorr_interlace_new(a2,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(LL,hpi,wx,pp) schedule(static)
  do i=1,LL(1)+2,2
    ri=i/2
    a1=hpi(1)*ri/2
    wx(i/2+1,1:2)=wcorr_interlace_new(a1,pp)
  enddo
  !$omp end parallel do

  !$omp parallel do default(private) shared(vc,LL,wx,wy,wz,npeff) schedule(static)
  do k=1,LL(3)
    do j=1,LL(2)
      do i=1,LL(1)+2,2
        d2= wx(i/2+1,1)*wy(j,1)*wz(k,1) + wx(i/2+1,1)*wy(j,2)*wz(k,2) + wx(i/2+1,2)*wy(j,1)*wz(k,2) + wx(i/2+1,2)*wy(j,2)*wz(k,1)
        vc(i:i+1,j,k)=d2/npeff
      enddo
    enddo
  enddo
  !$omp end parallel do
endsubroutine getshotnoise_interlace_new

function wcorr_interlace_new(aa,pp)
  implicit none
  real(4)::wcorr_interlace_new(2)
  real(4)::aa
  integer(4)::pp
  real(4)::sinx,cosx
  sinx=sin(aa)
  cosx=cos(aa)

  if (pp.eq.2) then
    wcorr_interlace_new(1)=cosx**2
    wcorr_interlace_new(2)=sinx**2
  !elseif (pp.eq.3) then
  !  wcorr_interlace_new(1)=coshx**4
  elseif (pp.eq.4) then
    wcorr_interlace_new(1)=(1-2./3*sinx**2)*cosx**4
    wcorr_interlace_new(2)=(1-2./3*cosx**2)*sinx**4
  !elseif (pp.eq.5) then
  !  wcorr_interlace_new(1)=(1-1./3*sinhx**2)*coshx**6
  elseif (pp.eq.6) then
    wcorr_interlace_new(1)=(1-sinx**2+2./15*sinx**4)*cosx**6
    wcorr_interlace_new(2)=(1-cosx**2+2./15*cosx**4)*sinx**6
  !elseif (pp.eq.7) then
  !  wcorr_interlace_new(1)=(1-2./3*sinhx**2+2./45*sinhx**4)*coshx**8
  elseif(pp.eq.8) then
    wcorr_interlace_new(1)=(1-4./3*sinx**2+2./5*sinx**4-4./315*sinx**6)*cosx**8
    wcorr_interlace_new(2)=(1-4./3*cosx**2+2./5*cosx**4-4./315*cosx**6)*sinx**8
  else
    wcorr_interlace_new(1:2)=1
  endif
endfunction wcorr_interlace_new

subroutine calcpkell(vc,LL,pk8,kbin,filename)
  implicit none
  integer(4)::LL(3)
  real(4)::vc(LL(1)+2,LL(2),LL(3))
  integer(4)::kbin
  character(*)::filename
  real(8)::pk8(9,kbin)
  real(4)::ri,rj,rk,kmode,mu
  real(4)::pl0,pl2,pl4
  integer(4)::bin
  real(4)::dr
  real(4)::factor

  pk8=0.
  write(*,*) 'calculating power spectrum'

  !$omp parallel do default(private) shared(vc,LL,kbasic,kmin,dk,kbin,logbin) &
  !$omp reduction(+:pk8) schedule(guided)
  do k=1,LL(3)
    if (k.le.LL(3)/2+1) then
      rk=k-1
    else
      rk=k-1-LL(3)
    endif
    rk=kbasic*rk
    do j=1,LL(2)
      if (j.le.LL(2)/2+1) then
        rj=j-1
      else
        rj=j-1-LL(2)
      endif
      rj=kbasic*rj
      do i=1,LL(1)+2,2
        ri=i/2
        ri=kbasic*ri

        if (i.eq.1.and.j.eq.1.and.k.eq.1) cycle
        if (i.eq.1) then
          factor=1
        else
          factor=2
        endif
        kmode=sqrt(ri*ri+rj*rj+rk*rk)
        mu=rk/kmode
        if (logbin) kmode=alog10(kmode)
        bin=floor((kmode-kmin)/dk)+1
        if (bin.lt.1.or.bin.gt.kbin) cycle
        pl0=1.*vc(i,j,k)
        pl2=0.5*(3*mu**2-1)*vc(i,j,k)
        pl4=0.125*(35*mu**4-30*mu**2+3)*vc(i,j,k)
        pk8(1,bin)=pk8(1,bin)+kmode*factor
        pk8(2,bin)=pk8(2,bin)+kmode**2*factor
        pk8(3,bin)=pk8(3,bin)+pl0*factor
        pk8(4,bin)=pk8(4,bin)+pl0**2*factor
        pk8(5,bin)=pk8(5,bin)+pl2*factor
        pk8(6,bin)=pk8(6,bin)+pl2**2*factor
        pk8(7,bin)=pk8(7,bin)+pl4*factor
        pk8(8,bin)=pk8(8,bin)+pl4**2*factor
        pk8(9,bin)=pk8(9,bin)+1.*factor
      enddo
    enddo
  enddo
  !$omp end parallel do

  pk8(1,:)=pk8(1,:)/pk8(9,:)
  pk8(2,:)=pk8(2,:)/pk8(9,:)-pk8(1,:)**2
  pk8(3,:)=pk8(3,:)/pk8(9,:)
  pk8(4,:)=pk8(4,:)/pk8(9,:)-pk8(3,:)**2
  pk8(5,:)=pk8(5,:)/pk8(9,:)
  pk8(6,:)=pk8(6,:)/pk8(9,:)-pk8(5,:)**2
  pk8(7,:)=pk8(7,:)/pk8(9,:)
  pk8(8,:)=pk8(8,:)/pk8(9,:)-pk8(7,:)**2

  if (usecenter) then
    do i=1,kbin
      pk8(1,i)=kmin+(i-0.5)*dk
    enddo
    if (logbin) pk8(1,:)=10.**pk8(1,:)
    pk8(2,:)=0.
  else
    if (logbin) then
      pk8(1,:)=10.**pk8(1,:)
      pk8(2,:)=(10.**sqrt(pk8(2,:))-1.)*pk8(1,:)
    endif
  endif

  write(*,*) 'writing: ',trim(filename)
  open(32,file=filename,form='formatted',status='replace')
  do i=1,kbin
    write(32,'(9e14.6)') pk8(1:9:2,i)
    if (debug) write(*,*) real(pk8(1:9:2,i))
  enddo
  close(32)
endsubroutine calcpkell

endprogram outputwincorr

max_stars_repo_name: ColdThunder/FFTps

max_issues_repo_name: ColdThunder/FFTps

max_forks_repo_name: ColdThunder/FFTps


text: Describe Users/ScarlettW here.


text:
program main
    use mpi

    integer :: my_rank,petot,ierr

    call mpi_init(ierr)
    call mpi_comm_rank(mpi_comm_world,my_rank,ierr)
    call mpi_comm_size(mpi_comm_world,petot,ierr)
    print *, "My rank is :",my_rank,"Total process is : ",petot
    call mpi_finalize(ierr)

end program main

max_stars_repo_name: haruka-soybean/xFortran

max_issues_repo_name: haruka-soybean/xFortran

max_forks_repo_name: haruka-soybean/xFortran