#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 ifFor 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 ifFor 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 ifThe 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()