function ind2pos(i,n,L)
integer, parameter :: dp=kind(1.d0)
integer :: i, n
real(dp) :: ind2pos, L
ind2pos = ((i-1)-(n-1)/ 2.0)*L/(n-1)
end function ind2pos

subroutine save_gnuplot(temp,n,L,filename,ix_start,iy_start,ix_size, iy_size,cartesianComm,nprocs,rank)
use mpi
! save the temperature distribution
! the ascii format is suitable for splot gnuplot function
integer, parameter                 :: dp=kind(1.d0)
integer :: ix,iy,n, ix_size, iy_size, ix_start, iy_start
integer :: ix_size_r, iy_size_r, ix_start_r, iy_start_r, nprocs,iproc,cartesianComm
real(dp), dimension(0:ix_size+1,0:iy_size+1) :: temp
real(dp), dimension(:,:), allocatable :: temp_glob
real(dp) :: ind2pos, L
character(len=*) :: filename
integer :: rank
integer :: mystatus(MPI_STATUS_SIZE)

allocate(temp_glob(n,n))

if(rank == 0) then
    temp_glob(ix_start+1:ix_start+ix_size,iy_start+1:iy_start+iy_size) = temp(1:ix_size,1:iy_size)
    do iproc = 1,nprocs-1
        call MPI_Recv(ix_start_r, 1, MPI_INTEGER, iproc, 400, cartesianComm, MPI_STATUS_IGNORE, ierr)
        call MPI_Recv(iy_start_r, 1, MPI_INTEGER, iproc, 401, cartesianComm, MPI_STATUS_IGNORE, ierr)
        call MPI_Recv(ix_size_r, 1, MPI_INTEGER, iproc, 402, cartesianComm, MPI_STATUS_IGNORE, ierr)
        call MPI_Recv(iy_size_r, 1, MPI_INTEGER, iproc, 403, cartesianComm, MPI_STATUS_IGNORE, ierr)
        call MPI_Recv(temp_glob(ix_start_r+1:ix_start_r+ix_size_r,iy_start_r+1:iy_start_r+iy_size_r), ix_size_r*iy_size_r, &
                      MPI_DOUBLE, iproc, 404, cartesianComm, MPI_STATUS_IGNORE, ierr)
    enddo
    open(unit=20,file=filename,form='formatted')
    do iy=1,n
        do ix=1,n
            write(20,*) ind2pos(ix,n,L),ind2pos(iy,n,L),temp_glob(ix,iy)
        enddo
        write(20,*)
    enddo
    close(20)
else
    call MPI_Send(ix_start, 1, MPI_INTEGER, 0, 400, cartesianComm, ierr)
    call MPI_Send(iy_start, 1, MPI_INTEGER, 0, 401, cartesianComm, ierr)
    call MPI_Send(ix_size,  1, MPI_INTEGER, 0, 402, cartesianComm, ierr)
    call MPI_Send(iy_size,  1, MPI_INTEGER, 0, 403, cartesianComm, ierr)
    call MPI_Send(temp(1:ix_size,1:iy_size), ix_size*iy_size, MPI_DOUBLE, 0, 404, cartesianComm, ierr)
endif

deallocate(temp_glob)

end subroutine save_gnuplot

subroutine init_field(temp,n,L,ix_start,iy_start,ix_size,iy_size)
! initialize the T field
integer, parameter :: dp=kind(1.d0)
integer :: ix,iy,n,ix_start, iy_start, ix_size,iy_size
real(dp), dimension(0:ix_size+1,0:iy_size+1) :: temp
real(dp) :: ind2pos, L, x, y
real(dp), parameter :: sigma = 0.1d0
real(dp), parameter :: tmax = 100.d0

do iy=0,iy_size+1
    do ix=0,ix_size+1
        x = ind2pos(ix+ix_start,n,L)
        y = ind2pos(iy+iy_start,n,L)
        temp(ix,iy) = tmax*exp(-(x**2+y**2)/(2.0*sigma**2))
    enddo
enddo
end subroutine init_field