! -----------------------------------------------------------------------------
! Frame   : Procedures and program used in JGR data report 
!           "Use of twenty years CLUSTER/FGM data to observe the mean behavior 
!           of the magnetic field and current density of Earth's magnetosphere" 
!           By P. Robert and M. Dunlop, October 2021
!
! Library : read_3D_grid.f90
! Object  : routine for reading a 3D grid
! Author  : P. Robert, LPP-ScientiDev, 2020-Oct.2021
! Mail    : patrick.robert@lpp.polytechnique.fr
! -----------------------------------------------------------------------------

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
  
  module dim_grid
  
  real(kind=4), dimension(:,:,:), allocatable   :: grid_Vx, grid_Vy, grid_Vz, grid_V0
  real(kind=4), dimension(:,:,:), allocatable   :: gri_div, gri_dsc, gri_abj
  integer     , dimension(:,:,:), allocatable   :: nval
  integer     , dimension(:,:),   allocatable   :: nval_2D
  integer     , dimension(:,:),   allocatable   :: imagexz,imagexy,imageyz
  
  real(kind=4), dimension(:,:),   allocatable   :: grid_2D
  real(kind=4), dimension(:),     allocatable   :: gri_Px, gri_Py, gri_Pz

  integer :: nx,ny,nz,nbvec
  real(kind=4) :: x1,y1,z1,x2,y2,z2,dx,dy,dz
  
  end module dim_grid

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
  
  subroutine read_3D_grid(ifc,grid_name)
  
  use dim_grid
  
  character(len=*) :: grid_name
  
  x1=-20.
  x2= 20.
  y1=-20.
  y2= 20.
  z1=-20.
  z2= 20.

  iosta=0
  n=0
  
  print*, '--------------------------------'
  print*, 'reading ',trim(grid_name)
  close(ifc)
  open(ifc,file=grid_name)
  read(ifc,*) nx,ny,nz, nbvec
  print*, 'nx,ny,nz,nbvec=',nx,ny,nz,nbvec,' ie ntot=',nx*ny*nz
  
  if(nbvec /= 1 .and. nbvec /= 4 .and. nbvec /=7) stop 'nbvec only 4 or 7'
  
  dx=(x2-x1)/float(nx)
  dy=(y2-y1)/float(ny)
  dz=(z2-z1)/float(nz)
  
  allocate(nval(nx,ny,nz))
  nval(:,:,:)=0
  
! B or J
  if(nbvec == 4 .or. nbvec == 7) then 
     allocate(grid_Vx(nx,ny,nz))
     allocate(grid_Vy(nx,ny,nz))
     allocate(grid_Vz(nx,ny,nz))
     allocate(grid_V0(nx,ny,nz))

     grid_Vx(:,:,:)=0.
     grid_Vy(:,:,:)=0.
     grid_Vz(:,:,:)=0.
     grid_V0(:,:,:)=0.
  endif
  
! J only 

  if(nbvec == 7) then
     allocate(gri_div(nx,ny,nz))
     allocate(gri_dsc(nx,ny,nz))
     allocate(gri_abj(nx,ny,nz))
     
     gri_div(:,:,:)=0.
     gri_dsc(:,:,:)=0.
     gri_abj(:,:,:)=0.
  endif

  ntot=0
  nvt =0
  ipcp=0
  iosta2=0
  
  write(*,*)
  write(*,'("reading...  ",i9," points")') nx*ny*nz

  do i=1,nx
  ipc=int(25.*float(i)/float(nx))
  if(ipc /= ipcp) write(*,'("reading...  ",i3," %")') 4*ipc
  ipcp=ipc
  
  do j=1,ny
     do k=1,nz
        if(nbvec == 4) read(ifc,*,iostat=iosta2,end=10) grid_Vx(i,j,k),grid_Vy(i,j,k),grid_Vz(i,j,k), grid_V0(i,j,k),nval(i,j,k) 
        if(nbvec == 7) read(ifc,*,iostat=iosta2,end=10) grid_Vx(i,j,k),grid_Vy(i,j,k),grid_Vz(i,j,k), grid_V0(i,j,k), &
                                          gri_div(i,j,k),gri_dsc(i,j,k),gri_abj(i,j,k), nval(i,j,k)                 
        if(nbvec == 1) read(ifc,*,iostat=iosta2,end=10) grid_V0(i,j,k),nval(i,j,k) 
        
        if(iosta2 /= 0) then
                 print*,'iosta2=',iosta2
                 print*,'record #',ntot,' cannot be read'
                 go to 10
        endif
         
         ntot=ntot+1
        nvt =nvt +nval(i,j,k)
     enddo
  enddo
  enddo

 10 continue
 
  ipc=100
  write(*,'("read !!!     ",i3," %")') ipc
  write(*,*)
  
  close(ifc)


  print*, trim(grid_name),' read, nvt,ntot=',nvt,ntot
  print*, '-----------------------------------------------------'
  
  return
  end 

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

  subroutine alloc_grid1D()
  
  use dim_grid

  allocate(gri_Px(nx))
  allocate(gri_Py(ny))
  allocate(gri_Pz(nz))

  dx=(x2-x1)/float(nx)
  dy=(y2-y1)/float(ny)
  dz=(z2-z1)/float(nz)

  do i=1,nx
     gri_Px(i)=float(i-1)*dx +x1 +0.5*dx
  enddo

  do j=1,ny
     gri_Py(j)=float(j-1)*dy +y1 +0.5*dy
  enddo

  do k=1,nz
     gri_Pz(k)=float(k-1)*dz +z1 +0.5*dz 
  enddo
!-----------------------------------------------------

  return
  end 
  
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX 

  subroutine give_model_B(X,Y,Z,VX,VY,VZ,iout)
      
  use dim_grid  
  
  common /Rweigh/ Rint
  
  ! print*,dx,dy,dz
  
  ix=int((x-x1)/dx) +1
  jy=int((y-y1)/dy) +1
  kz=int((z-z1)/dz) +1
  
  if(ix < 1 .or. jy < 1 .or. kz < 1) then
     iout=1
     return
  endif
 
  if(ix > nx .or. jy > ny .or. kz > nz) then
     iout=1
     return
  endif  
  
  
! limites a l'interieur de Rint  
    
  ixr1=int((x-Rint-x1)/dx) +1
  jyr1=int((y-Rint-y1)/dy) +1
  kzr1=int((z-Rint-z1)/dz) +1
    
  ixr2=int((x+Rint-x1)/dx) +1
  jyr2=int((y+Rint-y1)/dy) +1
  kzr2=int((z+Rint-z1)/dz) +1
  
  ixr1=max(min(ixr1,nx),1)
  jyr1=max(min(jyr1,ny),1)
  kzr1=max(min(kzr1,nz),1)
  
  ixr2=max(min(ixr2,nx),1)
  jyr2=max(min(jyr2,ny),1)
  kzr2=max(min(kzr2,nz),1)

  
! --------------------
!    Pondération de Franke-Little
!    R. Franke, « Scattered data Interpolation: Tests of some methods », 
!    Mathematics of Computation, vol. 38,‎ 1982, p. 181–200
  
! ---------------------------------------  
! X
! --------------------------------------- 

  dxb=2.*Rint
  dyb=2.*Rint
  dzb=2.*Rint
  
  somvx=0.
  somvy=0.
  somvz=0.
  
  somwx=0.
  somwy=0.
  somwz=0.
  
  do i=ixr1,ixr2
     do j=jyr1,jyr2
       do k=kzr1,kzr2
          if(nval(i,j,k) /= 0) then
          
           dxb=abs(gri_Px(i)-x)
           dyb=abs(gri_Py(j)-y)
           dzb=abs(gri_Pz(k)-z)
           
           wx=max(1. -dxb/Rint,0.)
           wy=max(1. -dyb/Rint,0.)
           wz=max(1. -dzb/Rint,0.)
           
           somvx=somvx +grid_Vx(i,j,k)*wx
           somvy=somvy +grid_Vy(i,j,k)*wy
           somvz=somvz +grid_Vz(i,j,k)*wz
           
           somwx=somwx +wx
           somwy=somwy +wy
           somwz=somwz +wz
          endif
       enddo
     enddo
  enddo
  
  if(somwx > 0.) then
     vx=somvx/somwx
     ioutx=0
     else
     vx=0.
     ioutx=1
  endif
  
  if(somwy > 0.) then
     vy=somvy/somwy
     iouty=0
     else
     vy=0.
     iouty=1
  endif
  
  if(somwz > 0.) then
     vz=somvz/somwz
     ioutz=0
     else
     vz=0.
     ioutz=1
  endif  
  
  if(ioutx ==1 .or. iouty == 1 .or. ioutz == 1) then
     iout= 1
     else
     iout=0
   endif     
 

  return
  end
  
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
