program mpizetasum
** mpi control declarations  
      include 'mpif.h'
      integer numprocs, numslaves, slave, myrank, ierr, serr, rerr, tag
      integer master /0/   
      integer status(MPI_STATUS_SIZE)
** other declarations    
      integer indices(2),integer n,termsperproc
      double precision zetasum, partialsum  
** universal parameter and data initalizations 
      numterms = 1000000
      s = 2
** mpi driver
      tag =1
      call mpi_init(ierr)
      call mpi_comm_rank(MPI_COMM_WORLD, myrank, ierr)
      call mpi_comm_size(MPI_COMM_WORLD, numprocs, ierr)
      numslaves = numprocs-1
      if (myrank .EQ. master) then
         goto 200 
       else
         goto 300
       endif
** mpi cleanup
100   if (myrank .EQ. master) then
        print *, "zetasum = ", zetasum
        endif
      call mpi_finalize(ierr)
      goto 400
** master node code
200   termsperproc = numterms/numslaves
      lastslaveterm = termsperproc*numslaves
      do slave = 1,numslaves
         indices(1) = termsperproc*(slave-1)+1
         indices(2) = termsperproc*slave
         call mpi_send(indices,2,MPI_INTEGER,slave,tag,MPI_COMM_WORLD,serr)    
         enddo  
      do n = lastslaveterm+1, numterms
         zetasum = zetasum + 1/n**s
         enddo
      print *, "tailsum = ", zetasum   
      do slave = 1,numslaves
        call mpi_recv(partialsum,1,MPI_DOUBLE_PRECISION,slave,tag,MPI_COMM_WORLD,status,rerr)
        zetasum=zetasum+partialsum 
        print 250, slave, partialsum,zetasum
        enddo
250   FORMAT('slave = ',I3,' partial sum = ',F18.15 ,' zetasum = ',F18.15)              
      goto 100   
** slave node code
300   call mpi_recv(indices,2,MPI_INTEGER,master,tag,MPI_COMM_WORLD,status,rerr)  
      partialsum = 0
      do n = indices(1), indices(2)
         partialsum = partialsum + 1/n**s
         enddo   
      call mpi_send(partialsum,1,MPI_DOUBLE_PRECISION,master,tag,MPI_COMM_WORLD,serr)     
      goto 100            
400   end