text:
*dk,refinefc
      subroutine refine_face()
C
C     CHANGE HISTORY -
C
C        $Log: refine_face.f,v $
C        Revision 2.00  2007/11/09 20:04:01  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.7   05 May 2000 15:33:58   dcg
CPVCS    refresh mbndry value after cmo_newlen
CPVCS
CPVCS       Rev 1.6   Fri Jan 22 15:46:46 1999   dcg
CPVCS    Define the tolerance xst2 to be 1.d-9
CPVCS
CPVCS       Rev 1.5   Mon Apr 14 16:59:24 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.4   11/07/95 17:24:40   dcg
CPVCS    change flag to 2 in mmgetblk calls
CPVCS
CPVCS       Rev 1.3   09/29/95 09:14:10   het
CPVCS    Put in added attributes inheritance
CPVCS
CPVCS       Rev 1.2   06/05/95 10:36:54   het
CPVCS    Make changes for hybrid_grids
CPVCS
CPVCS       Rev 1.1   05/26/95 13:17:52   het
CPVCS    Replace subroutine parameter list with subroutine calles.
CPVCS
CPVCS       Rev 1.0   11/10/94 12:18:00   pvcs
CPVCS    Original version.
C
      implicit none
C
      integer nplen,nvalues
      parameter (nplen=10000000)
      parameter (nvalues=3)
C
      pointer (ipitet, itet)
      pointer (ipitet, itet1)
      pointer (ipjtet, jtet)
      pointer (ipjtet, jtet1)
      integer itet(4,*), jtet(4,*)
      integer itet1(*), jtet1(*)

      pointer (ipitetclr, itetclr)
      pointer (ipitettyp, itettyp)
      pointer (ipitetoff, itetoff)
      pointer (ipjtetoff, jtetoff)
      integer itetclr(*), itettyp(*),
     *        itetoff(*), jtetoff(*)
C
      pointer (ipitflag, itflag)
      integer itflag(*)

      pointer (ipitetnn, itetnn)
      pointer (ipitetnn1, itetnn1)
      pointer (ipitetnn2, itetnn2)
      integer itetnn(4,*), itetnn1(4,*), itetnn2(4,*)

      pointer (ipkfix, kfix)
      pointer (ipkfix, kfix1)
      integer kfix(4,*), kfix1(*)
C
      pointer (iplist_sink, list_sink)
      pointer (iplist_source, list_source)
      integer list_sink(*), list_source(nvalues,*)
C
      integer ierror,ier,ierrw,ics,icscode
      integer npoints,length,icmotype,ntets,mbndry,nen,nef,
     *        nface1,icount,ibound,i,i1,i2,i3,i4,i5,
     *        it,nedge,j,j1,j2,j3,j4,l,l1,l2,l3,
     *        ifc,jt,jf,if,iclrt1,iclrt2,
     *        ifaceiter,nface1_save
      integer lenitetclr,lenitettyp,lenitetoff,lenjtetoff,
     *        lenitet,lenjtet,lenxic,lenyic,lenzic
      integer nadd1,iedgeiter,npointsnew,ntetsnew,nedge_save
      integer k,kf,ke,kt,irefine,iface,nelementsmm,
     *        ntetsinc,nnodesmm,npointsinc,inc,inc1,inc2,
     *        kflast,ktlast,kelast,jcount,idnum,idum,
     *        itstart,itlast,ifstart,ielast,iestart,iflast

      real*8 xa,ya,za,xb,yb,zb,xc,yc,zc,xn,yn,zn,
     *       xd,yd,zd,xl,yl,zl,xv,yv,zv
      real*8 xdotmin,xdot,xdotl,ds23,ds2i,ds3i
      real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,xint,yint,zint
      real*8 xxsmall,q,vol,dist
C
      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      real*8 xic(*), yic(*), zic(*)

      pointer (ipxfix, xfix)
      real*8 xfix(4,*)
      pointer (ipxweight_source, xweight_source)
      real*8 xweight_source(nvalues,*)

c  a tolerance
      real*8 xst2
      data xst2/1.0d-9/

      integer itetface0(4), itetface1(4,4)
      data itetface0 / 3, 3, 3, 3 /
      data itetface1 / 2, 3, 4, 1,
     *                 1, 4, 3, 2,
     *                 1, 2, 4, 3,
     *                 1, 3, 2, 4 /

      integer itetface2(3,3,4)
      data itetface2 / 3, 4, 2,
     *                         4, 2, 3,
     *                         2, 3, 4,
     *                         4, 3, 1,
     *                         3, 1, 4,
     *                         1, 4, 3,
     *                         2, 4, 1,
     *                         4, 1, 2,
     *                         1, 2, 4,
     *                         3, 2, 1,
     *                         2, 1, 3,
     *                         1, 3, 2 /

      integer itetface3(2,3,4)
      data itetface3 / 2, 1,
     *                 3, 1,
     *                 4, 1,
     *                 1, 1,
     *                 4, 3,
     *                 3, 2,
     *                 1, 2,
     *                 2, 3,
     *                 4, 2,
     *                 1, 3,
     *                 3, 3,
     *                 2, 2 /

      real*8 a,b,c,d,e,f,crosx,crosy,crosz
      crosx(a,b,c,d,e,f)=b*f-c*e
      crosy(a,b,c,d,e,f)=c*d-a*f
      crosz(a,b,c,d,e,f)=a*e-b*d

      real*8 crosx1,crosy1,crosz1,volume
      crosx1(i,j,k)=(yic(j)-yic(i))*(zic(k)-zic(i))-
     *              (yic(k)-yic(i))*(zic(j)-zic(i))
      crosy1(i,j,k)=(xic(k)-xic(i))*(zic(j)-zic(i))-
     *              (xic(j)-xic(i))*(zic(k)-zic(i))
      crosz1(i,j,k)=(xic(j)-xic(i))*(yic(k)-yic(i))-
     *              (xic(k)-xic(i))*(yic(j)-yic(i))
      volume(i1,i2,i3,i4)=(xic(i4)-xic(i1))*crosx1(i1,i2,i3)+
     *                    (yic(i4)-yic(i1))*crosy1(i1,i2,i3)+
     *                    (zic(i4)-zic(i1))*crosz1(i1,i2,i3)
C
      character*132 logmess
      character*32 cmo, cmolength
      character*32 isubname, iblknam, iprtnam

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C BEGIN begin
C
      isubname='refine_face'
C
      call cmo_get_name(cmo,ierror)
      call cmo_get_info('nnodes',cmo,npoints,length,icmotype,ierror)
      call cmo_get_info('nelements',cmo,ntets,length,icmotype,ierror)
      call cmo_get_info('mbndry',cmo,mbndry,length,icmotype,ierror)
      call cmo_get_info('nodes_per_element',cmo,
     *                  nen,length,icmotype,ierror)
      call cmo_get_info('faces_per_element',cmo,
     *                  nef,length,icmotype,ierror)
      call cmo_get_info('xic',cmo,ipxic,lenxic,icmotype,ierror)
      call cmo_get_info('yic',cmo,ipyic,lenyic,icmotype,ierror)
      call cmo_get_info('zic',cmo,ipzic,lenzic,icmotype,ierror)
      call cmo_get_info('itetclr',cmo,ipitetclr,lenitetclr,icmotype,ier)
      call cmo_get_info('itettyp',cmo,ipitettyp,lenitettyp,icmotype,ier)
      call cmo_get_info('itetoff',cmo,ipitetoff,lenitetoff,icmotype,ier)
      call cmo_get_info('jtetoff',cmo,ipjtetoff,lenjtetoff,icmotype,ier)
      call cmo_get_info('itet',cmo,ipitet,lenitet,icmotype,ierror)
      call cmo_get_info('jtet',cmo,ipjtet,lenjtet,icmotype,ierror)
C
      xxsmall=1.0d-30
      length=ntets
      call mmgetblk("itflag",isubname,ipitflag,length,1,icscode)
      length=4*ntets
      call mmgetblk("itetnn" ,isubname,ipitetnn ,length,1,icscode)
      call mmgetblk("itetnn1",isubname,ipitetnn1,length,1,icscode)
      call mmgetblk("itetnn2",isubname,ipitetnn2,length,1,icscode)
      call mmgetblk("kfix",isubname,ipkfix,length,2,icscode)
      call mmgetblk("xfix",isubname,ipxfix,length,2,icscode)
C
      do it=1,ntets
         do i=1,4
            itetnn(i,it)=itet(i,it)
            if(jtet(i,it).ge.mbndry) then
               itetnn1(i,it)=0
               itetnn2(i,it)=0
            else
               itetnn1(i,it)=1+(jtet(i,it)-1)/4
               itetnn2(i,it)=jtet(i,it)-4*(itetnn1(i,it)-1)
            endif
         enddo
      enddo
      do it=1,ntets
         itflag(it)=0
      enddo
      irefine=0
      npointsnew=npoints
      ntetsnew=ntets
      do it=1,ntets
         do j=1,4
            kfix(j,it)=0
         enddo
      enddo
      do 10 i=1,ntets
         xa=xic(itet(1,i))
         ya=yic(itet(1,i))
         za=zic(itet(1,i))
         xb=xic(itet(2,i))-xa
         yb=yic(itet(2,i))-ya
         zb=zic(itet(2,i))-za
         xc=xic(itet(3,i))-xa
         yc=yic(itet(3,i))-ya
         zc=zic(itet(3,i))-za
         xd=xic(itet(4,i))-xa
         yd=yic(itet(4,i))-ya
         zd=zic(itet(4,i))-za
         xn=crosx(xb,yb,zb,xc,yc,zc)
         yn=crosy(xb,yb,zb,xc,yc,zc)
         zn=crosz(xb,yb,zb,xc,yc,zc)
         x2=crosx(xn,yn,zn,xb,yb,zb)
         y2=crosy(xn,yn,zn,xb,yb,zb)
         z2=crosz(xn,yn,zn,xb,yb,zb)
         q=-0.5*(xc*xb+yc*yb+zc*zb-xc*xc-yc*yc-zc*zc)/
     *          (x2*xc+y2*yc+z2*zc+xxsmall)
         xl=q*x2+0.5*xb
         yl=q*y2+0.5*yb
         zl=q*z2+0.5*zb
         d=-0.5*(xd*xd+yd*yd+zd*zd)
         q=-(xd*xl+yd*yl+zd*zl+d)/(xd*xn+yd*yn+zd*zn+xxsmall)
         vol=xn*xd+yn*yd+zn*zd
         xv=q*xn+xl+xa
         yv=q*yn+yl+ya
         zv=q*zn+zl+za
         dist=(xv-xa)**2+(yv-ya)**2+(zv-za)**2
C********j1=cvmgn(jtet(1,i),1,jtet(1,i)-mbndry)
C********j1=cvmgm(j1,j1-mbndry,j1-mbndry)
         if((jtet(1,i)-mbndry).ne.0) then
            j1=jtet(1,i)
         else
            j1=1
         endif
         if((j1-mbndry).lt.0) then
            j1=j1
         else
            j1=j1-mbndry
         endif
         xfix(1,i)=(xic(itet1(j1))-xv)**2+
     *             (yic(itet1(j1))-yv)**2+
     *             (zic(itet1(j1))-zv)**2
C********j1=cvmgn(jtet(2,i),1,jtet(2,i)-mbndry)
C********j1=cvmgm(j1,j1-mbndry,j1-mbndry)
         if((jtet(2,i)-mbndry).ne.0) then
            j1=jtet(2,i)
         else
            j1=1
         endif
         if((j1-mbndry).lt.0) then
            j1=j1
         else
            j1=j1-mbndry
         endif
         xfix(2,i)=(xic(itet1(j1))-xv)**2+
     *             (yic(itet1(j1))-yv)**2+
     *             (zic(itet1(j1))-zv)**2
C********j1=cvmgn(jtet(3,i),1,jtet(3,i)-mbndry)
C********j1=cvmgm(j1,j1-mbndry,j1-mbndry)
         if((jtet(3,i)-mbndry).ne.0) then
            j1=jtet(3,i)
         else
            j1=1
         endif
         if((j1-mbndry).lt.0) then
            j1=j1
         else
            j1=j1-mbndry
         endif
         xfix(3,i)=(xic(itet1(j1))-xv)**2+
     *             (yic(itet1(j1))-yv)**2+
     *             (zic(itet1(j1))-zv)**2
C********j1=cvmgn(jtet(4,i),1,jtet(4,i)-mbndry)
C********j1=cvmgm(j1,j1-mbndry,j1-mbndry)
         if((jtet(4,i)-mbndry).ne.0) then
            j1=jtet(4,i)
         else
            j1=1
         endif
         if((j1-mbndry).lt.0) then
            j1=j1
         else
            j1=j1-mbndry
         endif
         xfix(4,i)=(xic(itet1(j1))-xv)**2+
     *             (yic(itet1(j1))-yv)**2+
     *             (zic(itet1(j1))-zv)**2
         ifc=4*(i-1)
C********kfix(1,i)=cvmgp(ifc+1,0,dist-xfix(1,i)-xst2*dist)
C********kfix(2,i)=cvmgp(ifc+2,0,dist-xfix(2,i)-xst2*dist)
C********kfix(3,i)=cvmgp(ifc+3,0,dist-xfix(3,i)-xst2*dist)
C********kfix(4,i)=cvmgp(ifc+4,0,dist-xfix(4,i)-xst2*dist)
C********kfix(1,i)=cvmgm(kfix(1,i),-i,-vol)
         if(dist-xfix(1,i)-xst2*dist .ge. 0.0) then
            kfix(1,i)=ifc+1
         else
            kfix(1,i)=0
         endif
         if(dist-xfix(2,i)-xst2*dist .ge. 0.0) then
            kfix(2,i)=ifc+2
         else
            kfix(2,i)=0
         endif
         if(dist-xfix(3,i)-xst2*dist .ge. 0.0) then
            kfix(3,i)=ifc+3
         else
            kfix(3,i)=0
         endif
         if(dist-xfix(4,i)-xst2*dist .ge. 0.0) then
            kfix(4,i)=ifc+4
         else
            kfix(4,i)=0
         endif
         if(-vol.lt.0.0) then
            kfix(1,i)=kfix(1,i)
         else
            kfix(1,i)=-i
         endif
         do k=1,4
            if(jtet(k,i).gt.0.and.jtet(k,i).lt.mbndry) then
               jt=1+(jtet(k,i)-1)/4
               jf=jtet(k,i)-4*(jt-1)
               iclrt1=itetclr(i)
               iclrt2=itetclr(jt)
               if(iclrt1.ne.iclrt2) then
                  if(iclrt1.ne.0.and.iclrt2.ne.0) then
C*****               kfix(k,i)=0
                  endif
               endif
             endif
          enddo
  10  continue

C compiler: WARNING This statement can not be reached.
C     nface1=0
C     ^

      nface1=0
      do it=1,ntets
         do j=1,4
            jt=itetnn1(j,it)
            jf=itetnn2(j,it)
            i1=itet(jf,jt)
            i2=itet(itetface1(1,jf),jt)
            i3=itet(itetface1(2,jf),jt)
            i4=itet(itetface1(3,jf),jt)
            if(jt.gt.0.and.jt.le.ntets) then
               if(itetclr(it).ne.itetclr(jt)) then
                  if(kfix(j,it).gt.0) then
                     if(nface1.gt.0) then
                        do l=1,nface1
                           kt=1+(kfix1(l)-1)/4
                           kf=kfix1(l)-4*(kt-1)
                           j1=itet(kf,kt)
                           j2=itet(itetface1(1,kf),kt)
                           j3=itet(itetface1(2,kf),kt)
                           j4=itet(itetface1(3,kf),kt)
                           if((i2.eq.j2.and.i3.eq.j3.and.i4.eq.j4).or.
     *                        (i2.eq.j3.and.i3.eq.j4.and.i4.eq.j2).or.
     *                        (i2.eq.j4.and.i3.eq.j2.and.i4.eq.j3)) then
                              goto 5
                           endif
                        enddo
                     endif
                     nface1=nface1+1
                     kfix1(nface1)=kfix(j,it)
   5                 continue
                  endif
               endif
            endif
         enddo
      enddo
      goto 7

C  From compiler
C  WARNING: This statement can not be reached.
C      nface1=0
C      ^

      nface1=0
      do it=1,ntets
         do j=1,4
            kfix(j,it)=4*(it-1)+j
            jt=itetnn1(j,it)
            jf=itetnn2(j,it)
            i1=itet(jf,jt)
            i2=itet(itetface1(1,jf),jt)
            i3=itet(itetface1(2,jf),jt)
            i4=itet(itetface1(3,jf),jt)
            if(jt.gt.0.and.jt.le.ntets) then
               if(nface1.gt.0) then
                  do l=1,nface1
                     kt=1+(kfix1(l)-1)/4
                     kf=kfix1(l)-4*(kt-1)
                     j1=itet(kf,kt)
                     j2=itet(itetface1(1,kf),kt)
                     j3=itet(itetface1(2,kf),kt)
                     j4=itet(itetface1(3,kf),kt)
                     if((i2.eq.j2.and.i3.eq.j3.and.i4.eq.j4).or.
     *                  (i2.eq.j3.and.i3.eq.j4.and.i4.eq.j2).or.
     *                  (i2.eq.j4.and.i3.eq.j2.and.i4.eq.j3)) then
                        goto 6
                     endif
                  enddo
               endif
               nface1=nface1+1
               kfix1(nface1)=kfix(j,it)
   6           continue
            endif
         enddo
      enddo
   7  continue
      ifaceiter=0
C
      length=nface1
      call mmgetblk('list_sink',isubname,iplist_sink,length,1,icscode)
      length=nvalues*nface1
      call mmgetblk('list_source',isubname,iplist_source,length,1,
     *              icscode)
      call mmgetblk('xweight_source',isubname,ipxweight_source,length,2,
     *              icscode)
C
      nadd1=0
  11  continue
      do it=1,ntets
         do i=1,4
            itetnn(i,it)=itet(i,it)
            if(jtet(i,it).ge.mbndry) then
               itetnn1(i,it)=0
               itetnn2(i,it)=0
            else
               itetnn1(i,it)=1+(jtet(i,it)-1)/4
               itetnn2(i,it)=jtet(i,it)-4*(itetnn1(i,it)-1)
            endif
         enddo
      enddo
      do it=1,ntets
         itflag(it)=0
      enddo
      npointsnew=npoints
      ntetsnew=ntets
        ifaceiter=ifaceiter+1
        write(logmess,'(a,2i10)') "Face iteration: ",ifaceiter,nface1
        call writloga('default',0,logmess,0,ierrw)
      do it=1,ntets
         itflag(it)=0
      enddo
      write(logmess,'(a,i10)')
     *   "Interface face that need reconnecting: ",nface1
      call writloga('default',0,logmess,0,ierrw)
      if(nface1.ne.0) then
         nface1_save=0
         do iface=1,nface1
            it=1+(kfix1(iface)-1)/4
            if=kfix1(iface)-4*(it-1)
            i1=itet(if,it)
            i2=itet(itetface1(1,if),it)
            i3=itet(itetface1(2,if),it)
            i4=itet(itetface1(3,if),it)
            jt=1+(jtet(if,it)-1)/4
            jf=jtet(if,it)-4*(jt-1)
            i5=itet(jf,jt)
            if(it.lt.0.or.jt.lt.0) then
               write(logmess,'(a,2i10)') "Error1: ",it,jt
               call writloga('default',0,logmess,0,ierrw)
               stop
            endif
            if(it.gt.ntets.or.jt.gt.ntets) then
               goto 20
            endif
            if(itflag(it).eq.0.and.itflag(jt).eq.0) then
               irefine=irefine+1
               call mmfindbk('xic',cmo,ipxic,length,icscode)
               if((npointsnew+1).gt.length) then
                  npointsinc=npointsnew+1000
                  call cmo_set_info('nnodes',cmo,npointsinc,1,1,ierror)
                  call mmgetlen(ipitetclr,nelementsmm,icscode)
                  call cmo_set_info('nelements',cmo,nelementsmm,1,1,ier)
                  call cmo_newlen(cmo,ierror)
                  call cmo_get_intinfo('mbndry',cmo,mbndry,length,
     *                    icmotype,ierror)
                  call cmo_get_info('xic',cmo,ipxic,lenxic,icmotype,ier)
                  call cmo_get_info('yic',cmo,ipyic,lenyic,icmotype,ier)
                  call cmo_get_info('zic',cmo,ipzic,lenzic,icmotype,ier)
               endif
               npointsnew=npointsnew+1
C
               nadd1=nadd1+1
               list_sink(nadd1)=npointsnew
               list_source(1,nadd1)=i2
               list_source(2,nadd1)=i3
               list_source(3,nadd1)=i4
               xweight_source(1,nadd1)=1.0
               xweight_source(2,nadd1)=1.0
               xweight_source(3,nadd1)=1.0
C
               call mmgetlen(ipitflag,length,icscode)
               if((ntetsnew+4).gt.length) then
                  inc=1000
                  ntetsinc=ntetsnew+inc
                  call cmo_set_info('nelements',cmo,ntetsinc,1,1,ierror)
                  call mmfindbk('xic',cmo,ipxic,nnodesmm,icscode)
                  call cmo_set_info('nnodes',cmo,nnodesmm,1,1,ierror)
                  call cmo_newlen(cmo,ierror)
                  call cmo_get_intinfo('mbndry',cmo,mbndry,length,
     *                    icmotype,ierror)
                  call cmo_get_info('itetclr',cmo,
     *                              ipitetclr,lenitetclr,icmotype,ier)
                  call cmo_get_info('itettyp',cmo,
     *                              ipitettyp,lenitettyp,icmotype,ier)
                  call cmo_get_info('itetoff',cmo,
     *                              ipitetoff,lenitetoff,icmotype,ier)
                  call cmo_get_info('jtetoff',cmo,
     *                              ipjtetoff,lenjtetoff,icmotype,ier)
                  call cmo_get_info('itet',cmo,
     *                              ipitet,lenitet,icmotype,ierror)
                  call cmo_get_info('jtet',cmo,
     *                              ipjtet,lenjtet,icmotype,ierror)
                  call mmgetnam(ipitflag,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipitflag,inc,
     *                          ics)
                  do idum=ntetsnew+1,ntetsnew+inc
                     itflag(idum)=0
                  enddo
                  inc1=nen*inc
                  call mmgetnam(ipitetnn,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipitetnn,inc1,
     *                          ics)
                  inc2=nef*inc
                  call mmgetnam(ipitetnn1,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipitetnn1,inc2,
     *                          ics)
                  call mmgetnam(ipitetnn2,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipitetnn2,inc2,
     *                          ics)
                  call mmgetnam(ipkfix,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipkfix,inc2,
     *                          ics)
                  call mmgetnam(ipxfix,iblknam,iprtnam,ics)
                  call mmincblk(iblknam,iprtnam,ipxfix,inc2,
     *                          ics)
               endif
               itflag(it)=1
               itetnn(1,it)=i1
               itetnn(2,it)=i2
               itetnn(3,it)=i3
               itetnn(4,it)=npointsnew
               itetnn1(1,it)=-1
               itetnn1(2,it)=-1
               itetnn1(3,it)=-1
               itetnn1(4,it)=-1
               itetnn2(1,it)=-1
               itetnn2(2,it)=-1
               itetnn2(3,it)=-1
               itetnn2(4,it)=-1
               ntetsnew=ntetsnew+1
                  itetclr(ntetsnew)=itetclr(it)
                  itettyp(ntetsnew)=itettyp(it)
                  itetoff(ntetsnew)=nen*(ntetsnew-1)
                  jtetoff(ntetsnew)=nef*(ntetsnew-1)
                  itetnn(1,ntetsnew)=i1
                  itetnn(2,ntetsnew)=i4
                  itetnn(3,ntetsnew)=i2
                  itetnn(4,ntetsnew)=npointsnew
                  itetnn1(1,ntetsnew)=-1
                  itetnn1(2,ntetsnew)=-1
                  itetnn1(3,ntetsnew)=-1
                  itetnn1(4,ntetsnew)=-1
                  itetnn2(1,ntetsnew)=-1
                  itetnn2(2,ntetsnew)=-1
                  itetnn2(3,ntetsnew)=-1
                  itetnn2(4,ntetsnew)=-1
               ntetsnew=ntetsnew+1
                  itetclr(ntetsnew)=itetclr(jt)
                  itettyp(ntetsnew)=itettyp(jt)
                  itetoff(ntetsnew)=nen*(ntetsnew-1)
                  jtetoff(ntetsnew)=nef*(ntetsnew-1)
                  itetnn(1,ntetsnew)=i5
                  itetnn(2,ntetsnew)=i2
                  itetnn(3,ntetsnew)=i4
                  itetnn(4,ntetsnew)=npointsnew
                  itetnn1(1,ntetsnew)=-1
                  itetnn1(2,ntetsnew)=-1
                  itetnn1(3,ntetsnew)=-1
                  itetnn1(4,ntetsnew)=-1
                  itetnn2(1,ntetsnew)=-1
                  itetnn2(2,ntetsnew)=-1
                  itetnn2(3,ntetsnew)=-1
                  itetnn2(4,ntetsnew)=-1
               ntetsnew=ntetsnew+1
                  itetclr(ntetsnew)=itetclr(jt)
                  itettyp(ntetsnew)=itettyp(jt)
                  itetoff(ntetsnew)=nen*(ntetsnew-1)
                  jtetoff(ntetsnew)=nef*(ntetsnew-1)
                  itetnn(1,ntetsnew)=i5
                  itetnn(2,ntetsnew)=i4
                  itetnn(3,ntetsnew)=i3
                  itetnn(4,ntetsnew)=npointsnew
                  itetnn1(1,ntetsnew)=-1
                  itetnn1(2,ntetsnew)=-1
                  itetnn1(3,ntetsnew)=-1
                  itetnn1(4,ntetsnew)=-1
                  itetnn2(1,ntetsnew)=-1
                  itetnn2(2,ntetsnew)=-1
                  itetnn2(3,ntetsnew)=-1
                  itetnn2(4,ntetsnew)=-1
               ntetsnew=ntetsnew+1
                  itetclr(ntetsnew)=itetclr(jt)
                  itettyp(ntetsnew)=itettyp(jt)
                  itetoff(ntetsnew)=nen*(ntetsnew-1)
                  jtetoff(ntetsnew)=nef*(ntetsnew-1)
                  itetnn(1,ntetsnew)=i5
                  itetnn(2,ntetsnew)=i3
                  itetnn(3,ntetsnew)=i2
                  itetnn(4,ntetsnew)=npointsnew
                  itetnn1(1,ntetsnew)=-1
                  itetnn1(2,ntetsnew)=-1
                  itetnn1(3,ntetsnew)=-1
                  itetnn1(4,ntetsnew)=-1
                  itetnn2(1,ntetsnew)=-1
                  itetnn2(2,ntetsnew)=-1
                  itetnn2(3,ntetsnew)=-1
                  itetnn2(4,ntetsnew)=-1
               itflag(jt)=1
               itetclr(jt)=itetclr(it)
               itettyp(jt)=itettyp(it)
               itetoff(jt)=nen*(jt-1)
               jtetoff(jt)=nef*(jt-1)
               itetnn(1,jt)=i1
               itetnn(2,jt)=i3
               itetnn(3,jt)=i4
               itetnn(4,jt)=npointsnew
               itetnn1(1,jt)=-1
               itetnn1(2,jt)=-1
               itetnn1(3,jt)=-1
               itetnn1(4,jt)=-1
               itetnn2(1,jt)=-1
               itetnn2(2,jt)=-1
               itetnn2(3,jt)=-1
               itetnn2(4,jt)=-1
            else
               nface1_save=nface1_save+1
               kfix1(nface1_save)=kfix1(iface)
            endif
   20       continue
         enddo
      endif
      write(logmess,'(a,i10,a,i10)')
     *   "Face-refined tets: old=",ntets," new=",ntetsnew
      call writloga('default',0,logmess,0,ierrw)
      if(ntetsnew.gt.ntets) then
         do it=1,ntets
            if(itflag(it).ne.0) then
               do i=1,4
                  kt=1+(jtet(i,it)-1)/4
                  kf=jtet(i,it)-4*(kt-1)
                  if(kt.le.ntets) then
                     itetnn1(kf,kt)=-1
                     itetnn2(kf,kt)=-1
                  endif
                  itetnn1(i,it)=-1
                  itetnn2(i,it)=-1
               enddo
            endif
         enddo
         npoints=npointsnew
         ntets=ntetsnew
         call cmo_set_info('nnodes',cmo,npoints,1,1,ierror)
         call cmo_set_info('nelements',cmo,ntets,1,1,ierror)
         call geniee(itetnn,itetnn1,itetnn2,4,4,ntets,npoints,
     *        3,npoints,ntets)
         do it=1,ntets
            do i=1,4
               itet(i,it)=itetnn(i,it)
               if(itetnn1(i,it).gt.0.and.itetnn1(i,it).le.ntets) then
                  jtet(i,it)=4*(itetnn1(i,it)-1)+itetnn2(i,it)
               else
                  jtet(i,it)=mbndry
               endif
            enddo
         enddo
         if(nface1_save.gt.0) then
            nface1=nface1_save
            goto 11
         endif
      endif
      goto 9999
 9999 continue
C
      call cmo_get_name(cmo,ierror)
C
      cmolength='nnodes'
      call cmo_interpolate(cmo,cmo,
     *                     cmolength,
     *                     nadd1,nvalues,
     *                     list_sink,list_source,xweight_source,
     *                     ierror)
      if(ierror.ne.0) call x3d_error(isubname,'cmo_interpolate',ier)
C
      call cmo_set_info('nnodes',cmo,npoints,1,1,ierror)
      call cmo_set_info('nelements',cmo,ntets,1,1,ierror)
C
      call mmrelprt(isubname,icscode)
C
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


! { dg-do compile }
!
! PR fortran/55763
!
! Contributed by Harald Anlauf
!

module gfcbug122
  implicit none
  type myobj
     class(*), allocatable :: x
   contains
     procedure :: print
  end type myobj
contains
  subroutine print(this)
    class(myobj) :: this
    select type (this)
    type is (integer) ! { dg-error "Unexpected intrinsic type 'INTEGER'" }
    type is (real) ! { dg-error "Unexpected intrinsic type 'REAL'" }
    type is (complex) ! { dg-error "Unexpected intrinsic type 'COMPLEX'" }
    type is (character(len=*)) ! { dg-error "Unexpected intrinsic type 'CHARACTER'" }
    end select
  end subroutine print
end module gfcbug122


text:
c     end ranfkbinary

c==============================================================================

      subroutine ranfmodmult( A, B, C )

c------------------------------------------------------------------------------
c
c     Ths routine calculates the product of the first two arguments.  All three
c     arguments are represented as arrays of 12-bit integers, each making up
c     the digits of one radix-4096 number.  The multiplication is done
c     piecemeal.
c
c     It returns the product in the third argument.
c
c------------------------------------------------------------------------------

      include 'pranf.inc'

      integer A, B, C, j1, j2, j3, j4, k1, k2, k3, k4
      integer a1,a2,a3,b1,b2,b3

      dimension A( IRandNumSize ),
     1     B( IRandNumSize ),
     2     C( IRandNumSize )

c------------------------------------------------------------------------------

c     write( *, 10 )
 10   format( ' ranfmodmult', /, '------------' )

c     j1 = A( 1 ) * B( 1 )
c     j2 = A( 1 ) * B( 2 ) + A( 2 ) * B( 1 )
c     j3 = A( 1 ) * B( 3 ) + A( 2 ) * B( 2 ) + A( 3 ) * B( 1 )
c     j4= A( 1 ) * B( 4 ) + A( 2 ) * B( 3 ) + A( 3 ) * B( 2 ) + A( 4 ) * B( 1 )
      a1 = A(1)
      a2 = A(2)
      a3 = A(3)
      b1 = B(1)
      b2 = B(2)
      b3 = B(3)
      j1 = a1 * b1
      j2 = a1 * b2 + a2 * b1
      j3 = a1 * b3 + a2 * b2 + a3 * b1
      j4 = a1 * B( 4 ) + a2 * b3 + a3 * b2 + A( 4 ) * b1

      k1 = j1
      k2 = j2 + k1 / 4096
      k3 = j3 + k2 / 4096
      k4 = j4 + k3 / 4096

      C( 1 ) = mod( k1, 4096 )
      C( 2 ) = mod( k2, 4096 )
      C( 3 ) = mod( k3, 4096 )
      C( 4 ) = mod( k4, 4096 )

c     write( *, 20 )
 20   format( '  leaving ranfmodmult' )

      return
      end

ext: f

max_stars_repo_name: gmlueck/llvm-test-suite

max_issues_repo_name: gmlueck/llvm-test-suite

max_forks_repo_name: gmlueck/llvm-test-suite


text:
SUBROUTINE rod_km(km,ea,length)
!
! This subroutine forms the stiffness matrix of a 1-d "rod" element.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::ea,length
 REAL(iwp),INTENT(OUT)::km(:,:)
 REAL(iwp)::one=1.0_iwp
 km(1,1)=one
 km(2,2)=one
 km(1,2)=-one
 km(2,1)=-one
 km=km*ea/length
RETURN
END SUBROUTINE rod_km

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
program test_tail
use iso_fortran_env,only:int8,int16,int32,int64,real32,real64,real128
use mod_testing,only:assert,initialize_tests,report_tests
use mod_functional

implicit none

logical,dimension(:),allocatable :: tests
logical :: test_failed
integer :: n,norder,ntests
integer,parameter :: stdout = 6

complex(kind=real32),dimension(:),allocatable :: c_r4
complex(kind=real64),dimension(:),allocatable :: c_r8
complex(kind=real128),dimension(:),allocatable :: c_r16

c_r4 = [(1,2),(2,4)]
c_r8 = c_r4
c_r16 = c_r4

n = 1
ntests = 13
call initialize_tests(tests,ntests)

tests(n) = assert(all(tail([1_int8,2_int8]) == [2_int8]),'tail, int8')
n = n + 1

tests(n) = assert(all(tail([1_int16,2_int16]) == [2_int16]),'tail, int16')
n = n + 1

tests(n) = assert(all(tail([1_int32,2_int32]) == [2_int32]),'tail, int32')
n = n + 1

tests(n) = assert(all(tail([1_int64,2_int64]) == [2_int64]),'tail, int64')
n = n + 1

tests(n) = assert(all(tail([1._real32,2._real32]) == [2._real32]),'tail, real32')
n = n + 1

tests(n) = assert(all(tail([1._real64,2._real64]) == [2._real64]),'tail, real64')
n = n + 1

tests(n) = assert(all(tail([1._real128,2._real128]) == [2._real128]),'tail, real128')
n = n + 1

tests(n) = assert(all(tail(c_r4) == [c_r4(2)]),'tail, complex real32')
n = n + 1

tests(n) = assert(all(tail(c_r8) == [c_r8(2)]),'tail, complex real64')
n = n + 1

tests(n) = assert(all(tail(c_r16) == [c_r16(2)]),'tail, complex real128')
n = n + 1

tests(n) = assert(size(tail([1._real32])) == 0,'size(tail([1])) == 0')
n = n + 1

tests(n) = assert(size(tail(tail([1._real32]))) == 0,'size(tail(tail([1]))) == 0')
n = n + 1

tests(n) = assert(all(tail([1,2]) == .tail.[1,2]),'tail operator, .tail.x')
n = n + 1

test_failed = .false.
call report_tests(tests,test_failed)
if(test_failed)stop 1

endprogram test_tail

max_stars_repo_name: Gorkowski/functional-fortran

max_issues_repo_name: Gorkowski/functional-fortran

max_forks_repo_name: Gorkowski/functional-fortran


module OpenSpecModule
!  Code in this file defines values for OPEN-statement specifiers.  Some
!  of the values are extensions to ANSI Fortran 90 and 95.  One of the
!  specifiers is not included in ANSI FORTRAN 77. The included
!  specifiers are ACCESS, FORM and ACTION.
!
      CHARACTER(len=20) :: ACCESS,FORM,ACTION(2)
!
!
!  Specifiers for OPEN statements for unformatted files, which are
!  sometimes compiler specific.
!  The included specifiers are ACCESS and FORM.
!
!  ACCESS specifier --
!
!    Standard Fortran -- Use unless there is a reason to do otherwise.
!      DATA ACCESS/'SEQUENTIAL'/
       DATA ACCESS/'STREAM'/
!
!
!  FORM specifier --
!
!    Standard Fortran, which results in vendor dependent (non-portable)
!    files.  Use unless there is a reason to do otherwise.
      DATA FORM/'UNFORMATTED'/
!
!    Non-standard Fortran that causes code compiled by Compaq (Digital)
!    Fortran on personal computers to use unstructured non-formatted
!    files.  This may make it possible for the non-formatted files used
!    by MODFLOW to be used with programs that are compiled by other
!    compilers.
!      DATA FORM/'BINARY'/
!
!
!  OPEN-statement specifiers related to file-sharing.
!
!  ACTION specifier --
!
!    Standard FORTRAN 77 -- Eliminate the ACTION= specifier from all
!    OPEN statements in the source-code files.
!
!    Standard Fortran 90 and 95 -- Use unless there is a reason to do
!    otherwise.
      DATA (ACTION(IACT),IACT=1,2)/'READ','READWRITE'/
!
!    Non-standard Fortran that causes code compiled by the Lahey LF90
!    compiler to create files that can be shared.  For use when parallel
!    processing is used or to enable an editor to view output files
!    while the program is running.
!      DATA (ACTION(I),I=1,2)/'READ,DENYWRITE','READWRITE,DENYNONE'/
!
end module OpenSpecModule


program test_filter
    use fplus_filter
    implicit none

    ! local variables
    integer :: i, f
    real (kind=8) :: pi

    ! array for test-data
    integer, parameter :: n = 10000
    real (kind=8), dimension(n) :: tdata

    ! fill the array with data
    print*, "4th Order Butterworth filter with cutoff periods 100 and 200:"
    pi = ACOS(-1.)
    do f = 1, 100
        do i = 1, n
            tdata(i) = sin(pi/2*i/real(f))
        end do
        ! perform the filtering
        call BWBPASS(tdata,n,100.0_8,200.0_8,2)
        print*, "Period: ", f*4, "Sum: ", sum(abs(tdata))
    end do

end program


text:
OPEN(3,FILE='FSEQ')
!  FSEQ is a formatted sequential fife by default.
WRITE(3,'(A,I3)') 'RECORD',1
WRITE(3, '()')
WRITE(3, '(A11) ') 'The 3rd one'
CLOSE(3)
END

max_stars_repo_name: quchunguang/test


Module declare

   IMPLICIT NONE

   INTEGER(KIND=8)::TIMESTART,TIMERUN,count_rate,cr,jstep,jstart(1000)
   REAL(KIND=8)::LASTTIME,rate

   INTEGER NATOMOLD,NATOM,NATOMNEW,NATOMMAX,NSYN,NSYNMAX,NPG,NPGMAX,NBONDGLY,NGLYMAX,NBONDPEP,NPEPMAX,NLOOP,NLOOPMAX
   INTEGER NPGOLD,NPGNEW,CAPLEN,NCAP,NBONDCAP,NBONDGLYOLD!,CAPBOUND1,CAPBOUND2
   INTEGER SYNRATIO,NATOMCAP,NB1,NB2,JRESET,LENMAX,newlengthave,length
   INTEGER NBOND,NBONDDEL,NBONDGLYSTART,NSIG
   INTEGER NL,NLOOPDEL,NB,PGCAP1,PGCAP2
   INTEGER PGLENMAX,LOOPLENMAX,JLYT,JRECOMP
   INTEGER TIMES(10),YEARS,MONS,DAYS,HOURS,MINS,SECS
   INTEGER N,NTHREADS!,OMP_GET_NUM_THREADS
   INTEGER NSTART,NR,JFORCE1,JFORCE2,SLOWLOAD
   INTEGER NTIP,ntip1,ntip2,NSEC,NPICK,ILOOP,JLINK,NNEW,JRFORCE,NRAND
   INTEGER NS,NS1,NPEPTOTAL,NATOMSTART,NPGSTART,NGLYNEW,ns_ter,jsyndir,na_red
   INTEGER NATOMDEL,OLDNATOMDEL,NDEL,N1,N2,NA,nstore,npress
   INTEGER JR,NTOTAL,JFILE,NFRAME,JL,JS,NTRAP,JCHANGE,JSCALE,A22,JSKIP,SKIP
   INTEGER IPG,JUNIT,NSYNREP,NTRANS,NHOLD,NRELAX,J1,J2,JEXIT,JPRINT,JRANDS,N_EDHOLD,noldbond,jpgnew,noldmax
   integer oldtypmax,jredist,jrestart,jsteric,ipg1,ipg2,jfast1,jfast2,jprevent,joldbond,nnewhoop
   integer pgbelowbondmax,ntemlink,ntemlinkmax,nprint

   INTEGER(kind=8),DIMENSION(:),ALLOCATABLE::jtrap,J_EDHOLD
   INTEGER,DIMENSION(:),ALLOCATABLE:: EDOSIG,TP2SIG,EDCAP,JDEACT,COUPLING
   INTEGER,DIMENSION(:,:),ALLOCATABLE::PGID,newpgid,EDHOLD,SIGBOND,EDPEP,JTERM
   INTEGER,DIMENSION(:,:),ALLOCATABLE::NEWPGIDtemp,TRANSGLY,GLYTIP,GLYSEC,oldglytip,oldglysec,temlink
   INTEGER,DIMENSION(:),ALLOCATABLE::PGLEN,newpglentemp,PGDIR,DNOR,ATOR,PEPDIR,glyid
   INTEGER,DIMENSION(:),ALLOCATABLE::NEWPGLEN,NEWPGTYP,NEWPGDIR,ATORNEW,DNORNEW,PEPDIRNEW
   INTEGER,DIMENSION(:,:),ALLOCATABLE::BOND,BONDGLY,BONDPEP,oldbond
   INTEGER,DIMENSION(:,:),ALLOCATABLE::NEWBONDGLY,NEWBONDPEP,SIGCROSS,ADDSIG
   INTEGER,DIMENSION(:),ALLOCATABLE::BONTYP,oldtyp,pgbelowbond
   INTEGER,DIMENSION(:),ALLOCATABLE::NEWBONTYP,SYNMARK,JR_ENDO_ON,JR_ENDO_OFF,SIGCLEAVE,EDLOCKIN
   INTEGER,DIMENSION(:,:),ALLOCATABLE::LOOP,NEWLOOP,PARTNER
   INTEGER,DIMENSION(:,:),ALLOCATABLE::CRLKAGE,NETWORK,SYNTHESIS,SYNPG
   INTEGER,DIMENSION(:),ALLOCATABLE::LOOPLEN,NEWLOOPLEN,LOOPTYP,NEWLOOPTYP
   INTEGER,DIMENSION(:),ALLOCATABLE:: SYNDIR,SYNLOOP,GTLOAD,GTATRANS
   INTEGER,DIMENSION(:),ALLOCATABLE:: NEINUM,GLYNUM
   INTEGER,DIMENSION(:,:,:),ALLOCATABLE::NEIBOND,NEIGLY,TPPEP

!   INTEGER,DIMENSION(:),ALLOCATABLE:: temlen,temtyp
!   INTEGER,DIMENSION(:,:),ALLOCATABLE::temid
!   integer ntem

   DOUBLE PRECISION PI!= 3.141592653589793239D0
   DOUBLE PRECISION L_G,K_G,L_P,K_P,THETA_0,KTHETA,PRES,P_COUPL,P_UCOUP,MreBlen0
   DOUBLE PRECISION P_ACTI,P_DEACT,PTERM0,pterm,P_GLYIN,P_PEPIN,P_PEPOUT,P_EDHOLD,PCLEAVE_SPON,PEDCAP_SPON,PEDCAP_INACT
   DOUBLE PRECISION PMATURE,PLYT,GROWTH,P_GLYOUT,PBEND,P_TRANSFAST,P_TRANSSLOW,prmbond,wthick,pthick,kaxis
   DOUBLE PRECISION p_trans,lambda

   DOUBLE PRECISION sigma,kmidcell,fmid,fconstrict,pconstrict0!,pconstrict(1000),pconstrictold(1000)
   DOUBLE PRECISION delx,invdelx,dpcons,xpbp,ftotal,lfinal,lcons

   DOUBLE PRECISION XCAP1,XCAP2,R,GAM,FAMAX,DREACT,DREACT2,DIST
   DOUBLE PRECISION FPMAX1,FPMAX2,FPMAX3,FPMAX4,FPMAX5,FPMAX6
   DOUBLE PRECISION DX,DY,DZ,FPBPMAX,INVMPBP,AREA,ORAD
   DOUBLE PRECISION LSTART,LREP,LSWITCH,KSWITCH,MSWITCH
   DOUBLE PRECISION KGTASE,KTPASE,LTPASE,KPAIR,LPAIR,KGTTP,LGTTP13,LGTTP2,l_edcap
   DOUBLE PRECISION KEDASE,LEDASE,KGTED,LGTED,KWALL,KLEAD,DELTA,INVDELTA,BETA,KSUR,KSIDE,LSIDE,ktped,ltped

   real(kind=8)::ratio(1000),newratio

   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:: X,Y,Z,XEDASEOLD,YEDASEOLD,ZEDASEOLD
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:: XNEW,YNEW,ZNEW,XLEAD,YLEAD,ZLEAD,XEDASE,YEDASE,ZEDASE
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FX,FY,FZ,FAMAG,FXGLY,FYGLY,FZGLY,FXSUR,FYSUR,FZSUR
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FXPEP,FYPEP,FZPEP,FXTHETA,FYTHETA,FZTHETA
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FXREP,FYREP,FZREP
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FXPRES,FYPRES,FZPRES,FXSYN,FYSYN,FZSYN
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FXCAPBI1,FXCAPBI2,FXCAPMON,FXCAPEDO
   DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE::XGTASE,YGTASE,ZGTASE,FXGTASE,FYGTASE,FZGTASE,LGTASE
   DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE::FXTPASE,FYTPASE,FZTPASE,XTPASE,YTPASE,ZTPASE
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::FXEDASE,FYEDASE,FZEDASE,FMAGPBP
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::R1,R2,EDRAD,BIRAD1,BIRAD2,MORAD,RANDS
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::RXGTASE,RYGTASE,RZGTASE,RXTPASE,RYTPASE,RZTPASE
   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::RXEDASE,RYEDASE,RZEDASE,fycon,fzcon

   DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE:: XCEN,YCEN,ZCEN,UX,UY,UZ,MreBlen,P_CURV,P_UCURV
   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)::YCENBI1,ZCENBI1,YCENBI2,ZCENBI2
   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)::YCENMON,ZCENMON,YCENEDO,ZCENEDO,randstore

   DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:)::xsyn,pconstrict,pconstrictold

   REAL,DIMENSION(:),ALLOCATABLE::XW,YW,ZW
   character (len=512) command

end module


!---------------------------------------------------!
! Copyright (c) 2016 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
!-------10--------20--------30--------40--------50--------60--------70--------80--------90
subroutine write_td_results
  use global_variables
  implicit none
  integer :: it

  open(20,file="td_result.out")
  do it = 0,Nt_iter
    write(20,"(999e26.16e3)")it*dt,dipole_t(it),dipole_t(it)-dipole_t(0) &
      ,quadrupole_t(it),quadrupole_t(it)-quadrupole_t(0),norm_t(it)
  end do
  close(20)

  return
end subroutine write_td_results


text:
      MODULE MOD_DECSTE

      CONTAINS

      SUBROUTINE DECSTE(LSRAT,LSPOP,JOBMAX,EPSILON,REDUCE,IHIST,IFRRA,ITORA,LSEXPO,
     $                  IFLUX,IDAT,LEVELPL,N,IPLOTF,NEWWRC,
     $                  NGAMR,NGAML,AGAMR,AGAML,DELTAC,LINE,MAXIND,
     $                  TPLOT,Y0,TEFFE,GRAD,ALDMDT,VINF,BET,PROLIB,LBLANK)

      USE COMMON_BLOCK

!     DECODING INPUT OPTIONS FOR STEAL.FOR

      IMPLICIT REAL*8(A - H, O - Z)

      CHARACTER KARTE*80
      DIMENSION LEVELPL(N)
      INTEGER   NGAMR(10),NGAML(10)
      DIMENSION AGAMR(10),AGAML(10)
      LOGICAL   LINE (MAXIND),PROLIB,NOTEMP,TPLOT
C***  COMMON / COMPLOT / TRANSFERS THE MAXIMUM OF THE Y-AXIS TO SUBR. PLOTT
      COMMON / COMPLOT / YMAX

C***  NO PRINT OPTIONS SET (DEFAULT)
      TEFFE=-1.
      LSRAT=-1
      LSPOP=-1
      LSEXPO=-1
      IHIST=-1
      IFLUX=-1
      IDAT=-1
      LBLANK=0
C***  OLDSTART=.F.:  POPNUMBERS ARE CALCULATED IN SUBR. POPZERO (DEFAULT)
C***  OLDSTART=.T.:  POPNUMBERS FROM OLD MODEL FILE
      OLDSTART=.FALSE.
C***  JOBMAX EXCEEDED (DEFAULT)
      JOBMAX=-1
C***  DEMANDED CONVERGENCE ACCURACY (DEFAULT)
      EPSILON=.005
C*** NO REDUCTION OF CORRECTIONS (DEFAULT)
      REDUCE = 1.
C***  NO POPNUMBER PLOT (DEFAULT)
      LEVELPL(1)=0
      NPLOT=0
C***  NO FLUX PLOT (DEFAULT)
      IPLOTF=-1
      PROLIB=.FALSE.
C***  NUMBER OF REPEAT ITERATIONS BETWEEN TWO WRCONT-JOBS (DEFAULT)
      NEWWRC=11
C***  NO SCHARMER AMPLIFICATION OF CORRECTIONS
      DO 21 I=1,10
      NGAMR(I)=999999999
      NGAML(I)=999999999
      AGAMR(I)=0.
      AGAML(I)=0.
   21 CONTINUE
      NGAMR(1)=0
      NGAML(1)=0
      MGC=1
      MGR=1
      MGL=1
      DELTAC=1.0
      IPRICC=0
      IPRILC=0
C***  LINES THAT HAVE NOT BEEN TREATED IN THE RADIATION TRANSFER CALCULATIONS,
C***  I.E. WHICH ARE NOT QUOTED IN A 'LINE' OPTION CARD, ARE NOT SCHARMER-
C*** AMPLIFIED.
      DO 32 IN=1,MAXIND
      LINE(IN)=.FALSE.
   32 CONTINUE
C***  DISTANCE MODULUS (DEFAULT VALUE)
      Y0=.0
C***  NO SUPPRESSION OF TEMPERATURE CORRECTIONS (DEFAULT)
c*** when using the temperature correction procedurte:
c    this option should be switched off within steal
      NOTEMP=.TRUE.
C*** NO TEMPERATURE PLOT (DEFAILT)
      TPLOT=.FALSE.

      OPEN (1,FILE='CARDS')
      REWIND 1
    8 READ (1,4,END=66) KARTE
    4 FORMAT (A)

      IF ( KARTE(:10) .EQ. 'LINE BLANK') THEN
C                           ==========
            IF (KARTE(17:18).EQ.'OF') THEN
               LBLANK=0
               PRINT *,' STEAL: LINE BLANKETING OPTION OFF'
            ELSE IF (KARTE(17:18).EQ.'ON') THEN
               LBLANK=1
               PRINT *,' STEAL: LINE BLANKETING OPTION ON'
            ELSE IF (KARTE(17:18).EQ.'PR') THEN
               LBLANK=2
               PRINT *,' STEAL: LINE BLANKETING OPTION ON (print on)'
            ELSE IF (KARTE(17:18).EQ.'IN') THEN
               LBLANK=-1
               PRINT *,' STEAL: LB TABLE INPUT'
            ELSE
               PRINT *,' NO VALID OPTION'
            ENDIF
            GOTO 8
            ENDIF

      IF (KARTE(:6) .EQ. 'NEWWRC' ) THEN
C                         ======
            DECODE (80,9,KARTE) XL
    9       FORMAT (7X,F10.0)
            NEWWRC=NINT(XL)
            GOTO 8
            ENDIF
      IF (KARTE(:9) .EQ. 'PLOT FLUX') THEN
C                         =========
            IPLOTF=1
            GOTO 8
            ENDIF
      IF ( KARTE(:6) .EQ. 'PROLIB' ) THEN
C                          ======
         DECODE(80,13,KARTE) TEFFE,GRAD,ALDMDT,VINF,BET,Y0
   13    FORMAT(10X,6F10.3)
            IPLOTF=1
            PROLIB=.TRUE.
            GOTO 8
            ENDIF
      IF (KARTE(:10) .EQ. 'PRINT RATE' ) THEN
                          !===========
            DECODE (80,1,KARTE) XL,XFRRA,XTORA

!0        1        2         3         4         5
!23456789 12345678 123456789 123456789 123456789 1234567890
!RINT RATE IIIIIIIII         IIIIIIIIII        IIIIIIIIII
    1       FORMAT (11X,F10.0,10X,F10.0,8X,F10.0)
            LSRAT=NINT(XL)
            IFRRA=NINT(XFRRA)
            ITORA=NINT(XTORA)
            IF (LSRAT.EQ.0) LSRAT=1
            IF (IFRRA.GT.ITORA) ITORA=IFRRA
            PRINT '("DECSTE: RATE =",I5,":",I5,":",I5)', IFRRA,LSRAT,ITORA
            GOTO 8
            ENDIF
      IF (KARTE(:10) .EQ. 'PRINT HIST' ) THEN
C                          ==========
            IHIST=1
            GOTO 8
            ENDIF
      IF (KARTE(:10) .EQ. 'PRINT FLUX' ) THEN
C                          ==========
            IFLUX=1
            GOTO 8
            ENDIF
      IF (KARTE(:11) .EQ. 'PRINT DATOM' ) THEN
C                          ===========
            IDAT=1
            GOTO 8
            ENDIF
      IF (KARTE(:9) .EQ. 'PRINT POP' ) THEN
C                         =========
            DECODE (80,3,KARTE) XL
    3       FORMAT (9X,F10.0)
            LSPOP=NINT(XL)
            IF (LSPOP.EQ.0) LSPOP=1
            PRINT *,' PRINT POP ',LSPOP,' DECODED'
            GOTO 8
            ENDIF
      IF (KARTE(:10) .EQ. 'PRINT EXPO' ) THEN
C                          ==========
            DECODE (80,1,KARTE) XL
            LSEXPO=NINT(XL)
            IF (LSEXPO .EQ. 0) LSEXPO=1
            PRINT *,' PRINT EXPO ',LSEXPO,' DECODED'
            GOTO 8
            ENDIF
      IF (KARTE(:7) .EQ. 'JOBMAX=' ) THEN
C                         =======
            DECODE (80,11,KARTE) XL
   11       FORMAT (7X,F10.9)
            JOBMAX=NINT(XL)
!      IF (JOBMAX .GE. 1000) THEN
!            PRINT *,'JOBMAX .GT. 1000: INVALID OPTION'
!            STOP 'ERROR'
!            ENDIF
            GOTO 8
            ENDIF
      IF (KARTE(:7) .EQ. 'EPSILON' ) THEN
C                         =======
            DECODE (80,14,KARTE) EPSILON
   14       FORMAT (8X,F10.0)
            GOTO 8
            ENDIF
      IF (KARTE(:6) .EQ. 'REDUCE' ) THEN
C                         ======
            DECODE (80,17,KARTE) REDUCE
   17       FORMAT (7X,F10.0)
            IF (REDUCE.EQ..0) REDUCE=.5
            GOTO 8
            ENDIF
      IF (KARTE(:8) .EQ. 'PLOT POP') THEN
C                         ========
            IF (NPLOT+1 .EQ. N) GOTO 8
            DECODE (30,5,KARTE) L1,L2,L3
    5 FORMAT (9X,3(I2,1X))
            IF (L1 .LE. 0) GOTO 8
            NPLOT=NPLOT+1
            ENCODE (8,2,LEVELPL(NPLOT)) L1,L2,L3
    2       FORMAT (3I2)
            LEVELPL(NPLOT+1)=0
            GOTO 8
            ENDIF
      IF (KARTE(:11) .EQ. 'PRINT CCORE') THEN
C                          ===========
            IPRICC=1
            GOTO 8
            ENDIF
      IF (KARTE(:11) .EQ. 'PRINT LCORE') THEN
C                          ===========
            DECODE (80,7,KARTE) XLC
    7       FORMAT (11X,F10.0)
            IPRILC=NINT(XLC)
            GOTO 8
            ENDIF
      IF (KARTE(:6) .EQ. 'GAMMA=') THEN
C                         ======
            MGC=MGC+1
            MGR=MGR+1
            MGL=MGL+1
            IF (MGC .GT. 10 .OR. MGR .GT. 10 .OR. MGL .GT. 10) THEN
               write (6,*) ('TOO MANY GAMMA OPTIONS')
               STOP 'ERROR'
               ENDIF
            DECODE (80,12,KARTE) GAMMA, ANG
   12       FORMAT (6X,F11.0,9X,F10.0)
            NGAMR(MGR)=ANG
            NGAML(MGL)=ANG
            AGAMR(MGR)=GAMMA
            AGAML(MGL)=GAMMA
            GOTO 8
            ENDIF
C      IF (KARTE(:6) .EQ. 'GAMMAC' ) THEN
C                         ======
C            MGC=MGC+1
C            IF (MGC .GT. 10) STOP 'ERROR'
C            DECODE (80,6,KARTE) GAMMAC,ANGC
C            GOTO 8
C            ENDIF
    6       FORMAT (7X,F10.0,9X,F10.0)
      IF (KARTE(:6) .EQ. 'GAMMAR' ) THEN
C                         ======
            MGR=MGR+1
            IF (MGR .GT. 10) STOP 'ERROR'
            DECODE (80,6,KARTE) GAMMAR,ANGR
            NGAMR(MGR)=ANGR
            AGAMR(MGR)=GAMMAR
            GOTO 8
            ENDIF
      IF (KARTE(:6) .EQ. 'GAMMAL' ) THEN
C                         ======
            MGL=MGL+1
            IF (MGL .GT. 10) STOP 'ERROR'
            DECODE (80,6,KARTE) GAMMAL,ANGL
            NGAML(MGL)=ANGL
            AGAML(MGL)=GAMMAL
            GOTO 8
            ENDIF
      IF (KARTE(:6) .EQ. 'DELTAC') THEN
C                         ======
            DECODE (80,10,KARTE) DELTAC
   10       FORMAT (7X,F10.0)
            GOTO 8
            ENDIF
      IF ( KARTE(:4) .EQ. 'LINE' ) THEN
C                          ====
            DECODE (80,41,KARTE) IND1
   41       FORMAT(4X,I3)
            IF (KARTE(9:10) .EQ. 'TO') THEN
                  DECODE (80,43,KARTE) IND2
   43             FORMAT(15X,I3)
                  ELSE
                  IND2=IND1
                  ENDIF
            IF (IND2 .GE. IND1) THEN
                  INC=1
                  ELSE
                  INC=-1
                  ENDIF
            DO 45 IND=IND1,IND2,INC
            IF (IND .GT. MAXIND) THEN
               PRINT *,' ERROR STOP IN DECSTE:  IND  .GT. MAXIND'
               write (6,*) (' IND  .GT. MAXIND')
               STOP 'ERROR'
               ENDIF
            LINE(IND)=.TRUE.
   45       CONTINUE
            GOTO 8
            ENDIF
      IF (KARTE(:8) .EQ. 'OLDSTART') THEN
C                         ========
            OLDSTART=.TRUE.
            GOTO 8
            ENDIF
      IF (KARTE(:8) .EQ. 'NO TEMPE') THEN
C                         ========
         NOTEMP=.TRUE.
         GOTO 8
         ENDIF
      IF (KARTE(:6).EQ.'PLOT T') THEN
C                       ======
            DECODE (80,15,KARTE) YMAX
   15       FORMAT (22X,F6.0)
            TPLOT=.TRUE.
         GOTO 8
            ENDIF

      GOTO 8

   66 CONTINUE
      CLOSE (1)
      RETURN

      END subroutine
      end module

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


text:
C************************************************************************
C
C                           HISTOGRAMMING PACKAGE
C                                 M. MANGANO
C
C************************************************************************
C
C**********************************************************************
C    SIMPLE HISTOGRAMMING PACKAGE --  SIMPLIFIED VERSION OF HBOOK
C    BY Michelangelo Mangano    NOVEMBER 1988
C    LAST REVISED NOVEMBER 9, 1988
C    LAST REVISED JUNE 12, 1989  (ADD SCATTER PLOTS)
C    LAST REVISED oct 1990 (Add multi-plots on one page, routines MULTITOP,
C         			MTFILL,...)
C    LAST REVISED Jun 2003 by FM: passed to double precision AND
C    included the date for linux
C**********************************************************************
C
C Fills up to 500 histograms with up to 100 bins.
C Gives a data file (to be specified in the calling program by assigning
C a file name to unit 98) and a topdrawer file (to be specified in the
C calling program by assigning a file name to unit 99).
C
C INITIALIZATION:
C Call once INIHIST; this just resets a few counters and logicals
C Call MBOOK(N,'TITLE',DEL,XMIN,XMAX) for each histogram to be booked.
C N (an integer) is the label of the histogram;
C 'TITLE' is the name of the histogram (no more then 100 characters);
C DEL (real*8) is the bin size;
C XMIN (real*8) is the lower limit of the first bin;
C XMAX (real*8)is the upper limit of the last  bin
C Example:
C      call mbook(2,'pt distribution',1.,10,70)
C This call initializes histogram number 2, called 'pt distribution';
C The bin size will be 1. (possibly GeV, if that's what you want), the
C first bin being  10<x<11. and the last one being 69.<x<70
C
C FILLING:
C When it's time, call MFILL(N,X,Y); this will add Y (real*8) to the bin
C in which X (real*8) happens to be, within histogram N.
C
C PLAYING AROUND:
C At the end of the day you may want to sum, divide, cancel, etc.etc.
C various histograms (bin by bin). Then you call MOPERA(I,'O',J,K,X,Y).
C The 1-character string O can take the following values:
C +  : sums       X*(hist I) with Y*(hist J) and puts the result in hist K;
C -  : subtracts  X*(hist I) with Y*(hist J) and puts the result in hist K;
C *  : multiplies X*(hist I) with Y*(hist J) and puts the result in hist K;
C /  : divides    X*(hist I) with Y*(hist J) and puts the result in hist K;
C F  : multiplies hist I by the factor X, and puts the result in hist K;
C R  : takes the square root of  hist  I, and puts the result in hist K;if
C      the value at a given bin is less than or equal to 0, puts 0 in K
C S  : takes the square      of  hist  I, and puts the result in hist K;
C L  : takes the log_10 of  hist  I, and puts the result in hist K; if the
C      value at a given bin is less than or equal to 0, puts 0 in K
C M  : statistical analysis; if I contains the weights (let's say WGT),
C      J contains variable times weight (F*WGT) and K contains the
C      variable squared times the weight (F**2*WGT), then, after using 'M',
C      J will contain the average value of the variable <F> and K will
C      contain the sigma of the average: sigma=sqrt(<F**2>-<F>**2).
C      If WGT=1. for all the entries, then it is enough to put I=J, and
C      it is not necessary to book a hist with the weights.
C V  : estimates errors for vegas evaluation of differential distributions.
C      Fill I with the values of
C      the functions do integrate times the Vegas weight (fun*wgt); fill
C      J with fun**2*wgt; then K will contain an estimate of the error
C      of the integration. Putting X=1/(#of iterations) performs the
C      avegare over the iterations, and gives the right normalization to
C      the differential distribution, I, and to the errors, K. J stays the same.
C
C FINAL ACCOUNTING:
C Now we can finalize our histograms; MFINAL(N) will calculate the integral
C of the histogram N, the mean value of the X variable and its RMS.
C If we now want to renormalize the hist's, we can call MNORM(N,X), which
C will normalize the integral to X  -- CAUTION: do not call MNORM before
C MFINAL, it will blow up.
C
C OUTPUT:
C To get a .dat file containing the values of the histograms, together with
C some information (like integral, mean values, etc.etc.) call MPRINT(N),
C for each hist N that you want in the .dat file. Before the call to MPRINT
C you want to open unit 98 and give it a name:
C     OPEN(UNIT=98,NAME='NAME.DAT',STATUS='NEW')
C If you want a topdrawer file with a plot of the hist values, call
C MTOP(N,M,'X','Y','SCALE'). The points of the plot will be taken from histogram
C N, the error bars from histogram M. 'SCALE', character*(*), determines
C the scale for y, logarithmic or linear (SCALE=LOG,LIN).
C If you do not want error bars, keep
C a histogram of zeros, or just call a hist that had not been booked.
C X will appear as a 'bottom title', and Y will appear as a 'left title'.
C The top title is by default the name of the histogram itself.
C A little box below the plot will contain some information on the plot
C itself. Before calling MTOP,
C     OPEN(UNIT=99,NAME='NAME.TOP',STATUS='NEW')
C--------------------------------------------------------------------------
C
C  COMMON/HISTO/  Histogram N
C
C   BOOK(N),      Three-letter character-string: 'NO' if histogram was not
C		  Booked, 'YES' otherwise
C   TITLE(N),     Title of the histogram
C
C   HMIN(N),      Min value of x range
C   HMAX(N),      Max value of x range
C   HDEL(N),      Bin width
C   NBIN(N),      Total number of bins
C   USCORE(N),    Total integral of underscores with x < HMIN(N)
C   OSCORE(N),    Total integral of onderscores with x > HMAX(N)
C   IUSCORE(N),   Number of entries with x < HMIN(N)
C   IOSCORE(N),   Number of entries with x > HMAX(N)
C   IENT(N),      Total number of entries within x range HMIN(N)<x<HMAX(N)
C   HINT(N),      Integral of the histogram within HMIN(N)<x<HMAX(N)
C   HAVG(N),      Average value of x, weighted over the x range of the histo
C   HSIG(N),      Quadratic dispersion of x around the average
C   HIST(N,L),    Value of bin L-th
C   XHIS(N,L),    Central x value of bin L-th
C   IHIS(N,L),    Number of entries within bin L-th
C   NHIST         Total number of booked histograms
C

      SUBROUTINE INIHIST
C*************************************************************
c     initialization
      IMPLICIT NONE
      include 'dbook.inc'
c     APPLgrid commons
      include "reweight_appl.inc"
      include "appl_common.inc"
      integer iappl
      common /for_applgrid/ iappl
c
c     LOCAL
c
      INTEGER I
C
C     START
C
      NHIST=0
      DO 1, I=1,NPLOTS
   1  BOOK(I)=' NO'
C     Initialize the number of bins of the aMCfast grids
      if(iappl.ne.0) appl_obs_nbins = 0
      END

      SUBROUTINE MBOOK(N,TIT,DEL,XMIN,XMAX)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
      CHARACTER*(*) TIT
      REAL*8 DEL,XMIN,XMAX
C
C     GLOBAL
C
      include 'dbook.inc'
C
C     LOCAL
C
      INTEGER I
Cq
C     START
C
      NHIST = MAX(N,NHIST)
      IF(BOOK(N)(1:1).EQ.'Y') THEN
	 CALL MWARN('MBOOK')
         WRITE(*,*) 'Histogram',N,TITLE(N),' already in use. '
         WRITE(*,*) 'superseded by ',TIT
      ENDIF
      BOOK(N) = 'YES'
      TITLE(N) = ' '//TIT
1     HDEL(N) = DEL
      NBIN(N) = INT((XMAX-XMIN)/(DEL*0.999999d0))
      IF(NBIN(N).GT.100) THEN
	WRITE(*,*) 'TOO MANY BINS (',NBIN(N),') REQUIRED IN HIST ',N
	WRITE(*,*) 'RE-ENTER BIN SIZE DELTA (OLD BIN = ',DEL,' ):'
	READ(*,*) DEL
	GO TO 1
      ENDIF
      HMIN(N) = XMIN
      HMAX(N) = NBIN(N)*DEL+XMIN
      IF(abs(HMAX(N)-XMAX).gt.0.001d0*DEL) THEN
	 CALL MWARN('MBOOK')
         WRITE(*,*)
     #'Histogram ', TIT, ' Change of upper limit:',xmax,'-->',HMAX(N)
      ENDIF
      IENT(N) = 0
      IUSCORE(N) = 0
      IOSCORE(N) = 0
      USCORE(N) = 0
      OSCORE(N) = 0
      HAVG(N) = 0
      HINT(N) = 0
      HSIG(N) = 0
      DO I=1,NBIN(N)
         XHIS(N,I)=HMIN(N)+HDEL(N)*(DFLOAT(I)-0.5d0)
         IHIS(N,I)=0
         HIST(N,I)=0
      ENDDO
      END

      SUBROUTINE MFILL4(N,X,Y)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
      REAL*8 X,Y
C
C     LOCAL
C
      REAL*8 XI
      INTEGER I
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
c
c     START
c

      IF(X.LT.HMIN(N)) THEN
         USCORE(N) = USCORE(N) + Y
         IUSCORE(N) = IUSCORE(N) + 1
      ELSEIF(X.GT.HMAX(N)) THEN
         OSCORE(N) = OSCORE(N) + Y
         IOSCORE(N) = IOSCORE(N) + 1
      ELSE
         XI=((X-HMIN(N))/HDEL(N))+1
         I=INT(XI)
         IENT(N)=IENT(N)+1
         IHIS(N,I)=IHIS(N,I)+1
         HIST(N,I)=HIST(N,I)+Y
      ENDIF
      END

      SUBROUTINE MINTEG(NIN,NOUT,IDIR,IPOW)
C If IPOW=1 performs the integral of the distribution contained in histogram
C NIN up to the value specified by the abscissa (if IDIR=1) or from this
C value on (if IDIR=-1). The resulting integral distribution is put into
C NOUT, which is automatically booked if NOUT.ne.NIN .  Choosing IPOW=2
C the routine will return the square root of the integral of the squares,
C as is required, for example, for the propagation of the mean quadratic error
C of a given distribution. Overscores and underscores are included.
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER NIN,NOUT,IDIR,IPOW
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
c
c     LOCAL
c
      INTEGER M,I,L
      CHARACTER*14  C
      DIMENSION C(2)
      DATA C/' INTEG BELOW X',' INTEG ABOVE X'/
C
C     EXTERNAL
C
      REAL*8 SUMPOW
C
C     START
C
      IF(BOOK(NIN)(1:1).NE.'Y') RETURN
      M = NBIN(NIN)
      I = (IDIR + 3)/2
      IF(NOUT.NE.NIN.AND.BOOK(NOUT)(1:1).NE.'Y') THEN
      	CALL MBOOK(NOUT,TITLE(NIN)//C(I),
     &                HDEL(NIN),HMIN(NIN),HMAX(NIN))
      ENDIF
      IF(IDIR.EQ.1) THEN
         HIST(NOUT,1) = SUMPOW(HIST(NIN,1),USCORE(NIN),IPOW)
         IHIS(NOUT,1) = IHIS(NIN,1) + IUSCORE(NIN)
         XHIS(NOUT,1) = XHIS(NIN,1) + HDEL(NIN)/2
         DO L=2,M
            HIST(NOUT,L) = SUMPOW(HIST(NIN,L),HIST(NOUT,L-1),IPOW)
            IHIS(NOUT,L) = IHIS(NIN,L) + IHIS(NOUT,L-1)
            XHIS(NOUT,L) = XHIS(NIN,L) + HDEL(NIN)/2
         ENDDO
         OSCORE(NOUT) = SUMPOW(OSCORE(NIN),HIST(NIN,M),IPOW)
         IOSCORE(NOUT) = IOSCORE(NIN) + IHIS(NIN,M)
      ELSEIF(IDIR.EQ.-1) THEN
         HIST(NOUT,M) = SUMPOW(HIST(NIN,M),OSCORE(NIN),IPOW)
         IHIS(NOUT,M) = IHIS(NIN,M) + IOSCORE(NIN)
         XHIS(NOUT,M) = XHIS(NIN,M) - HDEL(NIN)/2
         DO L=M-1,1,-1
            HIST(NOUT,L) = SUMPOW(HIST(NIN,L),HIST(NOUT,L+1),IPOW)
            IHIS(NOUT,L) = IHIS(NIN,L) + IHIS(NOUT,L+1)
            XHIS(NOUT,L) = XHIS(NIN,L) - HDEL(NIN)/2
         ENDDO
         USCORE(NOUT) = SUMPOW(USCORE(NIN),HIST(NIN,1),IPOW)
         IUSCORE(NOUT) = IUSCORE(NIN)+IHIS(NIN,1)
      ELSE
         CALL MWARN('MINTEG')
         WRITE(*,*) 'Wrong idir in minteg: OPERATION NOT PERFORMED'
         RETURN
      ENDIF
      END

      REAL*8 FUNCTION SUMPOW(X,Y,IPOW)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      REAL*8 X,Y
      INTEGER IPOW
C
C     START
C
      IF(IPOW.EQ.1) THEN
         SUMPOW = X + Y
      ELSEIF(IPOW.EQ.2) THEN
         SUMPOW = DSQRT(X**2+Y**2)
      ELSEIF(IPOW.EQ.0) THEN
         CALL MWARN('SUMPOW')
         WRITE(*,*)'Error: IPOW=0 not allowed in SUMPOW'
      ELSE
         SUMPOW = (X**IPOW+Y**IPOW)**(1./IPOW)
      ENDIF
      END

      SUBROUTINE MOPERA(I,OPER,J,K,X,Y)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER I,J,K
      CHARACTER OPER*1
      REAL*8 X,Y
C
C     LOCAL
C
      REAL*8 XXX,XSUM,XSUMSQ,XNORM,XAVG,XSQAVG
      INTEGER L
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     START
C
      IF( (BOOK(I)(1:1).NE.'Y'.AND.BOOK(I).NE.'NST') .OR.
     &    (BOOK(J)(1:1).NE.'Y'.AND.BOOK(J).NE.'NST') ) RETURN
      IF(NBIN(I).NE.NBIN(J).AND.(OPER.EQ.'+'.OR.OPER.EQ.'-'.OR.OPER.EQ.
     &    '*'.OR.OPER.EQ.'/'.OR.OPER.EQ.'M'.OR.OPER.EQ.'A')) THEN
	  CALL MWARN('MOPERA')
          WRITE(*,*) I,J
  20      FORMAT(' ****** INCOMPATIBLE OPERATION HIST ',I2,' &',I2,
     &    '*******'/)
          RETURN
      ENDIF
      IF(BOOK(K)(1:1).NE.'Y')
     &  CALL MBOOK(K,TITLE(I),HDEL(I),HMIN(I),HMAX(I))
      IF(OPER.EQ.'E') THEN
c If I contains the accumulated weights, J the accumulated squares of the
c weights and IHIS(J,1) the number of accumulated entries, 'E' will add
c the average value of I to K and will put in J the quadratic dispersion.
         IF(IHIS(J,1).NE.0) THEN
            XXX = 1./IHIS(J,1)
         ELSE
            XXX = 0
         ENDIF
         DO L=1,NBIN(I)
            XSUM   = HIST(I,L)
            XSUMSQ = HIST(J,L)
            HIST(K,L)=HIST(K,L) + XXX*HIST(I,L)
            HIST(J,L)=XXX*DSQRT(ABS(XSUMSQ-XSUM**2*XXX))
         ENDDO
         IENT(K)=IENT(K)+IENT(I)
         XSUM = USCORE(I)
         XSUMSQ = USCORE(J)
         USCORE(K) = USCORE(K)+XXX*XSUM
         USCORE(J) = XXX*DSQRT(ABS(XSUMSQ-XSUM**2*XXX))
         XSUM = OSCORE(I)
         XSUMSQ = OSCORE(J)
         OSCORE(K) = OSCORE(K)+XXX*XSUM
         OSCORE(J) = XXX*DSQRT(ABS(XSUMSQ-XSUM**2*XXX))
      ELSEIF(OPER.EQ.'Q') THEN
         DO L=1,NBIN(I)
            HIST(K,L) = SQRT(HIST(J,L)**2+HIST(I,L)**2)
         ENDDO
         USCORE(K) = SQRT(USCORE(J)**2+USCORE(I)**2)
         OSCORE(K) = SQRT(OSCORE(J)**2+OSCORE(I)**2)
      ELSEIF(OPER.EQ.'A') THEN
         DO L=1,NBIN(I)
            HIST(J,L) = HIST(J,L) + HIST(I,L)
            IHIS(J,L) = IHIS(J,L) + IHIS(I,L)
            HIST(K,L) = HIST(K,L) + HIST(I,L)**2
            IHIS(K,L) = IHIS(K,L) + 1
            HIST(I,L) = 0
            IHIS(I,L) = 0
         ENDDO
         IENT(J) = IENT(J)+IENT(I)
         IUSCORE(J) = IUSCORE(J) + IUSCORE(I)
         USCORE(J) = USCORE(J) + USCORE(I)
         IOSCORE(J) = IOSCORE(J) + IOSCORE(I)
         OSCORE(J) = OSCORE(J) + OSCORE(I)
         IENT(K) = IENT(K)+1
         IUSCORE(K) = IUSCORE(K) + 1
         USCORE(K) = USCORE(K) + USCORE(I)**2
         IOSCORE(K) = IOSCORE(K) + 1
         OSCORE(K) = OSCORE(K) + OSCORE(I)**2
         IENT(I) = 0
         IUSCORE(I) = 0
         IOSCORE(I) = 0
         USCORE(I) = 0
         OSCORE(I) = 0
      ELSEIF(OPER.EQ.'X') THEN
         DO L=1,NBIN(I)
            HIST(I,L) = 0
            IHIS(I,L) = 0
         ENDDO
         IENT(I) = 0
         IUSCORE(I) = 0
         IOSCORE(I) = 0
         USCORE(I) = 0
         OSCORE(I) = 0
      ELSE
        DO L=1,NBIN(I)
      	IF(OPER.EQ.'+') THEN
       	  HIST(K,L)=X*HIST(I,L) + Y*HIST(J,L)
      	ELSEIF(OPER.EQ.'-') THEN
      	  HIST(K,L)=X*HIST(I,L) - Y*HIST(J,L)
      	ELSEIF(OPER.EQ.'*') THEN
      	  HIST(K,L)=X*HIST(I,L) * Y*HIST(J,L)
      	ELSEIF(OPER.EQ.'/') THEN
          IF(Y.EQ.0..OR.HIST(J,L).EQ.0.) THEN
            HIST(K,L)=0.
          ELSE
            HIST(K,L)=X*HIST(I,L) / (Y*HIST(J,L))
          ENDIF
       	ELSEIF(OPER.EQ.'F') THEN
      	  HIST(K,L)=X*HIST(I,L)
      	ELSEIF(OPER.EQ.'R') THEN
          IF(HIST(I,L).GT.0.) THEN
            HIST(K,L)=X*DSQRT(HIST(I,L))
          ELSE
            HIST(K,L)=0.
          ENDIF
      	ELSEIF(OPER.EQ.'S') THEN
          HIST(K,L)=X*HIST(I,L)**2
      	ELSEIF(OPER.EQ.'L') THEN
          IF(HIST(I,L).EQ.0..OR.J.EQ.0.) THEN
             HIST(K,L)=0.
           ELSE
             HIST(K,L)=X*LOG10(Y*HIST(I,L))
           ENDIF
      	ELSEIF(OPER.EQ.'M') THEN
           IF(I.NE.J) XNORM=HIST(I,L)
           IF(I.EQ.J) XNORM=DFLOAT(IHIS(J,L))
           IF(XNORM.NE.0.) THEN
             XAVG=HIST(J,L)/XNORM
             HIST(K,L)=
     &       DSQRT(ABS(-XAVG**2+HIST(K,L)/XNORM)/DFLOAT(IHIS(I,L)))
             HIST(J,L)=XAVG
           ELSE
             HIST(K,L)=0.
             HIST(J,L)=0.
           ENDIF
      	ELSEIF(OPER.EQ.'V') THEN
           XAVG=HIST(I,L)*X
           XSQAVG=HIST(J,L)*X
           XNORM=DFLOAT(IHIS(I,L))*X
           IF(XNORM.NE.0.) THEN
              HIST(K,L)=DSQRT(ABS(XSQAVG-XAVG**2)/XNORM)
              HIST(I,L)=XAVG
           ELSE
              HIST(K,L)=0.
           ENDIF
      	ELSE
	 CALL MWARN('MOPERA')
         WRITE(*,*) OPER
   5     FORMAT(' ****** OPERATION ="',A1,'" UNKNOWN ********'/)
         RETURN
        ENDIF
        ENDDO
      	IF(OPER.EQ.'+') THEN
       	  USCORE(K)=X*USCORE(I) + Y*USCORE(J)
       	  OSCORE(K)=X*OSCORE(I) + Y*OSCORE(J)
      	ELSEIF(OPER.EQ.'-') THEN
      	  USCORE(K)=X*USCORE(I) - Y*USCORE(J)
      	  OSCORE(K)=X*OSCORE(I) - Y*OSCORE(J)
      	ELSEIF(OPER.EQ.'*') THEN
      	  USCORE(K)=X*USCORE(I) * Y*USCORE(J)
      	  OSCORE(K)=X*OSCORE(I) * Y*OSCORE(J)
      	ELSEIF(OPER.EQ.'/') THEN
          IF(Y.EQ.0..OR.USCORE(J).EQ.0.) THEN
            USCORE(K)=0.
          ELSE
            USCORE(K)=X*USCORE(I) / (Y*USCORE(J))
          ENDIF
          IF(Y.EQ.0..OR.OSCORE(J).EQ.0.) THEN
            OSCORE(K)=0.
          ELSE
            OSCORE(K)=X*OSCORE(I) / (Y*OSCORE(J))
          ENDIF
       	ELSEIF(OPER.EQ.'F') THEN
      	  USCORE(K)=X*USCORE(I)
      	  OSCORE(K)=X*OSCORE(I)
      	ELSEIF(OPER.EQ.'R') THEN
          IF(USCORE(I).GT.0.) THEN
            USCORE(K)=X*DSQRT(USCORE(I))
          ELSE
            USCORE(K)=0.
          ENDIF
          IF(OSCORE(I).GT.0.) THEN
            OSCORE(K)=X*DSQRT(OSCORE(I))
          ELSE
            OSCORE(K)=0.
          ENDIF
      	ELSEIF(OPER.EQ.'S') THEN
          USCORE(K)=X*USCORE(I)**2
          OSCORE(K)=X*OSCORE(I)**2
      	ELSEIF(OPER.EQ.'L') THEN
          IF(USCORE(I).EQ.0..OR.J.EQ.0.) THEN
             USCORE(K)=0.
           ELSE
             USCORE(K)=X*LOG10(Y*USCORE(I))
           ENDIF
          IF(OSCORE(I).EQ.0..OR.J.EQ.0.) THEN
             OSCORE(K)=0.
           ELSE
             OSCORE(K)=X*LOG10(Y*OSCORE(I))
           ENDIF
        ENDIF
      ENDIF
      RETURN
      END

      SUBROUTINE MZERO(N)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
C
C     LOCAL
C
      INTEGER I
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     START
C
      BOOK(N)='RES'
      IENT(N)=0
      IUSCORE(N)=0
      IOSCORE(N)=0
      HAVG(N)=0.
      HINT(N)=0.
      DO 1 I=1,NBIN(N)
   1  HIST(N,I)=0.
      END

      SUBROUTINE MRESET(N)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     START
C
      BOOK(N)='RES'
      END

      SUBROUTINE PUTTAG(J,NAME)
C*************************************************************
c Per marcare un istogramma
      IMPLICIT NONE
C
C     ARGUMENTS
C
      CHARACTER*(*) NAME
      INTEGER J
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     LOCAL
C
      CHARACTER*(*) TAG
C
C     START
C
      BOOK(J) = NAME
      RETURN
      ENTRY GETTAG(J,TAG)
      TAG = BOOK(J)
      END

      SUBROUTINE MFINAL(N)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
C
C     LOCAL
C
      INTEGER J,IF
      REAL*8 AVG,XIN,SIG,X
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     START
C
      IF(BOOK(N)(1:1).NE.'Y') RETURN
      AVG=0
      XIN=0
      SIG=0
      IF=0
      DO J=1,NBIN(N)
         X=HIST(N,J)
 	 AVG=AVG+X*XHIS(N,J)
         XIN=XIN+X
	 IF(X.NE.0) IF=1
      ENDDO
      IF(XIN.EQ.0) GO TO 10
      AVG = AVG/XIN
      DO J=1,NBIN(N)
         SIG=HIST(N,J)*(XHIS(N,J)-AVG)**2+SIG
      ENDDO
      SIG=DSQRT(ABS(SIG/XIN))
 10   CONTINUE
      HINT(N) = XIN
      HAVG(N) = AVG
      HSIG(N) = SIG
      END

      SUBROUTINE MNORM(N,X)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
      REAL*8  X
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     LOCAL
C
      INTEGER I
      REAL*8 Y
C
C     START
C
      IF(BOOK(N)(:1).NE.'Y')RETURN
      IF(HINT(N).EQ.0.) THEN
	CALL MWARN('MNORM')
	WRITE(*,*)' INTEGRAL HIST ',N,' IS ZERO: CANNOT RENORMALIZE'
	RETURN
      ELSE
	Y=X/HINT(N)
      ENDIF
      DO 1, I=1,NBIN(N)
    1 HIST(N,I)=HIST(N,I)*Y
      HINT(N)=X
      OSCORE(N)=OSCORE(N)*Y
      USCORE(N)=USCORE(N)*Y
      END

      SUBROUTINE MPRINT(N)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N
C
C     LOCAL
C
      INTEGER INI,J
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     START
C
      DATA INI/0/
      IF(INI.EQ.0) THEN
c     CALL IDATE(IMON,IDAY,IYEAR)
c     CALL TIME(CTIME)
      INI=1
      ENDIF
      IF(BOOK(N)(:1).NE.'Y') RETURN
C      WRITE(98,7) N,IYEAR,IMON,IDAY,CTIME(1:5)
      WRITE(98,*) TITLE(N)
      DO 1 J=1,NBIN(N)
      IF(HIST(N,J).EQ.0.) GO TO 1
      WRITE(98,'(3X,F10.4,2X,E15.4)')
     &                            XHIS(N,J),HIST(N,J)
    1 CONTINUE
      WRITE(98,15) HAVG(N),HSIG(N),HINT(N)
      WRITE(98,20) IENT(N),IUSCORE(N),IOSCORE(N)
C    7 FORMAT(4X,'HIST = ',I3,'   19',I2,'-',I2,'-',I2,1X,A5/)
   10 FORMAT(4X,2E10.3)
   15 FORMAT(/' AVG =',E10.3,4X,' RMS =',E10.3,' INTEGRAL =',E10.3,/)
   20 FORMAT(' ENTRIES=',I5,4X,'UNDERSCORE=',I5,4x,'OVERSCORE=',I5,//)
      END

      SUBROUTINE MTOP4(N,M,BTIT,LTIT,SCALE)

C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER N,M
      CHARACTER*(*) LTIT,BTIT
      CHARACTER*(*) SCALE
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     LOCAL
C
      INTEGER INI,J
      real*8 maxbin
C
C     START
C
      DATA INI/0/
      IF(INI.EQ.0) THEN
c      CALL IDATE(IMON,IDAY,IYEAR)
c      CALL TIME(CTIME)
      INI=1
      ENDIF
      IF(BOOK(N)(:1).NE.'Y') RETURN
cRF
      IF (N.eq.1)WRITE(99,'(A)')' SET DEVICE POSTSCRIPT ORIENT=3'

      WRITE(99,100) TITLE(N),BTIT,LTIT,SCALE,HMIN(N),HMAX(N)
  100 FORMAT( /1x,
     &' SET INTENSITY 4'/,1X,
     &' SET WINDOW Y 2.5 TO 9.'/,1X,
     &' SET WINDOW X 2.5 TO 12.'/,1X,
     &' SET FONT DUPLEX '/1X,
     &' SET SYMBOL 5O SIZE 1.8'/,1X,
     &' TITLE TOP ','"',A20,'"',/1X,
     &' TITLE BOTTOM ','"',A20,'"',/1X,
     &' TITLE LEFT ','"',A20,'"',/1X,
     &' SET SCALE Y ',A,/1X,
     &' (SET TICKS TOP OFF)   '/1x,
     &' SET LIMITS X ',F10.5,' ',F10.5,/1X,
     &' SET ORDER X Y DY ')
      maxbin=0d0
      DO J=1,NBIN(N)
         maxbin=max(maxbin,abs(HIST(N,J)))
      ENDDO
      DO 1 J=1,NBIN(N)
         WRITE(99,'(3X,F10.4,2(2X,E15.4))')
     &        XHIS(N,J),HIST(N,J),HIST(M,J)
    1 CONTINUE
      WRITE(99,200)
  200 FORMAT('   HISTO')
  300 FORMAT( /1x,
     &' (INFO-BOX'/,1X,
     &' BOX 7.3 1.3 SIZE 9.5 0.7'/,1X,
     &' SET WINDOW Y 0. TO 2.'/,1X,
     &' SET TITLE SIZE -1.5'/1X,
     &' SET FONT DUPLEX '/1X,
     &' TITLE 3.8 1.4 "X-sect =',1PE10.3,'(pb)   AVG =',1PE10.3,
     &             '   RMS =',1PE10.3,'"',/1X,
     &' TITLE 3.8 1.1 "Entries =',I8,2x,'Entries =',I8,2X
     &                ,'Undersc =',I6,2X,'Oversc =',I6,'"',/1X,
     &' SET TITLE SIZE -2')
      WRITE(99,400)
  400 FORMAT('   NEW PLOT')
      END

      SUBROUTINE MULTITOP4(NH,NE,N,M,BTIT,LTIT,SCA)
C*************************************************************
      IMPLICIT NONE
C
C     ARGUMENTS
C
      INTEGER NH,NE,N,M
      CHARACTER*(*) LTIT,BTIT
      CHARACTER*(*) SCA
C
C     GLOBAL
C
      INCLUDE 'dbook.inc'
C
C     LOCAL
C
      REAL*8 YTIT,XTIT,FEXP,FMAX,FMIN,X,XMX,FMX,FMN
      REAL*8 YU,XU,YL,XL,YTIT0,YD,XTIT0,TITX,TITS,TICS,SRED,XD
      INTEGER IBIN,IPNS,J,I,NOLD,IFRAME,IFRMAX,IP,NS,MOLD
      INTEGER INI
C
C     START
C
      CHARACTER SCALE*3
      CHARACTER*7  PLOT(4)
      DATA PLOT/'SOLID','DASHES','DOTS','DOTDASH'/
C  PLOT SIZE, CORNERS
      REAL*8 WIDTH,HEIGHT,XCORN,YCORN
      DATA WIDTH,HEIGHT/11.5,8.5/,XCORN,YCORN/1.5,1./
C  PLOT VERSUS TEXT FRACTION
      REAL*8 XPFRAC,YPFRAC,XTFRAC,YTFRAC
      DATA XPFRAC,YPFRAC/0.75,0.75/,XTFRAC,YTFRAC/0.25,0.25/
C  DEFAULT SIZES
      REAL*8 TIC0,LAB0,TIT0,LABS
      DATA TIT0,LAB0,TIC0/3d0,3d0,0.06d0/
      DATA INI/0/
      IF(INI.EQ.0) THEN
c      CALL IDATE(IMON,IDAY,IYEAR)
c      CALL TIME(CTIME)
      IFRAME=0
C      WRITE(99,71) IYEAR,IMON,IDAY,CTIME(1:5)
C   71 FORMAT(4X,' (   19',I2,' -',I2,' -',I2,1X,A5/)
      INI=1
      ENDIF
      IF(SCA.EQ.'REF') THEN
	IFRAME=0
	RETURN
      ENDIF
      IF(BOOK(NH)(:1).NE.'Y') RETURN
      IFRMAX=N*M
      IFRAME=IFRAME+1
      IF(IFRAME.GT.IFRMAX.OR.N.NE.NOLD.OR.M.NE.MOLD) THEN
      	IFRAME=1
        WRITE(99,202)
C        WRITE(99,1) IMON,IDAY,CTIME(1:5)
C  1     FORMAT(' SET FONT DUPLEX',/,'  SET TITLE SIZE 2',/,
C     +      ' TITLE 12.8 9 ANGLE -90 ','" MLM   ',I2,'-',I2,1X,A5,'"')
      ENDIF
      IF(IFRAME.EQ.1) THEN
    	I=1
	J=1
      ELSEIF(IFRAME.LE.IFRMAX) THEN
	IF(I.LE.N) I=I+1
        IF(I.GT.N) THEN
		I=1
		J=J+1
	ENDIF
      ENDIF
      IF(N.EQ.NOLD) GO TO 10
      NS=N-1
      XD=WIDTH/DFLOAT(N)
      SRED=DSQRT(DFLOAT(N*M))
      TITS=TIT0/SRED
      LABS=LAB0/SRED
      TICS=TIC0/SRED
      XTIT0=0.55*XPFRAC*XD
      NOLD=N
10    IF(M.EQ.MOLD) GO TO 20
      YD=HEIGHT/DFLOAT(M)
      YTIT0=0.06*YD
      MOLD=M
20    CONTINUE
      XL=(I-1)*XD + XCORN
      YL=(M-J)*YD + YCORN
      XU=XL+XD*XPFRAC
      YU=YL+YD*YPFRAC
      IP=0
      XMX=0.
      DO IBIN=1,NBIN(NH)
        X=HIST(NH,IBIN)
        IF(X.NE.0.) THEN
           IF(XMX.EQ.0.) THEN
              FMX = X + HIST(NE,IBIN)
              FMN = X - HIST(NE,IBIN)
           ELSE
              FMX=MAX(FMX,X + HIST(NE,IBIN))
              FMN=MIN(FMN,X - HIST(NE,IBIN))
           ENDIF
        ENDIF
        XMX=MAX(XMX,ABS(X)+ HIST(NE,IBIN))
      ENDDO
      SCALE=SCA
50    IF(SCALE.EQ.'LIN') THEN
        IF(FMN.GE.0.)   FMIN=0.
        IF(FMN.LT.0.)   FMIN=FMN*1.3
        IF(FMX.GT.0.)   FMAX=FMX*1.3
        IF(FMX.LT.0.)   FMAX=0.
      ELSEIF(SCALE.EQ.'LOG') THEN
cRF
c$$$        IF(FMN.LE.0.) THEN
c$$$                SCALE='LIN'
c$$$                GO TO 50
c$$$        ENDIF
        FMAX=10.**( AINT(LOG10(ABS(FMX))+1000001) - 1000000 )
        FMIN=10.**( AINT(LOG10(ABS(FMN))+1000000) - 1000000 )
      ENDIF
      WRITE(99,100) TITS,LABS,TICS,XL,XU,YL,YU
100   FORMAT(2X,'( SET FONT DUPLEX',/,
     *       2X,'SET TITLE SIZE ',F8.4,/,
     *       2X,'SET LABEL SIZE ',F8.4,/,
     *       2X,'SET TICKS TOP OFF SIZE ',F8.4,/,
     *       2X,'SET WINDOW X ',F8.4,' TO ',F8.4,/,
     *       2X,'SET WINDOW Y ',F8.4,' TO ',F8.4)
      XTIT=XL+XTIT0
      YTIT=YU+YTIT0
      WRITE(99,101) XL,YTIT,TITLE(NH)(1:40)
101   FORMAT('  TITLE ',2(F8.4,1X),'"',A,'"')
      YTIT=YTIT-2.*YTIT0
      WRITE(99,102) XTIT,YTIT,HINT(NH)
102   FORMAT('  TITLE ',2(F8.4,1X),'" INT=',1PE10.3,'"')
      YTIT=YTIT-YTIT0
      WRITE(99,103) XTIT,YTIT,IENT(NH)
103   FORMAT('  TITLE ',2(F8.4,1X),'" ENT=',I9,'"')
      YTIT=YTIT-YTIT0
      IF(USCORE(NH).NE.0.) THEN
        WRITE(99,104) XTIT,YTIT,USCORE(NH)
104     FORMAT('  TITLE ',2(F8.4,1X),'" UFL=',1PE10.3,'"')
        YTIT=YTIT-YTIT0
      ENDIF
      IF(OSCORE(NH).NE.0.) THEN
        WRITE(99,105) XTIT,YTIT,OSCORE(NH)
105     FORMAT('  TITLE ',2(F8.4,1X),'" OFL=',1PE10.3,'"')
        YTIT=YTIT-YTIT0
      ENDIF
      WRITE(99,106) XTIT,YTIT,XU,YTIT,XTIT,YTIT,XTIT,YU
106   FORMAT(2X,'SET ORD X Y ',/,2(F8.4,1X),/,2(F8.4,1X),/,
     *       2X,'JOIN TEXT',/,
     *       2X,2(F8.4,1X),/,2(F8.4,1X),/,
     *       2X,'JOIN TEXT')
      WRITE(99,108) TITS*1.5
108   FORMAT(2X,'SET TITLE SIZE ',F8.4)
      WRITE(99,107) BTIT,XL-0.75*XD*XTFRAC,YL+(YU-YL)/3.,LTIT,SCALE,
     * HMIN(NH),HMAX(NH),FMIN,FMAX
107   FORMAT(
     &' TITLE BOTTOM ','"',A,'"',/1X,
     &' TITLE ',f10.5,f10.5,' ANGLE 90 ','"',A,'"',/1X,
     &' SET SCALE Y ',A,/1X,
     &' SET TICKS TOP OFF   '/1x,
     &' SET LIMITS X ',F10.5,' ',F10.5,/1X,
     &' SET LIMITS Y ',1PE10.3,' ',1PE10.3,/1X,
     &' SET ORDER X Y DY')
C
C  END HEADER , FILL TOPDRAWER WITH DATA
C
      ENTRY MTFILL(NH,NE,N,M,BTIT,LTIT,SCA)
      IP=IP+1
      IF(IP.GT.4) IP=1
      WRITE(99,110) TITLE(NH),HINT(NH),IENT(NH)
110   FORMAT(' ( ',A,/,' ( INT=',1PE10.3,'  ENTRIES=',I12)
      DO 200 IBIN=1,NBIN(NH)
      WRITE(99,'(3X,F10.4,2(2X,E15.4))')
     &          XHIS(NH,IBIN),HIST(NH,IBIN),HIST(NE,IBIN)
200   CONTINUE
      WRITE(99,201)  PLOT(IP)
      IF(BOOK(NE).NE.'NO')   WRITE(99,*)  '  PLOT'
201   FORMAT(2X,'HIST ',A)
202   FORMAT('   NEW PLOT',/,/)
203   RETURN
      END

      SUBROUTINE MWARN(ROUT)
C*************************************************************
      CHARACTER*(*) ROUT
      WRITE(*,*) '***********************************************'
      WRITE(*,*) '***** WARNING CALLED FROM ROUTINE ',ROUT,':'
      END

C*******************************************************************
C     END OF THE HISTOGRAMMING PACKAGE
C*******************************************************************

      subroutine mfill(n,var,www)
      implicit none
      integer n,n_by4
      double precision var,www
c     APPLgrid commons
      include "reweight_appl.inc"
      include "appl_common.inc"
      integer iappl
      common /for_applgrid/ iappl
      integer j
      if(iappl.ne.0)then
         do j=1,nh_obs
            if(n.eq.ih_obs(j))then
               appl_obs_num   = j
               appl_obs_histo = var
c     Fill the reference APPLgrid histograms
               call APPL_fill_ref
c     Fill the APPLgrid files
               call APPL_fill
            endif
         enddo
      endif
      n_by4=4*(n-1)+1
      call mfill4(n_by4,var,www)
      return
      end

      subroutine bookup(n,string,del,xl,xu)
      implicit none
      integer i,n,n_by4
      character*(*) string
      double precision del,xl,xu
c     APPLgrid commons
      include "reweight_appl.inc"
      include "appl_common.inc"
      integer iappl
      common /for_applgrid/ iappl
c     Initialize the grids only if the switch "iappl" is different from zero
c     and if the title string containes the word "central" and does not contain
c     the word "Born".
      if(iappl.ne.0.and.index(string,"central").ne.0.and.
     1                  index(string,"Born").eq.0)then
c     Observable parameters
c     Compute number of bins and edges only if they have not been given by the user.
         if(appl_obs_nbins.eq.0)then
            appl_obs_nbins = int( ( xu - xl ) / del / 0.999999d0 )
c     compute bin edges
            do i=0,appl_obs_nbins
               appl_obs_bins(i) = xl + i * del
            enddo
         endif
         appl_obs_min = appl_obs_bins(0)
         appl_obs_max = appl_obs_bins(appl_obs_nbins)
         if(abs(appl_obs_max-xu).gt.0.00000001d0)then
            write(*,*) 'APPLgrid Histogram: ',
     1                 'Change of the upper limit:',xu,'-->',
     2                  appl_obs_max
         endif
c     Initialize APPLgrid routines
         call APPL_init
c     Keep track of the position of this histogram
         nh_obs = nh_obs + 1
         ih_obs(nh_obs) = n
c     Reset number of bins to zero
         appl_obs_nbins = 0
      endif
      n_by4=4*(n-1)+1
      call bookup4(n_by4,string,del,xl,xu)
      return
      end

      subroutine MTOP(n,string1,string2,string3)
      implicit none
      integer n,n_by4,m_by4
      character*(*) string1,string2,string3
c     APPLgrid commons
      include "reweight_appl.inc"
      include "appl_common.inc"
      integer iappl
      common /for_applgrid/ iappl
      integer j
      if(iappl.ne.0)then
         do j=1,nh_obs
            if(n.eq.ih_obs(j))then
               appl_obs_num = j
               call APPL_fill_ref_out
               call APPL_term
            endif
         enddo
      endif
      n_by4=4*(n-1)+1
      m_by4=n_by4+3
c write the 'n' plots with the 'n+3' error bars
      call mtop4(n_by4,m_by4,string1,string2,string3)
      return
      end

      subroutine multitop(n,lr,lh,string1,string2,string3)
      implicit none
      integer n,n_by4,lr,lh,m_by4
      character*(*) string1,string2,string3
c     APPLgrid commons
      include "reweight_appl.inc"
      include "appl_common.inc"
      integer iappl
      common /for_applgrid/ iappl
      integer j
      if(iappl.ne.0)then
         do j=1,nh_obs
            if(n.eq.ih_obs(j))then
               appl_obs_num = j
               call APPL_fill_ref_out
               call APPL_term
            endif
         enddo
      endif
      n_by4=4*(n-1)+1
      m_by4=n_by4+3
c write the 'n' plots with the 'n+3' error bars
      call multitop4(n_by4,m_by4,lr,lh,string1,string2,string3)
      return
      end

      subroutine bookup4(n,string,del,xl,xu)
      implicit none
      integer n
      character*(*) string
      double precision del,xl,xu
c
c Per ogni istogramma da fare, ne sono richiesti quattro
c In n si accumulano i valori in outfun.
c A ogni iterazione di vegas l'istogramma n viene sommato a n+1,
c e il quadrato del suo valore viene sommato a n+2,
c L'istogramma n viene anche usato alla fine per combinare
c i totali dei vari contributi sig0,sig2, etc. mentre
c l'istogramma n+3 viene usato per combinare gli
c errori dei vari contributi sig0,sig2, etc.
c Non si vuole che n e n+3 vengano salvati o riesumati.
c Cambiando il tag in N,N+3, questo non avviene (si guardi
c in mbook e save/restart e anche mclear in questo programma.
c
      call mbook(n,  string,del,xl,xu)
      call mbook(n+1,'tmp ',del,xl,xu)
      call mbook(n+2,'tmp square',del,xl,xu)
      call mbook(n+3,'error ',del,xl,xu)
      call puttag(n,'YST')
      call puttag(n+3,'NST')
      return
      end

      subroutine accum(inclde)
      implicit real * 8 (a-h,o-z)
      include 'dbook.inc'
      PARAMETER (NMB=NPLOTS)
      character * 3 tag
      logical inclde
c
c     Accumula i valori e i valori al quadrato per l'analisi statistica,
c     e svuota l'istogramma di accumulo.
c
      do j=1,nmb-3
         call gettag(j,tag)
         if(tag.eq.'YST') then
             if (inclde) then
c Sum the results in histos 'j' onto 'j+1' and the squares to 'j+2'
c This also empties histograms 'j'
                call mopera(j,'A',j+1,j+2,dum,dum)
             else
c Only empty the histograms 'j'.
                call mopera(j,'X',j,j,dum,dum)
             endif
         endif
      enddo
      return
      end

      subroutine mclear
      implicit real * 8 (a-h,o-z)
      character * 70 files(20), filn
      data nfil/0/
c
c sum up files
      call sumfil(files,nfil)
      nfil=0
      return
      entry addfil(filn)
c adds file filn to the list of save files.
      nfil = nfil+1
      if(nfil.gt.20) then
         write(*,*) 'mclear: too many files.'
         stop
      endif
      files(nfil) = filn
      end

      subroutine sumfil(files,nfil)
      implicit real * 8 (a-h,o-z)
      include 'dbook.inc'
      PARAMETER (NMB=NPLOTS)
      external dummyinit
      character * 3 tag
      character * 70 files(*)
c
      do j=1,nfil
         call resume(files(j),dummyinit,'NO')
c
c completa l'analisi statistica
c
         do k=1,nmb-3
            call gettag(k,tag)
            if(tag.eq.'YST') then
               call addup(k)
            endif
         enddo
      enddo
      end

      subroutine dummyinit()
      return
      end

      subroutine addup(j)
      implicit none
      integer j
      real*8 dum
c
c accumula j+1 riscalato in j e pone in j+2 la stima dell'errore
c
      call mopera(j+1,'E',j+2,j,dum,dum)
c
c accumula l'errore in quadratura
c
      call mopera(j+2,'Q',j+3,j+3,dum,dum)
      end

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


module fields_arrays
  use constants, only: run_name_size
  implicit none

  private

  public :: phi, apar, bpar, phinew, aparnew, bparnew
  public :: apar_ext, time_field, response_file

  !Main fields
  complex, dimension (:,:,:), allocatable :: phi,    apar,    bpar
  complex, dimension (:,:,:), allocatable :: phinew, aparnew, bparnew
  !For antenna
  complex, dimension (:,:,:), allocatable :: apar_ext
  ! (-ntgrid:ntgrid,ntheta0,naky) replicated

  !Timing data
  real :: time_field(2)=0.

  !For response data
  character(run_name_size) :: response_file
end module fields_arrays


text:
	write (*, 10000)
10000	format (' Entered dummy main ')
	end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
        program reverse

        real*8 x(500), y(500)

        open(10,file='lev.bot')

        i = 0

 10     continue
          i = i + 1
          read(10,*,end=20) x(i), y(i)
          goto 10
 20     continue
        close(10)

        num = i - 1

        open(10,file='new.bot')

        do i = 1, num
          write(10,*) x(num-i+1), y(num-i+1)
        end do
        close(10)

        stop
        end


text:
MODULE AccretionShockUtilitiesModule

  USE KindModule, ONLY: &
    DP, &
    Zero, &
    Half, &
    One, &
    Two, &
    Three, &
    Five, &
    TwoPi
  USE ProgramHeaderModule, ONLY: &
    nDOFX, &
    nNodes, &
    nNodesX, &
    nDimsX
  USE MeshModule, ONLY: &
    MeshX, &
    NodeCoordinate
  USE GeometryFieldsModule, ONLY: &
    iGF_Gm_dd_11, &
    iGF_Gm_dd_22, &
    iGF_Gm_dd_33
  USE FluidFieldsModule, ONLY: &
    iCF_D, &
    iCF_S1, &
    iCF_S2, &
    iCF_S3, &
    iCF_E, &
    iCF_Ne, &
    iPF_D, &
    iAF_P
  USE EquationOfStateModule_IDEAL, ONLY: &
    Gamma_IDEAL
  USE EquationOfStateModule, ONLY: &
    ComputePressureFromPrimitive
  USE Euler_BoundaryConditionsModule, ONLY: &
    ApplyBoundaryConditions_Euler
  USE Euler_UtilitiesModule_Relativistic, ONLY: &
    ComputePrimitive_Euler_Relativistic
  USE UnitsModule, ONLY: &
    Erg, &
    Gram, &
    Centimeter
  USE QuadratureModule, ONLY: &
    GetQuadrature
  USE Euler_ErrorModule, ONLY: &
    DescribeError_Euler

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: WriteInitialConditionsToFile
  PUBLIC :: ComputeAccretionShockDiagnostics

CONTAINS

  SUBROUTINE WriteInitialConditionsToFile &
    ( iX_B0, iX_E0, iX_B1, iX_E1, uGF, uCF, FileName )

    INTEGER ,           INTENT(in)    :: &
      iX_B0(3), iX_E0(3), iX_B1(3), iX_E1(3)
    REAL(DP),           INTENT(in)    :: &
      uGF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    REAL(DP),           INTENT(inout) :: &
      uCF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    CHARACTER(LEN=128),   INTENT(in)  :: &
      FileName

    INTEGER  :: iNX, iX1
    REAL(DP) :: D(nNodesX(1)), V1(nNodesX(1)), V2(nNodesX(1)), V3(nNodesX(1)), &
                E(nNodesX(1)), Ne(nNodesX(1)), P (nNodesX(1))
    INTEGER :: iErr(nNodesX(1),iX_B1(1):iX_E1(1))
    CHARACTER(LEN=16) :: FMT

    IF( nDimsX .GT. 1 )THEN

      WRITE(*,*) 'WriteInitialConditionsToFile not implemented for nDimsX > 1.'
      RETURN

    END IF

    WRITE(FMT,'(A3,I3.3,A10)') '(SP', nDOFX, 'ES25.16E3)'

    OPEN( 101, FILE = TRIM( FileName ) // '_D.dat' )
    OPEN( 102, FILE = TRIM( FileName ) // '_V.dat' )
    OPEN( 103, FILE = TRIM( FileName ) // '_P.dat' )

    WRITE(101,'(A16)') TRIM( FMT )
    WRITE(102,'(A16)') TRIM( FMT )
    WRITE(103,'(A16)') TRIM( FMT )

    CALL ApplyBoundaryConditions_Euler( iX_B0, iX_E0, iX_B1, iX_E1, uCF )

    DO iX1 = iX_B1(1), iX_E1(1)

      iErr(:,iX1) = 0

      CALL ComputePrimitive_Euler_Relativistic &
             ( uCF(:,iX1,1,1,iCF_D ), &
               uCF(:,iX1,1,1,iCF_S1), &
               uCF(:,iX1,1,1,iCF_S2), &
               uCF(:,iX1,1,1,iCF_S3), &
               uCF(:,iX1,1,1,iCF_E ), &
               uCF(:,iX1,1,1,iCF_Ne), &
               D, V1, V2, V3, E, Ne, &
               uGF(:,iX1,1,1,iGF_Gm_dd_11), &
               uGF(:,iX1,1,1,iGF_Gm_dd_22), &
               uGF(:,iX1,1,1,iGF_Gm_dd_33), &
               iErr(:,iX1) )

      CALL ComputePressureFromPrimitive &
             ( D, E, Ne, P )

      WRITE(101,FMT) D
      WRITE(102,FMT) V1
      WRITE(103,FMT) P

    END DO

    CLOSE( 103 )
    CLOSE( 102 )
    CLOSE( 101 )

    IF( ANY( iErr .NE. 0 ) )THEN

      DO iX1 = iX_B1(1), iX_E1(1)
      DO iNX = 1, nNodesX(1)

        CALL DescribeError_Euler( iErr(iNX,iX1) )

      END DO
      END DO

    END IF

  END SUBROUTINE WriteInitialConditionsToFile

  SUBROUTINE ComputeAccretionShockDiagnostics &
    ( iX_B0, iX_E0, iX_B1, iX_E1, uPF, uAF, Power_Legendre )

    INTEGER , INTENT(in)    :: iX_B0(3), iX_E0(3), iX_B1(3), iX_E1(3)
    REAL(DP), INTENT(in)    :: uPF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    REAL(DP), INTENT(in)    :: uAF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    REAL(DP), INTENT(inout) :: Power_Legendre(0:)

    CALL ComputePowerInLegendreModes &
           ( iX_B0, iX_E0, iX_B1, iX_E1, uPF, uAF, Power_Legendre )

  END SUBROUTINE ComputeAccretionShockDiagnostics

  SUBROUTINE ComputePowerInLegendreModes &
    ( iX_B0, iX_E0, iX_B1, iX_E1, uPF, uAF, Power )

    ! --- Compute power in Legendre modes a la
    !     Blondin & Mezzacappa (2006), ApJ, 642, 401 ---

    INTEGER , INTENT(in)    :: iX_B0(3), iX_E0(3), iX_B1(3), iX_E1(3)
    REAL(DP), INTENT(in)    :: uPF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    REAL(DP), INTENT(in)    :: uAF(1:,iX_B1(1):,iX_B1(2):,iX_B1(3):,1:)
    REAL(DP), INTENT(inout) :: Power(0:)

    INTEGER  :: iX1, iX2, iX3, iNX, iNX1, iNX2, nX(3), nDOF(3), &
                NodesX2(nNodes), iX, iNX2_L, iNX2_U, iNX1_L, iNX1_U
    REAL(DP) :: dX1, dX2
    REAL(DP) :: xQ(nNodes), wQ(nNodes)
    REAL(DP) :: G(0:2,nNodes,iX_B0(1):iX_E0(1))
    REAL(DP) :: Field(nDOFX,iX_B0(1):iX_E0(1), &
                            iX_B0(2):iX_E0(2), &
                            iX_B0(3):iX_E0(3))
    REAL(DP) :: FieldUnits
    REAL(DP), ALLOCATABLE :: X1(:), X2(:), P0(:), P1(:), P2(:)

    nX   = iX_E0 - iX_B0 + 1
    nDOF = nNodesX * nX

    FieldUnits = ( Erg / Centimeter**3 ) &
                   / ( Gram / Centimeter**3)**( Gamma_IDEAL )

    CALL GetQuadrature( nNodes, xQ, wQ )

    ALLOCATE( X1(iX_B0(1):iX_B0(1)+nDOF(1)-1) )
    ALLOCATE( X2(iX_B0(2):iX_B0(2)+nDOF(2)-1) )
    ALLOCATE( P0(iX_B0(2):iX_B0(2)+nDOF(2)-1) )
    ALLOCATE( P1(iX_B0(2):iX_B0(2)+nDOF(2)-1) )
    ALLOCATE( P2(iX_B0(2):iX_B0(2)+nDOF(2)-1) )

    ! --- Populate X1 array ---

    iX = 0

    DO iX1 = iX_B1(2), iX_E1(2)

      DO iNX1 = 1, nNodesX(1)

        iX = iX + 1

        X1(iX) = NodeCoordinate( MeshX(1), iX1, iNX1 ) / Centimeter

      END DO

    END DO

    ! --- Populate X2 array ---

    iX = 0

    DO iX2 = iX_B0(2), iX_E0(2)

      DO iNX2 = 1, nNodesX(2)

        iX = iX + 1

        X2(iX) = NodeCoordinate( MeshX(2), iX2, iNX2 )

      END DO

    END DO

    CALL ComputeLegendrePolynomials( X2, P0, P1, P2 )

    ! --- Define field used for computing power (e.g., entropy) ---

    DO iX3 = iX_B0(3), iX_E0(3)
    DO iX2 = iX_B0(2), iX_E0(2)
    DO iX1 = iX_B0(1), iX_E0(1)
    DO iNX = 1, nDOFX

      ! --- Entropy ---

      Field(iNX,iX1,iX2,iX3) &
        = LOG10( ( uAF(iNX,iX1,iX2,iX3,iAF_P) &
                     / uPF(iNX,iX1,iX2,iX3,iPF_D)**( Gamma_IDEAL ) ) &
                     / FieldUnits )

    END DO
    END DO
    END DO
    END DO

    Power = Zero
    G     = Zero

    ! --- Loop over radii ---

    DO iX1 = iX_B0(1), iX_E0(1)

      dX1 = MeshX(1) % Width(iX1) / Centimeter

      DO iNX1 = 1, nNodes

        CALL ComputeIntegrationNodesX2( iNX1, NodesX2 )

        ! --- For each radius, compute moments (G functions, Eq. (7)) ---

        DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)

          dX2 = MeshX(2) % Width(iX2)

          iNX2_L = iX_B0(2) + nNodesX(2) * ( iX2 - iX_B0(2) )
          iNX2_U = iNX2_L + nNodesX(2) - 1

          G(0,iNX1,iX1) &
            = G(0,iNX1,iX1) &
                + dX2 * SUM( wQ * Field(NodesX2,iX1,iX2,iX3) &
                               * P0(iNX2_L:iNX2_U) * SIN( X2(iNX2_L:iNX2_U) ) )

          G(1,iNX1,iX1) &
            = G(1,iNX1,iX1) &
                + dX2 * SUM( wQ * Field(NodesX2,iX1,iX2,iX3) &
                               * P1(iNX2_L:iNX2_U) * SIN( X2(iNX2_L:iNX2_U) ) )

          G(2,iNX1,iX1) &
            = G(2,iNX1,iX1) &
                + dX2 * SUM( wQ * Field(NodesX2,iX1,iX2,iX3) &
                               * P2(iNX2_L:iNX2_U) * SIN( X2(iNX2_L:iNX2_U) ) )

        END DO ! iX2
        END DO ! iX3

      END DO ! iNX1

      ! --- Compute powers (Eq. (8)) ---

      iNX1_L = iX_B0(1) + nNodesX(1) * ( iX1 - iX_B0(1) )
      iNX1_U = iNX1_L + nNodesX(1) - 1

      Power(0) &
        = Power(0) &
            + TwoPi * dX1 * SUM( wQ * G(0,:,iX1)**2 * X1(iNX1_L:iNX1_U)**2 )

      Power(1) &
        = Power(1) &
            + TwoPi * dX1 * SUM( wQ * G(1,:,iX1)**2 * X1(iNX1_L:iNX1_U)**2 )

      Power(2) &
        = Power(2) &
            + TwoPi * dX1 * SUM( wQ * G(2,:,iX1)**2 * X1(iNX1_L:iNX1_U)**2 )

    END DO ! X1

    DEALLOCATE( P2 )
    DEALLOCATE( P1 )
    DEALLOCATE( P0 )
    DEALLOCATE( X2 )
    DEALLOCATE( X1 )

  END SUBROUTINE ComputePowerInLegendreModes

  SUBROUTINE ComputeLegendrePolynomials( X2, P0, P1, P2 )

    REAL(DP), INTENT(in)  :: X2(:)
    REAL(DP), INTENT(out) :: P0(:), P1(:), P2(:)

    REAL(DP) :: X(SIZE(X2))

    X = COS( X2 )

    P0 = SQRT( One   / Two )
    P1 = SQRT( Three / Two ) * X
    P2 = SQRT( Five  / Two ) * Half * ( Three * X**2 - One )

  END SUBROUTINE ComputeLegendrePolynomials

  SUBROUTINE ComputeIntegrationNodesX2( iNX1, NodesX2 )

    INTEGER, INTENT(in)  :: iNX1
    INTEGER, INTENT(out) :: NodesX2(:)

    IF     ( nNodes .EQ. 1 )THEN

      NodesX2 = [ 1 ]

    ELSE IF( nNodes .EQ. 2 )THEN

      IF( iNX1 .EQ. 1 )THEN

        NodesX2 = [ 1, 3 ]

      ELSE

        NodesX2 = [ 2, 4 ]

      END IF

    ELSE IF( nNodes .EQ. 3 )THEN

      IF( iNX1 .EQ. 1 )THEN

        NodesX2 = [ 1, 4, 7  ]

      ELSE IF( iNX1 .EQ. 2 )THEN

        NodesX2 = [ 2, 5, 9 ]

      ELSE

        NodesX2 = [ 3, 6, 9 ]

      END IF

    ELSE

      PRINT*, 'Not implemented for nNodes > 3. Stopping...'
      STOP

    END IF

  END SUBROUTINE ComputeIntegrationNodesX2

END MODULE AccretionShockUtilitiesModule

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


program functionsTest
    use nt_TypesModule
    use nt_NeuronModule
    use sft_AssertModule
    use sft_SuiteModule
    implicit none

    type(sft_Suite) :: testSuite

    call sft_init(testSuite)

    call sft_summary(testSuite)

end program functionsTest