! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module cover_module
! ------------------------------------------------------------------------------
! Code area 8: indexing
! -------------------------------------------------------------------------------

!!****h* Conquest/cover_module
!!  NAME
!!   cover_module
!!  PURPOSE
!!   This module holds variables and routines related to covering
!!   sets.  A covering set is (at the moment) the smallest orthorhombic
!!   collection of small groups that will fully enclose a primary set
!!   out to a specified radius.  It's used for searching for neighbours
!!   and communicating details about small groups between processors.
!!
!!   More details (exhaustive details !) can be found in the Conquest
!!   notes "Matrix multiplication in Conquest: Practical details" 
!!   which are in the Conquest documentation repository (stored
!!   in the directory TechnicalMatMult/)
!!  USES
!!   GenComms, basic_types, datatypes, global_module
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   19/02/00 by D.R.Bowler
!!    Implementing changes for transposes (adding ncover_rem and
!!    routines to find iprim_part)
!!   03/03/00 by DRB
!!    Finishing up transposes
!!   29/03/00 by DRB 
!!    Creating a generalised cover type and altering make_gcs to use it
!!   11/04/00 by DRB
!!    Correcting cover type
!!   19/04/00 by DRB
!!    Adding general group and primary set types to make GCS totally 
!!    general
!!   20/04/00 by DRB
!!    Removed type declaration and placed in basic_types module
!!    Continued generalisation of make_cs
!!    Removed GOTOs from indexx (replaced with do...while)
!!   21/06/2001 dave
!!    Added ROBODoc headers and removed stop statements and MPI_Abort
!!    statements throughout
!!   19/04/2002 drb
!!    Corrected reference to notes and added USE field for whole module
!!   17/06/2002 dave
!!    Added RCS Id tag and improved headers a little
!!   2008/02/04 17:12 dave
!!    Changes for output to file not stdoutx
!!   2008/05/16 ast
!!    Added some timers
!!   2014/09/15 18:30 lat
!!    fixed call start/stop_timer to timer_module (not timer_stdlocks_module !)
!!   2016/01/28 16:42 dave
!!    Changed ewald_CS to ion_ion_CS
!!   2019/11/04 15:12 dave
!!    Replace indexx from Numerical Recipes with generic heapsort
!!  SOURCE
!!
module cover_module

  ! Module use
  use datatypes
  use basic_types  
  use global_module,          only: io_lun
  use timer_module,           only: start_timer, stop_timer
  use timer_stdclocks_module, only: tmr_std_indexing, tmr_std_allocation


  implicit none

  save

  ! DCS = domain CS; BCS = bundle CS
  type(cover_set) :: DCS_blocks
  type(cover_set) :: DCS_parts
  type(cover_set) :: BCS_blocks
  type(cover_set) :: BCS_parts
  type(cover_set) :: ion_ion_CS
  type(cover_set) :: D2_CS ! for DFT-D2

!!***
contains

!!****f* cover_module/make_cs *
!!
!!  NAME 
!!   make_cs - makes a covering set
!!  USAGE
!! 
!!  PURPOSE
!!   sbrt make_cs: create a covering set out of groups given by the 
!!   group_set variable, for a cutoff gr_rcut around the primary set 
!!   defined by the widths (nw_x,y,z) and left-spans (nlx,y,z).
!!   N.B. The widths and left spans must have been converted into units
!!   of the CS groups !
!!   The origin is also passed; this is a problem is we're making a
!!   CS for one group type out of the other.  It MUST be the group 
!!   nearest the origin partition of the primary set.
!!   The results are put into the derived type set.
!!   The logical variable members decides whether or not the 
!!   x,y,zcover variables are created (NOT for the cross-group CS)
!!
!!   The notes above are a little out of date - the widths and left-spans
!!   are now part of the primary set type, and conversion is done by 
!!   convert_primary (below)
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added some timers
!!   2013/07/01 M.Arita
!!    Added allocation of ig_cover
!!   2013/08/21 M.Arita
!!    Initialised and calculated ig_cover for loading matrices
!!  SOURCE
!!
  subroutine make_cs(myid, gr_rcut, set, groups, prim, mx_mcell, &
                     x_mem_cell, y_mem_cell, z_mem_cell)

    ! Module usage
    use global_module
    use basic_types
    use GenComms, ONLY: cq_abort, inode, ionode
    use memory_module, ONLY: reg_alloc_mem, reg_dealloc_mem, type_int,&
         type_dbl
    use functions, ONLY: heapsort_integer_index

    implicit none

    ! Passed variables
    type(cover_set) :: set
    type(group_set), target   :: groups
    type(primary_set), target :: prim
    integer, intent(in) :: myid
    integer, OPTIONAL :: mx_mcell
    real(double), dimension(:), intent(IN), OPTIONAL :: x_mem_cell
    real(double), dimension(:), intent(IN), OPTIONAL :: y_mem_cell
    real(double), dimension(:), intent(IN), OPTIONAL :: z_mem_cell
    real(double) :: gr_rcut ! Cutoff radius for CS

    ! Local variables
    ! Number of periodic images of a group
    integer :: nrepx(groups%mx_gedge)
    integer :: nrepy(groups%mx_gedge)
    integer :: nrepz(groups%mx_gedge)
    ! x,y,z numbering of CS groups (for CC labels)
    integer, allocatable, dimension(:)  :: nx_in_cover
    integer, allocatable, dimension(:)  :: ny_in_cover
    integer, allocatable, dimension(:)  :: nz_in_cover
    ! Variables for irreducible CS
    integer :: ind_min(groups%mx_gcell)
    integer :: ngcx_min(groups%mx_gcell)
    integer :: ngcy_min(groups%mx_gcell)
    integer :: ngcz_min(groups%mx_gcell)
    integer :: min_sort(groups%mx_gcell)
    ! Converted units of primary cell size
    integer :: nx_o,ny_o,nz_o
    ! General scalars
    integer :: ng_in_cell,nmodx,nmody,nmodz
    integer :: noccx,nremx,minx,ngcx,noccy,nremy,miny,ngcy,noccz,nremz
    integer :: minz,ngcz,ng_in_min,ind,nqx,nqy,nqz,ind_qart,ino,ind_cover
    integer :: nrx,nry,nrz,nsx,nsy,nsz,ni,nnd,stat
    integer :: nm_in_cover
    real(double) :: dcellx,dcelly,dcellz,xadd,yadd,zadd
    logical :: members

    call start_timer(tmr_std_indexing)
    ! This determines whether or not we find information about all 
    ! the members of the groups in the covering set
    members = (PRESENT(mx_mcell).AND.PRESENT(x_mem_cell).AND.&
         PRESENT(y_mem_cell).AND.PRESENT(z_mem_cell))  
    nnd = myid+1
    ! Start the derived type - assign useful variables
    set%rcut = gr_rcut     ! Cutoff around primary set
    set%groups => groups   ! Groups out of which to build CS
    set%prim => prim       ! Primary set around which to build CS
    ! Find the origin, width and left span of the CS
    call convert_primary(groups,prim,set,nx_o,ny_o,nz_o)
    ! This gives the total number of groups
    set%ng_cover=set%ncoverx*set%ncovery*set%ncoverz
    set%mx_gcover = set%ng_cover
    call allocate_cs(set,numprocs,members)
    call start_timer(tmr_std_allocation)
    allocate(nx_in_cover(set%ng_cover),ny_in_cover(set%ng_cover),nz_in_cover(set%ng_cover),STAT=stat)
    if(stat/=0) call cq_abort("Error allocating nx_in_cover: ",set%ng_cover,stat)
    call reg_alloc_mem(area_index, 3*set%ng_cover,type_int)
    call stop_timer(tmr_std_allocation)
    !if(set%ng_cover>set%mx_gcover) then
    !   call cq_abort('make_cs: too many groups in GCS',&
    !        set%ng_cover,set%mx_gcover)
    !endif
    ! Conversion factors from unit cell lengths->groups
    dcellx=rcellx/real(groups%ngcellx,double)
    dcelly=rcelly/real(groups%ngcelly,double)
    dcellz=rcellz/real(groups%ngcellz,double)
    ! Fully explained in notes mentioned above
    nmodx=((groups%ngcellx+set%nspanlx-1)/groups%ngcellx)*groups%ngcellx
    nmody=((groups%ngcelly+set%nspanly-1)/groups%ngcelly)*groups%ngcelly
    nmodz=((groups%ngcellz+set%nspanlz-1)/groups%ngcellz)*groups%ngcellz
    ! create apparatus for periodic boundary conditions ---------------
    ! ... x-direction
    noccx=set%ncoverx/groups%ngcellx
    nremx=set%ncoverx-noccx*groups%ngcellx
    minx=min(set%ncoverx,groups%ngcellx)
    if(minx>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in x-edge',minx)
    endif
    do ngcx=1,minx
       if(ngcx<=nremx) then
          nrepx(ngcx)=noccx+1
       else
          nrepx(ngcx)=noccx
       endif
    enddo
    ! ... y-direction
    noccy=set%ncovery/groups%ngcelly
    nremy=set%ncovery-noccy*groups%ngcelly
    miny=min(set%ncovery,groups%ngcelly)
    if(miny>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in y-edge',miny)
    endif
    do ngcy=1,miny
       if(ngcy<=nremy) then
          nrepy(ngcy)=noccy+1
       else
          nrepy(ngcy)=noccy
       endif
    enddo
    ! ... z-direction
    noccz=set%ncoverz/groups%ngcellz
    nremz=set%ncoverz-noccz*groups%ngcellz
    minz=min(set%ncoverz,groups%ngcellz)
    if(minz>groups%mx_gedge) then
       call cq_abort('make_cs: too many groups in z-edge',minz)
    endif
    do ngcz=1,minz
       if(ngcz<=nremz) then
          nrepz(ngcz)=noccz+1
       else
          nrepz(ngcz)=noccz
       endif
    enddo
    ! go over groups in GCS periodic-irreducible set, calculating
    ! simulation-cell (node-order, home-start) label of each 
    ng_in_cell = groups%ngcellx*groups%ngcelly*groups%ngcellz
    ng_in_min = minx*miny*minz
    ind=0
    do ngcx=1,minx
       do ngcy=1,miny
          do ngcz=1,minz
             ind=ind+1
             nqx=1+mod(nx_o+ngcx-set%nspanlx-2+nmodx,groups%ngcellx)
             nqy=1+mod(ny_o+ngcy-set%nspanly-2+nmody,groups%ngcelly)
             nqz=1+mod(nz_o+ngcz-set%nspanlz-2+nmodz,groups%ngcellz)
             ind_qart=(nqx-1)*groups%ngcelly*groups%ngcellz+(nqy-1)*groups%ngcellz+nqz
             ino=groups%inv_ngnode(ind_qart)
             ind_min(ind)=1+mod(ino-groups%inode_beg(nnd)+ ng_in_cell,ng_in_cell)
             ngcx_min(ind)=ngcx
             ngcy_min(ind)=ngcy
             ngcz_min(ind)=ngcz
          enddo
       enddo
    enddo
    ! sort minimum CS by nodes 
    call heapsort_integer_index(ng_in_min,ind_min,min_sort)

    ! go over all GCS groups in node-periodic-grouped order 
    ind_cover=0
    do ind=1,ng_in_min
       ngcx=ngcx_min(min_sort(ind))
       ngcy=ngcy_min(min_sort(ind))
       ngcz=ngcz_min(min_sort(ind))
       do nrx=1,nrepx(ngcx)
          do nry=1,nrepy(ngcy)
             do nrz=1,nrepz(ngcz)
                ind_cover=ind_cover+1
                nx_in_cover(ind_cover)=ngcx-1-set%nspanlx+(nrx-1)*groups%ngcellx
                ny_in_cover(ind_cover)=ngcy-1-set%nspanly+(nry-1)*groups%ngcelly
                nz_in_cover(ind_cover)=ngcz-1-set%nspanlz+(nrz-1)*groups%ngcellz
             enddo
          enddo
       enddo
    enddo
    ! make covering set for primary set 
    !if(members) set%icover_ibeg(1)=1
    if(members) nm_in_cover = 0
    do ind_cover=1,set%ng_cover
       nsx=nx_in_cover(ind_cover)
       nsy=ny_in_cover(ind_cover)
       nsz=nz_in_cover(ind_cover)
       nqx=1+mod(nx_o+nsx+nmodx-1,groups%ngcellx)
       nqy=1+mod(ny_o+nsy+nmody-1,groups%ngcelly)
       nqz=1+mod(nz_o+nsz+nmodz-1,groups%ngcellz)
       xadd=(nx_o+nsx-nqx)*dcellx
       yadd=(ny_o+nsy-nqy)*dcelly
       zadd=(nz_o+nsz-nqz)*dcellz
       ind_qart=(nqx-1)*groups%ngcelly*groups%ngcellz+&
            (nqy-1)*groups%ngcellz+nqz
       set%lab_cell(ind_cover)=ind_qart
       set%lab_cover(ind_cover)=&
            (nsx+set%nspanlx)*set%ncovery*set%ncoverz+&
            (nsy+set%nspanly)*set%ncoverz+(nsz+set%nspanlz)+1
       !!  TM  temporary ????
       set%inv_lab_cover(set%lab_cover(ind_cover))=ind_cover
       !!  TM  temporary ????
       if(members) nm_in_cover = nm_in_cover + groups%nm_group(ind_qart)
       !if(members.AND.groups%nm_group(ind_qart)>0) then
       !   set%n_ing_cover(ind_cover)=groups%nm_group(ind_qart)
       !   if(set%icover_ibeg(ind_cover)+groups%nm_group(ind_qart)-1> &
       !        set%mx_mcover) then
       !      call cq_abort('make_cs: xcover dimension exceeded',& 
       !           set%icover_ibeg(ind_cover)+groups%nm_group(ind_qart)-1, &
       !           set%mx_mcover)
       !   endif
       !   if(groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1> &
       !        mx_mcell) then
       !      call cq_abort('make_cs: x_atom_cover dim. exceeded', &
       !           groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1)
       !   endif
       !   do ni=1,groups%nm_group(ind_qart)
       !      set%xcover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           x_mem_cell(groups%icell_beg(ind_qart)+ni-1)+xadd
       !      set%ycover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           y_mem_cell(groups%icell_beg(ind_qart)+ni-1)+yadd
       !      set%zcover(set%icover_ibeg(ind_cover)+ni-1)= &
       !           z_mem_cell(groups%icell_beg(ind_qart)+ni-1)+zadd
       !   enddo
       !   if(ind_cover.lt.set%ng_cover) then
       !      set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
       !           groups%nm_group(ind_qart)
       !   endif
       !else if(members) then
       !   if(ind_cover < set%ng_cover) then   !TM 26/Jun/2003
       !      set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
       !           groups%nm_group(ind_qart)
       !   endif  ! (ind_cover < set%ng_cover)  !TM 26/Jun/2003
       !   set%n_ing_cover(ind_cover)=0
       !endif ! members.AND.nm_group>0
    enddo
    if(members) then ! Now generate member information
       if(inode==ionode.AND.iprint_index>1) write(io_lun,*) 'Members in covering set: ',nm_in_cover
       if(members) set%icover_ibeg(1)=1
       set%mx_mcover = nm_in_cover
       call start_timer(tmr_std_allocation)
       allocate(set%xcover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members x: ',set%mx_mcover,stat)
       allocate(set%ycover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members y: ',set%mx_mcover,stat)
       allocate(set%zcover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating cover set members z: ',set%mx_mcover,stat)
       allocate(set%ig_cover(set%mx_mcover),STAT=stat)
       if(stat/=0) call cq_abort('Error allocating ig_cover:',set%mx_mcover,stat)
       call reg_alloc_mem(area_index, 3*set%mx_mcover,type_dbl)
       call stop_timer(tmr_std_allocation)
       do ind_cover=1,set%ng_cover
          nsx=nx_in_cover(ind_cover)
          nsy=ny_in_cover(ind_cover)
          nsz=nz_in_cover(ind_cover)
          nqx=1+mod(nx_o+nsx+nmodx-1,groups%ngcellx)
          nqy=1+mod(ny_o+nsy+nmody-1,groups%ngcelly)
          nqz=1+mod(nz_o+nsz+nmodz-1,groups%ngcellz)
          xadd=(nx_o+nsx-nqx)*dcellx
          yadd=(ny_o+nsy-nqy)*dcelly
          zadd=(nz_o+nsz-nqz)*dcellz
          ind_qart= set%lab_cell(ind_cover)
          if(groups%nm_group(ind_qart)>0) then
             set%n_ing_cover(ind_cover)=groups%nm_group(ind_qart)
             if(groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1> mx_mcell) then
                call cq_abort('make_cs: x_atom_cover dim. exceeded', &
                     groups%icell_beg(ind_qart)+groups%nm_group(ind_qart)-1)
             endif
             do ni=1,groups%nm_group(ind_qart)
                set%xcover(set%icover_ibeg(ind_cover)+ni-1)= &
                     x_mem_cell(groups%icell_beg(ind_qart)+ni-1)+xadd
                set%ycover(set%icover_ibeg(ind_cover)+ni-1)= &
                     y_mem_cell(groups%icell_beg(ind_qart)+ni-1)+yadd
                set%zcover(set%icover_ibeg(ind_cover)+ni-1)= &
                     z_mem_cell(groups%icell_beg(ind_qart)+ni-1)+zadd
                set%ig_cover(set%icover_ibeg(ind_cover)+ni-1)= &
                     id_glob(groups%icell_beg(ind_qart)+ni-1)
             enddo
             if(ind_cover.lt.set%ng_cover) then
                set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
                     groups%nm_group(ind_qart)
             endif
          else 
             if(ind_cover < set%ng_cover) then   !TM 26/Jun/2003
                set%icover_ibeg(ind_cover+1)=set%icover_ibeg(ind_cover)+ &
                     groups%nm_group(ind_qart)
             endif  ! (ind_cover < set%ng_cover)  !TM 26/Jun/2003
             set%n_ing_cover(ind_cover)=0
          endif ! nm_group>0
       enddo
    end if
    !if(members) then ! Added by TM 17/07/00
    !   do np=1,set%ng_cover  ! Loop over partitions in GCS
    !      if(set%n_ing_cover(np).gt.0) then  ! Are there atoms ?
    !         if(set%icover_ibeg(np)+set%n_ing_cover(np)-1.gt.&
    !              set%mx_mcover) then
    !            call cq_abort('make_cs: member index error', &
    !                 set%icover_ibeg(np)+set%n_ing_cover(np)-1)
    !         endif
    !      endif
    !   enddo
    !endif ! (members)
    call start_timer(tmr_std_allocation)
    deallocate(nx_in_cover,ny_in_cover,nz_in_cover,STAT=stat)
    if(stat/=0) call cq_abort("Error deallocating nx_in_cover: ",set%ng_cover,stat)
    call reg_dealloc_mem(area_index, 3*set%ng_cover,type_int)
    call stop_timer(tmr_std_allocation)
    call stop_timer(tmr_std_indexing)
    return
  end subroutine make_cs
!!***

!!****f* cover_module/make_iprim *
!!
!!  NAME 
!!   make_iprim
!!  USAGE
!! 
!!  PURPOSE
!!   sbrt make_iprim: creates iprim_group, which gives the covering set CC 
!!   label for each primary set group - I'm pretty sure that this is 
!!   right, but it's a really hard problem
!!
!!   Definitely right - it was broken a while back, but now fixed 
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added timer
!!  SOURCE
!!
  subroutine make_iprim(set,prim)

    use basic_types
    use GenComms, ONLY: cq_abort
    use memory_module, ONLY: reg_alloc_mem, type_int
    use global_module, ONLY: area_index

    implicit none

    ! Passed variables
    type(cover_set) :: set
    type(primary_set) :: prim

    ! Local variables
    integer :: pr,nsx,nsy,nsz,ind_group,i,j,stat

    call start_timer(tmr_std_indexing)
    if(.NOT.ASSOCIATED(set%iprim_group)) then
       call start_timer(tmr_std_allocation)
       allocate(set%iprim_group(prim%mx_iprim),STAT=stat)
       if(stat/=0) call cq_abort('make_iprim: error allocating memory')
       call reg_alloc_mem(area_index,prim%mx_iprim,type_int)
       call stop_timer(tmr_std_allocation)
    endif
    pr = 0 ! Indexes primary members
    do i=1,prim%groups_on_node
       nsx = prim%idisp_primx(i)
       nsy = prim%idisp_primy(i)
       nsz = prim%idisp_primz(i)
       ind_group = (nsx+set%nspanlx)*set%ncovery*set%ncoverz+&
            (nsy+set%nspanly)*set%ncoverz+(nsz+set%nspanlz)+1
       if(prim%nm_nodgroup(i) > 0) then
          do j=1,prim%nm_nodgroup(i)
             pr = pr+1
             if(pr>prim%mx_iprim) call cq_abort('make_iprim: error in pr index',pr)
             set%iprim_group(pr) = ind_group
          enddo
       endif ! (prim%nm_nodgroup(i) > 0) then
    enddo
    call stop_timer(tmr_std_indexing)
  end subroutine make_iprim
!!***

!!****f* cover_module/send_ncover *
!!
!!  NAME 
!!   send_ncover
!!  USAGE
!! 
!!  PURPOSE
!!   Once make_cs has been run on every processor, this routine
!!   sends sizes of the CS to all other processors (and gathers theirs)
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    ROBODoc header and changed to gsum
!!   2008/08/12 10:49 dave
!!    Tidying: removed unnecessary variable ncovers and allocate/deallocate calls
!!  SOURCE
!!
! =====================================================================
!   sbrt send_ncover: Called once at the start of a run (after make_gcs), 
!   to ensure that all nodes know all other nodes GCS sizes (ncoverx)
!   for transpose work
! ---------------------------------------------------------------------
  subroutine send_ncover(set,nnd)

    ! Module declaration
    use global_module, ONLY: numprocs
    use GenComms, ONLY: gsum, cq_abort

    implicit none

    ! Passed variables
    type(cover_set) :: set
    integer :: nnd ! Node no. starting from one

    set%ncover_rem = 0
    set%ncover_rem(3*(nnd-1)+1) = set%ncoverx
    set%ncover_rem(3*(nnd-1)+2) = set%ncovery
    set%ncover_rem(3*(nnd-1)+3) = set%ncoverz
    call gsum(set%ncover_rem, 3*numprocs)
    return
  end subroutine send_ncover
!!***

!!****f* cover_module/convert_primary *
!!
!!  NAME 
!!   convert_primary
!!  USAGE
!! 
!!  PURPOSE
!!   Given a primary set and a cutoff, creates the widths and left-spans
!!   of the CS in the specified groups and returns the origin of the CS
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99 ?
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and used cq_abort
!!  SOURCE
!!
  subroutine convert_primary(groups,prim,set,ncx_o,ncy_o,ncz_o)

    use datatypes
    use global_module
    use basic_types
    use numbers, ONLY: very_small
    use GenComms, ONLY: cq_abort

    implicit none

    ! Passed variables
    type(group_set)   :: groups
    type(primary_set) :: prim
    type(cover_set)   :: set
    integer :: ncx_o,ncy_o,ncz_o                ! Origin for cover set

    ! Local variables
    real(double) :: ro_x,ro_y,ro_z              ! Origin in reals
    real(double) :: ro_cx,ro_cy,ro_cz           ! CS origin in reals
    real(double) :: dx,dy,dz,dcx,dcy,dcz        ! Sizes of groups in reals
    real(double) :: lhx,lhy,lhz,lhcx,lhcy,lhcz  ! Left hand corners
    real(double) :: rhx,rhy,rhz,rhcx,rhcy,rhcz  ! Right hand corners
    real(double) :: cslx,csly,cslz              ! Covering set left

    if(.NOT.associated(prim%groups)) then
       call cq_abort('convert_primary: primary set groups pointer not assoc')
    endif
    ! First, find sizes of groups in terms of unit cell
    dx =  rcellx/real(prim%groups%ngcellx,double)   ! Check for a more portable way !
    dy =  rcelly/real(prim%groups%ngcelly,double)
    dz =  rcellz/real(prim%groups%ngcellz,double)
    dcx = rcellx/real(groups%ngcellx,double)
    dcy = rcelly/real(groups%ngcelly,double)
    dcz = rcellz/real(groups%ngcellz,double)
    ! Convert origin of prim to reals (add very_small to prevent ambiguity)
    ro_x = (real(prim%nx_origin-1,double))*dx + very_small 
    ro_y = (real(prim%ny_origin-1,double))*dy + very_small
    ro_z = (real(prim%nz_origin-1,double))*dz + very_small
    ! Convert this to number of CS groups
    ro_cx = ro_x/dcx
    ro_cy = ro_y/dcy
    ro_cz = ro_z/dcz
    ! Round and offset (partition 1,1,1 has LH corner at 0.0,0.0,0.0)
    ! Accept implicit type conversion
    ncx_o = 1+anint(ro_cx)  
    ncy_o = 1+anint(ro_cy)
    ncz_o = 1+anint(ro_cz)
    set%nx_origin = ncx_o
    set%ny_origin = ncy_o
    set%nz_origin = ncz_o
    ! Find CS origin in reals
    ro_cx = (real(ncx_o-1,double))*dcx
    ro_cy = (real(ncy_o-1,double))*dcy
    ro_cz = (real(ncz_o-1,double))*dcz
    ! Start by finding the left and right hand corners of the CS
    lhx = dx*(prim%nx_origin-1-prim%nleftx)-set%rcut
    lhy = dy*(prim%ny_origin-1-prim%nlefty)-set%rcut
    lhz = dz*(prim%nz_origin-1-prim%nleftz)-set%rcut
    rhx = dx*(prim%nx_origin-1-prim%nleftx+prim%nw_primx)+set%rcut  
    rhy = dy*(prim%ny_origin-1-prim%nlefty+prim%nw_primy)+set%rcut
    rhz = dz*(prim%nz_origin-1-prim%nleftz+prim%nw_primz)+set%rcut
    ! The cover set left span must be offset from the CS origin
    cslx = (ro_cx - lhx)/dcx
    csly = (ro_cy - lhy)/dcy
    cslz = (ro_cz - lhz)/dcz
    ! Convert to integers
    set%nspanlx = ceiling(cslx)
    set%nspanly = ceiling(csly)
    set%nspanlz = ceiling(cslz)
    ! The width of the CS must be the distance from the LH corner to the 
    ! RH corner - this way we pick the SMALLEST CS needed.
    set%ncoverx = ceiling(rhx/dcx)-(ncx_o-1-set%nspanlx)
    set%ncovery = ceiling(rhy/dcy)-(ncy_o-1-set%nspanly)
    set%ncoverz = ceiling(rhz/dcz)-(ncz_o-1-set%nspanlz)
    ! Calculate the corners of the CS for checks
    lhcx = dcx*(ncx_o-1-set%nspanlx)
    lhcy = dcy*(ncy_o-1-set%nspanly)
    lhcz = dcz*(ncz_o-1-set%nspanlz)
    rhcx = dcx*(ncx_o-1-set%nspanlx+set%ncoverx)
    rhcy = dcy*(ncy_o-1-set%nspanly+set%ncovery)
    rhcz = dcz*(ncz_o-1-set%nspanlz+set%ncoverz)
    ! First check that the LH and RH corners are large enough
    ! LH
    if(lhx+very_small<lhcx) then
       write(io_lun,*) 'CS too small; adjusting nspanlx',lhx,lhcx
       set%nspanlx = set%nspanlx+ceiling((lhcx-lhx)/dcx)
    endif
    if(lhy+very_small<lhcy) then
       write(io_lun,*) 'CS too small; adjusting nspanly',lhy,lhcy
       set%nspanly = set%nspanly+ceiling((lhcy-lhy)/dcy)
    endif
    if(lhz+very_small<lhcz) then
       write(io_lun,*) 'CS too small; adjusting nspanlz',lhz,lhcz
       set%nspanlz = set%nspanlz+ceiling((lhcz-lhz)/dcz)
    endif
    ! RH
    if(rhx>rhcx+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncoverx',rhx,rhcx
       set%ncoverx = set%ncoverx+ceiling((rhx-rhcx)/dcx)
    endif
    if(rhy>rhcy+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncovery',rhy,rhcy
       set%ncovery = set%ncovery+ceiling((rhy-rhcy)/dcy)
    endif
    if(rhz>rhcz+very_small) then
       write(io_lun,*) 'CS too small; adjusting ncoverz',rhz,rhcz
       set%ncoverz = set%ncoverz+ceiling((rhz-rhcz)/dcz)
    endif
    ! Now check that the CS isn't too large
    ! LH
    if(lhx-lhcx>dcx) then
       write(io_lun,*) 'CS too big; adjusting nspanlx',lhx,lhcx
       set%nspanlx = set%nspanlx-floor((lhx-lhcx)/dcx)
    endif
    if(lhy-lhcy>dcy) then
       write(io_lun,*) 'CS too big; adjusting nspanly',lhy,lhcy
       set%nspanly = set%nspanly-floor((lhy-lhcy)/dcy)
    endif
    if(lhz-lhcz>dcz) then
       write(io_lun,*) 'CS too big; adjusting nspanlz',lhz,lhcz
       set%nspanlz = set%nspanlz-floor((lhz-lhcz)/dcz)
    endif
    ! RH
    if(rhcx-rhx>dcx) then
       write(io_lun,*) 'CS too big; adjusting ncoverx',rhx,rhcx
       set%ncoverx = set%ncoverx-floor((rhx-rhcx)/dcx)
    endif
    if(rhcy-rhy>dcy) then
       write(io_lun,*) 'CS too big; adjusting ncovery',rhy,rhcy
       set%ncovery = set%ncovery-floor((rhy-rhcy)/dcy)
    endif
    if(rhcz-rhz>dcz) then
       write(io_lun,*) 'CS too big; adjusting ncoverz',rhz,rhcz
       set%ncoverz = set%ncoverz-floor((rhz-rhcz)/dcz)
    endif
    return
  end subroutine convert_primary
!!***

!!****f* cover_module/allocate_cs *
!!
!!  NAME 
!!   allocate_cs
!!  USAGE
!! 
!!  PURPOSE
!!   allocates memory to the CS derived type
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave 
!!    Added ROBODoc header and used cq_abort
!!   2008/05/16 ast
!!    Added timer
!!  SOURCE
!!
  subroutine allocate_cs(set,mx_node,members)

    ! Module usage
    use basic_types
    use GenComms, ONLY: cq_abort
    use global_module, ONLY: area_index
    use memory_module, ONLY: reg_alloc_mem, type_int

    implicit none

    ! Passed variables
    integer :: mx_node
    type(cover_set) :: set
    logical :: members

    ! Local variables
    integer :: stat

    call start_timer(tmr_std_allocation)
    allocate(set%lab_cell(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(1)')
    endif
    allocate(set%lab_cover(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(2)')
    endif
    !TM temporary ???
    allocate(set%inv_lab_cover(set%mx_gcover),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(3)')
    endif
    !TM temporary ???
    allocate(set%ncover_rem(3*mx_node),STAT=stat)
    if(stat/=0) then
       call cq_abort('allocate_cs: error(4)')
    endif
    call reg_alloc_mem(area_index,3*set%mx_gcover+3*mx_node,type_int)
    ! Be safe with pointer
    nullify(set%iprim_group)
    if(members) then
       allocate(set%n_ing_cover(set%mx_gcover),STAT=stat)
       if(stat/=0) then
          call cq_abort('allocate_cs: error(5)')
       endif
       allocate(set%icover_ibeg(set%mx_gcover),STAT=stat)
       if(stat/=0) then
          call cq_abort('allocate_cs: error(6)')
       endif
       call reg_alloc_mem(area_index,2*set%mx_gcover,type_int)
       ! These will be allocated later !
       !allocate(set%xcover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(7)')
       !endif
       !allocate(set%ycover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(8)')
       !endif
       !allocate(set%zcover(set%mx_mcover),STAT=stat)
       !if(stat/=0) then
       !   call cq_abort('allocate_cs: error(9)')
       !endif
    else
       nullify(set%n_ing_cover,set%icover_ibeg, &
            set%xcover,set%ycover,set%zcover)
    endif
    call stop_timer(tmr_std_allocation)
    return
  end subroutine allocate_cs
!!***

!!****f* cover_module/deallocate_cs *
!!
!!  NAME 
!!   deallocate_cs
!!  USAGE
!! 
!!  PURPOSE
!!   Deallocates memory from CS derived type
!!  INPUTS
!! 
!! 
!!  USES
!! 
!!  AUTHOR
!!   D.R.Bowler
!!  CREATION DATE
!!   24/11/99
!!  MODIFICATION HISTORY
!!   21/06/2001 dave
!!    Added ROBODoc header and cq_abort
!!   2008/05/16 ast
!!    Added timer
!!   2018/07/04 08:47 dave
!!    Added members logical to be consistent with allocate
!!   2018/09/06 10:35 dave
!!    Tidying and bug fix: members was bracketing ALL variables and
!!    stat wasn't initialised  
!!  SOURCE
!!
  subroutine deallocate_cs(set,members)

    ! Module usage
    use basic_types
    use GenComms, ONLY: cq_abort
    use global_module, ONLY: area_index
    use memory_module, ONLY: reg_dealloc_mem, type_int

    implicit none

    ! Passed variables
    type(cover_set) :: set
    logical :: members

    ! Local variables
    integer :: stat

    stat = 0
    call start_timer(tmr_std_allocation)
    call reg_dealloc_mem(area_index,2*set%mx_gcover,type_int)
    call reg_dealloc_mem(area_index,3*set%mx_gcover,type_int)
    call reg_dealloc_mem(area_index,size(set%ncover_rem),type_int)
    if(ASSOCIATED(set%iprim_group)) then
       call reg_dealloc_mem(area_index,size(set%iprim_group),type_int)
       deallocate(set%iprim_group,STAT=stat)
       if(stat/=0) call cq_abort('deallocate_cs: error(1)')
    endif
    deallocate(set%ncover_rem, set%inv_lab_cover, set%lab_cover,set%lab_cell, STAT=stat)
    if(stat/=0) call cq_abort('deallocate_cs: error(2)')
    if(members) then
       deallocate(set%zcover,set%ycover,set%xcover, set%icover_ibeg,set%n_ing_cover, STAT=stat)
       if(stat/=0) call cq_abort('deallocate_cs: error(3)')
    endif
    call stop_timer(tmr_std_allocation)
    return
  end subroutine deallocate_cs
!!***
  
end module cover_module
