PHOS/reconstruction/clusterlib.f

   1         SUBROUTINE CLUST(iwl,gamthr)
   2 C
   3 C      9.12.92  Kolosov V.N.
   4 C        GLASS CLUSTERS SEARCH
   5 C
   6         common /error/ier
   7
   8         include 'comgl.for'
   9         include 'comalwl.for'
  10         include 'comwlgen.for'
  11
  12         integer mark(nglal),icl(ialw)
  13
  14
  15         miwl=1
  16         mawl=nwall
  17 c       write (*,*) ' nwall ',nwall
  18         if(iwl.gt.0.and.iwl.le.nwall) then
  19                 miwl=iwl
  20                 mawl=iwl
  21         endif
  22
  23         do iw=miwl,mawl
  24         ip=ipwgl(iw)
  25         ng=nglw(iw)
  26 c       call wlf(wl,ng,mgl(ip+1),egl(ip+1))
  27         ncls(iw)=0
  28         ipwlcl(iw)=ntclst
  29
  30         nc=ntclst
  31         do k=ip+1,ip+ng
  32         m=mgl(k)
  33         e=egl(k)
  34 c       write (*,*) ' m,e ',m,e
  35         if(icl(m).eq.0.and.e.gt.gamthr) then
  36 c New cluster
  37                 if(nc.eq.nmxcl) then
  38                 write (*,*) ' max cluster err !!! ',nmxcl
  39                 ier=36
  40                 return
  41                 endif
  42                 nc=nc+1
  43 co      write (*,*) ' new cluster ',nc,ipcl
  44                 icl(m)=nc
  45                 ipclst(nc)=ipcl
  46                 ncl=1
  47                 ncld=1
  48                 DO WHILE (NCLD.gt.0)
  49                 NCLD=0
  50                 DO 2 K1=ip+1,ip+NG
  51 co      write (*,*) ' neib? #,m,icl',k1-ip,mgl(k1),icl(mgl(k1))
  52                 IF(MARK(K1).NE.0) GOTO 2
  53                 M1=MGl(K1)
  54                 IF(ICL(M1).NE.NC) GOTO 2
  55 co      write (*,*) ' search ns,wl',ns(m1),wl(m1)
  56                         DO L=1,NS(M1)
  57                         M1S=KS(L,M1)
  58 co      write (*,*) ' sosed ',l,m1s,wl(m1s)
  59                         IF(WL(M1S).NE.0..AND.ICL(M1S).EQ.0) THEN
  60                                 ICL(M1S)=NC
  61                                 NCLD=NCLD+1
  62                         ENDIF
  63                         ENDDO
  64                 MARK(K1)=NC
  65 c Save here
  66                 if(ipcl.eq.nglal) then
  67                 write (*,*) ' over cluster buff !!! ',ipcl
  68                 ier=37
  69                 return
  70                 endif
  71                 eclst(ipcl+1)=egl(k1)
  72                 mclst(ipcl+1)=mgl(k1)
  73                 ipcl=ipcl+1
  74 2               CONTINUE
  75                 ncl=ncl+ncld
  76                 ENDDO
  77                 lgclst(nc)=ncl
  78 C   END  NEW  CLUSTER
  79         endif
  80         enddo
  81         ncls(iw)=nc-ntclst
  82         ntclst=ntclst+nc
  83 c  zero marks
  84         DO K=ip+1,ip+NG
  85         m=mgl(k)
  86         MARK(K)=0
  87         icl(m)=0
  88         ENDDO
  89 c New wall
  90         enddo
  91
  92         return
  93
  94         entry zclust
  95
  96         ipcl=0
  97         ntclst=0
  98         do i=1,nalw
  99         ncls(i)=0
 100         enddo
 101         do i=1,ialw
 102         icl(i)=0
 103         enddo
 104         do i=1,nglal
 105         mark(i)=0
 106         enddo
 107
 108         return
 109         end
 110
 111         SUBROUTINE cCLUST(iwl)
 112         common /error/ier
 113
 114         include 'comgl.for'
 115         include 'comalwl.for'
 116         include 'comwlgen.for'
 117
 118         real xt(3),xmi(3),xma(3)
 119
 120         miwl=1
 121         mawl=nwall
 122         if(iwl.gt.0.and.iwl.le.nwall) then
 123                 miwl=iwl
 124                 mawl=iwl
 125         endif
 126         do iw=miwl,mawl
 127 c new wall
 128         ipwcc=ipwlcl(iw)
 129         etwc=0.
 130         nck(iw)=ncls(iw)
 131
 132         do i=1,ncls(iw)
 133 c new cluster
 134         lc=lgclst(i+ipwcc)
 135         ipgcc=ipclst(i+ipwcc)
 136         et=0.
 137         emax=0.
 138         call vzero(xt,3)
 139         call ucopy(x(1,mclst(ipgcc+1)),xmi,3)
 140         call ucopy(x(1,mclst(ipgcc+1)),xma,3)
 141                 do ni=1,lc
 142 c new glass in cluster
 143                 nic=ipgcc+ni
 144                 mgc=mclst(nic)
 145                 egc=eclst(nic)
 146                 if(emax.lt.egc) then
 147                         emax=egc
 148                         mmax=mgc
 149                 endif
 150                 if(xmi(1).gt.x(1,mgc)) xmi(1)=x(1,mgc)
 151                 if(xmi(2).gt.x(2,mgc)) xmi(2)=x(2,mgc)
 152                 if(xmi(3).gt.x(3,mgc)) xmi(3)=x(3,mgc)
 153                 if(xma(1).lt.x(1,mgc)) xma(1)=x(1,mgc)
 154                 if(xma(2).lt.x(2,mgc)) xma(2)=x(2,mgc)
 155                 if(xma(3).lt.x(3,mgc)) xma(3)=x(3,mgc)
 156
 157                 xt(1)=xt(1)+x(1,mgc)*egc
 158                 xt(2)=xt(2)+x(2,mgc)*egc
 159                 xt(3)=xt(3)+x(3,mgc)*egc
 160                 et=et+egc
 161 c end glass in cluster
 162                 enddo
 163         if(et.eq.0.) then
 164                 write (*,*) ' error!! energy cluster=0 '
 165                 write (*,*) ' wall,nc,lc,ipgcc ',iw,i,lc,ipgcc
 166                 ier=33
 167                 return
 168         endif
 169         xt(1)=xt(1)/et
 170         xt(2)=xt(2)/et
 171         xt(3)=xt(3)/et
 172
 173         eck(i,iw)=et
 174         emxck(i,iw)=emax
 175         call ucopy(xt,xavck(1,i,iw),3)
 176         call ucopy(xmi,xmick(1,i,iw),3)
 177         call ucopy(xma,xmack(1,i,iw),3)
 178
 179         etwc=etwc+et
 180 c end cluster
 181         enddo
 182         etck(iw)=etwc
 183 c end wall
 184         enddo
 185
 186         return
 187         end