  program plot_3D_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 : plot_3D_grid.f90
! Object  : plot slice of 3D grid in XZ, XY and YZ plane
! Input   : path of a 3D grid file
!           origine x for YZ plane, y for XZ plane and z for XY plane (RE)
!           optional plot of low values (1=plot, 0=noplot)
!           width of the slice for the 3 planes (RE)
! Output  : PostScript file 
! Author  : P. Robert, LPP-ScientiDev, 2020-Oct.2021
! Mail    : patrick.robert@lpp.polytechnique.fr
! -----------------------------------------------------------------------------
  use dim_grid ! to have access to grid_Vx(nx,ny,nz) etc.

  character(len=255) :: grid_path,grid_name,psfile
  character(len=1)   :: BJ,pm,ckp
  character(len=10)  :: coxg,coyg,cozg
  character(len= 8)  :: ctolx,ctoly,ctolz
  
!  real(kind=4), dimension(:,:,:), allocatable   :: grid_nbp

  print*,'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
  print*
  print*, 'grid 3D  path? ex: grids_results/grid_Bxyz_080080080_p001.dat'
  read (*,'(a)') grid_path
  print*, grid_path
  
  print*, 'ox,oy,oz for cut plane ? (ex: 0. 0. 0.)'
  read *,  oxg,oyg,ozg
  print*,  oxg,oyg,ozg
  
  print*, 'tolerence x y z ? (ex: 2 RE)'
  read *, tolx,toly,tolz
  print*, tolx,toly,tolz
  
  print*,'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'

  write(coxg,'("x=",f5.1)') oxg  
  write(coyg,'("y=",f5.1)') oyg
  write(cozg,'("z=",f5.1)') ozg
  
  write(ctolx,'(" tol=",f3.1)') tolx
  write(ctoly,'(" tol=",f3.1)') toly
  write(ctolz,'(" tol=",f3.1)') tolz
  
  i=index(grid_path,'/',back=.true.)
  if(i>0) then
      grid_name=grid_path(i+1:)
          else
      grid_name=grid_path
  endif
  
  print*, 'grid_name=',trim(grid_name)
  read(grid_name,'(5x,a1,14x,a1,i2,1x,a3)') BJ,pm,idta,ckp
  dta=float(idta)
  if(pm == 'm') dta=-dta

  print*, '------------------------------------------------------------------'
  print*,'grid_name= ',trim(grid_name),'  dta=',dta,' kp=',ckp
  print*, '------------------------------------------------------------------'


! ---------------
! read grid data
! ---------------
  
  call read_3D_grid(1,grid_path)

! ---------------
! plot grid data
! ---------------
  
  print*,'dx,dy,dz=',dx,dy,dz

  itx=int(tolx/dx +0.5)
  ity=int(toly/dy +0.5)
  itz=int(tolz/dz +0.5)
  
  print*, 'itx,ity,itz=',itx,ity,itz

! ==============================================================   
!  Module pour le vecteur de la grille 
! ==============================================================   

  print*, 'plot du module de '//BJ//' dans le plan xz at oy=',oyg 
  i=index(grid_name,'.')
  psfile=grid_name(1:i-1)//'_XZ_mod.ps'
  call plot_mod_XZ(psfile,grid_V0,coyg//trim(grid_name)//ctoly,BJ,oyg,ity)
  
  print*, 'plot du module de '//BJ//' dans le plan xy at ozg=',ozg
  psfile=grid_name(1:i-1)//'_XY_mod.ps'
  call plot_mod_XY(psfile,grid_V0,cozg//trim(grid_name)//ctolz,BJ,ozg,itz)

  print*, 'plot du module de '//BJ//' dans le plan yz at oxg=',oxg
  psfile=grid_name(1:i-1)//'_YZ_mod.ps'
  call plot_mod_YZ(psfile,grid_V0,cozg//trim(grid_name)//ctolx,BJ,oxg,itx)

! ============================================================== 
! Nombre de points de la grille
! ============================================================== 
  
  ! allocate(grid_nbp(nx,ny,nz))
  
  ! grid_nbp(:,:,:)=float(nval(:,:,:))
  
  ! print*, 'plot du  nombre de points dans le plan xz at oy=',oyg 
  ! i=index(grid_name,'.')
  ! psfile='grid_Npnt'//grid_name(10:i-1)//'_XZ_Nbp.ps'
  ! call plot_mod_XZ(psfile,grid_V0,coyg//trim(grid_name)//ctoly,BJ,oyg,ity)
  
  ! print*, 'plot du  nombre de points dans le plan xy at oz=',ozg
  ! psfile='grid_Npnt'//grid_name(10:i-1)//'_XY_Nbp.ps'
  ! call plot_mod_XY(psfile,grid_V0,cozg//trim(grid_name)//ctolz,BJ,ozg,itz)
  
  ! print*, 'plot du  nombre de points dans le plan yz at ox=',oxg
  ! psfile='grid_Npnt'//grid_name(10:i-1)//'_YZ_Nbp.ps'
  ! call plot_mod_YZ(psfile,grid_V0,coxg//trim(grid_name)//ctolx,BJ,oxg,itx)
  
! ============================================================== 
! Si J plot en plus de div/curl et (J,B) angle
! ============================================================== 

  if(BJ /= 'J') go to 10
  if(BJ /= 'A') go to 10 ! provisoire pour eviter ces plots
  print*, 'nbvec=',nbvec
  if(nbvec /= 7) stop 'J but nbvec /= 7, Dsc et AnJB not plotted'
  
  print*, 'plot du module de div/curl dans le plan xz at oyg=',oyg
  psfile='grid_DisC'//grid_name(10:i-1)//'_XZ_mod.ps'
  call plot_mod_XZ(psfile,gri_dsc,coyg//trim(grid_name)//ctoly,'R',oyg,ity)
  
  print*, 'plot du module de div/curl dans le plan xy at ozg=',ozg
  psfile='grid_DisC'//grid_name(10:i-1)//'_XY_mod.ps'
  call plot_mod_XY(psfile,gri_dsc,cozg//trim(grid_name)//ctolz,'R',ozg,itz)

  print*, 'plot du module de div/curl dans le plan yz at oxg=',oxg
  psfile='grid_DisC'//grid_name(10:i-1)//'_YZ_mod.ps'
  call plot_mod_YZ(psfile,gri_dsc,coxg//trim(grid_name)//ctolx,'R',oxg,itx)
  
  
  
  print*, 'plot du module de (J,B) dans le plan xz at oyg=',oyg  
  psfile='grid_AnJB'//grid_name(10:i-1)//'_XZ_mod.ps'
  call plot_mod_XZ(psfile,gri_abj,coyg//trim(grid_name)//ctoly,'A',oyg,ity)
  
  print*, 'plot du module de (J,B) dans le plan xy at ozg=',ozg
  psfile='grid_AnJB'//grid_name(10:i-1)//'_XY_mod.ps'
  call plot_mod_XY(psfile,gri_abj,cozg//trim(grid_name)//ctolz,'A',ozg,itz)

  print*, 'plot du module de (J,B) dans le plan yz at oxg=',oxg
  psfile='grid_AnJB'//grid_name(10:i-1)//'_YZ_mod.ps'
  call plot_mod_YZ(psfile,gri_abj,coxg//trim(grid_name)//ctolx,'A',oxg,itx)
  
! ============================================================== 
10 continue
  print*
  print*, "    plot_3D_grid.exe  : NORMAL TERMINATION"
         stop "plot_3D_grid.exe  : NORMAL TERMINATION"

  end 

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
  subroutine plot_mod_XZ(psfile,grid,grid_name,BJ,oyg,ity)
  
  use dim_grid
  
  real(kind=4), dimension(nx,ny,nz) :: grid
  character(len=10) ::  colmap
  character(len=1)  :: BJ
  character(len=*)  :: grid_name,psfile
  character(len=80) :: label
  
  if(allocated(grid_2D)) deallocate(grid_2D)
  if(allocated(nval_2D)) deallocate(nval_2D)
  if(allocated(imagexz)) deallocate(imagexz)
  
  allocate(grid_2D(nx,nz))
  allocate(nval_2D(nx,nz))
  allocate(imagexz(nx,nz))

  
  grid_2D(:,:)=0.
  nval_2D(:,:)=0
  
  orix=1.5
  oriy=4.5
  sizx=18.
  sizy=18.

  tc=0.4 
  
  call dopegra_(7,psfile)
  call dfigsiz_(sizx,sizy)
  call dfigori_(orix,oriy)
  call dfiglim_(x1,x2,z1,z2)
  call dcolmap_('spectro064')
  call gcolmap_(colmap,nlevel)

! moyenne autour de Y=oyg

  noyg=max(min(int((oyg-y1)/dy),ny),1) 

  iy1=max(min(noyg -ity,ny),1)
  iy2=max(min(noyg +ity,ny),1)
  
  print*,'iy1,iy2=',iy1,iy2
  
  do ix=1,nx
       Pxi=float(ix-1)*dx +x1 +0.5*dx
  do iz=1,nz
       Pzi=float(iz-1)*dz +z1 +0.5*dz
    ! on reduit l'intervalle d'integration pres de la Terre
     r=sqrt(Pxi**2 +Pzi**2)
     idy=ity
     if(r < 3.) idy=-noyg
    iy1=max(min(noyg -idy,ny),1)
    iy2=max(min(noyg +idy,ny),1)
  do iy=iy1,iy2
     if(nval(ix,iy,iz) /= 0 .and. abs(grid(ix,iy,iz)) < 1.e30) then
        grid_2D(ix,iz)= grid_2D(ix,iz) +grid(ix,iy,iz)
        nval_2D(ix,iz)= nval_2D(ix,iz) +1
     endif
  enddo
  enddo
  enddo

! normalisation

  do ix=1,nx
  do iz=1,nz
  if(nval_2D(ix,iz) /= 0) then
                grid_2D(ix,iz)=grid_2D(ix,iz)/float(nval_2D(ix,iz))
                          else
                grid_2D(ix,iz)=-1.
  endif
  enddo
  enddo
     
  print*, 'legende'
  
  if(BJ == 'B') label='B0 log(nT)'
  if(BJ == 'J') label='J0 log(nA/m2)'
  if(BJ == 'D') label='|Div(B)| log(nT/m)'
  if(BJ == 'R') label='|Div/curl|'
  if(BJ == 'A') label='(J,B) angle degrees)'
  if(BJ == 'N') label='counts log'
  
  ilog=1
  if(BJ == 'A') ilog=0
  if(BJ == 'R') ilog=0
  
  print*, 'BJ=',BJ,'   ilog=',ilog,'   label=',trim(label)
  
  call ppagcha_(orix,oriy+sizy+3.*tc,-1,tc,tc,0.,trim(label))
  call ppagcha_(orix,oriy+sizy+tc,-1,tc,tc,0.,trim(grid_name))
  call plot_grid_mod(grid_2D,nx,nz,imagexz,ilog,'XZ')
 
  call dclogra_
  
  return
  end

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
  subroutine plot_mod_YZ(psfile,grid,grid_name,BJ,oxg,itx)
  
  use dim_grid
  
  real(kind=4), dimension(nx,ny,nz) :: grid
  character(len=10) ::  colmap
  character(len=1)  :: BJ
  character(len=*)  :: grid_name,psfile
  character(len=80) :: label
  
  if(allocated(grid_2D)) deallocate(grid_2D)
  if(allocated(nval_2D)) deallocate(nval_2D)
  if(allocated(imageyz)) deallocate(imageyz)
  
  allocate(grid_2D(ny,nz))
  allocate(nval_2D(ny,nz))
  allocate(imageyz(ny,nz))

  
  grid_2D(:,:)=0.
  nval_2D(:,:)=0
  
  orix=1.
  oriy=4.5
  sizx=18.
  sizy=18.

  tc=0.4 
  
  call dopegra_(7,psfile)
  call dfigsiz_(sizx,sizy)
  call dfigori_(orix,oriy)
  call dfiglim_(y1,y2,z1,z2)
  call dcolmap_('spectro064')
  call gcolmap_(colmap,nlevel)

! moyenne autour de X=oxg

  noxg=max(min(int((oxg-x1)/dx),nx),1) 

  ix1=max(min(noxg -itx,nx),1)
  ix2=max(min(noxg +itx,nx),1)
  
  print*,'ix1,ix2=',ix1,ix2
  
  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
    ! on reduit l'intervalle d'integration pres de la Terre
     r=sqrt(Pyi**2 +Pzi**2)
     idx=itx
     if(r < 3.) idx=-noxg
    ix1=max(min(noxg -idx,nx),1)
    ix2=max(min(noxg +idx,nx),1)
  do ix=ix1,ix2
     if(nval(ix,iy,iz) /= 0 .and. abs(grid(ix,iy,iz)) < 1.e30) then
        grid_2D(iy,iz)= grid_2D(iy,iz) +grid(ix,iy,iz)
        nval_2D(iy,iz)= nval_2D(iy,iz) +1
     endif
  enddo
  enddo
  enddo

! normalisation

  do iy=1,ny
  do iz=1,nz
  if(nval_2D(iy,iz) /= 0) then
                grid_2D(iy,iz)=grid_2D(iy,iz)/float(nval_2D(iy,iz))
                          else
                grid_2D(iy,iz)=-1.
  endif
  enddo
  enddo
     
  print*, 'legende'
  
  if(BJ == 'B') label='B0 log(nT)'
  if(BJ == 'J') label='J0 log(nA/m2)'
  if(BJ == 'D') label='|Div(B)| log(nT/m)'
  if(BJ == 'R') label='|Div/curl|'
  if(BJ == 'A') label='(J,B) angle degrees)'
  if(BJ == 'N') label='counts log'
  
  ilog=1
  if(BJ == 'A') ilog=0
  if(BJ == 'R') ilog=0
  
  print*, 'BJ=',BJ,'   ilog=',ilog,'   label=',trim(label)
  
  call ppagcha_(orix,oriy+sizy+3.*tc,-1,tc,tc,0.,trim(label))
  call ppagcha_(orix,oriy+sizy+tc,-1,tc,tc,0.,trim(grid_name))
  call plot_grid_mod(grid_2D,ny,nz,imageyz,ilog,'YZ')
 
  call dclogra_
  
  return
  end

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
  subroutine plot_mod_XY(psfile, grid,grid_name,BJ,ozg,itz)
  
  use dim_grid
  
  real(kind=4), dimension(nx,ny,nz) :: grid
  character(len=10) ::  colmap
  character(len=1)  :: BJ
  character(len=*)  :: grid_name,psfile
  character(len=80) :: label

  if(allocated(grid_2D)) deallocate(grid_2D)
  if(allocated(nval_2D)) deallocate(nval_2D)
  if(allocated(imagexy)) deallocate(imagexy)
  
  allocate(grid_2D(nx,ny))
  allocate(nval_2D(nx,ny))
  allocate(imagexy(nx,ny))
  
  grid_2D(:,:)=0.
  nval_2D(:,:)=0
  
  orix=1.5
  oriy=4.5
  sizx=18.
  sizy=18.

  tc=0.4
  
  call dopegra_(7,psfile)
  call dfigsiz_(sizx,sizy)
  call dfigori_(orix,oriy)
  call dfiglim_(x1,x2,y1,y2)
  call dcolmap_('spectro064')
  call gcolmap_(colmap,nlevel) 

! moyenne autour de Z=ozg

  nozg=max(min(int((ozg-z1)/dz),nz),1) 

  iz1=max(min(nozg -itz,nz),1)
  iz2=max(min(nozg +itz,nz),1)

  do ix=1,nx
     Pxi=float(ix-1)*dx +x1 +0.5*dx
  do iy=1,ny
     Pyi=float(iy-1)*dy +y1 +0.5*dy
    ! on reduit l'intervalle d'integration pres de la Terre
     r=sqrt(Pxi**2 +Pyi**2)
     idz=itz
     if(r < 3.) idz=-nozg
    iz1=max(min(nozg -idz,nz),1)
    iz2=max(min(nozg +idz,nz),1)

  do iz=iz1,iz2
     if(nval(ix,iy,iz) /= 0 .and. abs(grid(ix,iy,iz)) < 1.e30) then
        grid_2D(ix,iy)= grid_2D(ix,iy) +grid(ix,iy,iz)
        nval_2D(ix,iy)= nval_2D(ix,iy) +1
     endif
  enddo
  enddo
  enddo

! normalisation

  do ix=1,nx
  do iy=1,ny
  if(nval_2D(ix,iy) /= 0) then
                grid_2D(ix,iy)=grid_2D(ix,iy)/float(nval_2D(ix,iy))
                          else
                grid_2D(ix,iy)=-1.
  endif
  enddo
  enddo
  
  print*, 'legende'
  
  if(BJ == 'B') label='B0 log(nT)'
  if(BJ == 'J') label='J0 log(nA/m2)'
  if(BJ == 'D') label='|Div(B)| log(nT/m)'
  if(BJ == 'R') label='|Div/curl|'
  if(BJ == 'A') label='(J,B) angle degrees)'
  if(BJ == 'N') label='counts log'
  
  ilog=1
  if(BJ == 'A') ilog=0
  if(BJ == 'R') ilog=0
  
  print*, 'BJ=',BJ,'   ilog=',ilog,'   label=',trim(label)
  
  call ppagcha_(orix,oriy+sizy+3.*tc,-1,tc,tc,0.,trim(label))
  call ppagcha_(orix,oriy+sizy+tc,-1,tc,tc,0.,trim(grid_name))
  call plot_grid_mod(grid_2D,nx,ny,imagexy,ilog,'XY')

  call dclogra_
  
  return
  end

! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
 
  subroutine plot_grid_mod(grid,iresx,iresy,image,ilog,XX)

  real(kind=4) :: grid(iresx,iresy)
  integer :: image(iresx,iresy)!,nvalgr(iresx,iresy)
  character(len=2) :: XX
  real (kind=4), dimension(1072) :: bowshox,bowshoz
  real(kind=4) :: romagx(201),romagy(201)
  character(len=10) ::  colmap

  call gfigsiz_(sizx,sizy)
  call gfigori_(orix,oriy)
  call gfiglimx(x1,x2)
  call gfiglimy(y1,y2)
  call dcolmap_('spectro064')
  call gcolmap_(colmap,nlevel)

! passage en log (deviens -8 si egal zero et -10 si negatif)

  if(ilog == 1 ) then
       print*, 'data in log scale'
       valzer=1.e-8
       valneg=1.e-10
       call cscalog_(grid,iresx,iresy,valzer,valneg)

       pcmin=10.
       pcmax=1.
                 else
       pcmin=0.
       pcmax=0.
  endif

  print*, 'calcul des min et max de la grid'
  
  smin=minval(grid)
  smax=maxval(grid)

!  call  cscasli_(grid,iresx,iresy,smin,smax,pcmin,pcmax)

  print*
  print*, 'Image: valmin,valmax=',smin,smax
  print*
  print*,'imposed values?'
  read(*,*,end=99)  vmin,vmax
  print*, vmin,vmax

  if((vmax-vmin) > 0.) then
           smin=vmin
           smax=vmax
  endif

  print*, 'valeurs prises:',smin,smax

  swhi= smax
  sbla= smin
  
  print*, 'calcul de l''image avec smin, smax=',smin,smax
  call cscaima_(grid,iresx,iresy,smin,smax,image,sbla,swhi)

! on met les valeurs faibles en blanc a la place de noir
 
  call crgbico_(255,255,255,iblan)
  call crgbico_(  0,  0,  0,inoir)
  
  do i=1,iresx
  do j=1,iresy
     if(image(i,j) == inoir) image(i,j)=iblan
  enddo
  enddo
  
  print*, 'plot de l''image'
  call ppagima_(orix,oriy,-1,sizx,sizy,0.,image,iresx,iresy)

  if(XX == 'YZ') go to 10

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

  print*,'bow shock'
  call cobowsho2(bowshox,bowshoz,nbow)
  call dlincol_('m')
  call dlinwid_(5)
  call pfigcur_(bowshox,bowshoz,nbow)

  print*,'magnetopause PR'
  call comagpau_pr(romagx,romagy,nmagpau)
  call dlincol_('b')
  call dlinwid_(7)
  call pfigcur_(romagx,romagy,201)
  call dlincol_('n')
!-----------------------------------------------------
  
10 continue

! graduations

  call dlincol_('n')
  call dlinwid_(3)
  call pfigfra_
  call pfiggrax(0.,5.,1.,'(f4.0)')
  call pfiggray(0.,5.,1.,'(f4.0)')

  call dlabsiz_(0.4,0.4)
  call pfiglabx(xx(1:1)//' GSM')
  call pfiglaby(xx(2:2)//' GSM')

! plot de la Terre

  call dlincol_('c')
  call dfilzon_
  call pfigcir_(0.,0.,1.,90.,270.,5.)
  call pfilzon_

  call dlincol_('y')
  call dfilzon_
  call pfigcir_(0.,0.,1.,-90.,90.,5.)
  call pfilzon_

  call dlincol_('n')
  call pfigcir_(0.,0.,1.,0.,360.,5.)
  call pfiglin_(0.,-1.,0.,1.)
  call dlincol_('n')
  call pfigfra_

! plot de la mire (redefini une figure)

  sx=sizx/3.
  sy=0.8
  ox=orix+sizx-sx
  oy=oriy+sizy+1.
  tgra=0.30
  call plotpalhsbx(ox,oy,sx,sy,vmin,vmax,tgra)

  print*, 'nb level=',nlevel
  
  99 continue

  return
  end
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!
  subroutine plotpalhsbx(ox,oy,sx,sy,xmin,xmax,tgra)
!
! --------------------------------------------------------------------+--
! plot d'une palette HSB
! --------------------------------------------------------------------+--
!
  character*8 forx
!
! *** palette hsb avec cadre
!
  call  ppagcmah(ox,oy,sx,sy,1)
  call  ppagrec_(ox,oy,sx,sy)
!
! *** graduation de l'axe avec les bornes exactes
!
! *   calcul des intervalles des graduations
!
  call  cfiggra_(xmin,xmax,x1a,x2a,bgx,sgx,forx)
!
! *   definition de la figure incluse correspondant aux grad. rondes
!
  call  dfigori_(ox,oy)
  call  dfigsiz_(sx,sy)
  call  dfiglimx(xmin,xmax)
!
! *   trace des graduations
!
  call  dstipos_(' o','ii')
  call  dgrapos_('t','l')
  call  dstisiz_(sy/5.,sy/10.)
  call  dgrasiz_(tgra,tgra)
  call  pfiggrax(x1a+bgx,bgx,sgx,forx)
!
  return
  END
!
! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
!


