geogrid/src/parallel_module.F

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! MODULE PARALLEL_MODULE
   3 !
   4 ! This module provides routines for parallelizing.
   5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   6 module parallel_module
   7
   8 #ifdef _MPI
   9 include 'mpif.h'
  10 #endif
  11
  12    integer, parameter :: IO_NODE = 0
  13
  14    integer, parameter :: HALO_WIDTH = 3
  15
  16    integer, pointer, dimension(:,:) :: processors, &
  17                                        proc_minx, proc_maxx, &
  18                                        proc_miny, proc_maxy
  19    integer :: nprocs, &
  20               my_proc_id, &
  21               nproc_x, nproc_y, &
  22               my_x, my_y, &
  23               my_minx, my_miny, my_maxx, my_maxy, &
  24               comm
  25
  26
  27    contains
  28
  29
  30    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  31    ! Name: parallel_start
  32    !
  33    ! Purpose: For MPI, the purpose of this routine is to basically set up
  34    !   a communicator for a rectangular mesh, and determine how many processors
  35    !   in the x and y directions there will be.
  36    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  37    subroutine parallel_start()
  38
  39       implicit none
  40
  41       ! Arguments
  42
  43       ! Local variables
  44 #ifdef _MPI
  45       integer :: mpi_rank, mpi_size
  46       integer :: mpi_ierr
  47       integer, dimension(2) :: dims, coords
  48       integer :: rectangle, myleft, myright, mytop, mybottom
  49       integer :: mini, m, n
  50       logical, dimension(2) :: periods
  51
  52       ! Find out our rank and the total number of processors
  53       call MPI_Init(mpi_ierr)
  54       call MPI_Comm_rank(MPI_COMM_WORLD, mpi_rank, mpi_ierr)
  55       call MPI_Comm_size(MPI_COMM_WORLD, mpi_size, mpi_ierr)
  56
  57       comm = MPI_COMM_WORLD
  58
  59       nprocs = mpi_size
  60       my_proc_id = mpi_rank
  61
  62       ! Code from RSL to get number of procs in m and n directions
  63       mini = 2*nprocs
  64       nproc_x = 1
  65       nproc_y = nprocs
  66       do m = 1, nprocs
  67         if ( mod( nprocs, m ) == 0 ) then
  68           n = nprocs / m
  69           if ( abs(m-n) < mini  ) then
  70             mini = abs(m-n)
  71             nproc_x = m
  72             nproc_y = n
  73           end if
  74         end if
  75       end do
  76
  77       ! Calculate which patch current processor will work on
  78       my_x = mod(mpi_rank,nproc_x)
  79       my_y = mpi_rank / nproc_x
  80
  81 #else
  82       comm = 0
  83       my_proc_id = IO_NODE
  84       nprocs = 1
  85       my_x = 0
  86       my_y = 0
  87       nproc_x = 1
  88       nproc_y = 1
  89 #endif
  90
  91       nullify(processors)
  92       nullify(proc_minx)
  93       nullify(proc_maxx)
  94       nullify(proc_miny)
  95       nullify(proc_maxy)
  96
  97    end subroutine parallel_start
  98
  99
 100    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 101    ! Name: parallel_get_tile_dims
 102    !
 103    ! Purpose: To compute the starting and ending indices of the patch that the
 104    !   calling processor is to work on. When there are multiple processors,
 105    !   appropriate data structures describing the range of indices being
 106    !   worked on by other processors are also allocated and filled
 107    !   (processors, proc_minx, proc_maxx, proc_miny, proc_maxy).
 108    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 109    subroutine parallel_get_tile_dims(idim, jdim)
 110
 111       implicit none
 112
 113       ! Arguments
 114       integer, intent(in) :: idim, jdim
 115
 116       ! Local variables
 117 #ifdef _MPI
 118       integer :: i, j, ix, iy, px, py
 119       integer, dimension(2) :: buffer
 120       integer :: mpi_ierr
 121       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 122
 123       !
 124       ! Determine starting and ending grid points in x and y direction that we will work on
 125       !
 126       ! NOTE:
 127       ! For now, copy code from RSL_LITE's module_dm.F until build mechanism to link
 128       ! WRF and WPS code is worked out more.
 129       ! Eventually, it would probably be best to use module_dm code without copying
 130       !
 131       my_minx = -1
 132       j = 1
 133       do i = 1, idim
 134          call task_for_point(i, j, 1, idim, 1, jdim, nproc_x, nproc_y, px, py)
 135          if ( px == my_x ) then
 136             my_maxx = i
 137             if ( my_minx == -1 ) my_minx = i
 138          end if
 139       end do
 140
 141       my_miny = -1
 142       i = 1
 143       do j = 1, jdim
 144          call task_for_point(i, j, 1, idim, 1, jdim, nproc_x, nproc_y, px, py)
 145          if ( py == my_y ) then
 146             my_maxy = j
 147             if ( my_miny == -1 ) my_miny = j
 148          end if
 149       end do
 150
 151       ! Create space to hold information about which other processors are
 152       !   working on which parts of the domain
 153       allocate(processors(0:nproc_x-1, 0:nproc_y-1))
 154       allocate(proc_minx(0:nproc_x-1, 0:nproc_y-1))
 155       allocate(proc_miny(0:nproc_x-1, 0:nproc_y-1))
 156       allocate(proc_maxx(0:nproc_x-1, 0:nproc_y-1))
 157       allocate(proc_maxy(0:nproc_x-1, 0:nproc_y-1))
 158
 159       ! Exchange information with other processors
 160       if (my_proc_id == IO_NODE) then
 161          processors(my_x, my_y) = my_proc_id
 162          do i=1,nprocs-1
 163             call MPI_Recv(buffer, 2, MPI_INTEGER, i, MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 164             processors(buffer(1), buffer(2)) = mpi_stat(MPI_SOURCE)
 165          end do
 166       else
 167          buffer(1) = my_x
 168          buffer(2) = my_y
 169          call MPI_Send(buffer, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 170       end if
 171
 172       do ix=0,nproc_x-1
 173          do iy=0,nproc_y-1
 174             call parallel_bcast_int(processors(ix,iy), IO_NODE)
 175          end do
 176       end do
 177
 178       proc_minx(my_x, my_y) = my_minx
 179       proc_maxx(my_x, my_y) = my_maxx
 180       proc_miny(my_x, my_y) = my_miny
 181       proc_maxy(my_x, my_y) = my_maxy
 182
 183       do ix=0,nproc_x-1
 184          do iy=0,nproc_y-1
 185             call parallel_bcast_int(proc_minx(ix,iy), processors(ix,iy))
 186             call parallel_bcast_int(proc_maxx(ix,iy), processors(ix,iy))
 187             call parallel_bcast_int(proc_miny(ix,iy), processors(ix,iy))
 188             call parallel_bcast_int(proc_maxy(ix,iy), processors(ix,iy))
 189          end do
 190       end do
 191
 192 #else
 193       allocate(processors(0:nproc_x-1, 0:nproc_y-1))
 194       allocate(proc_minx(0:nproc_x-1, 0:nproc_y-1))
 195       allocate(proc_miny(0:nproc_x-1, 0:nproc_y-1))
 196       allocate(proc_maxx(0:nproc_x-1, 0:nproc_y-1))
 197       allocate(proc_maxy(0:nproc_x-1, 0:nproc_y-1))
 198
 199       processors(0,0) = IO_NODE
 200       proc_minx(0,0) = 1
 201       proc_miny(0,0) = 1
 202       proc_maxx(0,0) = idim
 203       proc_maxy(0,0) = jdim
 204       my_minx = 1
 205       my_maxx = idim
 206       my_miny = 1
 207       my_maxy = jdim
 208
 209 #endif
 210
 211    end subroutine parallel_get_tile_dims
 212
 213
 214    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 215    ! Copied from RSL_LITE's task_for_point.c until a good way can be found to
 216    !    get the build mechanism to use the original code in RSL_LITE.
 217    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 218    subroutine task_for_point(i_p, j_p, ids_p, ide_p, jds_p, jde_p, npx, npy, px, py)
 219
 220       implicit none
 221
 222       ! Arguments
 223       integer, intent(in) :: i_p, j_p, ids_p, ide_p, jds_p, jde_p, npx, npy
 224       integer, intent(out) :: px, py
 225
 226       ! Local variables
 227       integer :: a, b, rem, idim, jdim, i, j, ids, jds, ide, jde
 228
 229       i = i_p - 1
 230       j = j_p - 1
 231       ids = ids_p - 1
 232       jds = jds_p - 1
 233       ide = ide_p - 1
 234       jde = jde_p - 1
 235
 236       idim = ide-ids+1
 237       jdim = jde-jds+1
 238
 239       i = max(i,ids)
 240       i = min(i,ide)
 241       rem = mod(idim, npx)
 242       a = ( rem / 2 ) * ( (idim / npx) + 1 )
 243       b = a + ( npx - rem ) * ( idim / npx )
 244       if ( i-ids < a ) then
 245          px = (i-ids) / ( (idim / npx) + 1 )
 246       else if ( i-ids < b ) then
 247          px = ( a / ( (idim / npx) + 1 ) ) + (i-a-ids) / ( ( b - a ) / ( npx - rem ) )
 248       else
 249          px = ( a / ( (idim / npx) + 1 ) ) + (b-a-ids) / ( ( b - a ) / ( npx - rem ) ) + &
 250                                              (i-b-ids) / ( ( idim / npx ) + 1 )
 251       end if
 252
 253       j = max(j,jds)
 254       j = min(j,jde)
 255       rem = mod(jdim, npy)
 256       a = ( rem / 2 ) * ( (jdim / npy) + 1 )
 257       b = a + ( npy - rem ) * ( jdim / npy )
 258       if ( j-jds < a ) then
 259          py = (j-jds) / ( (jdim / npy) + 1 )
 260       else if ( j-jds < b ) then
 261          py = ( a / ( (jdim / npy) + 1 ) ) + (j-a-jds) / ( ( b - a ) / ( npy - rem ) )
 262       else
 263          py = ( a / ( (jdim / npy) + 1 ) ) + (b-a-jds) / ( ( b - a ) / ( npy - rem ) ) + &
 264                                              (j-b-jds) / ( ( jdim / npy ) + 1 )
 265       end if
 266
 267    end subroutine task_for_point
 268
 269
 270    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 271    ! Name: gather_whole_field_i
 272    !
 273    ! Purpose:
 274    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 275    subroutine gather_whole_field_i(patch_array, ms1, me1, ms2, me2, ms3, me3, &
 276                                  ps1, pe1, ps2, pe2, ps3, pe3, &
 277                                  domain_array, ds1, de1, ds2, de2, ds3, de3)
 278
 279       implicit none
 280
 281       ! Arguments
 282       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 283                              ms1, me1, ms2, me2, ms3, me3, &
 284                              ds1, de1, ds2, de2, ds3, de3
 285       integer, dimension(ms1:me1,ms2:me2,ms3:me3), intent(in) :: patch_array
 286       integer, dimension(ds1:de1,ds2:de2,ds3:de3), intent(inout) :: domain_array
 287
 288       ! Local variables
 289 #ifdef _MPI
 290       integer :: i, ii, j, jj, kk
 291       integer, dimension(2) :: idims, jdims
 292       integer :: mpi_ierr
 293       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 294
 295       if (my_proc_id == IO_NODE) then
 296
 297          do j=0,nproc_y-1
 298             do i=0,nproc_x-1
 299                if (processors(i,j) /= IO_NODE) then
 300                   call MPI_Recv(jdims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 301                   call MPI_Recv(idims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 302                   do kk=ds3,de3
 303 ! BUG: Check on mpi_stat and mpi_ierr
 304                      call MPI_Recv(domain_array(idims(1):idims(2),jdims(1):jdims(2),kk), &
 305                                                (idims(2)-idims(1)+1)*(jdims(2)-jdims(1)+1), &
 306                                    MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 307                   end do
 308                else
 309                   domain_array(ps1:pe1,ps2:pe2,ps3:pe3) = patch_array(ps1:pe1,ps2:pe2,ps3:pe3)
 310                end if
 311             end do
 312          end do
 313
 314       else
 315
 316          jdims(1) = ps2
 317          jdims(2) = pe2
 318          call MPI_Send(jdims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 319          idims(1) = ps1
 320          idims(2) = pe1
 321          call MPI_Send(idims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 322          do kk=ps3,pe3
 323             call MPI_Send(patch_array(ps1:pe1,ps2:pe2,kk), (pe1-ps1+1)*(pe2-ps2+1), MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 324 ! BUG: Check on mpi_ierr
 325          end do
 326      end if
 327 #endif
 328
 329    end subroutine gather_whole_field_i
 330
 331
 332    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 333    ! Name: gather_whole_field_r
 334    !
 335    ! Purpose:
 336    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 337    subroutine gather_whole_field_r(patch_array, ms1, me1, ms2, me2, ms3, me3, &
 338                                  ps1, pe1, ps2, pe2, ps3, pe3, &
 339                                  domain_array, ds1, de1, ds2, de2, ds3, de3)
 340
 341       implicit none
 342
 343       ! Arguments
 344       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 345                              ms1, me1, ms2, me2, ms3, me3, &
 346                              ds1, de1, ds2, de2, ds3, de3
 347       real, dimension(ms1:me1,ms2:me2,ms3:me3), intent(in) :: patch_array
 348       real, dimension(ds1:de1,ds2:de2,ds3:de3), intent(inout) :: domain_array
 349
 350       ! Local variables
 351 #ifdef _MPI
 352       integer :: i, ii, j, jj, kk
 353       integer, dimension(2) :: idims, jdims
 354       integer :: mpi_ierr
 355       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 356
 357       if (my_proc_id == IO_NODE) then
 358
 359          do j=0,nproc_y-1
 360             do i=0,nproc_x-1
 361                if (processors(i,j) /= IO_NODE) then
 362                   call MPI_Recv(jdims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 363                   call MPI_Recv(idims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 364                   do kk=ds3,de3
 365 ! BUG: Check on mpi_stat and mpi_ierr
 366                      call MPI_Recv(domain_array(idims(1):idims(2),jdims(1):jdims(2),kk), &
 367                                    (idims(2)-idims(1)+1)*(jdims(2)-jdims(1)+1), &
 368                                    MPI_REAL, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 369                   end do
 370                else
 371                   domain_array(ps1:pe1,ps2:pe2,ps3:pe3) = patch_array(ps1:pe1,ps2:pe2,ps3:pe3)
 372                end if
 373             end do
 374          end do
 375
 376       else
 377
 378          jdims(1) = ps2
 379          jdims(2) = pe2
 380          call MPI_Send(jdims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 381          idims(1) = ps1
 382          idims(2) = pe1
 383          call MPI_Send(idims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 384          do kk=ps3,pe3
 385             call MPI_Send(patch_array(ps1:pe1,ps2:pe2,kk), (pe1-ps1+1)*(pe2-ps2+1), MPI_REAL, 0, my_proc_id, comm, mpi_ierr)
 386 ! BUG: Check on mpi_ierr
 387          end do
 388       end if
 389 #endif
 390
 391    end subroutine gather_whole_field_r
 392
 393
 394    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 395    ! Name: scatter_whole_field_i
 396    !
 397    ! Purpose:
 398    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 399    subroutine scatter_whole_field_i(patch_array, ms1, me1, ms2, me2, ms3, me3, &
 400                                  ps1, pe1, ps2, pe2, ps3, pe3, &
 401                                  domain_array, ds1, de1, ds2, de2, ds3, de3)
 402
 403       implicit none
 404
 405       ! Arguments
 406       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 407                              ms1, me1, ms2, me2, ms3, me3, &
 408                              ds1, de1, ds2, de2, ds3, de3
 409       integer, dimension(ms1:me1,ms2:me2,ms3:me3), intent(inout) :: patch_array
 410       integer, dimension(ds1:de1,ds2:de2,ds3:de3), intent(in) :: domain_array
 411
 412       ! Local variables
 413 #ifdef _MPI
 414       integer :: i, ii, j, jj, kk
 415       integer, dimension(2) :: idims, jdims
 416       integer :: mpi_ierr
 417       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 418
 419       if (my_proc_id == IO_NODE) then
 420
 421          do j=0,nproc_y-1
 422             do i=0,nproc_x-1
 423                if (processors(i,j) /= IO_NODE) then
 424                   call MPI_Recv(jdims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 425                   call MPI_Recv(idims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 426                   do kk=ds3,de3
 427 ! BUG: Check on mpi_stat and mpi_ierr
 428                      call MPI_Send(domain_array(idims(1):idims(2),jdims(1):jdims(2),kk), &
 429                                                (idims(2)-idims(1)+1)*(jdims(2)-jdims(1)+1), &
 430                                    MPI_INTEGER, processors(i,j), my_proc_id, comm, mpi_ierr)
 431                   end do
 432                else
 433                   patch_array(ps1:pe1,ps2:pe2,ps3:pe3) = domain_array(ps1:pe1,ps2:pe2,ps3:pe3)
 434                end if
 435             end do
 436          end do
 437
 438       else
 439
 440          jdims(1) = ps2
 441          jdims(2) = pe2
 442          call MPI_Send(jdims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 443          idims(1) = ps1
 444          idims(2) = pe1
 445          call MPI_Send(idims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 446          do kk=ps3,pe3
 447             call MPI_Recv(patch_array(ps1:pe1,ps2:pe2,kk), (pe1-ps1+1)*(pe2-ps2+1), &
 448                           MPI_INTEGER, 0, MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 449 ! BUG: Check on mpi_ierr
 450          end do
 451      end if
 452 #endif
 453
 454    end subroutine scatter_whole_field_i
 455
 456
 457    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 458    ! Name: scatter_whole_field_r
 459    !
 460    ! Purpose:
 461    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 462    subroutine scatter_whole_field_r(patch_array, ms1, me1, ms2, me2, ms3, me3, &
 463                                  ps1, pe1, ps2, pe2, ps3, pe3, &
 464                                  domain_array, ds1, de1, ds2, de2, ds3, de3)
 465
 466       implicit none
 467
 468       ! Arguments
 469       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 470                              ms1, me1, ms2, me2, ms3, me3, &
 471                              ds1, de1, ds2, de2, ds3, de3
 472       real, dimension(ms1:me1,ms2:me2,ms3:me3), intent(inout) :: patch_array
 473       real, dimension(ds1:de1,ds2:de2,ds3:de3), intent(in) :: domain_array
 474
 475       ! Local variables
 476 #ifdef _MPI
 477       integer :: i, ii, j, jj, kk
 478       integer, dimension(2) :: idims, jdims
 479       integer :: mpi_ierr
 480       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 481
 482       if (my_proc_id == IO_NODE) then
 483
 484          do j=0,nproc_y-1
 485             do i=0,nproc_x-1
 486                if (processors(i,j) /= IO_NODE) then
 487                   call MPI_Recv(jdims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 488                   call MPI_Recv(idims, 2, MPI_INTEGER, processors(i,j), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 489                   do kk=ds3,de3
 490 ! BUG: Check on mpi_stat and mpi_ierr
 491                      call MPI_Send(domain_array(idims(1):idims(2),jdims(1):jdims(2),kk), &
 492                                                (idims(2)-idims(1)+1)*(jdims(2)-jdims(1)+1), &
 493                                    MPI_REAL, processors(i,j), my_proc_id, comm, mpi_ierr)
 494                   end do
 495                else
 496                   patch_array(ps1:pe1,ps2:pe2,ps3:pe3) = domain_array(ps1:pe1,ps2:pe2,ps3:pe3)
 497                end if
 498             end do
 499          end do
 500
 501       else
 502
 503          jdims(1) = ps2
 504          jdims(2) = pe2
 505          call MPI_Send(jdims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 506          idims(1) = ps1
 507          idims(2) = pe1
 508          call MPI_Send(idims, 2, MPI_INTEGER, 0, my_proc_id, comm, mpi_ierr)
 509          do kk=ps3,pe3
 510             call MPI_Recv(patch_array(ps1:pe1,ps2:pe2,kk), (pe1-ps1+1)*(pe2-ps2+1), &
 511                           MPI_REAL, 0, MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 512 ! BUG: Check on mpi_ierr
 513          end do
 514      end if
 515 #endif
 516
 517    end subroutine scatter_whole_field_r
 518
 519
 520    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 521    ! Name: exchange_halo_r
 522    !
 523    ! Purpose:
 524    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 525    subroutine exchange_halo_r(patch_array, &
 526                               ms1, me1, ms2, me2, ms3, me3, &
 527                               ps1, pe1, ps2, pe2, ps3, pe3)
 528
 529       implicit none
 530
 531       ! Arguments
 532       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 533                              ms1, me1, ms2, me2, ms3, me3
 534       real, dimension(ms1:me1,ms2:me2,ms3:me3), intent(inout) :: patch_array
 535
 536       ! Local variables
 537 #ifdef _MPI
 538       integer :: jj, kk
 539       integer :: mpi_ierr
 540       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 541
 542       !
 543       ! Get left edge of halo
 544       !
 545       if (my_x /= (nproc_x - 1)) then
 546          do kk=ps3,pe3
 547             do jj=ms2,me2
 548                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_REAL, &
 549                              processors(my_x+1,my_y), my_proc_id, comm, mpi_ierr)
 550             end do
 551          end do
 552       end if
 553       if (my_x /= 0) then
 554          do kk=ps3,pe3
 555             do jj=ms2,me2
 556                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 557                              processors(my_x-1,my_y), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 558             end do
 559          end do
 560       end if
 561
 562       !
 563       ! Get right edge of halo
 564       !
 565       if (my_x /= 0) then
 566          do kk=ps3,pe3
 567             do jj=ms2,me2
 568                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 569                              processors(my_x-1,my_y), my_proc_id, comm, mpi_ierr)
 570             end do
 571          end do
 572       end if
 573       if (my_x /= (nproc_x - 1)) then
 574          do kk=ps3,pe3
 575             do jj=ms2,me2
 576                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_REAL, &
 577                              processors(my_x+1,my_y), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 578             end do
 579          end do
 580       end if
 581
 582       !
 583       ! Get bottom edge of halo
 584       !
 585       if (my_y /= (nproc_y - 1)) then
 586          do kk=ps3,pe3
 587             do jj=pe2-HALO_WIDTH+1,pe2
 588                call MPI_Send(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_REAL, &
 589                              processors(my_x,my_y+1), my_proc_id, comm, mpi_ierr)
 590             end do
 591          end do
 592       end if
 593       if (my_y /= 0) then
 594          do kk=ps3,pe3
 595             do jj=ms2,ms2+HALO_WIDTH-1
 596                call MPI_Recv(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_REAL, &
 597                              processors(my_x,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 598             end do
 599          end do
 600       end if
 601
 602       !
 603       ! Get top edge of halo
 604       !
 605       if (my_y /= 0) then
 606          do kk=ps3,pe3
 607             do jj=ps2,ps2+HALO_WIDTH-1
 608                call MPI_Send(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_REAL, &
 609                              processors(my_x,my_y-1), my_proc_id, comm, mpi_ierr)
 610             end do
 611          end do
 612       end if
 613       if (my_y /= (nproc_y - 1)) then
 614          do kk=ps3,pe3
 615             do jj=me2-HALO_WIDTH+1,me2
 616                call MPI_Recv(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_REAL, &
 617                              processors(my_x,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 618             end do
 619          end do
 620       end if
 621
 622       !
 623       ! Get lower-right corner of halo
 624       !
 625       if (my_y /= (nproc_y - 1) .and. my_x /= 0) then
 626          do kk=ps3,pe3
 627             do jj=pe2-HALO_WIDTH+1,pe2
 628                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 629                              processors(my_x-1,my_y+1), my_proc_id, comm, mpi_ierr)
 630             end do
 631          end do
 632       end if
 633       if (my_y /= 0 .and. my_x /= (nproc_x - 1)) then
 634          do kk=ps3,pe3
 635             do jj=ms2,ms2+HALO_WIDTH-1
 636                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_REAL, &
 637                              processors(my_x+1,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 638             end do
 639          end do
 640       end if
 641
 642       !
 643       ! Get upper-left corner of halo
 644       !
 645       if (my_y /= 0 .and. my_x /= (nproc_x - 1)) then
 646          do kk=ps3,pe3
 647             do jj=ps2,ps2+HALO_WIDTH-1
 648                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_REAL, &
 649                              processors(my_x+1,my_y-1), my_proc_id, comm, mpi_ierr)
 650             end do
 651          end do
 652       end if
 653       if (my_y /= (nproc_y - 1) .and. my_x /= 0) then
 654          do kk=ps3,pe3
 655             do jj=me2-HALO_WIDTH+1,me2
 656                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 657                              processors(my_x-1,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 658             end do
 659          end do
 660       end if
 661
 662       !
 663       ! Get upper-right corner of halo
 664       !
 665       if (my_y /= 0 .and. my_x /= 0) then
 666          do kk=ps3,pe3
 667             do jj=ps2,ps2+HALO_WIDTH-1
 668                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 669                              processors(my_x-1,my_y-1), my_proc_id, comm, mpi_ierr)
 670             end do
 671          end do
 672       end if
 673       if (my_y /= (nproc_y - 1) .and. my_x /= (nproc_x - 1)) then
 674          do kk=ps3,pe3
 675             do jj=me2-HALO_WIDTH+1,me2
 676                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_REAL, &
 677                              processors(my_x+1,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 678             end do
 679          end do
 680       end if
 681
 682       !
 683       ! Get lower-left corner of halo
 684       !
 685       if (my_y /= (nproc_y - 1) .and. my_x /= (nproc_x - 1)) then
 686          do kk=ps3,pe3
 687             do jj=pe2-HALO_WIDTH+1,pe2
 688                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_REAL, &
 689                              processors(my_x+1,my_y+1), my_proc_id, comm, mpi_ierr)
 690             end do
 691          end do
 692       end if
 693       if (my_y /= 0 .and. my_x /= 0) then
 694          do kk=ps3,pe3
 695             do jj=ms2,ms2+HALO_WIDTH-1
 696                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_REAL, &
 697                              processors(my_x-1,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 698             end do
 699          end do
 700       end if
 701 #endif
 702
 703    end subroutine exchange_halo_r
 704
 705
 706    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 707    ! Name: exchange_halo_i
 708    !
 709    ! Purpose:
 710    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 711    subroutine exchange_halo_i(patch_array, &
 712                               ms1, me1, ms2, me2, ms3, me3, &
 713                               ps1, pe1, ps2, pe2, ps3, pe3)
 714
 715       implicit none
 716
 717       ! Arguments
 718       integer, intent(in) :: ps1, pe1, ps2, pe2, ps3, pe3, &
 719                              ms1, me1, ms2, me2, ms3, me3
 720       integer, dimension(ms1:me1,ms2:me2,ms3:me3), intent(inout) :: patch_array
 721
 722       ! Local variables
 723 #ifdef _MPI
 724       integer :: jj, kk
 725       integer :: mpi_ierr
 726       integer, dimension(MPI_STATUS_SIZE) :: mpi_stat
 727
 728       !
 729       ! Get left edge of halo
 730       !
 731       if (my_x /= (nproc_x - 1)) then
 732          do kk=ps3,pe3
 733             do jj=ms2,me2
 734                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 735                              processors(my_x+1,my_y), my_proc_id, comm, mpi_ierr)
 736             end do
 737          end do
 738       end if
 739       if (my_x /= 0) then
 740          do kk=ps3,pe3
 741             do jj=ms2,me2
 742                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 743                              processors(my_x-1,my_y), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 744             end do
 745          end do
 746       end if
 747
 748       !
 749       ! Get right edge of halo
 750       !
 751       if (my_x /= 0) then
 752          do kk=ps3,pe3
 753             do jj=ms2,me2
 754                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 755                              processors(my_x-1,my_y), my_proc_id, comm, mpi_ierr)
 756             end do
 757          end do
 758       end if
 759       if (my_x /= (nproc_x - 1)) then
 760          do kk=ps3,pe3
 761             do jj=ms2,me2
 762                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 763                              processors(my_x+1,my_y), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 764             end do
 765          end do
 766       end if
 767
 768       !
 769       ! Get bottom edge of halo
 770       !
 771       if (my_y /= (nproc_y - 1)) then
 772          do kk=ps3,pe3
 773             do jj=pe2-HALO_WIDTH+1,pe2
 774                call MPI_Send(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_INTEGER, &
 775                              processors(my_x,my_y+1), my_proc_id, comm, mpi_ierr)
 776             end do
 777          end do
 778       end if
 779       if (my_y /= 0) then
 780          do kk=ps3,pe3
 781             do jj=ms2,ms2+HALO_WIDTH-1
 782                call MPI_Recv(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_INTEGER, &
 783                              processors(my_x,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 784             end do
 785          end do
 786       end if
 787
 788       !
 789       ! Get top edge of halo
 790       !
 791       if (my_y /= 0) then
 792          do kk=ps3,pe3
 793             do jj=ps2,ps2+HALO_WIDTH-1
 794                call MPI_Send(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_INTEGER, &
 795                              processors(my_x,my_y-1), my_proc_id, comm, mpi_ierr)
 796             end do
 797          end do
 798       end if
 799       if (my_y /= (nproc_y - 1)) then
 800          do kk=ps3,pe3
 801             do jj=me2-HALO_WIDTH+1,me2
 802                call MPI_Recv(patch_array(ms1:me1,jj,kk), (me1-ms1+1), MPI_INTEGER, &
 803                              processors(my_x,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 804             end do
 805          end do
 806       end if
 807
 808       !
 809       ! Get lower-right corner of halo
 810       !
 811       if (my_y /= (nproc_y - 1) .and. my_x /= 0) then
 812          do kk=ps3,pe3
 813             do jj=pe2-HALO_WIDTH+1,pe2
 814                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 815                              processors(my_x-1,my_y+1), my_proc_id, comm, mpi_ierr)
 816             end do
 817          end do
 818       end if
 819       if (my_y /= 0 .and. my_x /= (nproc_x - 1)) then
 820          do kk=ps3,pe3
 821             do jj=ms2,ms2+HALO_WIDTH-1
 822                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 823                              processors(my_x+1,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 824             end do
 825          end do
 826       end if
 827
 828       !
 829       ! Get upper-left corner of halo
 830       !
 831       if (my_y /= 0 .and. my_x /= (nproc_x - 1)) then
 832          do kk=ps3,pe3
 833             do jj=ps2,ps2+HALO_WIDTH-1
 834                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 835                              processors(my_x+1,my_y-1), my_proc_id, comm, mpi_ierr)
 836             end do
 837          end do
 838       end if
 839       if (my_y /= (nproc_y - 1) .and. my_x /= 0) then
 840          do kk=ps3,pe3
 841             do jj=me2-HALO_WIDTH+1,me2
 842                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 843                              processors(my_x-1,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 844             end do
 845          end do
 846       end if
 847
 848       !
 849       ! Get upper-right corner of halo
 850       !
 851       if (my_y /= 0 .and. my_x /= 0) then
 852          do kk=ps3,pe3
 853             do jj=ps2,ps2+HALO_WIDTH-1
 854                call MPI_Send(patch_array(ps1:ps1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 855                              processors(my_x-1,my_y-1), my_proc_id, comm, mpi_ierr)
 856             end do
 857          end do
 858       end if
 859       if (my_y /= (nproc_y - 1) .and. my_x /= (nproc_x - 1)) then
 860          do kk=ps3,pe3
 861             do jj=me2-HALO_WIDTH+1,me2
 862                call MPI_Recv(patch_array(me1-HALO_WIDTH+1:me1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 863                              processors(my_x+1,my_y+1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 864             end do
 865          end do
 866       end if
 867
 868       !
 869       ! Get lower-left corner of halo
 870       !
 871       if (my_y /= (nproc_y - 1) .and. my_x /= (nproc_x - 1)) then
 872          do kk=ps3,pe3
 873             do jj=pe2-HALO_WIDTH+1,pe2
 874                call MPI_Send(patch_array(pe1-HALO_WIDTH+1:pe1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 875                              processors(my_x+1,my_y+1), my_proc_id, comm, mpi_ierr)
 876             end do
 877          end do
 878       end if
 879       if (my_y /= 0 .and. my_x /= 0) then
 880          do kk=ps3,pe3
 881             do jj=ms2,ms2+HALO_WIDTH-1
 882                call MPI_Recv(patch_array(ms1:ms1+HALO_WIDTH-1,jj,kk), HALO_WIDTH, MPI_INTEGER, &
 883                              processors(my_x-1,my_y-1), MPI_ANY_TAG, comm, mpi_stat, mpi_ierr)
 884             end do
 885          end do
 886       end if
 887 #endif
 888
 889    end subroutine exchange_halo_i
 890
 891
 892    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 893    ! Name: parallel_bcast_logical
 894    !
 895    ! Purpose:
 896    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 897    subroutine parallel_bcast_logical(lval)
 898
 899       implicit none
 900
 901       ! Argument
 902       logical, intent(inout) :: lval
 903
 904       ! Local variables
 905 #ifdef _MPI
 906       integer :: mpi_ierr
 907
 908       call MPI_Bcast(lval, 1, MPI_LOGICAL, IO_NODE, comm, mpi_ierr)
 909 #endif
 910
 911    end subroutine parallel_bcast_logical
 912
 913
 914    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 915    ! Name: parallel_bcast_int
 916    !
 917    ! Purpose:
 918    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 919    subroutine parallel_bcast_int(ival, from_whom)
 920
 921       implicit none
 922
 923       ! Argument
 924       integer, intent(inout) :: ival
 925       integer, intent(in), optional :: from_whom
 926
 927       ! Local variables
 928 #ifdef _MPI
 929       integer :: mpi_ierr
 930
 931       if (present(from_whom)) then
 932          call MPI_Bcast(ival, 1, MPI_INTEGER, from_whom, comm, mpi_ierr)
 933       else
 934          call MPI_Bcast(ival, 1, MPI_INTEGER, IO_NODE, comm, mpi_ierr)
 935       end if
 936 #endif
 937
 938    end subroutine parallel_bcast_int
 939
 940
 941    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 942    ! Name: parallel_bcast_real
 943    !
 944    ! Purpose:
 945    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 946    subroutine parallel_bcast_real(rval, from_whom)
 947
 948       implicit none
 949
 950       ! Argument
 951       real, intent(inout) :: rval
 952       integer, intent(in), optional :: from_whom
 953
 954       ! Local variables
 955 #ifdef _MPI
 956       integer :: mpi_ierr
 957
 958       if (present(from_whom)) then
 959          call MPI_Bcast(rval, 1, MPI_REAL, from_whom, comm, mpi_ierr)
 960       else
 961          call MPI_Bcast(rval, 1, MPI_REAL, IO_NODE, comm, mpi_ierr)
 962       end if
 963 #endif
 964
 965    end subroutine parallel_bcast_real
 966
 967
 968    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 969    ! Name: parallel_bcast_char
 970    !
 971    ! Purpose:
 972    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 973    subroutine parallel_bcast_char(cval, n, from_whom)
 974
 975       implicit none
 976
 977       ! Argument
 978       integer, intent(in) :: n
 979       character (len=n), intent(inout) :: cval
 980       integer, intent(in), optional :: from_whom
 981
 982       ! Local variables
 983 #ifdef _MPI
 984       integer :: mpi_ierr
 985
 986       if (present(from_whom)) then
 987          call MPI_Bcast(cval, n, MPI_CHARACTER, from_whom, comm, mpi_ierr)
 988       else
 989          call MPI_Bcast(cval, n, MPI_CHARACTER, IO_NODE, comm, mpi_ierr)
 990       end if
 991 #endif
 992
 993    end subroutine parallel_bcast_char
 994
 995
 996    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 997    ! Name: parallel_finish
 998    !
 999    ! Purpose: Free up, deallocate, and for MPI, finalize.
1000    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1001    subroutine parallel_finish()
1002
1003       implicit none
1004
1005       ! Arguments
1006
1007       ! Local variables
1008 #ifdef _MPI
1009       integer :: mpi_ierr
1010
1011       call MPI_Finalize(mpi_ierr)
1012 #endif
1013
1014       if (associated(processors)) deallocate(processors)
1015       if (associated(proc_minx)) deallocate(proc_minx)
1016       if (associated(proc_maxx)) deallocate(proc_maxx)
1017       if (associated(proc_miny)) deallocate(proc_miny)
1018       if (associated(proc_maxy)) deallocate(proc_maxy)
1019
1020    end subroutine parallel_finish
1021
1022
1023    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1024    ! Name: parallel_abort
1025    !
1026    ! Purpose: Terminate everything
1027    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1028    subroutine parallel_abort()
1029
1030       implicit none
1031
1032       ! Arguments
1033
1034       ! Local variables
1035 #ifdef _MPI
1036       integer :: mpi_ierr, mpi_errcode
1037
1038       call MPI_Abort(MPI_COMM_WORLD, mpi_errcode, mpi_ierr)
1039 #endif
1040
1041       stop
1042
1043    end subroutine parallel_abort
1044
1045 end module parallel_module