module grids_eqdata
use precision
use grids_setup
implicit none

public
 

real(lk),dimension(:),allocatable :: psiEq0,meshni,meshne,meshti,meshte,&

                                     bmeshr_unit,bmeshz_unit,bmesh,dpsidr,dpsidz,&

                                     meshni1d,meshne1d,meshti1d,meshte1d,&

                                     kapanir,kapaner,kapatir,kapater,&

                                     kapaniz,kapanez,kapatiz,kapatez,&

                                     kapanip,kapanep,kapatip,kapatep,&

                                     jmesh,gupp,guzz,guss,gups,volume_analytic,&
! diagnostic quantities below

                                     db0mag_dr_bsp,db0mag_dz_bsp, &

                                     db0r_dr_bsp,db0r_dz_bsp, &
                                     
                                     db0z_dr_bsp,db0z_dz_bsp, &

                                     bmesh_bsp, bmeshr_bsp, bmeshz_bsp, &

                                     dbr_value_dr_psisp, dbr_value_dz_psisp, &
       
                                     dbz_value_dr_psisp, dbz_value_dz_psisp, &
               
                                     db0mag_dr_psisp, db0mag_dz_psisp, &
  
                                     grad_b_scale, curl_b_scale, curl_Bmag_scale  

! Er grids
real(lk), dimension(:), allocatable :: phip00


save
#ifndef GPU_UM
    !$acc declare create(meshni, meshne, meshti, meshte,phip00)
    !$acc declare create(meshni1d, meshti1d, meshne1d, meshte1d)
    !$acc declare create(kapanip, kapanep, kapatip, kapatep)
#endif

 
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine frc_mesh_equilibrium_field_initial
    use global_parameters,only:gtcout,mype
    use spline_function, only: spline2d,spline1d
    use equilibrium, only: ni_Eq, ne_Eq, ti_Eq, te_Eq,&
                           lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: sfcore, sfsol
    use mapping_splines, only: lsp_er, dsp_er, phip00_sp
    use read_in_parameters, only: iflux_sol,uni_profile,construct_nT,n_coef,t_coef,n0_1d,t0_1d,BC_type
!    use cylindrical_function, only: radial_decay_1d, radial_decay_2d

    integer :: i, ij, ii

!!! define 1d profile for temperature and density
    real :: t_curve_fit, n_curve_fit


!================initialize plasma profile=================== 

    allocate(meshni1d(mpsi), meshne1d(mpsi), meshti1d(mpsi), meshte1d(mpsi))


    if(construct_nT==1)then

        if(sfcore==1)then
            do i=1,mpsi_core
                meshte1d(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
   

                t_curve_fit = 1.0_lk+t_coef(1)*(tanh((t_coef(2)-psimesh_core(i))/t_coef(3))-1.0_lk)
                meshti1d(i)= t0_1d*t_curve_fit;
    
                n_curve_fit = 1.0_lk+n_coef(1)*(tanh((n_coef(2)-psimesh_core(i))/n_coef(3))-1.0_lk)
                meshni1d(i)= n0_1d*n_curve_fit;
                meshne1d(i)= n0_1d*n_curve_fit;
    
            enddo
        endif

        if(sfsol==1)then
            do i=1,mpsi_sol
                meshte1d(mpsi_core+i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
    
                t_curve_fit = 1.0_lk+t_coef(1)*(tanh((t_coef(2)-psimesh_sol(i))/t_coef(3))-1.0_lk)
                meshti1d(mpsi_core+i)= t0_1d*t_curve_fit;
    
                n_curve_fit = 1.0_lk+n_coef(1)*(tanh((n_coef(2)-psimesh_sol(i))/n_coef(3))-1.0_lk)
                meshni1d(mpsi_core+i)= n0_1d*n_curve_fit;
                meshne1d(mpsi_core+i)= n0_1d*n_curve_fit;
    
            enddo
        endif

    else
        if(sfcore==1)then
            do i=1,mpsi_core
                meshni1d(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
                meshne1d(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
                meshti1d(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
                meshte1d(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
            enddo
        endif
       
        if(sfsol==1)then
            do i=1,mpsi_sol
                meshni1d(mpsi_core+i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
                meshne1d(mpsi_core+i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
                meshti1d(mpsi_core+i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
                meshte1d(mpsi_core+i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
            enddo
        endif
    
    endif

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    allocate(meshni(mgrid),meshti(mgrid),meshte(mgrid),meshne(mgrid))
    allocate(psiEq0(mgrid),bmeshr_unit(mgrid),bmeshz_unit(mgrid),bmesh(mgrid),&
               dpsidr(mgrid),dpsidz(mgrid))
    allocate(kapanip(mpsi),kapanep(mpsi),kapatip(mpsi),kapatep(mpsi))
    allocate(kapanir(mgrid),kapaniz(mgrid),kapatir(mgrid),kapatiz(mgrid),&
            kapaner(mgrid),kapanez(mgrid),kapater(mgrid),kapatez(mgrid))
   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      Core region EQ quantities on simulation grid 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    if(sfcore==1)CALL grids_eq_inicore


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      SOL region EQ quantities on simulation grid 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if(sfsol==1)CALL grids_eq_inisol


!LOAD==1 uniform profile simulation
    if(uni_profile==1)then

! Keep original physical value of gradients when we increase the machine size
        kapanip(:) = kapanip(:)*meshni1d(:)/meshni1d(iflux)
        kapanep(:) = kapanep(:)*meshne1d(:)/meshne1d(iflux)
        kapatip(:) = kapatip(:)*meshti1d(:)/meshti1d(iflux)
        kapatep(:) = kapatep(:)*meshte1d(:)/meshte1d(iflux)

! Calculating the equilibrium kappa for turbulence drive
        kapanir(:)=kapanir(:)*meshni(:)/meshni1d(iflux)
        kapaniz(:)=kapaniz(:)*meshni(:)/meshni1d(iflux)
        kapaner(:)=kapaner(:)*meshne(:)/meshne1d(iflux)
        kapanez(:)=kapanez(:)*meshne(:)/meshne1d(iflux)
        kapatir(:)=kapatir(:)*meshti(:)/meshti1d(iflux)
        kapatiz(:)=kapatiz(:)*meshti(:)/meshti1d(iflux)
        kapater(:)=kapater(:)*meshte(:)/meshte1d(iflux)
        kapatez(:)=kapatez(:)*meshte(:)/meshte1d(iflux)

        meshni1d(:)=meshni1d(iflux)
        meshne1d(:)=meshne1d(iflux)
        meshti1d(:)=meshti1d(iflux)
        meshte1d(:)=meshte1d(iflux)
        meshni(:)=meshni1d(iflux)
        meshne(:)=meshne1d(iflux)
        meshti(:)=meshti1d(iflux)
        meshte(:)=meshte1d(iflux)
    endif

    if (BC_type==1)then
        call radial_decay_1d(kapanip)
        call radial_decay_1d(kapanep)
        call radial_decay_1d(kapatip)
        call radial_decay_1d(kapatep)
    
        call radial_decay_2d(kapanir)
        call radial_decay_2d(kapaniz)
        call radial_decay_2d(kapatir)
        call radial_decay_2d(kapatiz)
        call radial_decay_2d(kapaner)
        call radial_decay_2d(kapanez)
        call radial_decay_2d(kapater)
        call radial_decay_2d(kapatez)
    endif
#ifndef GPU_UM
    !$acc update device(meshni, meshne, meshti, meshte)
    !$acc update device(meshni1d, meshti1d, meshne1d, meshte1d)
    !$acc update device(kapanip, kapanep, kapatip, kapatep)
#endif
! gauge tensor elements calculation
    CALL setgeo
    
!==================== initialize Er =========================     
    allocate(phip00(mpsi))
    if(sfcore==1)then
        do i=1,mpsi_core
            phip00(i)=spline1d(psimesh_core(i),0,lsp_er,dsp_er,phip00_sp,0)
        enddo
    endif
   
    if(sfsol==1)then
        do i=1,mpsi_sol
            phip00(mpsi_core+i)=spline1d(psimesh_sol(i),0,lsp_er,dsp_er,phip00_sp,0)
        enddo
    endif
    !!! Rescale the phip00 profile based on the diagnosis point R and B, so the shear at 
    !!! diag point is omega_s(ii) = gamma_0*gam_ratio
    ii=mgrid_core+igrid_sol(iflux_sol)+int(0.5_lk*real(mpara_sol(iflux_sol)))  ! middle radial grid index
    ! Notice Er shear at ii is: omega_s(ii) = d2phi_dp2(ii)*R(ii)^2*B(ii)
    phip00(:) = phip00(:)*6.0_lk/Rsol(ii-mgrid_core)**2/bmesh(ii)
#ifndef GPU_UM
    !$acc update device(phip00)
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    TESTING parallel derivative by calculating b0\cdot\grad B_0   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    CALL magnetic_field_scale_length

! output key information of plasma profiles
    if(mype==0)then
        
        call eq_dataout_mype0

        if(sfcore==1)then
            write(gtcout,*)""
            write(gtcout,*)"Key information of CORE plasma profile at outer midplane:"
            write(gtcout,*)"parallel diag index/parallel grid number=",1,"/",mpara_core(1)
            write(gtcout,110)
            do i=1,mpsi_core
                ij=igrid_core(i) + 1
                write(gtcout,120)i,rg1d_core(i),psimesh_core(i),&
                                 kapanir(ij),kapatir(ij),kapaner(ij),kapater(ij)
            enddo
        endif
        if(sfsol==1)then
            write(gtcout,*)""
            write(gtcout,*)"Key information of SOL plasma profile at outer midplane:"
            write(gtcout,*)"parallel diag index/parallel grid number=",&
                            int(0.5_lk*real(mpara_sol(2))) + 1,"/",mpara_sol(2)
            write(gtcout,115)
            do i=1,mpsi_sol
                ij=mgrid_core + igrid_sol(i) + int(0.5_lk*real(mpara_sol(i))) + 1
                write(gtcout,120)i,rg1d_sol(i),psimesh_sol(i),&
                                 kapanir(ij),kapatir(ij),kapaner(ij),kapater(ij)
            enddo
        endif
    endif
110 format (1X,T4,"i",T8,"rg1d_core",T24,"psimesh_core",T40,"kapanir",T56,"kapatir",T72,"kapaner",T88,"kapater")
115 format (1X,T4,"i",T8,"rg1d_sol",T24,"psimesh_sol",T40,"kapanir",T56,"kapatir",T72,"kapaner",T88,"kapater")
120 format (1X,i4,7g16.8)

    call mode_parameter_cal

end subroutine frc_mesh_equilibrium_field_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine mode_parameter_cal
    use global_parameters, only: ulength, utime, rho0, pi, pi2, eden0, etemp0, r0, nmode, &
                                 gtcout, mype
    use read_in_parameters, only: aelectron,qion,aion, iflux_sol

    use simulation_domain, only: psi0_sol
    use grids_setup, only: psimesh_sol,iflux
    use spline_function, only: spline1d
    use mapping_splines, only: phip00_sp, lsp_er, dsp_er
    integer :: i,j,n,m,kr_number,ii
    real(lk) :: omega_temp_star,omega_ele_star, rhos, kperp2_rhos2, rhoi, kpara
    real(lk) :: gam_zeta,gam_r,temp,emfactor,canon
    real(lk) :: Er_shear,psifrc

! output the key parameters of selected modes
    if(mype==0)then

        temp=(etemp0/eden0)*7.43e2*7.43e2/(r0*r0)
        emfactor=3.0*1e10*utime/ulength
        canon=rho0*rho0/temp/aelectron/emfactor/emfactor
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        n=nmode   ! zeta mode number  (determined by toroidaln in gtc.in)
        m=1       ! Z mode number
        kr_number=1


        ii=mgrid_core+igrid_sol(iflux_sol)+int(0.5_lk*real(mpara_sol(iflux_sol)))  ! middle radial point value
        gam_zeta=real(n)/Rsol(ii-mgrid_core)
        gam_r=real(kr_number)*pi/(rg1d_sol(mpsi_sol)-rg1d_sol(1))
        kpara=2.0_lk*pi*real(m)/totals_sol(mpsi_sol)
        omega_ele_star=-meshte(ii)/bmesh(ii)*kapaner(ii)*gam_zeta/qion
        omega_temp_star=meshti(ii)*kapatir(ii)*gam_zeta/qion/bmesh(ii)
        rhos=sqrt(meshte(ii)/bmesh(ii)/bmesh(ii)*(aion/qion/qion))
        rhoi=sqrt(meshti(ii)/bmesh(ii)/bmesh(ii))*(sqrt(aion)/qion)

        ! Er shear at the diagnosis point at middle plane
        psifrc = psimesh_sol(iflux_sol)
        Er_shear = spline1d(psifrc-psi0_sol,1,lsp_er,dsp_er,phip00_sp,0)
        Er_shear = Er_shear*6.0_lk

        write(gtcout,*)'Debye length =',sqrt(temp),'r0'
        write(gtcout,*)'kzeta =',gam_zeta/(ulength*0.01_lk),'m^-1'
        write(gtcout,*)'v_thi =',rhoi*bmesh(ii)*(ulength*0.01_lk)/utime,'m/s'
        write(gtcout,*)"totals=",totals_sol(iflux_sol),'r0'
        write(gtcout,*)"vth_i/totals=",rhoi/utime/totals_sol(iflux_sol),'Hz'
        write(gtcout,*)'kpara =',kpara/(ulength*0.01_lk), 'm^-1'
        write(gtcout,*)"kapanir =",kapanir(ii)
        write(gtcout,*)"kapatir =",kapatir(ii)
        write(gtcout,*)"eta_i =",kapatir(ii)/kapanir(ii)
        write(gtcout,*)"meshte =",meshte(ii)
        write(gtcout,*)"meshti =",meshti(ii)
        write(gtcout,*)"omega_temp_star =",omega_temp_star,'=',omega_temp_star/(pi2*utime),'Hz'
        write(gtcout,*)"omega_ele_star=",omega_ele_star,'=',omega_ele_star/(pi2*utime),'Hz'
        write(gtcout,*)"rhos=",rhos,'r0'
        write(gtcout,*)"rhoi=",rhoi,'r0'
        write(gtcout,*)"kperp_rhos=",gam_zeta*rhos
        write(gtcout,*)"kperp_rhoi=",gam_zeta*rhoi
        write(gtcout,*)"bmesh =", bmesh(ii)
        write(gtcout,*)"Rsol(iflux_sol) =", Rsol(igrid_sol(iflux_sol)+int(0.5_lk*real(mpara_sol(iflux_sol))))
        write(gtcout,*)"Er shear =", Er_shear,'omega_cp (i.e. 1.0/utime)'
        call FLUSH(gtcout)
    endif


end subroutine mode_parameter_cal
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function diff_1d_profile(n,f,x)result(diff_f)
    integer, intent(in) :: n
    real(lk), intent(in) :: f(n),x(n)
    real(lk) :: diff_f(n)
!    real(lk), intent(out) :: diff_f(n)

    integer :: i

    do i=2,n-1
        diff_f(i) = (f(i+1)-f(i-1))/(x(i+1)-x(i-1))
    enddo
    diff_f(1) = (f(2)-f(1))/(x(2)-x(1))
    diff_f(n) = (f(n)-f(n-1))/(x(n)-x(n-1))
end function diff_1d_profile


subroutine grids_eq_inicore
    use read_in_parameters,only: min_r,uni_profile,construct_nT
    use global_parameters,only:gtcout,size_scaling,mype
    use spline_function, only: spline2d
    use equilibrium, only: ni_Eq, ne_Eq, ti_Eq, te_Eq, psirz_Eq, &
                           lsr, lsz, spdR, spdZ, zmax

    integer :: i,j,ij,ii,mtest,icount,ierror
    real(lk) :: nitmp, netmp, titmp, tetmp

    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j
            meshni(ij)=meshni1d(i)
            meshne(ij)=meshne1d(i)
            meshti(ij)=meshti1d(i)
            meshte(ij)=meshte1d(i)
        enddo
    enddo

! Calculate bmesh and psi on the simulation grids_setup
!$omp parallel do private(i) 
    do i=1,mgrid_core
        psiEq0(i) = spline2d(Rcore(i),Zcore(i)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        dpsidr(i) = spline2d(Rcore(i),Zcore(i)+zmax,1,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        dpsidz(i) = spline2d(Rcore(i),Zcore(i)+zmax,2,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    enddo

!$omp parallel do private(i,j,ij) 
    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j
            bmeshr_unit(ij)=-dpsidz(ij)/max(Rcore(ij),min_r)
            bmeshz_unit(ij)=dpsidr(ij)/max(Rcore(ij),min_r)
            bmesh(ij)=sqrt(bmeshr_unit(ij)**2+bmeshz_unit(ij)**2)
            bmeshr_unit(ij)=bmeshr_unit(ij)/bmesh(ij)
            bmeshz_unit(ij)=bmeshz_unit(ij)/bmesh(ij)
        enddo
    enddo

! KEEP REALISTIC KAPPA, DON't USE equilibrium quantities based on simulation grids
! Calculate the \partial /\partial \psi by using the middle plane value
!$omp parallel do private(i,ij,nitmp,netmp,titmp,tetmp) 
    do i=1,mpsi_core
        ij=igrid_core(i)+1
        kapanir(ij)=spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
        kapaner(ij)=spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
        kapatir(ij)=spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
        kapater(ij)=spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,te_Eq,0,0)
  
        nitmp = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
        netmp = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
        titmp = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
        tetmp = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
 
        kapanip(i)=kapanir(ij)/dpsidr(ij)/nitmp
        kapanep(i)=kapaner(ij)/dpsidr(ij)/netmp
        kapatip(i)=kapatir(ij)/dpsidr(ij)/titmp
        kapatep(i)=kapater(ij)/dpsidr(ij)/tetmp
    enddo

!!! compute kappa value based on 1d profile:
    if (construct_nT==1) then
!    !!! modify the kappa value to achieve constant omega* theoretically:
!    !!!       omega* = (kappani_r+kappati_r)*Ti/B*(nmode/R)
        kapatip(1:mpsi_core)=diff_1d_profile(mpsi_core,log(meshti1d(1:mpsi_core)),psimesh_core)
        kapatep(1:mpsi_core)=diff_1d_profile(mpsi_core,log(meshte1d(1:mpsi_core)),psimesh_core)
        kapanip(1:mpsi_core)=diff_1d_profile(mpsi_core,log(meshni1d(1:mpsi_core)),psimesh_core)
        kapanep(1:mpsi_core)=diff_1d_profile(mpsi_core,log(meshne1d(1:mpsi_core)),psimesh_core)
    endif

! Keep original physical value of kappa when we increase the machine size
    do i=1,mpsi_core
        kapanip(i) = size_scaling*kapanip(i)!*1.6195_lk
        kapanep(i) = size_scaling*kapanep(i)!*1.6195_lk
        kapatip(i) = size_scaling*kapatip(i)!*0.8035_lk
        kapatep(i) = size_scaling*kapatep(i)!*0.8035_lk
    enddo

! Calculating the equilibrium kappa for turbulence drive
!$omp parallel do private(i,j,ij) 
    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j
            kapanir(ij)=kapanip(i)*dpsidr(ij)
            kapaniz(ij)=kapanip(i)*dpsidz(ij)
            kapatir(ij)=kapatip(i)*dpsidr(ij)
            kapatiz(ij)=kapatip(i)*dpsidz(ij)
            kapaner(ij)=kapanep(i)*dpsidr(ij)
            kapanez(ij)=kapanep(i)*dpsidz(ij)
            kapater(ij)=kapatep(i)*dpsidr(ij)
            kapatez(ij)=kapatep(i)*dpsidz(ij)
        enddo
    enddo

end subroutine grids_eq_inicore

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine grids_eq_inisol
    use read_in_parameters,only: min_r,uni_profile,construct_nT
    use global_parameters,only:gtcout,size_scaling,mype
    use spline_function, only: spline2d
    use equilibrium, only: ni_Eq, ne_Eq, ti_Eq, te_Eq, psirz_Eq,&
                           lsr, lsz, spdR, spdZ, zmax

    integer :: i,j,ij,ijsol,ii,mtest,icount,ierror
    real(lk) :: nitmp, netmp, titmp, tetmp

    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=mgrid_core+igrid_sol(i)+j
            meshni(ij)=meshni1d(mpsi_core+i)
            meshne(ij)=meshne1d(mpsi_core+i)
            meshti(ij)=meshti1d(mpsi_core+i)
            meshte(ij)=meshte1d(mpsi_core+i)
        enddo
    enddo

!$omp parallel do private(i) 
    do i=1,mgrid_sol
        psiEq0(mgrid_core+i) = spline2d(Rsol(i),Zsol(i)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        dpsidr(mgrid_core+i) = spline2d(Rsol(i),Zsol(i)+zmax,1,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        dpsidz(mgrid_core+i) = spline2d(Rsol(i),Zsol(i)+zmax,2,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    enddo

!$omp parallel do private(i,j,ij,ijsol) 
    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=mgrid_core+igrid_sol(i)+j
            ijsol=igrid_sol(i)+j
            bmeshr_unit(ij)=-dpsidz(ij)/max(Rsol(ijsol),min_r)
            bmeshz_unit(ij)=dpsidr(ij)/max(Rsol(ijsol),min_r)
            bmesh(ij)=sqrt(bmeshr_unit(ij)**2+bmeshz_unit(ij)**2)
            bmeshr_unit(ij)=bmeshr_unit(ij)/bmesh(ij)
            bmeshz_unit(ij)=bmeshz_unit(ij)/bmesh(ij)
        enddo
    enddo

! KEEP REALISTIC KAPPA, DON't USE equilibrium quantities on simulation grids
! Calculate the partial partial \psi by using the middle plane value

!$omp parallel do private(i,ij,ijsol,nitmp,netmp,titmp,tetmp) 
    do i=1,mpsi_sol
        ij=mgrid_core+igrid_sol(i)+int(0.5_lk*real(mpara_sol(i)+1))
        ijsol=igrid_sol(i)+int(0.5_lk*real(mpara_sol(i)+1))
        kapanir(ij)=spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
        kapaner(ij)=spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
        kapatir(ij)=spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
        kapater(ij)=spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,te_Eq,0,0)
 
        nitmp = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,ni_Eq,0,0)
        netmp = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,ne_Eq,0,0)
        titmp = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,ti_Eq,0,0)
        tetmp = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,te_Eq,0,0)
 
        kapanip(mpsi_core+i)=kapanir(ij)/dpsidr(ij)/nitmp
        kapanep(mpsi_core+i)=kapaner(ij)/dpsidr(ij)/netmp
        kapatip(mpsi_core+i)=kapatir(ij)/dpsidr(ij)/titmp
        kapatep(mpsi_core+i)=kapater(ij)/dpsidr(ij)/tetmp
    enddo


!!! compute kappa value based on 1d profile:
    if (construct_nT==1) then
!    !!! modify the kappa value to achieve constant omega* theoretically:
!    !!!       omega* = (kappani_r+kappati_r)*Ti/B*(nmode/R)
        kapatip(mpsi_core+1:mpsi)=diff_1d_profile(mpsi_sol,log(meshti1d(mpsi_core+1:mpsi)),psimesh_sol)
        kapatep(mpsi_core+1:mpsi)=diff_1d_profile(mpsi_sol,log(meshte1d(mpsi_core+1:mpsi)),psimesh_sol)
        kapanip(mpsi_core+1:mpsi)=diff_1d_profile(mpsi_sol,log(meshni1d(mpsi_core+1:mpsi)),psimesh_sol)
        kapanep(mpsi_core+1:mpsi)=diff_1d_profile(mpsi_sol,log(meshne1d(mpsi_core+1:mpsi)),psimesh_sol)
    endif


! Keep original physical value of kappa when we increase the machine size
    do i=1,mpsi_sol
        kapanip(mpsi_core+i) = size_scaling*kapanip(mpsi_core+i)!*1.6195_lk
        kapanep(mpsi_core+i) = size_scaling*kapanep(mpsi_core+i)!*1.6195_lk
        kapatip(mpsi_core+i) = size_scaling*kapatip(mpsi_core+i)!*0.8035_lk
        kapatep(mpsi_core+i) = size_scaling*kapatep(mpsi_core+i)!*0.8035_lk
    enddo

! Calculating the equilibrium kappa for turbulence drive
!$omp parallel do private(i,j,ij) 
    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=mgrid_core+igrid_sol(i)+j
            kapanir(ij)=kapanip(mpsi_core+i)*dpsidr(ij)
            kapaniz(ij)=kapanip(mpsi_core+i)*dpsidz(ij)
            kapatir(ij)=kapatip(mpsi_core+i)*dpsidr(ij)
            kapatiz(ij)=kapatip(mpsi_core+i)*dpsidz(ij)
            kapaner(ij)=kapanep(mpsi_core+i)*dpsidr(ij)
            kapanez(ij)=kapanep(mpsi_core+i)*dpsidz(ij)
            kapater(ij)=kapatep(mpsi_core+i)*dpsidr(ij)
            kapatez(ij)=kapatep(mpsi_core+i)*dpsidz(ij)
        enddo
    enddo

end subroutine grids_eq_inisol
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine setgeo
    use global_parameters,only:gtcout,mype
    use mapping_splines, only: sspcore_psi_theta, sspsol_psi_z, msp_psi_core, msp_psi_sol, &
                               lspcore,lstcore,del_psi_core,del_theta_core, &
                               lspsol,lszsol,del_psi_sol,del_z_sol
    use simulation_domain, only: sfcore, sfsol
    use spline_function, only: spline2d,sppsiR1d
    use cylindrical_function, only: geotheta,dgeotheta
    use equilibrium, only: psirz_Eq,&
                           lsr,lsz,spdR,spdZ, zmax, psi_sep, r_Eq, z_Eq
    use cell_volume, only: get_cell_volume_MC
    use read_in_parameters,only: min_r, Zsol_domain0

    integer i,j,ij,ijsol,mtest
    real(lk) rdum,zdum,rdum_inv,tdum,dtdr,dtdz,dsdp,dsdt,dsdz,&
             dpsidr,dpsidz,del_s,del_p, &
             jaco_frc,psidum,psidumtmp


    allocate(gupp(mgrid),guzz(mgrid),guss(mgrid),gups(mgrid),jmesh(mgrid),volume_analytic(mgrid),STAT=mtest)
    gupp=0.0_lk
    guzz=0.0_lk
    guss=0.0_lk
    gups=0.0_lk
    jmesh=0.0_lk
    volume_analytic=0.0_lk
    
    if(sfcore==1)then
        do i=1,mpsi_core
            del_p=0.5_lk*(psimesh_core(min(mpsi_core,i+1))-psimesh_core(i))+0.5_lk*(psimesh_core(i)-psimesh_core(max(1,i-1)))
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j

                rdum = max(min_r, abs(Rcore(ij)))
                rdum_inv = 1.0_lk/rdum
                zdum=Zcore(ij)
                tdum=geotheta(rdum,zdum)
      
                dtdr=dgeotheta(1,rdum,zdum)
                dtdz=dgeotheta(2,rdum,zdum)
                dsdp=spline2d(psimesh_core(i)-msp_psi_core(1),tdum,1,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
                dsdt=spline2d(psimesh_core(i)-msp_psi_core(1),tdum,2,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
                dpsidr=spline2d(rdum,zdum+zmax,1,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
                dpsidz=spline2d(rdum,zdum+zmax,2,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    
                jmesh(ij)=rdum_inv*dpsidr*dsdt*dtdz-rdum_inv*dpsidz*dsdt*dtdr
                gupp(ij)=dpsidr*dpsidr+dpsidz*dpsidz
                guzz(ij)=rdum_inv*rdum_inv
                guss(ij)=(dsdp*dpsidr+dsdt*dtdr)**2+(dsdp*dpsidz+dsdt*dtdz)**2
                gups(ij)=dsdp*(dpsidr*dpsidr+dpsidz*dpsidz)+dsdt*(dpsidr*dtdr+dpsidz*dtdz)
    
                if(j==1 .or. j==mpara_core(i))then
                    del_s=0.5_lk*(Score(igrid_core(i)+2)-Score(igrid_core(i)+1)) + 0.5_lk*(Score(igrid_core(i)+mpara_core(i))-Score(igrid_core(i)+mpara_core(i)-1))
                else
                    del_s=0.5_lk*(Score(ij+1)-Score(ij))+0.5_lk*(Score(ij)-Score(ij-1))
                endif

                jaco_frc = 1.0_lk/jmesh(ij)
                volume_analytic(ij)=jaco_frc*del_p*del_s
            enddo
        enddo    
    endif
    
    
    call corepolsmooth(jmesh)
    call corepolsmooth(guss)
    call corepolsmooth(gups)
    
    
    if(sfsol==1)then
        do i=1,mpsi_sol
            del_p=0.5_lk*(psimesh_sol(min(mpsi_sol,i+1))-psimesh_sol(i))+0.5_lk*(psimesh_sol(i)-psimesh_sol(max(1,i-1)))

            do j=1,mpara_sol(i)
                ijsol=igrid_sol(i)+j
                ij=mgrid_core+ijsol
            
                if(psimesh_sol(i) .gt. 1.01_lk*psi_sep)then
                    psidum = psimesh_sol(i)
                    rdum = max(min_r, Rsol(ijsol))
                else
! use psi in the cell center for separatrix to calculate jacobian        
                    psidum = psimesh_sol(i) + del_p
                    call sppsiR1d(0,lsr,lsz,spdR,spdZ,r_Eq,z_Eq,psirz_Eq,Zsol(ijsol)+zmax,psidum,rdum,psidumtmp)
                endif
                rdum_inv = 1.0_lk/rdum
        
                dsdp=spline2d(psidum-msp_psi_sol(1),Zsol(ijsol)-Zsol_domain0,1,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
                dsdz=spline2d(psidum-msp_psi_sol(1),Zsol(ijsol)-Zsol_domain0,2,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
                dpsidr=spline2d(rdum,Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
                dpsidz=spline2d(rdum,Zsol(ijsol)+zmax,2,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

                jmesh(ij)=rdum_inv*dpsidr*dsdz
                gupp(ij)=dpsidr*dpsidr + dpsidz*dpsidz
                guzz(ij)=rdum_inv*rdum_inv
                guss(ij)=(dsdp*dpsidr)**2 + (dsdp*dpsidz+dsdz)**2
                gups(ij)=dsdp*gupp(ij) + dpsidz*dsdz
        
                if(j==1)then
                    del_s=0.5_lk*(Ssol(ijsol+1)-Ssol(ijsol))
                elseif(j==mpara_sol(i))then
                    del_s=0.5_lk*(Ssol(ijsol)-Ssol(ijsol-1))
                else
                    del_s=0.5_lk*(Ssol(ijsol+1)-Ssol(ijsol))+0.5_lk*(Ssol(ijsol)-Ssol(ijsol-1))
                endif
          
                jaco_frc = 1.0_lk/jmesh(ij)
                volume_analytic(ij) = jaco_frc*del_p*del_s
            enddo
        enddo
    endif

    call get_cell_volume_MC

end subroutine setgeo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine magnetic_field_scale_length
    use global_parameters,only:gtcout,mype,min_b
    use spline_function, only: spline2d
    use equilibrium, only: B0Mag_Eq, B0R_Eq, B0Z_Eq, psirz_Eq,&
                           lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: sfcore, sfsol
    use grids_fields,only: gradient
    use read_in_parameters, only: min_r

    integer :: i,j,ij,ijsol,ii,mtest,icount,ierror
    real(lk) rdum, b_inv

    allocate(bmesh_bsp(mgrid),bmeshr_bsp(mgrid),bmeshz_bsp(mgrid),&
             db0mag_dr_bsp(mgrid),db0mag_dz_bsp(mgrid),&
             db0r_dr_bsp(mgrid),db0r_dz_bsp(mgrid),&
             db0z_dr_bsp(mgrid),db0z_dz_bsp(mgrid))

    allocate(dbr_value_dr_psisp(mgrid), dbr_value_dz_psisp(mgrid), &
             dbz_value_dr_psisp(mgrid), dbz_value_dz_psisp(mgrid), &
             db0mag_dr_psisp(mgrid), db0mag_dz_psisp(mgrid))   

    allocate(grad_b_scale(mgrid),curl_b_scale(mgrid), curl_Bmag_scale(mgrid))

    if(sfcore==1)then
!$omp parallel do private(i,j,ij,rdum, b_inv) 
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                rdum = max(min_r, Rcore(ij))
                b_inv = 1.0_lk/max(min_b,bmesh(ij))

                bmesh_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)
                bmeshr_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                bmeshz_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
                db0mag_dr_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)
                db0mag_dz_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,2,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)

                db0r_dr_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                db0r_dz_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,2,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                db0z_dr_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,1,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
                db0z_dz_bsp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,2,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)

                dbr_value_dr_psisp(ij) = -spline2d(Rcore(ij),Zcore(ij)+zmax,21,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbr_value_dz_psisp(ij) = -spline2d(Rcore(ij),Zcore(ij)+zmax,22,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbz_value_dr_psisp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,11,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbz_value_dz_psisp(ij) = spline2d(Rcore(ij),Zcore(ij)+zmax,12,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                db0mag_dr_psisp(ij) = bmeshr_unit(ij)*dbr_value_dr_psisp(ij) + &
                                      bmeshz_unit(ij)*dbz_value_dr_psisp(ij)
                db0mag_dz_psisp(ij) = bmeshr_unit(ij)*dbr_value_dz_psisp(ij) + &
                                      bmeshz_unit(ij)*dbz_value_dz_psisp(ij)
                grad_b_scale(ij) = (bmeshz_unit(ij)*db0mag_dr_psisp(ij)-bmeshr_unit(ij)*db0mag_dz_psisp(ij))*b_inv
                curl_b_scale(ij) = grad_b_scale(ij) + (dbr_value_dz_psisp(ij)-dbz_value_dr_psisp(ij))*b_inv
                curl_Bmag_scale(ij) =  (dbr_value_dz_psisp(ij)-dbz_value_dr_psisp(ij))*b_inv
            enddo
        enddo
    endif
    
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,ijsol,rdum,b_inv) 
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                ijsol=igrid_sol(i)+j
                rdum = max(min_r, Rsol(ijsol))
                b_inv = 1.0_lk/max(min_b,bmesh(ij))

                bmesh_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)
                bmeshr_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                bmeshz_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,0,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
                db0mag_dr_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)
                db0mag_dz_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,2,lsr,lsz,spdR,spdZ,B0Mag_Eq,0,0)

                db0r_dr_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                db0r_dz_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,2,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
                db0z_dr_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,1,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
                db0z_dz_bsp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,2,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)

                dbr_value_dr_psisp(ij) = -spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,21,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbr_value_dz_psisp(ij) = -spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,22,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbz_value_dr_psisp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,11,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                dbz_value_dz_psisp(ij) = spline2d(Rsol(ijsol),Zsol(ijsol)+zmax,12,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)/rdum
                db0mag_dr_psisp(ij) = bmeshr_unit(ij)*dbr_value_dr_psisp(ij) + &
                                      bmeshz_unit(ij)*dbz_value_dr_psisp(ij)
                db0mag_dz_psisp(ij) = bmeshr_unit(ij)*dbr_value_dz_psisp(ij) + &
                                      bmeshz_unit(ij)*dbz_value_dz_psisp(ij)
                grad_b_scale(ij) = (bmeshz_unit(ij)*db0mag_dr_psisp(ij)-bmeshr_unit(ij)*db0mag_dz_psisp(ij))*b_inv
                curl_b_scale(ij) = grad_b_scale(ij) + (dbr_value_dz_psisp(ij)-dbz_value_dr_psisp(ij))*b_inv
                curl_Bmag_scale(ij) =  (dbr_value_dz_psisp(ij)-dbz_value_dr_psisp(ij))*b_inv
            enddo
        enddo
    endif

endsubroutine magnetic_field_scale_length
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine eq_dataout_mype0
    use read_in_parameters, only: numerfrc
    implicit none
  
    integer igrid_eq, i, j, ij

    igrid_eq=100
    open(igrid_eq,file='grids_eq.out',status='replace')
    if(numerfrc==1)then
        write(igrid_eq,101)mpsi_core,mpara_core(1),mgrid_core
        write(igrid_eq,101)mpsi_sol,mpara_sol(1),mgrid_sol

        write(igrid_eq,102)psimesh_core,psimesh_sol
        write(igrid_eq,102)kapanip,kapatip,kapanep,kapatep

        write(igrid_eq,102)psiEq0, dpsidr, dpsidz 
        write(igrid_eq,102)bmeshr_unit, bmeshz_unit, bmesh

        write(igrid_eq,102)bmeshr_bsp, bmeshz_bsp, bmesh_bsp
        write(igrid_eq,102)db0mag_dr_bsp, db0mag_dz_bsp
        write(igrid_eq,102)db0r_dr_bsp, db0r_dz_bsp, db0z_dr_bsp, db0z_dz_bsp

        write(igrid_eq,102)grad_b_scale, curl_b_scale, curl_Bmag_scale
        write(igrid_eq,102)meshni, meshti, meshne, meshte 
        write(igrid_eq,102)kapanir, kapaniz, kapatir, kapatiz 
        write(igrid_eq,102)kapaner, kapanez, kapater, kapatez 
        
    elseif(numerfrc==2)then
        write(igrid_eq,101)mpsi_core,mpara_core(1),mgrid_core
       
        write(igrid_eq,102)psimesh_core
        write(igrid_eq,102)kapanip,kapatip,kapanep,kapatep

        write(igrid_eq,102)psiEq0, dpsidr, dpsidz
        write(igrid_eq,102)bmeshr_unit, bmeshz_unit, bmesh
        write(igrid_eq,102)grad_b_scale, curl_b_scale, curl_Bmag_scale
        write(igrid_eq,102)meshni, meshti, meshne, meshte 
        write(igrid_eq,102)kapanir, kapaniz, kapatir, kapatiz 
        write(igrid_eq,102)kapaner, kapanez, kapater, kapatez 
                  
    elseif(numerfrc==3 .or. numerfrc==4)then
        write(igrid_eq,101)mpsi_sol,mpara_sol(1),mgrid_sol

        write(igrid_eq,102)psimesh_sol
        write(igrid_eq,102)kapanip,kapatip,kapanep,kapatep

        write(igrid_eq,102)psiEq0, dpsidr, dpsidz
        write(igrid_eq,102)bmeshr_unit, bmeshz_unit, bmesh
        write(igrid_eq,102)grad_b_scale, curl_b_scale, curl_Bmag_scale
        write(igrid_eq,102)meshni, meshti, meshne, meshte 
        write(igrid_eq,102)kapanir, kapaniz, kapatir, kapatiz 
        write(igrid_eq,102)kapaner, kapanez, kapater, kapatez 
           
    endif                 

    close(igrid_eq)

101 format(i6)
102 format(e13.6)

end subroutine eq_dataout_mype0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine corepolsmooth(farray)
    use global_parameters
    use grids_setup
    use simulation_domain, only: sfcore
  
    integer :: i,j,ij,ijn,ij1,ij2,ij3,ij4,icount
    real(lk),dimension(mgrid_core),intent(inout) :: farray(mgrid_core)
    real(lk),dimension(mgrid_core) :: farraytmp

!$omp parallel do private(i)
    do i=1,mgrid_core
        farraytmp(i)=farray(i)
    enddo
  
!$omp parallel do private(i)
    do i=1,mpsi_core
        farraytmp(igrid_core(i)+1)=0.5_lk*(farray(igrid_core(i)+1)+farray(igrid_core(i)+mpara_core(i)))
        farraytmp(igrid_core(i)+mpara_core(i))=0.5_lk*(farray(igrid_core(i)+1)+farray(igrid_core(i)+mpara_core(i)))

        farray(igrid_core(i)+1)=farraytmp(igrid_core(i)+1)
        farray(igrid_core(i)+mpara_core(i))=farraytmp(igrid_core(i)+mpara_core(i)) 
    enddo

! Constant deltas grid!
! poloidal smoothing (-0.0625 0.25 0.625 0.25 -0.0625)
!$omp parallel do private(i,j,ij,ij1,ij2,ij3,ij4)
    do i=1,mpsi_core
        do j=3,mpara_core(i)-2
            ij=igrid_core(i)+j
            ij1=ij+1
            ij2=ij+2
            ij3=ij-1
            ij4=ij-2

            farraytmp(ij)=0.625_lk*farray(ij)+0.25_lk*(farray(ij1)+farray(ij3))-&
                          0.0625_lk*(farray(ij2)+farray(ij4))
        enddo
    enddo

!$omp parallel do private(i,ij,ijn)
    do i=1,mpsi_core
        ij=igrid_core(i)+1
        ijn=igrid_core(i)+mpara_core(i)

        farraytmp(ij)=0.625_lk*farray(ij)+0.25_lk*(farray(ij+1)+farray(ijn-1))-&
                      0.0625_lk*(farray(ij+2)+farray(ijn-2))
        farraytmp(ij+1)=0.625_lk*farray(ij+1)+0.25_lk*(farray(ij+2)+farray(ij))-&
                        0.0625_lk*(farray(ij+3)+farray(ijn-1))
        farraytmp(ijn)=0.625_lk*farray(ijn)+0.25_lk*(farray(ij+1)+farray(ijn-1))-&
                       0.0625_lk*(farray(ij+2)+farray(ijn-2))
        farraytmp(ijn-1)=0.625_lk*farray(ijn-1)+0.25_lk*(farray(ijn)+farray(ijn-2))-&
                         0.0625_lk*(farray(ij+1)+farray(ijn-3))
    enddo
  
!$omp parallel do private(i)
    do i=1,mgrid_core
       farray(i)=farraytmp(i)
    enddo
end subroutine corepolsmooth
!================================================================  
end module grids_eqdata
