initialize_mpi Subroutine

  subroutine initialize_mpi(params)
    !! @note Subroutine that initialize MPI communications.@endnote
    !! Through this subroutine the default MPI communicator MPI_COMM_WORLD 
    !! is initialized. Also, a Cartesian
    TYPE(KORC_PARAMS), INTENT(INOUT) 	:: params
    !! Core KORC simulation parameters.
    INTEGER 				:: mpierr
    !! MPI error status.
    INTEGER, PARAMETER 			:: NDIMS = 1
    !! Number of dimensions of non-standard topology. 
    !! NDIMS=1 for a 1-D MPI topology, NDIMS=2 for a 2-D MPI topology, 
    !! and NDIMS=3 for a 3-D MPI topology.
    INTEGER, DIMENSION(:), ALLOCATABLE 	:: DIMS
    !! Dimension of the non-standard MPI topology params::mpi_params::mpi_topo. 
    !! This is equal to the number of MPI processes in KORC.
    LOGICAL 				:: all_mpis_initialized = .FALSE.
    !! Flag to determine if all the MPI processes were initialized correctly.
    LOGICAL 				:: mpi_process_initialized = .FALSE.
    !! Flag to determine if a given MPI process was initialized correctly.
    LOGICAL, PARAMETER 			:: REORDER = .FALSE.
    !! Flag to determine if the new MPI topology params::mpi_params::mpi_topo 
    !! needs to be re-ordered.
    LOGICAL, DIMENSION(:), ALLOCATABLE 	:: PERIODS !< Something here
    !! Array of logicals determining what dimensions of the new MPI 
    !! topology params::mpi_params::mpi_topo are periodic (T) or not (F).
    INTEGER :: ii
    !! Variable to iterate over different MPI processes.
    LOGICAL :: mpiinit = .FALSE.

    !call MPI_INITIALIZED(mpiinit,mpierr)
    !write(6,*) 'initialized after',mpiinit
    
    
    call MPI_INIT(mpierr)

    !write(6,*) 'mpi_init error code',mpierr
    
    !call MPI_INITIALIZED(mpiinit,mpierr)
    !write(6,*) 'initialized after',mpiinit
    

    if (mpierr .NE. MPI_SUCCESS) then
       write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
       write(6,'(/," ERROR: Initializing MPI. Aborting... ")')
       write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
       call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
    end if
    
    call MPI_INITIALIZED(mpi_process_initialized,mpierr)

    !write(6,*) 'initialized after',mpi_process_initialized
    
    call MPI_REDUCE(mpi_process_initialized,all_mpis_initialized,1, &
         MPI_LOGICAL,MPI_LAND,0,MPI_COMM_WORLD,mpierr)

    !write(6,*) 'made it here 2'
    
    call MPI_BCAST(all_mpis_initialized,1, &
         MPI_LOGICAL,0,MPI_COMM_WORLD,mpierr)

    call MPI_COMM_SIZE(MPI_COMM_WORLD, params%mpi_params%nmpi, mpierr)
    if (mpierr .NE. MPI_SUCCESS) then
       write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
       write(6,'(/," ERROR: Obtaining size of communicator. Aborting... ")')
       write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
       call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
    end if

    ! * * * Getting the rank of the MPI process in the WORLD COMMON communicator * * * !
    call MPI_COMM_RANK(MPI_COMM_WORLD, params%mpi_params%rank, mpierr)
    if (mpierr .NE. MPI_SUCCESS) then
       write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
       write(6,'(/," ERROR: Obtaining MPI rank. Aborting... ")')
       write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
       call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
    end if
    
    ! * * * Here a Cartesian topology for MPI is created * * * !
    ALLOCATE(DIMS(NDIMS))
    ALLOCATE(PERIODS(NDIMS))

    ! This loop isn't necessary but helps to do things more general in the future
    do ii=1_idef,NDIMS
       DIMS(ii) = params%mpi_params%nmpi
       PERIODS(ii) = .TRUE.
    end do

    ! * * * Here a periodic topology for MPI is created * * * !
    call MPI_CART_CREATE(MPI_COMM_WORLD, NDIMS, DIMS, PERIODS, REORDER, &
         params%mpi_params%mpi_topo, mpierr)
    if (mpierr .NE. MPI_SUCCESS) then
       write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
       write(6,'(/," ERROR: Creating new MPI topology. Aborting... ")')
       write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
       call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
    end if

    ! * * * Getting the rank of the MPI process in the new topology * * * !
    call MPI_COMM_RANK(params%mpi_params%mpi_topo, params%mpi_params%rank_topo, mpierr)
    if (mpierr .NE. MPI_SUCCESS) then
       write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
       write(6,'(/," ERROR: Obtaining new MPI topology ranks. Aborting... ")')
       write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
       call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
    end if

    DEALLOCATE(DIMS)
    DEALLOCATE(PERIODS)


    if (all_mpis_initialized) then

       call MPI_BARRIER(MPI_COMM_WORLD,mpierr)

       call set_paths(params)
       
       if (params%mpi_params%rank.EQ.0) then
          write(output_unit_write,'(/,"* * * * * COMMUNICATIONS  * * * * *")')
          write(output_unit_write,'(/,"  MPI communications initialized!  ")')
          write(output_unit_write,'(/,"  Number of MPI processes: ",I5)') params%mpi_params%nmpi
          write(output_unit_write,'(/,"* * * * * * * * * * * * * * * * * *")')
       end if

    else
       if (params%mpi_params%rank.EQ.0) then
          write(6,'(/,"* * * * * * * COMMUNICATIONS * * * * * * *")')
          write(6,'(/," ERROR: MPI not initialized. Aborting... ")')
          write(6,'(/,"* * * * * * * * * ** * * * * * * * * * * *")')
          call MPI_ABORT(MPI_COMM_WORLD, -10, mpierr)
       end if
    end if
 
    !...
  end subroutine initialize_mpi