!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006, 2010, 2015, 2017
!        Salvatore Filippone    Cranfield University
!        Alfredo Buttari        CNRS-IRIT, Toulouse
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
! package: psb_i_base_vect_mod
!
! This module contains the definition of the psb_i_base_vect type which
! is a container for dense vectors.
!  This is encapsulated instead of being just a simple array to allow for
!  more complicated situations, such as GPU programming, where the memory
!  area we are interested in is not easily accessible from the host/Fortran
!  side. It is also meant to be encapsulated in an outer type, to allow
!  runtime switching as per the STATE design pattern, similar to the
!  sparse matrix types.
!
!
module psb_i_base_vect_mod

  use psb_const_mod
  use psb_error_mod
  use psb_realloc_mod

  !> \namespace  psb_base_mod  \class psb_i_base_vect_type
  !! The psb_i_base_vect_type
  !! defines a middle level  integer(psb_ipk_) encapsulated dense vector.
  !! The encapsulation is needed, in place of a simple array, to allow
  !! for complicated situations, such as GPU programming, where the memory
  !!  area we are interested in is not easily accessible from the host/Fortran
  !!  side. It is also meant to be encapsulated in an outer type, to allow
  !!  runtime switching as per the STATE design pattern, similar to the
  !!  sparse matrix types.
  !!
  type psb_i_base_vect_type
    !> Values.
    integer(psb_ipk_), allocatable :: v(:)
    integer(psb_ipk_), allocatable :: combuf(:)
    integer(psb_mpik_), allocatable :: comid(:,:)
  contains
    !
    !  Constructors/allocators
    !
    procedure, pass(x) :: bld_x    => i_base_bld_x
    procedure, pass(x) :: bld_n    => i_base_bld_n
    generic, public    :: bld      => bld_x, bld_n
    procedure, pass(x) :: all      => i_base_all
    procedure, pass(x) :: mold     => i_base_mold
    !
    ! Insert/set. Assembly and free.
    ! Assembly does almost nothing here, but is important
    ! in derived classes.
    !
    procedure, pass(x) :: ins_a    => i_base_ins_a
    procedure, pass(x) :: ins_v    => i_base_ins_v
    generic, public    :: ins      => ins_a, ins_v
    procedure, pass(x) :: zero     => i_base_zero
    procedure, pass(x) :: asb      => i_base_asb
    procedure, pass(x) :: free     => i_base_free
    !
    ! Sync: centerpiece of handling of external storage.
    ! Any derived class having extra storage upon sync
    ! will guarantee that both fortran/host side and
    ! external side contain the same data. The base
    ! version is only a placeholder.
    !
    procedure, pass(x) :: sync     => i_base_sync
    procedure, pass(x) :: is_host  => i_base_is_host
    procedure, pass(x) :: is_dev   => i_base_is_dev
    procedure, pass(x) :: is_sync  => i_base_is_sync
    procedure, pass(x) :: set_host => i_base_set_host
    procedure, pass(x) :: set_dev  => i_base_set_dev
    procedure, pass(x) :: set_sync => i_base_set_sync

    !
    ! These are for handling gather/scatter in new
    ! comm internals implementation.
    !
    procedure, nopass  :: use_buffer   => i_base_use_buffer
    procedure, pass(x) :: new_buffer   => i_base_new_buffer
    procedure, nopass  :: device_wait  => i_base_device_wait
    procedure, pass(x) :: maybe_free_buffer  => i_base_maybe_free_buffer
    procedure, pass(x) :: free_buffer  => i_base_free_buffer
    procedure, pass(x) :: new_comid    => i_base_new_comid
    procedure, pass(x) :: free_comid   => i_base_free_comid

    !
    ! Basic info
    procedure, pass(x) :: get_nrows => i_base_get_nrows
    procedure, pass(x) :: sizeof    => i_base_sizeof
    procedure, nopass  :: get_fmt   => i_base_get_fmt
    !
    ! Set/get data from/to an external array; also
    ! overload assignment.
    !
    procedure, pass(x) :: get_vect => i_base_get_vect
    procedure, pass(x) :: set_scal => i_base_set_scal
    procedure, pass(x) :: set_vect => i_base_set_vect
    generic, public    :: set      => set_vect, set_scal
    !
    ! Gather/scatter. These are needed for MPI interfacing.
    ! May have to be reworked.
    !
    procedure, pass(x) :: gthab    => i_base_gthab
    procedure, pass(x) :: gthzv    => i_base_gthzv
    procedure, pass(x) :: gthzv_x  => i_base_gthzv_x
    procedure, pass(x) :: gthzbuf  => i_base_gthzbuf
    generic, public    :: gth      => gthab, gthzv, gthzv_x, gthzbuf
    procedure, pass(y) :: sctb     => i_base_sctb
    procedure, pass(y) :: sctb_x   => i_base_sctb_x
    procedure, pass(y) :: sctb_buf => i_base_sctb_buf
    generic, public    :: sct      => sctb, sctb_x, sctb_buf

  end type psb_i_base_vect_type

  public  :: psb_i_base_vect
  private :: constructor, size_const
  interface psb_i_base_vect
    module procedure constructor, size_const
  end interface psb_i_base_vect

contains

  !
  ! Constructors.
  !

  !> Function  constructor:
  !! \brief     Constructor from an array
  !!  \param   x(:)  input array to be copied
  !!
  function constructor(x) result(this)
    integer(psb_ipk_)   :: x(:)
    type(psb_i_base_vect_type) :: this
    integer(psb_ipk_) :: info

    this%v = x
    call this%asb(size(x,kind=psb_ipk_),info)
  end function constructor

  !> Function  constructor:
  !! \brief     Constructor from size
  !!  \param    n   Size of vector to be built.
  !!
  function size_const(n) result(this)
    integer(psb_ipk_), intent(in) :: n
    type(psb_i_base_vect_type) :: this
    integer(psb_ipk_) :: info

    call this%asb(n,info)

  end function size_const

  !
  ! Build from a sample
  !

  !> Function  bld_x:
  !! \memberof  psb_i_base_vect_type
  !! \brief     Build method from an array
  !!  \param   x(:)  input array to be copied
  !!
  subroutine i_base_bld_x(x,this)
    use psb_realloc_mod
    integer(psb_ipk_), intent(in) :: this(:)
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(size(this),x%v,info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
      return
    end if
    x%v(:)  = this(:)

  end subroutine i_base_bld_x

  !
  ! Create with size, but no initialization
  !

  !> Function  bld_n:
  !! \memberof  psb_i_base_vect_type
  !! \brief     Build method with size (uninitialized data)
  !!  \param    n    size to be allocated.
  !!
  subroutine i_base_bld_n(x,n)
    use psb_realloc_mod
    integer(psb_ipk_), intent(in) :: n
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(n,x%v,info)
    call x%asb(n,info)

  end subroutine i_base_bld_n

  !> Function  base_all:
  !! \memberof  psb_i_base_vect_type
  !! \brief     Build method with size (uninitialized data) and
  !!            allocation return code.
  !!  \param    n    size to be allocated.
  !!  \param info  return code
  !!
  subroutine i_base_all(n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    integer(psb_ipk_), intent(in)               :: n
    class(psb_i_base_vect_type), intent(out)    :: x
    integer(psb_ipk_), intent(out)              :: info

    call psb_realloc(n,x%v,info)

  end subroutine i_base_all

  !> Function  base_mold:
  !! \memberof  psb_i_base_vect_type
  !! \brief     Mold method: return a variable with the same dynamic type
  !!  \param    y returned variable
  !!  \param info  return code
  !!
  subroutine i_base_mold(x, y, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(in)   :: x
    class(psb_i_base_vect_type), intent(out), allocatable :: y
    integer(psb_ipk_), intent(out)              :: info

    allocate(psb_i_base_vect_type :: y, stat=info)

  end subroutine i_base_mold

  !
  ! Insert a bunch of values at specified positions.
  !
  !> Function  base_ins:
  !! \memberof  psb_i_base_vect_type
  !! \brief Insert coefficients.
  !!
  !!
  !!         Given  a list of N pairs
  !!           (IRL(i),VAL(i))
  !!         record a new coefficient in X such that
  !!            X(IRL(1:N)) = VAL(1:N).
  !!
  !!         - the update operation will perform either
  !!               X(IRL(1:n)) = VAL(1:N)
  !!           or
  !!               X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N)
  !!           according to the value of DUPLICATE.
  !!
  !!
  !!  \param n     number of pairs in input
  !!  \param irl(:)  the input row indices
  !!  \param val(:)  the input coefficients
  !!  \param dupl    how to treat duplicate entries
  !!  \param info  return code
  !!
  !
  subroutine i_base_ins_a(n,irl,val,dupl,x,info)
    use psi_serial_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in)               :: n, dupl
    integer(psb_ipk_), intent(in)               :: irl(:)
    integer(psb_ipk_), intent(in)        :: val(:)
    integer(psb_ipk_), intent(out)              :: info

    integer(psb_ipk_) :: i, isz

    info = 0
    if (psb_errstatus_fatal()) return

    if (.not.allocated(x%v)) then
      info = psb_err_invalid_vect_state_
    else if (n > min(size(irl),size(val))) then
      info = psb_err_invalid_input_

    else
      isz = size(x%v)
      select case(dupl)
      case(psb_dupl_ovwrt_)
        do i = 1, n
          !loop over all val's rows

          ! row actual block row
          if ((1 <= irl(i)).and.(irl(i) <= isz)) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i)) = val(i)
          end if
        enddo

      case(psb_dupl_add_)

        do i = 1, n
          !loop over all val's rows
          if ((1 <= irl(i)).and.(irl(i) <= isz)) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i)) = x%v(irl(i)) +  val(i)
          end if
        enddo

      case default
        info = 321
        ! !$      call psb_errpush(info,name)
        ! !$      goto 9999
      end select
    end if
    call x%set_host()
    if (info /= 0) then
      call psb_errpush(info,'base_vect_ins')
      return
    end if

  end subroutine i_base_ins_a

  subroutine i_base_ins_v(n,irl,val,dupl,x,info)
    use psi_serial_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in)               :: n, dupl
    class(psb_i_base_vect_type), intent(inout)  :: irl
    class(psb_i_base_vect_type), intent(inout)  :: val
    integer(psb_ipk_), intent(out)              :: info

    integer(psb_ipk_) :: i, isz

    info = 0
    if (psb_errstatus_fatal()) return

    if (irl%is_dev()) call irl%sync()
    if (val%is_dev()) call val%sync()
    if (x%is_dev())   call x%sync()
    call x%ins(n,irl%v,val%v,dupl,info)

    if (info /= 0) then
      call psb_errpush(info,'base_vect_ins')
      return
    end if

  end subroutine i_base_ins_v

  !
  !> Function  base_zero
  !! \memberof  psb_i_base_vect_type
  !! \brief Zero out contents
  !!
  !
  subroutine i_base_zero(x)
    use psi_serial_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout)    :: x

    if (allocated(x%v)) x%v=izero
    call x%set_host()
  end subroutine i_base_zero

  !
  ! Assembly.
  ! For derived classes: after this the vector
  ! storage is supposed to be in sync.
  !
  !> Function  base_asb:
  !! \memberof  psb_i_base_vect_type
  !! \brief Assemble vector: reallocate as necessary.
  !!
  !!  \param n     final size
  !!  \param info  return code
  !!
  !

  subroutine i_base_asb(n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    integer(psb_ipk_), intent(in)              :: n
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    info = 0
    if (x%get_nrows() < n) &
         & call psb_realloc(n,x%v,info)
    if (info /= 0) &
         & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
    call x%sync()
  end subroutine i_base_asb

  !
  !> Function  base_free:
  !! \memberof  psb_i_base_vect_type
  !! \brief Free vector
  !!
  !!  \param info  return code
  !!
  !
  subroutine i_base_free(x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(out)              :: info

    info = 0
    if (allocated(x%v)) deallocate(x%v, stat=info)
    if (info == 0) call x%free_buffer(info)
    if (info == 0) call x%free_comid(info)
    if (info /= 0) call &
         & psb_errpush(psb_err_alloc_dealloc_,'vect_free')

  end subroutine i_base_free

  !
  !> Function  base_free_buffer:
  !! \memberof  psb_i_base_vect_type
  !! \brief Free aux buffer
  !!
  !!  \param info  return code
  !!
  !
  subroutine i_base_free_buffer(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    if (allocated(x%combuf)) &
         &  deallocate(x%combuf,stat=info)
  end subroutine i_base_free_buffer

  !
  !> Function  base_maybe_free_buffer:
  !! \memberof  psb_i_base_vect_type
  !! \brief Conditionally Free aux buffer.
  !!        In some derived classes, e.g. GPU,
  !!        does not really frees to avoid  runtime
  !!        costs
  !!
  !!  \param info  return code
  !!
  !
  subroutine i_base_maybe_free_buffer(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    info = 0
    if (psb_get_maybe_free_buffer())&
         &  call x%free_buffer(info)

  end subroutine i_base_maybe_free_buffer

  !
  !> Function  base_free_comid:
  !! \memberof  psb_i_base_vect_type
  !! \brief Free aux MPI communication id buffer
  !!
  !!  \param info  return code
  !!
  !
  subroutine i_base_free_comid(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    if (allocated(x%comid)) &
         &  deallocate(x%comid,stat=info)
  end subroutine i_base_free_comid

  !
  ! The base version of SYNC & friends does nothing, it's just
  ! a placeholder.
  !
  !
  !> Function  base_sync:
  !! \memberof  psb_i_base_vect_type
  !! \brief Sync: base version is a no-op.
  !!
  !
  subroutine i_base_sync(x)
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x

  end subroutine i_base_sync

  !
  !> Function  base_set_host:
  !! \memberof  psb_i_base_vect_type
  !! \brief Set_host: base version is a no-op.
  !!
  !
  subroutine i_base_set_host(x)
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x

  end subroutine i_base_set_host

  !
  !> Function  base_set_dev:
  !! \memberof  psb_i_base_vect_type
  !! \brief Set_dev: base version is a no-op.
  !!
  !
  subroutine i_base_set_dev(x)
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x

  end subroutine i_base_set_dev

  !
  !> Function  base_set_sync:
  !! \memberof  psb_i_base_vect_type
  !! \brief Set_sync: base version is a no-op.
  !!
  !
  subroutine i_base_set_sync(x)
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x

  end subroutine i_base_set_sync

  !
  !> Function  base_is_dev:
  !! \memberof  psb_i_base_vect_type
  !! \brief Is  vector on external device    .
  !!
  !
  function i_base_is_dev(x) result(res)
    implicit none
    class(psb_i_base_vect_type), intent(in) :: x
    logical  :: res

    res = .false.
  end function i_base_is_dev

  !
  !> Function  base_is_host
  !! \memberof  psb_i_base_vect_type
  !! \brief Is  vector on standard memory    .
  !!
  !
  function i_base_is_host(x) result(res)
    implicit none
    class(psb_i_base_vect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function i_base_is_host

  !
  !> Function  base_is_sync
  !! \memberof  psb_i_base_vect_type
  !! \brief Is  vector on sync               .
  !!
  !
  function i_base_is_sync(x) result(res)
    implicit none
    class(psb_i_base_vect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function i_base_is_sync

  !
  ! Size info.
  !
  !
  !> Function  base_get_nrows
  !! \memberof  psb_i_base_vect_type
  !! \brief  Number of entries
  !!
  !
  function i_base_get_nrows(x) result(res)
    implicit none
    class(psb_i_base_vect_type), intent(in) :: x
    integer(psb_ipk_) :: res

    res = 0
    if (allocated(x%v)) res = size(x%v)

  end function i_base_get_nrows

  !
  !> Function  base_get_sizeof
  !! \memberof  psb_i_base_vect_type
  !! \brief  Size in bytes
  !!
  !
  function i_base_sizeof(x) result(res)
    implicit none
    class(psb_i_base_vect_type), intent(in) :: x
    integer(psb_long_int_k_) :: res

    ! Force 8-byte integers.
    res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows()

  end function i_base_sizeof

  !
  !> Function  base_get_fmt
  !! \memberof  psb_i_base_vect_type
  !! \brief  Format
  !!
  !
  function i_base_get_fmt() result(res)
    implicit none
    character(len=5) :: res
    res = 'BASE'
  end function i_base_get_fmt

  !
  !
  !
  !> Function  base_get_vect
  !! \memberof  psb_i_base_vect_type
  !! \brief  Extract a copy of the contents
  !!
  !
  function  i_base_get_vect(x) result(res)
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), allocatable                 :: res(:)
    integer(psb_ipk_) :: info

    if (.not.allocated(x%v)) return
    if (.not.x%is_host()) call x%sync()
    allocate(res(x%get_nrows()),stat=info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
      return
    end if
    res(:) = x%v(:)
  end function i_base_get_vect

  !
  ! Reset all values
  !
  !
  !> Function  base_set_scal
  !! \memberof  psb_i_base_vect_type
  !! \brief  Set all entries
  !! \param val   The value to set
  !!
  subroutine i_base_set_scal(x,val,first,last)
    class(psb_i_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in) :: val
    integer(psb_ipk_), optional :: first, last

    integer(psb_ipk_) :: info, first_, last_

    first_=1
    last_=size(x%v)
    if (present(first)) first_ = max(1,first)
    if (present(last))  last_  = min(last,last_)

    if (x%is_dev()) call x%sync()
    x%v(first_:last_) = val
    call x%set_host()

  end subroutine i_base_set_scal

  !
  !> Function  base_set_vect
  !! \memberof  psb_i_base_vect_type
  !! \brief  Set all entries
  !! \param val(:)  The vector to be copied in
  !!
  subroutine i_base_set_vect(x,val,first,last)
    class(psb_i_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in) :: val(:)
    integer(psb_ipk_), optional :: first, last

    integer(psb_ipk_) :: info, first_, last_, nr

    first_=1
    last_=min(psb_size(x%v),size(val))
    if (present(first)) first_ = max(1,first)
    if (present(last))  last_  = min(last,last_)

    if (allocated(x%v)) then
      if (x%is_dev()) call x%sync()
      x%v(first_:last_) = val(1:last_-first_+1)
    else
      x%v = val
    end if
    call x%set_host()

  end subroutine i_base_set_vect

  !
  ! Gather: Y = beta * Y + alpha * X(IDX(:))
  !
  !
  !> Function  base_gthab
  !! \memberof  psb_i_base_vect_type
  !! \brief gather into an array
  !!    Y = beta * Y + alpha * X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  !! \param alpha
  !! \param beta
  subroutine i_base_gthab(n,idx,alpha,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) :: alpha, beta, y(:)
    class(psb_i_base_vect_type) :: x

    if (x%is_dev()) call x%sync()
    call psi_gth(n,idx,alpha,x%v,beta,y)

  end subroutine i_base_gthab
  !
  ! shortcut alpha=1 beta=0
  !
  !> Function  base_gthzv
  !! \memberof  psb_i_base_vect_type
  !! \brief gather into an array special alpha=1 beta=0
  !!    Y = X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  subroutine i_base_gthzv_x(i,n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i,n
    class(psb_i_base_vect_type) :: idx
    integer(psb_ipk_) ::  y(:)
    class(psb_i_base_vect_type) :: x

    if (idx%is_dev()) call idx%sync()
    call x%gth(n,idx%v(i:),y)

  end subroutine i_base_gthzv_x

  !
  ! New comm internals impl.
  !
  subroutine i_base_gthzbuf(i,n,idx,x)
    use psi_serial_mod
    integer(psb_ipk_) :: i,n
    class(psb_i_base_vect_type) :: idx
    class(psb_i_base_vect_type) :: x

    if (.not.allocated(x%combuf)) then
      call psb_errpush(psb_err_alloc_dealloc_,'gthzbuf')
      return
    end if
    if (idx%is_dev()) call idx%sync()
    if (x%is_dev()) call x%sync()
    call x%gth(n,idx%v(i:),x%combuf(i:))

  end subroutine i_base_gthzbuf
  !
  !> Function  base_device_wait:
  !! \memberof  psb_i_base_vect_type
  !! \brief device_wait: base version is a no-op.
  !!
  !
  subroutine i_base_device_wait()
    implicit none

  end subroutine i_base_device_wait

  function i_base_use_buffer() result(res)
    logical :: res

    res = .true.
  end function i_base_use_buffer

  subroutine i_base_new_buffer(n,x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)              :: n
    integer(psb_ipk_), intent(out)             :: info

    call psb_realloc(n,x%combuf,info)
  end subroutine i_base_new_buffer

  subroutine i_base_new_comid(n,x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)              :: n
    integer(psb_ipk_), intent(out)             :: info

    call psb_realloc(n,2_psb_ipk_,x%comid,info)
  end subroutine i_base_new_comid

  !
  ! shortcut alpha=1 beta=0
  !
  !> Function  base_gthzv
  !! \memberof  psb_i_base_vect_type
  !! \brief gather into an array special alpha=1 beta=0
  !!    Y = X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  subroutine i_base_gthzv(n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) ::  y(:)
    class(psb_i_base_vect_type) :: x

    if (x%is_dev()) call x%sync()
    call psi_gth(n,idx,x%v,y)

  end subroutine i_base_gthzv

  !
  ! Scatter:
  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
  !
  !
  !> Function  base_sctb
  !! \memberof  psb_i_base_vect_type
  !! \brief scatter into a class(base_vect)
  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  !! \param beta
  !! \param x(:)
  subroutine i_base_sctb(n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) :: beta, x(:)
    class(psb_i_base_vect_type) :: y

    if (y%is_dev()) call y%sync()
    call psi_sct(n,idx,x,beta,y%v)
    call y%set_host()

  end subroutine i_base_sctb

  subroutine i_base_sctb_x(i,n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i, n
    class(psb_i_base_vect_type) :: idx
    integer(psb_ipk_) :: beta, x(:)
    class(psb_i_base_vect_type) :: y

    if (idx%is_dev()) call idx%sync()
    call y%sct(n,idx%v(i:),x,beta)
    call y%set_host()

  end subroutine i_base_sctb_x

  subroutine i_base_sctb_buf(i,n,idx,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i, n
    class(psb_i_base_vect_type) :: idx
    integer(psb_ipk_) :: beta
    class(psb_i_base_vect_type) :: y

    if (.not.allocated(y%combuf)) then
      call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
      return
    end if
    if (y%is_dev()) call y%sync()
    if (idx%is_dev()) call idx%sync()
    call y%sct(n,idx%v(i:),y%combuf(i:),beta)
    call y%set_host()

  end subroutine i_base_sctb_buf

end module psb_i_base_vect_mod

module psb_i_base_multivect_mod

  use psb_const_mod
  use psb_error_mod
  use psb_realloc_mod
  use psb_i_base_vect_mod

  !> \namespace  psb_base_mod  \class psb_i_base_vect_type
  !! The psb_i_base_vect_type
  !! defines a middle level  integer(psb_ipk_) encapsulated dense vector.
  !! The encapsulation is needed, in place of a simple array, to allow
  !! for complicated situations, such as GPU programming, where the memory
  !!  area we are interested in is not easily accessible from the host/Fortran
  !!  side. It is also meant to be encapsulated in an outer type, to allow
  !!  runtime switching as per the STATE design pattern, similar to the
  !!  sparse matrix types.
  !!
  private
  public  :: psb_i_base_multivect, psb_i_base_multivect_type

  type psb_i_base_multivect_type
    !> Values.
    integer(psb_ipk_), allocatable :: v(:,:)
    integer(psb_ipk_), allocatable :: combuf(:)
    integer(psb_mpik_), allocatable :: comid(:,:)
  contains
    !
    !  Constructors/allocators
    !
    procedure, pass(x) :: bld_x    => i_base_mlv_bld_x
    procedure, pass(x) :: bld_n    => i_base_mlv_bld_n
    generic, public    :: bld      => bld_x, bld_n
    procedure, pass(x) :: all      => i_base_mlv_all
    procedure, pass(x) :: mold     => i_base_mlv_mold
    !
    ! Insert/set. Assembly and free.
    ! Assembly does almost nothing here, but is important
    ! in derived classes.
    !
    procedure, pass(x) :: ins      => i_base_mlv_ins
    procedure, pass(x) :: zero     => i_base_mlv_zero
    procedure, pass(x) :: asb      => i_base_mlv_asb
    procedure, pass(x) :: free     => i_base_mlv_free
    !
    ! Sync: centerpiece of handling of external storage.
    ! Any derived class having extra storage upon sync
    ! will guarantee that both fortran/host side and
    ! external side contain the same data. The base
    ! version is only a placeholder.
    !
    procedure, pass(x) :: sync     => i_base_mlv_sync
    procedure, pass(x) :: is_host  => i_base_mlv_is_host
    procedure, pass(x) :: is_dev   => i_base_mlv_is_dev
    procedure, pass(x) :: is_sync  => i_base_mlv_is_sync
    procedure, pass(x) :: set_host => i_base_mlv_set_host
    procedure, pass(x) :: set_dev  => i_base_mlv_set_dev
    procedure, pass(x) :: set_sync => i_base_mlv_set_sync

    !
    ! Basic info
    procedure, pass(x) :: get_nrows => i_base_mlv_get_nrows
    procedure, pass(x) :: get_ncols => i_base_mlv_get_ncols
    procedure, pass(x) :: sizeof    => i_base_mlv_sizeof
    procedure, nopass  :: get_fmt   => i_base_mlv_get_fmt
    !
    ! Set/get data from/to an external array; also
    ! overload assignment.
    !
    procedure, pass(x) :: get_vect => i_base_mlv_get_vect
    procedure, pass(x) :: set_scal => i_base_mlv_set_scal
    procedure, pass(x) :: set_vect => i_base_mlv_set_vect
    generic, public    :: set      => set_vect, set_scal

    !
    ! These are for handling gather/scatter in new
    ! comm internals implementation.
    !
    procedure, nopass  :: use_buffer   => i_base_mlv_use_buffer
    procedure, pass(x) :: new_buffer   => i_base_mlv_new_buffer
    procedure, nopass  :: device_wait  => i_base_mlv_device_wait
    procedure, pass(x) :: maybe_free_buffer  => i_base_mlv_maybe_free_buffer
    procedure, pass(x) :: free_buffer  => i_base_mlv_free_buffer
    procedure, pass(x) :: new_comid    => i_base_mlv_new_comid
    procedure, pass(x) :: free_comid   => i_base_mlv_free_comid

    !
    ! Gather/scatter. These are needed for MPI interfacing.
    ! May have to be reworked.
    !
    procedure, pass(x) :: gthab    => i_base_mlv_gthab
    procedure, pass(x) :: gthzv    => i_base_mlv_gthzv
    procedure, pass(x) :: gthzm    => i_base_mlv_gthzm
    procedure, pass(x) :: gthzv_x  => i_base_mlv_gthzv_x
    procedure, pass(x) :: gthzbuf  => i_base_mlv_gthzbuf
    generic, public    :: gth      => gthab, gthzv, gthzm, gthzv_x, gthzbuf
    procedure, pass(y) :: sctb     => i_base_mlv_sctb
    procedure, pass(y) :: sctbr2   => i_base_mlv_sctbr2
    procedure, pass(y) :: sctb_x   => i_base_mlv_sctb_x
    procedure, pass(y) :: sctb_buf => i_base_mlv_sctb_buf
    generic, public    :: sct      => sctb, sctbr2, sctb_x, sctb_buf
  end type psb_i_base_multivect_type

  interface psb_i_base_multivect
    module procedure constructor, size_const
  end interface psb_i_base_multivect

contains

  !
  ! Constructors.
  !

  !> Function  constructor:
  !! \brief     Constructor from an array
  !!  \param   x(:)  input array to be copied
  !!
  function constructor(x) result(this)
    integer(psb_ipk_)   :: x(:,:)
    type(psb_i_base_multivect_type) :: this
    integer(psb_ipk_) :: info

    this%v = x
    call this%asb(size(x,dim=1,kind=psb_ipk_),size(x,dim=2,kind=psb_ipk_),info)
  end function constructor

  !> Function  constructor:
  !! \brief     Constructor from size
  !!  \param    n   Size of vector to be built.
  !!
  function size_const(m,n) result(this)
    integer(psb_ipk_), intent(in) :: m,n
    type(psb_i_base_multivect_type) :: this
    integer(psb_ipk_) :: info

    call this%asb(m,n,info)

  end function size_const

  !
  ! Build from a sample
  !

  !> Function  bld_x:
  !! \memberof  psb_i_base_multivect_type
  !! \brief     Build method from an array
  !!  \param   x(:)  input array to be copied
  !!
  subroutine i_base_mlv_bld_x(x,this)
    use psb_realloc_mod
    integer(psb_ipk_), intent(in) :: this(:,:)
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(size(this,1),size(this,2),x%v,info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
      return
    end if
    x%v(:,:)  = this(:,:)

  end subroutine i_base_mlv_bld_x

  !
  ! Create with size, but no initialization
  !

  !> Function  bld_n:
  !! \memberof  psb_i_base_multivect_type
  !! \brief     Build method with size (uninitialized data)
  !!  \param    n    size to be allocated.
  !!
  subroutine i_base_mlv_bld_n(x,m,n)
    use psb_realloc_mod
    integer(psb_ipk_), intent(in) :: m,n
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(m,n,x%v,info)
    call x%asb(m,n,info)

  end subroutine i_base_mlv_bld_n

  !> Function  base_mlv_all:
  !! \memberof  psb_i_base_multivect_type
  !! \brief     Build method with size (uninitialized data) and
  !!            allocation return code.
  !!  \param    n    size to be allocated.
  !!  \param info  return code
  !!
  subroutine i_base_mlv_all(m,n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    integer(psb_ipk_), intent(in)               :: m,n
    class(psb_i_base_multivect_type), intent(out) :: x
    integer(psb_ipk_), intent(out)              :: info

    call psb_realloc(m,n,x%v,info)

  end subroutine i_base_mlv_all

  !> Function  base_mlv_mold:
  !! \memberof  psb_i_base_multivect_type
  !! \brief     Mold method: return a variable with the same dynamic type
  !!  \param    y returned variable
  !!  \param info  return code
  !!
  subroutine i_base_mlv_mold(x, y, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(in)   :: x
    class(psb_i_base_multivect_type), intent(out), allocatable :: y
    integer(psb_ipk_), intent(out)              :: info

    allocate(psb_i_base_multivect_type :: y, stat=info)

  end subroutine i_base_mlv_mold

  !
  ! Insert a bunch of values at specified positions.
  !
  !> Function  base_mlv_ins:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Insert coefficients.
  !!
  !!
  !!         Given  a list of N pairs
  !!           (IRL(i),VAL(i))
  !!         record a new coefficient in X such that
  !!            X(IRL(1:N)) = VAL(1:N).
  !!
  !!         - the update operation will perform either
  !!               X(IRL(1:n)) = VAL(1:N)
  !!           or
  !!               X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N)
  !!           according to the value of DUPLICATE.
  !!
  !!
  !!  \param n     number of pairs in input
  !!  \param irl(:)  the input row indices
  !!  \param val(:)  the input coefficients
  !!  \param dupl    how to treat duplicate entries
  !!  \param info  return code
  !!
  !
  subroutine i_base_mlv_ins(n,irl,val,dupl,x,info)
    use psi_serial_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in)               :: n, dupl
    integer(psb_ipk_), intent(in)               :: irl(:)
    integer(psb_ipk_), intent(in)        :: val(:,:)
    integer(psb_ipk_), intent(out)              :: info

    integer(psb_ipk_) :: i, isz

    info = 0
    if (psb_errstatus_fatal()) return

    if (.not.allocated(x%v)) then
      info = psb_err_invalid_vect_state_
    else if (n > min(size(irl),size(val))) then
      info = psb_err_invalid_input_

    else
      isz = size(x%v,1)
      select case(dupl)
      case(psb_dupl_ovwrt_)
        do i = 1, n
          !loop over all val's rows

          ! row actual block row
          if ((1 <= irl(i)).and.(irl(i) <= isz)) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i),:) = val(i,:)
          end if
        enddo

      case(psb_dupl_add_)

        do i = 1, n
          !loop over all val's rows
          if ((1 <= irl(i)).and.(irl(i) <= isz)) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i),:) = x%v(irl(i),:) +  val(i,:)
          end if
        enddo

      case default
        info = 321
        ! !$      call psb_errpush(info,name)
        ! !$      goto 9999
      end select
    end if
    if (info /= 0) then
      call psb_errpush(info,'base_mlv_vect_ins')
      return
    end if

  end subroutine i_base_mlv_ins

  !
  !> Function  base_mlv_zero
  !! \memberof  psb_i_base_multivect_type
  !! \brief Zero out contents
  !!
  !
  subroutine i_base_mlv_zero(x)
    use psi_serial_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout)    :: x

    if (allocated(x%v)) x%v=izero

  end subroutine i_base_mlv_zero

  !
  ! Assembly.
  ! For derived classes: after this the vector
  ! storage is supposed to be in sync.
  !
  !> Function  base_mlv_asb:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Assemble vector: reallocate as necessary.
  !!
  !!  \param n     final size
  !!  \param info  return code
  !!
  !

  subroutine i_base_mlv_asb(m,n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    integer(psb_ipk_), intent(in)              :: m,n
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    if ((x%get_nrows() < m).or.(x%get_ncols()<n)) &
         & call psb_realloc(m,n,x%v,info)
    if (info /= 0) &
         & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')

  end subroutine i_base_mlv_asb

  !
  !> Function  base_mlv_free:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Free vector
  !!
  !!  \param info  return code
  !!
  !
  subroutine i_base_mlv_free(x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(out)              :: info

    info = 0
    if (allocated(x%v)) deallocate(x%v, stat=info)
    if (info /= 0) call &
         & psb_errpush(psb_err_alloc_dealloc_,'vect_free')

  end subroutine i_base_mlv_free

  !
  ! The base version of SYNC & friends does nothing, it's just
  ! a placeholder.
  !
  !
  !> Function  base_mlv_sync:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Sync: base version is a no-op.
  !!
  !
  subroutine i_base_mlv_sync(x)
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x

  end subroutine i_base_mlv_sync

  !
  !> Function  base_mlv_set_host:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Set_host: base version is a no-op.
  !!
  !
  subroutine i_base_mlv_set_host(x)
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x

  end subroutine i_base_mlv_set_host

  !
  !> Function  base_mlv_set_dev:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Set_dev: base version is a no-op.
  !!
  !
  subroutine i_base_mlv_set_dev(x)
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x

  end subroutine i_base_mlv_set_dev

  !
  !> Function  base_mlv_set_sync:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Set_sync: base version is a no-op.
  !!
  !
  subroutine i_base_mlv_set_sync(x)
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x

  end subroutine i_base_mlv_set_sync

  !
  !> Function  base_mlv_is_dev:
  !! \memberof  psb_i_base_multivect_type
  !! \brief Is  vector on external device    .
  !!
  !
  function i_base_mlv_is_dev(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    logical  :: res

    res = .false.
  end function i_base_mlv_is_dev

  !
  !> Function  base_mlv_is_host
  !! \memberof  psb_i_base_multivect_type
  !! \brief Is  vector on standard memory    .
  !!
  !
  function i_base_mlv_is_host(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function i_base_mlv_is_host

  !
  !> Function  base_mlv_is_sync
  !! \memberof  psb_i_base_multivect_type
  !! \brief Is  vector on sync               .
  !!
  !
  function i_base_mlv_is_sync(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function i_base_mlv_is_sync

  !
  ! Size info.
  !
  !
  !> Function  base_mlv_get_nrows
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Number of entries
  !!
  !
  function i_base_mlv_get_nrows(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    integer(psb_ipk_) :: res

    res = 0
    if (allocated(x%v)) res = size(x%v,1)

  end function i_base_mlv_get_nrows

  function i_base_mlv_get_ncols(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    integer(psb_ipk_) :: res

    res = 0
    if (allocated(x%v)) res = size(x%v,2)

  end function i_base_mlv_get_ncols

  !
  !> Function  base_mlv_get_sizeof
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Size in bytesa
  !!
  !
  function i_base_mlv_sizeof(x) result(res)
    implicit none
    class(psb_i_base_multivect_type), intent(in) :: x
    integer(psb_long_int_k_) :: res

    ! Force 8-byte integers.
    res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()

  end function i_base_mlv_sizeof

  !
  !> Function  base_mlv_get_fmt
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Format
  !!
  !
  function i_base_mlv_get_fmt() result(res)
    implicit none
    character(len=5) :: res
    res = 'BASE'
  end function i_base_mlv_get_fmt

  !
  !
  !
  !> Function  base_mlv_get_vect
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Extract a copy of the contents
  !!
  !
  function  i_base_mlv_get_vect(x) result(res)
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), allocatable                 :: res(:,:)
    integer(psb_ipk_) :: info,m,n
    m = x%get_nrows()
    n = x%get_ncols()
    if (.not.allocated(x%v)) return
    call x%sync()
    allocate(res(m,n),stat=info)
    if (info /= 0) then
      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
      return
    end if
    res(1:m,1:n) = x%v(1:m,1:n)
  end function i_base_mlv_get_vect

  !
  ! Reset all values
  !
  !
  !> Function  base_mlv_set_scal
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Set all entries
  !! \param val   The value to set
  !!
  subroutine i_base_mlv_set_scal(x,val)
    class(psb_i_base_multivect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in) :: val

    integer(psb_ipk_) :: info
    x%v = val

  end subroutine i_base_mlv_set_scal

  !
  !> Function  base_mlv_set_vect
  !! \memberof  psb_i_base_multivect_type
  !! \brief  Set all entries
  !! \param val(:)  The vector to be copied in
  !!
  subroutine i_base_mlv_set_vect(x,val)
    class(psb_i_base_multivect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in) :: val(:,:)
    integer(psb_ipk_) :: nr
    integer(psb_ipk_) :: info

    if (allocated(x%v)) then
      nr = min(size(x%v,1),size(val,1))
      nc = min(size(x%v,2),size(val,2))

      x%v(1:nr,1:nc) = val(1:nr,1:nc)
    else
      x%v = val
    end if

  end subroutine i_base_mlv_set_vect

  function i_base_mlv_use_buffer() result(res)
    logical :: res

    res = .true.
  end function i_base_mlv_use_buffer

  subroutine i_base_mlv_new_buffer(n,x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)              :: n
    integer(psb_ipk_), intent(out)             :: info

    integer(psb_ipk_)               :: nc
    nc = x%get_ncols()
    call psb_realloc(n*nc,x%combuf,info)
  end subroutine i_base_mlv_new_buffer

  subroutine i_base_mlv_new_comid(n,x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)              :: n
    integer(psb_ipk_), intent(out)             :: info

    call psb_realloc(n,2_psb_ipk_,x%comid,info)
  end subroutine i_base_mlv_new_comid

  subroutine i_base_mlv_maybe_free_buffer(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    info = 0
    if (psb_get_maybe_free_buffer())&
         &  call x%free_buffer(info)

  end subroutine i_base_mlv_maybe_free_buffer

  subroutine i_base_mlv_free_buffer(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    if (allocated(x%combuf)) &
         &  deallocate(x%combuf,stat=info)
  end subroutine i_base_mlv_free_buffer

  subroutine i_base_mlv_free_comid(x,info)
    use psb_realloc_mod
    implicit none
    class(psb_i_base_multivect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info

    if (allocated(x%comid)) &
         &  deallocate(x%comid,stat=info)
  end subroutine i_base_mlv_free_comid

  !
  ! Gather: Y = beta * Y + alpha * X(IDX(:))
  !
  !
  !> Function  base_mlv_gthab
  !! \memberof  psb_i_base_multivect_type
  !! \brief gather into an array
  !!    Y = beta * Y + alpha * X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  !! \param alpha
  !! \param beta
  subroutine i_base_mlv_gthab(n,idx,alpha,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) :: alpha, beta, y(:)
    class(psb_i_base_multivect_type) :: x
    integer(psb_ipk_) :: nc

    if (x%is_dev()) call x%sync()
    if (.not.allocated(x%v)) then
      return
    end if
    nc = psb_size(x%v,2_psb_ipk_)
    call psi_gth(n,nc,idx,alpha,x%v,beta,y)

  end subroutine i_base_mlv_gthab
  !
  ! shortcut alpha=1 beta=0
  !
  !> Function  base_mlv_gthzv
  !! \memberof  psb_i_base_multivect_type
  !! \brief gather into an array special alpha=1 beta=0
  !!    Y = X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  subroutine i_base_mlv_gthzv_x(i,n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i,n
    class(psb_i_base_vect_type) :: idx
    integer(psb_ipk_) ::  y(:)
    class(psb_i_base_multivect_type) :: x

    if (x%is_dev()) call x%sync()
    call x%gth(n,idx%v(i:),y)

  end subroutine i_base_mlv_gthzv_x

  !
  ! shortcut alpha=1 beta=0
  !
  !> Function  base_mlv_gthzv
  !! \memberof  psb_i_base_multivect_type
  !! \brief gather into an array special alpha=1 beta=0
  !!    Y = X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  subroutine i_base_mlv_gthzv(n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) ::  y(:)
    class(psb_i_base_multivect_type) :: x
    integer(psb_ipk_) :: nc

    if (x%is_dev()) call x%sync()
    if (.not.allocated(x%v)) then
      return
    end if
    nc = psb_size(x%v,2_psb_ipk_)

    call psi_gth(n,nc,idx,x%v,y)

  end subroutine i_base_mlv_gthzv
  !
  ! shortcut alpha=1 beta=0
  !
  !> Function  base_mlv_gthzv
  !! \memberof  psb_i_base_multivect_type
  !! \brief gather into an array special alpha=1 beta=0
  !!    Y = X(IDX(:))
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  subroutine i_base_mlv_gthzm(n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) ::  y(:,:)
    class(psb_i_base_multivect_type) :: x
    integer(psb_ipk_) :: nc

    if (x%is_dev()) call x%sync()
    if (.not.allocated(x%v)) then
      return
    end if
    nc = psb_size(x%v,2_psb_ipk_)

    call psi_gth(n,nc,idx,x%v,y)

  end subroutine i_base_mlv_gthzm

  !
  ! New comm internals impl.
  !
  subroutine i_base_mlv_gthzbuf(i,ixb,n,idx,x)
    use psi_serial_mod
    integer(psb_ipk_) :: i, ixb, n
    class(psb_i_base_vect_type) :: idx
    class(psb_i_base_multivect_type) :: x
    integer(psb_ipk_) :: nc

    if (.not.allocated(x%combuf)) then
      call psb_errpush(psb_err_alloc_dealloc_,'gthzbuf')
      return
    end if
    if (idx%is_dev()) call idx%sync()
    if (x%is_dev()) call x%sync()
    nc = x%get_ncols()
    call x%gth(n,idx%v(i:),x%combuf(ixb:))

  end subroutine i_base_mlv_gthzbuf

  !
  ! Scatter:
  ! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
  !
  !
  !> Function  base_mlv_sctb
  !! \memberof  psb_i_base_multivect_type
  !! \brief scatter into a class(base_mlv_vect)
  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
  !! \param n  how many entries to consider
  !! \param idx(:) indices
  !! \param beta
  !! \param x(:)
  subroutine i_base_mlv_sctb(n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) :: beta, x(:)
    class(psb_i_base_multivect_type) :: y
    integer(psb_ipk_) :: nc

    if (y%is_dev()) call y%sync()
    nc = psb_size(y%v,2_psb_ipk_)
    call psi_sct(n,nc,idx,x,beta,y%v)
    call y%set_host()

  end subroutine i_base_mlv_sctb

  subroutine i_base_mlv_sctbr2(n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    integer(psb_ipk_) :: beta, x(:,:)
    class(psb_i_base_multivect_type) :: y
    integer(psb_ipk_) :: nc

    if (y%is_dev()) call y%sync()
    nc = y%get_ncols()
    call psi_sct(n,nc,idx,x,beta,y%v)
    call y%set_host()

  end subroutine i_base_mlv_sctbr2

  subroutine i_base_mlv_sctb_x(i,n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i, n
    class(psb_i_base_vect_type) :: idx
    integer( psb_ipk_) :: beta, x(:)
    class(psb_i_base_multivect_type) :: y

    call y%sct(n,idx%v(i:),x,beta)

  end subroutine i_base_mlv_sctb_x

  subroutine i_base_mlv_sctb_buf(i,iyb,n,idx,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: i, iyb, n
    class(psb_i_base_vect_type) :: idx
    integer(psb_ipk_) :: beta
    class(psb_i_base_multivect_type) :: y
    integer(psb_ipk_) :: nc

    if (.not.allocated(y%combuf)) then
      call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
      return
    end if
    if (y%is_dev()) call y%sync()
    if (idx%is_dev()) call idx%sync()
    nc = y%get_ncols()
    call y%sct(n,idx%v(i:),y%combuf(iyb:),beta)
    call y%set_host()

  end subroutine i_base_mlv_sctb_buf

  !
  !> Function  base_device_wait:
  !! \memberof  psb_i_base_vect_type
  !! \brief device_wait: base version is a no-op.
  !!
  !
  subroutine i_base_mlv_device_wait()
    implicit none

  end subroutine i_base_mlv_device_wait

end module psb_i_base_multivect_mod


text:
subroutine case_DoTimeStep(a)
   use typre
   use Mod_masterVariables
   use Mod_GeneralCase
   implicit none
   class(GeneralCase), target :: a

   type (masterVariables), pointer :: m => NULL()
   m => a%caseVars%masterVars

   !Iterate the case
   m%kfl_gocou = 1
   m%iiter = 1

   !This is only local coupling, not case coupling
   !-----local to the container------
   local_coupling: do while (m%kfl_gocou==1)

       !Mesh Projections for FMALE
       call a%MeshProjections
       call a%Doiter
       call a%LocalCouplingConvergence

   end do local_coupling

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
  subroutine radar_simulator( &
    hp, &
    nprof,ngate, &
    undef, &
    hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
    p_matrix,t_matrix,rh_matrix, &
    Ze_non,Ze_ray,a_to_vol,g_to_vol,dBZe, &
    g_to_vol_in,g_to_vol_out)

  use m_mrgrnk
  use array_lib
  use math_lib
  use optics_lib
  use radar_simulator_types
  use scale_LUTs_io
  implicit none

! Purpose:
!
!   Simulates a vertical profile of radar reflectivity
!   Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand.
!   but has been substantially modified since that time by
!   Laura Fowler and Roger Marchand (see modifications below).
!
! Inputs:
!
!   [hp]              structure that defines hydrometeor types and other radar properties
!
!   [nprof]           number of hydrometeor profiles
!   [ngate]           number of vertical layers
!
!   [undef]           missing data value
!   (The following 5 arrays must be in order from closest to the radar
!    to farthest...)
!
!   [hgt_matrix]      height of hydrometeors (km)
!   [p_matrix]        pressure profile (hPa)
!   [t_matrix]        temperature profile (K)
!   [rh_matrix]       relative humidity profile (%) -- only needed if gaseous aborption calculated.
!
!   [hm_matrix]       table of hydrometeor mixing rations (g/kg)
!   [re_matrix]       table of hydrometeor effective radii.  0 ==> use defaults. (units=microns)
!   [Np_matrix]	      table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
!
! Outputs:
!
!   [Ze_non]          radar reflectivity without attenuation (dBZ)
!   [Ze_ray]          Rayleigh reflectivity (dBZ)
!   [h_atten_to_vol]  attenuation by hydromets, radar to vol (dB)
!   [g_atten_to_vol]  gaseous atteunation, radar to vol (dB)
!   [dBZe]            effective radar reflectivity factor (dBZ)
!
! Optional:
!   [g_to_vol_in]     integrated atten due to gases, r>v (dB).
!                     If present then is used as gaseous absorption, independently of the
!                     value in use_gas_abs
!   [g_to_vol_out]    integrated atten due to gases, r>v (dB).
!                     If present then gaseous absorption for each profile is returned here.
!
! Created:
!   11/28/2005  John Haynes (haynes@atmos.colostate.edu)
!
! Modified:
!   09/2006  placed into subroutine form (Roger Marchand,JMH)
!   08/2007  added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
!   01/2008  'Do while' to determine if hydrometeor(s) present in volume
!             changed for vectorization purposes (A. Bodas-Salcedo)
!
!   07/2010  V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT)
!  ... All hydrometeor and radar simulator properties now included in hp structure
!  ... hp structure should be initialized by call to radar_simulator_init prior
!  ... to calling this subroutine.
!     Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added
!  ... Changes implement by Roj Marchand following work by Laura Fowler
!
!   10/2011  Modified ngate loop to go in either direction depending on flag
!     hp%radar_at_layer_one.  This affects the direction in which attenuation is summed.
!
!     Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple
!     summation. (Roger Marchand)
!
!
! ----- INPUTS -----

  logical, parameter  ::  DO_LUT_TEST = .false.
  logical, parameter  ::  DO_NP_TEST = .false.

  type(class_param), intent(inout) :: hp

  integer, intent(in) ::  nprof,ngate

  real undef
  real*8, dimension(nprof,ngate), intent(in) :: &
    hgt_matrix, p_matrix,t_matrix,rh_matrix

  real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
  real*8, dimension(hp%nhclass,nprof,ngate), intent(in)    :: Np_matrix

! ----- OUTPUTS -----
  real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
       g_to_vol,dBZe,a_to_vol

! ----- OPTIONAL -----
  real*8, optional, dimension(nprof,ngate) :: &
  g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
                           ! the same gaseous absorption in different calls. Optional to allow compatibility
                           ! with original version. A. Bodas April 2008.
!  real*8, dimension(nprof,ngate) :: kr_matrix

! ----- INTERNAL -----

  real, parameter :: one_third = 1.0/3.0
  real*8 :: t_kelvin
  integer :: &
  phase, & ! 0=liquid, 1=ice
  ns       ! number of discrete drop sizes

  logical :: hydro	! true=hydrometeor in vol, false=none
  real*8 :: &
  rho_a, &   ! air density (kg m^-3)
  gases      ! function: 2-way gas atten (dB/km)

  real*8, dimension(:), allocatable :: &
  Di, Deq, &   ! discrete drop sizes (um)
  Ni, &        ! discrete concentrations (cm^-3 um^-1)
  rhoi         ! discrete densities (kg m^-3)

  real*8, dimension(nprof, ngate) :: &
  z_vol, &      ! effective reflectivity factor (mm^6/m^3)
  z_ray, &                      ! reflectivity factor, Rayleigh only (mm^6/m^3)
  kr_vol, &     ! attenuation coefficient hydro (dB/km)
  g_vol         ! attenuation coefficient gases (dB/km)

  integer,parameter :: KR8 = selected_real_kind(15,300)
  real*8, parameter :: xx = -1.0_KR8
  real*8,  dimension(:), allocatable :: xxa
  real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
  real*8 :: half_a_atten_current,half_a_atten_above
  real*8 :: half_g_atten_current,half_g_atten_above
  integer*4 :: tp, i, j, k, pr, itt, iff

  real*8    step,base, Np
  integer*4 iRe_type,n,max_bin

  integer   start_gate,end_gate,d_gate

  logical :: g_to_vol_in_present, g_to_vol_out_present

  ! Logicals to avoid calling present within the loops
  g_to_vol_in_present  = present(g_to_vol_in)
  g_to_vol_out_present = present(g_to_vol_out)

  !
  ! load scaling matricies from disk -- but only the first time this subroutine is called
  !
  if(hp%load_scale_LUTs) then
    call load_scale_LUTs(hp)
    hp%load_scale_LUTs=.false.
    hp%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
  endif

  pi = acos(-1.0)

!   ----- Initialisation -----
  g_to_vol = 0.0
  a_to_vol = 0.0
  z_vol    = 0.0
  z_ray    = 0.0
  kr_vol   = 0.0

!   // loop over each range gate (ngate) ... starting with layer closest to the radar !
  if(hp%radar_at_layer_one) then
    start_gate=1
    end_gate=ngate
    d_gate=1
  else
    start_gate=ngate
    end_gate=1
    d_gate=-1
  endif
  do k=start_gate,end_gate,d_gate
  ! // loop over each profile (nprof)
    do pr=1,nprof
      t_kelvin = t_matrix(pr,k)
!     :: determine if hydrometeor(s) present in volume
      hydro = .false.
      do j=1,hp%nhclass
        if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
          hydro = .true.
          exit
        endif
      enddo

!     :: if there is hydrometeor in the volume
      if (hydro) then

        rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin))
!       :: loop over hydrometeor type
        do tp=1,hp%nhclass
          if (hm_matrix(tp,pr,k) <= 1E-12) cycle
          phase = hp%phase(tp)
          if (phase==0) then
            itt = infind(hp%mt_ttl,t_kelvin)
          else
            itt = infind(hp%mt_tti,t_kelvin)
          endif
          if (re_matrix(tp,pr,k).eq.0) then
            call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, &
              hp%dtype(tp),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
              hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),Re)
            re_matrix(tp,pr,k)=Re
          else
            if (Np_matrix(tp,pr,k)>0) then
               print *, 'Warning: Re and Np set for the same ', &
                        'volume & hydrometeor type.  Np is being ignored.'
            endif
            Re = re_matrix(tp,pr,k)
          endif

          iRe_type=1
          if(Re.gt.0) then
            ! determine index in to scale LUT
            !
            ! distance between Re points (defined by "base" and "step") for
            ! each interval of size Re_BIN_LENGTH
            ! Integer asignment, avoids calling floor intrinsic
            n=Re/Re_BIN_LENGTH
            if (n>=Re_MAX_BIN) n=Re_MAX_BIN-1
            step=hp%step_list(n+1)
            base=hp%base_list(n+1)
            iRe_type=Re/step
            if (iRe_type.lt.1) iRe_type=1

            Re=step*(iRe_type+0.5)	! set value of Re to closest value allowed in LUT.
            iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)

            ! make sure iRe_type is within bounds
            if (iRe_type.ge.nRe_types) then
!               write(*,*) 'Warning: size of Re exceed value permitted ', &
!                    'in Look-Up Table (LUT).  Will calculate. '
               ! no scaling allowed
               iRe_type=nRe_types
               hp%Z_scale_flag(tp,itt,iRe_type)=.false.
            else
               ! set value in re_matrix to closest values in LUT
              if (.not. DO_LUT_TEST) re_matrix(tp,pr,k)=Re
            endif
          endif
          ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
          ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
          if( (.not. hp%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST)  then
            ! :: create a distribution of hydrometeors within volume
            select case(hp%dtype(tp))
              case(4)
                ns = 1
                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
                Di = hp%p1(tp)
                Ni = 0.
              case default
                ns = nd   ! constant defined in radar_simulator_types.f90
                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
                Di = hp%D
                Ni = 0.
            end select
            call dsd(hm_matrix(tp,pr,k),re_matrix(tp,pr,k),Np_matrix(tp,pr,k), &
                     Di,Ni,ns,hp%dtype(tp),rho_a,t_kelvin, &
                     hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
                     hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp))

            ! calculate particle density
            if (phase == 1) then
              if (hp%rho(tp) < 0) then
                ! Use equivalent volume spheres.
                hp%rho_eff(tp,1:ns,iRe_type) = 917  				! solid ice == equivalent volume approach
                Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6
                ! alternative is to comment out above two lines and use the following block
                ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density
                !
                ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3)   !MG Mie approach

                ! as the particle size gets small it is possible that the mass to size relationship of
                ! (given by power law in hclass.data) can produce impossible results
                ! where the mass is larger than a solid sphere of ice.
                ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
                ! do i=1,ns
                ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
                ! hp%rho_eff(tp,i,iRe_type) = 917
                ! endif
                ! enddo
              else
                ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
                hp%rho_eff(tp,1:ns,iRe_type) = 917
                Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0))
                ! alternative ... coment out above two lines and use the following for MG-Mie
                ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp)   !MG Mie approach
              endif
            else
              ! I assume here that water phase droplets are spheres.
              ! hp%rho should be ~ 1000  or hp%apm=524 .and. hp%bpm=3
              Deq = Di
            endif

            ! calculate effective reflectivity factor of volume
            xxa = -9.9
            rhoi = hp%rho_eff(tp,1:ns,iRe_type)
            call zeff(hp%freq,Deq,Ni,ns,hp%k2,t_kelvin,phase,hp%do_ray, &
                      ze,zr,kr,xxa,xxa,rhoi)

            ! test code ... compare Np value input to routine with sum of DSD
            ! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation
            ! not just the DSD representation given by Ni
            if(Np_matrix(tp,pr,k)>0 .and. DO_NP_TEST ) then
              Np = path_integral(Ni,Di,1,ns-1)/rho_a*1E6
              ! Note: Representation is not great or small Re < 2
              if( (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)>0.1 ) then
                write(*,*) 'Error: Np input does not match sum(N)'
                write(*,*) tp,pr,k,Re,Ni(1),Ni(ns),10*log10(ze)
                write(*,*) Np_matrix(tp,pr,k),Np,(Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)
                write(*,*)
              endif
            endif

            deallocate(Di,Ni,rhoi,xxa,Deq)

            ! LUT test code
            ! This segment of code compares full calculation to scaling result
            if ( hp%Z_scale_flag(tp,itt,iRe_type) .and. DO_LUT_TEST )  then
              scale_factor=rho_a*hm_matrix(tp,pr,k)
              ! if more than 2 dBZe difference print error message/parameters.
              if ( abs(10*log10(ze) - 10*log10(hp%Ze_scaled(tp,itt,iRe_type) * &
                   scale_factor)) > 2 ) then
                write(*,*) 'Roj Error: ',tp,itt,iRe_type,hp%Z_scale_flag(tp,itt,iRe_type),n,step,base
                write(*,*) 10*log10(ze),10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor)
                write(*,*) hp%Ze_scaled(tp,itt,iRe_type),scale_factor
                write(*,*) re_matrix(tp,pr,k),Re
                write(*,*)
              endif
            endif

          else ! can use z scaling
            scale_factor=rho_a*hm_matrix(tp,pr,k)
            zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor
            ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
            kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor
          endif  ! end z_scaling

          kr_vol(pr,k) = kr_vol(pr,k) + kr
          z_vol(pr,k)  = z_vol(pr,k)  + ze
          z_ray(pr,k)  = z_ray(pr,k)  + zr

          ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
          if ( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then
            if (iRe_type>1) then
              scale_factor=rho_a*hm_matrix(tp,pr,k)
              hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
              hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
              hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
              hp%Z_scale_flag(tp,itt,iRe_type) = .true.
              hp%Z_scale_added_flag(tp,itt,iRe_type)=.true.
            endif
          endif

        enddo	! end loop of tp (hydrometeor type)

      else
!     :: volume is hydrometeor-free
        kr_vol(pr,k) = 0
        z_vol(pr,k)  = undef
        z_ray(pr,k)  = undef
      endif

      !     :: attenuation due to hydrometeors between radar and volume
      !
      ! NOTE old scheme integrates attenuation only for the layers ABOVE
      ! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
      ! which may be a problem.   ROJ
      ! in the new scheme I assign half the attenuation to the current layer
      if(d_gate==1) then
        ! dheight calcuations assumes hgt_matrix points are the cell mid-points.
        if (k>2) then
          ! add to previous value to half of above layer + half of current layer
          a_to_vol(pr,k)=  a_to_vol(pr,k-1) + &
             (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
        else
          a_to_vol(pr,k)=  kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
        endif
      else   ! d_gate==-1
        if(k<ngate) then
          ! add to previous value half of above layer + half of current layer
          a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
              (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
        else
          a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
        endif
      endif

      !     :: attenuation due to gaseous absorption between radar and volume
      if (g_to_vol_in_present) then
        g_to_vol(pr,k) = g_to_vol_in(pr,k)
      else
        if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .and. (pr == 1)) ) then
          g_vol(pr,k) = gases(p_matrix(pr,k),t_kelvin,rh_matrix(pr,k),hp%freq)
          if (d_gate==1) then
            if (k>1) then
              ! add to previous value to half of above layer + half of current layer
              g_to_vol(pr,k) =  g_to_vol(pr,k-1) + &
                  0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
            else
              g_to_vol(pr,k)=  0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
            endif
          else   ! d_gate==-1
            if (k<ngate) then
              ! add to previous value to half of above layer + half of current layer
              g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
                 0.5*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
            else
              g_to_vol(pr,k)= 0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
            endif
          endif
        elseif(hp%use_gas_abs == 2) then
          ! using value calculated for the first column
          g_to_vol(pr,k) = g_to_vol(1,k)
        elseif (hp%use_gas_abs == 0) then
          g_to_vol(pr,k) = 0
        endif
      endif

      ! Compute Rayleigh reflectivity, and full, attenuated reflectivity
      if ((hp%do_ray == 1) .and. (z_ray(pr,k) > 0)) then
        Ze_ray(pr,k) = 10*log10(z_ray(pr,k))
      else
        Ze_ray(pr,k) = undef
      endif
      if (z_vol(pr,k) > 0) then
        Ze_non(pr,k) = 10*log10(z_vol(pr,k))
        dBZe(pr,k) = Ze_non(pr,k)-a_to_vol(pr,k)-g_to_vol(pr,k)
      else
        dBZe(pr,k) = undef
        Ze_non(pr,k) = undef
      endif

    enddo   ! end loop over pr (profile)

  enddo ! end loop of k (range gate)

  ! Output array with gaseous absorption
  if (g_to_vol_out_present) g_to_vol_out = g_to_vol

  ! save any updates made
  if (hp%update_scale_LUTs) call save_scale_LUTs(hp)

  end subroutine radar_simulator

max_stars_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_issues_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_forks_repo_licenses: zlib-acknowledgement RSA-MD FTL


text:
      subroutine higgsw(br)
      implicit none
      include 'constants.f'
      include 'ewcouple.f'
      include 'masses.f'
      double precision wff,mfsq,br
c-----approximate form for the width of the standard model higgs
c-----valid for low masses
      wff(mfsq)=dsqrt(2d0)/8d0/pi*gf*hmass*mfsq
     & *(1d0-4d0*mfsq/hmass**2)**1.5d0

      hwidth=3d0*(wff(mbsq)+wff(mcsq))+wff(mtausq)
      write(6,*) 'hmass,hwidth',hmass,hwidth
      write(6,*) 'mtausq,mcsq,mbsq',mtausq,mcsq,mbsq
      write(6,*)
      br=3d0*wff(mbsq)/hwidth
      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen