#include "macrodef.fh"

#ifdef Vertex
#ifdef Cell
#error Both Cell and Vertex are defined
#endif    
!---------------------------------------------------------------------------------------------------
! copy a point of data into data target for vertext center code
!---------------------------------------------------------------------------------------------------
  subroutine pointcopy(wei,llbout,uubout,ext_out,data_out,xx,yy,zz,dv)
  implicit none
  integer,intent(in) :: wei
  integer,dimension(3),intent(in) ::ext_out
  real*8,dimension(3) :: llbout,uubout
  real*8,dimension(ext_out(1),ext_out(2),ext_out(3)),intent(inout)::data_out
  real*8,intent(in) :: xx,yy,zz,dv

  real*8,dimension(3) :: ho
  integer :: i,j,k

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::pointcopy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

!!!   
  if(any(ext_out == 1))then
       write(*,*)"fmisc.f90::pointcopy: meets iolated points for out data"
       write(*,*) llbout,uubout
       stop
  else
    ho = (uubout-llbout)/(ext_out-1)
  endif    
  i = idint((xx-llbout(1))/ho(1)+0.4)+1
  j = idint((yy-llbout(2))/ho(2)+0.4)+1
  k = idint((zz-llbout(3))/ho(3)+0.4)+1

  if(i<1 .or. i>ext_out(1) .or. &
     j<1 .or. j>ext_out(2) .or. &
     k<1 .or. k>ext_out(3) )then
     write(*,*)"i,j,k = ",i,j,k
     write(*,*)"ext = ",ext_out
     stop
  endif
  if(dabs(llbout(1)+(i-1)*ho(1)-xx)>ho(1)/2 .or. &
     dabs(llbout(2)+(j-1)*ho(2)-yy)>ho(2)/2 .or. &
     dabs(llbout(3)+(k-1)*ho(3)-zz)>ho(3)/2 )then    
     write(*,*)"fmisc.f90::pointcopy: llbout = ",llbout
     write(*,*)"fmisc.f90::pointcopy: ho = ",ho
     write(*,*)"fmisc.f90::pointcopy: x,y,z = ",llbout(1)+(i-1)*ho(1),llbout(2)+(j-1)*ho(2),llbout(3)+(k-1)*ho(3)
     write(*,*)"fmisc.f90::pointcopy: point = ",xx,yy,zz
     stop
   endif

  data_out(i,j,k)=dv

  return

  end subroutine pointcopy
!---------------------------------------------------------------------------------------------------
! copy a part of data from data source, for vertex center code
!---------------------------------------------------------------------------------------------------
  subroutine copy(wei,llbout,uubout,ext_out,data_out,llbin,uubin,ext_in,data_in,lcopy,ucopy)
  implicit none
  integer,intent(in) :: wei
  integer,dimension(3),intent(in) ::ext_out,ext_in
  real*8,dimension(3),intent(in) :: lcopy,ucopy
  real*8,dimension(3) :: llbout,uubout,llbin,uubin
  real*8,dimension(ext_out(1),ext_out(2),ext_out(3)),intent(inout)::data_out
  real*8,dimension(ext_in(1),ext_in(2),ext_in(3)),intent(in)::data_in

  real*8,dimension(3) :: ho,hi
  integer,dimension(3) :: illo,iuuo,illi,iuui

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::copy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

!!!   
  if(any(ext_out == 1))then
    if(any(ext_in == 1))then
       write(*,*)"fmisc.f90::copy: meets iolated points for both in and out data"
       write(*,*) llbin,uubin
       write(*,*) llbout,uubout
       stop
     else
       hi = (uubin-llbin)/(ext_in-1)
       ho = hi
    endif
  else
    ho = (uubout-llbout)/(ext_out-1)
    if(any(ext_in == 1))then
      hi = ho
    else
      hi = (uubin-llbin)/(ext_in-1)
      if(any(abs(hi-ho) > min(hi,ho)/2))then
         write(*,*)"fmisc.f90::copy: meets copy reqest for different numerical grid"
         write(*,*)hi,ho
         stop
      endif
    endif
  endif    
  illo = idint((lcopy-llbout)/ho+0.4)+1
  iuuo = ext_out - idint((uubout-ucopy)/ho+0.4)
  illi = idint((lcopy-llbin)/hi+0.4)+1
  iuui = ext_in - idint((uubin-ucopy)/hi+0.4)

  if(any(llbout-lcopy>ho/2) .or. any(ucopy-uubout>ho/2))then
     write(*,*)"fmisc.f90::copy: llbout = ",llbout
     write(*,*)"fmisc.f90::copy: uubout = ",uubout
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     write(*,*)"fmisc.f90::copy: ho = ",ho
     write(*,*)llbout-lcopy,ucopy-uubout
     stop
  elseif(any(llbin -lcopy>hi/2) .or. any(ucopy-uubin >hi/2))then
     write(*,*)"fmisc.f90::copy: llbin = ",llbin
     write(*,*)"fmisc.f90::copy: uubin = ",uubin
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     stop
  elseif(any(illo<1) .or. any(illi<1) .or. any(illo-iuuo>0) .or. any(illi-iuui>0) .or. &
         any(iuui-ext_in>0) .or. any(iuuo-ext_out>0))then
     write(*,*)"fmisc.f90::copy: illi = ",illi
     write(*,*)"fmisc.f90::copy: iuui = ",iuui
     write(*,*)"fmisc.f90::copy: illo = ",illo
     write(*,*)"fmisc.f90::copy: iuuo = ",iuuo     
     write(*,*)"fmisc.f90::copy: llbout = ",llbout
     write(*,*)"fmisc.f90::copy: uubout = ",uubout
     write(*,*)"fmisc.f90::copy: llbin = ",llbin
     write(*,*)"fmisc.f90::copy: uubin = ",uubin
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     stop
   endif

  data_out(illo(1):iuuo(1),illo(2):iuuo(2),illo(3):iuuo(3))=data_in(illi(1):iuui(1),illi(2):iuui(2),illi(3):iuui(3))

  return

  end subroutine copy
!-----------------------------------------------------------------------------------------------------------------
! three dimensional interpolation for vertex center grid structure  
  subroutine global_interp(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry,ORDN
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA

!~~~~~~> Other parameters:

  integer :: j,m,imin,jmin,kmin
  integer,dimension(3) :: cxB,cxT,cxI,cmin,cmax
  real*8,dimension(3) :: cx
  real*8, dimension(1:ORDN) :: x1a
  real*8, dimension(1:ORDN,1:ORDN,1:ORDN) :: ya
  integer, parameter :: NO_SYMM = 0, EQUATORIAL = 1, OCTANT = 2
  real*8 :: dX,dY,dZ,ddy
  real*8, parameter :: ONE=1.d0
  logical::decide3d

  imin = lbound(f,1)
  jmin = lbound(f,2)
  kmin = lbound(f,3)

  dX = X(imin+1)-X(imin)
  dY = Y(jmin+1)-Y(jmin)
  dZ = Z(kmin+1)-Z(kmin)

  forall( j = 1:ordn ) x1a(j) = ( j - 1 )* ONE

  cxI(1) = idint((x1-X(1))/dX+0.4)+1
  cxI(2) = idint((y1-Y(1))/dY+0.4)+1
  cxI(3) = idint((z1-Z(1))/dZ+0.4)+1

  cxB = cxI - ORDN/2+1
  cxT = cxB + ORDN - 1
       
  cmin = 1
  cmax = ex
  if(Symmetry == OCTANT  .and.dabs(X(1))<dX) cmin(1) = -ORDN/2+2
  if(Symmetry == OCTANT  .and.dabs(Y(1))<dY) cmin(2) = -ORDN/2+2
  if(Symmetry /= NO_SYMM .and.dabs(Z(1))<dZ) cmin(3) = -ORDN/2+2
  do m =1,3
   if(cxB(m) < cmin(m))then
      cxB(m) = cmin(m)
      cxT(m) = cxB(m) + ORDN - 1
   endif
   if(cxT(m) > cmax(m))then
      cxT(m) = cmax(m)
      cxB(m) = cxT(m) + 1 - ORDN
   endif
 enddo
 if(cxB(1)>0)then
  cx(1) = (x1 - X(cxB(1)))/dX
 else
  cx(1) = (x1 + X(2-cxB(1)))/dX
 endif
 if(cxB(2)>0)then
  cx(2) = (y1 - Y(cxB(2)))/dY
 else
  cx(2) = (y1 + Y(2-cxB(2)))/dY
 endif
 if(cxB(3)>0)then
  cx(3) = (z1 - Z(cxB(3)))/dZ
 else
  cx(3) = (z1 + Z(2-cxB(3)))/dZ
 endif

  if(decide3d(ex,f,f,cxB,cxT,SoA,ya,ORDN,Symmetry))then
     write(*,*)"global_interp position: ",x1,y1,z1
     write(*,*)"data range: ",X(1),X(ex(1)),Y(1),Y(ex(2)),Z(1),Z(ex(3))
     stop
  endif
  call polin3(x1a,x1a,x1a,ya,cx(1),cx(2),cx(3),f_int,ddy,ORDN)

  return

  end subroutine global_interp
!----------------------------------------------------------------
! decide which 3d data to be used does not surport PI-Symmetry yet 
!----------------------------------------------------------------
  function decide3d(ex,f,fpi,cxB,cxT,SoA,ya,ORDN,Symmetry)  result(gont)
  implicit none

  integer,                                 intent(in) :: ORDN,Symmetry
  integer,dimension(1:3)                 , intent(in) :: ex,cxB,cxT
  real*8, dimension(1:3)                 , intent(in) :: SoA
  real*8, dimension(ex(1),ex(2),ex(3))   , intent(in) :: f,fpi
  real*8, dimension(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3)), intent(out):: ya
  logical::gont

  integer,dimension(1:3) :: fmin1,fmin2,fmax1,fmax2
  integer::i,j,k,m

  gont=.false.
  do m=1,3
! check cxB and cxT are NaN or not  
    if(.not.(iabs(cxB(m)).ge.0)) gont=.true.
    if(.not.(iabs(cxT(m)).ge.0)) gont=.true.
    fmin1(m) = max(1,cxB(m))
    fmax1(m) = cxT(m)
    fmin2(m) = cxB(m)
    fmax2(m) = min(0,cxT(m))
    if((fmin1(m).le.fmax1(m)).and.(  fmin1(m)<1.or.  fmax1(m)>ex(m)))gont=.true.
    if((fmin2(m).le.fmax2(m)).and.(2-fmax2(m)<1.or.2-fmin2(m)>ex(m)))gont=.true.
  enddo
!sanity check
  if(gont)then
          write(*,*)"error in decide3d"
          write(*,*)((fmin1.le.fmax1).and.(  fmin1<1.or.  fmax1>ex))
          write(*,*)((fmin2.le.fmax2).and.(2-fmax2<1.or.2-fmin2>ex))
          write(*,*)"cxB, cxT and data shape:"
          write(*,*)cxB,cxT,ex
          write(*,*)"resulted fmin1, fmax1 and fmin2, fmax2:"
          write(*,*)fmin1,fmax1,fmin2,fmax2
  else

  do k=fmin1(3),fmax1(3)
     do j=fmin1(2),fmax1(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,j,k)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(2-i,j,k)*SoA(1)
        enddo
     enddo
     do j=fmin2(2),fmax2(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,2-j,k)*SoA(2)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(2-i,2-j,k)*SoA(1)*SoA(2)
        enddo
     enddo
  enddo
           
  do k=fmin2(3),fmax2(3)
     do j=fmin1(2),fmax1(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,j,2-k)*SoA(3)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(2-i,j,2-k)*SoA(1)*SoA(3)
        enddo
     enddo
     do j=fmin2(2),fmax2(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,2-j,2-k)*SoA(2)*SoA(3)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(2-i,2-j,2-k)*SoA(1)*SoA(2)*SoA(3)
        enddo
     enddo
  enddo

  endif

  end function decide3d

!---------------------------------------------------------------------------------------
subroutine symmetry_bd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1),-ord+1:extc(2),-ord+1:extc(3)),intent(out):: funcc
  real*8, dimension(1:3), intent(in) :: SoA

  integer::i

  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
   enddo
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+2,1:extc(3))*SoA(2)
   enddo
   do i=0,ord-1
      funcc(:,:,-i) = funcc(:,:,i+2)*SoA(3)
   enddo

end subroutine symmetry_bd

subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,-ord+1:extc(3)+ord),intent(out):: funcc
  real*8, dimension(1:3), intent(in) :: SoA

  integer::i

  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
   enddo
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+2,1:extc(3))*SoA(2)
      funcc(:,extc(2)+1+i,1:extc(3)) = funcc(:,extc(2)-1-i,1:extc(3))*SoA(2)
   enddo
   do i=0,ord-1
      funcc(:,:,-i) = funcc(:,:,i+2)*SoA(3)
      funcc(:,:,extc(3)+1+i) = funcc(:,:,extc(3)-1-i)*SoA(3)
   enddo

end subroutine symmetry_tbd

subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,extc(3)),intent(out):: funcc
  real*8, dimension(2), intent(in) :: SoA

  integer::i

  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1)
   enddo
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+2,1:extc(3))*SoA(2)
      funcc(:,extc(2)+1+i,1:extc(3)) = funcc(:,extc(2)-1-i,1:extc(3))*SoA(2)
   enddo

end subroutine symmetry_stbd

subroutine symmetry_sntbd(ord,extc,func,funcc,SoA,actd)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord,actd
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,extc(3)),intent(out):: funcc
  real*8, intent(in) :: SoA

  integer::i

  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
  if(actd==0)then
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA
   enddo
  elseif(actd==1)then
   do i=0,ord-1
      funcc(1:extc(1),-i,1:extc(3)) = funcc(1:extc(1),i+2,1:extc(3))*SoA
      funcc(1:extc(1),extc(2)+1+i,1:extc(3)) = funcc(1:extc(1),extc(2)-1-i,1:extc(3))*SoA
   enddo
  else
    write(*,*)"symmetry_sntbd: not recognized actd = ",actd
  endif

end subroutine symmetry_sntbd


subroutine d2dump(wei,llb,uub,ext,data_in,data_out,gord,SoA)
  implicit none
  integer,             intent(in) :: wei,gord
  integer,dimension(3),intent(in) :: ext
  real*8, dimension(3),intent(in) :: SoA
  real*8, dimension(3) :: llb,uub
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in)   ::data_in
  real*8, dimension(ext(1),ext(2)),       intent(inout)::data_out

  real*8  :: dZ
  integer :: i,j,k

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::copy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

  dZ = (uub(3)-llb(3))/(ext(3)-1)
  k = idint((0-llb(3))/dZ+0.4)+1
     
  if(k < 1)then
      write(*,*) "d2dump: something must be wrong"
      return
  endif

  data_out(i,j) = data_in(i,j,k)

end subroutine d2dump

#else
#ifdef Cell
!subroutine interp_2 support cell center only
!-----------------------------------------------------------------------------
!
! Interpolate function f using weights Delx, Dely and Delz
!
!-----------------------------------------------------------------------------

  subroutine interp_2(ex,f,f_int,il,iu,jl,ju,kl,ku,Dx,Dy,Dz,&
                                     ordn,SoA,symmetry)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  integer,                             intent(in) :: il,iu,jl,ju,kl,ku,ordn
  real*8,                              intent(in) :: Dx,Dy,Dz,SoA(3)

!~~~~~~> Other parameters:

  integer :: j,imin,jmin,kmin
  real*8, dimension(1:ordn) :: x1a
  real*8, dimension(1:ordn,1:ordn,1:ordn) :: ya
  real*8, parameter :: ONE=1.d0
  integer, parameter :: NO_SYMM = 0, EQUATORIAL = 1, OCTANT = 2
  real*8 :: ddy,symX,symY,symZ

  symX = SoA(1)
  symY = SoA(2)
  symZ = SoA(3)
 
  imin = lbound(f,1)
  jmin = lbound(f,2)
  kmin = lbound(f,3)

  forall( j = 1:ordn ) x1a(j) = ( j - 1 )* ONE

  ya(2:ordn,2:ordn,2:ordn) = f(il+1:iu,jl+1:ju,kl+1:ku)
   
  if( il < imin .and. symmetry < OCTANT ) then
      write(*,*) 'Error in interp_2!!!'
      stop
  endif
  if( il < imin ) then
    ya(1,2:ordn,2:ordn) = f(imin,jl+1:ju,kl+1:ku)* symX
  else
    ya(1,2:ordn,2:ordn) = f(il  ,jl+1:ju,kl+1:ku)
  endif
  
  if( jl < jmin .and. symmetry < OCTANT ) then
     write(*,*) 'Error in interp_2!!!'
     stop
  endif
   
  if( jl < jmin ) then
    ya(2:ordn,1,2:ordn) = f(il+1:iu,jmin,kl+1:ku)* symY
  else
    ya(2:ordn,1,2:ordn) = f(il+1:iu,jl,kl+1:ku)
  endif

  if( kl < kmin .and. symmetry < EQUATORIAL ) then
    write(*,*)  'Error in interp_2!!!'
    stop
  endif

  if( kl < kmin ) then
   ya(2:ordn,2:ordn,1) = f(il+1:iu,jl+1:ju,kmin)* symZ
  else
   ya(2:ordn,2:ordn,1) = f(il+1:iu,jl+1:ju,kl  )
  endif

  if( il < imin .and. jl < jmin ) then
   ya(1,1,2:ordn) = f(imin,jmin,kl+1:ku)* symX * symY
  else if( il >= imin .and. jl < jmin ) then
   ya(1,1,2:ordn) = f(il,jmin,kl+1:ku)* symY
  else if( il < imin .and. jl >= jmin ) then
   ya(1,1,2:ordn) = f(imin,jl,kl+1:ku)* symX
  else 
   ya(1,1,2:ordn) = f(il,jl,kl+1:ku)
  endif
   
  if( il < imin .and. kl < kmin ) then
   ya(1,2:ordn,1) = f(imin,jl+1:ju,kmin)* symX * symZ
  else if( il >= imin .and. kl < kmin ) then
   ya(1,2:ordn,1) = f(il,jl+1:ju,kmin)* symZ
  else if( il < imin .and. kl >= kmin ) then
   ya(1,2:ordn,1) = f(imin,jl+1:ju,kl)* symX
  else
   ya(1,2:ordn,1) = f(il,jl+1:ju,kl)
  endif

  if( jl < jmin .and. kl < kmin ) then
    ya(2:ordn,1,1) = f(il+1:iu,jmin,kmin)* symY * symZ
  else if( jl >= jmin .and. kl < kmin ) then
    ya(2:ordn,1,1) =   f(il+1:iu,jl,kmin)* symZ
  else if( jl < jmin .and. kl >= kmin ) then
    ya(2:ordn,1,1) =   f(il+1:iu,jmin,kl)* symY
  else
   ya(2:ordn,1,1) = f(il+1:iu,jl,kl)
  endif

  if( il < imin ) then
   if( jl < jmin .and. kl < kmin) then
    ya(1,1,1) = f(imin,jmin,kmin)* symX * symY * symZ
   else if( jl >= jmin .and. kl < kmin ) then
    ya(1,1,1) = f(imin,jl,kmin)* symX * symZ
   else if( jl < jmin .and. kl >= kmin ) then
    ya(1,1,1) = f(imin,jmin,kl)* symX * symY
   else
    ya(1,1,1) = f(imin,jl,kl)* symX
   endif
  else
   if( jl < jmin .and. kl < kmin) then
    ya(1,1,1) = f(il,jmin,kmin)* symY * symZ
   else if( jl >= jmin .and. kl < kmin ) then
    ya(1,1,1) = f(il,jl,kmin)* symZ
   else if( jl < jmin .and. kl >= kmin ) then
    ya(1,1,1) = f(il,jmin,kl)* symY
   else
    ya(1,1,1) = f(il,jl,kl)
   endif
  endif
  
  call polin3(x1a,x1a,x1a,ya,Dx,Dy,Dz,f_int,ddy,ordn)

  if(.not.(dabs(f_int).ge.0))then
    write(*,*)"find nan in interp_2:",f_int,"inputs are:"
!    write(*,*)ya
!    write(*,*)"-----------------------------------------"
!    write(*,*)f(il:iu,jl:ju,kl:ku)
    write(*,*)Dx,Dy,Dz,symx,symy,symz,ordn
    write(*,*)il,iu,jl,ju,kl,ku,ex,symmetry
  endif

  return

  end subroutine interp_2
!---------------------------------------------------------------------------------------------------
! copy a point of data into data target for vertext center code
!---------------------------------------------------------------------------------------------------
  subroutine pointcopy(wei,llbout,uubout,ext_out,data_out,xx,yy,zz,dv)
  implicit none
  integer,intent(in) :: wei
  integer,dimension(3),intent(in) ::ext_out
  real*8,dimension(3) :: llbout,uubout
  real*8,dimension(ext_out(1),ext_out(2),ext_out(3)),intent(inout)::data_out
  real*8,intent(in) :: xx,yy,zz,dv

  real*8,dimension(3) :: ho
  integer :: i,j,k

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::pointcopy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

!!!   
  ho = (uubout-llbout)/ext_out
  i = idint((xx-llbout(1))/ho(1)+0.4)+1
  j = idint((yy-llbout(2))/ho(2)+0.4)+1
  k = idint((zz-llbout(3))/ho(3)+0.4)+1

  if(i<1 .or. i>ext_out(1) .or. &
     j<1 .or. j>ext_out(2) .or. &
     k<1 .or. k>ext_out(3) )then
     write(*,*)"i,j,k = ",i,j,k
     write(*,*)"ext = ",ext_out
     stop
  endif
  if(dabs(llbout(1)+(i-0.5)*ho(1)-xx)>ho(1)/2 .or. &
     dabs(llbout(2)+(j-0.5)*ho(2)-yy)>ho(2)/2 .or. &
     dabs(llbout(3)+(k-0.5)*ho(3)-zz)>ho(3)/2 )then    
     write(*,*)"fmisc.f90::pointcopy: llbout = ",llbout
     write(*,*)"fmisc.f90::pointcopy: ho = ",ho
     write(*,*)"fmisc.f90::pointcopy: x,y,z = ",llbout(1)+(i-0.5)*ho(1),llbout(2)+(j-0.5)*ho(2),llbout(3)+(k-0.5)*ho(3)
     write(*,*)"fmisc.f90::pointcopy: point = ",xx,yy,zz
     stop
   endif

  data_out(i,j,k)=dv

  return

  end subroutine pointcopy
!---------------------------------------------------------------------------------------------------
! copy a part of data from data source, for cell center code
!---------------------------------------------------------------------------------------------------
  subroutine copy(wei,llbout,uubout,ext_out,data_out,llbin,uubin,ext_in,data_in,lcopy,ucopy)
  implicit none
  integer,intent(in) :: wei
  integer,dimension(3),intent(in) ::ext_out,ext_in
  real*8,dimension(3),intent(in) :: lcopy,ucopy
  real*8,dimension(3) :: llbout,uubout,llbin,uubin
  real*8,dimension(ext_out(1),ext_out(2),ext_out(3)),intent(inout)::data_out
  real*8,dimension(ext_in(1),ext_in(2),ext_in(3)),intent(in)::data_in

  real*8,dimension(3) :: ho,hi
  integer,dimension(3) :: illo,iuuo,illi,iuui

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::copy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

!!!   
  ho = (uubout-llbout)/ext_out
  hi = (uubin-llbin)/ext_in
  illo = idint((lcopy-llbout)/ho+0.4)+1
  iuuo = ext_out - idint((uubout-ucopy)/ho+0.4)
  illi = idint((lcopy-llbin)/hi+0.4)+1
  iuui = ext_in - idint((uubin-ucopy)/hi+0.4)

  if(any(llbout-lcopy>ho/2) .or. any(ucopy-uubout>ho/2))then
     write(*,*)"fmisc.f90::copy: llbout = ",llbout
     write(*,*)"fmisc.f90::copy: uubout = ",uubout
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     write(*,*)"fmisc.f90::copy: ho = ",ho
     write(*,*)llbout-lcopy,ucopy-uubout
     stop
  elseif(any(llbin -lcopy>hi/2) .or. any(ucopy-uubin >hi/2))then
     write(*,*)"fmisc.f90::copy: llbin = ",llbin
     write(*,*)"fmisc.f90::copy: uubin = ",uubin
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     stop
  elseif(any(illo<1) .or. any(illi<1) .or. any(illo-iuuo>0) .or. any(illi-iuui>0) .or. &
         any(iuui-ext_in>0) .or. any(iuuo-ext_out>0))then
     write(*,*)"fmisc.f90::copy: illi = ",illi
     write(*,*)"fmisc.f90::copy: iuui = ",iuui
     write(*,*)"fmisc.f90::copy: illo = ",illo
     write(*,*)"fmisc.f90::copy: iuuo = ",iuuo     
     write(*,*)"fmisc.f90::copy: llbout = ",llbout
     write(*,*)"fmisc.f90::copy: uubout = ",uubout
     write(*,*)"fmisc.f90::copy: llbin = ",llbin
     write(*,*)"fmisc.f90::copy: uubin = ",uubin
     write(*,*)"fmisc.f90::copy: llbcopy = ",lcopy
     write(*,*)"fmisc.f90::copy: uubcopy = ",ucopy
     stop
   endif

  data_out(illo(1):iuuo(1),illo(2):iuuo(2),illo(3):iuuo(3))=data_in(illi(1):iuui(1),illi(2):iuui(2),illi(3):iuui(3))

  return

  end subroutine copy
!--------------------------------------------------------------------------
! three dimensional interpolation for cell center grid structure  
  subroutine global_interp(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry,ORDN
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA

!~~~~~~> Other parameters:

  integer :: j,m,imin,jmin,kmin
  integer,dimension(3) :: cxB,cxT,cxI,cmin,cmax
  real*8,dimension(3) :: cx
  real*8, dimension(1:ORDN) :: x1a
  real*8, dimension(1:ORDN,1:ORDN,1:ORDN) :: ya
  integer, parameter :: NO_SYMM = 0, EQUATORIAL = 1, OCTANT = 2
  real*8 :: dX,dY,dZ,ddy
  real*8, parameter :: ONE=1.d0
  logical::decide3d

  imin = lbound(f,1)
  jmin = lbound(f,2)
  kmin = lbound(f,3)

  dX = X(imin+1)-X(imin)
  dY = Y(jmin+1)-Y(jmin)
  dZ = Z(kmin+1)-Z(kmin)

  forall( j = 1:ordn ) x1a(j) = ( j - 1 )* ONE

  cxI(1) = idint((x1-X(1))/dX+0.4)+1
  cxI(2) = idint((y1-Y(1))/dY+0.4)+1
  cxI(3) = idint((z1-Z(1))/dZ+0.4)+1

  cxB = cxI - ORDN/2+1
  cxT = cxB + ORDN - 1
       
  cmin = 1
  cmax = ex
  if(Symmetry == OCTANT  .and.dabs(X(1))<dX) cmin(1) = -ORDN/2+1
  if(Symmetry == OCTANT  .and.dabs(Y(1))<dY) cmin(2) = -ORDN/2+1
  if(Symmetry /= NO_SYMM .and.dabs(Z(1))<dZ) cmin(3) = -ORDN/2+1

  do m =1,3
   if(cxB(m) < cmin(m))then
      cxB(m) = cmin(m)
      cxT(m) = cxB(m) + ORDN - 1
   endif
   if(cxT(m) > cmax(m))then
      cxT(m) = cmax(m)
      cxB(m) = cxT(m) + 1 - ORDN
   endif
 enddo
 if(cxB(1)>0)then
  cx(1) = (x1 - X(cxB(1)))/dX
 else
  cx(1) = (x1 + X(1-cxB(1)))/dX
 endif
 if(cxB(2)>0)then
  cx(2) = (y1 - Y(cxB(2)))/dY
 else
  cx(2) = (y1 + Y(1-cxB(2)))/dY
 endif
 if(cxB(3)>0)then
  cx(3) = (z1 - Z(cxB(3)))/dZ
 else
  cx(3) = (z1 + Z(1-cxB(3)))/dZ
 endif

  if(decide3d(ex,f,f,cxB,cxT,SoA,ya,ORDN,Symmetry))then
     write(*,*)"global_interp position: ",x1,y1,z1
     write(*,*)"data range: ",X(1),X(ex(1)),Y(1),Y(ex(2)),Z(1),Z(ex(3))
     stop
  endif
  call polin3(x1a,x1a,x1a,ya,cx(1),cx(2),cx(3),f_int,ddy,ORDN)

  return

  end subroutine global_interp
!----------------------------------------------------------------
! decide which 3d data to be used does not surport PI-Symmetry yet 
!----------------------------------------------------------------
  function decide3d(ex,f,fpi,cxB,cxT,SoA,ya,ORDN,Symmetry)  result(gont)
  implicit none

  integer,                                 intent(in) :: ORDN,Symmetry
  integer,dimension(1:3)                 , intent(in) :: ex,cxB,cxT
  real*8, dimension(1:3)                 , intent(in) :: SoA
  real*8, dimension(ex(1),ex(2),ex(3))   , intent(in) :: f,fpi
  real*8, dimension(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3)), intent(out):: ya
  logical::gont

  integer,dimension(1:3) :: fmin1,fmin2,fmax1,fmax2
  integer::i,j,k,m

  gont=.false.
  do m=1,3
! check cxB and cxT are NaN or not  
    if(.not.(iabs(cxB(m)).ge.0)) gont=.true.
    if(.not.(iabs(cxT(m)).ge.0)) gont=.true.
    fmin1(m) = max(1,cxB(m))
    fmax1(m) = cxT(m)
    fmin2(m) = cxB(m)
    fmax2(m) = min(0,cxT(m))
    if((fmin1(m).le.fmax1(m)).and.(  fmin1(m)<1.or.  fmax1(m)>ex(m)))gont=.true.
    if((fmin2(m).le.fmax2(m)).and.(1-fmax2(m)<1.or.1-fmin2(m)>ex(m)))gont=.true.
  enddo
!sanity check
  if(gont)then
          write(*,*)"error in decide3d"
          write(*,*)((fmin1.le.fmax1).and.(  fmin1<1.or.  fmax1>ex))
          write(*,*)((fmin2.le.fmax2).and.(1-fmax2<1.or.1-fmin2>ex))
          write(*,*)"cxB, cxT and data shape:"
          write(*,*)cxB,cxT,ex
          write(*,*)"resulted fmin1, fmax1 and fmin2, fmax2:"
          write(*,*)fmin1,fmax1,fmin2,fmax2
  else

  do k=fmin1(3),fmax1(3)
     do j=fmin1(2),fmax1(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,j,k)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(1-i,j,k)*SoA(1)
        enddo
     enddo
     do j=fmin2(2),fmax2(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,1-j,k)*SoA(2)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(1-i,1-j,k)*SoA(1)*SoA(2)
        enddo
     enddo
  enddo
           
  do k=fmin2(3),fmax2(3)
     do j=fmin1(2),fmax1(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,j,1-k)*SoA(3)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(1-i,j,1-k)*SoA(1)*SoA(3)
        enddo
     enddo
     do j=fmin2(2),fmax2(2)
        do i=fmin1(1),fmax1(1)
           ya(i,j,k) = f(i,1-j,1-k)*SoA(2)*SoA(3)
        enddo
        do i=fmin2(1),fmax2(1)
           ya(i,j,k) = f(1-i,1-j,1-k)*SoA(1)*SoA(2)*SoA(3)
        enddo
     enddo
  enddo

  endif

  end function decide3d

!---------------------------------------------------------------------------------------
subroutine symmetry_bd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1),-ord+1:extc(2),-ord+1:extc(3)),intent(out):: funcc
  real*8, dimension(1:3), intent(in) :: SoA

  integer::i

!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
   enddo
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
   enddo
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
   do i=0,ord-1
      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
   enddo

end subroutine symmetry_bd

subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,-ord+1:extc(3)+ord),intent(out):: funcc
  real*8, dimension(1:3), intent(in) :: SoA

  integer::i

  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
   enddo
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
      funcc(:,extc(2)+1+i,1:extc(3)) = funcc(:,extc(2)-i,1:extc(3))*SoA(2)
   enddo
   do i=0,ord-1
      funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
      funcc(:,:,extc(3)+1+i) = funcc(:,:,extc(3)-i)*SoA(3)
   enddo

end subroutine symmetry_tbd

subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,extc(3)),intent(out):: funcc
  real*8, dimension(2), intent(in) :: SoA

  integer::i

  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
   enddo
   do i=0,ord-1
      funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
      funcc(:,extc(2)+1+i,1:extc(3)) = funcc(:,extc(2)-i,1:extc(3))*SoA(2)
   enddo

end subroutine symmetry_stbd

subroutine symmetry_sntbd(ord,extc,func,funcc,SoA,actd)
  implicit none

!~~~~~~> input arguments
  integer,intent(in) :: ord,actd
  integer,dimension(3),   intent(in) :: extc
  real*8, dimension(extc(1),extc(2),extc(3)),intent(in ):: func
  real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,extc(3)),intent(out):: funcc
  real*8, intent(in) :: SoA

  integer::i

  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
  if(actd==0)then
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA
      funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA
   enddo
  elseif(actd==1)then
   do i=0,ord-1
      funcc(1:extc(1),-i,1:extc(3)) = funcc(1:extc(1),i+1,1:extc(3))*SoA
      funcc(1:extc(1),extc(2)+1+i,1:extc(3)) = funcc(1:extc(1),extc(2)-i,1:extc(3))*SoA
   enddo
  else
    write(*,*)"symmetry_sntbd: not recognized actd = ",actd
  endif

end subroutine symmetry_sntbd

subroutine d2dump(wei,llb,uub,ext,data_in,data_out,gord,SoA)
  implicit none
  integer,intent(in) :: wei,gord
  integer,dimension(3),intent(in) ::ext
  real*8,dimension(3),intent(in) :: SoA
  real*8,dimension(3) :: llb,uub
  real*8,dimension(ext(1),ext(2),ext(3)),intent(in)::data_in
  real*8,dimension(ext(1),ext(2)),intent(inout)::data_out

  real*8 :: dZ
  integer :: i,j,k

!sanity check
  if(wei.ne.3)then
     write(*,*)"fmisc.f90::copy: this routine only surport 3 dimension"
     write(*,*)"dim = ",wei
     stop
  endif

  dZ = (uub(3)-llb(3))/ext(3)
  k = idint((0-llb(3))/dZ+0.4)+1

  select case (gord)
  case (2)
     if(k > 2)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5625d0*(data_in(i,j,k)+data_in(i,j,k-1))-0.0625d0*(data_in(i,j,k+1)+data_in(i,j,k-2))
         enddo
         enddo
     else if(k == 1)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5625d0*(data_in(i,j,k)+SoA(3)*data_in(i,j,k))-0.0625d0*(data_in(i,j,k+1)+SoA(3)*data_in(i,j,k+1))
         enddo
         enddo
     else
      write(*,*) "d2dump: something must be wrong, k = ",k
      return
     endif
  case (3)
     if(k > 3)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5859375d0*(data_in(i,j,k)+data_in(i,j,k-1))   &
                           -0.9765625d-1*(data_in(i,j,k+1)+data_in(i,j,k-2)) &
                           +0.1171875d-1*(data_in(i,j,k+2)+data_in(i,j,k-3))
         enddo
         enddo
     else if(k == 1)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5859375d0*(data_in(i,j,k)+SoA(3)*data_in(i,j,k))   &
                           -0.9765625d-1*(data_in(i,j,k+1)+SoA(3)*data_in(i,j,k+1)) &
                           +0.1171875d-1*(data_in(i,j,k+2)+SoA(3)*data_in(i,j,k+2))
         enddo
         enddo
     else
      write(*,*) "d2dump: something must be wrong, k = ",k
      return
     endif
  case (4)
     if(k > 4)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5981445312d0*(data_in(i,j,k)+data_in(i,j,k-1))   &
                            -0.1196289063d0*(data_in(i,j,k+1)+data_in(i,j,k-2)) &
                           +0.2392578125d-1*(data_in(i,j,k+2)+data_in(i,j,k-3)) &
                           -0.2441406250d-2*(data_in(i,j,k+3)+data_in(i,j,k-4))
         enddo
         enddo
     else if(k == 1)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.5981445312d0*(data_in(i,j,k)+SoA(3)*data_in(i,j,k))   &
                            -0.1196289063d0*(data_in(i,j,k+1)+SoA(3)*data_in(i,j,k+1)) &
                           +0.2392578125d-1*(data_in(i,j,k+2)+SoA(3)*data_in(i,j,k+2)) &
                           -0.2441406250d-2*(data_in(i,j,k+3)+SoA(3)*data_in(i,j,k+3))
         enddo
         enddo
     else
      write(*,*) "d2dump: something must be wrong, k = ",k
      return
     endif
  case (5)
     if(k > 5)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.6056213378d0*(data_in(i,j,k)+data_in(i,j,k-1))   &
                            -0.1345825196d0*(data_in(i,j,k+1)+data_in(i,j,k-2)) &
                           +0.3460693359d-1*(data_in(i,j,k+2)+data_in(i,j,k-3)) &
                           -0.6179809571d-2*(data_in(i,j,k+3)+data_in(i,j,k-4)) &
                           +0.5340576171d-3*(data_in(i,j,k+4)+data_in(i,j,k-5))
         enddo
         enddo
     else if(k == 1)then
         do i=1,ext(1)
         do j=1,ext(2)
             data_out(i,j) = 0.6056213378d0*(data_in(i,j,k)+SoA(3)*data_in(i,j,k))   &
                            -0.1345825196d0*(data_in(i,j,k+1)+SoA(3)*data_in(i,j,k+1)) &
                           +0.3460693359d-1*(data_in(i,j,k+2)+SoA(3)*data_in(i,j,k+2)) &
                           -0.6179809571d-2*(data_in(i,j,k+3)+SoA(3)*data_in(i,j,k+3)) &
                           +0.5340576171d-3*(data_in(i,j,k+4)+SoA(3)*data_in(i,j,k+4))
         enddo
         enddo
     else
      write(*,*) "d2dump: something must be wrong, k = ",k
      return
     endif
  case default
    write(*,*) "d2dump: not recognized ord = ",gord
    return
  end select

end subroutine d2dump

#else
#error Not define Vertex nor Cell
#endif  
#endif
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! common code for cell and vertex
!------------------------------------------------------------------------------
! Lagrangian polynomial interpolation
!------------------------------------------------------------------------------

!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
  subroutine polint6_neville(xa, ya, x, y, dy)
  implicit none

  real*8, dimension(6), intent(in) :: xa, ya
  real*8, intent(in) :: x
  real*8, intent(out) :: y, dy

  integer :: i, m, ns, n_m
  real*8, dimension(6) :: c, d, ho
  real*8 :: dif, dift, hp, h, den_val

  c = ya
  d = ya
  ho = xa - x

  ns = 1
  dif = abs(x - xa(1))

  do i = 2, 6
    dift = abs(x - xa(i))
    if (dift < dif) then
      ns = i
      dif = dift
    end if
  end do

  y = ya(ns)
  ns = ns - 1

  do m = 1, 5
    n_m = 6 - m
    do i = 1, n_m
      hp = ho(i)
      h  = ho(i+m)
      den_val = hp - h

      if (den_val == 0.0d0) then
        write(*,*) 'failure in polint for point',x
        write(*,*) 'with input points: ',xa
        stop
      end if

      den_val = (c(i+1) - d(i)) / den_val

      d(i) = h * den_val
      c(i) = hp * den_val
    end do

    if (2 * ns < n_m) then
      dy = c(ns + 1)
    else
      dy = d(ns)
      ns = ns - 1
    end if
    y = y + dy
  end do

  return
  end subroutine polint6_neville

!DIR$ ATTRIBUTES FORCEINLINE :: polint
  subroutine polint(xa, ya, x, y, dy, ordn)
  implicit none

  integer, intent(in) :: ordn
  real*8, dimension(ordn), intent(in) :: xa, ya
  real*8, intent(in) :: x
  real*8, intent(out) :: y, dy

  integer :: i, m, ns, n_m
  real*8, dimension(ordn) :: c, d, ho
  real*8 :: dif, dift, hp, h, den_val

  if (ordn == 6) then
    call polint6_neville(xa, ya, x, y, dy)
    return
  end if

  c = ya
  d = ya
  ho = xa - x

  ns = 1
  dif = abs(x - xa(1))

  do i = 2, ordn
    dift = abs(x - xa(i))
    if (dift < dif) then
      ns = i
      dif = dift
    end if
  end do

  y = ya(ns)
  ns = ns - 1

  do m = 1, ordn - 1
    n_m = ordn - m
    do i = 1, n_m
      hp = ho(i)
      h  = ho(i+m)
      den_val = hp - h

      if (den_val == 0.0d0) then
        write(*,*) 'failure in polint for point',x
        write(*,*) 'with input points: ',xa
        stop
      end if

      den_val = (c(i+1) - d(i)) / den_val

      d(i) = h * den_val
      c(i) = hp * den_val
    end do

    if (2 * ns < n_m) then
      dy = c(ns + 1)
    else
      dy = d(ns)
      ns = ns - 1
    end if
    y = y + dy
  end do

  return
  end subroutine polint
!------------------------------------------------------------------------------
! Compute Lagrange interpolation basis weights for one target point.
!------------------------------------------------------------------------------
!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
  subroutine polint_lagrange_weights(xa, x, w, ordn)
  implicit none

  integer, intent(in) :: ordn
  real*8, dimension(1:ordn), intent(in) :: xa
  real*8, intent(in) :: x
  real*8, dimension(1:ordn), intent(out) :: w

  integer :: i, j
  real*8 :: num, den, dx

  do i = 1, ordn
    num = 1.d0
    den = 1.d0
    do j = 1, ordn
      if (j /= i) then
        dx = xa(i) - xa(j)
        if (dx == 0.0d0) then
          write(*,*) 'failure in polint for point',x
          write(*,*) 'with input points: ',xa
          stop
        end if
        num = num * (x - xa(j))
        den = den * dx
      end if
    end do
    w(i) = num / den
  end do

  return
  end subroutine polint_lagrange_weights
!------------------------------------------------------------------------------
!
! interpolation in 2 dimensions, follow yx order
!
!------------------------------------------------------------------------------
  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
  implicit none

  integer,intent(in) :: ordn
  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
  real*8, intent(in) :: x1,x2
  real*8, intent(out) :: y,dy

#ifdef POLINT_LEGACY_ORDER
  integer  :: i,m
  real*8, dimension(ordn) :: ymtmp
  real*8, dimension(ordn) :: yntmp

  m=size(x1a)
  do i=1,m
    yntmp=ya(i,:)
    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
  end do
  call polint(x1a,ymtmp,x1,y,dy,ordn)
#else
  integer  :: j
  real*8, dimension(ordn) :: ymtmp
  real*8 :: dy_temp

  do j=1,ordn
    call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
  end do
  call polint(x2a, ymtmp, x2, y, dy, ordn)
#endif

  return
  end subroutine polin2
!------------------------------------------------------------------------------
!
! interpolation in 3 dimensions, follow zyx order
!
!------------------------------------------------------------------------------
  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
  implicit none

  integer,intent(in) :: ordn
  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
  real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
  real*8, intent(in) :: x1,x2,x3
  real*8, intent(out) :: y,dy

#ifdef POLINT_LEGACY_ORDER
  integer  :: i,j,m,n
  real*8, dimension(ordn,ordn) :: yatmp
  real*8, dimension(ordn) :: ymtmp
  real*8, dimension(ordn) :: yntmp
  real*8, dimension(ordn) :: yqtmp

  m=size(x1a)
  n=size(x2a)
  do i=1,m
   do j=1,n
    yqtmp=ya(i,j,:)
    call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
   end do
    yntmp=yatmp(i,:)
    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
  end do
  call polint(x1a,ymtmp,x1,y,dy,ordn)
#else
  integer  :: i, j, k
  real*8, dimension(ordn) :: w1, w2
  real*8, dimension(ordn) :: ymtmp
  real*8 :: yx_sum, x_sum

  call polint_lagrange_weights(x1a, x1, w1, ordn)
  call polint_lagrange_weights(x2a, x2, w2, ordn)

  do k = 1, ordn
    yx_sum = 0.d0
    do j = 1, ordn
      x_sum = 0.d0
      do i = 1, ordn
        x_sum = x_sum + w1(i) * ya(i,j,k)
      end do
      yx_sum = yx_sum + w2(j) * x_sum
    end do
    ymtmp(k) = yx_sum
  end do

  call polint(x3a, ymtmp, x3, y, dy, ordn)
#endif

  return
  end subroutine polin3
!--------------------------------------------------------------------------------------
! calculate L2norm
  subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                          f,f_out,gw)

  implicit none
!~~~~~~> Input parameters:
  integer,intent(in ):: ex(1:3)
  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
  integer,intent(in)::gw
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8, intent(out) :: f_out
!~~~~~~> Other variables:

  real*8, parameter :: ZEO = 0.D0
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
  integer::i,j,k,n_elements
  real*8, dimension(:), allocatable :: f_flat
  real*8, external :: DDOT

  dX = X(2) - X(1)
  dY = Y(2) - Y(1)
  dZ = Z(2) - Z(1)

! for ghost zone
   imin = gw+1
   jmin = gw+1
   kmin = gw+1

   imax = ex(1) - gw
   jmax = ex(2) - gw
   kmax = ex(3) - gw

!for patch boundary (i.e., not ghost boundary)

if(dabs(X(ex(1))-xmax) < dX) imax = ex(1)
if(dabs(Y(ex(2))-ymax) < dY) jmax = ex(2)
if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)
if(dabs(X(1)-xmin) < dX) imin = 1
if(dabs(Y(1)-ymin) < dY) jmin = 1
if(dabs(Z(1)-zmin) < dZ) kmin = 1

! Optimized with oneMKL BLAS DDOT for dot product
n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(n_elements))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
deallocate(f_flat)

f_out = f_out*dX*dY*dZ

  return

  end subroutine l2normhelper
!--------------------------------------------------------------------------------------
! calculate L2norm especially for shell Blocks
  subroutine l2normhelper_sh(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                          f,f_out,gw,ogw,Symmetry)

  implicit none
!~~~~~~> Input parameters:
  integer,intent(in ):: ex(1:3),Symmetry
  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
  integer,intent(in)::gw,ogw
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8, intent(out) :: f_out
!~~~~~~> Other variables:

  real*8, parameter :: ZEO = 0.D0
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
  integer::i,j,k,n_elements
  real*8, dimension(:), allocatable :: f_flat
  real*8, external :: DDOT

  real*8 :: PIo4

  PIo4 = dacos(-1.d0)/4.d0

  dX = X(2) - X(1)
  dY = Y(2) - Y(1)
  dZ = Z(2) - Z(1)

! for ghost zone
   imin = gw+1
   jmin = gw+1
   kmin = gw+1

   imax = ex(1) - gw
   jmax = ex(2) - gw
   kmax = ex(3) - gw

!for patch boundary (i.e., not ghost boundary)

if(dabs(X(ex(1))-xmax) < dX)then
  if(X(ex(1))-PIo4 > dX)then
    imax = ex(1)-ogw             ! for overlap zone
  else
    imax = ex(1)
  endif
endif
if(dabs(Y(ex(2))-ymax) < dY)then
  if(Y(ex(2))-PIo4 > dY)then
    jmax = ex(2)-ogw             ! for overlap zone
  else
    jmax = ex(2)
  endif
endif
if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)

if(dabs(X(1)-xmin) < dX)then
  if(X(1)+PIo4 < dX)then
    imin = 1+ogw             ! for overlap zone
  else
    imin = 1
  endif
endif
if(dabs(Y(1)-ymin) < dY)then
  if(Y(1)+PIo4 < dY)then
    jmin = 1+ogw             ! for overlap zone
  else
    jmin = 1
  endif
endif
if(dabs(Z(1)-zmin) < dZ) kmin = 1

!for Symmetry ghost points
if(Symmetry==1)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
  if(dabs(ymax-gw*dY)<dY.and.Y(ex(2))>0.d0) jmax = ex(2)-gw
endif
if(Symmetry==2)then
  if(dabs(xmin+gw*dX)<dX.and.X(1)<0.d0) imin = gw+1
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
endif

! Optimized with oneMKL BLAS DDOT for dot product
n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(n_elements))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
deallocate(f_flat)

f_out = f_out*dX*dY*dZ

  return

  end subroutine l2normhelper_sh
!--------------------------------------------------------------------------------------
! calculate L2norm especially for shell Blocks
! use root mean sqrt method
  subroutine l2normhelper_sh_rms(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                          f,f_out,gw,ogw,Symmetry,Nout)

  implicit none
!~~~~~~> Input parameters:
  integer,intent(in ):: ex(1:3),Symmetry
  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
  integer,intent(in)::gw,ogw
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8, intent(out) :: f_out
  integer,intent(out) :: Nout
!~~~~~~> Other variables:

  real*8, parameter :: ZEO = 0.D0
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
  integer::i,j,k
  real*8, dimension(:), allocatable :: f_flat
  real*8, external :: DDOT

  real*8 :: PIo4

  PIo4 = dacos(-1.d0)/4.d0

  dX = X(2) - X(1)
  dY = Y(2) - Y(1)
  dZ = Z(2) - Z(1)

! for ghost zone
   imin = gw+1
   jmin = gw+1
   kmin = gw+1

   imax = ex(1) - gw
   jmax = ex(2) - gw
   kmax = ex(3) - gw

!for patch boundary (i.e., not ghost boundary)

if(dabs(X(ex(1))-xmax) < dX)then
  if(X(ex(1))-PIo4 > dX)then
    imax = ex(1)-ogw             ! for overlap zone
  else
    imax = ex(1)
  endif
endif
if(dabs(Y(ex(2))-ymax) < dY)then
  if(Y(ex(2))-PIo4 > dY)then
    jmax = ex(2)-ogw             ! for overlap zone
  else
    jmax = ex(2)
  endif
endif
if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)

if(dabs(X(1)-xmin) < dX)then
  if(X(1)+PIo4 < dX)then
    imin = 1+ogw             ! for overlap zone
  else
    imin = 1
  endif
endif
if(dabs(Y(1)-ymin) < dY)then
  if(Y(1)+PIo4 < dY)then
    jmin = 1+ogw             ! for overlap zone
  else
    jmin = 1
  endif
endif
if(dabs(Z(1)-zmin) < dZ) kmin = 1

!for Symmetry ghost points
if(Symmetry==1)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
  if(dabs(ymax-gw*dY)<dY.and.Y(ex(2))>0.d0) jmax = ex(2)-gw
endif
if(Symmetry==2)then
  if(dabs(xmin+gw*dX)<dX.and.X(1)<0.d0) imin = gw+1
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
endif

! Optimized with oneMKL BLAS DDOT for dot product
Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(Nout))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [Nout])
f_out = DDOT(Nout, f_flat, 1, f_flat, 1)
deallocate(f_flat)

  return

  end subroutine l2normhelper_sh_rms
! locating the position of NaN
  subroutine ScaneNaN(ext,X,Y,Z,f)
  implicit none
  integer,dimension(3),intent(in) :: ext
  real*8,dimension(ext(1)) :: X
  real*8,dimension(ext(2)) :: Y
  real*8,dimension(ext(3)) :: Z
  real*8,dimension(ext(1),ext(2),ext(3)) :: f

  integer :: i,j,k

  do k=1,ext(3)
  do j=1,ext(2)
  do i=1,ext(1)
     if(abs(f(i,j,k)) .ne. abs(f(i,j,k))) write(*,*)X(i),Y(j),Z(k),f(i,j,k)
  enddo
  enddo
  enddo

  end subroutine ScaneNaN
! fortran version writefile
  subroutine fwritefile(time,nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax,NN,filename,data_out)
  implicit none

  real*8,intent(in) :: time,xmin,xmax,ymin,ymax,zmin,zmax
  integer,intent(in) :: nx,ny,nz,NN
  real*8,dimension(nx,ny,nz),intent(in) :: data_out
  Character(Len=NN) :: filename

 
!  Open( 12 , File = filename,form='BINARY', access="SEQUENTIAL",status="replace",action='Write')
  Open( 12 , File = filename,form='UNFORMATTED', access="DIRECT",status="replace",action='Write')
  Write( 12 ) time,nx,ny,nz,xmin,xmax,ymin,ymax,zmin,zmax,data_out

  Close( 12 )

  end subroutine fwritefile
!--------------------------------------------------------------------------
!
! average for interface 
!
!--------------------------------------------------------------------------
  subroutine average(ext,f1,f2,fout)
  implicit none
  integer,dimension(3),   intent(in) :: ext
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in):: f1,f2
  real*8, dimension(ext(1),ext(2),ext(3)),intent(out):: fout

  real*8,parameter::HLF=0.5d0

  fout = HLF*(f1+f2)

  return

  end subroutine average
!-----------------------------------------------------------------------------  
  subroutine average3(ext,f1,f2,fout)
  implicit none
  integer,dimension(3),   intent(in) :: ext
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in):: f1,f2
  real*8, dimension(ext(1),ext(2),ext(3)),intent(out):: fout
! f1 ----------                ^
!                fout ------p  | t
! f2 ----------                |
! 2 points, 1st order interpolation
! 1   2
! f2  f1
! *---*--> t
!    ^
! f=3/4*f_1 + 1/4*f_2

  real*8,parameter::C1=0.75d0,C2=0.25d0

  fout = C1*f1+C2*f2

  return

  end subroutine average3
!-----------------------------------------------------------------------------  
  subroutine average2(ext,f1,f2,f3,fout)
  implicit none
  integer,dimension(3),   intent(in) :: ext
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in):: f1,f2,f3
  real*8, dimension(ext(1),ext(2),ext(3)),intent(out):: fout
! f1 ----------                ^
!                fout ------   |
! f2 ----------                | t
!                              |
! f3 ----------                |
! 3 points, 2nd order interpolation
! 1   2   3
! f3  f2  f1
! *---*---*--> t
!       ^
! f=3/8*f_1 + 3/4*f_2 - 1/8*f_3

  real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
  integer :: i,j,k

  do concurrent (k=1:ext(3), j=1:ext(2), i=1:ext(1))
    fout(i,j,k) = C1*f1(i,j,k)+C2*f2(i,j,k)+C3*f3(i,j,k)
  end do

  return

  end subroutine average2
!-----------------------------------------------------------------------------  
  subroutine average2p(ext,f1,f2,f3,fout)
  implicit none
  integer,dimension(3),   intent(in) :: ext
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in):: f1,f2,f3
  real*8, dimension(ext(1),ext(2),ext(3)),intent(out):: fout
! f1 ----------                ^
!                fout ------p  |
! f2 ----------                | t
!                              |
! f3 ----------                |
! 3 points, 2nd order interpolation
! 1   2   3
! f3  f2  f1
! *---*---*--> t
!        ^
! f=21/32*f_1 + 7/16*f_2 - 3/32*f_3

  real*8,parameter::C1=5.d0/3.2d1,C2=1.5d1/1.6d1,C3=-3.d0/3.2d1

  fout = C1*f1+C2*f2+C3*f3

  return

  end subroutine average2p
!-----------------------------------------------------------------------------  
  subroutine average2m(ext,f1,f2,f3,fout)
  implicit none
  integer,dimension(3),   intent(in) :: ext
  real*8, dimension(ext(1),ext(2),ext(3)),intent(in):: f1,f2,f3
  real*8, dimension(ext(1),ext(2),ext(3)),intent(out):: fout
! f1 ----------                ^
!                fout ------m  |
! f2 ----------                | t
!                              |
! f3 ----------                |
! 3 points, 2nd order interpolation
! 1   2   3
! f3  f2  f1
! *---*---*--> t
!      ^
! f=5/32*f_1 + 15/16*f_2 - 3/32*f_3

  real*8,parameter::C1=5.d0/3.2d1,C2=1.5d1/1.6d1,C3=-3.d0/3.2d1

  fout = C1*f1+C2*f2+C3*f3

  return

  end subroutine average2m
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  subroutine lowerboundset(ex,chi0,TINNY)
  implicit none

!~~~~~~% Input parameters:

  integer ,intent(in):: ex(1:3)
  real*8  ,intent(in):: TINNY
  real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) ::chi0

  where(chi0 < TINNY) chi0 = TINNY

  return   

  end subroutine lowerboundset
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!global interpolation with given index and coeffients
  subroutine global_interpind(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry,inds,coef,sst)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry,ORDN,sst
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA
  integer,dimension(3),                intent(in) :: inds
  real*8, dimension(3*ORDN),           intent(in) :: coef

!~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,-ORDN+1:ex(3)+ORDN) :: fh
  integer :: m
  integer,dimension(3) :: cxB,cxT
  real*8, dimension(ORDN,ORDN,ORDN) :: ya
  real*8, dimension(ORDN,ORDN) :: tmp2
  real*8, dimension(ORDN) :: tmp1
  real*8, dimension(3) :: SoAh
  real*8, external :: DDOT

! +1 because c++ gives 0 for first point
  cxB = inds+1  
  cxT = cxB + ORDN - 1
      
  if(all(cxB>0).and.all(cxT<ex+1))then
     ya=f(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3))
  elseif(any(cxB<-ORDN+1).or.any(cxT>ex+ORDN))then
     write(*,*)"error in global_interpind, cxB = ",cxB
     write(*,*)"                           cxT = ",cxT 
     write(*,*)"                           ext = ",ex
     stop
  else
     if(sst==-1)then
        SoAh = SoA
        if(any(cxT>ex)) write(*,*)"error global_interpind sst =",sst
     elseif(sst==0.or.sst==1)then
        SoAh = SoA
        SoAh(3) = 0
        if(cxB(3)<1.or.cxT(3)>ex(3)) write(*,*)"error global_interpind sst =",sst
     elseif(sst==2.or.sst==3)then
        SoAh(1) = SoA(2)
        SoAh(2) = SoA(3)
        SoAh(3) = 0
        if(cxB(3)<1.or.cxT(3)>ex(3)) write(*,*)"error global_interpind sst =",sst
     elseif(sst==4.or.sst==5)then
        SoAh(1) = SoA(1)
        SoAh(2) = SoA(3)
        SoAh(3) = 0
        if(cxB(3)<1.or.cxT(3)>ex(3)) write(*,*)"error global_interpind sst =",sst,cxB(3),cxT(3)
     endif
     call symmetry_tbd(ORDN,ex,f,fh,SoAh)
     ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3))
  endif 

  ! Optimized with BLAS operations for better performance
  ! First dimension: z-direction weighted sum
  tmp2=0
  do m=1,ORDN
    tmp2 = tmp2 + coef(2*ORDN+m)*ya(:,:,m)
  enddo

  ! Second dimension: y-direction weighted sum
  tmp1=0
  do m=1,ORDN
    tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
  enddo

  ! Third dimension: x-direction weighted sum using BLAS DDOT
  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)

  return

  end subroutine global_interpind
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!global interpolation with given index and coeffients
! special for shell to shell
  subroutine global_interpind2d(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry,inds,coef,sst)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry,ORDN,sst
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA
  integer,dimension(3),                intent(in) :: inds
  real*8, dimension(2*ORDN),           intent(in) :: coef

!~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,ex(3)) :: fh
  integer :: m
  integer,dimension(2) :: cxB,cxT
  real*8, dimension(ORDN,ORDN) :: ya
  real*8, dimension(ORDN) :: tmp1
  real*8, dimension(2) :: SoAh
  real*8, external :: DDOT

! +1 because c++ gives 0 for first point
  cxB = inds(1:2)+1  
  cxT = cxB + ORDN - 1
      
  if(all(cxB>0).and.all(cxT<ex(1:2)+1))then
     ya=f(cxB(1):cxT(1),cxB(2):cxT(2),inds(3))
  elseif(any(cxB<-ORDN+1).or.any(cxT>ex(1:2)+ORDN))then
     write(*,*)"error in global_interpind2d, cxB = ",cxB
     write(*,*)"                             cxT = ",cxT 
     write(*,*)"                             ext = ",ex(1:2)
     stop
  else
     if(sst==-1)then
       write(*,*)"error in global_interpind2d, sst = ",sst
       stop
     elseif(sst==0.or.sst==1)then
        SoAh = SoA(1:2)
     elseif(sst==2.or.sst==3)then
        SoAh(1) = SoA(2)
        SoAh(2) = SoA(3)
     elseif(sst==4.or.sst==5)then
        SoAh(1) = SoA(1)
        SoAh(2) = SoA(3)
     endif
     call symmetry_stbd(ORDN,ex,f,fh,SoAh)
     ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),inds(3))
  endif 

  ! Optimized with BLAS operations
  tmp1=0
  do m=1,ORDN
    tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
  enddo

  ! Use BLAS DDOT for final weighted sum
  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)

  return

  end subroutine global_interpind2d
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!global interpolation with given index and coeffients
! special for shell to shell
! dumyd refer to source
  subroutine global_interpind1d(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry,indsi,coef,sst,dumyd)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3),symmetry,ORDN,sst,dumyd
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA
  integer,dimension(3),                intent(in) :: indsi
  real*8, dimension(ORDN),             intent(in) :: coef

!~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,ex(3)) :: fh
  integer :: m
  integer :: cxB,cxT
  real*8, dimension(ORDN) :: ya
  real*8 :: SoAh
  integer,dimension(3) :: inds
  real*8, external :: DDOT

! +1 because c++ gives 0 for first point
  inds = indsi + 1
  cxB = inds(1)  
  cxT = cxB + ORDN - 1
     
! active is rho  
  if(dumyd==1)then

  if(cxB>0.and.cxT<ex(1)+1)then
     ya=f(cxB:cxT,inds(2),inds(3))
  elseif(cxB<-ORDN+1.or.cxT>ex(1)+ORDN)then
     write(*,*)"error in global_interpind1d, cxB = ",cxB
     write(*,*)"                             cxT = ",cxT 
     write(*,*)"                             ext = ",ex(1)
     stop
  else
     if(sst==-1)then
       write(*,*)"error in global_interpind1d, sst = ",sst
       stop
     elseif(sst==0.or.sst==1)then
        SoAh = SoA(1)
     elseif(sst==2.or.sst==3)then
        SoAh = SoA(2)
     elseif(sst==4.or.sst==5)then
        SoAh = SoA(1)
     endif
     call symmetry_sntbd(ORDN,ex,f,fh,SoAh,1-dumyd)
     ya=fh(cxB:cxT,inds(2),inds(3))
  endif 

! active is sigma
  elseif(dumyd==0)then

  if(cxB>0.and.cxT<ex(2)+1)then
     ya=f(inds(2),cxB:cxT,inds(3))
  elseif(cxB<-ORDN+1.or.cxT>ex(2)+ORDN)then
     write(*,*)"error in global_interpind1d, cxB = ",cxB
     write(*,*)"                             cxT = ",cxT 
     write(*,*)"                             ext = ",ex(2)
     stop
  else
     if(sst==-1)then
       write(*,*)"error in global_interpind1d, sst = ",sst
       stop
     elseif(sst==0.or.sst==1)then
        SoAh = SoA(2)
     elseif(sst==2.or.sst==3)then
        SoAh = SoA(3)
     elseif(sst==4.or.sst==5)then
        SoAh = SoA(3)
     endif
     call symmetry_sntbd(ORDN,ex,f,fh,SoAh,1-dumyd)
     ya=fh(inds(2),cxB:cxT,inds(3))
  endif 

  else
          write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
  endif

  ! Optimized with BLAS DDOT for weighted sum
  f_int = DDOT(ORDN, coef, 1, ya, 1)

  return

  end subroutine global_interpind1d
!-----------------------------------------------------------------------------------------------------------------
! three dimensional interpolation for both vertex and cell center grid structure  
! for distinguishing shell and Cartesian
  subroutine global_interp_ss(ex,X,Y,Z,f,f_int,x1,y1,z1,ORDN,SoA,symmetry,sst)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3), symmetry,ORDN,sst
  real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1,z1
  real*8, dimension(3),                intent(in) :: SoA

!~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,-ORDN+1:ex(3)+ORDN) :: fh
  real*8, dimension(3) :: SoAh
  integer :: j,m,imin,jmin,kmin
  integer,dimension(3) :: cxB,cxT,cxI,cmin,cmax
  real*8,dimension(3) :: cx
  real*8, dimension(1:ORDN) :: x1a
  real*8, dimension(1:ORDN,1:ORDN,1:ORDN) :: ya
  integer, parameter :: NO_SYMM = 0, EQUATORIAL = 1, OCTANT = 2
  real*8 :: dX,dY,dZ,ddy
  real*8, parameter :: ONE=1.d0

  imin = lbound(f,1)
  jmin = lbound(f,2)
  kmin = lbound(f,3)

  dX = X(imin+1)-X(imin)
  dY = Y(jmin+1)-Y(jmin)
  dZ = Z(kmin+1)-Z(kmin)

  forall( j = 1:ordn ) x1a(j) = ( j - 1 )* ONE

  cxI(1) = idint((x1-X(1))/dX+0.4)+1
  cxI(2) = idint((y1-Y(1))/dY+0.4)+1
  cxI(3) = idint((z1-Z(1))/dZ+0.4)+1

  cxB = cxI - ORDN/2+1
  cxT = cxB + ORDN - 1

  cmin = 1
  cmax = ex

  if(sst==-1)then
     SoAh = SoA
     cmin = -ORDN+1
  elseif(sst==0.or.sst==1)then
     SoAh = SoA
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  elseif(sst==2.or.sst==3)then
     SoAh(1) = SoA(2)
     SoAh(2) = SoA(3)
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  elseif(sst==4.or.sst==5)then
     SoAh(1) = SoA(1)
     SoAh(2) = SoA(3)
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  endif
  do m =1,3
   if(cxB(m) < cmin(m))then
      cxB(m) = cmin(m)
      cxT(m) = cxB(m) + ORDN - 1
   endif
   if(cxT(m) > cmax(m))then
      cxT(m) = cmax(m)
      cxB(m) = cxT(m) + 1 - ORDN
   endif
 enddo
 cx(1) = (x1 - X(1))/dX-cxB(1)+1
 cx(2) = (y1 - Y(1))/dY-cxB(2)+1
 cx(3) = (z1 - Z(1))/dZ-cxB(3)+1

  call symmetry_tbd(ORDN,ex,f,fh,SoAh)
  ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3))

  call polin3(x1a,x1a,x1a,ya,cx(1),cx(2),cx(3),f_int,ddy,ORDN)

  return

  end subroutine global_interp_ss
!-----------------------------------------------------------------------------------------------------------------
! two dimensional interpolation for both vertex and cell center grid structure  
! for distinguishing shell and Cartesian
  subroutine global_interp_ss_2d(ex,X,Y,indZ,f,f_int,x1,y1,ORDN,SoA,symmetry,sst)
  implicit none

!~~~~~~> Input parameters:

  integer,                             intent(in) :: ex(1:3),indZ,symmetry,ORDN,sst
  real*8,intent(in) :: X(ex(1)),Y(ex(2))
  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f
  real*8,                              intent(out):: f_int
  real*8,                              intent(in) :: x1,y1
  real*8, dimension(3),                intent(in) :: SoA

!~~~~~~> Other parameters:

  real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,-ORDN+1:ex(3)+ORDN) :: fh
  real*8, dimension(3) :: SoAh
  integer :: j,m,imin,jmin,kmin
  integer,dimension(2) :: cxB,cxT,cxI,cmin,cmax
  real*8,dimension(2) :: cx
  real*8, dimension(1:ORDN) :: x1a
  real*8, dimension(1:ORDN,1:ORDN) :: ya
  integer, parameter :: NO_SYMM = 0, EQUATORIAL = 1, OCTANT = 2
  real*8 :: dX,dY,ddy
  real*8, parameter :: ONE=1.d0

! sanity check
  if(indZ < 1 .or. indZ > ex(3))then
      write(*,*)"error in global_interp_ss_2d, ext = ",ex(3),"ind = ",indZ
      return
  endif

  imin = lbound(f,1)
  jmin = lbound(f,2)
  kmin = lbound(f,3)

  dX = X(imin+1)-X(imin)
  dY = Y(jmin+1)-Y(jmin)

  forall( j = 1:ordn ) x1a(j) = ( j - 1 )* ONE

  cxI(1) = idint((x1-X(1))/dX+0.4)+1
  cxI(2) = idint((y1-Y(1))/dY+0.4)+1

  cxB = cxI - ORDN/2+1
  cxT = cxB + ORDN - 1

  cmin = 1
  cmax = ex(1:2)

  if(sst==-1)then
     SoAh = SoA
     cmin = -ORDN+1
  elseif(sst==0.or.sst==1)then
     SoAh = SoA
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  elseif(sst==2.or.sst==3)then
     SoAh(1) = SoA(2)
     SoAh(2) = SoA(3)
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  elseif(sst==4.or.sst==5)then
     SoAh(1) = SoA(1)
     SoAh(2) = SoA(3)
     SoAh(3) = 0
     cmin(1:2) = -ORDN+1
     cmax(1:2) = ex(1:2)+ORDN
  endif
  do m =1,2
   if(cxB(m) < cmin(m))then
      cxB(m) = cmin(m)
      cxT(m) = cxB(m) + ORDN - 1
   endif
   if(cxT(m) > cmax(m))then
      cxT(m) = cmax(m)
      cxB(m) = cxT(m) + 1 - ORDN
   endif
 enddo
 cx(1) = (x1 - X(1))/dX-cxB(1)+1
 cx(2) = (y1 - Y(1))/dY-cxB(2)+1

  call symmetry_tbd(ORDN,ex,f,fh,SoAh)
  ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),indZ)

  call polin2(x1a,x1a,ya,cx(1),cx(2),f_int,ddy,ORDN)

  return

  end subroutine global_interp_ss_2d
!------------------------------------------
!fortran version of Wigner d function
!Eq.(42) of PRD 77, 024027 (2008)
!we consider only theta in [0,pi]
!------------------------------------------
  function fWigner_d_function(l,m,s,costheta) result(gont)
  implicit none
  integer,intent(in) :: l,m,s
  real*8,intent(in) :: costheta

  real*8 :: gont

  integer :: t,C1,C2
  real*8 :: ffact,vv,sinht,cosht

  C1=max(0,m-s)
  C2=min(l+m,l-s)
  vv=0
  sinht=dsqrt((1.d0-costheta)/2.d0)
  cosht=dsqrt((1.d0+costheta)/2.d0);
  if(C1/2*2==C1)then
    do t=C1,C2,2
       vv=vv+cosht**(2*l+m-s-2*t)*sinht**(2*t+s-m)/(ffact(l+m-t)*ffact(l-s-t)*ffact(t)*ffact(t+s-m))
    enddo
    do t=C1+1,C2,2
       vv=vv-cosht**(2*l+m-s-2*t)*sinht**(2*t+s-m)/(ffact(l+m-t)*ffact(l-s-t)*ffact(t)*ffact(t+s-m))
    enddo
  else
    do t=C1,C2,2
       vv=vv-cosht**(2*l+m-s-2*t)*sinht**(2*t+s-m)/(ffact(l+m-t)*ffact(l-s-t)*ffact(t)*ffact(t+s-m))
    enddo
    do t=C1+1,C2,2
       vv=vv+cosht**(2*l+m-s-2*t)*sinht**(2*t+s-m)/(ffact(l+m-t)*ffact(l-s-t)*ffact(t)*ffact(t+s-m))
    enddo
  endif
  
  gont = vv*dsqrt(ffact(l+m)*ffact(l-m)*ffact(l+s)*ffact(l-s))

  return

  end function fWigner_d_function
!----------------------------------
! Optimized factorial function using lookup table for small N
! and log-gamma for large N to avoid overflow
  function ffact(N) result(gont)
  implicit none
  integer,intent(in) :: N

  real*8 :: gont
  integer :: i

  ! Lookup table for factorials 0! to 20! (precomputed)
  real*8, parameter, dimension(0:20) :: fact_table = [ &
    1.d0, 1.d0, 2.d0, 6.d0, 24.d0, 120.d0, 720.d0, 5040.d0, 40320.d0, &
    362880.d0, 3628800.d0, 39916800.d0, 479001600.d0, 6227020800.d0, &
    87178291200.d0, 1307674368000.d0, 20922789888000.d0, &
    355687428096000.d0, 6402373705728000.d0, 121645100408832000.d0, &
    2432902008176640000.d0 ]

! sanity check
  if(N < 0)then
     write(*,*) "ffact: error input for factorial"
     gont = 1.d0
     return
  endif

  ! Use lookup table for small N (fast path)
  if(N <= 20)then
     gont = fact_table(N)
  else
     ! Use log-gamma function for large N: N! = exp(log_gamma(N+1))
     ! This avoids overflow and is computed efficiently
     gont = exp(log_gamma(dble(N+1)))
  endif

  return

  end function ffact
!---------------------------  
!Eq.(41) of PRD 77, 024027 (2008)
!----------------------------------
  function Yslm(s,l,m,the,phi) result(gont)
  implicit none
  integer,intent(in) :: s,l,m
  real*8,intent(in) :: the,phi

  double complex :: gont

  real*8 :: fWigner_d_function,PI,rp

  PI = dacos(-1.d0)

  rp = fWigner_d_function(l,m,s,dcos(the))
  rp = rp*dsqrt((2*l+1.d0)/4.d0/PI)
  if(s/2*2.ne.s) rp = -rp

  gont = dcmplx(dcos(m*phi),dsin(m*phi))

  gont = rp*gont

  return

  end function Yslm
!------------------------------------------------------------------------------------  
subroutine set_value(ext,data_out,rr)

  IMPLICIT NONE

  integer, intent(in) :: ext(3)
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(out) :: data_out
  REAL*8, intent(in) :: rr

  data_out = rr

  return
end subroutine set_value
subroutine add_value(ext,data_out,rr)

  IMPLICIT NONE

  integer, intent(in) :: ext(3)
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(inout) :: data_out
  REAL*8, intent(in) :: rr

  data_out = data_out + rr

  return
end subroutine add_value
! copy array2 to array1  
subroutine array_copy(ext,data1,data2)

  IMPLICIT NONE

  integer, intent(in) :: ext(3)
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(out) :: data1
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(in) :: data2

  data1 = data2

  return
  end subroutine array_copy
! add array2 to array1  
subroutine array_add(ext,data1,data2)

  IMPLICIT NONE

  integer, intent(in) :: ext(3)
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(inout) :: data1
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(in) :: data2

  data1 = data1 + data2

  return
  end subroutine array_add
! subtract array2 from array1  
subroutine array_subtract(ext,data1,data2)

  IMPLICIT NONE

  integer, intent(in) :: ext(3)
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(inout) :: data1
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(in) :: data2

  data1 = data1 - data2

  return
  end subroutine array_subtract
! find out the maximum
subroutine find_maximum(ext,X,Y,Z,fun,val,pos,llb,uub)

  implicit none

  integer,intent(in) :: ext(3),llb(3),uub(3)
  real*8 :: X(ext(1)),Y(ext(2)),Z(ext(3))
  REAL*8, DIMENSION(ext(1),ext(2),ext(3)), intent(in) :: fun
  real*8,intent(out) :: val,pos(3)

  integer :: i,j,k,ii,jj,kk
  real*8 :: tmp

  tmp = 0.d0

  ii=1
  jj=1
  kk=1

  do k=llb(3)+1,ext(3)-uub(3)
  do j=llb(2)+1,ext(2)-uub(2)
  do i=llb(1)+1,ext(1)-uub(1)
      if(dabs(fun(i,j,k)) > tmp)then
         tmp = dabs(fun(i,j,k))
         ii = i
         jj = j
         kk = k
      endif
  enddo
  enddo
  enddo

  pos(1) = X(ii)
  pos(2) = Y(jj)
  pos(3) = Z(kk)
  val = tmp

  return

end subroutine