1      subroutine smd_vlist_init_system()
2      implicit none
3#include "errquit.fh"
4#include "inp.fh"
5#include "mafdecls.fh"
6#include "rtdb.fh"
7#include "util.fh"
8#include "global.fh"
9c
10      character*32 sp_vlist,sp_exlist,sp_coords
11      character*32 tag,pname
12      logical result
13      logical oexlist
14
15      pname = "smd_vlist_init_system"
16c
17      tag = "coordinates"
18      call smd_system_get_component(sp_coords,tag,result)
19      if(.not.result)
20     >  call errquit(
21     >       pname//'no component '//tag,0,0)
22
23      oexlist = .true.
24      tag = "excl_list"
25      call smd_system_get_component(sp_exlist,tag,result)
26      if(.not.result) oexlist = .false.
27
28      tag = "verlet_list"
29      call smd_system_get_component(sp_vlist,tag,result)
30      if(.not.result)
31     >  call errquit(
32     >       pname//'no component '//tag,0,0)
33
34
35      call smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coords)
36c
37      return
38      end
39
40      subroutine smd_vlist_init(oexlist,sp_vlist,sp_exlist,sp_coord)
41      implicit none
42#include "errquit.fh"
43#include "inp.fh"
44#include "mafdecls.fh"
45#include "rtdb.fh"
46#include "util.fh"
47#include "global.fh"
48c
49      character*(*) sp_vlist
50      character*(*) sp_exlist
51      character*(*) sp_coord
52      logical oexlist
53c
54      character*32 pname
55      character*80 tag
56      integer np,nl
57      integer i_p,i_l,i_c
58      integer h_l
59      integer i_list,i_clist
60      integer i
61      integer i_xp,i_xl
62      integer i_cl,h_cl
63      integer i_ct,h_ct
64      double precision rc2
65      integer nlb
66      logical result
67c
68      pname = "smd_vlist_init"
69c
70c      write(*,*) "in "//pname
71c
72c     get coordinate information
73c     --------------------------
74      tag = "coords"
75      call smd_data_get_index(sp_coord,tag,i_c,result)
76      if(.not. result)
77     >  call errquit(
78     >       pname//'error getting index for'//tag,0, RTDB_ERR)
79
80      call smd_data_get_size(sp_coord,tag,np,result)
81      if(.not. result)
82     >  call errquit(
83     >       pname//'error getting size for'//tag,0, RTDB_ERR)
84      np = np/3
85c
86c     get excluded list information
87c     -----------------------------
88      if(oexlist) then
89      tag = "exlist:pointer"
90      call smd_data_get_index(sp_exlist,tag,i_xp,result)
91      if(.not. result)
92     >  call errquit(
93     >       pname//'error getting index for'//tag,0, RTDB_ERR)
94      tag = "exlist:list"
95      call smd_data_get_index(sp_exlist,tag,i_xl,result)
96      if(.not. result)
97     >  call errquit(
98     >       pname//'error getting index for'//tag,0, RTDB_ERR)
99
100      end if
101c
102c     gestimate the size of pair list
103c     ------------------------------
104      nl =  min( 7*np*200, ma_inquire_avail(MT_INT))
105      nl = nl/7
106c
107c     create pointer memory
108c     ---------------------
109      call smd_namespace_create(sp_vlist)
110      tag = "vlist:pointer"
111      call smd_data_create(sp_vlist,"vlist:pointer",np,MT_INT)
112      call smd_data_get_index(sp_vlist,tag,i_p,result)
113      if(.not. result)
114     >  call errquit(
115     >       pname//'error getting index for'//tag,0, RTDB_ERR)
116
117c
118c    create temporary scratch array for list since
119c    we do not know the size yet
120c    ---------------------------------------------
121      if(.not.ma_push_get(mt_int,nl,'tmp l',h_l,i_l))
122     + call errquit(pname//'Failed to allocate memory for tmp l',
123     + nl, MA_ERR)
124
125      if(.not.ma_push_get(mt_dbl,3*nl,'tmp cl',h_cl,i_cl))
126     + call errquit(pname//'Failed to allocate memory for tmp l',
127     + nl, MA_ERR)
128
129      if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct))
130     + call errquit(pname//'Failed to allocate memory for tmp l',
131     + nl, MA_ERR)
132
133
134      call smd_cutoff_get_rcut_verlet(rc2)
135      rc2=rc2*rc2
136
137      if(oexlist) then
138      call smd_vlist_set(np,
139     +                   nl,
140     +                   rc2,
141     +                   dbl_mb(i_c),
142     +                   int_mb(i_xp),
143     +                   int_mb(i_xl),
144     +                   int_mb(i_p),
145     +                   int_mb(i_l),
146     +                   dbl_mb(i_cl),
147     +                   dbl_mb(i_ct),
148     +                   result)
149
150      else
151
152      call smd_vlist_set1(np,
153     +                   nl,
154     +                   rc2,
155     +                   dbl_mb(i_c),
156     +                   int_mb(i_p),
157     +                   int_mb(i_l),
158     +                   dbl_mb(i_cl),
159     +                   dbl_mb(i_ct),
160     +                   result)
161
162
163      end if
164c
165c     create list memory
166c     nl now contains the actual size
167c     we will buffer it though to allow for possible expansion
168c     --------------------------------------------------------
169      nlb = 500
170      tag = "vlist:list"
171      call smd_data_create(sp_vlist,tag,nl+nlb,MT_INT)
172      call smd_data_get_index(sp_vlist,tag,i_list,result)
173      if(.not. result)
174     >  call errquit(
175     >       pname//'error getting index for'//tag,0, RTDB_ERR)
176
177
178      do i=1,nl
179       int_mb(i_list+i-1) = int_mb(i_l+i-1)
180      end do
181
182      tag = "vlist:distances"
183      call smd_data_create(sp_vlist,tag,3*(nl+nlb),MT_DBL)
184      call smd_data_get_index(sp_vlist,tag,i_clist,result)
185      if(.not. result)
186     >  call errquit(
187     >       pname//'error getting index for'//tag,0, RTDB_ERR)
188
189      tag = "vlist:displacement"
190      call smd_data_create(sp_vlist,tag,3*np,MT_DBL)
191
192      do i=1,3*nl
193       dbl_mb(i_clist+i-1) = dbl_mb(i_cl+i-1)
194      end do
195
196      if(.not.ma_pop_stack(h_ct))
197     & call errquit(pname//'Failed to deallocate stack h_l',nl,
198     &       MA_ERR)
199
200
201      if(.not.ma_pop_stack(h_cl))
202     & call errquit(pname//'Failed to deallocate stack h_l',nl,
203     &       MA_ERR)
204
205
206      if(.not.ma_pop_stack(h_l))
207     & call errquit(pname//'Failed to deallocate stack h_l',nl,
208     &       MA_ERR)
209
210
211      return
212      end
213c
214      subroutine smd_vlist_update(olist,ocoord)
215      implicit none
216#include "errquit.fh"
217#include "inp.fh"
218#include "mafdecls.fh"
219#include "rtdb.fh"
220#include "util.fh"
221#include "global.fh"
222c
223      logical olist
224      logical ocoord
225c
226      character*32 sp_vlist
227      character*32 sp_exlist
228      character*32 sp_coord
229c
230      character*32 pname
231      character*80 tag
232      integer np,nl
233      integer i_p,i_l,i_c
234      integer h_l
235      integer i_list,i_clist
236      integer i
237      integer i_xp,i_xl
238      double precision rc2
239      integer nlb
240      logical result
241      integer i_ct,h_ct
242      logical oexlist
243c
244      pname = "smd_vlist_update"
245c
246c      write(*,*) "in "//pname
247c
248      if((.not.olist).and.(.not.ocoord)) return
249c
250c     get components
251c     --------------
252      tag = "coordinates"
253      call smd_system_get_component(sp_coord,tag,result)
254      if(.not.result)
255     >  call errquit(
256     >       pname//'no component '//tag,0,0)
257
258      oexlist = .true.
259      tag = "excl_list"
260      call smd_system_get_component(sp_exlist,tag,result)
261      if(.not.result) oexlist =.false.
262
263      tag = "verlet_list"
264      call smd_system_get_component(sp_vlist,tag,result)
265      if(.not.result)
266     >  call errquit(
267     >       pname//'no component '//tag,0,0)
268
269c
270      if(.not.olist) then
271        call smd_vlist_update_coord(sp_vlist,sp_coord)
272        goto 200
273      end if
274c
275c     get coordinate information
276c     --------------------------
277      tag = "coords"
278      call smd_data_get_index(sp_coord,tag,i_c,result)
279      if(.not. result)
280     >  call errquit(
281     >       pname//'error getting index for'//tag,0, RTDB_ERR)
282
283      call smd_data_get_size(sp_coord,tag,np,result)
284      if(.not. result)
285     >  call errquit(
286     >       pname//'error getting size for'//tag,0, RTDB_ERR)
287      np = np/3
288c
289c     get excluded list information
290c     -----------------------------
291      if(oexlist) then
292      tag = "exlist:pointer"
293      call smd_data_get_index(sp_exlist,tag,i_xp,result)
294      if(.not. result)
295     >  call errquit(
296     >       pname//'error getting index for'//tag,0, RTDB_ERR)
297      tag = "exlist:list"
298      call smd_data_get_index(sp_exlist,tag,i_xl,result)
299      if(.not. result)
300     >  call errquit(
301     >       pname//'error getting index for'//tag,0, RTDB_ERR)
302
303      end if
304c
305c     get verlet list information
306c     ---------------------------
307      tag = "vlist:pointer"
308      call smd_data_get_index(sp_vlist,tag,i_p,result)
309      if(.not. result)
310     >  call errquit(
311     >       pname//'error getting index for'//tag,0, RTDB_ERR)
312
313      tag = "vlist:list"
314      call smd_data_get_index(sp_vlist,tag,i_list,result)
315      if(.not. result)
316     >  call errquit(
317     >       pname//'error getting index for'//tag,0, RTDB_ERR)
318      call smd_data_get_size(sp_vlist,tag,nl,result)
319      if(.not. result)
320     >  call errquit(
321     >       pname//'error getting size for'//tag,0, RTDB_ERR)
322
323      tag = "vlist:distances"
324      call smd_data_get_index(sp_vlist,tag,i_clist,result)
325      if(.not. result)
326     >  call errquit(
327     >       pname//'error getting index for'//tag,0, RTDB_ERR)
328c
329c    create temporary scratch array for list
330c    ---------------------------------------------
331
332      if(.not.ma_push_get(mt_dbl,3*np,'tmp',h_ct,i_ct))
333     + call errquit(pname//'Failed to allocate memory for tmp l',
334     + nl, MA_ERR)
335
336
337      call smd_cutoff_get_rcut_verlet(rc2)
338      rc2=rc2*rc2
339
340      if(oexlist) then
341      call smd_vlist_set(np,
342     +                   nl,
343     +                   rc2,
344     +                   dbl_mb(i_c),
345     +                   int_mb(i_xp),
346     +                   int_mb(i_xl),
347     +                   int_mb(i_p),
348     +                   int_mb(i_list),
349     +                   dbl_mb(i_clist),
350     +                   dbl_mb(i_ct),
351     +                   result)
352
353      else
354
355      call smd_vlist_set1(np,
356     +                   nl,
357     +                   rc2,
358     +                   dbl_mb(i_c),
359     +                   int_mb(i_p),
360     +                   int_mb(i_list),
361     +                   dbl_mb(i_clist),
362     +                   dbl_mb(i_ct),
363     +                   result)
364
365
366      end if
367
368      if(.not.ma_pop_stack(h_ct))
369     & call errquit(pname//'Failed to deallocate stack h_l',nl,
370     &       MA_ERR)
371
372
373200   continue
374      return
375      end
376c
377      subroutine smd_vlist_update_coord(sp_vlist,sp_coord)
378      implicit none
379#include "errquit.fh"
380#include "inp.fh"
381#include "mafdecls.fh"
382#include "rtdb.fh"
383#include "util.fh"
384#include "global.fh"
385c
386      character*(*) sp_vlist
387      character*(*) sp_coord
388c
389      character*32 pname
390      character*80 tag
391      integer np,nl
392      integer i_p,i_l,i_c
393      integer i_list,i_clist
394      integer i_cl
395      logical result
396c
397      pname = "smd_vlist_init"
398c
399c      write(*,*) "in "//pname
400c
401c     get coordinate information
402c     --------------------------
403      tag = "coords"
404      call smd_data_get_index(sp_coord,tag,i_c,result)
405      if(.not. result)
406     >  call errquit(
407     >       pname//'error getting index for'//tag,0, RTDB_ERR)
408
409      call smd_data_get_size(sp_coord,tag,np,result)
410      if(.not. result)
411     >  call errquit(
412     >       pname//'error getting size for'//tag,0, RTDB_ERR)
413      np = np/3
414c
415c     get list information
416c     --------------------
417      tag = "vlist:pointer"
418      call smd_data_get_index(sp_vlist,tag,i_p,result)
419      if(.not. result)
420     >  call errquit(
421     >       pname//'error getting index for'//tag,0, RTDB_ERR)
422
423      tag = "vlist:list"
424      call smd_data_get_index(sp_vlist,tag,i_list,result)
425      if(.not. result)
426     >  call errquit(
427     >       pname//'error getting index for'//tag,0, RTDB_ERR)
428
429
430      tag = "vlist:distances"
431      call smd_data_get_index(sp_vlist,tag,i_clist,result)
432      if(.not. result)
433     >  call errquit(
434     >       pname//'error getting index for'//tag,0, RTDB_ERR)
435
436      nl = int_mb(i_p+np-1) - int_mb(i_p)
437      call smd_vlist_update_coord0(np,nl,
438     +                   dbl_mb(i_c),
439     +                   int_mb(i_p),
440     +                   int_mb(i_l),
441     +                   dbl_mb(i_cl))
442
443      call smd_lat_rebox(nl,dbl_mb(i_clist))
444
445      return
446      end
447c
448      subroutine smd_vlist_update_coord0(np,nl,c,point,
449     >                         list,cl)
450c
451c     update coordinates in verlet pair list
452c     point(i) is an index to list() array
453c     that contains all the pairs of atom i
454c     In other words all the atoms paired with atom i
455c     are contained in list(point(i)),...,list(point(i+1)-1)
456c     np     [in]  number of atoms (which is also size of pointer array)
457c     nl     [in]  list size
458c     c      [in]  coordinates
459c     point  [in] verlet pointer
460c     list   [in] verlet list
461c     cl     [out] list of vectors (ri-rj)
462      implicit none
463#include "errquit.fh"
464      integer np,nl
465      double precision c(np,3)
466      integer point(np)
467      integer list(nl)
468      double precision cl(nl,3)
469c
470      integer i,j
471      integer nlist
472      integer jnab,jbeg,jend
473
474      character*30 pname
475
476      pname = "smd_vlist_update_coord0"
477
478      nlist=0
479      do i=1,np-1
480       jbeg=point(i)
481       jend=point(i+1)-1
482
483       if(jbeg.le.jend)then
484
485        do jnab=jbeg,jend
486
487         j=list(jnab)
488
489         nlist = nlist + 1
490         if(nlist.gt.nl)
491     >       call errquit(
492     >       pname//'out of bounds',0, RTDB_ERR)
493
494         cl(nlist,1)=c(i,1)-c(j,1)
495         cl(nlist,2)=c(i,2)-c(j,2)
496         cl(nlist,3)=c(i,3)-c(j,3)
497
498        enddo
499
500       end if
501      enddo
502
503
504200   continue
505      return
506
507      END
508c
509      subroutine smd_vlist_set(np,nl,vcutsq,c,xp,xl,point,
510     >                         list,cl,ct,result)
511c
512c     constructs verlet pairt list
513c     point(i) is an index to list() array
514c     that contains all the pairs of atom i
515c     In other words all the atoms paired with atom i
516c     are contained in list(point(i)),...,list(point(i+1)-1)
517c     np     [in]  number of atoms (which is also size of pointer array)
518c     nl     [inout] size of verlet list
519c     c      [in]  coordinates
520c     xp     [in]  excluded pointer
521c     xl     [in]  excluded list
522c     point  [out] verlet pointer
523c     list   [out] verlet list
524c     cl     [out] list of vectors (ri-rj)
525c     result [out] status of subroutine
526      implicit none
527#include "errquit.fh"
528      integer np
529      integer nl
530      double precision vcutsq
531      double precision c(np,3)
532      integer xp(np)
533      integer xl(*)
534      integer point(np)
535      integer list(nl)
536      double precision cl(nl,3)
537      double precision ct(np,3)
538      logical result
539c
540      integer i,j,k
541      integer nlist
542      integer eatm
543      double precision rij(3),rijsq,cc(1,3)
544
545      character*30 pname
546
547      pname = "smd_vlist_set"
548
549      result = .false.
550      nlist=0
551
552c      write(*,*) "vcutsq=",vcutsq
553      do i=1,np-1
554
555       point(i)=nlist+1
556       if(xp(i).ne.xp(i+1))eatm=xp(i)
557
558       k = 0
559       do j=i+1,np
560
561         k = k +1
562         ct(k,1)=c(i,1)-c(j,1)
563         ct(k,2)=c(i,2)-c(j,2)
564         ct(k,3)=c(i,3)-c(j,3)
565
566       end do
567
568
569       call smd_lat_rebox(np,ct)
570
571
572       k = 0
573       do j=i+1,np
574
575        k = k + 1
576        if((xp(i).ne.xp(i+1)).and.(xl(eatm).eq.j))then
577
578         eatm=min(eatm+1,(xp(i+1)-1))
579
580        else
581
582         rij(1)=ct(k,1)
583         rij(2)=ct(k,2)
584         rij(3)=ct(k,3)
585
586         rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)
587
588
589         if(rijsq.lt.vcutsq)then
590
591          nlist=nlist+1
592
593          if(nlist.gt.nl)then
594           result = .false.
595           goto 200
596          endif
597
598          list(nlist)=j
599          cl(nlist,1)=rij(1)
600          cl(nlist,2)=rij(2)
601          cl(nlist,3)=rij(3)
602
603         endif
604
605        endif
606
607       enddo
608
609      enddo
610
611      point(np)=nlist+1
612
613      nl = nlist
614
615      result = .true.
616
617200   continue
618      return
619
620      END
621
622      subroutine smd_vlist_set1(np,nl,vcutsq,c,point,
623     >                         list,cl,ct,result)
624c
625c     constructs verlet pairt list
626c     point(i) is an index to list() array
627c     that contains all the pairs of atom i
628c     In other words all the atoms paired with atom i
629c     are contained in list(point(i)),...,list(point(i+1)-1)
630c     np     [in]  number of atoms (which is also size of pointer array)
631c     nl     [inout] size of verlet list
632c     c      [in]  coordinates
633c     xp     [in]  excluded pointer
634c     xl     [in]  excluded list
635c     point  [out] verlet pointer
636c     list   [out] verlet list
637c     cl     [out] list of vectors (ri-rj)
638c     result [out] status of subroutine
639      implicit none
640#include "errquit.fh"
641      integer np
642      integer nl
643      double precision vcutsq
644      double precision c(np,3)
645      integer point(np)
646      integer list(nl)
647      double precision cl(nl,3)
648      double precision ct(np,3)
649      logical result
650c
651      integer i,j,k
652      integer nlist
653      integer eatm
654      double precision rij(3),rijsq,cc(1,3)
655
656      character*30 pname
657
658      pname = "smd_vlist_set"
659
660      result = .false.
661      nlist=0
662
663c      write(*,*) "vcutsq=",vcutsq
664      do i=1,np-1
665
666       point(i)=nlist+1
667
668       k = 0
669       do j=i+1,np
670
671         k = k +1
672         ct(k,1)=c(i,1)-c(j,1)
673         ct(k,2)=c(i,2)-c(j,2)
674         ct(k,3)=c(i,3)-c(j,3)
675
676       end do
677
678
679       call smd_lat_rebox(np,ct)
680
681
682       k = 0
683       do j=i+1,np
684
685        k = k + 1
686
687         rij(1)=ct(k,1)
688         rij(2)=ct(k,2)
689         rij(3)=ct(k,3)
690
691         rijsq=rij(1)*rij(1)+rij(2)*rij(2)+rij(3)*rij(3)
692
693
694         if(rijsq.lt.vcutsq)then
695
696          nlist=nlist+1
697
698          if(nlist.gt.nl)then
699           result = .false.
700           goto 200
701          endif
702
703          list(nlist)=j
704          cl(nlist,1)=rij(1)
705          cl(nlist,2)=rij(2)
706          cl(nlist,3)=rij(3)
707
708         endif
709
710       enddo
711
712      enddo
713
714      point(np)=nlist+1
715
716      nl = nlist
717
718      result = .true.
719
720200   continue
721      return
722
723      END
724
725      SUBROUTINE smd_vlist_test(lupdate)
726      implicit none
727#include "smd_system.fh"
728#include "mafdecls.fh"
729#include "errquit.fh"
730c
731      logical lupdate
732c
733      integer na
734      integer i_v,i_vd
735      double precision t
736      character*72 sp_vel
737      character*72 sp_vlist
738      logical result
739      character*72 tag
740      character*30 pname
741      double precision vcut,rcut
742c
743      pname = "smd_vlist_test"
744c
745c     get velocity aray
746c     -----------------
747      tag = "velocity"
748      call smd_system_get_component(sp_vel,tag,result)
749      if(.not.result)
750     >  call errquit(
751     >       pname//'no component '//tag,0,0)
752
753      tag = "vel"
754      call smd_data_get_index(sp_vel,tag,i_v,result)
755      if(.not. result)
756     >  call errquit(
757     >       pname//'error getting index for'//tag,0, 0)
758      call  smd_data_get_size(sp_vel,tag,na,result)
759      if(.not. result)
760     >  call errquit(
761     >       pname//'error getting index for'//tag,0, 0)
762      na = na/3
763c
764c     get time step
765c     ------------
766      if(.not.smd_system_tstep(t))
767     >  call errquit(
768     >       pname//'no time step ',0,0)
769
770c
771c     get verlet displacement array
772c     -----------------------------
773      tag = "verlet_list"
774      call smd_system_get_component(sp_vlist,tag,result)
775      if(.not.result)
776     >  call errquit(
777     >       pname//'no component '//tag,0,0)
778
779      tag = "vlist:displacement"
780      call smd_data_get_index(sp_vlist,tag,i_vd,result)
781      if(.not. result)
782     >  call errquit(
783     >       pname//'error getting index for'//tag,0, 0)
784
785       call smd_cutoff_get_rcut(rcut)
786       call smd_cutoff_get_rcut_verlet(vcut)
787
788       call smd_vlist_test0(na,t,rcut,vcut,
789     >     dbl_mb(i_vd),dbl_mb(i_v),lupdate)
790
791      return
792
793      END
794
795      SUBROUTINE smd_vlist_test0(natms,tstep,rcut,vcut,ivv,vvv,lupdate)
796
797      implicit none
798c
799      integer natms
800      double precision rcut,vcut
801      double precision tstep
802      logical lupdate
803      double precision ivv(natms,3)
804      double precision vvv(natms,3)
805c
806      integer i,exceed
807
808      double precision  tstepsq
809      double precision  dispmax,dispxsq,dispysq,dispzsq,disprsq
810
811      logical lnew
812
813      data lnew/.true./
814
815      save lnew
816
817      tstepsq=tstep**2
818
819      if(lnew)then
820
821       lupdate=.true.
822       lnew=.false.
823
824       do i=1,natms
825
826        ivv(i,1)=0.0
827        ivv(i,2)=0.0
828        ivv(i,3)=0.0
829
830       enddo
831
832      else
833
834       lupdate=.false.
835
836       dispmax=((vcut-rcut)/2.0)**2
837
838       do i=1,natms
839
840        ivv(i,1)=ivv(i,1)+vvv(i,1)
841        ivv(i,2)=ivv(i,2)+vvv(i,2)
842        ivv(i,3)=ivv(i,3)+vvv(i,3)
843
844       enddo
845
846       exceed=0
847
848       do i=1,natms
849
850        dispxsq=ivv(i,1)**2
851        dispysq=ivv(i,2)**2
852        dispzsq=ivv(i,3)**2
853        disprsq=tstepsq*(dispxsq+dispysq+dispzsq)
854        if(disprsq.gt.dispmax) then
855          exceed=exceed+1
856          write(*,*) "verlet update disp",disprsq,dispmax
857        end if
858        if(exceed.ge.2)lupdate=.true.
859
860       enddo
861
862       if(lupdate)then
863
864        do i=1,natms
865
866         ivv(i,1)=0
867         ivv(i,2)=0
868         ivv(i,3)=0
869
870        enddo
871
872       endif
873
874      endif
875
876      return
877
878      END
879c $Id$
880