! This file is part of GTC version 4.4 
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without 
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice, 
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice, 
!    this list of conditions and the following disclaimer in the documentation 
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be 
!    used to endorse or promote products derived from this software without 
!    specific prior written permission.

subroutine eqplot
  use global_parameters
  use field_array
  use equilibrium
  use spline_function, only: sptem, tem_dp, spden, den_dp, spzeff, sptrot, sper, spq, &
       dq_dp, spgpsi, sppressure, sprpsi, sptorpsi, sprgpsi, sppsitor, sppsirg, splcos, &
       splsin,spjac, currenti, delb, spx, spz, spb, spf
#ifdef _SDP
  use sdp, only: sdp_equilibrium_output, sdp_diag
#endif
  implicit none
  
  integer,parameter :: mskip=4
  integer i,j,ieq,nplot,nrad
  real(lk) delpsi,dum,pdum,tdum,deltheta,datap(lsp),data1d(lsp),&
    datax(mpsi/mskip+1,lst),dataz(mpsi/mskip+1,lst),data2d(mpsi/mskip+1,lst,5)

#ifdef ADIOS2
  call eqplot_adios2
#endif
! #else

  ieq=120 !open equilibiurm plot output file
  open(ieq,file='equilibrium.out',status='replace')
101 format(i6)
102 format(e16.8)

!# of 1D radial plots
  nplot=29
!  nrad=mpsi+1
  nrad=lsp
  write(ieq,101)nplot,nrad

!0: radial axis using poloidal flux function
!  delpsi=(psimesh(mpsi)-psimesh(0))/real(nrad-1)
  delpsi=psiw/real(nrad-1)
  do i=1,nrad
!     datap(i)=psimesh(0)+delpsi*real(i-1)
     datap(i)=delpsi*real(i-1)
  enddo
  write(ieq,102)datap 
  
!1: sqaure-root of normalized toroidal flux function
  data1d=0.0
  write(ieq,102)data1d 

!2: minor radius
  data1d=0.0
  write(ieq,102)data1d 

!3: major radius
  data1d=0.0
  write(ieq,102)data1d 

#ifdef _FRC
!4: Te
  do i=1,nrad
     data1d(i)=sptem(datap(i),tfepp)/(rho0*rho0)
  enddo
  write(ieq,102)data1d 

!5: -d(ln(Te))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tfepp)/sptem(dum,tfepp)
  enddo
  write(ieq,102)data1d 

!6: ne
  do i=1,nrad
     data1d(i)=spden(datap(i),nfepp)
  enddo
  write(ieq,102)data1d 

!7: -d(ln(ne))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nfepp)/spden(dum,nfepp)
  enddo
  write(ieq,102)data1d 
#else
!4: Te
  do i=1,nrad
     data1d(i)=sptem(datap(i),tepp)/tepp(1,1)
  enddo
  write(ieq,102)data1d 

!5: -d(ln(Te))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tepp)/sptem(dum,tepp)
  enddo
  write(ieq,102)data1d 

!6: ne
  do i=1,nrad
     data1d(i)=spden(datap(i),nepp)/nepp(1,1)
  enddo
  write(ieq,102)data1d 

!7: -d(ln(ne))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nepp)/spden(dum,nepp)
  enddo
  write(ieq,102)data1d 
#endif
!8: ti
  do i=1,nrad
     data1d(i)=sptem(datap(i),tipp)/tepp(1,1)
  enddo
  write(ieq,102)data1d 

!9: -d(ln(ti))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tipp)/sptem(dum,tipp)
  enddo
  write(ieq,102)data1d

!10: ni
  do i=1,nrad
     data1d(i)=spden(datap(i),nipp)/nepp(1,1)
  enddo
  write(ieq,102)data1d 

!11: -d(ln(ni))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nipp)/spden(dum,nipp)
  enddo
  write(ieq,102)data1d 

#ifndef _FRC
!12: tf
  do i=1,nrad
     data1d(i)=sptem(datap(i),tfpp)/tepp(1,1)
  enddo
  write(ieq,102)data1d 

!13: -d(ln(tf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tfpp)/sptem(dum,tfpp)
  enddo
  write(ieq,102)data1d

!14: nf
  do i=1,nrad
     data1d(i)=spden(datap(i),nfpp)/nepp(1,1)
  enddo
  write(ieq,102)data1d 

!15: -d(ln(nf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nfpp)/spden(dum,nfpp)
  enddo
  write(ieq,102)data1d 
#else
!12: tf
  do i=1,nrad
     data1d(i)=sptem(datap(i),tepp)/(rho0*rho0)
  enddo
  write(ieq,102)data1d 

!13: -d(ln(tf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tepp)/sptem(dum,tepp)
  enddo
  write(ieq,102)data1d

!14: nf
  do i=1,nrad
     data1d(i)=spden(datap(i),nepp)
  enddo
  write(ieq,102)data1d 

!15: -d(ln(nf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nepp)/spden(dum,nepp)
  enddo
  write(ieq,102)data1d 
#endif
!16: zeff
  do i=1,nrad
     data1d(i)=spzeff(datap(i))
  enddo
  write(ieq,102)data1d 

!17: toroidal rotation
  do i=1,nrad
     data1d(i)=sptrot(datap(i))
  enddo
  write(ieq,102)data1d 

!18: radial electric field
  do i=1,nrad
     data1d(i)=sper(datap(i))
  enddo
  write(ieq,102)data1d 

!19: q profile
  do i=1,nrad
     data1d(i)=spq(datap(i))
  enddo
  write(ieq,102)data1d 

!20: d(ln(q))/dpsi
  do i=1,nrad
     dum=datap(i)
     data1d(i)=dq_dp(dum)/spq(dum)
  enddo
  write(ieq,102)data1d

!21: gcurrent profile
  do i=1,nrad
     data1d(i)=spgpsi(datap(i))
  enddo
  write(ieq,102)data1d 

!22: pressure profile
  do i=1,nrad
     data1d(i)=sppressure(datap(i))
  enddo
  write(ieq,102)data1d 

!23: minor radius
  do i=1,nrad
     data1d(i)=sprpsi(datap(i))
  enddo
  write(ieq,102)data1d 

!24: toroidal flux
  do i=1,nrad
     data1d(i)=sptorpsi(datap(i))
  enddo
  write(ieq,102)data1d

!25: radial grid: rgpsi
  do i=1,nrad
     data1d(i)=sprgpsi(datap(i))
  enddo
  write(ieq,102)data1d

!26: inverse of spline torpsi: psitor
  tdum=0.0
  pdum=sptorpsi(psiw)
  dum=(pdum-tdum)/real(nrad-1)
  do i=1,nrad
     pdum=tdum+dum*real(i-1)
     data1d(i)=sppsitor(pdum)
  enddo
  write(ieq,102)data1d

!27: inverse of spline rgpsi: psirg
  tdum=0.0
  pdum=sprgpsi(psiw)
  dum=(pdum-tdum)/real(nrad-1)
  do i=1,nrad
     pdum=tdum+dum*real(i-1)
     data1d(i)=sppsirg(pdum)
  enddo
  write(ieq,102)data1d

!28: error of spline cos in [0, pi/2]
  do i=1,nrad
     dum=0.25*spdtheta*real(i-1)
     data1d(i)=splcos(dum)-cos(dum)
  enddo
  write(ieq,102)data1d

!29: error of spline sin in [0, pi/2]
  do i=1,nrad
     dum=0.25*spdtheta*real(i-1)
     data1d(i)=splsin(dum)-sin(dum)
  enddo
  write(ieq,102)data1d

!# of 2D contour plots on poloidal plane
  nplot=5
  write(ieq,101)nplot,mpsi/mskip+1,lst
  deltheta=2.0*pi/real(lst-1)

  do i=1,mpsi/mskip+1
     pdum=psimesh((i-1)*mskip)
     do j=1,lst
        tdum=deltheta*real(j-1)
#ifdef _TOROIDAL3D
        datax(i,j)=spx(pdum,tdum,0.0_lk)
        dataz(i,j)=spz(pdum,tdum,0.0_lk)
        data2d(i,j,1)=spb(pdum,tdum,0.0_lk)
        data2d(i,j,4)=spf(pdum,tdum,0.0_lk)
#else
        datax(i,j)=spx(pdum,tdum)
        dataz(i,j)=spz(pdum,tdum)
        data2d(i,j,1)=spb(pdum,tdum)
        data2d(i,j,4)=spf(pdum,tdum)
#endif
        data2d(i,j,2)=spjac(pdum,tdum)
        !javierhn: currenti in Boozer coordinates is going to be always 
        ! constant in a flux surface. Not really useful here. It should be
        ! deleted. Keep it due the idl script
        data2d(i,j,3)=currenti(pdum,tdum)
        !javierhn: delb is no longer used and it has been set to zero
        data2d(i,j,5)=delb(pdum,tdum)
     enddo
  enddo

!0-1: mesh points on (X,Z)
  write(ieq,102)datax,dataz 

!2: b-field
  write(ieq,102)data2d(:,:,1)

!3: Jacobian
  write(ieq,102)data2d(:,:,2)

!4: icurrent
  write(ieq,102)data2d(:,:,3)

!5: zeta (magnetic coord.) to phi (cylindrical coord.)
  write(ieq,102)data2d(:,:,4)

!6: delb
  write(ieq,102)data2d(:,:,5)

close(ieq)

! End ifdef for ADIOS2
! #endif

#ifdef _SDP
   call sdp_equilibrium_output
   call sdp_diag
#endif

end subroutine eqplot


!-----------------------------------------------------------------------------
#ifdef ADIOS2
subroutine eqplot_adios2
  use global_parameters
  use field_array
  use equilibrium
  use spline_function, only: sptem, tem_dp, spden, den_dp, spzeff, sptrot, sper, spq, &
       dq_dp, spgpsi, sppressure, sprpsi, sptorpsi, sprgpsi, sppsitor, sppsirg, splcos, &
       splsin,spjac, currenti, delb, spx, spz, spb, spf
#ifdef _SDP
  use sdp, only: sdp_equilibrium_output, sdp_diag
#endif
  use adios2
  use adios2_common
  implicit none

  integer,parameter :: mskip=4
  integer i,j,ieq,nplot,nrad, ierr
  real(lk) delpsi,dum,pdum,tdum,deltheta,datap(lsp),data1d(lsp),&
    datax(mpsi/mskip+1,lst),dataz(mpsi/mskip+1,lst),data2d(mpsi/mskip+1,lst,5)

  type(adios2_variable) :: var_psi
  type(adios2_variable) :: var_lsp
  type(adios2_variable) :: var_dpsi
  type(adios2_variable) :: var_Te0
  type(adios2_variable) :: var_ne0
  type(adios2_variable) :: var_dlgne_dpsi_neg
  type(adios2_variable) :: var_dlgTe_dpsi_neg
  type(adios2_variable) :: var_Ti0
  type(adios2_variable) :: var_dlgTi_dpsi_neg
  type(adios2_variable) :: var_ni0
  type(adios2_variable) :: var_dlgni_dpsi_neg
  type(adios2_variable) :: var_Tf0
  type(adios2_variable) :: var_dlgtf_dpsi_neg
  type(adios2_variable) :: var_nf0
  type(adios2_variable) :: var_dlgnf_dpsi_neg
  type(adios2_variable) :: var_Zeff
  type(adios2_variable) :: var_u_tor
  type(adios2_variable) :: var_Er
  type(adios2_variable) :: var_q
  type(adios2_variable) :: var_dlgq_dpsi
  type(adios2_variable) :: var_g
  type(adios2_variable) :: var_P
  type(adios2_variable) :: var_r_1d
  type(adios2_variable) :: var_torpsi
  type(adios2_variable) :: var_r_1dg
  type(adios2_variable) :: var_dtor
  type(adios2_variable) :: var_psi_tor
  type(adios2_variable) :: var_drg
  type(adios2_variable) :: var_psi_rg
  type(adios2_variable) :: var_dtheta
  type(adios2_variable) :: var_cos_err
  type(adios2_variable) :: var_sin_err
  type(adios2_variable) :: var_npsi
  type(adios2_variable) :: var_ntheta
  type(adios2_variable) :: var_R
  type(adios2_variable) :: var_Z
  type(adios2_variable) :: var_B
  type(adios2_variable) :: var_Jacobian
  type(adios2_variable) :: var_I
  type(adios2_variable) :: var_nu
  type(adios2_variable) :: var_delta

  if (.not. adios2_eqplot_io_initialized) then
    call adios2_declare_io (adios2_eqplot_io, adios2_obj, adios2_eqplot_io_id, ierr)

    call adios2_define_variable (var_psi, adios2_eqplot_io, "psi", adios2_type_real_lk, 1, &
      1_8 * shape(datap), 1_8 * (/0/), 1_8 * shape(datap), .false., ierr)
    
    call adios2_define_variable (var_lsp, adios2_eqplot_io, "lsp", adios2_type_integer4, ierr)

    call adios2_define_variable (var_dpsi, adios2_eqplot_io, "dpsi", adios2_type_real_lk, ierr)

    call adios2_define_variable (var_dtor, adios2_eqplot_io, "dtor", adios2_type_real_lk, ierr)

    call adios2_define_variable (var_drg, adios2_eqplot_io, "drg", adios2_type_real_lk, ierr)

    call adios2_define_variable (var_dtheta, adios2_eqplot_io, "dtheta", adios2_type_real_lk, ierr)

    call adios2_define_variable (var_npsi, adios2_eqplot_io, "npsi", adios2_type_integer4, ierr)

    call adios2_define_variable (var_ntheta, adios2_eqplot_io, "ntheta", adios2_type_integer4, ierr)

    call adios2_define_variable (var_Te0, adios2_eqplot_io, "Te0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgTe_dpsi_neg, adios2_eqplot_io, "dlgTe_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_ne0, adios2_eqplot_io, "ne0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgne_dpsi_neg, adios2_eqplot_io, "dlgne_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_Ti0, adios2_eqplot_io, "Ti0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgTi_dpsi_neg, adios2_eqplot_io, "dlgTi_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_ni0, adios2_eqplot_io, "ni0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgni_dpsi_neg, adios2_eqplot_io, "dlgni_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_Tf0, adios2_eqplot_io, "Tf0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgtf_dpsi_neg, adios2_eqplot_io, "dlgtf_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_nf0, adios2_eqplot_io, "nf0", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgnf_dpsi_neg, adios2_eqplot_io, "dlgnf_dpsi_neg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_Zeff, adios2_eqplot_io, "Zeff", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_u_tor, adios2_eqplot_io, "u_tor", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_Er, adios2_eqplot_io, "Er", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_q, adios2_eqplot_io, "q", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_dlgq_dpsi, adios2_eqplot_io, "dlgq_dpsi", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_g, adios2_eqplot_io, "g", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_P, adios2_eqplot_io, "P", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_r_1d, adios2_eqplot_io, "r_1d", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_torpsi, adios2_eqplot_io, "torpsi", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_r_1dg, adios2_eqplot_io, "rg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_psi_tor, adios2_eqplot_io, "psi_tor", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_psi_rg, adios2_eqplot_io, "psi_rg", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_cos_err, adios2_eqplot_io, "cos_err", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_sin_err, adios2_eqplot_io, "sin_err", adios2_type_real_lk, 1, &
      1_8 * shape(data1d), 1_8 * (/0/), 1_8 * shape(data1d), .false., ierr)

    call adios2_define_variable (var_R, adios2_eqplot_io, "R", adios2_type_real_lk, 2, &
      1_8 * shape(datax), 1_8 * (/0,0/), 1_8 * shape(datax), .false., ierr)

    call adios2_define_variable (var_Z, adios2_eqplot_io, "Z", adios2_type_real_lk, 2, &
      1_8 * shape(dataz), 1_8 * (/0,0/), 1_8 * shape(dataz), .false., ierr)

    call adios2_define_variable (var_B, adios2_eqplot_io, "B", adios2_type_real_lk, 2, &
      1_8 * shape(data2d(:,:,1)), 1_8 * (/0,0/), 1_8 * shape(data2d(:,:,1)), .false., ierr)

    call adios2_define_variable (var_Jacobian, adios2_eqplot_io, "Jacobian", adios2_type_real_lk, 2, &
      1_8 * shape(data2d(:,:,2)), 1_8 * (/0,0/), 1_8 * shape(data2d(:,:,2)), .false., ierr)

    call adios2_define_variable (var_I, adios2_eqplot_io, "I", adios2_type_real_lk, 2, &
      1_8 * shape(data2d(:,:,3)), 1_8 * (/0,0/), 1_8 * shape(data2d(:,:,3)), .false., ierr)

    call adios2_define_variable (var_nu, adios2_eqplot_io, "nu", adios2_type_real_lk, 2, &
      1_8 * shape(data2d(:,:,4)), 1_8 * (/0,0/), 1_8 * shape(data2d(:,:,4)), .false., ierr)

    call adios2_define_variable (var_delta, adios2_eqplot_io, "delta", adios2_type_real_lk, 2, &
      1_8 * shape(data2d(:,:,5)), 1_8 * (/0,0/), 1_8 * shape(data2d(:,:,5)), .false., ierr)

    call adios2_open (adios2_eqplot_e, adios2_eqplot_io, "equilibrium.bp", &
      adios2_mode_write, mpi_comm_self, ierr)

    adios2_eqplot_io_initialized = .true.
  endif ! if .not. adios2_eqplot_io_initialized

  call adios2_begin_step(adios2_eqplot_e, adios2_step_mode_append, ierr)

  !# of 1D radial plots
  nplot=29
  nrad=lsp

  !0: radial axis using poloidal flux function
  ! delpsi=(psimesh(mpsi)-psimesh(0))/real(nrad-1)
  delpsi=psiw/real(nrad-1)
  do i=1,nrad
     datap(i)=delpsi*real(i-1)
  enddo

  call adios2_put(adios2_eqplot_e, var_lsp, nrad, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_dpsi, delpsi, adios2_mode_sync, ierr)
  
  call adios2_put(adios2_eqplot_e, var_psi, datap, adios2_mode_sync, ierr)

#ifdef _FRC
  !4: Te
  do i=1,nrad
     data1d(i)=sptem(datap(i),tfepp)/(rho0*rho0)
  enddo
  call adios2_put(adios2_eqplot_e, var_Te0, data1d, adios2_mode_sync, ierr)
  
  !5: -d(ln(Te))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tfepp)/sptem(dum,tfepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgTe_dpsi_neg, data1d, adios2_mode_sync, ierr)
  
  !6: ne
  do i=1,nrad
     data1d(i)=spden(datap(i),nfepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_ne0, data1d, adios2_mode_sync, ierr)
  
  !7: -d(ln(ne))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nfepp)/spden(dum,nfepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgne_dpsi_neg, data1d, adios2_mode_sync, ierr)
#else
  !4: Te
  do i=1,nrad
     data1d(i)=sptem(datap(i),tepp)/tepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_Te0, data1d, adios2_mode_sync, ierr)
  
  !5: -d(ln(Te))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tepp)/sptem(dum,tepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgTe_dpsi_neg, data1d, adios2_mode_sync, ierr)
  
  !6: ne
  do i=1,nrad
     data1d(i)=spden(datap(i),nepp)/nepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_ne0, data1d, adios2_mode_sync, ierr)
  
  !7: -d(ln(ne))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nepp)/spden(dum,nepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgne_dpsi_neg, data1d, adios2_mode_sync, ierr)
#endif
  
  !8: ti
  do i=1,nrad
     data1d(i)=sptem(datap(i),tipp)/tepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_Ti0, data1d, adios2_mode_sync, ierr)

  !9: -d(ln(ti))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tipp)/sptem(dum,tipp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgTi_dpsi_neg, data1d, adios2_mode_sync, ierr)

  !10: ni
  do i=1,nrad
     data1d(i)=spden(datap(i),nipp)/nepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_ni0, data1d, adios2_mode_sync, ierr)

  !11: -d(ln(ni))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nipp)/spden(dum,nipp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgni_dpsi_neg, data1d, adios2_mode_sync, ierr)

#ifndef _FRC
  !12: tf
  do i=1,nrad
     data1d(i)=sptem(datap(i),tfpp)/tepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_Tf0, data1d, adios2_mode_sync, ierr)

  !13: -d(ln(tf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tfpp)/sptem(dum,tfpp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgtf_dpsi_neg, data1d, adios2_mode_sync, ierr)

  !14: nf
  do i=1,nrad
     data1d(i)=spden(datap(i),nfpp)/nepp(1,1)
  enddo
  call adios2_put(adios2_eqplot_e, var_nf0, data1d, adios2_mode_sync, ierr)

  !15: -d(ln(nf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nfpp)/spden(dum,nfpp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgnf_dpsi_neg, data1d, adios2_mode_sync, ierr)
#else
  !12: tf
  do i=1,nrad
     data1d(i)=sptem(datap(i),tepp)/(rho0*rho0)
  enddo
  call adios2_put(adios2_eqplot_e, var_Tf0, data1d, adios2_mode_sync, ierr)

  !13: -d(ln(tf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-tem_dp(dum,tepp)/sptem(dum,tepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgtf_dpsi_neg, data1d, adios2_mode_sync, ierr)

  !14: nf
  do i=1,nrad
     data1d(i)=spden(datap(i),nepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_nf0, data1d, adios2_mode_sync, ierr)

  !15: -d(ln(nf))/dr
  do i=1,nrad
     dum=datap(i)
     data1d(i)=-den_dp(dum,nepp)/spden(dum,nepp)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgnf_dpsi_neg, data1d, adios2_mode_sync, ierr)
#endif
  !16: zeff
  do i=1,nrad
     data1d(i)=spzeff(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_Zeff, data1d, adios2_mode_sync, ierr)

  !17: toroidal rotation
  do i=1,nrad
     data1d(i)=sptrot(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_u_tor, data1d, adios2_mode_sync, ierr)

  !18: radial electric field
  do i=1,nrad
     data1d(i)=sper(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_Er, data1d, adios2_mode_sync, ierr)

  !19: q profile
  do i=1,nrad
     data1d(i)=spq(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_q, data1d, adios2_mode_sync, ierr)

  !20: d(ln(q))/dpsi
  do i=1,nrad
     dum=datap(i)
     data1d(i)=dq_dp(dum)/spq(dum)
  enddo
  call adios2_put(adios2_eqplot_e, var_dlgq_dpsi, data1d, adios2_mode_sync, ierr)

  !21: gcurrent profile
  do i=1,nrad
     data1d(i)=spgpsi(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_g, data1d, adios2_mode_sync, ierr)

  !22: pressure profile
  do i=1,nrad
     data1d(i)=sppressure(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_P, data1d, adios2_mode_sync, ierr)

  !23: minor radius
  do i=1,nrad
     data1d(i)=sprpsi(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_r_1d, data1d, adios2_mode_sync, ierr)

  !24: toroidal flux
  do i=1,nrad
     data1d(i)=sptorpsi(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_torpsi, data1d, adios2_mode_sync, ierr)

  !25: radial grid: rgpsi
  do i=1,nrad
     data1d(i)=sprgpsi(datap(i))
  enddo
  call adios2_put(adios2_eqplot_e, var_r_1dg, data1d, adios2_mode_sync, ierr)

  !26: inverse of spline torpsi: psitor
  tdum=0.0
  pdum=sptorpsi(psiw)
  dum=(pdum-tdum)/real(nrad-1)
  do i=1,nrad
     pdum=tdum+dum*real(i-1)
     data1d(i)=sppsitor(pdum)
  enddo
  call adios2_put(adios2_eqplot_e, var_dtor, dum, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_psi_tor, data1d, adios2_mode_sync, ierr)

  !27: inverse of spline rgpsi: psirg
  tdum=0.0
  pdum=sprgpsi(psiw)
  dum=(pdum-tdum)/real(nrad-1)
  do i=1,nrad
     pdum=tdum+dum*real(i-1)
     data1d(i)=sppsirg(pdum)
  enddo
  call adios2_put(adios2_eqplot_e, var_drg, dum, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_psi_rg, data1d, adios2_mode_sync, ierr)

  !28: error of spline cos in [0, pi/2]
  do i=1,nrad
     dum=0.25*spdtheta*real(i-1)
     data1d(i)=splcos(dum)-cos(dum)
  enddo
  call adios2_put(adios2_eqplot_e, var_dtheta, dum, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_cos_err, data1d, adios2_mode_sync, ierr)

  !29: error of spline sin in [0, pi/2]
  do i=1,nrad
     dum=0.25*spdtheta*real(i-1)
     data1d(i)=splsin(dum)-sin(dum)
  enddo
  call adios2_put(adios2_eqplot_e, var_sin_err, data1d, adios2_mode_sync, ierr)

  !# of 2D contour plots on poloidal plane
  nplot=5
  deltheta=2.0*pi/real(lst-1)

  ! If this is changed to adios2_mode_deferred, you cannot do mpsi/mskip+1
  ! You have to store it in a temp variable and pass it to adios2_put
  call adios2_put(adios2_eqplot_e, var_npsi, mpsi/mskip+1, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_ntheta, lst, adios2_mode_sync, ierr)

  do i=1,mpsi/mskip+1
     pdum=psimesh((i-1)*mskip)
     do j=1,lst
        tdum=deltheta*real(j-1)
#ifdef _TOROIDAL3D
        datax(i,j)=spx(pdum,tdum,0.0_lk)
        dataz(i,j)=spz(pdum,tdum,0.0_lk)
        data2d(i,j,1)=spb(pdum,tdum,0.0_lk)
        data2d(i,j,4)=spf(pdum,tdum,0.0_lk)
#else
        datax(i,j)=spx(pdum,tdum)
        dataz(i,j)=spz(pdum,tdum)
        data2d(i,j,1)=spb(pdum,tdum)
        data2d(i,j,4)=spf(pdum,tdum)
#endif
        data2d(i,j,2)=spjac(pdum,tdum)
        !javierhn: currenti in Boozer coordinates is going to be always 
        ! constant in a flux surface. Not really useful here. It should be
        ! deleted. Keep it due the idl script
        data2d(i,j,3)=currenti(pdum,tdum)
        !javierhn: delb is no longer used and it has been set to zero
        data2d(i,j,5)=delb(pdum,tdum)
     enddo
  enddo

  !0-1: mesh points on (X,Z)
  call adios2_put(adios2_eqplot_e, var_R, datax, adios2_mode_sync, ierr)
  call adios2_put(adios2_eqplot_e, var_Z, dataz, adios2_mode_sync, ierr)

  !2: b-field
  call adios2_put(adios2_eqplot_e, var_B, data2d(:,:,1), adios2_mode_sync, ierr)

  !3: Jacobian
  call adios2_put(adios2_eqplot_e, var_Jacobian, data2d(:,:,2), adios2_mode_sync, ierr)

  !4: icurrent
  call adios2_put(adios2_eqplot_e, var_I, data2d(:,:,3), adios2_mode_sync, ierr)

  !5: zeta (magnetic coord.) to phi (cylindrical coord.)
  call adios2_put(adios2_eqplot_e, var_nu, data2d(:,:,4), adios2_mode_sync, ierr)

  !6: delb
  call adios2_put(adios2_eqplot_e, var_delta, data2d(:,:,5), adios2_mode_sync, ierr)

  call adios2_end_step(adios2_eqplot_e, ierr)
end subroutine eqplot_adios2
!end ifdef ADIOS2
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
