Serial

The first example is from the SEPlib program Nmo3d. In SEPlib there is a module that helps initialize the parallel environment. Our only ifdef will be in whether to include the MPI version of the parallel environment or a dummy serial version. We also need to initialize our parallel environment. The optional arguments signify the input and output to the program. If provided to the sep_par_start routine, intag will be set to "in" for a serial run and "intag" for a parallel run. The variable outtag will be set to "out" or "outtag" and will be initialized with a copy of the intag history file.

#ifdef SEP_MPI
  use sep_par_mpi
#else
  use sep_par
#endif
  call sep_par_start(intag=intag,outtag=outtag)

Initializing the output dataset is unaffected. If we are running in a parallel environment and the input is sectioned, the output tag will be sectioned. Otherwise a standard SEPlib file will be created.

  call init_sep3d(in,out,"OUTPUT",ctag=outtag)  !Initialize SEP

It is often useful to know the thread number for the current process. The below code section bases its decision on how often to print progress on whether it is the master or slave process.

  if(sep_thread_num()==0) then
    call from_param("pct_master",pct_print,2.)
  else
    call from_param("pct_slave",pct_print,10.)
  end if

For NMO our data can be gridded in an arbitrary style. We simply need to loop through all of the data. The init_loop_calc function figures out looping parameters based on the dimensionality of the data in%ndims, the dimensions of the data in%n, and the maximum amount of data we want to read ntr*in%n(1)*in%n(2). The identifier "MAIN" just gives a name for the current loop. The do_sep_loop will increment windowing parameters fwind and nwind until we have looped through entire dataset, at which point the return value will change.

  if(0/=init_loop_calc(in%ndims,in%n,"MAIN",ntr*in%n(1)*in%n(2))) &
     call seperr("trouble initing loop")
   do while(0==do_sep_loop("MAIN",nwind,fwind))

To read the data (whether it is a regular or irregular cube) we first request the number of traces that fall in our current window. In this case, the library will return nh, the number of traces that the current thread locally owns in the given window.

   call sep3d_grab_headers(intag,in,nh,fwind=fwind(2:),nwind=nwind(2:))

If the thread owns any of data in the window, read it.

   if (nh.ne.0) then
     if (.not. sep3d_read_data(intag,in,input(:,:nh))) then
          call erexit("trouble reading data")
        end if

For a regular dataset, which elements of the current window are owned by the current thread is important. We can request to know the coordinates of the traces. The returned 2-D array is of dimensions (ndim-1 by nh), where ndim is the dimensions of the data and nh is the number of traces in the current block.

 if(.not. sep3d_grab_coords(in,coords)) &
           call seperr("trouble grabbing coords")

Writing is again automatic. We will automatically write just the local portion of our data.

 if(.not. sep3d_write_data(outtag,out,output,nwind=nwind,fwind=fwind)) then
    call seperr("trouble writing out ")
 end if

The only remaining step we need to do is to make sure that the total number of traces is correctly calculated and make any calls needed to end the current parallel environment. The library keeps track of how many traces it has written to a given file. The sep3d_update_ntraces call tells the master thread how many traces each slave thread has written. The total number is then written to the output history (and possibly header format file).

 if(.not. sep3d_update_ntraces(out)) call seperr("trouble updating ntraces")
   call sep3d_rite_num_traces(outtag,out)
 call sep_par_stop()