  program c_curl_from_B3D_grid
  
! -----------------------------------------------------------------------------
! 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
!
! Program : c_curl_from_B3D_grid.f90
! Object  : Compute J 3D grid from a B 3D grid
! Input   : path of a B 3D grid file
! Output  : J 3D grid file
! Author  : P. Robert, LPP-ScientiDev, 2020-Oct.2021
! Mail    : patrick.robert@lpp.polytechnique.fr
! -----------------------------------------------------------------------------

  use dim_grid

  integer     , dimension(:,:,:), allocatable   :: nvalJ
  
  real satmag(3,4), satpos(3,4)

  character(len=255) :: grid_path,grid_name,datafileJ
  character(len=1)  :: BJ,pm
  character(len=3)  :: cdta,ckp
  
  common /vecjfa/  vjx134 , vjx142 , vjx123 , vjx432 , &
                   vjy134 , vjy142 , vjy123 , vjy432 , &
                   vjz134 , vjz142 , vjz123 , vjz432 

  print*,'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
  print*
  print*, 'grid 3D  path? ex: grid_B_mIGRF/grid_Bxyz_160160160_p00_AKp.dat'
  read (*,'(a)') grid_path
  print*, grid_path

  i=index(grid_path,'/',back=.true.)
  if(i>0) then
      grid_name=grid_path(i+1:)
          else
      grid_name=grid_path
  endif

  print*, '------------------------------------------------------------------'
  print*,'grid_name= ',trim(grid_name)
  print*, '------------------------------------------------------------------'

  read(grid_name,'(5x,a1,14x,a1,i2,1x,a3)')BJ,pm,idta,ckp
  dta=float(idta)
  if(pm == 'm') dta=-dta
  cdta=grid_name(21:23)
  
! ---------------
! read grid data
! ---------------
  
  call read_3D_grid(1,grid_path)
  
  print*, 'x1,x2,y1,y2,z1,z2=',x1,x2,y1,y2,z1,z2
  print*, 'nx,ny,nz,dx,dy,dz=',nx,ny,nz,dx,dy,dz
  
! ============================================================== 
! calcul de J
! ============================================================== 

  allocate(nvalJ(nx,ny,nz))
  
  nvalj(:,:,:)=0
  
  write(datafileJ, '(a,3i3.3,5a)') 'grid_Jxyz_',nx,ny,nz,'_',cdta,'_',ckp,'.dat'
  
  print*
  print*, 'creation du fichier ',trim(datafileJ)
  print*, '------------------------------------------'
  print*
  
  open(13,file=datafileJ)

  write(13,'(i8,1x,i8,1x,i8,1x,i1)') nx,ny,nz,7 

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

  ipcp=0

ipcp=0
            
  do ix=1,nx
       ipc=int(50.*float(ix)/float(nx))
       if(ipc /= ipcp) write(*,'("writing... ",i3," %")') 2*ipc
       ipcp=ipc
       
       Pxi=float(ix-1)*dx +x1 +0.5*dx
  do iy=1,ny
       Pyi=float(iy-1)*dy +y1 +0.5*dy
  do iz=1,nz
       Pzi=float(iz-1)*dz +z1 +0.5*dz
     
     
     if(ix == nx .or. iy == ny .or. iz == nz) then
                nvalj(ix,iy,iz)=0
                go to 10
     endif
       
     if(nval(ix  ,iy  ,iz  ) /= 0 .and. &
        nval(ix+1,iy  ,iz  ) /= 0 .and. &
        nval(ix  ,iy+1,iz  ) /= 0 .and. &
        nval(ix  ,iy  ,iz+1) /= 0 ) then
        
        nvalj(ix,iy,iz)=1 
        
        satpos(1,1)=Pxi
        satpos(2,1)=Pyi
        satpos(3,1)=Pzi
        
        satpos(1,2)=Pxi +dx
        satpos(2,2)=Pyi
        satpos(3,2)=Pzi
        
        satpos(1,3)=Pxi
        satpos(2,3)=Pyi +dy
        satpos(3,3)=Pzi
        
        satpos(1,4)=Pxi
        satpos(2,4)=Pyi
        satpos(3,4)=Pzi +dz
        
        satmag(1,1)=grid_Vx(ix,iy,iz)
        satmag(2,1)=grid_Vy(ix,iy,iz)
        satmag(3,1)=grid_Vz(ix,iy,iz)
        
        satmag(1,2)=grid_Vx(ix+1,iy,iz)
        satmag(2,2)=grid_Vy(ix+1,iy,iz)
        satmag(3,2)=grid_Vz(ix+1,iy,iz)
        
        satmag(1,3)=grid_Vx(ix,iy+1,iz)
        satmag(2,3)=grid_Vy(ix,iy+1,iz)
        satmag(3,3)=grid_Vz(ix,iy+1,iz)
        
        satmag(1,4)=grid_Vx(ix,iy,iz+1)
        satmag(2,4)=grid_Vy(ix,iy,iz+1)
        satmag(3,4)=grid_Vz(ix,iy,iz+1)

        satpos(:,:)=satpos(:,:)*6372.2
        
        call cocurldiv(satpos,satmag,vjx,vjy,vjz,curl,div,ajb,elong,plana)
        
        vj0=sqrt(vjx**2 +vjy**2 +vjz**2)
      
        dsc=99.9999
        if(abs(curl) > 1.e-20) dsc=div/curl
            
             else

        nvalj(ix,iy,iz)=0
        vjx=1.e30
        vjy=1.e30
        vjz=1.e30
        vj0=1.e30
        div=1.e30
        ajb=1.e30
        dsc=1.e30
        
     endif        
 10  continue       
     write(13,'(7E13.5,i8)') vjx,vjy,vjz,vj0,abs(div),abs(dsc),ajb, nvalj(ix,iy,iz)

  enddo
  enddo
  enddo

  close(13)

! ============================================================== 

  print*
  print*, "    c_curl_from_B3D_grid.exe  : NORMAL TERMINATION"
         stop "c_curl_from_B3D_grid.exe  : NORMAL TERMINATION"

  end 

!     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 

