ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/nano_mpi/src/xyz_io_module.F90
Revision: 4
Committed: Mon Jun 10 17:18:36 2002 UTC (22 years, 1 month ago) by chuckv
File size: 25818 byte(s)
Log Message:
Import Root

File Contents

# User Rev Content
1 chuckv 4 module xyz_io
2     use definitions, ONLY : DP,machdep_lnblnk,machdep_flush
3     use parameter
4     use simulation
5     use model_module
6     use status
7     use file_units, ONLY : next_unit
8     #ifdef MPI
9     use mpi_module
10     #endif
11     implicit none
12     PRIVATE
13    
14    
15     !various common formats
16     character(len=*), parameter :: xyz_line1 = "(i11)"
17     character(len=*), parameter :: xyz_line2 = "(i11,3es12.3)"
18     character(len=*), parameter :: xyz_coord = "(a3,6es12.3,i3)"
19    
20    
21    
22     ! atype interface
23     character(len=2) :: charident
24     character(len=2) :: atid
25    
26     integer :: dump_unit
27     integer :: eor_unit
28     integer :: force_unit
29    
30     #ifdef MPI
31     integer, dimension(mpi_status_size) :: istatus
32     #endif
33    
34     public :: open_xyz_files
35     public :: read_initial_config
36     public :: dump_config
37     public :: write_final_config
38     public :: dump_forces
39     #ifdef MPI
40     public :: mpi_dump_q
41     #endif
42    
43    
44     contains
45    
46     !! open input and output files for simulation
47     subroutine open_xyz_files(this_error)
48     character(len=80) :: dump_file, msg
49     character(len=80) :: force_file
50     character(len=*), parameter :: routine_name = "OPEN_FILES"
51     integer, optional :: this_error
52     integer :: err
53     #ifndef MPI
54     integer, parameter :: node = 0
55     #endif
56    
57     if (present(this_error)) this_error = 0
58    
59     if (node == 0) then
60    
61     dump_unit = next_unit()
62    
63     dump_file = file_prefix(1:machdep_lnblnk(file_prefix)) // '.dump'
64     open(file=dump_file, unit=dump_unit, status='replace', iostat=err, &
65     action='write')
66     if (err > 0) then
67     write(msg,*) "Error opening dump file", trim(dump_file)
68     call info(routine_name, msg)
69     if (present(this_error)) this_error = -1
70     end if
71     write(msg,*) 'configurations will dump to: ', trim(dump_file)
72     call info(routine_name, msg)
73    
74    
75    
76     !!$ force_file = file_prefix(1:machdep_lnblnk(file_prefix)) // '.force'
77     !!$ open(file=force_file, unit=force_unit, status='replace', iostat= err, &
78     !!$ action='write')
79     !!$
80     !!$ if (err > 0) then
81     !!$ write(msg,*) "Error opening force file", trim(force_file)
82     !!$ call info(routine_name, msg)
83     !!$ end if
84     !!$ write(msg,*) 'Force dump will be in: ', trim(force_file)
85     !!$ call info(routine_name, msg)
86    
87     end if
88     end subroutine open_xyz_files
89    
90    
91    
92     subroutine read_initial_config(this_error)
93     ! Reads In the current configuration in XYZ file format. The next
94     ! frame in XYZ is simply appended to the end of the current frame.
95    
96     !standard io
97     integer :: fu, i, a, atom, nats, id, junk,j
98     integer :: natoms_tmp, id_tmp
99     integer, intent(out), optional :: this_error
100     character(len=2), allocatable, dimension(:) :: atid_tmp, atid_array
101     character(len=120) :: msg
102     integer :: read_error
103     integer :: xyz_input_unit
104     integer :: err
105     character(len=80) :: input_file
106    
107     !mpi stuff
108     integer :: iproc, ncount
109     #ifndef MPI
110     integer,parameter :: node = 0
111     #endif
112    
113     natoms_tmp = 0
114    
115     if (present(this_error)) this_error = 0
116    
117     call info("READ_INITIAL_CONFIG","Reading old trajectory file")
118     #ifdef MPI
119     if (node == 0) then
120     #endif
121     input_file = file_prefix(1:machdep_lnblnk(file_prefix)) // '.in'
122     xyz_input_unit = next_unit()
123     open(file=input_file, unit=xyz_input_unit, status='old',iostat=err, &
124     action='read')
125     if (err > 0 ) then
126     write(msg,*) "Error opening config file", trim(input_file)
127     call info("READ_INITIAL_CONFIG", msg)
128     if (present(this_error)) this_error = -1
129     end if
130    
131     write(msg,*) 'reading initial config from: ', trim(input_file)
132     call info("READ_INITIAL_CONFIG", msg)
133     #ifdef MPI
134     endif
135     #endif
136    
137    
138     #ifndef MPI
139    
140     read(unit=xyz_input_unit,fmt=*,iostat=read_error) natoms_tmp
141     read(xyz_input_unit,*)
142    
143     if (read_error /= 0) then
144     write(msg,*) "Error reading particle configuration: ", i, "Atom: ",atom
145     call info("READ_INITIAL_CONFIG",msg)
146     if (present(this_error)) this_error = -1
147     return
148     else if (natoms /= natoms_tmp) then
149     write(msg,*) 'Atom number mismatch read: ',natoms_tmp, &
150     'Allocated: ', natoms
151     call info('READ_INITIAL_CONFIG', msg)
152     if (present(this_error)) this_error = -1
153     return
154     end if
155    
156    
157    
158     ! this block is only true for monatomic systems
159     !! fix me for general case
160     nmol = natoms
161     do i = 1, nmol
162     na(i) = 1
163     atom_index(i,1) = i
164     enddo
165    
166     do i = 1, nmol
167     do a = 1, na(i)
168     atom = atom_index(i,a)
169     read(unit=xyz_input_unit,fmt=*,iostat=read_error) atid, q(1,atom), q(2,atom), q(3,atom), &
170     v(1,atom), v(2,atom), v(3,atom)
171     if (read_error /= 0) then
172     write(msg,*) "Error reading particle configuration: ", i, "Atom: ",atom
173     call info("READ_INITIAL_CONFIG",msg)
174     if (present(this_error)) this_error = -1
175     return
176     end if
177     id = atype_atoi(atid)
178     ident(atom) = id
179     ! this line only works for monatomic systems:
180     model(i) = atype_identifier(id)
181     mass(atom) = atype_mass(id)
182     end do
183     end do
184    
185     #else
186     if (node == 0) then
187     allocate(atid_tmp(max_local))
188     else
189     allocate(atid_array(max_local))
190     end if
191    
192    
193     if (node == 0) then
194     read(xyz_input_unit,*) natoms_tmp
195     read(xyz_input_unit,*)
196     end if
197    
198     call mpi_bcast(natoms_tmp,1,mpi_integer,0,mpi_comm_world,mpi_err)
199     if (mpi_err /= 0) then
200     call error("READ_INITIAL_CONFIG","Error bcast natoms_tmp")
201     endif
202    
203     if (natoms_tmp /= natoms) then
204     write(msg,*) 'Atom number mismatch read: ',natoms_tmp, &
205     'Allocated: ', natoms
206     call error('READ_INITIAL_CONFIG', msg)
207     end if
208    
209     call mpi_bcast(natoms,1,mpi_integer,0,mpi_comm_world,mpi_err)
210     if (mpi_err /= 0) then
211     call error("READ_INITIAL_CONFIG","Error bcast natoms")
212     endif
213    
214     ! this block is only true for monatomic systems
215     nmol = natoms
216    
217     do i = 1, maxmol_local
218     na(i) = 1
219     atom_index(i,1) = i
220     end do
221    
222     if (node == 0) then
223     do iproc = 0, nprocs - 1
224    
225     ncount = nint(float(iproc+1)/nprocs*nmol) - &
226     nint(float(iproc)/nprocs*nmol)
227     do i = 1,ncount
228     do a = 1, na(i)
229     atom = atom_index(i,a)
230     read(xyz_input_unit,*) atid_tmp(i), q_tmp(1,atom), q_tmp(2,atom), &
231     q_tmp(3,atom), &
232     v_tmp(1,atom), v_tmp(2,atom), v_tmp(3,atom)
233     end do
234     end do
235    
236     if (iproc.eq.0) then
237     do i = 1,ncount
238     do a = 1, na(i)
239     atom = atom_index(i,a)
240     q(1,atom) = q_tmp(1,atom)
241     q(2,atom) = q_tmp(2,atom)
242     q(3,atom) = q_tmp(3,atom)
243     v(1,atom) = v_tmp(1,atom)
244     v(2,atom) = v_tmp(2,atom)
245     v(3,atom) = v_tmp(3,atom)
246     id = atype_atoi(atid_tmp(i))
247     ident(atom) = id
248     ! this line only works for monatomic systems:
249     model(i) = atype_identifier(id)
250     mass(atom) = atype_mass(id)
251     end do
252     enddo
253     else
254    
255     call mpi_send(atid_tmp,2*ncount,mpi_character,iproc,0, &
256     mpi_comm_world,mpi_err)
257     if (mpi_err /= 0) then
258     call error("READ_INITIAL_CONFIG","Error mpi_send atid_tmp")
259     endif
260     call mpi_send(q_tmp,3*ncount,mpi_double_precision, &
261     iproc,0,mpi_comm_world,mpi_err)
262     if (mpi_err /= 0) then
263     call error("READ_INITIAL_CONFIG","Error mpi_send q_tmp")
264     endif
265     call mpi_send(v_tmp,3*ncount,mpi_double_precision, &
266     iproc,0,mpi_comm_world,mpi_err)
267     if (mpi_err /= 0) then
268     call error("READ_INITIAL_CONFIG","Error mpi_send v_tmp")
269     endif
270    
271     endif
272    
273     end do
274     close(xyz_input_unit)
275     else
276     call mpi_recv(atid_array,2*max_local,mpi_character,0,0, &
277     mpi_comm_world,istatus,mpi_err)
278     if (mpi_err /= 0) then
279     call error("READ_INITIAL_CONFIG","Error reading atid_array")
280     endif
281     call mpi_recv(q,3*max_local, &
282     mpi_double_precision,0,0,mpi_comm_world,istatus,mpi_err)
283     if (mpi_err /= 0) then
284     call error("READ_INITIAL_CONFIG","Error reading q_array")
285     endif
286     call mpi_recv(v,3*max_local, &
287     mpi_double_precision,0,0,mpi_comm_world,istatus,mpi_err)
288     if (mpi_err /= 0) then
289     call error("READ_INITIAL_CONFIG","Error reading v_array")
290     endif
291    
292     do j = 1,nlocal
293     ! array error here.....
294     id = atype_atoi(atid_array(j))
295     ident(j) = id
296     ! this line only works for monatomic systems:
297     model(j) = atype_identifier(id)
298     mass(j) = atype_mass(id)
299     end do
300    
301     end if
302     call mpi_barrier(mpi_comm_world,mpi_err)
303    
304    
305    
306     if (node == 0) then
307     deallocate(atid_tmp)
308     close(xyz_input_unit)
309     else
310     deallocate(atid_array)
311     end if
312     #endif
313     end subroutine read_initial_config
314    
315     subroutine dump_config(time)
316     include 'headers/atom.h'
317    
318     real(kind=DP), intent(in) :: time
319     #ifdef MPI
320     integer, dimension(maxmol_local) :: na_tmp
321     integer, dimension(maxmol_local,3) :: atom_index_tmp
322     integer, dimension(max_local) :: ident_tmp
323     #endif
324     integer ::i,j, iproc, a
325     integer :: atom
326     integer :: icnt, i_a_cnt
327     character(len=10) :: chartmp
328    
329     #ifndef MPI
330     write(dump_unit,*) natoms
331     write(dump_unit,*) time
332    
333     do i = 1, natoms
334     charident = atype_identifier(ident(i)) // ' '
335     write(unit=dump_unit,fmt=xyz_coord) charident, (q(j,i),j=1,3), &
336     (v(j,i),j=1,3), ident(i)
337     end do
338    
339     call machdep_flush(dump_unit)
340     #else
341     if (node /= 0) then
342     call mpi_send(na,nmol_local,mpi_integer,0,0, &
343     mpi_comm_world,mpi_err)
344     if (mpi_err /= 0) then
345     call error("DUMP_CONFIG","Error sending na")
346     endif
347     call mpi_send(atom_index,3*nmol_local,mpi_integer,0,0, &
348     mpi_comm_world,mpi_err)
349     if (mpi_err /= 0) then
350     call error("DUMP_CONFIG","Error sending atom_index")
351     endif
352     call mpi_send(ident,nlocal,mpi_integer,0,0, &
353     mpi_comm_world,mpi_err)
354     if (mpi_err /= 0) then
355     call error("DUMP_CONFIG","Error sending ident")
356     endif
357     call mpi_send(q,3*nlocal,mpi_double_precision,0,0, &
358     mpi_comm_world,mpi_err)
359     if (mpi_err /= 0) then
360     call error("DUMP_CONFIG","Error sending q")
361     endif
362     call mpi_send(v,3*nlocal,mpi_double_precision,0,0, &
363     mpi_comm_world,mpi_err)
364     if (mpi_err /= 0) then
365     call error("DUMP_CONFIG","Error sending v")
366     endif
367    
368     else
369    
370    
371     write(dump_unit,*) natoms
372     write(dump_unit,*) time
373    
374    
375     do i = 1, nlocal
376     charident = atype_identifier(ident(i)) // ' '
377     write(unit=dump_unit,fmt=xyz_coord) charident, (q(j,i),j=1,3), &
378     (v(j,i),j=1,3), ident(i)
379     end do
380    
381     do iproc = 1,nprocs - 1
382     icnt = 0
383     call mpi_recv(na_tmp, maxmol_local, mpi_integer,iproc,0, &
384     mpi_comm_world,istatus,mpi_err)
385    
386     call mpi_get_count(istatus,mpi_integer,icnt,mpi_err)
387    
388     call mpi_recv(atom_index_tmp,3*icnt,mpi_integer,iproc,0, &
389     mpi_comm_world,istatus,mpi_err)
390     if (mpi_err /= 0) then
391     call error("DUMP_CONFIG","Error recieving atom_index")
392     endif
393    
394     call mpi_recv(ident_tmp,max_local,mpi_integer,iproc,0, &
395     mpi_comm_world,istatus,mpi_err)
396     if (mpi_err /= 0) then
397     call error("DUMP_CONFIG","Error recieving ident_tmp")
398     endif
399    
400     call mpi_get_count(istatus,mpi_integer,i_a_cnt,mpi_err)
401    
402     call mpi_recv(q_tmp,3*i_a_cnt,mpi_double_precision,iproc,0, &
403     mpi_comm_world,istatus,mpi_err)
404     if (mpi_err /= 0) then
405     call error("DUMP_CONFIG","Error recieving q_tmp")
406     endif
407    
408     call mpi_recv(v_tmp,3*i_a_cnt,mpi_double_precision,iproc,0, &
409     mpi_comm_world,istatus,mpi_err)
410     if (mpi_err /= 0) then
411     call error("DUMP_CONFIG","Error recieving v_tmp")
412     endif
413    
414    
415     do i = 1,icnt
416     do a = 1, na_tmp(i)
417     atom = atom_index_tmp(i,a)
418     atid = atype_identifier(ident_tmp(atom)) // ' '
419     write(dump_unit,fmt=xyz_coord) atid, q_tmp(1,atom), &
420     q_tmp(2,atom), q_tmp(3,atom), v_tmp(1,atom), &
421     v_tmp(2,atom), v_tmp(3,atom), ident_tmp(atom)
422     end do
423     enddo
424    
425     end do
426     endif
427    
428     call mpi_barrier(mpi_comm_world,mpi_err)
429    
430     if (node == 0) &
431     call machdep_flush(dump_unit)
432     #endif
433     end subroutine dump_config
434    
435    
436     subroutine write_final_config()
437     ! Reads In the current configuration in XYZ file format. The next
438     ! frame in XYZ is simply appended to the end of the current frame.
439    
440     include 'headers/atom.h'
441    
442     character(len = 80) :: eor_file
443     character(len = 120) :: msg
444     integer err
445    
446     #ifdef MPI
447     integer, dimension(maxmol_local) :: na_tmp
448     integer, dimension(maxmol_local,3) :: atom_index_tmp
449     integer, dimension(max_local) :: ident_tmp
450     #endif
451     integer ::i,j, iproc, a
452     integer :: atom
453     integer :: icnt, i_a_cnt
454    
455     character(len=10) :: chartmp
456    
457    
458    
459    
460     #ifdef MPI
461     if (node == 0) then
462     #endif
463     eor_unit = next_unit()
464    
465     eor_file = file_prefix(1:machdep_lnblnk(file_prefix)) // '.eor'
466     open(file=eor_file, unit=eor_unit, status='replace', iostat= err, &
467     action='write')
468    
469     if (err > 0) then
470     write(msg,*) "Error opening eor file", trim(eor_file)
471     call info("DUMP_FINAL_CONFIG", msg)
472     end if
473     write(msg,*) 'final configuration will be in: ', trim(eor_file)
474     call info("DUMP_FINAL_CONFIG", msg)
475    
476     #ifdef MPI
477     endif
478     #endif
479    
480    
481    
482    
483    
484    
485    
486    
487     #ifndef MPI
488    
489    
490     write(eor_unit,fmt=xyz_line1) natoms
491     write(eor_unit,*)
492    
493     do i = 1, natoms
494     do a = 1, na(i)
495     atom = atom_index(i,a)
496     atid = atype_identifier(ident(atom)) // ' '
497     write(eor_unit,fmt=xyz_coord) atid, q(1,atom), q(2,atom), q(3,atom), &
498     v(1,atom), v(2,atom), v(3,atom), ident(i)
499     end do
500     end do
501    
502     close(eor_unit)
503     #else
504    
505     if (node /= 0) then
506    
507     call mpi_send(na,nmol_local,mpi_integer,0,0, &
508     mpi_comm_world,mpi_err)
509     if (mpi_err /= 0) then
510     call error("WRITE_FINALCONFIG","Error sending na")
511     endif
512     call mpi_send(atom_index,3*nmol_local,mpi_integer,0,0, &
513     mpi_comm_world,mpi_err)
514     if (mpi_err /= 0) then
515     call error("WRITE_FINALCONFIG","Error sending atom_index")
516     endif
517     call mpi_send(ident,nlocal,mpi_integer,0,0, &
518     mpi_comm_world,mpi_err)
519     if (mpi_err /= 0) then
520     call error("WRITE_FINALCONFIG","Error sending ident")
521     endif
522     call mpi_send(q,3*nlocal,mpi_double_precision,0,0, &
523     mpi_comm_world,mpi_err)
524     if (mpi_err /= 0) then
525     call error("WRITE_FINALCONFIG","Error sending q")
526     endif
527     call mpi_send(v,3*nlocal,mpi_double_precision,0,0, &
528     mpi_comm_world,mpi_err)
529     if (mpi_err /= 0) then
530     call error("WRITE_FINALCONFIG","Error sending v")
531     endif
532    
533     else
534    
535    
536     write(eor_unit,*) natoms
537     write(eor_unit,*)
538    
539    
540     do i = 1, nlocal
541     charident = atype_identifier(ident(i)) // ' '
542     write(unit=eor_unit,fmt=xyz_coord) charident, (q(j,i),j=1,3), &
543     (v(j,i),j=1,3), ident(i)
544     end do
545    
546     do iproc = 1,nprocs - 1
547     icnt = 0
548    
549     call mpi_recv(na_tmp, maxmol_local, mpi_integer,iproc,0, &
550     mpi_comm_world,istatus,mpi_err)
551     if (mpi_err /= 0) then
552     call error("WRITE_FINALCONFIG","Error receiving na_tmp")
553     endif
554    
555     call mpi_get_count(istatus,mpi_integer,icnt,mpi_err)
556    
557     call mpi_recv(atom_index_tmp,3*icnt,mpi_integer,iproc,0, &
558     mpi_comm_world,istatus,mpi_err)
559     if (mpi_err /= 0) then
560     call error("WRITE_FINALCONFIG","Error receiving atom_index_tmp")
561     endif
562    
563     call mpi_recv(ident_tmp,max_local,mpi_integer,iproc,0, &
564     mpi_comm_world,istatus,mpi_err)
565     if (mpi_err /= 0) then
566     call error("WRITE_FINAL_CONFIG","Error receiving ident_tmp")
567     endif
568    
569     call mpi_get_count(istatus,mpi_integer,i_a_cnt,mpi_err)
570     call mpi_recv(q_tmp,3*i_a_cnt,mpi_double_precision,iproc,0, &
571     mpi_comm_world,istatus,mpi_err)
572     if (mpi_err /= 0) then
573     call error("WRITE_FINAL_CONFIG","Error receiving q_tmp")
574     endif
575     call mpi_recv(v_tmp,3*i_a_cnt,mpi_double_precision,iproc,0, &
576     mpi_comm_world,istatus,mpi_err)
577     if (mpi_err /= 0) then
578     call error("WRITE_FINAL_CONFIG","Error receiving v_tmp")
579     endif
580    
581    
582     do i = 1,icnt
583     do a = 1, na_tmp(i)
584     atom = atom_index_tmp(i,a)
585     atid = atype_identifier(ident_tmp(atom)) // ' '
586     write(eor_unit,fmt=xyz_coord) atid, q_tmp(1,atom), &
587     q_tmp(2,atom), q_tmp(3,atom), v_tmp(1,atom), &
588     v_tmp(2,atom), v_tmp(3,atom), ident_tmp(atom)
589     end do
590     enddo
591    
592     end do
593     endif
594    
595     call mpi_barrier(mpi_comm_world,mpi_err)
596    
597     if (node == 0) then
598     call machdep_flush(eor_unit)
599     close(eor_unit)
600     end if
601     #endif
602     end subroutine write_final_config
603    
604    
605     subroutine dump_forces(this_time)
606     ! Reads In the current configuration in XYZ file format. The next
607     ! frame in XYZ is simply appended to the end of the current frame.
608    
609     ! include 'headers/fileio.h'
610     ! include 'headers/atom.h'
611    
612     real( kind=DP ), intent(in) :: this_time
613     integer ::i,j, iproc, a
614     #ifdef MPI
615     integer, dimension(max_local) :: tag_tmp
616     #endif
617     integer :: icnt
618     character(Len=*), parameter :: forcefmt = '(I3,3es18.8)'
619    
620     #ifdef MPI
621     if (node /= 0) then
622    
623     call mpi_send(tag_local,nlocal,mpi_integer,0,0, &
624     mpi_comm_world,mpi_err)
625    
626     call mpi_send(f,3*nlocal,mpi_double_precision,0,0, &
627     mpi_comm_world,mpi_err)
628    
629     call mpi_send(rho,nlocal,mpi_double_precision,0,0, &
630     mpi_comm_world,mpi_err)
631    
632     ! call mpi_send(frho,nlocal,mpi_double_precision,0,0, &
633     ! mpi_comm_world,mpi_err)
634    
635     else
636    
637     write(force_unit,'(es12.4)') this_time
638    
639    
640     do i = 1, nlocal
641     write(unit=force_unit,fmt=*) tag_local(i),rho(i),f(1,i), &
642     f(2,i),f(3,i)
643     end do
644    
645     do iproc = 1,nprocs - 1
646     icnt = 0
647     f_tmp = 0.0d0
648    
649    
650     call mpi_recv(tag_tmp,max_local,mpi_integer,iproc,0, &
651     mpi_comm_world,istatus,mpi_err)
652     call mpi_get_count(istatus,mpi_integer,icnt,mpi_err)
653     call mpi_recv(f_tmp,3*icnt,mpi_double_precision,iproc,0, &
654     mpi_comm_world,istatus,mpi_err)
655     call mpi_recv(rho_tmp,icnt,mpi_double_precision,iproc,0, &
656     mpi_comm_world,istatus,mpi_err)
657     ! call mpi_recv(frho_tmp,icnt,mpi_double_precision,iproc,0, &
658     ! mpi_comm_world,istatus,mpi_err)
659    
660    
661     do i = 1,icnt
662     write(unit=force_unit,fmt=*) tag_tmp(i),rho_tmp(i), &
663     f_tmp(1,i),f_tmp(2,i),f_tmp(3,i)
664     enddo
665    
666     end do
667     endif
668    
669     call mpi_barrier(mpi_comm_world,mpi_err)
670    
671     if (node == 0) then
672     call machdep_flush(dump_unit)
673     end if
674     #else
675     write(force_unit,'(es12.4)') this_time
676     do i = 1, natoms
677     write(unit=force_unit,fmt=*) i,rho(i),f(1,i), &
678     f(2,i),f(3,i)
679     end do
680    
681     #endif
682    
683     end subroutine dump_forces
684    
685    
686     !! remove me
687     #ifdef MPI
688     subroutine mpi_dump_q()
689     integer, parameter :: q_row_unit = 110
690     integer, parameter :: q_col_unit = 111
691     integer :: processor
692     integer :: i, err
693     integer :: nrow_count,ncol_count
694     character(len=60) :: q_row_file,q_col_file,q_file
695     character(len=5) :: char_node
696     real(kind=8), dimension(3,max_row) :: q_row_tmp
697     real(kind=8), dimension(3,max_col) :: q_col_tmp
698     integer, dimension(max_row) :: tag_row_tmp,ident_row_tmp
699     integer, dimension(max_col) :: tag_col_tmp,ident_col_tmp
700     integer :: nrow_tmp, ncol_tmp, nlocal_tmp
701    
702     if (node /= 0 ) then
703     call mpi_send(nlocal,1,mpi_integer,0,0,mpi_comm_world,mpi_err)
704     call mpi_send(nrow,1,mpi_integer,0,0,mpi_comm_world,mpi_err)
705     call mpi_send(ncol,1,mpi_integer,0,0,mpi_comm_world,mpi_err)
706     call mpi_send(ident_row,nrow,mpi_integer,0,0, &
707     mpi_comm_world,mpi_err)
708     call mpi_send(ident_col,ncol,mpi_integer,0,0, &
709     mpi_comm_world,mpi_err)
710     call mpi_send(q_row,3*nrow,mpi_double_precision,0,0, &
711     mpi_comm_world,mpi_err)
712     call mpi_send(q_col,3*ncol,mpi_double_precision,0,0, &
713     mpi_comm_world,mpi_err)
714     call mpi_send(tag_row,nrow,mpi_integer,0,0, &
715     mpi_comm_world,mpi_err)
716     call mpi_send(tag_col,ncol,mpi_integer,0,0, &
717     mpi_comm_world,mpi_err)
718     else
719    
720     !! node 0 writes its arrays first
721     write(char_node,"(a1,i1,a1)") "_",node,"_"
722     char_node = trim(char_node)
723     q_row_file = file_prefix(1:machdep_lnblnk(file_prefix)) // &
724     char_node(1:machdep_lnblnk(char_node)) // &
725     '.row'
726     open(file=q_row_file, unit=q_row_unit, status='replace', iostat= err, &
727     action='write')
728     write(unit=q_row_unit,fmt="(i5,i5)") node,nrow
729     do i = 1, nrow
730     write(unit=q_row_unit,fmt="(i5,3es12.4,i4)") tag_row(i),q_row(1,i), &
731     q_row(2,i), q_row(3,i), ident_row(i)
732     end do
733    
734     close(unit=q_row_unit)
735    
736     q_col_file = file_prefix(1:machdep_lnblnk(file_prefix)) // &
737     char_node(1:machdep_lnblnk(char_node)) // &
738     '.col'
739     open(file=q_col_file, unit=q_col_unit, status='replace', iostat= err, &
740     action='write')
741     write(unit=q_col_unit,fmt="(i5,i5)") node,ncol
742     do i = 1, ncol
743     write(unit=q_col_unit,fmt="(i5,3es12.4,i4)") tag_col(i),q_col(1,i), q_col(2,i), &
744     q_col(3,i), ident_col(i)
745     end do
746    
747     close(unit=q_col_unit)
748     write(*,*) node, nlocal, nrow, ncol
749     do processor = 1, nprocs - 1
750     !! recv data from other procs
751     call mpi_recv(nlocal_tmp,1,mpi_integer,processor,0,mpi_comm_world,&
752     istatus,mpi_err)
753     call mpi_recv(nrow_tmp,1,mpi_integer,processor,0,mpi_comm_world,&
754     istatus,mpi_err)
755     call mpi_recv(ncol_tmp,1,mpi_integer,processor,0,mpi_comm_world,&
756     istatus,mpi_err)
757     write(*,*) processor, nlocal_tmp, nrow_tmp, ncol_tmp
758     call mpi_recv(ident_row_tmp,nrow_tmp,mpi_integer,processor,0, &
759     mpi_comm_world,istatus,mpi_err)
760     write(*,*) "Past ident row"
761     ! call mpi_get_count(istatus,mpi_integer, nrow_count,mpi_err)
762     call mpi_recv(ident_col_tmp,ncol_tmp,mpi_integer,processor,0, &
763     mpi_comm_world,istatus,mpi_err)
764     ! call mpi_get_count(istatus,mpi_integer,ncol_count,mpi_err)
765     write(*,*) "Past ident_col "
766     call mpi_recv(q_row_tmp,3*nrow_tmp,mpi_double_precision,processor,0,&
767     mpi_comm_world,istatus,mpi_err)
768     write(*,*) "Past q_row"
769     call mpi_recv(q_col_tmp,3*ncol_tmp,mpi_double_precision,processor,0,&
770     mpi_comm_world,istatus,mpi_err)
771     write(*,*) "Past q_col"
772     call mpi_recv(tag_row_tmp,nrow_tmp,mpi_integer,processor,0, &
773     mpi_comm_world,istatus,mpi_err)
774     write(*,*) "Past tag row"
775     call mpi_recv(tag_col_tmp,ncol_tmp,mpi_integer,processor,0, &
776     mpi_comm_world,istatus,mpi_err)
777     write(*,*) "Past tag col"
778     write(char_node,"(a1,i1,a1)") "_",processor,"_"
779     char_node = trim(char_node)
780     q_row_file = file_prefix(1:machdep_lnblnk(file_prefix)) // &
781     char_node(1:machdep_lnblnk(char_node)) // '.row'
782     open(file=q_row_file, unit=q_row_unit, status='replace', iostat= err, &
783     action='write')
784     write(unit=q_row_unit,fmt="(i5,i5)") node,nrow
785     do i = 1, nrow_tmp
786     write(unit=q_row_unit,fmt="(i5,3es12.4,i4)") tag_row_tmp(i),q_row_tmp(1,i), &
787     q_row_tmp(2,i),q_row_tmp(3,i), ident_row_tmp(i)
788     end do
789    
790     close(unit=q_row_unit)
791    
792     q_col_file = file_prefix(1:machdep_lnblnk(file_prefix)) // &
793     char_node(1:machdep_lnblnk(char_node)) // &
794     '.col'
795     open(file=q_col_file, unit=q_col_unit, status='replace', iostat= err, &
796     action='write')
797     write(unit=q_col_unit,fmt="(i5,i5)") node,ncol
798     do i = 1, ncol_tmp
799     write(unit=q_col_unit,fmt="(i5,3es12.4,i4)") tag_col_tmp(i),q_col_tmp(1,i), &
800     q_col_tmp(2,i), q_col_tmp(3,i), ident_col_tmp(i)
801     end do
802    
803     close(unit=q_col_unit)
804     enddo
805    
806     endif
807     call mpi_barrier(mpi_comm_world,mpi_err)
808     end subroutine mpi_dump_q
809    
810     #endif
811    
812     end module xyz_io