1      subroutine sp_start(lout,lfntop,filtop,lfnrst,filrst,
2     + lsyn,filsyn,nsyn,
3     + rsht,rlng,rinp,rsgmi,
4     + npxi,npyi,npzi,nbxi,nbyi,nbzi,
5     + npbt,nbxt,boxt,vlatt,lpbc,
6     + nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai,
7     + ldb,lbp,fld,lpol,lfre,temp,tempw,temps,lqd,iguidi,
8     + ldbg,idbg,prjct,mbbi,nseri,isld,irset,ictrl,nseqi,
9     + i_lseqi,ndumsi,nbgeti,npreci,madbxi)
10c
11c
12c $Id$
13c
14      implicit none
15c
16#include "sp_common.fh"
17#include "mafdecls.fh"
18#include "msgids.fh"
19#include "global.fh"
20c
21      logical sp_rdrest,sp_diffbb
22      external sp_rdrest,sp_diffbb
23c
24      integer lfntop,lfnrst,lsyn,nsyn,ldbg,idbg,lout
25      character*255 filtop,filrst,filsyn
26      real*8 rsht,rlng,rinp,boxt(3),vlatt(3,3),rtemp(3)
27      real*8 temp,tempw,temps,fld,rsgmi
28      integer itemp(2)
29      integer npxi,npyi,npzi,nbxi,nbyi,nbzi,isld,irset,nseqi,i_lseqi
30      integer nwmi,mwmi,nwai,mwai,nsfi,msfi,nsmi,msmi,nsai,msai
31      integer lenscr,ldb,lbp,npbt,nbxt,iguidi,mbbi,nseri,ictrl,ndumsi
32      integer nbgeti,npreci,madbxi
33      logical lpol,lfre,ignore,lpbc,lqd
34      character*80 prjct
35c
36      integer i,j,l_f,i_f
37c
38      if(.not.ma_verify_allocator_stuff()) then
39      call md_abort('ERROR IN MEMORY USE AT START SP_START',0)
40      endif
41c
42      me=ga_nodeid()
43      np=ga_nnodes()
44c
45      project=prjct
46c
47      idebug=idbg
48      lfndbg=ldbg
49      icntrl=ictrl
50      lfnout=lout
51c
52      lfnsyn=lsyn
53      nfsync=nsyn
54      isload=isld
55      ireset=irset
56      nbget=nbgeti
57      nprec=npreci
58      ibget=nbget
59      madbox=madbxi
60c
61      rshort=rsht
62      rlong=max(rsht,rlng)
63      rbox=rinp
64      ntwin=1
65      if(rlng.gt.rsht) ntwin=2
66      ltwin=rlng.gt.rsht
67      lnode0=lqmd
68      lqmd=lqd
69c
70      loadb=ldb
71      lbpair=lbp
72      factld=fld
73c
74      lpola=lpol
75      lfree=lfre
76c
77      nbx=nbxi
78      nby=nbyi
79      nbz=nbzi
80      nbxin=nbxi
81      nbyin=nbyi
82      nbzin=nbzi
83c
84      npx=npxi
85      npy=npyi
86      npz=npzi
87c
88      mbbl=0
89      nbbdif=-1
90c
91      mwm=mwmi
92      msa=max(msai,msmi)
93      mbbreq=mbbi
94      mbblp=0
95      nserie=nseri
96c
97      iguide=iguidi
98c
99      npack=0
100      npackw=0
101c
102      rsgm=rsgmi
103c
104      call sp_nrnode()
105c
106      call sp_dimens(lfnrst,filrst)
107c
108      call sp_alloc()
109c
110      call sp_decomp(int_mb(i_iown),dbl_mb(i_boxs),
111     + int_mb(i_buren),int_mb(i_bindex))
112c
113      call sp_initip(int_mb(i_iown),int_mb(i_ipl))
114c
115      call sp_numbb(int_mb(i_iown),dbl_mb(i_boxs))
116c
117      if(sp_diffbb(dbl_mb(i_boxs),int_mb(i_rng)))
118     + call sp_listbb(int_mb(i_iown),dbl_mb(i_boxs),int_mb(i_bb))
119c
120      call sp_alloc2()
121c
122      if(ireset.eq.0) then
123      ignore=sp_rdrest(lfnrst,filrst,dbl_mb(i_boxs))
124      endif
125      if(np.gt.1) then
126      call ga_brdcst(msp_02,nsm,ma_sizeof(mt_int,1,mt_byte),0)
127      endif
128      msm=max(1,nsm)
129      if(.not.ma_push_get(mt_dbl,msm*3,'xscr',l_xscr,i_xscr))
130     + call md_abort('Failed to allocate xscr',0)
131c
132      lenscr=3*max(mwm*mwa,msa)
133      if(.not.ma_push_get(mt_dbl,lenscr,'x',l_x,i_x))
134     + call md_abort('Failed to allocate x',0)
135      if(.not.ma_push_get(mt_dbl,lenscr,'v',l_v,i_v))
136     + call md_abort('Failed to allocate v',0)
137      if(.not.ma_push_get(mt_dbl,lenscr,'f',l_f,i_f))
138     + call md_abort('Failed to allocate f',0)
139      if(.not.ma_push_get(mt_dbl,3*mwm*mwa,'r',l_r,i_r))
140     + call md_abort('Failed to allocate r',0)
141      lenscr=max(miw2*mwm,mis2*msa)
142      if(.not.ma_push_get(mt_int,lenscr,'i',l_i,i_i))
143     + call md_abort('Failed to allocate i',0)
144      lenscr=ma_inquire_avail(mt_byte)/
145     + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+
146     + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1
147      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx))
148     + call md_abort('Failed to allocate bx',0)
149      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv))
150     + call md_abort('Failed to allocate bv',0)
151      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf))
152     + call md_abort('Failed to allocate bf',0)
153      if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br))
154     + call md_abort('Failed to allocate br',0)
155      if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi))
156     + call md_abort('Failed to allocate bi',0)
157      if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n))
158     + call md_abort('Failed to allocate n',0)
159c
160      call ga_sync()
161c
162      if(.not.ma_verify_allocator_stuff()) then
163      call md_abort('ERROR IN MEMORY USE BEFORE RDRST IN SP_START',0)
164      endif
165c
166      call sp_rdrst(lfnrst,filrst,lfntop,filtop,
167     + temp,tempw,temps,int_mb(i_ipl),
168     + dbl_mb(i_x),dbl_mb(i_v),dbl_mb(i_f),dbl_mb(i_r),int_mb(i_i),
169     + dbl_mb(i_x),dbl_mb(i_v),dbl_mb(i_f),dbl_mb(i_xscr),int_mb(i_i),
170     + dbl_mb(i_bx),dbl_mb(i_bv),dbl_mb(i_bf),dbl_mb(i_br),
171     + int_mb(i_bi),lenscr,
172     + dbl_mb(i_bx),dbl_mb(i_bv),dbl_mb(i_bf),int_mb(i_bi),
173     + lenscr,int_mb(i_n),
174     + int_mb(i_iown),dbl_mb(i_boxs),int_mb(i_lseq),int_mb(i_sndx))
175c
176      if(.not.ma_verify_allocator_stuff()) then
177      call md_abort('ERROR IN MEMORY USE AFTER RDRST IN SP_START',0)
178      endif
179c
180      if(.not.ma_pop_stack(l_n))
181     + call md_abort('Failed to deallocate n',0)
182      if(.not.ma_pop_stack(l_bi))
183     + call md_abort('Failed to deallocate bi',0)
184      if(.not.ma_pop_stack(l_br))
185     + call md_abort('Failed to deallocate br',0)
186      if(.not.ma_pop_stack(l_bf))
187     + call md_abort('Failed to deallocate bf',0)
188      if(.not.ma_pop_stack(l_bv))
189     + call md_abort('Failed to deallocate bv',0)
190      if(.not.ma_pop_stack(l_bx))
191     + call md_abort('Failed to deallocate bx',0)
192      if(.not.ma_pop_stack(l_i))
193     + call md_abort('Failed to deallocate i',0)
194      if(.not.ma_pop_stack(l_r))
195     + call md_abort('Failed to deallocate r',0)
196      if(.not.ma_pop_stack(l_f))
197     + call md_abort('Failed to deallocate f',0)
198      if(.not.ma_pop_stack(l_v))
199     + call md_abort('Failed to deallocate v',0)
200      if(.not.ma_pop_stack(l_x))
201     + call md_abort('Failed to deallocate x',0)
202c
203      if(np.gt.1) then
204      rtemp(1)=temp
205      rtemp(2)=tempw
206      rtemp(3)=temps
207      call ga_brdcst(msp_01,rtemp,ma_sizeof(mt_dbl,3,mt_byte),0)
208      temp=rtemp(1)
209      tempw=rtemp(2)
210      temps=rtemp(3)
211      itemp(1)=nsm
212      itemp(2)=nsf
213      call ga_brdcst(msp_02,itemp,ma_sizeof(mt_int,2,mt_byte),0)
214      nsm=itemp(1)
215      nsf=itemp(2)
216      endif
217c
218      nwmi=nwm
219      nwai=nwa
220      nsmi=nsm
221      nsai=nsa
222c
223      mwmi=mwm
224      mwai=mwa
225      msai=msa
226c
227      npbt=npbtyp
228      nbxt=nbxtyp
229c
230      do 1 i=1,3
231      boxt(i)=box(i)
232      do 2 j=1,3
233      vlatt(i,j)=vlat(i,j)
234    2 continue
235    1 continue
236c
237      nsfi=nsf
238      msfi=max(1,nsf)
239      msmi=max(1,nsm)
240c
241      third=1.0d0/3.0d0
242      nldup=-2
243      ipairf=-1
244      ipairt=-1
245      lpipo=.false.
246c
247      lpbc=npbtyp.gt.0
248c
249      if(me.eq.0) then
250      if(nfsync.gt.0) then
251      open(unit=lfnsyn,file=filsyn(1:index(filsyn,' ')-1),
252     + status='unknown')
253      write(lfnsyn,3000) np
254 3000 format(i5)
255      endif
256      endif
257c
258      nseqi=nseq
259      i_lseqi=i_lseq
260c
261      ndumsi=ndums
262c
263      if(.not.ma_verify_allocator_stuff()) then
264      call md_abort('ERROR IN MEMORY USE AT EXIT SP_START',0)
265      endif
266c
267      return
268      end
269      subroutine sp_dimens(lfnrst,filrst)
270c
271      implicit none
272c
273#include "sp_common.fh"
274#include "msgids.fh"
275#include "mafdecls.fh"
276#include "global.fh"
277#include "geom.fh"
278#include "util.fh"
279c
280      integer lfnrst
281      character*255 filrst
282c
283      character*1 cdum
284      integer i,j,jdum,ibx,iby,ibz,itemp(11)
285      real*8 rtemp(5),rsgmr
286c
287      mbbl=0
288c
289      if(me.eq.0) then
290c
291      open(unit=lfnrst,file=filrst(1:index(filrst,' ')-1),
292     + status='old',form='formatted',err=9999)
293      rewind(lfnrst)
294c
295      do 2 i=1,3
296      read(lfnrst,1001) cdum
297 1001 format(a1)
298    2 continue
299      read(lfnrst,1006) nhist
300 1006 format(32x,i5)
301      do 6 i=1,nhist
302      read(lfnrst,1007) hist(i)
303 1007 format(a)
304    6 continue
305      read(lfnrst,1002) npbtyp,nbxtyp,rsgmr
306 1002 format(2i5,f12.6)
307      if(rsgm.lt.0.0d0) rsgm=rsgmr
308      read(lfnrst,1004) ((vlat(i,j),j=1,3),i=1,3)
309 1004 format(3f12.6)
310      box(1)=vlat(1,1)
311      box(2)=vlat(2,2)
312      box(3)=vlat(3,3)
313      read(lfnrst,1003) jdum
314 1003 format(40x,i5)
315      read(lfnrst,1001) cdum
316      if(jdum.ne.0) then
317      read(lfnrst,1001) cdum
318      endif
319      read(lfnrst,1005) nwm,nwa,nsm,nsa,nwmc,nsf,nseq
320 1005 format(7i10)
321      close(unit=lfnrst,status='keep')
322c
323      bsize=max(rshort+half*rsgm,half*(rlong+half*rsgm),rbox)
324c
325      if(util_print('distribution',print_default)) then
326      if(me.eq.0) write(lfnout,2001)
327 2001 format(/,' Distribution information',/)
328      if(me.eq.0) write(lfnout,2002) rshort,rsgm,rlong,rbox,box,bsize
329 2002 format(' Short range cutoff',t35,f12.6,/,
330     + ' Segment size',t35,f12.6,/,' Long range cutoff',t35,f12.6,/,
331     + ' Box size rbox',t35,f12.6,//,' Box dimension',t35,3f12.6,//,
332     + ' Initial cell size',t35,f12.6,/)
333      endif
334c
335      if(nbx*nby*nbz.lt.np) then
336      nbx=int(box(1)/bsize)
337      nby=int(box(2)/bsize)
338      nbz=int(box(3)/bsize)
339      endif
340c
341      nbx=max(1,nbx,npx)
342      nby=max(1,nby,npy)
343      nbz=max(1,nbz,npz)
344c
345      if(util_print('distribution',print_default)) then
346      if(me.eq.0) then
347      write(lfnout,2003) nbx,nby,nbz
348 2003 format(' Initial cell distribution',t35,3i5)
349      endif
350      endif
351c
352      if(nbxin.eq.0) then
353      nred(1)=0
354      nred(2)=0
355      nred(3)=0
356      if(nbx.gt.npx.and.mod(nbx,npx).gt.0) then
357      nbx=(nbx/npx)*npx
358      nred(1)=nbx
359      endif
360      endif
361      if(nbyin.eq.0) then
362      if(nby.gt.npy.and.mod(nby,npy).gt.0) then
363      nby=(nby/npy)*npy
364      nred(2)=nby
365      endif
366      endif
367      if(nbzin.eq.0) then
368      if(nbz.gt.npz.and.mod(nbz,npz).gt.0) then
369      nbz=(nbz/npz)*npz
370      nred(3)=nbz
371      endif
372      endif
373c
374      if(util_print('distribution',print_default)) then
375      if(me.eq.0) then
376      write(lfnout,2004) nbx,nby,nbz
377 2004 format(' Final cell distribution',t35,3i5,/)
378      endif
379      endif
380c
381      bxmin=bsize/dble(int((dble(nbx)*bsize)/box(1))+1)
382      bymin=bsize/dble(int((dble(nby)*bsize)/box(2))+1)
383      bzmin=bsize/dble(int((dble(nbz)*bsize)/box(3))+1)
384c
385      if(util_print('distribution',print_default)) then
386      if(me.eq.0) then
387      write(lfnout,2005) bxmin,bymin,bzmin
388 2005 format(' Minimum cell size',t35,3f12.6)
389      endif
390      endif
391c
392      endif
393c
394      if(np.gt.1) then
395      itemp(1)=nwm
396      itemp(2)=nwa
397      itemp(3)=nsm
398      itemp(4)=nsa
399      itemp(5)=nbx
400      itemp(6)=nby
401      itemp(7)=nbz
402      itemp(8)=npbtyp
403      itemp(9)=nbxtyp
404      itemp(10)=nsf
405      itemp(11)=nseq
406      rtemp(1)=bsize
407      rtemp(2)=1.001d0*bxmin
408      rtemp(3)=1.001d0*bymin
409      rtemp(4)=1.001d0*bzmin
410      rtemp(5)=rsgm
411      call ga_brdcst(msp_03,itemp,ma_sizeof(mt_int,11,mt_byte),0)
412      call ga_brdcst(msp_04,rtemp,ma_sizeof(mt_dbl,5,mt_byte),0)
413      call ga_brdcst(msp_05,box,ma_sizeof(mt_dbl,3,mt_byte),0)
414      call ga_brdcst(msp_06,vlat,ma_sizeof(mt_dbl,9,mt_byte),0)
415      nwm=itemp(1)
416      nwa=itemp(2)
417      nsm=itemp(3)
418      nsa=itemp(4)
419      nbx=itemp(5)
420      nby=itemp(6)
421      nbz=itemp(7)
422      npbtyp=itemp(8)
423      nbxtyp=itemp(9)
424      nsf=itemp(10)
425      nseq=itemp(11)
426      bsize=rtemp(1)
427      bxmin=rtemp(2)
428      bymin=rtemp(3)
429      bzmin=rtemp(4)
430      rsgm=rtemp(5)
431      endif
432c
433      rbbl=rlong+half*rsgm
434c
435      boxh(1)=half*box(1)
436      boxh(2)=half*box(2)
437      boxh(3)=half*box(3)
438c
439      maxbox=max(nbx,nby,nbz)
440      nbtot=nbx*nby*nbz
441      lpbc0=nbx.eq.1.or.nby.eq.1.or.nbz.eq.1.or.
442     + npx.eq.1.or.npy.eq.1.or.npz.eq.1.or.lpbc9
443c
444      mbox=30
445      do 3 ibx=1,nbx
446      do 4 iby=1,nby
447      do 5 ibz=1,nbz
448      if(me.eq.((((ibz-1)*npz)/nbz)*npy+(((iby-1)*npy)/nby))*npx
449     + +((ibx-1)*npx)/nbx) mbox=mbox+1
450    5 continue
451    4 continue
452    3 continue
453      mbxloc=mbox-30
454c
455      if(np.gt.1) call ga_igop(msp_07,mbox,1,'max')
456c
457      mwa=max(1,nwa)
458      msag=max(1,msa,(mbox-30+madbox)*((nwm*nwa+nsa)/nbtot+1)+1)
459      msa=max(1,msa,(mbox-30+madbox)*((nwm*nwa+nsa)/nbtot+1)+1)
460      mwmg=max(1,msag/mwa+1)
461      mwm=max(1,mwm,msa/mwa+1)
462c
463      msa=min(msa,2*nsa+1)
464      mwm=min(mwm,2*nwm+1)
465      msag=min(msag,2*nsa+1)
466      mwmg=min(mwmg,2*nwm+1)
467c
468      if(lnode0) then
469      msa=nsa+1
470      mwm=nwm+1
471      msag=nsa+1
472      mwmg=nwm+1
473      endif
474c
475      if(util_print('distribution',print_default)) then
476      if(me.eq.0) then
477      write(lfnout,2006) mbox-30,madbox,nbtot
478 2006 format(/,' ARRAY DIMENSION INFORMATION',//,
479     + ' Number cells per processor:  ',i7,/,
480     + ' Number of buffer cells:      ',i7,/,
481     + ' Total number of cells:       ',i7)
482      endif
483      endif
484c
485      return
486 9999 call md_abort('Failed to open restart file',0)
487      return
488      end
489      subroutine sp_alloc
490c
491      implicit none
492c
493#include "sp_common.fh"
494#include "mafdecls.fh"
495#include "global.fh"
496c
497      integer isize
498c
499      if(.not.ma_push_get(mt_int,np,'bindex',l_bindex,i_bindex))
500     + call md_abort('Failed to allocate bindex',0)
501      if(.not.ma_push_get(mt_int,54*np,'buren',l_buren,i_buren))
502     + call md_abort('Failed to allocate buren',0)
503      if(.not.ma_push_get(mt_int,3*maxbox,'owner',l_iown,i_iown))
504     + call md_abort('Failed to allocate owner',0)
505      if(.not.ma_push_get(mt_dbl,3*maxbox,'bxsiz',l_boxs,i_boxs))
506     + call md_abort('Failed to allocate bxsiz',0)
507      if(.not.ma_push_get(mt_int,6*maxbox,'ibxrg',l_boxr,i_boxr))
508     + call md_abort('Failed to allocate ibxrg',0)
509      if(.not.ma_push_get(mt_int,6*maxbox,'rng',l_rng,i_rng))
510     + call md_abort('Failed to allocate rng',me)
511c
512      if(.not.ma_push_get(mt_int,mip2*mbox,'ipl',l_ipl,i_ipl))
513     + call md_abort('Failed to allocate ipl',0)
514      if(.not.ma_push_get(mt_int,mip2*mbox,'jpl',l_jpl,i_jpl))
515     + call md_abort('Failed to allocate jpl',0)
516c
517      mseq=nseq
518c
519      if(.not.ma_push_get(mt_int,mseq,'lseq',l_lseq,i_lseq))
520     + call md_abort('Failed to allocate lseq',0)
521      if(.not.ma_push_get(mt_int,mseq,'sndx',l_sndx,i_sndx))
522     + call md_abort('Failed to allocate sndx',0)
523c
524      if(.not.ga_create(mt_int,np*mbox,mip2,'ip',mbox,mip2,ga_ip))
525     + call md_abort('Failed to create global array ip',0)
526c
527      return
528      end
529      subroutine sp_alloc2
530c
531      implicit none
532c
533#include "sp_common.fh"
534#include "mafdecls.fh"
535#include "global.fh"
536#include "util.fh"
537c
538      integer isize
539c
540      if(nbget.ne.0) then
541      msa=max(1,msa,(mbox-30+nrempr+madbox)*((nwm*nwa+nsa)/nbtot+1)+1)
542      mwm=max(1,mwm,msa/mwa+1)
543      msa=min(msa,2*nsa+1)
544      mwm=min(mwm,2*nwm+1)
545      endif
546c
547      if(util_print('distribution',print_default)) then
548      if(me.eq.0) then
549      if(nbget.ne.0) then
550      write(lfnout,2005) nrempr
551 2005 format(' Number of remote cell pairs: ',i7)
552      if(nbget.gt.0) then
553      write(lfnout,2006) nbget
554 2006 format(' Number of prefetch cells:    ',i7)
555      endif
556      endif
557      write(lfnout,2007) mwm,mwmg
558 2007 format(' Dimension solvent local:     ',i7,', global:',i7)
559      write(lfnout,2008) msa,msag
560 2008 format(' Dimension solute local:      ',i7,', global:',i7)
561      endif
562      endif
563c
564      if(.not.ga_create(mt_int,np*mwmg,miw2,'iw',mwmg,miw2,ga_iw))
565     + call md_abort('Failed to create global array iw',0)
566      isize=6+12*mwa
567      if(lpola) isize=6+18*mwa
568      if(lpola.and.lfree) isize=6+30*mwa
569      if(.not.ga_create(mt_dbl,np*mwmg,isize,'w',mwmg,isize,ga_w))
570     + call md_abort('Failed to create global array w',0)
571      if(.not.ga_create(mt_int,np*msag,mis2,'is',msag,mis2,ga_is))
572     + call md_abort('Failed to create global array is',0)
573      isize=39
574      if(lpola) isize=45
575      if(lpola.and.lfree) isize=57
576      if(.not.ga_create(mt_dbl,np*msag,isize,'s',msag,isize,ga_s))
577     + call md_abort('Failed to create global array s',0)
578c
579      if(.not.ga_create(mt_int,np*mwmg,1,'iwz',mwmg,1,ga_iwz))
580     + call md_abort('Failed to create global array iwz',0)
581      if(.not.ga_create(mt_int,np*msag,1,'isz',msag,1,ga_isz))
582     + call md_abort('Failed to create global array isz',0)
583c
584      return
585      end
586      subroutine sp_free()
587c
588      implicit none
589c
590#include "sp_common.fh"
591#include "mafdecls.fh"
592#include "global.fh"
593c
594      if(.not.ga_destroy(ga_isz))
595     + call md_abort('Error ga_destroy isz',ga_isz)
596      if(.not.ga_destroy(ga_iwz))
597     + call md_abort('Error ga_destroy iwz',ga_iwz)
598c
599      if(.not.ga_destroy(ga_s))
600     + call md_abort('Error ga_destroy s',ga_s)
601      if(.not.ga_destroy(ga_is))
602     + call md_abort('Error ga_destroy is',ga_is)
603      if(.not.ga_destroy(ga_w))
604     + call md_abort('Error ga_destroy w',ga_w)
605      if(.not.ga_destroy(ga_iw))
606     + call md_abort('Error ga_destroy iw',ga_iw)
607      if(.not.ga_destroy(ga_ip))
608     + call md_abort('Error ga_destroy iw',ga_ip)
609c
610      if(.not.ma_pop_stack(l_sndx))
611     + call md_abort('Failed to deallocate sndx',0)
612      if(.not.ma_pop_stack(l_lseq))
613     + call md_abort('Failed to deallocate lseq',0)
614      if(.not.ma_pop_stack(l_jpl))
615     + call md_abort('Failed to deallocate jpl',0)
616      if(.not.ma_pop_stack(l_ipl))
617     + call md_abort('Failed to deallocate ipl',0)
618      if(.not.ma_pop_stack(l_rng))
619     + call md_abort('Failed to deallocate rng',0)
620      if(.not.ma_pop_stack(l_boxr))
621     + call md_abort('Failed to deallocate boxr',0)
622      if(.not.ma_pop_stack(l_boxs))
623     + call md_abort('Failed to deallocate boxs',0)
624      if(.not.ma_pop_stack(l_iown))
625     + call md_abort('Failed to deallocate iown',0)
626      if(.not.ma_pop_stack(l_buren))
627     + call md_abort('Failed to deallocate buren',0)
628      if(.not.ma_pop_stack(l_bindex))
629     + call md_abort('Failed to deallocate bindex',0)
630c
631      return
632      end
633      subroutine sp_rdrst(lfnrst,filrst,lfntop,filtop,
634     + temp,tempw,temps,ipl,xw,vw,fw,xwcr,iwl,xs,vs,fs,xscr,isl,
635     + bxw,bvw,bfw,brw,ibw,nw,bxs,bvs,bfs,ibs,ns,ndx,
636     + ibownr,boxsiz,lseq,isndx)
637c
638      implicit none
639c
640#include "sp_common.fh"
641#include "mafdecls.fh"
642#include "msgids.fh"
643c
644      integer sp_btop
645      external sp_btop
646c
647      integer lfnrst,lfntop
648      character*255 filrst,filtop
649      integer nw,ns
650      integer ibownr(maxbox,3)
651      integer ipl(mbox,mip2),iwl(mwm,miw2),isl(msa,mis2)
652      real*8 xw(mwm,3,mwa),xs(msa,3),xwcr(mwm,3)
653      real*8 vw(mwm,3,mwa),vs(msa,3),xscr(msm,3)
654      real*8 fw(mwm,3,mwa),fs(msa,3)
655      real*8 bxw(nw,3,mwa),bxs(ns,3),brw(nw,3)
656      real*8 bvw(nw,3,mwa),bvs(ns,3)
657      real*8 bfw(nw,3,mwa),bfs(ns,3)
658      integer ibw(nw,2),ibs(ns,mis2)
659      real*8 boxsiz(maxbox,3)
660      integer ndx(nw),lseq(mseq),isndx(mseq)
661      real*8 temp,tempw,temps
662c
663      character*1 cdum
664      real*8 rdum,cgx,cgy,cgz
665      integer i,j,k,idum,jdum,kdum,number,ncyc,numw
666      integer icyc,ibx,iby,ibz,ipx,ipy,ipz,node,new,nold
667      integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs
668      integer ili,ihi,jli,jhi,ilp,ihp,jlp,jhp
669c
670      integer nat,nqt,naw,nbw,nhw,ndw,now,ntw,nnw,nsmr
671      real*8 boxi(3)
672      logical lforces
673      character*80 card
674c
675      integer l,m,ib(3),nbox,nb(3),joff,npars,icount,lasts
676      real*8 xtmin,xtmax,xtx,xt(3)
677c
678      if(.not.ma_verify_allocator_stuff()) then
679      call md_abort('ERROR IN MEMORY USE ENTERING SP_RDRST',0)
680      endif
681c
682      lpbc9=.true.
683      nsmr=0
684      icount=0
685      lasts=0
686      ndums=0
687c
688      if(me.eq.0) then
689c
690      open(unit=lfnrst,file=filrst(1:index(filrst,' ')-1),
691     + status='old',form='formatted',err=9999)
692      rewind(lfnrst)
693c
694      open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
695     + status='old',form='formatted',err=9899)
696      rewind(lfntop)
697c
698      write(lfnout,6000)
699 6000 format(/,' TOPOLOGY FILE INFORMATION',/)
700c
701      read(lfntop,2001,end=9897,err=9898) card
702 2001 format(a80)
703      write(lfnout,6001) card
704 6001 format(' Title',t15,a)
705      read(lfntop,2001,end=9897,err=9898) card
706      write(lfnout,6002) card
707 6002 format(t15,a)
708      read(lfntop,2001,end=9897,err=9898) card
709      write(lfnout,6002) card
710c
711      read(lfntop,2001,end=9897,err=9898) card
712      write(lfnout,6003) card(1:12),card(13:32),card(33:42)
713 6003 format(' Version',t15,a,/,' Date',t15,a,/,' Force field',t15,a)
714c
715      read(lfntop,2002,end=9897,err=9898) npars
716      read(lfntop,2002,end=9897,err=9898) nat
717      read(lfntop,2002,end=9897,err=9898) nqt
718      read(lfntop,2002,end=9897,err=9898) nseq
719 2002 format(i5)
720      read(lfntop,2001,end=9897,err=9898)
721      do 1103 k=1,npars
722      do 102 i=1,nat
723      read(lfntop,2001,end=9897,err=9898)
724      do 103 j=i,nat
725      read(lfntop,2001,end=9897,err=9898)
726  103 continue
727  102 continue
728 1103 continue
729      do 104 i=1,nqt*npars
730      read(lfntop,2001,end=9897,err=9898)
731  104 continue
732      do 4104 i=1,nseq
733      read(lfntop,2001,end=9897,err=9898)
734 4104 continue
735      read(lfntop,2003,end=9897,err=9898) naw,nbw,nhw,ndw,now,ntw,nnw
736 2003 format(5i7,2i10)
737      read(lfntop,2001,end=9897,err=9898)
738      do 105 i=1,naw
739      read(lfntop,2001,end=9897,err=9898)
740  105 continue
741      do 106 i=1,nbw*(npars+1)
742      read(lfntop,2001,end=9897,err=9898)
743  106 continue
744      do 107 i=1,nhw*(npars+1)
745      read(lfntop,2001,end=9897,err=9898)
746  107 continue
747      do 108 i=1,ndw*(npars+1)
748      read(lfntop,2001,end=9897,err=9898)
749  108 continue
750      do 109 i=1,now*(npars+1)
751      read(lfntop,2001,end=9897,err=9898)
752  109 continue
753      if(ntw.gt.0) then
754      read(lfntop,2004,end=9897,err=9898)
755      read(lfntop,2004,end=9897,err=9898)
756 2004 format(11i7)
757      endif
758      if(nnw.gt.0) then
759      read(lfntop,2005,end=9997,err=9998)
760      read(lfntop,2005,end=9997,err=9998)
761 2005 format(11i7)
762      endif
763      read(lfntop,2001,end=9897,err=9898)
764      do 204 i=1,npars
765      read(lfntop,2001,end=9897,err=9898)
766  204 continue
767c
768      write(lfnout,6100)
769 6100 format(/,' RESTART FILE INFORMATION',/)
770c
771      read(lfnrst,1001,end=9997,err=9998) card
772 1001 format(a80)
773      write(lfnout,6101) card
774 6101 format(' Title',t15,a)
775      read(lfnrst,1001,end=9997,err=9998) card
776      write(lfnout,6102) card
777 6102 format(t15,a)
778      read(lfnrst,1001,end=9997,err=9998) card
779      write(lfnout,6002) card
780c
781      read(lfnrst,1016) card(1:32),nhist,lforces
782 1016 format(a32,i5,4x,l1)
783      write(lfnout,6103) card(1:12),card(13:32)
784 6103 format(' Version',t15,a,/,' Date',t15,a)
785c
786      if(nhist.gt.0) write(lfnout,6104)
787 6104 format(/,' History',/)
788c
789      if(nhist.gt.0) then
790      do 21 i=1,nhist
791      read(lfnrst,1017) hist(i)
792 1017 format(a)
793      write(lfnout,6105) hist(i)
794 6105 format(1x,a)
795   21 continue
796      endif
797      if(nhist.lt.mxhist) then
798      nhist=nhist+1
799      else
800      do 22 i=1,nhist-1
801      hist(i)=hist(i+1)
802   22 continue
803      endif
804      do 23 i=1,80
805      hist(nhist)(i:i)=' '
806   23 continue
807      read(lfnrst,1002,end=9997,err=9998) npbtyp,nbxtyp
808 1002 format(i5,i5)
809      read(lfnrst,1003,end=9997,err=9998) ((vlat(i,j),j=1,3),i=1,3)
810 1003 format(3f12.6)
811      read(lfnrst,1004,end=9997,err=9998) jdum
812 1004 format(40x,i5)
813      read(lfnrst,1005,end=9997,err=9998) temp,tempw,temps
814 1005 format(3f12.6)
815      do 2 i=1,3
816      box(i)=vlat(i,i)
817      boxh(i)=half*box(i)
818      boxi(i)=one/box(i)
819      do 3 j=1,3
820      vlati(i,j)=vlat(i,j)
821    3 continue
822    2 continue
823      call matinv(vlati,3,3)
824      if(jdum.ne.0) then
825      read(lfnrst,1001,end=9997,err=9998) cdum
826      endif
827      read(lfnrst,1006,end=9997,err=9998) idum
828 1006 format(70x,i5)
829      if(idum.gt.0) then
830      read(lfnrst,1007,end=9997,err=9998) idum,jdum,kdum
831 1007 format(3i5)
832      read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,idum)
833      read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,jdum)
834      read(lfnrst,1008,end=9997,err=9998) (rdum,i=1,kdum)
835 1008 format(4e20.12)
836      endif
837c
838      if(nwm.gt.0) then
839      number=0
840      ncyc=nwm/nw+1
841      numw=nw
842      do 4 icyc=1,ncyc
843      if(nwm-number.lt.numw) numw=nwm-number
844      do 44 i=1,numw
845      read(lfnrst,1009,end=9997,err=9998)
846     + ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
847 1009 format(2x,6f13.8)
848      if(lforces) read(lfnrst,1109,end=9997,err=9998)
849     + ((bfw(i,j,k),j=1,3),k=1,nwa)
850 1109 format(2x,6e13.6)
851      read(lfnrst,1010,end=9997,err=9998) ibw(i,2),(brw(i,k),k=1,3)
852 1010 format(i1,1x,3f13.8)
853   44 continue
854      do 5 i=1,numw
855      cgx=zero
856      cgy=zero
857      cgz=zero
858      do 6 k=1,nwa
859      cgx=cgx+bxw(i,1,k)
860      cgy=cgy+bxw(i,2,k)
861      cgz=cgz+bxw(i,3,k)
862      if(.not.lforces) then
863      bfw(i,1,k)=zero
864      bfw(i,2,k)=zero
865      bfw(i,3,k)=zero
866      endif
867    6 continue
868      ibx=0
869      iby=0
870      ibz=0
871      if(nbxtyp.ne.1) then
872      xt(1)=cgx
873      xt(2)=cgy
874      xt(3)=cgz
875      else
876      xt(1)=box(1)*(vlati(1,1)*cgx+vlati(1,2)*cgy+vlati(1,3)*cgz)
877      xt(2)=box(2)*(vlati(2,1)*cgx+vlati(2,2)*cgy+vlati(2,3)*cgz)
878      xt(3)=box(3)*(vlati(3,1)*cgx+vlati(3,2)*cgy+vlati(3,3)*cgz)
879      endif
880      do 7 j=1,nbx-1
881      if(xt(1)/nwa+boxh(1).gt.boxsiz(j,1)) ibx=j
882    7 continue
883      do 8 j=1,nby-1
884      if(xt(2)/nwa+boxh(2).gt.boxsiz(j,2)) iby=j
885    8 continue
886      do 9 j=1,nbz-1
887      if(xt(3)/nwa+boxh(3).gt.boxsiz(j,3)) ibz=j
888    9 continue
889c
890      if(npbtyp.gt.0) then
891      if(ibx.ge.nbx) ibx=ibx-nbx
892      if(iby.ge.nby) iby=iby-nby
893      if(ibx.lt.0) ibx=ibx+nbx
894      if(iby.lt.0) iby=iby+nby
895      if(npbtyp.eq.1) then
896      if(ibz.ge.nbz) ibz=ibz-nbz
897      if(ibz.lt.0) ibz=ibz+nbz
898      else
899      if(ibz.ge.nbz) ibz=nbz-1
900      if(ibz.lt.0) ibz=0
901      endif
902      else
903      if(ibx.ge.nbx) ibx=nbx-1
904      if(iby.ge.nby) iby=nby-1
905      if(ibz.ge.nbz) ibz=nbz-1
906      if(ibx.lt.0) ibx=0
907      if(iby.lt.0) iby=0
908      if(ibz.lt.0) ibz=0
909      endif
910      ipx=ibownr(ibx+1,1)
911      ipy=ibownr(iby+1,2)
912      ipz=ibownr(ibz+1,3)
913c
914      ndx(i)=(ipz*npy+ipy)*npx+ipx
915      ibw(i,1)=(ibz*nby+iby)*nbx+ibx
916      if(lnode0) then
917      ndx(i)=0
918      ibw(i,1)=0
919      endif
920    5 continue
921      do 10 node=0,np-1
922      new=0
923      do 11 i=1,numw
924      if(ndx(i).eq.node) then
925      new=new+1
926      iwl(new,lwgmn)=number+i
927      iwl(new,lwnod)=node
928      do 124 j=1,3
929      do 12 k=1,nwa
930      xw(new,j,k)=bxw(i,j,k)
931      vw(new,j,k)=bvw(i,j,k)
932      fw(new,j,k)=bfw(i,j,k)
933   12 continue
934      if(nserie.eq.0) then
935      xwcr(new,j)=zero
936      else
937      xwcr(new,j)=brw(i,j)
938      endif
939  124 continue
940      iwl(new,lwbox)=ibw(i,1)
941      iwl(new,lwdyn)=5*ibw(i,2)
942      endif
943   11 continue
944c
945      if(new.gt.0) then
946      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
947      call ga_get(ga_ip,ilp,ilp,jlp,jhp,ipl,mbox)
948      nold=ipl(1,2)
949      if(nold+new.gt.mwm) call md_abort('Dimension mwm too small',0)
950      ipl(1,2)=ipl(1,2)+new
951      call ga_put(ga_ip,ilp,ilp,jlp,jhp,ipl,mbox)
952      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
953      ili=ili+nold
954      ihi=ili+new-1
955      call ga_put(ga_iw,ili,ihi,jli,jhi,iwl,mwm)
956      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
957      ilw=ilw+nold
958      ihw=ilw+new-1
959      call ga_put(ga_w,ilw,ihw,jlw,jlw+3*mwa-1,xw(1,1,1),mwm)
960      call ga_put(ga_w,ilw,ihw,jlw+3*mwa,jlw+6*mwa-1,vw(1,1,1),mwm)
961      call ga_put(ga_w,ilw,ihw,jlw+6*mwa,jlw+6*mwa+2,xwcr(1,1),mwm)
962      call ga_put(ga_w,ilw,ihw,jlw+6*mwa+3,jlw+9*mwa+2,fw(1,1,1),mwm)
963      endif
964c
965   10 continue
966      number=number+numw
967    4 continue
968      endif
969c
970      if(nsa.gt.0) then
971      nb(1)=nbx
972      nb(2)=nby
973      nb(3)=nbz
974      nbox=0
975      k=1
976      joff=0
977      do 13 i=1,nsa
978      if(k.gt.ns) call md_abort('Increase memory for sp_rdrst buffer',0)
979      read(lfnrst,1011,end=9997,err=9998)
980     + ibs(k,11),(bxs(k,j),j=1,3),(bvs(k,j),j=1,3),ibs(k,12)
981 1011 format(i1,1x,6f13.8,i5)
982      if(ibs(k,12).lt.0) ndums=ndums+1
983      if(lforces) then
984      read(lfnrst,1111,end=9997,err=9998) (bfs(k,j),j=1,3)
985 1111 format(2x,6e13.6)
986      else
987      bfs(k,1)=zero
988      bfs(k,2)=zero
989      bfs(k,3)=zero
990      endif
991      read(lfntop,2009,end=9897,err=9898) (ibs(k,j),j=1,10)
992 2009 format(16x,i3,4i7,5i5)
993c 2009 format(16x,10i5)
994      if(ibs(k,3).ne.lasts) then
995      icount=icount+1
996      lasts=ibs(k,3)
997      isndx(icount)=lasts
998      endif
999      ibs(k,3)=icount
1000      if(nsmr.lt.ibs(k,2)) nsmr=ibs(k,2)
1001      if(nserie.eq.0) then
1002      if(nsf.lt.ibs(k,1)) nsf=ibs(k,1)
1003      else
1004      if(nsf.lt.ibs(k,1))
1005     + call md_abort('Error in number of solute fractions',nsf)
1006      endif
1007c
1008c     if segment of this atom is different from the segment of previous atom
1009c     then distribute all previous atoms in the list
1010c
1011      if(k.gt.1) then
1012      if(ibs(k,3).ne.ibs(k-1,3)) then
1013      new=k-1
1014      goto 14
1015      endif
1016      endif
1017c
1018c     if this is the last atom distribute
1019c
1020      if(i.eq.nsa) then
1021      new=k
1022      goto 14
1023      endif
1024c
1025c     read next atom
1026c
1027      k=k+1
1028      goto 13
1029c
1030c     distribute atoms 1 through new
1031c
1032   14 continue
1033c
1034c     determine the center of geometry
1035c
1036      do 15 l=1,3
1037      xtmax=bxs(1,l)
1038      xtmin=bxs(1,l)
1039      do 16 j=1,new
1040      xtmax=max(xtmax,bxs(j,l))
1041      xtmin=min(xtmin,bxs(j,l))
1042   16 continue
1043      xtx=0.5d0*(xtmax+xtmin)
1044      if(npbtyp.ne.0) then
1045      if(abs(xtx).gt.boxh(l)) then
1046      xtx=xtx-nint(xtx*boxi(l))*box(l)
1047      endif
1048      endif
1049      ib(l)=0
1050      do 17 m=1,nb(l)-1
1051      if(xtx+boxh(l).gt.boxsiz(m,l)) ib(l)=m
1052   17 continue
1053   15 continue
1054c
1055c     periodic boundaries
1056c
1057      if(npbtyp.gt.0) then
1058      m=2
1059      if(npbtyp.eq.1) m=3
1060      do 18 l=1,m
1061      if(ib(l).ge.nb(l)) ib(l)=ib(l)-nb(l)
1062      if(ib(l).lt.0) ib(l)=ib(l)+nb(l)
1063   18 continue
1064      if(npbtyp.gt.1) then
1065      if(ib(3).ge.nb(3)) ib(3)=nb(3)-1
1066      if(ib(3).lt.0) ib(3)=0
1067      endif
1068      else
1069      do 19 l=1,3
1070      if(ib(l).ge.nb(l)) ib(l)=nb(l)-1
1071      if(ib(l).lt.0) ib(l)=0
1072   19 continue
1073      endif
1074c
1075c     determine owning node
1076c
1077      if(.not.lnode0) nbox=(ib(3)*nb(2)+ib(2))*nb(1)+ib(1)
1078      node=sp_btop(nbox,ibownr)
1079c
1080      do 120 j=1,new
1081      isl(j,lsgan)=joff+j
1082      isl(j,lsfrc)=ibs(j,1)
1083      isl(j,lsmol)=ibs(j,2)
1084      isl(j,lssgm)=ibs(j,3)
1085      isl(j,lsgrp)=ibs(j,4)
1086      isl(j,lspgr)=ibs(j,5)
1087      isl(j,lsatt)=ibs(j,6)
1088      isl(j,lsct1)=ibs(j,7)
1089      isl(j,lsct2)=ibs(j,8)
1090      isl(j,lsct3)=ibs(j,9)
1091      isl(j,lssss)=ibs(j,10)
1092      isl(j,lsdyn)=5*ibs(j,11)
1093      isl(j,lsbox)=nbox
1094      isl(j,lsnod)=node
1095      isl(j,lshop)=0
1096      if(ibs(j,12).gt.0) isl(j,lshop)=ibs(j,12)*2
1097      if(ibs(j,12).lt.0) isl(j,lshop)=(-ibs(j,12))*2+1
1098      do 121 l=1,3
1099      xs(j,l)=bxs(j,l)
1100      vs(j,l)=bvs(j,l)
1101      fs(j,l)=bfs(j,l)
1102  121 continue
1103  120 continue
1104      joff=joff+new
1105c
1106c     communicate data to node
1107c
1108      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
1109      call ga_get(ga_ip,ilp,ilp+1,jlp,jhp,ipl,mbox)
1110      nold=ipl(2,2)
1111      if(nold+new.gt.msa)
1112     + call md_abort('Dimension msa too small (1)',nold+new)
1113      ipl(2,2)=ipl(2,2)+new
1114      call ga_put(ga_ip,ilp,ilp+1,jlp,jhp,ipl,mbox)
1115      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
1116      ili=ili+nold
1117      ihi=ili+new-1
1118      call ga_put(ga_is,ili,ihi,jli,jhi,isl,msa)
1119      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
1120      ils=ils+nold
1121      ihs=ils+new-1
1122      call ga_put(ga_s,ils,ihs,jls,jls+2,xs(1,1),msa)
1123      call ga_put(ga_s,ils,ihs,jls+3,jls+5,vs(1,1),msa)
1124      call ga_put(ga_s,ils,ihs,jls+9,jls+11,fs(1,1),msa)
1125c
1126c     make first atom of next segment first in the list
1127c
1128      if(k.gt.new) then
1129      do 122 j=1,12
1130      ibs(1,j)=ibs(k,j)
1131  122 continue
1132      do 123 j=1,3
1133      bxs(1,j)=bxs(k,j)
1134      bvs(1,j)=bvs(k,j)
1135      bfs(1,j)=bfs(k,j)
1136  123 continue
1137      k=2
1138      if(i.eq.nsa) then
1139      new=1
1140      k=1
1141      goto 14
1142      endif
1143      endif
1144c
1145   13 continue
1146      endif
1147c
1148      if(nsm.gt.0) then
1149      do 31 i=1,nsm
1150      xscr(i,1)=zero
1151      xscr(i,2)=zero
1152      xscr(i,3)=zero
1153   31 continue
1154      do 32 i=1,nsm
1155      read(lfnrst,1012,end=99,err=99) (xscr(i,j),j=1,3)
1156 1012 format(2x,3f13.8)
1157      if(nserie.eq.0) then
1158      xscr(i,1)=zero
1159      xscr(i,2)=zero
1160      xscr(i,3)=zero
1161      endif
1162   32 continue
1163      endif
1164c
1165      if(.not.ma_verify_allocator_stuff()) then
1166      call md_abort('ERROR IN MEMORY IN SP_RDRST',1)
1167      endif
1168c
1169      if(nseq.gt.0) then
1170      read(lfnrst,1013) (lseq(i),i=1,nseq)
1171 1013 format(20i3)
1172      endif
1173c
1174      if(.not.ma_verify_allocator_stuff()) then
1175      call md_abort('ERROR IN MEMORY IN SP_RDRST',2)
1176      endif
1177c
1178   99 continue
1179c
1180      close(unit=lfnrst,status='keep')
1181      close(unit=lfntop,status='keep')
1182c
1183      endif
1184c
1185      if(nseq.gt.0) then
1186      call ga_brdcst(msp_22,lseq,nseq*ma_sizeof(mt_int,1,mt_byte),0)
1187      endif
1188      call ga_brdcst(msp_28,ndums,ma_sizeof(mt_int,1,mt_byte),0)
1189c
1190      if(.not.ma_verify_allocator_stuff()) then
1191      call md_abort('ERROR IN MEMORY USE EXITING SP_RDRST',0)
1192      endif
1193c
1194      return
1195 9897 call md_abort('EOF encountered on topology file',0)
1196 9898 call md_abort('Error reading topology file',1)
1197 9899 call md_abort('Error opening topology file',2)
1198 9997 call md_abort('EOF encountered on restart file',0)
1199 9998 call md_abort('Error reading restart file',0)
1200 9999 call md_abort('Error opening restart file',0)
1201      return
1202      end
1203      integer function sp_btop(ibox,ibownr)
1204c
1205      implicit none
1206c
1207#include "sp_common.fh"
1208c
1209      integer ibox,ibownr(maxbox,3)
1210      integer iboxx,iboxy,iboxz
1211c
1212      iboxx=mod(ibox,nbx)
1213      iboxy=mod((ibox-iboxx)/nbx,nby)
1214      iboxz=((ibox-iboxx)/nbx-iboxy)/nby
1215      sp_btop=(ibownr(iboxz+1,3)*npy+ibownr(iboxy+1,2))*npx+
1216     + ibownr(iboxx+1,1)
1217c
1218      return
1219      end
1220      subroutine sp_decomp(ibownr,boxsiz,iburen,ibindx)
1221c
1222      implicit none
1223c
1224#include "sp_common.fh"
1225#include "msgids.fh"
1226#include "global.fh"
1227c
1228      integer ibownr(maxbox,3)
1229      real*8 boxsiz(maxbox,3)
1230      integer iburen(np,27,2)
1231      integer ibindx(np)
1232c
1233      integer ibx,iby,ibz,i,j,ix,jx,iy,jy,iz,jz,jnode,nrnod,nbiown
1234c
1235c     check dimensions of ibownr
1236c
1237      if(nbx.gt.maxbox.or.nby.gt.maxbox.or.nbz.gt.maxbox)
1238     + call md_abort('Dimension maxbox too small',0)
1239c
1240c     determine the node dimension for each sub box
1241c
1242      do 1 ibx=1,nbx
1243      ibownr(ibx,1)=((ibx-1)*npx)/nbx
1244      boxsiz(ibx,1)=(box(1)*dble(ibx))/dble(nbx)
1245    1 continue
1246      do 2 iby=1,nby
1247      ibownr(iby,2)=((iby-1)*npy)/nby
1248      boxsiz(iby,2)=(box(2)*dble(iby))/dble(nby)
1249    2 continue
1250      do 3 ibz=1,nbz
1251      ibownr(ibz,3)=((ibz-1)*npz)/nbz
1252      boxsiz(ibz,3)=(box(3)*dble(ibz))/dble(nbz)
1253    3 continue
1254      if(iand(idebug,1).eq.1) then
1255      write(lfndbg,8003) nbx,nby,nbz,maxbox,npx,npy,npz
1256 8003 format('ibownr in sp_decomp',7i5)
1257      write(lfndbg,8004) (ibownr(ibx,1),ibx=1,nbx)
1258 8004 format('ibownr x',/,(20i5))
1259      write(lfndbg,8005) (ibownr(iby,2),iby=1,nby)
1260 8005 format('ibownr y',/,(20i5))
1261      write(lfndbg,8006) (ibownr(ibz,3),ibz=1,nbz)
1262 8006 format('ibownr z',/,(20i5))
1263      write(lfndbg,8002)
1264 8002 format('boxlist')
1265      endif
1266c
1267c     determine neighboring nodes and store in neighb(27,2)
1268c     such that:
1269c
1270c     neighb(n,1) is the n-th neighbor of node me
1271c     neighb(n,2) is the node of which node me is the n-th neighbor
1272c
1273c     a value of -1 indicates that such node does not exist
1274c
1275      do 4 i=1,27
1276      neighb(i,1)=-1
1277      neighb(i,2)=-1
1278    4 continue
1279c
1280      do 5 ix=1,3
1281      jx=mex+ix-2
1282      if(npbtyp.gt.0) then
1283      if(npx.gt.2.and.jx.lt.0) jx=jx+npx
1284      if(npx.gt.2.and.jx.ge.npx) jx=jx-npx
1285      endif
1286      if(jx.ge.0.and.jx.lt.npx) then
1287      do 6 iy=1,3
1288      jy=mey+iy-2
1289      if(npbtyp.gt.0) then
1290      if(npy.gt.2.and.jy.lt.0) jy=jy+npy
1291      if(npy.gt.2.and.jy.ge.npy) jy=jy-npy
1292      endif
1293      if(jy.ge.0.and.jy.lt.npy) then
1294      do 7 iz=1,3
1295      jz=mez+iz-2
1296      if(npbtyp.eq.1) then
1297      if(npz.gt.2.and.jz.lt.0) jz=jz+npz
1298      if(npz.gt.2.and.jz.ge.npz) jz=jz-npz
1299      endif
1300      if(jz.ge.0.and.jz.lt.npz) then
1301      jnode=((jz*npy)+jy)*npx+jx
1302      neighb(3*(3*(ix-1)+(iy-1))+iz,1)=jnode
1303      neighb(3*(3*(3-ix)+(3-iy))+4-iz,2)=jnode
1304      endif
1305    7 continue
1306      endif
1307    6 continue
1308      endif
1309    5 continue
1310c
1311      nbiown=30
1312      do 8 ibx=1,nbx
1313      do 9 iby=1,nby
1314      do 10 ibz=1,nbz
1315      nrnod=(ibownr(ibz,3)*npy+ibownr(iby,2))*npx+ibownr(ibx,1)
1316      if(me.eq.nrnod) nbiown=nbiown+1
1317   10 continue
1318    9 continue
1319    8 continue
1320c
1321      do 11 j=1,27
1322      do 12 i=1,np
1323      iburen(i,j,1)=0
1324      iburen(i,j,2)=0
1325   12 continue
1326   11 continue
1327      do 13 i=1,np
1328      ibindx(i)=0
1329   13 continue
1330      do 14 j=1,27
1331      iburen(me+1,j,1)=neighb(j,1)
1332      iburen(me+1,j,2)=neighb(j,2)
1333      if(neighb(j,1).ge.0) ibindx(neighb(j,1)+1)=j
1334   14 continue
1335c
1336      call ga_igop(msp_27,iburen,np*54,'+')
1337c
1338c     broadcast the maximum number of sub-boxes per node to all nodes
1339c
1340      if(nbiown.gt.mbox) call md_abort('Error in mbox',0)
1341c
1342      return
1343      end
1344      subroutine sp_nrnode
1345c
1346      implicit none
1347c
1348#include "sp_common.fh"
1349c
1350      integer ix,iy,iz,npt,i,j,k,l
1351c
1352c     this routine distributes the available processes in the
1353c     Cartesian directions
1354c
1355c     npx : number of processes in x direction
1356c     npy : number of processes in y direction
1357c     npz : number of processes in z direction
1358c
1359c     determine node dimensions
1360c
1361      if(npx*npy*npz.ne.np) then
1362      if(npx+npy+npz.gt.0)
1363     + call md_abort('Specified npx*npy*npz ne np',0)
1364c
1365      npt=0
1366      do 1 i=1,np
1367      do 2 j=i,np
1368      do 3 k=j,np
1369      l=i*j*k
1370      if(l.eq.np) then
1371      if(l.eq.npt) then
1372      if(k.gt.npz) goto 3
1373      if(i+j+k.lt.npx+npy+npz) then
1374      npx=i
1375      npy=j
1376      npz=k
1377      endif
1378      goto 3
1379      else
1380      npt=np
1381      npx=i
1382      npy=j
1383      npz=k
1384      endif
1385      endif
1386    3 continue
1387    2 continue
1388    1 continue
1389c
1390      if(npx*npy*npz.ne.np) call md_abort('nrnode: code error',0)
1391      endif
1392c
1393c     determine processor location of me
1394c
1395      do 4 ix=1,npx
1396      do 5 iy=1,npy
1397      do 6 iz=1,npz
1398      if(me.eq.((iz-1)*npy+(iy-1))*npx+(ix-1)) then
1399      mex=ix-1
1400      mey=iy-1
1401      mez=iz-1
1402      endif
1403    6 continue
1404    5 continue
1405    4 continue
1406c
1407      return
1408      end
1409      subroutine sp_initip(ibownr,ipl)
1410c
1411      implicit none
1412c
1413#include "sp_common.fh"
1414#include "global.fh"
1415c
1416      integer ibownr(maxbox,3),ipl(mbox,mip2)
1417      integer i,j,ibx,iby,ibz,nrbox,nrnod,nbiown
1418      integer il,ih,jl,jh
1419c
1420      do 1 j=1,mip2
1421      do 2 i=1,30
1422      ipl(i,j)=0
1423    2 continue
1424    1 continue
1425c
1426      if(iand(idebug,1).eq.1) then
1427      write(lfndbg,8003) nbx,nby,nbz,maxbox
1428 8003 format('ibownr in sp_initip',4i5)
1429      write(lfndbg,8004) (ibownr(ibx,1),ibx=1,nbx)
1430 8004 format('ibownr x',/,(20i5))
1431      write(lfndbg,8005) (ibownr(iby,2),iby=1,nby)
1432 8005 format('ibownr y',/,(20i5))
1433      write(lfndbg,8006) (ibownr(ibz,3),ibz=1,nbz)
1434 8006 format('ibownr z',/,(20i5))
1435      call util_flush(lfndbg)
1436      if(iand(idebug,2).eq.2) write(lfndbg,8002)
1437 8002 format('boxlist')
1438      endif
1439      nbiown=0
1440      do 3 ibx=1,nbx
1441      do 4 iby=1,nby
1442      do 5 ibz=1,nbz
1443      nrbox=((ibz-1)*nby+iby-1)*nbx+ibx-1
1444      nrnod=(ibownr(ibz,3)*npy+ibownr(iby,2))*npx+ibownr(ibx,1)
1445      if(iand(idebug,2).eq.2) then
1446      write(lfndbg,8001) ibx,iby,ibz,nrbox,nrnod
1447 8001 format(5i5)
1448      call util_flush(lfndbg)
1449      endif
1450      if(me.eq.nrnod) then
1451      nbiown=nbiown+1
1452      ipl(30+nbiown,1)=nrbox
1453      ipl(30+nbiown,2)=0
1454      ipl(30+nbiown,3)=0
1455      ipl(30+nbiown,4)=0
1456      ipl(30+nbiown,5)=0
1457      endif
1458    5 continue
1459    4 continue
1460    3 continue
1461      ipl(1,1)=nbiown
1462c
1463      call ga_distribution(ga_ip,me,il,ih,jl,jh)
1464      call ga_put(ga_ip,il,ih,jl,jh,ipl,mbox)
1465c
1466      if(iand(idebug,2).eq.2) then
1467      write(lfndbg,8000) (i,ipl(30+i,1),i=1,ipl(1,1))
1468 8000 format('ipl',/,(2i5))
1469      call util_flush(lfndbg)
1470      endif
1471c
1472      return
1473      end
1474      subroutine sp_update_i(numsa,isl,numwm,iwl)
1475c
1476      implicit none
1477c
1478#include "sp_common.fh"
1479#include "mafdecls.fh"
1480#include "global.fh"
1481#include "msgids.fh"
1482c
1483      integer numsa,isl(msa,mis2)
1484      integer numwm,iwl(mwm,miw2)
1485c
1486      call sp_upd_i(numsa,isl,int_mb(i_pack),numwm,iwl,int_mb(i_packw))
1487c
1488      return
1489      end
1490      subroutine sp_upd_i(numsa,isl,islp,numwm,iwl,iwlp)
1491c
1492      implicit none
1493c
1494#include "sp_common.fh"
1495#include "global.fh"
1496#include "msgids.fh"
1497c
1498      integer numsa,isl(msa,mis2)
1499      integer numwm,iwl(mwm,miw2)
1500      integer iwlp(mwm,npackw),islp(msa,npack)
1501c
1502      integer il,ih,jl,jh
1503c
1504      if(numsa.gt.0) then
1505      call ga_distribution(ga_is,me,il,ih,jl,jh)
1506      if(npack.eq.0) then
1507      call ga_put(ga_is,il,il+numsa-1,jl,jh,isl,msa)
1508      else
1509      call sp_pack(numsa,isl,islp)
1510      call ga_put(ga_is,il,il+numsa-1,jl,jl+npack-1,islp,msa)
1511      endif
1512      endif
1513c
1514      if(numwm.gt.0) then
1515      call ga_distribution(ga_iw,me,il,ih,jl,jh)
1516      if(npackw.eq.0) then
1517      call ga_put(ga_iw,il,il+numwm-1,jl,jh,iwl,mwm)
1518      else
1519      call sp_packw(numwm,iwl,iwlp)
1520      call ga_put(ga_iw,il,il+numwm-1,jl,jl+npackw-1,iwlp,mwm)
1521      endif
1522      endif
1523c
1524      call ga_sync()
1525c
1526      return
1527      end
1528      subroutine sp_numbb(ibownr,boxsiz)
1529c
1530      implicit none
1531c
1532#include "sp_common.fh"
1533#include "mafdecls.fh"
1534#include "msgids.fh"
1535#include "global.fh"
1536#include "util.fh"
1537c
1538      integer ibownr(maxbox,3)
1539      real*8 boxsiz(maxbox,3)
1540      logical lside,leven
1541c
1542      integer ibx,iby,ibz,ipx,ipy,ipz,ibox,inode
1543      integer jbx,jby,jbz,kbx,kby,kbz,jpx,jpy,jpz
1544      integer ilx,ihx,ily,ihy,ilz,ihz
1545      integer jbox,jnode,mbblb,i
1546      real*8 dx,dxtmp,dy,dytmp,dz,dztmp,dist2
1547      character*255 string
1548c
1549      if(iand(idebug,2).eq.2) then
1550      write(lfndbg,'(a,i6)') 'boxsiz in sp_numbb'
1551      write(lfndbg,'(6f12.6)') (boxsiz(i,1),i=1,nbx)
1552      write(lfndbg,'(6f12.6)') (boxsiz(i,2),i=1,nby)
1553      write(lfndbg,'(6f12.6)') (boxsiz(i,3),i=1,nbz)
1554      call util_flush(lfndbg)
1555      endif
1556c
1557c     determine size of box-box pairlist
1558c     ----------------------------------
1559c
1560      mbbl=0
1561      mbblb=0
1562      do 1 ibx=0,nbx-1
1563      ipx=ibownr(ibx+1,1)
1564      do 2 iby=0,nby-1
1565      ipy=ibownr(iby+1,2)
1566      do 3 ibz=0,nbz-1
1567      ipz=ibownr(ibz+1,3)
1568      ibox=(ibz*nby+iby)*nbx+ibx
1569      inode=(ipz*npy+ipy)*npx+ipx
1570      if(inode.eq.me) then
1571      do 4 jbx=0,nbx-1
1572      kbx=jbx-ibx
1573      jpx=ibownr(jbx+1,1)
1574c
1575      if(ibx.le.jbx) then
1576      ilx=ibx
1577      ihx=jbx
1578      else
1579      ilx=jbx
1580      ihx=ibx
1581      endif
1582c
1583      dx=zero
1584      if(ibx.ne.jbx) then
1585      dx=boxsiz(ihx,1)-boxsiz(ilx+1,1)
1586      if(npbtyp.gt.0) then
1587      dxtmp=zero
1588      if(ilx.gt.0) dxtmp=boxsiz(ilx,1)
1589      if(ihx.lt.nbx-1) dxtmp=dxtmp-boxsiz(ihx+1,1)+box(1)
1590      if(dxtmp.lt.dx) dx=dxtmp
1591      if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx
1592      if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx
1593      endif
1594      endif
1595c
1596      do 5 jby=0,nby-1
1597      kby=jby-iby
1598      jpy=ibownr(jby+1,2)
1599c
1600      if(iby.le.jby) then
1601      ily=iby
1602      ihy=jby
1603      else
1604      ily=jby
1605      ihy=iby
1606      endif
1607c
1608      dy=zero
1609      if(iby.ne.jby) then
1610      dy=boxsiz(ihy,2)-boxsiz(ily+1,2)
1611      if(npbtyp.gt.0) then
1612      dytmp=zero
1613      if(ily.gt.0) dytmp=boxsiz(ily,2)
1614      if(ihy.lt.nby-1) dytmp=dytmp-boxsiz(ihy+1,2)+box(2)
1615      if(dytmp.lt.dy) dy=dytmp
1616      if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby
1617      if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby
1618      endif
1619      endif
1620c
1621      do 6 jbz=0,nbz-1
1622      kbz=jbz-ibz
1623      jpz=ibownr(jbz+1,3)
1624c
1625      if(ibz.le.jbz) then
1626      ilz=ibz
1627      ihz=jbz
1628      else
1629      ilz=jbz
1630      ihz=ibz
1631      endif
1632c
1633      dz=zero
1634      if(ibz.ne.jbz) then
1635      dz=boxsiz(ihz,3)-boxsiz(ilz+1,3)
1636      if(npbtyp.eq.1) then
1637      dztmp=zero
1638      if(ilz.gt.0) dztmp=boxsiz(ilz,3)
1639      if(ihz.lt.nbz-1) dztmp=dztmp-boxsiz(ihz+1,3)+box(3)
1640      if(dztmp.lt.dz) dz=dztmp
1641      if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz
1642      if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz
1643      endif
1644      endif
1645c
1646      jbox=(jbz*nby+jby)*nbx+jbx
1647      jnode=(jpz*npy+jpy)*npx+jpx
1648c
1649c     determine orientation jbox in relation to ibox
1650c
1651c     lside is true if
1652c
1653c     i: 0  j: 0  k:  +
1654c     i: 0  j:  + k:-0+
1655c     i:  + j:-0+ k:-0+
1656c
1657      lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0)
1658     + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0
1659c
1660c     determine if ibox is identical to jbox
1661c
1662c     lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0
1663c
1664c     determine if difference in box numbers is even or odd
1665c
1666      leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox)
1667c
1668c     calculate the distance between the two boxes
1669c
1670      if(nbxtyp.eq.1) then
1671      dist2=
1672     + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+
1673     + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+
1674     + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2
1675      else
1676      dist2=dx*dx+dy*dy+dz*dz
1677      endif
1678c
1679c     keep half of the box pairs
1680c
1681c     this test also appears in sp_numbb
1682c     any changes need to be made in both routines
1683c
1684      if((inode.eq.jnode.and.ibox.ge.jbox).or.(inode.ne.jnode.and.
1685     + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then
1686c
1687c     keep only those within maximum cutoff distance
1688c
1689      if(rbbl*rbbl.gt.dist2) mbbl=mbbl+1
1690      endif
1691      if(rbbl*rbbl.gt.dist2) mbblb=mbblb+1
1692 6    continue
1693 5    continue
1694 4    continue
1695      endif
1696 3    continue
1697 2    continue
1698 1    continue
1699c
1700      mbbl=mbbl+1
1701      if(mbblb+1.gt.mbbl) mbbl=mbblb+1
1702c
1703      if(np.gt.1) call ga_igop(msp_20,mbbl,1,'max')
1704      if(me.eq.0) then
1705      if(util_print('distribution',print_debug)) then
1706      write(lfnout,2001) mbbl
1707 2001 format(' Dimension of the cell-cell list is ',i7)
1708      endif
1709      endif
1710c
1711c     determine size of box-box pairlist assuming minimum cell sizes
1712c     --------------------------------------------------------------
1713c
1714      mbbl=0
1715      mbblb=0
1716      do 11 ibx=0,nbx-1
1717      ipx=ibownr(ibx+1,1)
1718      do 12 iby=0,nby-1
1719      ipy=ibownr(iby+1,2)
1720      do 13 ibz=0,nbz-1
1721      ipz=ibownr(ibz+1,3)
1722      ibox=(ibz*nby+iby)*nbx+ibx
1723      inode=(ipz*npy+ipy)*npx+ipx
1724      if(inode.eq.me) then
1725      do 14 jbx=0,nbx-1
1726      kbx=jbx-ibx
1727      jpx=ibownr(jbx+1,1)
1728c
1729      if(ibx.le.jbx) then
1730      ilx=ibx
1731      ihx=jbx
1732      else
1733      ilx=jbx
1734      ihx=ibx
1735      endif
1736c
1737      dx=zero
1738      if(ibx.ne.jbx) then
1739      dx=dble(ihx-ilx-1)*bxmin
1740      if(npbtyp.gt.0) then
1741      dxtmp=zero
1742      if(ilx.gt.0) dxtmp=dble(ilx)*bxmin
1743      if(ihx.lt.nbx-1) dxtmp=dxtmp+dble(nbx-ihx-1)*bxmin
1744      if(dxtmp.lt.dx) dx=dxtmp
1745      if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx
1746      if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx
1747      endif
1748      endif
1749c
1750      do 15 jby=0,nby-1
1751      kby=jby-iby
1752      jpy=ibownr(jby+1,2)
1753c
1754      if(iby.le.jby) then
1755      ily=iby
1756      ihy=jby
1757      else
1758      ily=jby
1759      ihy=iby
1760      endif
1761c
1762      dy=zero
1763      if(iby.ne.jby) then
1764      dy=dble(ihy-ily-1)*bymin
1765      if(npbtyp.gt.0) then
1766      dytmp=zero
1767      if(ily.gt.0) dytmp=dble(ily)*bymin
1768      if(ihy.lt.nby-1) dytmp=dytmp+dble(nby-ihy-1)*bymin
1769      if(dytmp.lt.dy) dy=dytmp
1770      if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby
1771      if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby
1772      endif
1773      endif
1774c
1775      do 16 jbz=0,nbz-1
1776      kbz=jbz-ibz
1777      jpz=ibownr(jbz+1,3)
1778c
1779      if(ibz.le.jbz) then
1780      ilz=ibz
1781      ihz=jbz
1782      else
1783      ilz=jbz
1784      ihz=ibz
1785      endif
1786c
1787      dz=zero
1788      if(ibz.ne.jbz) then
1789      dz=dble(ihz-ilz-1)*bzmin
1790      if(npbtyp.eq.1) then
1791      dztmp=zero
1792      if(ilz.gt.0) dztmp=dble(ilz)*bzmin
1793      if(ihz.lt.nbz-1) dztmp=dztmp-dble(nbz-ihz-1)*bzmin
1794      if(dztmp.lt.dz) dz=dztmp
1795      if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz
1796      if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz
1797      endif
1798      endif
1799c
1800      jbox=(jbz*nby+jby)*nbx+jbx
1801      jnode=(jpz*npy+jpy)*npx+jpx
1802c
1803c     determine orientation jbox in relation to ibox
1804c
1805c     lside is true if
1806c
1807c     i: 0  j: 0  k:  +
1808c     i: 0  j:  + k:-0+
1809c     i:  + j:-0+ k:-0+
1810c
1811      lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0)
1812     + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0
1813c
1814c     determine if ibox is identical to jbox
1815c
1816c     lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0
1817c
1818c     determine if difference in box numbers is even or odd
1819c
1820      leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox)
1821c
1822c     calculate the distance between the two boxes
1823c
1824      if(nbxtyp.eq.1) then
1825      dist2=
1826     + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+
1827     + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+
1828     + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2
1829      else
1830      dist2=dx*dx+dy*dy+dz*dz
1831      endif
1832c
1833c     keep half of the box pairs
1834c
1835c     this test also appears in sp_numbb
1836c     any changes need to be made in both routines
1837c
1838      if((inode.eq.jnode.and.ibox.ge.jbox).or.(inode.ne.jnode.and.
1839     + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then
1840c
1841c     keep only those within maximum cutoff distance
1842c
1843      if(rbbl*rbbl.gt.dist2) mbbl=mbbl+1
1844      endif
1845      if(rbbl*rbbl.gt.dist2) mbblb=mbblb+1
1846   16 continue
1847   15 continue
1848   14 continue
1849      endif
1850   13 continue
1851   12 continue
1852   11 continue
1853c
1854      mbbl=mbbl+1
1855      if(mbblb+1.gt.mbbl) mbbl=mbblb+1
1856c
1857c
1858c     allocate memory for the box-box list
1859c
1860      if(np.gt.1) call ga_igop(msp_20,mbbl,1,'max')
1861c
1862      if(me.eq.0) then
1863      if(util_print('distribution',print_debug)) then
1864      write(lfnout,2002) mbbl
1865 2002 format(' Dimension of the cell-cell list is ',i7)
1866      endif
1867      endif
1868c
1869      if(mbblp.eq.0) then
1870      mbbl=max(mbbl,mbbreq)
1871      if(.not.ma_push_get(mt_int,mbb2*mbbl,'bb',l_bb,i_bb))
1872     + call md_abort('Failed to allocate memory for bb',0)
1873      mbblp=mbbl
1874      else
1875      if(mbbl.gt.mbblp) then
1876      write(string,1111) mbblp,mbbl,mbbreq
1877 1111 format('error: lbbl increase from ',i6,' to ',i6,'(',i6,')')
1878      call md_abort(string,me)
1879c      call md_abort('lbbl increased beyond allocated memory',me)
1880      endif
1881      endif
1882c
1883      return
1884      end
1885      subroutine sp_listbb(ibownr,boxsiz,lbbl)
1886c
1887      implicit none
1888c
1889#include "sp_common.fh"
1890#include "global.fh"
1891#include "msgids.fh"
1892c
1893      integer ibownr(maxbox,3),lbbl(mbbl,mbb2)
1894      real*8 boxsiz(maxbox,3)
1895      logical lside,leven
1896c
1897      integer ibx,iby,ibz,ipx,ipy,ipz,ibox,inode
1898      integer jbx,jby,jbz,kbx,kby,kbz,jpx,jpy,jpz
1899      integer ilx,ihx,ily,ihy,ilz,ihz
1900      integer i,j,jbox,jnode,ltemp
1901      real*8 dx,dxtmp,dy,dytmp,dz,dztmp,dist2
1902c
1903c     Construction of the box-box pairlist
1904c
1905      nbbl=0
1906      do 1 ibx=0,nbx-1
1907      ipx=ibownr(ibx+1,1)
1908      do 2 iby=0,nby-1
1909      ipy=ibownr(iby+1,2)
1910      do 3 ibz=0,nbz-1
1911      ipz=ibownr(ibz+1,3)
1912      ibox=(ibz*nby+iby)*nbx+ibx
1913      inode=(ipz*npy+ipy)*npx+ipx
1914      if(inode.eq.me) then
1915      do 4 jbx=0,nbx-1
1916      kbx=jbx-ibx
1917      jpx=ibownr(jbx+1,1)
1918c
1919      if(ibx.le.jbx) then
1920      ilx=ibx
1921      ihx=jbx
1922      else
1923      ilx=jbx
1924      ihx=ibx
1925      endif
1926c
1927      dx=zero
1928      if(ibx.ne.jbx) then
1929      dx=boxsiz(ihx,1)-boxsiz(ilx+1,1)
1930      if(npbtyp.gt.0) then
1931      dxtmp=zero
1932      if(ilx.gt.0) dxtmp=boxsiz(ilx,1)
1933      if(ihx.lt.nbx-1) dxtmp=dxtmp-boxsiz(ihx+1,1)+box(1)
1934      if(dxtmp.lt.dx) dx=dxtmp
1935      if(kbx.gt.0.and.kbx.gt.iabs(kbx-nbx)) kbx=kbx-nbx
1936      if(kbx.lt.0.and.-kbx.gt.iabs(kbx+nbx)) kbx=kbx+nbx
1937      endif
1938      endif
1939c
1940      do 5 jby=0,nby-1
1941      kby=jby-iby
1942      jpy=ibownr(jby+1,2)
1943c
1944      if(iby.le.jby) then
1945      ily=iby
1946      ihy=jby
1947      else
1948      ily=jby
1949      ihy=iby
1950      endif
1951c
1952      dy=zero
1953      if(iby.ne.jby) then
1954      dy=boxsiz(ihy,2)-boxsiz(ily+1,2)
1955      if(npbtyp.gt.0) then
1956      dytmp=zero
1957      if(ily.gt.0) dytmp=boxsiz(ily,2)
1958      if(ihy.lt.nby-1) dytmp=dytmp-boxsiz(ihy+1,2)+box(2)
1959      if(dytmp.lt.dy) dy=dytmp
1960      if(kby.gt.0.and.kby.gt.iabs(kby-nby)) kby=kby-nby
1961      if(kby.lt.0.and.-kby.gt.iabs(kby+nby)) kby=kby+nby
1962      endif
1963      endif
1964c
1965      do 6 jbz=0,nbz-1
1966      kbz=jbz-ibz
1967      jpz=ibownr(jbz+1,3)
1968c
1969      if(ibz.le.jbz) then
1970      ilz=ibz
1971      ihz=jbz
1972      else
1973      ilz=jbz
1974      ihz=ibz
1975      endif
1976c
1977      dz=zero
1978      if(ibz.ne.jbz) then
1979      dz=boxsiz(ihz,3)-boxsiz(ilz+1,3)
1980      if(npbtyp.eq.1) then
1981      dztmp=zero
1982      if(ilz.gt.0) dztmp=boxsiz(ilz,3)
1983      if(ihz.lt.nbz-1) dztmp=dztmp-boxsiz(ihz+1,3)+box(3)
1984      if(dztmp.lt.dz) dz=dztmp
1985      if(kbz.gt.0.and.kbz.gt.iabs(kbz-nbz)) kbz=kbz-nbz
1986      if(kbz.lt.0.and.-kbz.gt.iabs(kbz+nbz)) kbz=kbz+nbz
1987      endif
1988      endif
1989c
1990      jbox=(jbz*nby+jby)*nbx+jbx
1991      jnode=(jpz*npy+jpy)*npx+jpx
1992c
1993c     determine orientation jbox in relation to ibox
1994c
1995c     lside is true if
1996c
1997c     i: 0  j: 0  k:  +
1998c     i: 0  j:  + k:-0+
1999c     i:  + j:-0+ k:-0+
2000c
2001      lside=(kbx.eq.0.and.kby.eq.0.and.kbz.ge.0)
2002     + .or.(kbx.eq.0.and.kby.gt.0) .or. kbx.gt.0
2003c
2004c     determine if ibox is identical to jbox
2005c
2006c     lsame=kbx.eq.0.and.kby.eq.0.and.kbz.eq.0
2007c
2008c     determine if difference in box numbers is even or odd
2009c
2010      leven=2*(iabs(ibox-jbox)/2).eq.iabs(ibox-jbox)
2011c
2012c     calculate the distance between the two boxes
2013c
2014      if(nbxtyp.eq.1) then
2015      dist2=
2016     + (vlat(1,1)*dx/box(1)+vlat(1,2)*dy/box(2)+vlat(1,3)*dz/box(3))**2+
2017     + (vlat(2,1)*dx/box(1)+vlat(2,2)*dy/box(2)+vlat(2,3)*dz/box(3))**2+
2018     + (vlat(3,1)*dx/box(1)+vlat(3,2)*dy/box(2)+vlat(3,3)*dz/box(3))**2
2019      else
2020      dist2=dx*dx+dy*dy+dz*dz
2021      endif
2022c
2023c     keep half of the box pairs
2024c
2025c     this test also appears in sp_numbb
2026c     any changes need to be made in both routines
2027c
2028      if((inode.eq.jnode.and.ibox.ge.jbox).or. (inode.ne.jnode.and.
2029     + ((lside.and.leven).or.(.not.lside.and..not.leven)))) then
2030c
2031c     keep only those within maximum cutoff distance
2032c
2033      if(rbbl*rbbl.gt.dist2) then
2034      nbbl=nbbl+1
2035      if(nbbl.gt.mbbl) call md_abort('Box-box list too small',mbbl)
2036      lbbl(nbbl,1)=jnode
2037      lbbl(nbbl,2)=jbox
2038      lbbl(nbbl,3)=ibox
2039      lbbl(nbbl,4)=0
2040      endif
2041      endif
2042 6    continue
2043 5    continue
2044 4    continue
2045      endif
2046 3    continue
2047 2    continue
2048 1    continue
2049      npprev=0
2050c
2051      nbbloc=0
2052      do 7 i=1,nbbl-1
2053      do 8 j=i+1,nbbl
2054      if((lbbl(i,1).ne.me.and.lbbl(j,1).eq.me).or.
2055     + (lbbl(i,1).gt.lbbl(j,1).and.lbbl(i,1).ne.me).or.
2056     + (lbbl(i,1).eq.lbbl(j,1).and.lbbl(i,2).gt.lbbl(j,2)).or.
2057     + (lbbl(i,1).eq.lbbl(j,1).and.lbbl(i,2).eq.lbbl(j,2).and.
2058     + lbbl(i,3).gt.lbbl(j,3))) then
2059      ltemp=lbbl(i,1)
2060      lbbl(i,1)=lbbl(j,1)
2061      lbbl(j,1)=ltemp
2062      ltemp=lbbl(i,2)
2063      lbbl(i,2)=lbbl(j,2)
2064      lbbl(j,2)=ltemp
2065      ltemp=lbbl(i,3)
2066      lbbl(i,3)=lbbl(j,3)
2067      lbbl(j,3)=ltemp
2068      endif
2069    8 continue
2070      if(lbbl(i,1).eq.me) nbbloc=i
2071    7 continue
2072      if(lbbl(nbbl,1).eq.me) nbbloc=nbbl
2073c
2074      nrempr=0
2075      if(nbget.ne.0.and.np.gt.1) then
2076      nrempr=1
2077      do 9 i=2,nbbl
2078      if(lbbl(i,2).ne.lbbl(i-1,2)) nrempr=nrempr+1
2079    9 continue
2080      call ga_igop(msp_29,nrempr,1,'max')
2081      nrempr=min(2*nrempr,nrempr+25)
2082      endif
2083c
2084      if(iand(idebug,2).eq.2) then
2085      write(lfndbg,8000) (i,(lbbl(i,j),j=1,3),i=1,nbbl)
2086 8000 format('lbbl',/,(4i5))
2087      call util_flush(lfndbg)
2088      endif
2089c
2090      return
2091      end
2092      subroutine sp_finish()
2093c
2094      implicit none
2095c
2096#include "sp_common.fh"
2097#include "mafdecls.fh"
2098c
2099      if(.not.ma_pop_stack(l_xscr))
2100     + call md_abort('Failed to deallocate xscr',0)
2101      if(.not.ma_pop_stack(l_bb))
2102     + call md_abort('Failed to deallocate memory for bb',0)
2103c
2104      call sp_free()
2105c
2106      return
2107      end
2108      subroutine sp_initf(fw,fs,llng,iwz,isz,lpair)
2109c
2110      implicit none
2111c
2112#include "sp_common.fh"
2113#include "mafdecls.fh"
2114#include "global.fh"
2115c
2116      real*8 fw(mwm,3,mwa,2),fs(msa,3,2)
2117      integer iwz(mwm),isz(msa)
2118      logical llng,lpair
2119c
2120      integer i,j,k,l,m,il,ih,jl,jh
2121c
2122      llong=llng
2123c
2124      m=1
2125      if(llong) m=2
2126c
2127      do 1 l=1,m
2128      if(nwm.gt.0) then
2129      do 2 k=1,mwa
2130      do 3 j=1,3
2131      do 4 i=1,mwm
2132      fw(i,j,k,l)=zero
2133    4 continue
2134    3 continue
2135    2 continue
2136      endif
2137      if(nsa.gt.0) then
2138      do 5 j=1,3
2139      do 6 i=1,msa
2140      fs(i,j,l)=zero
2141    6 continue
2142    5 continue
2143      endif
2144    1 continue
2145c
2146      if(nwm.gt.0) then
2147      call ga_distribution(ga_w,me,il,ih,jl,jh)
2148      call ga_put(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm)
2149      if(llong) call ga_put(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2,
2150     + fw(1,1,1,2),mwm)
2151      endif
2152      if(nsa.gt.0) then
2153      call ga_distribution(ga_s,me,il,ih,jl,jh)
2154      call ga_put(ga_s,il,ih,jl+6,jl+8,fs,msa)
2155      if(llong) call ga_put(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa)
2156      endif
2157c
2158      if(lpair) then
2159      do 7 i=1,mwm
2160      iwz(i)=0
2161    7 continue
2162      do 8 i=1,msa
2163      isz(i)=0
2164    8 continue
2165      call ga_zero(ga_iwz)
2166      call ga_zero(ga_isz)
2167      endif
2168c
2169      return
2170      end
2171      subroutine sp_copyg(fw,fs)
2172c
2173      implicit none
2174c
2175#include "sp_common.fh"
2176#include "global.fh"
2177c
2178      real*8 fw(mwm,3,mwa),fs(msa,3)
2179c
2180      integer il,ih,jl,jh
2181c
2182      if(nwm.gt.0) then
2183      call ga_distribution(ga_w,me,il,ih,jl,jh)
2184      call ga_put(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm)
2185      endif
2186      if(nsa.gt.0) then
2187      call ga_distribution(ga_s,me,il,ih,jl,jh)
2188      call ga_put(ga_s,il,ih,jl+6,jl+8,fs,msa)
2189      endif
2190c
2191      return
2192      end
2193      subroutine sp_final(fw,fs,lpair,iwz,isz)
2194c
2195      implicit none
2196c
2197#include "sp_common.fh"
2198#include "mafdecls.fh"
2199#include "global.fh"
2200c
2201      real*8 fw(mwm,3,mwa,2),fs(msa,3,2)
2202      logical lpair
2203      integer iwz(mwm),isz(msa)
2204c
2205      integer i,il,ih,jl,jh
2206c
2207      if(np.gt.0) then
2208      if(nwm.gt.0) then
2209      call ga_distribution(ga_w,me,il,ih,jl,jh)
2210      call ga_acc(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm,one)
2211      if(llong) call ga_acc(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2,
2212     + fw(1,1,1,2),mwm,one)
2213      call ga_get(ga_w,il,ih,jl+6*mwa+3,jl+9*mwa+2,fw,mwm)
2214      if(ltwin) call ga_get(ga_w,il,ih,jl+9*mwa+3,jl+12*mwa+2,
2215     + fw(1,1,1,2),mwm)
2216      endif
2217      if(nsa.gt.0) then
2218      call ga_distribution(ga_s,me,il,ih,jl,jh)
2219      call ga_acc(ga_s,il,ih,jl+6,jl+8,fs,msa,one)
2220      if(llong) call ga_acc(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa,one)
2221      call ga_get(ga_s,il,ih,jl+6,jl+8,fs,msa)
2222      if(ltwin) call ga_get(ga_s,il,ih,jl+9,jl+11,fs(1,1,2),msa)
2223      endif
2224      endif
2225c
2226      if(lpair) then
2227      if(nwm.gt.0) then
2228      call ga_distribution(ga_iwz,me,il,ih,jl,jh)
2229      call ga_acc(ga_iwz,il,ih,1,1,iwz,mwm,1)
2230      call ga_get(ga_iwz,il,ih,1,1,iwz,mwm)
2231      do 1 i=1,nwmloc
2232      iwz(i)=min(1,iwz(i))
2233    1 continue
2234      endif
2235      if(nsa.gt.0) then
2236      call ga_distribution(ga_isz,me,il,ih,jl,jh)
2237      call ga_acc(ga_isz,il,ih,1,1,isz,msa,1)
2238      call ga_get(ga_isz,il,ih,1,1,isz,msa)
2239      do 2 i=1,nsaloc
2240      isz(i)=min(1,isz(i))
2241    2 continue
2242      endif
2243      endif
2244c
2245      return
2246      end
2247      subroutine sp_wrtrst(lfnrst,filrst,lveloc,
2248     + pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs,xscr,prjct,
2249     + lseq)
2250c
2251      implicit none
2252c
2253#include "sp_common.fh"
2254#include "mafdecls.fh"
2255#include "global.fh"
2256c
2257      integer lfnrst
2258      character*255 filrst
2259      logical lveloc
2260      integer iwl(mwm,miw2),isl(msa,mis2),lseq(mseq)
2261      real*8 pres,temp,tempw,temps
2262      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
2263      real*8 xs(msa,3),vs(msa,3),fs(msa,3),xscr(msm,3)
2264      character*80 prjct
2265c
2266      integer lenscr
2267c
2268      project=prjct
2269c
2270      lenscr=ma_inquire_avail(mt_byte)/
2271     + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+
2272     + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1
2273      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx))
2274     + call md_abort('Failed to allocate bx',0)
2275      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv))
2276     + call md_abort('Failed to allocate bv',0)
2277      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf))
2278     + call md_abort('Failed to allocate bf',0)
2279      if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br))
2280     + call md_abort('Failed to allocate br',0)
2281      if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi))
2282     + call md_abort('Failed to allocate bi',0)
2283      if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n))
2284     + call md_abort('Failed to allocate n',0)
2285c
2286      call sp_wtrst(lfnrst,filrst,lveloc,pres,temp,tempw,temps,
2287     + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs,
2288     + xscr,int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
2289     + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx),
2290     + dbl_mb(i_bv),dbl_mb(i_bf),lseq)
2291c
2292      if(.not.ma_pop_stack(l_n))
2293     + call md_abort('Failed to deallocate n',0)
2294      if(.not.ma_pop_stack(l_bi))
2295     + call md_abort('Failed to deallocate bi',0)
2296      if(.not.ma_pop_stack(l_br))
2297     + call md_abort('Failed to deallocate br',0)
2298      if(.not.ma_pop_stack(l_bf))
2299     + call md_abort('Failed to deallocate bf',0)
2300      if(.not.ma_pop_stack(l_bv))
2301     + call md_abort('Failed to deallocate bv',0)
2302      if(.not.ma_pop_stack(l_bx))
2303     + call md_abort('Failed to deallocate bx',0)
2304c
2305      return
2306      end
2307      subroutine sp_wtrst(lfnrst,filrst,lveloc,pres,temp,tempw,temps,
2308     + iwl,iwlp,xw,vw,fw,xwcr,isl,islp,xs,vs,fs,xscr,
2309     + ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs,lseq)
2310c
2311      implicit none
2312c
2313#include "sp_common.fh"
2314#include "mafdecls.fh"
2315#include "global.fh"
2316c
2317      integer lfnrst,nb
2318      character*255 filrst
2319      logical lveloc
2320      real*8 pres,temp,tempw,temps
2321      integer iwl(mwm,miw2),isl(msa,mis2),lseq(mseq)
2322      integer iwlp(mwm,npackw),islp(msa,npack)
2323      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
2324      real*8 xs(msa,3),vs(msa,3),fs(msa,3),xscr(msm,3)
2325      integer ipl(mbox,mip2),ibw(nb),ibs(nb,2)
2326      real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3)
2327      real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3)
2328c
2329      integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan
2330      integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw
2331      integer ils,ihs,jls,jhs
2332      character*10 rdate,rtime
2333      character*18 user
2334#ifdef USE_POSIXF
2335      integer ilen,ierror
2336#endif
2337      integer idyn,idynp,ihop
2338      logical lforces
2339c
2340      lforces=iguide.ne.0
2341c
2342      if(ga_nodeid().eq.0) then
2343c
2344      call swatch(rdate,rtime)
2345#ifdef USE_POSIXF
2346      call pxfgetlogin(user, ilen, ierror)
2347#else
2348      call getlog(user)
2349#endif
2350      if(user(18:18).ne.' ') user='                  '
2351c
2352      rewind(lfnrst)
2353      write(lfnrst,1000)
2354 1000 format('Restart file',/,' ',/,' ')
2355      write(lfnrst,1001) 4.2,rdate,rtime,nhist,lforces
2356 1001 format(f12.6,2a10,i5,4x,l1)
2357      hist(nhist)(1:18)=user
2358      hist(nhist)(19:28)=rdate
2359      hist(nhist)(29:48)=rtime
2360      hist(nhist)(49:108)=project(1:60)
2361      do 10 i=1,nhist
2362      write(lfnrst,1009) hist(i)
2363 1009 format(a)
2364   10 continue
2365      write(lfnrst,1002) npbtyp,nbxtyp,rsgm,((vlat(i,j),j=1,3),i=1,3)
2366 1002 format(2i5,f12.6,/,(3f12.6))
2367      write(lfnrst,1003) pres
2368 1003 format(1pe12.5)
2369      write(lfnrst,1004) temp,tempw,temps
2370 1004 format(3f12.6)
2371      write(lfnrst,1005) nwm,nwa,nsm,nsa,nwmc,nsf,nseq,0,0
2372 1005 format(7i10,2i5)
2373c
2374      if(nwm.gt.0) then
2375      number=0
2376      ncyc=nwm/nb+1
2377      numw=nb
2378      do 1 icyc=1,ncyc
2379      if(nwm-number.lt.numw) numw=nwm-number
2380c
2381c     begin test code 10/31/2001
2382c     initialize ibw to check that all atoms have been received
2383c
2384      do 1112 i=1,nb
2385      ibw(i)=-1
2386 1112 continue
2387c
2388c     end test code
2389c
2390      do 2 node=np-1,0,-1
2391      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
2392      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
2393      nwmn=ipl(1,2)
2394      if(nwmn.gt.0) then
2395      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
2396      if(npackw.eq.0) then
2397      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm)
2398      else
2399      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm)
2400      call sp_unpackw(nwmn,iwl,iwlp)
2401      endif
2402      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
2403      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
2404      if(lveloc)
2405     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
2406      if(lforces)
2407     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm)
2408      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa,jlw+6*mwa+2,xwcr,mwm)
2409      do 3 i=1,nwmn
2410      j=iwl(i,lwgmn)-number
2411      if(j.gt.0.and.j.le.numw) then
2412      do 4 k=1,nwa
2413      bxw(j,1,k)=xw(i,1,k)
2414      bxw(j,2,k)=xw(i,2,k)
2415      bxw(j,3,k)=xw(i,3,k)
2416      bvw(j,1,k)=vw(i,1,k)
2417      bvw(j,2,k)=vw(i,2,k)
2418      bvw(j,3,k)=vw(i,3,k)
2419      if(lforces) then
2420      bfw(j,1,k)=fw(i,1,k)
2421      bfw(j,2,k)=fw(i,2,k)
2422      bfw(j,3,k)=fw(i,3,k)
2423      endif
2424    4 continue
2425      brw(j,1)=xwcr(i,1)
2426      brw(j,2)=xwcr(i,2)
2427      brw(j,3)=xwcr(i,3)
2428      ibw(j)=iwl(i,lwdyn)
2429      endif
2430    3 continue
2431      endif
2432    2 continue
2433      do 5 i=1,numw
2434      if(lveloc) then
2435      write(lfnrst,1006) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
2436      else
2437      write(lfnrst,1006) ((bxw(i,j,k),j=1,3),(zero,j=1,3),k=1,nwa)
2438      endif
2439 1006 format(2x,6f13.8)
2440      if(lforces) write(lfnrst,1106) ((bfw(i,j,k),j=1,3),k=1,nwa)
2441 1106 format(2x,6e13.6)
2442      idyn=iand(ibw(i),12)/4
2443      idynp=iand(ibw(i),3)
2444      write(lfnrst,1007) idynp,idyn,(brw(i,k),k=1,3)
2445 1007 format(2i1,3f13.8)
2446c
2447c     begin test code 10/31/2001
2448c     check if al atoms have been received
2449c
2450      if(ibw(i).lt.0)
2451     + call md_abort('Missing solvent in wtrst',i)
2452c
2453c     end test code
2454c
2455    5 continue
2456      number=number+numw
2457    1 continue
2458      endif
2459c
2460      if(nsa.gt.0) then
2461      number=0
2462      ncyc=nsa/nb+1
2463      nums=nb
2464      do 6 icyc=1,ncyc
2465      if(nsa-number.lt.nums) nums=nsa-number
2466c
2467c     begin test code 10/31/2001
2468c     initialize ibw to check that all atoms have been received
2469c
2470      do 1117 i=1,nb
2471      ibs(i,1)=-1
2472      ibs(i,2)=0
2473 1117 continue
2474c
2475c     end test code
2476c
2477      do 7 node=np-1,0,-1
2478      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
2479      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
2480      nsan=ipl(2,2)
2481      if(nsan.gt.0) then
2482      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
2483      if(npack.eq.0) then
2484      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa)
2485      else
2486      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
2487      call sp_unpack(nsan,isl,islp)
2488      endif
2489      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
2490      call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
2491      if(lveloc) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
2492      if(lforces) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
2493      do 8 i=1,nsan
2494      j=isl(i,lsgan)-number
2495      if(j.gt.0.and.j.le.nums) then
2496      bxs(j,1)=xs(i,1)
2497      bxs(j,2)=xs(i,2)
2498      bxs(j,3)=xs(i,3)
2499      bvs(j,1)=vs(i,1)
2500      bvs(j,2)=vs(i,2)
2501      bvs(j,3)=vs(i,3)
2502      if(lforces) then
2503      bfs(j,1)=fs(i,1)
2504      bfs(j,2)=fs(i,2)
2505      bfs(j,3)=fs(i,3)
2506      endif
2507      ibs(j,1)=isl(i,lsdyn)
2508      ibs(j,2)=isl(i,lshop)
2509      endif
2510    8 continue
2511      endif
2512    7 continue
2513      do 9 i=1,nums
2514      idyn=iand(ibs(i,1),12)/4
2515      idynp=iand(ibs(i,1),3)
2516      ihop=ibs(i,2)
2517      if(iand(ihop,1).eq.1) then
2518      ihop=-(ihop/2)
2519      else
2520      ihop=ihop/2
2521      endif
2522      if(lveloc) then
2523      write(lfnrst,1008) idynp,idyn,(bxs(i,j),j=1,3),(bvs(i,j),j=1,3),
2524     + ihop
2525      else
2526      write(lfnrst,1008) idynp,idyn,(bxs(i,j),j=1,3),(zero,j=1,3),ihop
2527      endif
2528 1008 format(2i1,6f13.8,i5)
2529      if(lforces) write(lfnrst,1108) (bfs(i,j),j=1,3)
2530 1108 format(2x,3e13.6)
2531c
2532c     begin test code 10/31/2001
2533c     check if al atoms have been received
2534c
2535      if(ibs(i,1).lt.0)
2536     + call md_abort('Missing solute atom in wtrst',i)
2537c
2538c     end test code
2539c
2540    9 continue
2541      number=number+nums
2542    6 continue
2543      endif
2544c
2545      if(nsm.gt.0) then
2546      do 21 i=1,nsm
2547      write(lfnrst,1109) (xscr(i,j),j=1,3)
2548 1109 format(2x,3f13.8)
2549   21 continue
2550      endif
2551c
2552      if(nseq.gt.0) then
2553      write(lfnrst,1013) (lseq(i),i=1,nseq)
2554 1013 format(20i3)
2555      endif
2556c
2557      endif
2558c
2559      return
2560 9999 continue
2561      call md_abort('Failed to open restart file',me)
2562      return
2563      end
2564      subroutine sp_print()
2565c
2566      implicit none
2567c
2568#include "sp_common.fh"
2569#include "mafdecls.fh"
2570c
2571      integer i_lcnt,l_lcnt
2572c
2573      if(me.eq.0) then
2574      write(lfnout,1000)
2575 1000 format(/,' DOMAIN DECOMPOSITION',/)
2576c
2577      write(lfnout,1001) np,npx,npy,npz
2578 1001 format(' Processor count ',i5,' =',i5,' x',i5,' x',i5)
2579      write(lfnout,1002) nbx*nby*nbz,nbx,nby,nbz
2580 1002 format(' Cell count      ',i5,' =',i5,' x',i5,' x',i5)
2581c
2582      if(mod(nbx,npx)+mod(nby,npy)+mod(nbz,npz).ne.0) then
2583      write(lfnout,1003)
2584 1003 format(/,' WARNING: Inefficient distribution of cells over ',
2585     + 'processors')
2586      endif
2587c
2588      write(lfnout,1004) bxmin,bymin,bzmin
2589 1004 format(/,' Minimum cell size ',f12.6,2(' x',f12.6))
2590c
2591      if(nred(1)+nred(2)+nred(3).gt.0) then
2592      write(lfnout,1005)
2593 1005 format(/,' Warning: ',/)
2594      if(nred(1).gt.0) write(lfnout,1006) 'x',nred(1)
2595      if(nred(2).gt.0) write(lfnout,1006) 'y',nred(2)
2596      if(nred(3).gt.0) write(lfnout,1006) 'z',nred(3)
2597 1006 format(' Reduced number of cells in ',a,'-dimension: ',i5)
2598      endif
2599c
2600      if(nable.eq.1) write(lfnout,1007)
2601 1007 format(/,' Read previous box pair list',/)
2602c
2603      if(nable.eq.2) write(lfnout,1008)
2604 1008 format(/,' Unable to read previous box pair list',/)
2605c
2606      endif
2607c
2608      if(.not.ma_push_get(mt_int,3*np,'lcnt',l_lcnt,i_lcnt))
2609     + call md_abort('Failed to allocate memory for lcnt',0)
2610c
2611      call sp_prtcnt(int_mb(i_lcnt))
2612c
2613      if(.not.ma_pop_stack(l_lcnt))
2614     + call md_abort('Failed to deallocate lcnt',0)
2615c
2616      return
2617      end
2618      subroutine sp_prtcnt(lcnt)
2619c
2620      implicit none
2621c
2622#include "sp_common.fh"
2623#include "msgids.fh"
2624#include "global.fh"
2625c
2626      integer lcnt(np,3)
2627c
2628      integer i,j
2629c
2630      do 1 i=1,np
2631      lcnt(i,1)=0
2632      lcnt(i,2)=0
2633      lcnt(i,3)=0
2634    1 continue
2635c
2636      lcnt(me+1,1)=mbxloc
2637      lcnt(me+1,2)=nwmloc*nwa
2638      lcnt(me+1,3)=nsaloc
2639c
2640      if(np.gt.1) call ga_igop(msp_08,lcnt,3*np,'+')
2641c
2642      if(me.eq.0) then
2643      write(lfnout,1000)
2644 1000 format(/,' Initial distribution p:b(w+s)',/)
2645      write(lfnout,1001) (i-1,(lcnt(i,j),j=1,3),i=1,np)
2646 1001 format(4(3x,i4,':',i5,'(',i7,'+',i7,')'))
2647      write(lfnout,1002) mwm,msa
2648 1002 format(/,' Dimension workarrays solvent ',i6,/,
2649     + 22x,'solute  ',i6)
2650      endif
2651c
2652      return
2653      end
2654      subroutine sp_printf(filtop,lfntop,isl,xs,fs,
2655     + npener,esa)
2656c
2657      implicit none
2658c
2659#include "sp_common.fh"
2660#include "mafdecls.fh"
2661#include "msgids.fh"
2662#include "global.fh"
2663c
2664      integer lfntop
2665      character*70 filtop
2666      integer isl(msa,mis2),npener
2667      real*8 xs(msa,3),fs(msa,3,2),esa(nsa,2)
2668c
2669      integer lenscr
2670c
2671      lenscr=ma_inquire_avail(mt_byte)/
2672     + ((6*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+
2673     + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1
2674      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx))
2675     + call md_abort('Failed to allocate bx',0)
2676      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv))
2677     + call md_abort('Failed to allocate bv',0)
2678      if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br))
2679     + call md_abort('Failed to allocate br',0)
2680      if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi))
2681     + call md_abort('Failed to allocate bi',0)
2682      if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n))
2683     + call md_abort('Failed to allocate n',0)
2684c
2685      call sp_prt_s(filtop,lfntop,
2686     + int_mb(i_ipl),isl,int_mb(i_pack),xs,fs,
2687     + lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
2688     + npener,esa)
2689c
2690      if(.not.ma_pop_stack(l_n))
2691     + call md_abort('Failed to deallocate n',0)
2692      if(.not.ma_pop_stack(l_bi))
2693     + call md_abort('Failed to deallocate bi',0)
2694      if(.not.ma_pop_stack(l_br))
2695     + call md_abort('Failed to deallocate br',0)
2696      if(.not.ma_pop_stack(l_bv))
2697     + call md_abort('Failed to deallocate bv',0)
2698      if(.not.ma_pop_stack(l_bx))
2699     + call md_abort('Failed to deallocate bx',0)
2700c
2701      return
2702      end
2703      subroutine sp_prt_s(filtop,lfntop,
2704     + ipl,isl,islp,xs,fs,nb,ibs,bxs,bfs,
2705     + npener,esa)
2706c
2707      implicit none
2708c
2709#include "sp_common.fh"
2710#include "msgids.fh"
2711#include "global.fh"
2712c
2713      integer lfntop,nb,npener
2714      character*70 filtop
2715      integer ipl(mbox,mip2),isl(msa,mis2),islp(msa,npack),ibs(nb)
2716      real*8 xs(msa,3),fs(msa,3)
2717      real*8 bxs(nb,3),bfs(nb,3)
2718      real*8 esa(nsa,2)
2719c
2720      integer i,j,k,number,icyc,ncyc,node,nums,nsan,idyn,ism,iss
2721      integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ils,ihs,jls,jhs
2722      integer naw,nbw,nhw,ndw,now,ntw,nnw,nat,nqt,idum,npars
2723      character*1 cdum
2724      character*16 cat
2725c
2726      call ga_distribution(ga_ip,me,ilp,ihp,jlp,jhp)
2727      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
2728      nsan=ipl(2,2)
2729c
2730      if(nsan.gt.0) then
2731      call ga_distribution(ga_s,me,ils,ihs,jls,jhs)
2732      call ga_put(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
2733      endif
2734c
2735      call ga_sync()
2736c
2737      if(me.ne.0) return
2738c
2739      if(npener.le.0) then
2740      write(lfnout,1007)
2741 1007 format(//,' Solute coordinates and forces',//,
2742     + '   mol  segment     atom  dt     x       y       z ',
2743     + '             fx          fy          fz',/)
2744      else
2745      write(lfnout,1008)
2746 1008 format(//,' Solute coordinates, forces and energies',//,
2747     + '   mol  segment     atom  dt     x       y       z ',
2748     + '             fx          fy          fz        enb',
2749     + '          eb          ep',/)
2750      endif
2751c
2752      open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
2753     + status='old',form='formatted')
2754      rewind(lfntop)
2755      do 101 i=1,4
2756      read(lfntop,2001) cdum
2757 2001 format(a1)
2758  101 continue
2759      read(lfntop,2002) npars
2760      read(lfntop,2002) nat
2761      read(lfntop,2002) nqt
2762      read(lfntop,2002) nseq
2763 2002 format(i5)
2764      read(lfntop,2001) cdum
2765      do 1102 k=1,npars
2766      do 102 i=1,nat
2767      read(lfntop,2001) cdum
2768      do 103 j=i,nat
2769      read(lfntop,2001) cdum
2770  103 continue
2771  102 continue
2772 1102 continue
2773      do 104 i=1,nqt*npars
2774      read(lfntop,2001) cdum
2775  104 continue
2776      do 4104 i=1,nseq
2777      read(lfntop,2001) cdum
2778 4104 continue
2779      read(lfntop,2003) naw,nbw,nhw,ndw,now,ntw,nnw
2780 2003 format(5i7,2i10)
2781      read(lfntop,2001) cdum
2782      do 105 i=1,naw
2783      read(lfntop,2001) cdum
2784  105 continue
2785      do 106 i=1,nbw*(npars+1)
2786      read(lfntop,2001) cdum
2787  106 continue
2788      do 107 i=1,nhw*(npars+1)
2789      read(lfntop,2001) cdum
2790  107 continue
2791      do 108 i=1,ndw*(npars+1)
2792      read(lfntop,2001) cdum
2793  108 continue
2794      do 109 i=1,now*(npars+1)
2795      read(lfntop,2001) cdum
2796  109 continue
2797      if(ntw.gt.0) then
2798      read(lfntop,2004) (idum,i=1,ntw)
2799      read(lfntop,2004) (idum,i=1,ntw)
2800 2004 format(11i7)
2801      endif
2802      if(nnw.gt.0) then
2803      read(lfntop,2005) (idum,i=1,nnw)
2804      read(lfntop,2005) (idum,i=1,nnw)
2805 2005 format(11i7)
2806      endif
2807      read(lfntop,2001) cdum
2808      do 204 i=1,npars
2809      read(lfntop,2001) cdum
2810  204 continue
2811c
2812      if(nsa.gt.0) then
2813      number=0
2814      ncyc=nsa/nb+1
2815      nums=nb
2816      do 6 icyc=1,ncyc
2817      if(nsa-number.lt.nums) nums=nsa-number
2818      do 7 node=np-1,0,-1
2819      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
2820      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
2821      nsan=ipl(2,2)
2822      if(nsan.gt.0) then
2823      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
2824      if(npack.eq.0) then
2825      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa)
2826      else
2827      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
2828      call sp_unpack(nsan,isl,islp)
2829      endif
2830      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
2831      call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
2832      call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
2833      do 8 i=1,nsan
2834      j=isl(i,lsgan)-number
2835      if(j.gt.0.and.j.le.nums) then
2836      bxs(j,1)=xs(i,1)
2837      bxs(j,2)=xs(i,2)
2838      bxs(j,3)=xs(i,3)
2839      bfs(j,1)=fs(i,1)
2840      bfs(j,2)=fs(i,2)
2841      bfs(j,3)=fs(i,3)
2842      ibs(j)=isl(i,lsdyn)
2843      endif
2844    8 continue
2845      endif
2846    7 continue
2847      do 9 i=1,nums
2848      read(lfntop,2009) cat,ism,iss
2849 2009 format(a16,3x,2i7)
2850c 2009 format(a16,3x,2i7)
2851      idyn=iand(ibs(i),3)
2852      if(npener.le.0) then
2853      write(lfnout,1009) ism,iss,cat,idyn,
2854     + (bxs(i,j),j=1,3),(bfs(i,j),j=1,3)
2855 1009 format(2i5,':',a16,i1,1x,3f8.4,3x,3f12.3)
2856      else
2857      write(lfnout,1010) ism,iss,cat,idyn,
2858     + (bxs(i,j),j=1,3),(bfs(i,j),j=1,3),esa(i,1),esa(i,2),
2859     + esa(i,1)+esa(i,2)
2860 1010 format(2i5,':',a16,i1,1x,3f8.4,3x,3f12.3,3(2x,1pe10.3))
2861      endif
2862    9 continue
2863    6 continue
2864      endif
2865c
2866      close(unit=lfntop,status='keep')
2867c
2868      return
2869      end
2870      subroutine sp_wtrest(lfn)
2871c
2872      implicit none
2873c
2874#include "sp_common.fh"
2875#include "mafdecls.fh"
2876c
2877      integer lfn
2878c
2879      integer i_ltmp,l_ltmp
2880c
2881      if(.not.ma_push_get(mt_int,(mbbl+1)*mbb2,'ltmp',l_ltmp,i_ltmp))
2882     + call md_abort('Failed to allocate memory for ltmp',0)
2883      call sp_wrest(lfn,int_mb(i_bb),int_mb(i_ltmp),mbbl+1,
2884     + dbl_mb(i_boxs))
2885      if(.not.ma_pop_stack(l_ltmp))
2886     + call md_abort('Failed to deallocate ltmp',0)
2887c
2888      return
2889      end
2890      subroutine sp_wrest(lfn,lbbl,ltemp,mdim,boxsiz)
2891c
2892      implicit none
2893c
2894#include "sp_common.fh"
2895#include "msgids.fh"
2896#include "global.fh"
2897c
2898      integer lfn,mdim
2899      integer lbbl(mbbl,mbb2),ltemp(mdim,mbb2)
2900      real*8 boxsiz(maxbox,3)
2901c
2902      integer i,j,k
2903c
2904      if(me.eq.0) then
2905      write(lfn,1000)
2906 1000 format('restart space')
2907      write(lfn,1001) np,mbbl,npx,npy,npz,nbx,nby,nbz,np
2908 1001 format(9i7)
2909      write(lfn,1002) (boxsiz(i,1),i=1,nbx)
2910      write(lfn,1002) (boxsiz(i,2),i=1,nby)
2911      write(lfn,1002) (boxsiz(i,3),i=1,nbz)
2912 1002 format(4e20.12)
2913      endif
2914c
2915      do 1 i=1,np
2916      if(i.eq.me+1) then
2917      do 2 k=1,mbb2
2918      ltemp(1,k)=0
2919      do 3 j=1,nbbl
2920      ltemp(j+1,k)=lbbl(j,k)
2921    3 continue
2922    2 continue
2923      ltemp(1,1)=nbbl
2924      else
2925      do 4 k=1,mbb2
2926      do 5 j=1,mdim
2927      ltemp(j,k)=0
2928    5 continue
2929    4 continue
2930      endif
2931c
2932      call ga_igop(msp_09,ltemp,mdim*mbb2,'+')
2933c
2934      if(me.eq.0) then
2935      write(lfn,1003) i-1,ltemp(1,1)
2936 1003 format(2i7)
2937      do 6 j=1,ltemp(1,1)
2938      write(lfn,1004) (ltemp(j+1,k),k=1,4)
2939 1004 format(8i10)
2940    6 continue
2941      endif
2942c
2943    1 continue
2944c
2945      return
2946      end
2947      logical function sp_rdrest(lfn,fil,boxsiz)
2948c
2949      implicit none
2950c
2951#include "sp_common.fh"
2952#include "mafdecls.fh"
2953c
2954      integer lfn
2955      character*255 fil
2956      real*8 boxsiz(maxbox,3)
2957c
2958      logical sp_rrest
2959      external sp_rrest
2960c
2961      integer i_ltmp,l_ltmp
2962c
2963      if(.not.ma_push_get(mt_int,(mbbl+1)*mbb2,'ltmp',l_ltmp,i_ltmp))
2964     + call md_abort('Failed to allocate memory for ltmp',0)
2965      sp_rdrest=sp_rrest(lfn,fil,int_mb(i_bb),int_mb(i_ltmp),
2966     + mbbl+1,boxsiz)
2967      if(.not.ma_pop_stack(l_ltmp))
2968     + call md_abort('Failed to deallocate ltmp',0)
2969c
2970      return
2971      end
2972      logical function sp_rrest(lfn,fil,lbbl,ltemp,mdim,boxsiz)
2973c
2974      implicit none
2975c
2976#include "sp_common.fh"
2977#include "mafdecls.fh"
2978#include "msgids.fh"
2979c
2980      integer lfn,mdim
2981      character*255 fil
2982      integer lbbl(mbbl,mbb2),ltemp(mdim,mbb2)
2983c
2984      integer i,j,k,node,npp,nbl,nbytes,nbxs
2985      integer npxp,npyp,npzp,nbxp,nbyp,nbzp,mbblx
2986      character*13 string
2987      real*8 boxsiz(maxbox,3)
2988c
2989      if(me.eq.0) then
2990      open(unit=lfn,file=fil(1:index(fil,' ')-1),
2991     + status='old',form='formatted',err=9999)
2992      rewind(lfn)
2993c
2994    1 continue
2995      npp=0
2996      read(lfn,1000,end=9997) string
2997 1000 format(a13)
2998      if(string.ne.'restart space') goto 1
2999      read(lfn,1001) npp,mbblx,npxp,npyp,npzp,nbxp,nbyp,nbzp,nbxs
3000 1001 format(9i7)
3001      if(mbblp.lt.mbblx) npp=0
3002      if(npxp.ne.npx) npp=0
3003      if(npyp.ne.npy) npp=0
3004      if(npzp.ne.npz) npp=0
3005      if(nbxp.ne.nbx) npp=0
3006      if(nbyp.ne.nby) npp=0
3007      if(nbzp.ne.nbz) npp=0
3008      if(nbxs.gt.0.and.npp.gt.0) then
3009      read(lfn,1002) (boxsiz(i,1),i=1,nbxp)
3010      read(lfn,1002) (boxsiz(i,2),i=1,nbyp)
3011      read(lfn,1002) (boxsiz(i,3),i=1,nbzp)
3012 1002 format(4e20.12)
3013      endif
3014 9997 continue
3015      do 2 i=1,mbb2
3016      ltemp(1,i)=0
3017    2 continue
3018      endif
3019c
3020      nbytes=ma_sizeof(mt_int,1,mt_byte)
3021      call ga_brdcst(msp_10,npp,nbytes,0)
3022c
3023      if(np.ne.npp) goto 9998
3024c
3025      nbytes=ma_sizeof(mt_dbl,3*maxbox,mt_byte)
3026      call ga_brdcst(msp_10,boxsiz,nbytes,0)
3027c
3028      nbytes=mdim*mbb2*ma_sizeof(mt_int,1,mt_byte)
3029c
3030      do 3 i=1,np
3031c
3032      if(me.eq.0) then
3033      read(lfn,1003) node,nbl
3034 1003 format(2i7)
3035      ltemp(1,1)=node
3036      ltemp(1,2)=nbl
3037      do 4 j=1,nbl
3038      read(lfn,1004) (ltemp(j+1,k),k=1,4)
3039 1004 format(8i10)
3040      if(mbb2.gt.4) then
3041      do 44 k=5,mbb2
3042      ltemp(j+1,k)=0
3043   44 continue
3044      endif
3045    4 continue
3046      endif
3047c
3048      call ga_brdcst(msp_11,ltemp,nbytes,0)
3049c
3050      if(ltemp(1,1).eq.me) then
3051      nbbl=ltemp(1,2)
3052      do 5 k=1,mbb2
3053      do 6 j=1,nbbl
3054      lbbl(j,k)=ltemp(j+1,k)
3055    6 continue
3056    5 continue
3057      endif
3058c
3059    3 continue
3060c
3061      nable=1
3062      sp_rrest=.true.
3063      return
3064 9998 continue
3065      nable=2
3066      sp_rrest=.false.
3067      return
3068 9999 continue
3069      call md_abort('Failed to open restart file',0)
3070      sp_rrest=.false.
3071      return
3072      end
3073      subroutine sp_wrttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs,
3074     + stime,pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs)
3075c
3076      implicit none
3077c
3078#include "sp_common.fh"
3079#include "mafdecls.fh"
3080#include "global.fh"
3081c
3082      integer lfntrj
3083      logical lxw,lvw,lfw,lxs,lvs,lfs
3084      integer iwl(mwm,miw2),isl(msa,mis2)
3085      real*8 stime,pres,temp,tempw,temps
3086      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
3087      real*8 xs(msa,3),vs(msa,3),fs(msa,3)
3088c
3089      integer lenscr
3090c
3091      lenscr=ma_inquire_avail(mt_byte)/
3092     + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+
3093     + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1
3094      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx))
3095     + call md_abort('Failed to allocate bx',0)
3096      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv))
3097     + call md_abort('Failed to allocate bv',0)
3098      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf))
3099     + call md_abort('Failed to allocate bf',0)
3100      if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br))
3101     + call md_abort('Failed to allocate br',0)
3102      if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi))
3103     + call md_abort('Failed to allocate bi',0)
3104      if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n))
3105     + call md_abort('Failed to allocate n',0)
3106c
3107      call sp_wttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs,
3108     + stime,pres,temp,tempw,temps,
3109     + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs,
3110     + int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
3111     + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
3112     + dbl_mb(i_bf))
3113c
3114      if(.not.ma_pop_stack(l_n))
3115     + call md_abort('Failed to deallocate n',0)
3116      if(.not.ma_pop_stack(l_bi))
3117     + call md_abort('Failed to deallocate bi',0)
3118      if(.not.ma_pop_stack(l_br))
3119     + call md_abort('Failed to deallocate br',0)
3120      if(.not.ma_pop_stack(l_bf))
3121     + call md_abort('Failed to deallocate bf',0)
3122      if(.not.ma_pop_stack(l_bv))
3123     + call md_abort('Failed to deallocate bv',0)
3124      if(.not.ma_pop_stack(l_bx))
3125     + call md_abort('Failed to deallocate bx',0)
3126c
3127      return
3128      end
3129      subroutine sp_wttrj(lfntrj,lxw,lvw,lfw,lxs,lvs,lfs,
3130     + stime,pres,temp,tempw,temps,iwl,iwlp,xw,vw,fw,xwcr,
3131     + isl,islp,xs,vs,fs,ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs)
3132c
3133      implicit none
3134c
3135#include "sp_common.fh"
3136#include "mafdecls.fh"
3137#include "global.fh"
3138c
3139      integer lfntrj,nb
3140      logical lxw,lvw,lfw,lxs,lvs,lfs
3141      real*8 stime,pres,temp,tempw,temps
3142      integer iwl(mwm,miw2),isl(msa,mis2)
3143      integer iwlp(mwm,npackw),islp(msa,npack)
3144      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
3145      real*8 xs(msa,3),vs(msa,3),fs(msa,3)
3146      integer ipl(mbox,mip2),ibw(nb),ibs(nb)
3147      real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3)
3148      real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3)
3149c
3150      integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan
3151      integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw
3152      integer ils,ihs,jls,jhs
3153      character*10 rdate,rtime
3154c
3155      logical lpw,lps
3156c
3157      real*8 dumdst,dist
3158      integer nonh
3159c
3160      dumdst=0.02d0
3161c
3162      lpw=.false.
3163      lps=.false.
3164c
3165      if(me.eq.0) then
3166c
3167      call swatch(rdate,rtime)
3168c
3169      write(lfntrj,1000)
3170 1000 format('frame')
3171      write(lfntrj,1001) stime,temp,pres,rdate,rtime
3172 1001 format(2f12.6,1pe12.5,1x,2a10)
3173      write(lfntrj,1002) ((vlat(i,j),j=1,3),i=1,3)
3174 1002 format(3f12.6)
3175      write(lfntrj,1003) lxw,lvw,lfw,lpw,lxs,lvs,lfs,lps,nwm,nwa,nsa
3176 1003 format(8l1,3i10)
3177c
3178      if((lxw.or.lvw.or.lfw).and.nwm.gt.0) then
3179      number=0
3180      ncyc=nwm/nb+1
3181      numw=nb
3182      do 1 icyc=1,ncyc
3183      if(nwm-number.lt.numw) numw=nwm-number
3184c
3185c     begin test code 10/31/2001
3186c     initialize ibw to check that all atoms have been received
3187c
3188      do 1112 i=1,nb
3189      ibw(i)=-1
3190 1112 continue
3191c
3192c     end test code
3193c
3194      do 2 node=np-1,0,-1
3195      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
3196      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
3197      nwmn=ipl(1,2)
3198      if(nwmn.gt.0) then
3199      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
3200      if(npackw.eq.0) then
3201      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm)
3202      else
3203      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm)
3204      call sp_unpackw(nwmn,iwl,iwlp)
3205      endif
3206      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
3207      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
3208      if(lvw)
3209     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
3210      if(lfw)
3211     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm)
3212      do 3 i=1,nwmn
3213      j=iwl(i,lwgmn)-number
3214      if(j.gt.0.and.j.le.numw) then
3215      do 4 k=1,nwa
3216      bxw(j,1,k)=xw(i,1,k)
3217      bxw(j,2,k)=xw(i,2,k)
3218      bxw(j,3,k)=xw(i,3,k)
3219    4 continue
3220      if(lvw) then
3221      do 5 k=1,nwa
3222      bvw(j,1,k)=vw(i,1,k)
3223      bvw(j,2,k)=vw(i,2,k)
3224      bvw(j,3,k)=vw(i,3,k)
3225    5 continue
3226      endif
3227      if(lfw) then
3228      do 51 k=1,nwa
3229      bfw(j,1,k)=fw(i,1,k)
3230      bfw(j,2,k)=fw(i,2,k)
3231      bfw(j,3,k)=fw(i,3,k)
3232   51 continue
3233      endif
3234      ibw(j)=iwl(i,lwdyn)
3235      endif
3236    3 continue
3237      endif
3238    2 continue
3239c     begin solvent output
3240c     low precision
3241      if(nprec.eq.0) then
3242      if(lfw) then
3243      if(lvw) then
3244      do 61 i=1,numw
3245      write(lfntrj,1014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),
3246     + (bfw(i,j,k),j=1,3),k=1,nwa)
3247 1014 format(6f8.3,3f8.1)
3248      if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i)
3249   61 continue
3250      else
3251      do 71 i=1,numw
3252      write(lfntrj,1015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa)
3253 1015 format(3f8.3,3f8.1)
3254      if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i)
3255   71 continue
3256      endif
3257      else
3258      if(lvw) then
3259      do 6 i=1,numw
3260      write(lfntrj,1004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
3261 1004 format(6f8.3)
3262      if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i)
3263    6 continue
3264      else
3265      do 7 i=1,numw
3266      write(lfntrj,1005) ((bxw(i,j,k),j=1,3),k=1,nwa)
3267 1005 format(3f8.3)
3268      if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i)
3269    7 continue
3270      endif
3271      endif
3272c     high precision
3273      else
3274      if(lfw) then
3275      if(lvw) then
3276      do 261 i=1,numw
3277      write(lfntrj,2014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),
3278     + (bfw(i,j,k),j=1,3),k=1,nwa)
3279 2014 format(6e12.6,/,3e12.6)
3280      if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i)
3281  261 continue
3282      else
3283      do 271 i=1,numw
3284      write(lfntrj,2015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa)
3285 2015 format(3e12.6,/,3e12.6)
3286      if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i)
3287  271 continue
3288      endif
3289      else
3290      if(lvw) then
3291      do 26 i=1,numw
3292      write(lfntrj,2004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
3293 2004 format(6e12.6)
3294      if(ibw(i).lt.0) call md_abort('Missing solvent in wttrj',i)
3295   26 continue
3296      else
3297      do 27 i=1,numw
3298      write(lfntrj,2005) ((bxw(i,j,k),j=1,3),k=1,nwa)
3299 2005 format(3e12.6)
3300      if(ibw(i).lt.0) call md_abort('Missing solvent in wtrst',i)
3301   27 continue
3302      endif
3303      endif
3304      endif
3305c     end of solvent output
3306      number=number+numw
3307    1 continue
3308      endif
3309c
3310      if((lxs.or.lvs.or.lfs).and.nsa.gt.0) then
3311      number=0
3312      ncyc=nsa/nb+1
3313      nums=nb
3314      do 8 icyc=1,ncyc
3315      if(nsa-number.lt.nums) nums=nsa-number
3316c
3317c     begin test code 10/31/2001
3318c     initialize ibw to check that all atoms have been received
3319c
3320      do 1117 i=1,nb
3321      ibs(i)=-1
3322 1117 continue
3323c
3324c     end test code
3325c
3326      do 9 node=np-1,0,-1
3327      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
3328      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
3329      nsan=ipl(2,2)
3330      if(nsan.gt.0) then
3331      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
3332      if(npack.eq.0) then
3333      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa)
3334      else
3335      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
3336      call sp_unpack(nsan,isl,islp)
3337      endif
3338      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
3339      call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
3340      if(lvs) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
3341      if(lfs) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
3342      nonh=0
3343      do 10 i=1,nsan
3344      j=isl(i,lsgan)-number
3345      if(j.gt.0.and.j.le.nums) then
3346      bxs(j,1)=xs(i,1)
3347      bxs(j,2)=xs(i,2)
3348      bxs(j,3)=xs(i,3)
3349      if(isl(i,lshop).eq.0) nonh=i
3350      if(iand(isl(i,lshop),1).eq.1) then
3351      if(nonh.gt.0) then
3352      bxs(j,1)=xs(nonh,1)-xs(i,1)
3353      bxs(j,2)=xs(nonh,2)-xs(i,2)
3354      bxs(j,3)=xs(nonh,3)-xs(i,3)
3355      dist=sqrt(bxs(j,1)*bxs(j,1)+bxs(j,2)*bxs(j,2)+bxs(j,3)*bxs(j,3))
3356      dist=dumdst/dist
3357      bxs(j,1)=xs(nonh,1)-dist*bxs(j,1)
3358      bxs(j,2)=xs(nonh,2)-dist*bxs(j,2)
3359      bxs(j,3)=xs(nonh,3)-dist*bxs(j,3)
3360      endif
3361      endif
3362      if(lvs) then
3363      bvs(j,1)=vs(i,1)
3364      bvs(j,2)=vs(i,2)
3365      bvs(j,3)=vs(i,3)
3366      endif
3367      if(lfs) then
3368      bfs(j,1)=fs(i,1)
3369      bfs(j,2)=fs(i,2)
3370      bfs(j,3)=fs(i,3)
3371      endif
3372      ibs(j)=isl(i,lsdyn)
3373      endif
3374   10 continue
3375      endif
3376    9 continue
3377c     begin solute output
3378c     low precision
3379      if(nprec.eq.0) then
3380      if(lfs) then
3381      if(lvs) then
3382      do 111 i=1,nums
3383      write(lfntrj,1016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3),
3384     + (bfs(i,j),j=1,3)
3385 1016 format(6f8.3,3f8.1)
3386      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3387  111 continue
3388      else
3389      do 121 i=1,nums
3390      write(lfntrj,1017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3)
3391 1017 format(3f8.3,3f8.1)
3392      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3393  121 continue
3394      endif
3395      else
3396      if(lvs) then
3397      do 11 i=1,nums
3398      write(lfntrj,1006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3)
3399 1006 format(6f8.3)
3400      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3401   11 continue
3402      else
3403      do 12 i=1,nums
3404      write(lfntrj,1007) (bxs(i,j),j=1,3)
3405 1007 format(3f8.3)
3406      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3407   12 continue
3408      endif
3409      endif
3410      else
3411c     high precision
3412      if(lfs) then
3413      if(lvs) then
3414      do 2111 i=1,nums
3415      write(lfntrj,2016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3),
3416     + (bfs(i,j),j=1,3)
3417 2016 format(6e12.6,/,3e12.6)
3418      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3419 2111 continue
3420      else
3421      do 2121 i=1,nums
3422      write(lfntrj,2017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3)
3423 2017 format(3e12.6,/,3e12.6)
3424      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3425 2121 continue
3426      endif
3427      else
3428      if(lvs) then
3429      do 211 i=1,nums
3430      write(lfntrj,2006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3)
3431 2006 format(6e12.6)
3432      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3433  211 continue
3434      else
3435      do 212 i=1,nums
3436      write(lfntrj,2007) (bxs(i,j),j=1,3)
3437 2007 format(3e12.6)
3438      if(ibs(i).lt.0) call md_abort('Missing solute atom in wtrst',i)
3439  212 continue
3440      endif
3441      endif
3442      end if
3443c     end solute output
3444      number=number+nums
3445    8 continue
3446      endif
3447c
3448      endif
3449c
3450      call util_flush(lfntrj)
3451c
3452      return
3453      end
3454      logical function sp_skip(lfntri,nskip)
3455c
3456      implicit none
3457c
3458#include "sp_common.fh"
3459c
3460      integer lfntri,nskip
3461c
3462      integer i
3463      character*80 card
3464c
3465      if(me.eq.0) then
3466      do 1 i=1,nskip
3467    2 continue
3468      read(lfntri,1000,end=9999) card
3469 1000 format(a)
3470      if(card(1:5).ne.'frame') goto 2
3471    1 continue
3472      endif
3473c
3474      sp_skip=.true.
3475c
3476      return
3477c
3478 9999 continue
3479c
3480      sp_skip=.false.
3481c
3482      return
3483      end
3484      logical function sp_rdtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs,
3485     + stime,pres,temp,tempw,temps,iwl,xw,vw,fw,xwcr,isl,xs,vs,fs,
3486     + numwm,numsa)
3487c
3488      implicit none
3489c
3490#include "sp_common.fh"
3491#include "mafdecls.fh"
3492#include "global.fh"
3493c
3494      integer lfntri
3495      logical lxw,lvw,lfw,lxs,lvs,lfs
3496      integer numwm,numsa,iwl(mwm,miw2),isl(msa,mis2)
3497      real*8 stime,pres,temp,tempw,temps
3498      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
3499      real*8 xs(msa,3),vs(msa,3),fs(msa,3)
3500c
3501      integer lenscr
3502      logical done
3503c
3504      external sp_rtrj
3505      logical sp_rtrj
3506c
3507      lenscr=ma_inquire_avail(mt_byte)/
3508     + ((9*mwa+3)*ma_sizeof(mt_dbl,1,mt_byte)+
3509     + (mis2+4)*ma_sizeof(mt_int,1,mt_byte))-1
3510      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bx',l_bx,i_bx))
3511     + call md_abort('Failed to allocate bx',0)
3512      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bv',l_bv,i_bv))
3513     + call md_abort('Failed to allocate bv',0)
3514      if(.not.ma_push_get(mt_dbl,lenscr*3*mwa,'bf',l_bf,i_bf))
3515     + call md_abort('Failed to allocate bf',0)
3516      if(.not.ma_push_get(mt_dbl,lenscr*3,'br',l_br,i_br))
3517     + call md_abort('Failed to allocate br',0)
3518      if(.not.ma_push_get(mt_int,lenscr*max(mis2,2),'bi',l_bi,i_bi))
3519     + call md_abort('Failed to allocate bi',0)
3520      if(.not.ma_push_get(mt_int,lenscr,'n',l_n,i_n))
3521     + call md_abort('Failed to allocate n',0)
3522c
3523      done=sp_rtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs,
3524     + stime,pres,temp,tempw,temps,
3525     + iwl,int_mb(i_packw),xw,vw,fw,xwcr,isl,int_mb(i_pack),xs,vs,fs,
3526     + int_mb(i_ipl),lenscr,int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
3527     + dbl_mb(i_bf),dbl_mb(i_br),int_mb(i_bi),dbl_mb(i_bx),dbl_mb(i_bv),
3528     + dbl_mb(i_bf))
3529c
3530      if(.not.ma_pop_stack(l_n))
3531     + call md_abort('Failed to deallocate n',0)
3532      if(.not.ma_pop_stack(l_bi))
3533     + call md_abort('Failed to deallocate bi',0)
3534      if(.not.ma_pop_stack(l_br))
3535     + call md_abort('Failed to deallocate br',0)
3536      if(.not.ma_pop_stack(l_bf))
3537     + call md_abort('Failed to deallocate bf',0)
3538      if(.not.ma_pop_stack(l_bv))
3539     + call md_abort('Failed to deallocate bv',0)
3540      if(.not.ma_pop_stack(l_bx))
3541     + call md_abort('Failed to deallocate bx',0)
3542c
3543      call ga_sync()
3544c
3545      call sp_gagetixv(me,iwl,int_mb(i_packw),xw,xwcr,vw,numwm,
3546     + isl,int_mb(i_pack),xs,vs,numsa,int_mb(i_ipl))
3547c
3548      sp_rdtrj=done
3549c
3550      return
3551      end
3552      subroutine sp_gettp(temp,pres)
3553c
3554      implicit none
3555c
3556#include "sp_common.fh"
3557c
3558      real*8 temp,pres
3559c
3560      temp=tmp
3561      pres=prs
3562c
3563      return
3564      end
3565c
3566      logical function sp_rtrj(lfntri,lxw,lvw,lfw,lxs,lvs,lfs,
3567     + stime,pres,temp,tempw,temps,iwl,iwlp,xw,vw,fw,xwcr,
3568     + isl,islp,xs,vs,fs,ipl,nb,ibw,bxw,bvw,bfw,brw,ibs,bxs,bvs,bfs)
3569c
3570      implicit none
3571c
3572#include "sp_common.fh"
3573#include "mafdecls.fh"
3574#include "global.fh"
3575#include "msgids.fh"
3576c
3577      integer lfntri,nb
3578      logical lxw,lvw,lfw,lxs,lvs,lfs
3579      real*8 stime,pres,temp,tempw,temps
3580      integer iwl(mwm,miw2),isl(msa,mis2)
3581      integer iwlp(mwm,npackw),islp(msa,npack)
3582      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),fw(mwm,3,mwa),xwcr(mwm,3)
3583      real*8 xs(msa,3),vs(msa,3),fs(msa,3)
3584      integer ipl(mbox,mip2),ibw(nb),ibs(nb)
3585      real*8 bxw(nb,3,mwa),bvw(nb,3,mwa),bfw(nb,3,mwa),brw(nb,3)
3586      real*8 bxs(nb,3),bvs(nb,3),bfs(nb,3)
3587c
3588      integer i,j,k,node,ncyc,icyc,numw,nums,number,nwmn,nsan
3589      integer ilp,ihp,jlp,jhp,ili,ihi,jli,jhi,ilw,ihw,jlw,jhw
3590      integer ils,ihs,jls,jhs
3591      character*10 rdate,rtime
3592      character*80 card
3593c
3594      logical lpw,lps
3595c
3596      lpw=.false.
3597      lps=.false.
3598c
3599      if(me.eq.0) then
3600c
3601  100 continue
3602      read(lfntri,1000,end=9999) card
3603 1000 format(a)
3604      if(card(1:5).ne.'frame') goto 100
3605c
3606      read(lfntri,1001) stime,temp,pres,rdate,rtime
3607 1001 format(2f12.6,1pe12.5,1x,2a10)
3608      tmp=temp
3609      prs=pres
3610      read(lfntri,1002) ((vlat(i,j),j=1,3),i=1,3)
3611 1002 format(3f12.6)
3612      read(lfntri,1000,end=9999) card
3613c
3614      if(card(8:8).eq.'T'.or.card(8:8).eq.'F') then
3615      read(card,1003) lxw,lvw,lfw,lpw,lxs,lvs,lfs,lps,nwm,nwa,nsa
3616 1003 format(8l1,3i10)
3617      elseif(card(6:6).eq.'T'.or.card(6:6).eq.'F') then
3618      read(card,1023) lxw,lvw,lfw,lxs,lvs,lfs,nwm,nwa,nsa
3619 1023 format(6l1,3i10)
3620      lpw=.false.
3621      lps=.false.
3622      else
3623      read(card,1033) lxw,lvw,lxs,lvs,nwm,nwa,nsa
3624 1033 format(4l1,3i10)
3625      lpw=.false.
3626      lps=.false.
3627      lfw=.false.
3628      lfs=.false.
3629      endif
3630c
3631      if((lxw.or.lvw.or.lfw).and.nwm.gt.0) then
3632      number=0
3633      ncyc=nwm/nb+1
3634      numw=nb
3635      do 1 icyc=1,ncyc
3636      if(nwm-number.lt.numw) numw=nwm-number
3637c
3638      if(nprec.eq.0) then
3639      if(lfw) then
3640      if(lvw) then
3641      do 61 i=1,numw
3642      read(lfntri,1014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),
3643     + (bfw(i,j,k),j=1,3),k=1,nwa)
3644 1014 format(6f8.3,3f8.1)
3645   61 continue
3646      else
3647      do 71 i=1,numw
3648      read(lfntri,1015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa)
3649 1015 format(3f8.3,3f8.1)
3650   71 continue
3651      endif
3652      else
3653      if(lvw) then
3654      do 6 i=1,numw
3655      read(lfntri,1004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
3656 1004 format(6f8.3)
3657    6 continue
3658      else
3659      do 7 i=1,numw
3660      read(lfntri,1005) ((bxw(i,j,k),j=1,3),k=1,nwa)
3661 1005 format(3f8.3)
3662    7 continue
3663      endif
3664      endif
3665      else
3666      if(lfw) then
3667      if(lvw) then
3668      do 261 i=1,numw
3669      read(lfntri,2014) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),
3670     + (bfw(i,j,k),j=1,3),k=1,nwa)
3671 2014 format(6e12.6,/,3e12.6)
3672  261 continue
3673      else
3674      do 271 i=1,numw
3675      read(lfntri,2015) ((bxw(i,j,k),j=1,3),(bfw(i,j,k),j=1,3),k=1,nwa)
3676 2015 format(3e12.6,/,3e12.6)
3677  271 continue
3678      endif
3679      else
3680      if(lvw) then
3681      do 26 i=1,numw
3682      read(lfntri,2004) ((bxw(i,j,k),j=1,3),(bvw(i,j,k),j=1,3),k=1,nwa)
3683 2004 format(6e12.6)
3684   26 continue
3685      else
3686      do 27 i=1,numw
3687      read(lfntri,2005) ((bxw(i,j,k),j=1,3),k=1,nwa)
3688 2005 format(3e12.6)
3689   27 continue
3690      endif
3691      endif
3692      endif
3693c
3694      do 2 node=np-1,0,-1
3695      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
3696      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
3697      nwmn=ipl(1,2)
3698      if(nwmn.gt.0) then
3699      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
3700      if(npackw.eq.0) then
3701      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+lwdyn-1,iwl,mwm)
3702      else
3703      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm)
3704      call sp_unpackw(nwmn,iwl,iwlp)
3705      endif
3706      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
3707      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
3708      if(lvw)
3709     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
3710      if(lfw)
3711     + call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm)
3712      do 3 i=1,nwmn
3713      j=iwl(i,lwgmn)-number
3714      if(j.gt.0.and.j.le.numw) then
3715      do 4 k=1,nwa
3716      xw(i,1,k)=bxw(j,1,k)
3717      xw(i,2,k)=bxw(j,2,k)
3718      xw(i,3,k)=bxw(j,3,k)
3719    4 continue
3720      if(lvw) then
3721      do 5 k=1,nwa
3722      vw(i,1,k)=bvw(j,1,k)
3723      vw(i,2,k)=bvw(j,2,k)
3724      vw(i,3,k)=bvw(j,3,k)
3725    5 continue
3726      endif
3727      if(lfw) then
3728      do 51 k=1,nwa
3729      fw(i,1,k)=bfw(j,1,k)
3730      fw(i,2,k)=bfw(j,2,k)
3731      fw(i,3,k)=bfw(j,3,k)
3732   51 continue
3733      endif
3734      endif
3735    3 continue
3736      call ga_put(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
3737      if(lvw)
3738     + call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
3739      if(lfw)
3740     + call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+6*mwa+3,jlw+9*mwa+2,fw,mwm)
3741      endif
3742    2 continue
3743      number=number+numw
3744    1 continue
3745      endif
3746c
3747      if((lxs.or.lvs.or.lfs).and.nsa.gt.0) then
3748      number=0
3749      ncyc=nsa/nb+1
3750      nums=nb
3751      do 8 icyc=1,ncyc
3752      if(nsa-number.lt.nums) nums=nsa-number
3753c
3754      if(nprec.eq.0) then
3755      if(lfs) then
3756      if(lvs) then
3757      do 111 i=1,nums
3758      read(lfntri,1016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3),
3759     + (bfs(i,j),j=1,3)
3760 1016 format(6f8.3,3f8.1)
3761  111 continue
3762      else
3763      do 121 i=1,nums
3764      read(lfntri,1017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3)
3765 1017 format(3f8.3,3f8.1)
3766  121 continue
3767      endif
3768      else
3769      if(lvs) then
3770      do 11 i=1,nums
3771      read(lfntri,1006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3)
3772 1006 format(6f8.3)
3773   11 continue
3774      else
3775      do 12 i=1,nums
3776      read(lfntri,1007) (bxs(i,j),j=1,3)
3777 1007 format(3f8.3)
3778   12 continue
3779      endif
3780      endif
3781      else
3782      if(lfs) then
3783      if(lvs) then
3784      do 2111 i=1,nums
3785      read(lfntri,2016) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3),
3786     + (bfs(i,j),j=1,3)
3787 2016 format(6e12.6,/,3e12.6)
3788 2111 continue
3789      else
3790      do 2121 i=1,nums
3791      read(lfntri,2017) (bxs(i,j),j=1,3),(bfs(i,j),j=1,3)
3792 2017 format(3e12.6,/,3e12.6)
3793 2121 continue
3794      endif
3795      else
3796      if(lvs) then
3797      do 211 i=1,nums
3798      read(lfntri,2006) (bxs(i,j),j=1,3),(bvs(i,j),j=1,3)
3799 2006 format(6e12.6)
3800  211 continue
3801      else
3802      do 212 i=1,nums
3803      read(lfntri,2007) (bxs(i,j),j=1,3)
3804 2007 format(3e12.6)
3805  212 continue
3806      endif
3807      endif
3808      endif
3809c
3810      do 9 node=np-1,0,-1
3811      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
3812      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
3813      nsan=ipl(2,2)
3814      if(nsan.gt.0) then
3815      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
3816      if(npack.eq.0) then
3817      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+lsdyn-1,isl,msa)
3818      else
3819      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
3820      call sp_unpack(nsan,isl,islp)
3821      endif
3822      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
3823      call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
3824      if(lvs) call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
3825      if(lfs) call ga_get(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
3826      do 10 i=1,nsan
3827      j=isl(i,lsgan)-number
3828      if(j.gt.0.and.j.le.nums) then
3829      xs(i,1)=bxs(j,1)
3830      xs(i,2)=bxs(j,2)
3831      xs(i,3)=bxs(j,3)
3832      if(lvs) then
3833      vs(i,1)=bvs(j,1)
3834      vs(i,2)=bvs(j,2)
3835      vs(i,3)=bvs(j,3)
3836      endif
3837      if(lfs) then
3838      fs(i,1)=bfs(j,1)
3839      fs(i,2)=bfs(j,2)
3840      fs(i,3)=bfs(j,3)
3841      endif
3842      endif
3843   10 continue
3844      call ga_put(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
3845      if(lvs) call ga_put(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
3846      if(lfs) call ga_put(ga_s,ils,ils+nsan-1,jls+6,jls+8,fs,msa)
3847      endif
3848    9 continue
3849      number=number+nums
3850    8 continue
3851      endif
3852c
3853      endif
3854c
3855      call ga_brdcst(msp_15,vlat,
3856     + 9*ma_sizeof(mt_dbl,1,mt_byte),0)
3857c
3858      box(1)=vlat(1,1)
3859      box(2)=vlat(2,2)
3860      box(3)=vlat(3,3)
3861      boxh(1)=5.0d-1*box(1)
3862      boxh(2)=5.0d-1*box(2)
3863      boxh(3)=5.0d-1*box(3)
3864c
3865      sp_rtrj=.true.
3866c
3867      return
3868c
3869 9999 continue
3870c
3871      sp_rtrj=.false.
3872c
3873      return
3874c
3875      end
3876c
3877      subroutine sp_wrtmro(lfnmro,stime,pres,temp,tempw,temps,
3878     + iwl,xw,vw,xwcr,isl,xs,vs,prjct)
3879c
3880      implicit none
3881c
3882#include "sp_common.fh"
3883#include "mafdecls.fh"
3884#include "global.fh"
3885c
3886      integer lfnmro
3887      integer iwl(mwm,miw2),isl(msa,mis2)
3888      real*8 stime,pres,temp,tempw,temps
3889      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3)
3890      real*8 xs(msa,3),vs(msa,3)
3891      character*80 prjct
3892c
3893      project=prjct
3894c
3895      call sp_wtmro(lfnmro,stime,pres,temp,tempw,temps,
3896     + iwl,int_mb(i_packw),xw,vw,xwcr,isl,int_mb(i_pack),xs,vs,
3897     + int_mb(i_ipl))
3898c
3899      return
3900      end
3901      subroutine sp_wtmro(lfnmro,stime,pres,temp,tempw,temps,
3902     + iwl,iwlp,xw,vw,xwcr,isl,islp,xs,vs,ipl)
3903c
3904      implicit none
3905c
3906#include "sp_common.fh"
3907#include "mafdecls.fh"
3908#include "global.fh"
3909c
3910      integer lfnmro
3911      real*8 stime,pres,temp,tempw,temps
3912      integer iwl(mwm,miw2),isl(msa,mis2)
3913      integer iwlp(mwm,npackw),islp(msa,npack)
3914      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3)
3915      real*8 xs(msa,3),vs(msa,3)
3916      integer ipl(mbox,mip2)
3917c
3918      integer j,k,l,nwmn,nsan,node,ilp,ihp,jlp,jhp,ili,ihi,jli,jhi
3919      integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs
3920      character*10 rdate,rtime
3921      character*18 user
3922#ifdef USE_POSIXF
3923      integer ilen,ierror
3924#endif
3925c
3926      write(lfnmro) nwm,nwa,nsa,stime,temp,pres,vlat,nhist
3927c
3928      call swatch(rdate,rtime)
3929#ifdef USE_POSIXF
3930      call pxfgetlogin(user, ilen, ierror)
3931#else
3932      call getlog(user)
3933#endif
3934      if(user(18:18).ne.' ') user='                  '
3935      hist(nhist)(1:18)=user
3936      hist(nhist)(19:28)=rdate
3937      hist(nhist)(29:48)=rtime
3938      hist(nhist)(49:108)=project(1:60)
3939      write(lfnmro) (hist(j),j=1,nhist)
3940c
3941      do 1 node=np-1,0,-1
3942      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
3943      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
3944      write(lfnmro) ((ipl(j,k),j=1,mbox),k=1,mip2)
3945      nwmn=ipl(1,2)
3946      nsan=ipl(2,2)
3947      if(nwmn.gt.0) then
3948      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
3949      if(npackw.eq.0) then
3950      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jhi,iwl,mwm)
3951      else
3952      call ga_get(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm)
3953      call sp_unpackw(nwmn,iwl,iwlp)
3954      endif
3955      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
3956      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
3957      call ga_get(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
3958      write(lfnmro) ((iwl(j,k),j=1,nwmn),k=1,miw2)
3959      write(lfnmro) (((xw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa)
3960      write(lfnmro) (((vw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa)
3961      endif
3962      if(nsan.gt.0) then
3963      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
3964      if(npack.eq.0) then
3965      call ga_get(ga_is,ili,ili+nsan-1,jli,jhi,isl,msa)
3966      else
3967      call ga_get(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
3968      call sp_unpack(nsan,isl,islp)
3969      endif
3970      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
3971      call ga_get(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
3972      call ga_get(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
3973      write(lfnmro) ((isl(j,k),j=1,nsan),k=1,mis2)
3974      write(lfnmro) ((xs(j,k),j=1,nsan),k=1,3)
3975      write(lfnmro) ((vs(j,k),j=1,nsan),k=1,3)
3976      endif
3977    1 continue
3978c
3979      return
3980      end
3981      logical function sp_rdmri(lfnmri,stime,pres,temp,tempw,temps,
3982     + iwl,xw,vw,xwcr,isl,xs,vs)
3983c
3984      implicit none
3985c
3986#include "sp_common.fh"
3987#include "mafdecls.fh"
3988#include "global.fh"
3989c
3990      logical sp_rmri
3991      external sp_rmri
3992c
3993      integer lfnmri
3994      integer iwl(mwm,miw2),isl(msa,mis2)
3995      real*8 stime,pres,temp,tempw,temps
3996      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3)
3997      real*8 xs(msa,3),vs(msa,3)
3998c
3999      sp_rdmri=sp_rmri(lfnmri,stime,pres,temp,tempw,temps,
4000     + iwl,int_mb(i_packw),xw,vw,xwcr,isl,int_mb(i_pack),xs,vs,
4001     + int_mb(i_ipl))
4002c
4003      return
4004      end
4005      logical function sp_rmri(lfnmri,stime,pres,temp,tempw,temps,
4006     + iwl,iwlp,xw,vw,xwcr,isl,islp,xs,vs,ipl)
4007c
4008      implicit none
4009c
4010#include "sp_common.fh"
4011#include "mafdecls.fh"
4012#include "msgids.fh"
4013#include "global.fh"
4014c
4015      integer lfnmri
4016      real*8 stime,pres,temp,tempw,temps
4017      integer iwl(mwm,miw2),isl(msa,mis2)
4018      integer iwlp(mwm,npackw),islp(msa,npack)
4019      real*8 xw(mwm,3,mwa),vw(mwm,3,mwa),xwcr(mwm,3)
4020      real*8 xs(msa,3),vs(msa,3)
4021      integer ipl(mbox,mip2)
4022c
4023      integer i,j,k,l,nwmn,nsan,node,ilp,ihp,jlp,jhp,ili,ihi,jli,jhi
4024      integer ilw,ihw,jlw,jhw,ils,ihs,jls,jhs
4025      integer ltemp(3)
4026      real*8 rtemp(12)
4027c
4028      if(me.eq.0) then
4029      read(lfnmri,err=9,end=9) ltemp,rtemp,nhist
4030      if(nhist.gt.0) read(lfnmri,err=9,end=9) (hist(j),j=1,nhist)
4031      do 1 node=np-1,0,-1
4032      read(lfnmri,err=9) ((ipl(j,k),j=1,mbox),k=1,mip2)
4033      nwmn=ipl(1,2)
4034      nsan=ipl(2,2)
4035      call ga_distribution(ga_ip,node,ilp,ihp,jlp,jhp)
4036      call ga_put(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
4037      if(nwmn.gt.0) then
4038      read(lfnmri,err=9) ((iwl(j,k),j=1,nwmn),k=1,miw2)
4039      read(lfnmri,err=9) (((xw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa)
4040      read(lfnmri,err=9) (((vw(j,k,l),j=1,nwmn),k=1,3),l=1,nwa)
4041      call ga_distribution(ga_iw,node,ili,ihi,jli,jhi)
4042      if(npackw.eq.0) then
4043      call ga_put(ga_iw,ili,ili+nwmn-1,jli,jhi,iwl,mwm)
4044      else
4045      call sp_packw(nwmn,iwl,iwlp)
4046      call ga_put(ga_iw,ili,ili+nwmn-1,jli,jli+npackw-1,iwlp,mwm)
4047      endif
4048      call ga_distribution(ga_w,node,ilw,ihw,jlw,jhw)
4049      call ga_put(ga_w,ilw,ilw+nwmn-1,jlw,jlw+3*mwa-1,xw,mwm)
4050      call ga_put(ga_w,ilw,ilw+nwmn-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
4051      endif
4052      if(nsan.gt.0) then
4053      read(lfnmri,err=9) ((isl(j,k),j=1,nsan),k=1,mis2)
4054      read(lfnmri,err=9) ((xs(j,k),j=1,nsan),k=1,3)
4055      read(lfnmri,err=9) ((vs(j,k),j=1,nsan),k=1,3)
4056      call ga_distribution(ga_is,node,ili,ihi,jli,jhi)
4057      if(npack.eq.0) then
4058      call ga_put(ga_is,ili,ili+nsan-1,jli,jhi,isl,msa)
4059      else
4060      call sp_pack(nsan,isl,islp)
4061      call ga_put(ga_is,ili,ili+nsan-1,jli,jli+npack-1,islp,msa)
4062      endif
4063      call ga_distribution(ga_s,node,ils,ihs,jls,jhs)
4064      call ga_put(ga_s,ils,ils+nsan-1,jls,jls+2,xs,msa)
4065      call ga_put(ga_s,ils,ils+nsan-1,jls+3,jls+5,vs,msa)
4066      endif
4067    1 continue
4068      endif
4069c
4070      call ga_brdcst(msp_12,ltemp,3*ma_sizeof(mt_int,1,mt_byte),0)
4071      call ga_brdcst(msp_13,rtemp,12*ma_sizeof(mt_dbl,1,mt_byte),0)
4072      nwm=ltemp(1)
4073      nwa=ltemp(2)
4074      nsa=ltemp(3)
4075      stime=rtemp(1)
4076      temp=rtemp(2)
4077      pres=rtemp(3)
4078      k=3
4079      do 2 i=1,3
4080      do 3 j=1,3
4081      k=k+1
4082      vlat(i,j)=rtemp(k)
4083    3 continue
4084    2 continue
4085      call ga_distribution(ga_ip,me,ilp,ihp,jlp,jhp)
4086      call ga_get(ga_ip,ilp,ihp,jlp,jhp,ipl,mbox)
4087      nwmloc=ipl(1,2)
4088      nsaloc=ipl(2,2)
4089      if(nwmloc.gt.0) then
4090      call ga_distribution(ga_iw,me,ili,ihi,jli,jhi)
4091      if(npackw.eq.0) then
4092      call ga_get(ga_iw,ili,ili+nwmloc-1,jli,jhi,iwl,mwm)
4093      else
4094      call ga_get(ga_iw,ili,ili+nwmloc-1,jli,jli+npackw-1,iwlp,mwm)
4095      call sp_unpackw(nwmloc,iwl,iwlp)
4096      endif
4097      call ga_distribution(ga_w,me,ilw,ihw,jlw,jhw)
4098      call ga_get(ga_w,ilw,ilw+nwmloc-1,jlw,jlw+3*mwa-1,xw,mwm)
4099      call ga_get(ga_w,ilw,ilw+nwmloc-1,jlw+3*mwa,jlw+6*mwa-1,vw,mwm)
4100      endif
4101      if(nsaloc.gt.0) then
4102      call ga_distribution(ga_is,me,ili,ihi,jli,jhi)
4103      if(npack.eq.0) then
4104      call ga_get(ga_is,ili,ili+nsaloc-1,jli,jhi,isl,msa)
4105      else
4106      call ga_get(ga_is,ili,ili+nsaloc-1,jli,jli+npack-1,islp,msa)
4107      call sp_unpack(nsaloc,isl,islp)
4108      endif
4109      call ga_distribution(ga_s,me,ils,ihs,jls,jhs)
4110      call ga_get(ga_s,ils,ils+nsaloc-1,jls,jls+2,xs,msa)
4111      call ga_get(ga_s,ils,ils+nsaloc-1,jls+3,jls+5,vs,msa)
4112      endif
4113c
4114      sp_rmri=.true.
4115      return
4116c
4117    9 continue
4118      sp_rmri=.false.
4119      return
4120      end
4121      subroutine sp_fix()
4122c
4123      implicit none
4124c
4125      return
4126      end
4127