!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module cell_volume
use precision
implicit none

public get_cell_volume_MC, &
       volume
private

  
real(lk), dimension(:), allocatable :: volume
save

contains
!=====================================================
subroutine volume_init
    use grids_setup, only: mgrid
   
    allocate(volume(mgrid))
    volume=0.0_lk
end subroutine volume_init
!======================================================
subroutine get_cell_volume_MC
    use global_parameters
    use grids_setup
    use read_in_parameters, only: Zsol_domain0, Zsol_domain1
    use mapping_splines, only: sspsol_psi_z, sspcore_psi_theta, msp_psi_core, msp_psi_sol, &
                               lspsol,lszsol,del_psi_sol,del_z_sol, &
                               lspcore,lstcore,del_psi_core,del_theta_core
    use simulation_domain, only: R1sol, Z_corner_left, Z_corner_width, R1core, &
                                 psi0_sol, psi1_sol, psi0_core, psi1_core, &
                                 sfcore, sfsol
    use spline_function, only: spline2d
    use cylindrical_function, only: geotheta
    use equilibrium, only: psi_sep, lsr, lsz, spdR, spdZ, zmax, &
                             psirz_Eq
    implicit none

    integer :: m,numpart,mtest,icount,ierror,ipsi(1),ii,jj,ij,i,iz(1),s,&
               nrsp,nzsp,np,&
               jj0,jj1,ij0,ij1,iz0(1),iz1(1),is0(1),is1(1)
    real(lk) :: rdum,zdum,psidum,tdum,dpsi,dp1,sdum,ns,wt,dz,&
                ns0,ns1,dz0,dz1,&
                sdum0,sdum1,ds0,ds1,&
                dpr,dr2,dpz,dz2,psifrc,sign_psi,spdR_inv,spdZ_inv,&
                jacotmp(mgrid),adum(mgrid),weight,&
                wp1,wp0,wt10,wt00,wt11,wt01
    real(lk),dimension(:,:),allocatable :: ztest
    real(lk),dimension(:),allocatable :: wptest,wtest0,wtest1,&
                                         delscore,delscore_inv
    integer,dimension(:),allocatable :: jtest0,jtest1
  
    if(mype==0)write(gtcout,*)"numerical calculation of cell(jacobian) volume"

    call volume_init

    spdR_inv=1.0_lk/spdR
    spdZ_inv=1.0_lk/spdZ
    if(sfcore==1)then
        allocate(delscore(mpsi_core),delscore_inv(mpsi_core))
        delscore(:)=deltas_core1d(:)/totals_core(:)
        delscore_inv(:)=1.0_lk/delscore(:)
    endif
    numpart=10000000

    s=1
    do while(s < 10)
        s = s + 1
        allocate(ztest(2,numpart))
        ztest=0.0_lk
        call random_number(ztest(1,:)) !R direction
        call random_number(ztest(2,:)) !Z direction
        
        if(sfsol==1)then    
!!$omp parallel do private(m) 
            do m=1,numpart
                ztest(1,m)=R1sol*sqrt(ztest(1,m))
                ztest(2,m)=Zsol_domain0+abs(Zsol_domain1-Zsol_domain0)*ztest(2,m)
            enddo
        elseif(sfsol==0 .and. sfcore==1)then
!!$omp parallel do private(m) 
            do m=1,numpart
                ztest(1,m)=R1core*sqrt(ztest(1,m))
                ztest(2,m)=Z_corner_left+abs(Z_corner_width)*ztest(2,m)
            enddo

        else
            if(mype==0)write(gtcout,*)"wrong sfcore and sfsol settings"
            stop
        endif
 
        allocate(jtest0(numpart),jtest1(numpart),wptest(numpart),STAT=mtest)
        if(mtest /= 0) then
          print *,mype,'jtest0, jtest1 and wptest',mtest
        endif
        jtest0=0
        jtest1=0
        wptest=0.0_lk
        allocate(wtest0(numpart),wtest1(numpart),STAT=mtest)
        if(mtest /= 0) then
          print *,mype,'wtest0 and wtest1',mtest
        endif
        wtest0=0.0_lk
        wtest1=0.0_lk
     
!!$omp parallel do private(m,rdum,zdum,nrsp,dpr,dr2,nzsp,dpz,dz2,psifrc,sign_psi,np,tdum,sdum0,sdum1,&
!!$omp& jj0,jj1,ij0,ij1,iz0,iz1,dz0,dz1,ns0,ns1)
        do m=1,numpart
            rdum=ztest(1,m)
            zdum=ztest(2,m)
      
            nrsp=max(1,min(lsr-1,ceiling(rdum*spdR_inv)))
            dpr=rdum-spdR*real(nrsp-1)
            dr2=dpr*dpr
     
            nzsp=max(1,min(lsz-1,ceiling((zdum+zmax)*spdZ_inv)))
            dpz=zdum+zmax-spdZ*real(nzsp-1)
            dz2=dpz*dpz
     
            psifrc=psirz_Eq(1,nrsp,nzsp)     + psirz_Eq(2,nrsp,nzsp)*dpr     + psirz_Eq(3,nrsp,nzsp)*dr2 &
                   +psirz_Eq(4,nrsp,nzsp)*dpz + psirz_Eq(5,nrsp,nzsp)*dpr*dpz + psirz_Eq(6,nrsp,nzsp)*dr2*dpz &
                   +psirz_Eq(7,nrsp,nzsp)*dz2 + psirz_Eq(8,nrsp,nzsp)*dpr*dz2 + psirz_Eq(9,nrsp,nzsp)*dr2*dz2
     
     
            if(sfsol==1)then
        
                if(psifrc .ge. psi0_sol .and. psifrc .le. psi1_sol)then
        
                    if(psifrc .ge. psi1_sol)then
                        wptest(m)=1.0_lk
                        np=mpsi_sol-1
                    elseif(psifrc .le. psi0_sol)then
                        wptest(m)=0.0_lk
                        np=1
                    else
                        do i=1,mpsi_sol-1
                            sign_psi=(psifrc-psimesh_sol(i))*(psifrc-psimesh_sol(i+1))
                            if(sign_psi .le. 0.0_lk)then
                                np=i
                                wptest(m)=(psifrc-psimesh_sol(np))/(psimesh_sol(np+1)-psimesh_sol(np))
                                exit
                            endif
                        enddo
                    endif
        
                    sdum0=spline2d(psifrc-msp_psi_sol(1),zdum-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
                    sdum1=sdum0
            
                    is0=minloc(abs(sdum0-Ssol(igrid_sol(np)+1:igrid_sol(np)+mpara_sol(np))))
                    ds0=sdum0-Ssol(igrid_sol(np)+is0(1))
                    if(ds0>0.0_lk)then
                        jj0=min(is0(1),mpara_sol(np)-1)
                    else
                        jj0=max(is0(1)-1,1)
                    endif
            
                    is1=minloc(abs(sdum1-Ssol(igrid_sol(np+1)+1:igrid_sol(np+1)+mpara_sol(np+1))))
                    ds1=sdum1-Ssol(igrid_sol(np+1)+is1(1))
                    if(ds1>0.0_lk)then
                        jj1=min(is1(1),mpara_sol(np+1)-1)
                    else
                        jj1=max(is1(1)-1,1)
                    endif
                    ij0=igrid_sol(np)+jj0
                    ij1=igrid_sol(np+1)+jj1
            
                    wtest0(m)=(sdum0-Ssol(ij0))/(Ssol(ij0+1)-Ssol(ij0))
                    wtest1(m)=(sdum1-Ssol(ij1))/(Ssol(ij1+1)-Ssol(ij1))
                    jtest0(m)=mgrid_core+ij0
                    jtest1(m)=mgrid_core+ij1
            
                endif
            endif
        
            if(sfcore==1)then
                if(psifrc .ge. psi0_core .and. psifrc .le. psi1_core)then
            
                    if(psifrc .ge. psi1_core)then
                        wptest(m)=1.0_lk
                        np=mpsi_core-1
                    elseif(psifrc .le. psi0_core)then
                        wptest(m)=0.0_lk
                        np=1
                    else
                        do i=1,mpsi_core-1
                            sign_psi=(psifrc-psimesh_core(i))*(psifrc-psimesh_core(i+1))
                            if(sign_psi .le. 0.0_lk)then
                                np=i
                                wptest(m)=(psifrc-psimesh_core(np))/(psimesh_core(np+1)-psimesh_core(np))
                                exit
                            endif
                        enddo
                    endif
        
                    tdum=geotheta(ztest(1,m),ztest(2,m))
                    sdum0=spline2d(psifrc-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
                    sdum1=sdum0
             
                    ns0=sdum0*delscore_inv(np)
                    ns1=sdum1*delscore_inv(np+1)
             
                    jtest0(m)=igrid_core(np)+max(1,min(ceiling(ns0),mpara_core(np)-1))
                    jtest1(m)=igrid_core(np+1)+max(1,min(ceiling(ns1),mpara_core(np+1)-1))
                    if(ns0 .lt. real(mpara_core(np)-1))then
                        wtest0(m)=ns0-int(ns0)
                    else
                        wtest0(m)=1.0_lk
                    endif
             
                    if(ns1 .lt. real(mpara_core(np+1)-1))then
                        wtest1(m)=ns1-int(ns1)
                    else
                        wtest1(m)=1.0_lk
                    endif
                endif
            endif
        
        enddo
        if(mype==0)write(gtcout,*)"volume iteration(",s,")"
     
!#ifdef _OPENMP
!!$omp parallel private(jacotmp)
!        jacotmp=0.0_lk   ! Set array elements to zero
!#endif
!!$omp do private(m,rdum,zdum,nrsp,dpr,dr2,nzsp,dpz,dz2,psifrc,weight,wp1,wp0,wt10,wt00,wt11,wt01,ij)
        do m=1,numpart
     
            rdum=ztest(1,m)
            zdum=ztest(2,m)
     
            nrsp=max(1,min(lsr-1,ceiling(rdum*spdR_inv)))
            dpr=rdum-spdR*real(nrsp-1)
            dr2=dpr*dpr
     
            nzsp=max(1,min(lsz-1,ceiling((zdum+zmax)*spdZ_inv)))
            dpz=zdum+zmax-spdZ*real(nzsp-1)
            dz2=dpz*dpz
     
            psifrc=psirz_Eq(1,nrsp,nzsp)     + psirz_Eq(2,nrsp,nzsp)*dpr     + psirz_Eq(3,nrsp,nzsp)*dr2 &
                +psirz_Eq(4,nrsp,nzsp)*dpz + psirz_Eq(5,nrsp,nzsp)*dpr*dpz + psirz_Eq(6,nrsp,nzsp)*dr2*dpz &
                +psirz_Eq(7,nrsp,nzsp)*dz2 + psirz_Eq(8,nrsp,nzsp)*dpr*dz2 + psirz_Eq(9,nrsp,nzsp)*dr2*dz2
      
            weight=0.0_lk
            if(sfcore==1)then
                if(psifrc .ge. psi0_core .and. psifrc .le. psi1_core)then
                    weight=1.0_lk
                endif
            endif
      
            if(sfsol==1)then
                if(psifrc .ge. psi0_sol .and. psifrc .le. psi1_sol)then
                    weight=1.0_lk
                endif
            endif
     
     
            wp1=wptest(m)               !outer flux surface
            wp0=1.0_lk-wp1              !inner flux surface
     
            wt10=wp0*wtest0(m)       !upper Z grid on inner flux surface
            wt00=wp0-wt10            !lower Z grid on inner flux surface
     
            wt11=wp1*wtest1(m)       !upper Z grid on outer flux surface
            wt01=wp1-wt11            !lower Z grid on outer flux surface
     
!#ifdef _OPENMP
!     ! Use thread-private temporary array dnitmp to store the results
!            ij=max(jtest0(m),1)
!            jacotmp(ij) = jacotmp(ij) + wt00*weight
!     
!            ij=ij+1
!            jacotmp(ij) = jacotmp(ij) + wt10*weight
!     
!            ij=max(jtest1(m),1)
!            jacotmp(ij) = jacotmp(ij) + wt01*weight
!     
!            ij=ij+1
!            jacotmp(ij) = jacotmp(ij) + wt11*weight
!#else
     ! If no loop-level parallelism, write directly into array "densityi()".
            ij=max(jtest0(m),1)
            volume(ij) = volume(ij) + wt00*weight
     
            ij=ij+1
            volume(ij) = volume(ij) + wt10*weight
     
            ij=max(jtest1(m),1)
            volume(ij) = volume(ij) + wt01*weight
     
            ij=ij+1
            volume(ij) = volume(ij) + wt11*weight
!#endif
        enddo
!#ifdef _OPENMP
!!$omp critical
!        do ij=1,mgrid
!            volume(ij) = volume(ij) + jacotmp(ij)
!        enddo
!!$omp end critical
!#endif
!!$omp end parallel
     
        deallocate(ztest,jtest0,jtest1,wtest0,wtest1,wptest)
    enddo

! If we have a particle decomposition on the toroidal domains, do a reduce
! operation to add up all the contributions to charge density on the grid
    if(npartdom>1)then
!!$omp parallel do private(ij)
        do ij=1,mgrid
            jacotmp(ij)=volume(ij)
            volume(ij)=0.0_lk
        enddo
        call MPI_ALLREDUCE(jacotmp,volume,mgrid,mpi_Rsize,MPI_SUM,partd_comm,ierror)
    endif
    adum=0.0_lk
    call MPI_ALLREDUCE(volume,adum,mgrid,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
    volume=adum/mtoroidal

    if(sfcore==1)then
! Periodic boundary for core region
        do i=1,mpsi_core
            volume(igrid_core(i)+1)=volume(igrid_core(i)+1)+volume(igrid_core(i)+mpara_core(i))
            volume(igrid_core(i)+mpara_core(i))=volume(igrid_core(i)+1)
        enddo
        deallocate(delscore,delscore_inv)
    endif
end subroutine get_cell_volume_MC

end module cell_volume
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
