!  Xpect takes expectation values of operators defined in prop.
!  Version considers diagonal case only. Nikolai Zobov (2002)
      program leek
      call xpect3
      stop
      end
!**************************************************001
      subroutine xpect3
!
!     xpect3, and the various subroutines it calls, calculate the
!     expectation values of vib-rot eigenfunctions with respect
!     the properties defined by subroutine prop.
!
      implicit double precision (a-h,o-y), logical (z)
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      namelist/prt/ zprint,ztra,zstart,zfit,zform,&
                     iket,itra,itra0,ilev
      common /timing/ itime0
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      character(len=8) title(9)

      write(6,200)
200   format(//5x,'program xpect3 (version of 19 March 2003):',/)
!
!     read in logical control parameters zprint,ztra,zstart,zfit,zform
!
!     zfit = .true. expectation values for potential fitting
!                  (relative to the ground state)
!  *  zprint= .true. gives extra print out - useful for debugging.
!  *  ztra= .true. outputs the property data to strean itra
!  *  the default values of all three parameters are .false.
!     the variable "zform" controls whether the derivatives are outputted
!     formatted (true) or unformatted (false - default). (jhs)

      read(5,prt)
      read(5,100) title

100   format(9a8)
      write(6,202) title
202   format(/5x,9a8)

!     read control parameters
      call system_clock(itime0,irate2,imax2)
      call insize
!     now do the calculations
      call core
      call system_clock(itime2,irate2,imax2)
      itime=(itime2-itime0)/irate2
      write(6,1) itime
1     format(/i10,'secs cpu time used'/)
      stop
      end

!**************************************************002
      block data
      implicit double precision (a-h,o-y), logical (z)
!
!     block data stores the default values of the logical control
!     parameters and the ground zero energy.
!     it also stores the default values of the ket input stream, iket,
!     and the output stream itra.
!
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      data zncor/.false./, zprint/.false./,zfit/.true./,&
           ztra/.true./, zstart/.false./,&
           iket/11/,itra/12/,itra0/28/,ilev/57/,zform/.true./
      end

!**************************************************003
      subroutine insize
!
!     subroutine insize reads in the parameters which control the
!     overall size of the problem. some of these are passed from the
!     output streams of programs triatom and rotlev via the input
!     stream iket.
!     this means that many parameters are duplicated, giving the program
!     the opportunity to check that the bra and ket are compatible.
!     integration parameters for the two radial co-ordinates and the
!     angular functions must be set by the user, and are inputted on
!     stream 5 (five).
!
      implicit double precision (a-h,o-y),logical (z)
!
!     the following size parameters have these functions:
!
!     mbass: the total number of basis functions for the ket.
!     npnt1: number of gauss-laguerre dvr points for r1
!     npnt2: number of gauss-laguerre dvr points for r2
!     lmax : max number of angular basis functions in ket
!     lpot : max order of legendre polys in property expansion
!     npot : number of gauss-legendre integration points
!     npnt : max(npnt1,npnt2)
!     idia : = 1 scattering coordinates hetronuclear diatomic
!            = 2 scattering coordinates homonuclear diatomic
!            = -1 radau coordinates hetronuclear diatomic
!            = -2 radau coordinates homonuclear diatomic
!     ipar : parity of basis - if |idia|=2: ipar=0 for even & =1 for odd
!          : taken from ipar1 and ipar2
!     kmin : kmin=1 for sym. rotational basis, =0 for anti-sym.
!          : for non-coriolis calculations, kmin= k.
!     neval: number of eigenvalues supplied from rotlev or triatom
!     ncoord: number of vibrational coordinates explicitly considered
!     if (ncoord .ne. 3) some of the above are dummies, see below.
!     nprop: number of properties to be considered
!     iprop: vector with the information on which properties
!            will be considered in the run
!
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      parameter (mxprop=300)
      common /pot/ iprop(mxprop)
!     save masses & g in case they are needed in the properties routine
!
!     read in control parameters of problem:
!     ncoord = 2: atom-diatom problem with diatom rigid
!     ncoord = 3: full 3-d triatomic problem
!
!     read in of data for the wavefunction
!
      if (zfit) write(6,1050)
 1050 format(/5x,'Property analysis for potential energy fitting')
      if (.not. zfit) write(6,1051)
 1051 format(/5x,'Property analysis to obtain expectation values')
      write(6,1057) iket
 1057 format(//5x,'input ket wavefunctions read from stream ',&
             'iket  =',i4)
      if (ztra) then
           write(6,1058) itra
 1058     format(5x,'output transition data written to stream ',&
             'itra  =',i4)
           if (zform) write(6,1055)
 1055     format(5x,'file written formatted')
           if (.not.zform) write(6,1056)
 1056     format(5x,'file written unformatted')
           if (.not. zfit) then
              if (zstart) write(6,1052)
 1052        format(5x,'file written from the start')
              if (.not.zstart) write(6,1053)
 1053        format(5x,'data appended to the end of file')
           endif
      endif

      open(unit=iket,form='unformatted')
      read(iket) idia,ipar,lmax,nmax11,nmax12,jrot,kmin1,neval
      read(iket) zembed,zmors1,zmors2,xmass,g1,g2,zncor
      read(iket) re1,diss1,we1,re2,diss2,we2

      ncoord=3
      if (idia .gt. -2) then
         if (nmax11 .eq. 1) ncoord=2
         zbisc=.false.
      else
         zbisc=.true.
         zembed=.true.
      endif
      if (zfit .and. max(jrot,ipar) .gt. 0) write(6,1005) itra0
 1005 format(  5x,'Ground state data read from       stream ',&
             'itra0 =',i4)
       write(6,1006) ilev
 1006 format(  5x,'Ground state energy read from     stream ',&
             'ilev   ',i4)

!     set other parameters to maximum values

      npnt1= nmax11
      npnt2= nmax12
!
      npnt=max(npnt1,npnt2)
!
!     read in of integration parameters and states to be considered
!     zeros will give the default values
!
      read(5,108) lpot,npropin,nprt,nv1
 108  format(5i5)
 101  format(20i5)
      if (npropin .le. 0) npropin=1
      if (npropin .gt. mxprop) then
        write(6,9999) npropin,mxprop
9999    format(/,5x,'npropin =',i3,' requested, but mxprop =',i3,&
               /,5x,' ***** stop *****')
      endif
      if (nprt .lt. npropin) nprt=npropin
      if (lpot.eq.0) lpot = lmax+jrot+mod(lmax+jrot,2)
      read(5,101) (iprop(i),i=1,npropin)
      if (iprop(1) .le. 0) then
         do 160 i=1,npropin
         iprop(i)=i
160      continue
      endif

      write(6,201) npropin,nprt,lpot,(iprop(i),i=1,npropin)
201   format(/,5x,'Number of properties to be considered   =',i5,&
              /,5x,'Number of properties computed in PROP   =',i5,&
              /,5x,'Number of angular DVR points            =',i5,&
              /,5x,'properties to be considered are numbered:',&
              /,5x,(25i5))

      if (zprint) write(6,*) '     Debug printing requested'
!
!     are we doing a non coriolis problem?
!     set up jk for the two cases and total basis set sizes
!
      if (jrot.eq.0) kmin1= 1
      if (zncor) then
         jk1= 1
         jrot=abs(jrot)
      else
         jk1= jrot + kmin1
      endif
      if (idia .gt. -2) then
         nrade=npnt1*npnt2
         nrado=nrade
         mbass=nrade*jk1*lmax
      else
         nrade=npnt1*(npnt1+1)/2
         nrado=npnt1*(npnt1-1)/2
         jt=jk1/2
         mbass=(nrade+nrado)*jt
         if (2*jt .ne. jk1) then
           if (ipar .eq. 0) mbass=mbass+nrade
           if (ipar .ne. 0) mbass=mbass+nrado
         endif
         mbass=mbass*lmax
       endif

!     reset number of vib-rot functions to be considered

!     write out data from triatom/rotlev runs

      write(6,300) mbass,nmax11,nmax12,lmax,jrot,ipar,kmin1,neval,nv1
300   format(5x,'total number of basis functions   = ',i5,/&
             5x,'number of r1 radial dvr points    = ',i5,/&
             5x,'number of r2 radial dvr points    = ',i5,/&
             5x,'number of angular functions       = ',i5,/&
             5x,'total angular momentum            = ',i5,/&
             5x,'ipar, vibrational parity          = ',i5,/&
             5x,'kmin= 1-p, where j + p is parity, = ',i5,/&
             5x,'total vib-rot functions available = ',i5,/&
             5x,'vib-rot functions taken (0 = all) = ',i5,/)

!     space for angular integration
      if (jrot == 0) then
         npot=lmax
      else
         ltop = lpot
         if (idia .eq. 2 .and. mod(ltop,2) .ne. ipar) ltop=ltop+1
         npot=((ltop+2)/2)*2
         nn2= npot/2
      endif
      return
      end

!**************************************************004
      subroutine core
!
      implicit double precision (a-h,o-y), logical (z)
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev

      integer, dimension (jk1) :: nbass1,lmin1

!     generate the subindex arrays needed for trans

      call genind(nbass1,lmin1,mbass,jk1,nbmax1,iket)

      call x3main(nbass1,lmin1)

      return
      end

!**************************************************005
      subroutine genind(nbass,lmin,mbass,jk,nbmax,ivec)
!
!     this subroutine sets up the sub-index arrays needed for trans
!     and reads in the basis function labels for the vectors.
!     it calculates nbmax, the largest value of nbass, neeeded to
!     dimension the space needed for the d-coefficients.
!
      implicit double precision(a-h,o-y),logical (z)

      integer, dimension (jk) :: nbass,lmin,lbass

      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      if (zprint) write(6,205) jk,mbass
205   format(//,'   j + kmin =',i3,'   mbass=',i6,/)

!     read in the basis function labels
      read(ivec) mbass0,lmin,lbass,nbass
      if (mbass0.gt.mbass) goto 999

!     generate the sub-index arrays and  find nbmax

      nbmax=nbass(1)
      do 2 k=2,jk
      nbmax=max(nbmax,nbass(k))
2     continue
      if (zprint) then
         write(6,201)
201      format(//,5x,'indices generated by genind',/)
         write(6,202) nbmax
202      format(5x,'nbmax= ',i5,/,5x,'nbass, lmin and lbass follow',/)
         write(6,*) (nbass(k), k=1,jk),(lmin(k), k=1,jk),&
                    (lbass(k), k=1,jk)
      endif
      return
999   write(6,200) ivec,mbass,mbass0
200   format(//,5x,'basis function dimensions in error',&
              /,5x,'unit =',i3,' mbass =',i5,' expected, =',i5,' found')
      stop
      end

!**************************************************006
      subroutine x3main(nbass1,lmin1)

!     in this part, all data etc for the ket are labelled 1;
!     all data etc for the bra are labelled 2.
!
      implicit double precision (a-h,o-y), logical (z)
      common /timing/ itime0
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev

      integer, dimension (jk1) :: nbass1,lmin1
!      double precision, dimension(neval) :: eh
      double precision, allocatable :: eh(:)
!      double precision, dimension(npot) :: xd,wtd,xalf
      double precision, allocatable :: xd(:),wtd(:),xalf(:)
!      double precision, dimension(npropin) :: prin,trin
      double precision, allocatable :: prin(:),trin(:)
!     arrays for lagpt
!      double precision, dimension(nrade*npot,npropin) :: d0
      double precision, allocatable :: d0(:,:)
!      double precision, dimension(npnt1) :: r1
      double precision, allocatable :: r1(:)
!      double precision, dimension(npnt2) :: r2
      double precision, allocatable :: r2(:)
!     arrays for trans
!      double precision, dimension(neval,npropin) :: props
      double precision, allocatable :: props(:,:)
!      double precision, dimension(neval*nbmax1/lmax*npot) :: dc1
      double precision, allocatable :: dc1(:)
!     array for dsrd
!      integer, dimension (max(npot,lmax)) :: iv
      integer, allocatable  :: iv(:)
      double precision, allocatable :: xp(:)
      integer, allocatable          :: ip(:)
      data x0/0.0d0/,xp5/0.5d0/
      !
      !
      ! s.y. added the xp paremeters, reading them from disk.
      !
      allocate(xp(npropin),ip(npropin))
      open(unit=32,status='old',form='formatted',file='pot.fit')
      read(32,*) npropin_t
      !
      if (npropin_t/=npropin) then
        !
        write(6,"('x3main: inconcistent number of parameters in unit=32 '2i8)") npropin_t,npropin
        stop 'inconcistent npropin in unit=32'
        !
      endif
      !
      do i = 1,npropin
         read(32,*) ip(i), xp(i)
      end do
      !
      close(unit=32)
      !
      allocate(eh(neval))
      allocate(xd(npot),wtd(npot),xalf(npot))
      allocate(prin(npropin),trin(npropin))
      allocate(d0(nrade*npot,npropin))
      allocate(r1(npnt1))
      allocate(r2(npnt1))
      allocate(props(neval,npropin))
      allocate(dc1(neval*nbmax1/lmax*npot))
      allocate(iv(max(npot,lmax)))
!     read in radial dvr grid points
      !
      call getrow(r1,npnt1,iket)
      if (idia .gt. -2)  call getrow(r2,npnt2,iket)
      if (jrot==0) then
        if (idia > -2) then
          read(iket) xalf
          do i=1,lmax
            xalf(2*lmax+1-i) = -xalf(i)
          end do
          read(iket)
          read(iket)
        else
          read(iket) xalf
          read(iket)
          read(iket)
          read(iket)
          read(iket)
        endif
      endif

      if (kmin1==0.and.jrot==1) then
          read(iket)
          read(iket)
          read(iket)
          read(iket)
          read(iket)
      endif

!     read in of energies for the ket

      read(iket) neval
      read(iket) eh

! J = 0
      if (ipar==0.and.jrot==0.and.zprint) then
         write(6,201)
201     format(//5x,'energies for the ket in cm-1',/)
         write(6,203) (eh(i)*219474.63,i=1,neval-1)
203     format(5d24.12,/)
      end if

!     zero properties array
      props = x0

      if (jrot==0) then
!     J=0 case
          call dsrd(dc1,iket,mbass,nbass1(1),jk1,iv,ibass)

          if (idia==-2) then
            iadd=ipar
          else
            i0=npnt1
          endif
        ii=0
        do i2=1,npnt2
        if (idia==-2) then
          rr2=r1(i2)
          i0=i2-iadd
        else
          rr2=r2(i2)
        endif
        do i1=1,i0

        if (idia==-2) then
          do 20 j=1,npot
          if(iv(j)==0) go to 20
          ii=ii+1
          call prop(ip,xp,prin,r1(i1),rr2,xalf(j),npropin,nprt)
          call prop(ip,xp,trin,rr2,r1(i1),xalf(j),npropin,nprt)
          do n=1,npropin
            d0(ii,n)=xp5*(prin(n)+trin(n))
          end do
  20     continue
        else
          do 21 j=1,lmax
          if(iv(j)==0) go to 21
          ii=ii+1
          call prop(ip,xp,prin,r1(i1),rr2,xalf(j),npropin,nprt)
          if (idia == 2) then
             call prop(ip,xp,trin,r1(i1),rr2,-xalf(j),npropin,nprt)
             do n=1,npropin
               d0(ii,n)=xp5*(prin(n)+trin(n))
             end do
          else
             do n=1,npropin
               d0(ii,n)=prin(n)
             end do
          endif
  21      continue
        end if
        end do
        end do

!       call to trans
        write(6,603)
603    format(/,5x,'trans called after')
        call system_clock(itime2,irate2,imax2)
        itime=(itime2-itime0)/irate2
        write(6,111) itime
        call trans(props,d0,dc1,ibass)

      else
!       J > 0 case
!       first generate angular grid and then call props

! Kolya
!  print *,'x3main: '
!  print *, d0(1,1),r1,r2,xd,wtd

        call lagpt(ip,xp,d0,r1,r2,xd,wtd)

        do 10 k1= 1,jk1
        if(nbass1(k1).eq.0) goto 10
        ndbass = (nbass1(k1)*npot)/lmax
        ibass = nbass1(k1)*neval
        kz = k1 -  kmin1

        call dsrd(dc1,iket,mbass,ndbass,jk1,iv,ibass)
        call rest(dc1,ndbass,jk1,xd,wtd,kz,lmin1)

!     call to trans
        write(6,603)
        call system_clock(itime2,irate2,imax2)
        itime=(itime2-itime0)/irate2
        write(6,111) itime
        ibasst = ibass/neval
        call trans(props,d0,dc1,ndbass)
10      continue
      end if

!     end of property calculation

      write(6,604)
604   format(/,5x,'outp called after')
      call system_clock(itime2,irate2,imax2)
      itime=(itime2-itime0)/irate2
      write(6,111) itime
111    format(/i10,' secs cpu time used'/)

      call outpt2(props,eh)

      return
      end

!**************************************************011
      subroutine lagpt(ip,xp,d0,r1,r2,xd,wtd)

!     subroutine lagpt obtains values of the dipole at the radial
!     dvr points and angular integration points

      implicit double precision(a-h,o-y), logical (z)
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar

      double precision :: xp(npropin)
      integer          :: ip(npropin)
      !
      double precision, dimension(nrade*npot,npropin) :: d0
      double precision, dimension(npnt1) :: r1
      double precision, dimension(npnt2) :: r2
      double precision, dimension(npot) :: xd,wtd,b,c
      double precision, dimension(npropin) :: prin

      data x0/0.0d0/,toler/1.0d-8/,&
           x1/1.0d0/,x2/2.0d0/,x3/3.0d0/,x4/4.0d0/

!     set up points & weights for npot point angular integration

      xnu= x0
      alf= xnu
      bta= alf
      do i=2,npot
        xi= dble(i)
        b(i)= (alf+bta)*(bta-alf)/&
            ((alf+bta+x2*xi)*(alf+bta+x2*xi-x2))
        c(i)= x4*(xi-x1)*(alf+xi-x1)*(bta+xi-x1)*(alf+bta+xi-x1)/&
            ((alf+bta+x2*xi-x1)*(alf+bta+x2*xi-x2)*&
             (alf+bta+x2*xi-x2)*(alf+bta+x2*xi-x3))
      end do

      call jacobi(npot,nn2,xd,wtd,alf,bta,b,c,csa,tsa)

      write(6,1000) npot,0,(xd(i),wtd(i),i=1,nn2)
 1000 format(//,i8,' point gauss-associated legendre integration',&
                  ' with k =',i2,&
             //,5x,'integration points',11x,'weights',&
             /,(f23.15,d25.12))
      write(6,1010) csa,tsa
 1010 format(/,5x,'computed sum of weights',d26.15,&
             /,5x,'exact    sum of weights',d26.15//)
      if (abs((csa-tsa)/tsa) .gt. toler) goto 930
!     define other integration points
      do i=1,nn2
        j=i+nn2
        xd(j)=-xd(nn2+1-i)
        wtd(j)=wtd(nn2+1-i)
      end do
      do i = 1, 2*nn2
        xd(i) = -xd(i)
      end do

!     calculate properties at (r1,r2,cos\theta)

      i0=npnt1

      icall = 0
      ii = 0
      do i2=1,npnt2
         if (idia .eq. -2) then
           rr2=r1(i2)
           i0=i2
         else
           rr2=r2(i2)
         end if
        do i1=1,i0
          do j=1,npot
            ii=ii+1
            icall = icall + 1
            call prop(ip,xp,prin,r1(i1),rr2,xd(j),npropin,nprt)
            do n=1, npropin
              d0(ii,n) = prin(n)
            end do
          end do
        end do
      end do

      return
  930 write(6,940)
  940 format(/,5x,'gauss-legendre weights in error: adjust algorithm')
      stop
      end

      subroutine jacobi(nn,nn2,x,a,alf,bta,b,c,csa,tsa)

!     calculates zeros x(i) of the nn'th order jacobi polynomial
!     pn(alf,bta) for the segment (-1,1) & and corresponding weights
!     for gauss-jacobi integration. this routine uses a brute force
!     search for zeros due to GJ Harris (2001).
!     note that for our purposes, alf= bta= nu.

      implicit double precision(a-h,o-z)
      double precision, dimension(nn) :: x,a,b,c,xt
      data x0/0.0d0/,x1/1.0d0/,x2/2.0d0/,x6/6.0d0/,eps/1.0d-12/
      fn= dble(nn)
      csa= x0
      beta= x1
      if (alf.eq.1) beta = beta/x6
      cc= beta*x2**(alf+bta+x1)
!      tsa= cc/x2
      tsa= cc
      do 15 i=2,nn
        cc= cc*c(i)
   15 continue

! step through the jacobi polynomial and
! notes a first guess for the posititions of the zeros in the
! array xt(ii). Zeros are found by looking for a change in sign.

      ii=0
      pm1=x1
      xxx=x1
 30   continue
      call recur(p,dp,pn1,xxx,nn,alf,bta,b,c)

      if (pm1*p .lt. x0) then
         pm1 = -pm1
         ii = ii +1
         xt(ii)=xxx
      endif

      if (ii .eq. nn2) then
         do 40 i=1,nn2
         call recur(ptemp,dp,pn1,xt(i),nn,alf,bta,b,c)
40      continue
      else
         xxx=xxx-0.0001
         if (xxx .gt. -0.0002) goto 30
         write(6,*) "Incorrect number",ii-1," of zeros found in JACOBI"
         stop
      endif

       do 20 i=1,nn2
       call root(xt(i),nn,alf,bta,dpn,pn1,b,c,eps)
       x(i)= xt(i)
       a(i)= cc/(dpn*pn1)
       csa= csa + a(i) + a(i)
  20  continue
      if (2*nn2 .ne. nn) csa=csa-a(nn2)
      return
      end

!***********************************************019
      subroutine root(x,nn,alf,bta,dpn,pn1,b,c,eps)

!     improves the approximate root x; in addition obtains
!          dpn = derivative of p(n) at x
!          pn1 = value of p(n-1) at x.

      implicit double precision(a-h,o-z)
      double precision, dimension(nn) :: b,c
      iter= 0
1     iter= iter + 1
      call recur(p,dp,pn1,x,nn,alf,bta,b,c)
      d = p/dp
      x = x - d
      if (abs(d) .gt. eps .and. iter .lt. 10) goto 1
      dpn= dp
      return
      end

!**************************************************020
      subroutine recur(pn,dpn,pn1,x,nn,alf,bta,b,c)
      implicit double precision(a-h,o-z)
      double precision, dimension(nn) :: b,c
      data x0/0.0d0/,x1/1.0d0/,x2/2.0d0/
      p1= x1
      p= x + (alf-bta)/(alf + bta + x2)
      dp1= x0
      dp= x1
      do 1 j=2,nn
      q= (x - b(j))*p - c(j)*p1
      dq= (x - b(j))*dp + p - c(j)*dp1
      p1= p
      p= q
      dp1= dp
      dp= dq
1     continue
      pn= p
      dpn= dp
      pn1= p1
      return
      end

!*******************************************************020
      subroutine dsrd(d,ivec,mmbass,nbass,jk,iv,ibass)

!     subroutine to read d coefficients from dstore data
!     if read directly from dvr3d (i.e. j .le. 0), then transform to an
!     fbr in theta

      implicit double precision (a-h,o-y), logical(z)

      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      double precision, dimension(neval,nbass) :: d
      double precision, dimension(nbass) :: temp,plegd
      integer, dimension(200)  :: iv

      if (jk.gt.1) then
         !
         ! s.y. changed the order
         !
         nread = (nbass*lmax)/npot
         d = 0.0d0
         read(ivec) ((d(i,j),j=1,nread),i=1,neval)
!        orthonormality test
!	 sumd11 = 0.0d0
!	 do i=1,nread
!	   sumd11 = sumd11 + d(1,i)*d(1,i)
!	 end do
!	 write(6,*) 'dsrd: nread=', nread,'sumd11=', sumd11
      else
         rewind ivec
         do 30 i=1,6
         read(ivec)
30       continue
         if (.not. zbisc) read(ivec)
         if(idia.gt.-2) then
           rewind ivec
           do i=1,7
             read(ivec)
           end do
           read(ivec) kz,maxleg,nidvr,lincr
         else
           read(ivec) kz,maxleg,nidvr
         end if

         leng=(maxleg+1)*nidvr
         call getrow(plegd,leng,ivec)
         if (zbisc) then
           read(ivec) iang,ibass
           read(ivec) (iv(ii),ii=1,nidvr)
         else
            iang=nidvr
            ibass=mmbass
            do i=1,nidvr
              iv(i)=1
            end do
         end if
         read(ivec)
         read(ivec)
         do l=1, neval
           call getrow(temp,ibass,ivec)
           do i=1, ibass
             d(l,i)=temp(i)
           end do
         end do
      end if

      return
      end

!******************************************************020a
      subroutine rest(d,nbass,jk,xd,wtd,kz,lmin1)

!     rest of dsrd
!     subroutine to read d coefficients from dstore data
!     if read directly from dvr3d (i.e. j .le. 0), then transform to an
!     fbr in theta

      implicit double precision (a-h,o-y), logical(z)

      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      integer, dimension (jk) :: lmin1
      double precision, dimension(neval,nbass) :: d,temp
      double precision, dimension(npot,lmax*max(1,idia)) :: plegd
      double precision, dimension(npot) :: xd,wtd

      data x0/0.0d0/,x1/1.0d0/

      temp = d
      d = x0

      nrad = nbass/npot
      jdia=max(1,idia)
      nang = lmax*jdia
      nang2 = lmax
      npot2 = npot/jdia
!     calculate the legendre function
      if (idia>-2.and.jk==0) kz=0
      call asleg(plegd,d,nang-1,xd,npot,kz)
      if (idia==2) then
        do i=1,npot
          scale = sqrt(wtd(i))
          jj = lmin1(kz+1) - kz + 1
          do j=1,nang2
            plegd(i,j) =  scale*plegd(i,jj)
            jj=jj+jdia
          end do
        end do
      else
        do i=1,npot
          scale = sqrt(wtd(i))
          do j=1,lmax*max(1,idia)
            plegd(i,j) = scale*plegd(i,j)
          end do
          end do
        end if
!        orthonormality test
!         do k1=1,nang2
!         do k2=1,nang2
!         sum = x0
!         do j=1,npot
!         sum = sum + plegd(j,k1)*plegd(j,k2)
!         end do
!         write(6,*) 'rest: k1,k2 =',k1,k2,'sum plegd  = ', sum
!         end do
!         end do

!     evaluate wavefunction at angular integration points

      beta=x1
      d = x0
      ipt=1
      jpt=1

      do  mn=1,nrad
        call dgemm('n','t',neval,npot,nang2,beta,temp(1,jpt),neval,plegd,&
                npot,beta,d(1,ipt),neval)
        ipt = ipt + npot
        jpt = jpt + nang2
      end do

      return
      end

!******************************************************021
      subroutine asleg(pleg,pnorm,lmax,x,nn2,m)

!     calculate polynomials 1 to lmax at x = cos(theta) for m = 0 or 1,
!     using the routine of press et al, numerical recipes, p. 182,
!     for the polynomial part of associated legendre functions.
!     a factor of sin(theta)**m has NOT been removed from all functions.

      implicit double precision (a-h,o-z)

      double precision, dimension(nn2,0:lmax) :: pleg
      double precision, dimension(nn2) :: x
      double precision, dimension(0:lmax) :: pnorm

      data x1/1.0d0/,x2/2.0d0/

      if (m .lt. 0) goto 999
      do 10 i=1,nn2
      if (abs(x(i)) .gt. x1) goto 999
      pmm = x1
      fact = x1
      somx2=sqrt((x1-x(i))*(x1+x(i)))
      do j=1,m
        pmm = -pmm * fact * somx2
        fact = fact + x2
      end do
      pleg(i,0)= pmm
      pmmp1= x(i)*dble(m+m+1)*pmm
      pleg(i,1)= pmmp1
      ll=1

      do l= 2+m,lmax+m
        r2lm1 = dble(l+l-1)
        rlpmm1= dble(l+m-1)
        rlmm  = dble(l-m)
        pll= (x(i)*r2lm1*pmmp1 - rlpmm1*pmm)/rlmm
        pmm= pmmp1
        pmmp1= pll
        ll=ll+1
        pleg(i,ll)= pll
      end do
10    continue

!     set up the normalisation constants
!     (pnorm)**2 = (2j + 1)/2   *   (j - k)! / (j + k)!
      jj = -1
      do 13 j = m,lmax+m
      fact = x1
      do i = j-m+1,j+m
        facti = dble(i)
        fact = fact * facti
      end do
      jj = jj + 1
      pnorm(jj) = sqrt(dble(j+j+1) / (fact + fact))
   13 continue
!     now normalise the polynomials
      do jj=0,lmax
        do i=1,nn2
          pleg(i,jj) = pleg(i,jj) * pnorm(jj)
        end do
      end do
      return
999   write(6,200)
200   format(//5x,'improper argument in subroutine asleg'/)
      stop
      end

!*********************************************************022

      subroutine trans(t,d0,dc1,ndbass)

!     subroutine trans is the main working routine in program xpect3
!     it carries out the necessary angular integrations using 3-j
!     symbols as defined by brink and satchler "angular momentum".
!     trans is called by x3main for each k-k' overlap integral.

!     it uses a fast matrix multiplication routine mxmb.

      implicit double precision (a-h,o-y), logical (z)
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      double precision, dimension(neval,npropin) :: t
      double precision, dimension(npot*nrade,npropin) :: d0
      double precision, dimension(neval,ndbass) :: dc1

      nrad = ndbass/npot
      if (zprint) write(6,1) ndbass, nrad
1    format('  trans: ndbass, nrad = ',2i6)

!    take square of the wavefunction
     dc1=dc1**2

!    orthonormality test
!     sumdc111 = 0.0d0
!     do i=1,ndbass
!       sumdc111 = sumdc111 + dc1(1,i)
!     end do
!     write(6,*) 'trans: sum dc1(1,i) 1-1,',ndbass,'=',sumdc111

      if (jrot==0) then
        do i=1, neval
          do n=1, npropin
            do j=1, ndbass
              t(i,n) = t(i,n) + d0(j,n)*dc1(i,j)
            end do
          end do
        end do
      else

        n0 = npnt1
        i0 = 1
        if(nrad==nrade) i0 = 0
        do i=1, neval
          do n=1, npropin
            ic = 0
            ip = 0
            do n2 = 1,npnt2
              if(idia==-2) n0 = n2 - i0
              do n1 = 1,n0
                do j = 1,npot
                  ic = ic + 1
                  ip = ip + 1
                  t(i,n) = t(i,n) + d0(ip,n)*dc1(i,ic)
                end do
              end do
              if(i0==1) ip = ip + npot
            end do
          end do
        end do

      end if

      return
      end

!**************************************************026a
      subroutine outpt2(props,eh)
!
!     outpt2 outputs the data at the end of a property calculation
!
      implicit double precision(a-h,o-y), logical(z)
      common/dim/ ncoord,npnt,npnt1,npnt2,nrade,nrado,&
                  lpot,npot,nbmax1,mbass,kmin1,jk1,neval,&
                  nn2,lmax,npropin,nprt,jrot,idia,ipar
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      common /logic/ zmors1,zncor,zprint,ztra,zstart,zmors2,&
                     zfit,zform,iket,itra,itra0,ilev
      parameter (mxprop=100)
      common /pot/ iprop(mxprop)
      double precision, dimension (neval) :: eh
      double precision, dimension (neval,npropin) :: props,xprop
      double precision, dimension (npropin) :: pground

      character(len=2) s1
!     autocm converts atomic units to wavenumbers
      data autocm/ 2.19474624d+05/,x0/0.0d0/,s1/'  '/

      if(zncor) then
         j1= jk1
      else
         j1= jk1-kmin1
      endif

      write(6,203) jrot, kmin1, idia, ipar
203   format(//5x,'j1=',i4,'  kmin1=',i4,&
                 '  idia=',i4,'  ipar=',i4,//)

      ezero2=x0
      read(ilev,*,end=510,err=510)
      read(ilev,*,end=510,err=510)
      read(ilev,505,end=510,err=510) ezero2
  505 format(d20.12)
  510 continue
      ezero = ezero2 * autocm
      write(6,1005) ezero2
 1005 format(5x,'ground zero in hartrees ',d16.8)
      write(6,204) ezero
204   format(5x,'ground zero =',d16.8,' cm-1')

!
      if (ipar.eq.0.and.kmin1.eq.1) then
         s1 = 'A1'
      else if (ipar.eq.1.and.kmin1.eq.0) then
         s1 = 'A2'
      else if (ipar.eq.1.and.kmin1.eq.1) then
         s1 = 'B1'
      else if (ipar.eq.0.and.kmin1.eq.0) then
         s1 = 'B2'
      end if

      if (ztra .and. zform) then
          open(unit=itra,form='formatted')
      else if (ztra) then
          open(unit=itra,form='unformatted')
      endif
      write(6,305)
305  format(//,' J     ie       energy   property')

      if (ipar==0.and.jrot==0.and.zfit) then
! Kolya
        do ie=2,neval
          do n=1,npropin
            xprop(ie,n) = props(ie,n) - props(1,n)
          end do
        end do
        do n=1,npropin
          write(itra0,*) props(1,n)
        end do

        !
        ! s.y. change: now we write absolute values of derivatives, all of them!
        !
        do 4 ie=1,neval
        xeh=eh(ie)-ezero2
        xe1 = xeh*autocm
        if (ztra) write(itra,308) jrot,s1,ie,xe1,xeh
308     format(i2,1x,a2,i4,3x,f10.3,5(2x,d13.6)/&
            (17x,5(2x,d13.6)))
        do n=1, npropin
           !
           ! s.y. changed xprop to props
           !
           write(12,*) props(ie,n)
        end do
        write(6,308) jrot,s1,ie,xe1,(xprop(ie,n),n=1,npropin)
 4      continue
        !
        continue
      else if (zfit) then

!   Read in ground state shifts
         if (zform) then
            open(unit=itra0,form='formatted')
            read(itra0,*) pground
         else
            open(unit=itra0,form='unformatted')
            read(itra0)
            read(itra0) pground
         endif
         close (unit=itra0)

         do ie=1,neval
            do n=1,npropin
            xprop(ie,n) = props(ie,n) - pground(n)
            end do
         end do

         do ie=1,neval
           xeh= eh(ie)-ezero2
           xe1 = xeh*autocm
           if (zform) then
              if (ztra) write(itra,308) jrot,s1,ie,xe1,xeh
              do n = 1,npropin
                !
                ! s.y.
                !
                if(ztra) write(itra,*) props(ie,n)
                !
              end do
           else
              if (ztra) write(itra) jrot,s1,ie,xe1,xeh
              do n = 1,npropin
                !
                ! s.y.
                !
                if(ztra) write(itra,*) props(ie,n)
                !if(ztra .and. xprop(ie,n).ne.x0) write(itra) props(ie,n)
              end do
          endif
          write(6,308) jrot,s1,ie,xe1,(xprop(ie,n),n=1,npropin)
        end do

      else

        do ie=1,neval
          xeh=eh(ie)-ezero2
          xe1 = xeh*autocm
          write(6,308) jrot,s1,ie,xe1,(props(ie,n),n=1,npropin)
        enddo

!     write out to stream itra?

        if (.not. ztra) return

          if (idia .gt. 0) then
            isym=0
          else
            isym=abs(idia)
            if (ipar .eq. 1) isym=-isym
          endif

        if (zform) then
!         is this the first write-out?
          if(.not.zstart) then
10          read(itra,*,end=90)
            goto 10
90         continue
! *****   inclusion of the following card is machine dependent *****
            backspace itra
          endif
!
          write(itra,701) jrot,neval,ipar,kmin1
          write(itra,*) npropin,max(idia,0),isym,ezero
          write(itra,*) (iprop(n),n=1,npropin)
          write(itra,*) eh*autocm
701      format(i4,1x,i4,1x,i4,1x,i4)
          do n=1,npropin
            call outro2(props(1,n),neval,itra)
          end do

        else
!         is this the first write-out?
          if(.not.zstart) then
15        read(itra,end=95)
          goto 15
95        continue
! *****   inclusion of the following card is machine dependent *****
          backspace itra
        endif
!
        write(itra) jrot,neval,ipar,kmin1
        write(itra) npropin,max(idia,0),isym,ezero
        write(itra) (iprop(n),n=1,npropin)
        write(itra) eh*autocm
        do n=1,npropin
          call outrow(props(1,n),neval,itra)
        end do
      endif
      endif

      close (unit=itra)
      return
      end
!                                                **027
      subroutine getrow(row,nrow,iunit)
!
      implicit double precision (a-h,o-y)
      double precision, dimension(nrow) :: row
      read(iunit) row
      return
      end
!                                                **028
      subroutine outrow(row,nrow,iunit)
!
      implicit double precision (a-h,o-y)
      double precision, dimension(nrow) :: row
      write(iunit) row
      return
      end
!                                                **028a
      subroutine outro2(row,nrow,iunit)
!
      implicit double precision (a-h,o-y)
      double precision, dimension(nrow) :: row
      write(iunit,*) row
      return
      end
!                                                **029
      subroutine dprint(dz0,dz1,nrad,lpot,jdia,npropin)
!
!     subroutine to print out the property integrals calculated
!     in subroutine lagpt for debugging purposes.
!
      implicit double precision (a-h,o-y), logical (z)
      common /mass/ xmass(3),g1,g2,zembed,zbisc
      double precision, dimension (nrad,npropin) :: dz0
      double precision, dimension (nrad,npropin,lpot) :: dz1
      write(6,200)
200   format(//,5x,'p matrices calculated by subroutine lagpt',/)
      do 1 n=1,npropin
      do 11 i=1,nrad
      write(6,*) dz0(i,n)
      write(6,*) (dz1(i,n,k),k=jdia,lpot,jdia)
11    continue
 1    continue
      return
      end


text:
! Name = ARPIT KUMAR JAIN ROLLNO = 180122009
PROGRAM question4
    IMPLICIT NONE

    INTEGER :: n, factorial = 1, i
    PRINT *, "Enter number: "

    ! reading integer from screen
    READ *, n

    ! finding factorial by looping 1 to n and storing multiplication on factorial variable
    DO i = 1, n
        factorial = factorial * i
    ENDDO

    ! printing factorial value on screen
    print *, factorial

    STOP
END PROGRAM

max_stars_repo_path: 2. Evaluation-2/4.f90

max_issues_repo_path: 2. Evaluation-2/4.f90

max_forks_repo_path: 2. Evaluation-2/4.f90


text:
module simulation_mod
    use sim_setup_mod
    use atp_mod
    use mps_mod
    use atp_material_removal_mod
    use types_mod
    use convert_mpar_mod
    use output_mod
    use usr_interface_mod
    implicit none

    private ! by default
    public  :: setup_simulations
    public  :: simulate, simulate_gradient
    public  :: optanalyzer

    integer :: mpar_fid
    integer :: eext_fid
    integer :: evec_fid

    contains
! Main simulation routine
subroutine simulate(error, xvar, f_data, evec)
    implicit none
    double precision, intent(out)   :: error
    double precision, intent(in)    :: xvar(:)
    type(fdata_typ), intent(inout)  :: f_data
    double precision, optional, allocatable, intent(out) :: evec(:)
    type(evectmp_typ), allocatable  :: ev_tmp(:)
    integer                         :: k1, stype, nerr1, nerr2
    double precision, allocatable   :: mpar(:)
    double precision                :: error_tmp
    logical                         :: failed
    procedure(user_sim_template),pointer:: usr_sim ! Address to user subroutine

    failed = .false.

    ! Convert xvar to mpar, where mpar includes all material parameters
    allocate(mpar(size(f_data%glob%ipar_init)))
    call xvar_to_mpar(mpar, xvar, f_data%glob)

    ! Allocate ev_tmp if needed
    if (present(evec)) then
        allocate(ev_tmp(size(f_data%sim)))
        nerr1 = 0
    endif

    error = 0.d0
    do k1=1,size(f_data%sim)
        stype = f_data%sim(k1)%stype
        if      (stype == 0) then
            ! External simulation
            if (present(evec)) then
                call simulate_external(error_tmp, mpar, f_data, k1, ev_tmp(k1)%evec_tmp)
            else
                call simulate_external(error_tmp, mpar, f_data, k1)
            endif
        elseif  (stype == -1) then
            ! User simulation
            usr_sim => f_data%sim(k1)%usr_sim%sim_addr

            if (present(evec)) then
                call usr_sim(error_tmp, mpar, f_data%sim(k1)%usr_sim%user_data, ev_tmp(k1)%evec_tmp)
            else
                call usr_sim(error_tmp, mpar, f_data%sim(k1)%usr_sim%user_data)
            endif
        elseif  (stype == 1) then
            ! ATP simulation
            if (present(evec)) then
                call atp_simulate(error_tmp, mpar, f_data, k1, ev_tmp(k1)%evec_tmp)
            else
                call atp_simulate(error_tmp, mpar, f_data, k1)
            endif
        elseif  (stype == 2) then
            ! MP simulation (MPS)
            if (present(evec)) then
                call mp_simulate(error_tmp, mpar, f_data, k1,  ev_tmp(k1)%evec_tmp)
            else
                call mp_simulate(error_tmp, mpar, f_data, k1)
            endif
        elseif (stype == 11) then
            ! ATP machine ("material removal") (no error output for this stage)
            call atp_material_removal(error_tmp, mpar, f_data, k1)
            if (present(evec)) then
                allocate(ev_tmp(k1)%evec_tmp(1))
                ev_tmp(k1)%evec_tmp = 0.d0
            endif
        else
            call write_output('Simulation '//int2str(k1)//': stype = '//int2str(stype)//' is not supported', 'error', 'sim')
        endif
        ! Add current error to total error
        error = error + error_tmp

        if (present(evec)) then
            nerr1 = nerr1 + size(ev_tmp(k1)%evec_tmp)
        endif

        if (error>(huge(1.d0)/10)) then
            failed = .true.
            ! No point in keep simulating if error is already infinity for current parameters
            exit
        endif

    enddo

    if (present(evec)) then
        allocate(evec(nerr1))
        nerr1 = 1
        do k1=1,size(f_data%sim)
            nerr2 = (nerr1 - 1) + size(ev_tmp(k1)%evec_tmp)
            evec(nerr1:nerr2) = ev_tmp(k1)%evec_tmp
            nerr1 = nerr2 + 1
        enddo
    endif

end subroutine simulate

! Simulation routine if numerical gradient should be obtained for optimization
subroutine simulate_gradient(error, gradient, xvar, f_data, evec, evec_gradient, resnr)
    implicit none
    double precision, intent(out)   :: error, gradient(:)
    double precision, intent(in)    :: xvar(:)
    type(fdata_typ), intent(inout)  :: f_data
    double precision, optional, allocatable, intent(out) :: evec(:), evec_gradient(:,:)
    integer, optional               :: resnr  ! Save results for the base simulation (without pertubation?) if >0
                                        ! This options overrides whatever is in f_data%glob%resnr for the base sim
    double precision                :: error_pert, pertubation, xvar_pert(size(xvar))
    integer                         :: k1
    double precision, allocatable   :: evec_pert(:)
    integer                         :: tmp_resnr

    if (present(resnr)) then
        tmp_resnr = f_data%glob%resnr
        f_data%glob%resnr = resnr
    endif

    ! Get error at current xvar
    if (present(evec)) then
        call simulate(error, xvar, f_data, evec)
        if (present(evec_gradient)) then
            allocate(evec_gradient(size(evec), size(xvar)))
            allocate(evec_pert(size(evec)))
        endif
    else
        call simulate(error, xvar, f_data)
    endif

    if (present(resnr)) then
        f_data%glob%resnr = tmp_resnr
    endif

    xvar_pert = xvar
    pertubation = f_data%glob%num_grad_pert
    do k1=1,size(xvar)
        xvar_pert = xvar
        xvar_pert(k1) = xvar_pert(k1) + pertubation
        if (present(evec)) then
            call simulate(error_pert, xvar_pert, f_data, evec_pert)
            if (present(evec_gradient)) then
                evec_gradient(:,k1) = (evec_pert-evec)/pertubation
            endif
        else
            call simulate(error_pert, xvar_pert, f_data)
        endif
        gradient(k1) = (error_pert-error)/pertubation   !Numerical derivative
    enddo

end subroutine simulate_gradient

! Subroutine for setting up simulation (reading in expdata etc)
subroutine setup_simulations(f_data)
implicit none
    type(fdata_typ), intent(inout)  :: f_data
    integer                         :: k1, stype

    do k1=1,size(f_data%sim)
        stype = f_data%sim(k1)%stype
        if (stype==0) then      ! External simulation
            ! No setup
        elseif (stype==1) then  ! ATP simulation
            call setup_simulation(f_data, k1)
        elseif (stype==2) then  ! MPS
            call setup_simulation(f_data, k1)
        elseif (stype==11) then ! ATP material removal
            ! No setup
        endif
    enddo

end subroutine

! Subroutine for analyzing a found optimum
subroutine optanalyzer(xopt, f_data, error, dfdx, evec, corr)
    implicit none
    double precision                :: xopt(:), error
    type(fdata_typ)					:: f_data
    double precision, allocatable   :: dfdx(:), corr(:,:), evec(:), evec_gradient(:,:), evec_norm(:)
    integer                         :: k1, k2, nvar
    integer                         :: tmp_resnr

    nvar = size(xopt)

    ! allocate output depending on sizes of input
    allocate(dfdx(nvar), corr(nvar, nvar), evec_norm(nvar))

    ! Ensure that results are not written for each pertubation
    tmp_resnr = f_data%glob%resnr ! Save old setting
    f_data%glob%resnr = f_data%glob%opt_resnr

    ! Get numerical derivatives at optimum, asking for result output at the base simulation
	call simulate_gradient(error, dfdx, xopt, f_data, evec, evec_gradient, tmp_resnr)

    f_data%glob%resnr = tmp_resnr ! Reset old setting (incremented if a new result has been written)

    ! Calculate correlation
    do k1=1,nvar
        evec_norm(k1) = sqrt(sum(evec_gradient(:,k1)*evec_gradient(:,k1)))
    enddo

    corr = 1.d0 ! For diagonal
    do k1 = 1,nvar
        do k2 = (k1+1),nvar
            corr(k1,k2) = sum(evec_gradient(:,k1)*evec_gradient(:,k2))/(evec_norm(k1)*evec_norm(k2))
            corr(k2,k1) = corr(k1,k2)
        enddo
    enddo

end subroutine optanalyzer

! Subroutine for calling an external simulation (stype=0)
subroutine simulate_external(error, mpar, f_data, simnr, evec)
use compiler_mod
implicit none
    double precision                :: error, mpar(:), tmp
    type(fdata_typ)                    :: f_data
    integer                         :: simnr, status, nerr, k1
    character(len=(strl+10))        :: command  !Make slightly larger than the script field of f_data%sim(simnr)
    character(len=1)                :: eout_str !String for writing logical if all errors should be written
    double precision, allocatable, optional :: evec(:)
    character(len=strl)             :: mpar_file, efile, evecfile

    write(mpar_file, "(A)") 'mpar_'//int2str(simnr)//'.txt'
    write(efile, "(A)") 'eext_'//int2str(simnr)//'.txt'
    write(evecfile, "(A)") 'evec_'//int2str(simnr)//'.txt'

    if (present(evec)) then
        eout_str = '1'
    else
        eout_str = '0'
    endif

    ! Generate command line call with additional arguments:
    ! 1) eout_str (1 if evec is requested, 0 otherwise)
    ! 2) resnr (which result file number it currently is)
    ! 3) mpar_file (mpar_simnr.txt) gives the material parameters
    ! 4) efile (eext_simnr.txt) gives the file to which the scalar error should be written
    ! 5) evecfile (evec_simnr.txt) gives the file to which the error vector should be written (if requested)
    command = trim(f_data%sim(simnr)%ext_cmd%script) // ' ' // eout_str // ' ' // int2str(f_data%glob%resnr) // ' ' &
                                             // trim(mpar_file) // ' ' // trim(efile) // ' ' // trim(evecfile)

    if (f_data%glob%resnr > 0) then
        if (f_data%sim(simnr)%ext_cmd%logfile /= '') then
            command = trim(command) // ' >> ' // trim(f_data%sim(simnr)%ext_cmd%logfile) // ' 2>&1 '
        endif
    endif

    ! Write material parameters to default file 'mpar.txt'
    open(newunit=mpar_fid, file=trim(mpar_file), status='replace', action='write', IOSTAT=status)
    write(mpar_fid, "(ES20.10E3)") mpar
    close(mpar_fid)

    ! Make command line call
    status = system(command)   ! Option: call execute_command_line() or call system()

    if (status/=0) then
        call write_output('Error in external call "'//trim(command)//'"', 'error', 'sim', halt=.false.)
        call write_output('Output is '//int2str(status), 'error', 'sim', loc=.false.)
    endif

    ! Aquire error from external simulation
    open(newunit=eext_fid, file=efile, status='old', action='read', IOSTAT=status)
    if (status/=0) then
        close(eext_fid)
        call write_output('Cannot open the error file '//trim(efile)//' from the external script', 'error', 'sim', halt=.false.)
        call write_output('Please verify that this file has been created', 'error', 'sim', loc=.false.)
    endif

    rewind(eext_fid)    ! Ensure that file is at beginning before reading
    read(eext_fid, *, IOSTAT=status) error
    if (status/=0) then
        close(eext_fid)
        call write_output('Invalid format in the error file "'//trim(efile)//'"', 'error', 'sim', halt=.false.)
    endif

    close(eext_fid)

    ! Potentially also get error vector from external
    if (present(evec)) then
        open(newunit=evec_fid, file=evecfile, status='old', action='read', IOSTAT=status)
        if (status/=0) then
            close(evec_fid)
            call write_output('Cannot open the error file '//trim(evecfile)//' from the external script', 'error', 'sim', halt=.false.)
            call write_output('Please verify that this file has been created', 'error', 'sim', loc=.false.)
        endif

        rewind(evec_fid)    ! Ensure that file is at beginning before reading

        ! Get size of evec if first time:
        if (f_data%sim(simnr)%n_errors==-1) then
            nerr = -1 ! One will also be added after status=-1, hence start with -1
            do while(status==0)
                read(evec_fid, *, IOSTAT=status) tmp
                nerr=nerr+1
            enddo
            rewind(evec_fid)
        else
            nerr = f_data%sim(simnr)%n_errors
        endif
        allocate(evec(nerr))

        ! Read in error vector
        do k1=1,nerr
            read(evec_fid, *, IOSTAT=status) tmp
            if (status==0) then
                evec(k1) = tmp
            else
                call write_output('Problem reading error vector input file from external script on line '//int2str(k1), 'error', 'sim', halt=.false.)
                call write_output('Note that the same number of lines should exists every time', 'error', 'sim', loc=.false.)
            endif
        enddo

        ! Check that the entire file has been read:
        read(evec_fid, *, IOSTAT=status) tmp
        if (status==0) then
            call write_output('Not all lines in error vector input file read', 'warning', 'sim', halt=.false.)
            call write_output('Please ensure that the same number of lines exists every time', 'warning', 'sim', loc=.false.)
        endif

        close(evec_fid)
    endif

end subroutine simulate_external

end module simulation_mod

max_stars_repo_name: KnutAM/matmodfit

max_issues_repo_name: KnutAM/matmodfit

max_forks_repo_name: KnutAM/matmodfit


text:
fl

\ Run this script to initialize the file system

\
\
\ first partion has system information and partition files
\ it is a 32 Megabyte partition, should be more than
\ adequate for a lot of configuration files.
\
\ This information is generated automatically, after an sd_init
\ _sd_maxblock will be valid
\
[ifndef $C_a_doconl
    h52 wconstant $C_a_doconl
]
\
\ constant ( x -- ) skip blanks parse the next word and create a constant, allocate a long, 4 bytes
[ifndef constant
: constant
	lockdict create $C_a_doconl w, l, forthentry freedict
;
]

sd_init

1 constant _sysstart

\ _sd_CrEaTe ( n1 n2 -- ) n1 - starting block, n2 - last block + 1, CREATE a file system, WIPES OUT DATA
: _sd_CrEaTe
	sd_cogbuf h80 bounds
	do
		h20202020 i COG!
	loop
	h20202F01 dup sd_cogbuf COG!
	sd_cogbuf h20 + COG!
	h10000 sd_cogbuf h2A + COG!
	h81 sd_cogbuf h2B + COG!
	over sd_cogbuf h2E + COG!
	over sd_cogbuf h2F + COG!
	over sd_cogbuf h32 + COG!
	sd_cogbuf h33 + COG!
	dup h81 + sd_cogbuf h36 + COG!
	sd_lock	dup sd_blockwrite
	_sd_initdir
	sd_unlock
;

\ 32 Megabytes
h10000 _sysstart + constant _sysblocks

_sd_maxblock L@ constant _usrend

_sysstart _sysblocks _sd_CrEaTe
_sysstart _sysblocks + _usrend _sd_CrEaTe

_sysstart sd_mount

h2 fwrite .sdcardinfo

                      SDCard size:...

	_sd_maxblock L@ 2/ _nf C@++
	c" .sdcardinfo" sd_append
	c"  Kbytes~h0D~h0D" C@++ c" .sdcardinfo" sd_append

	c" File System: sys~h0D                      Start Block:" C@++ c" .sdcardinfo" sd_append

	_sysstart _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D                        End Block:" C@++ c" .sdcardinfo" sd_append

	_sysstart _sysblocks + _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D        File System Size (blocks):" C@++ c" .sdcardinfo" sd_append

	_sysblocks _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D        File System Size  (bytes):" C@++ c" .sdcardinfo" sd_append

	_sysblocks 2/ _nf C@++
	c" .sdcardinfo" sd_append

	c"  Kbytes~h0D~h0D" C@++ c" .sdcardinfo" sd_append

	c" File System: usr~h0D                      Start Block:" C@++ c" .sdcardinfo" sd_append

	_sysstart _sysblocks + _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D                        End Block:" C@++ c" .sdcardinfo" sd_append
	_usrend _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D        File System Size (blocks):" C@++ c" .sdcardinfo" sd_append
	_usrend _sysstart _sysblocks + - _nf C@++ c" .sdcardinfo" sd_append

	c" ~h0D        File System Size  (bytes):" C@++ c" .sdcardinfo" sd_append

	_usrend _sysstart _sysblocks + - 2/ _nf C@++
	c" .sdcardinfo" sd_append
	c"  Kbytes~h0D~h0D" C@++ c" .sdcardinfo" sd_append

h100 fwrite sdboot.f

hA state orC! cr

c" sdboot.f  -  initializing~h0D~h0D" .cstr
1 sd_mount

\ mountsys ( -- ) mount the system disk
[ifndef mountsys
: mountsys sd_init ...

	_sysstart <# # # # # # # # # # # #> C@++ c" sdboot.f" sd_append

	c"  sd_mount ; ~h0D]~h0D" C@++  c" sdboot.f"  sd_append

	c" \ mountusr ( -- ) mount the user disk~h0D[ifndef mountusr~h0D: mountusr sd_init " C@++ c" sdboot.f" sd_append

	_sysstart _sysblocks + <# #s #> C@++ c" sdboot.f" sd_append

	c"  sd_mount ; ~h0D]~h0D"  C@++ c" sdboot.f" sd_append

	c" c~h22 sdboot.f  -  Loading usrboot.f~h7Eh0D~h7Eh0D~h22 .cstr~h0D~h0D" C@++ c" sdboot.f" sd_append

	c" fload usrboot.f~h0D" C@++  c" sdboot.f"  sd_append

	c" c~h22 sdboot.f  -  DONE - PropForth Loaded ~h7Eh0D~h7Eh0D~h22 .cstr~h0D~h0D" C@++ c" sdboot.f" sd_append
	c" hA state andnC!~h0D~h0D" C@++ c" sdboot.f" sd_append

100 fwrite usrboot.f
version W@  .cstr cr

c" usrboot.f  -  initializing~h0D~h0D" .cstr

1 sd_mount

fread .sdcardinfo

c" usrboot.f  -  DONE~h0D~h0D" .cstr
...

mkdir .partion-sys/

_sysstart _sysblocks + sd_mount

mkdir .partion-usr/

_sysstart sd_mount

forget _sysstart

ext: f

max_stars_repo_path: PropForthV5.5-20130317/CurrentRelease/Examples/IMU/GY80/GY80 datalogger and Graph/00-sdfsInitScript.f


text:
C     Exemple de la section 3.1.4 du rapport sur les dependances

C     Version F90 avec motifs. La semantique est un peu differente
C     puisqu'il y a un effet de bord dans f() et que f() est appele
C     moins souvent que dans la version f77.

C     A compiler avec f90

      program cas314_f90

      real X(0:159), Y(0:11)

      do i1 = 0, 39
            X(4*i1:4*i1+3) = f()
      enddo

      print *, X

      do i2 = 0, 11
         Y(i2) = g(X(i2+10:i2+19))
      enddo

c     Valeurs attendues: 72, avec comme increments la suite (+3, +3, +2, +2)*
C     72.0 75.0 78.0 80.0 82.0 85.0 88.0 90.0 92.0 95.0 98.0 100.0
      print *, Y

      end

      function f()
c      static c
      data c /3./
      c = c + 1.
      f = c
      end

      function g(z)
      real z(0:9)

      print *, 'z(0)=', z(0)

      g = 0.
      do i = 0, 9
         g = g + z(i)
      enddo
      end


! { dg-do link }
! PR fortran/40011
!
! Contributed by Joost VandeVondele
!
!
! Before no "one" function was generated with -fwhole-file.
!
!
SUBROUTINE one ( )
END SUBROUTINE one

SUBROUTINE two ( )
END SUBROUTINE two

MODULE mod
CONTAINS
  SUBROUTINE three ( )
    CALL two ( )
  END SUBROUTINE three
  SUBROUTINE four ( )
      CALL one ( )
  END SUBROUTINE four
END MODULE mod
END


text:
      subroutine fdate(str)
C
C#######################################################################
C
C        $Log: lagrit_fdate.f,v $
C        Revision 2.00  2007/11/05 19:46:00  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.21   02 Oct 2007 12:40:28   spchu
CPVCS    original version
C
C#######################################################################
C
      character *(*) str
      character*8 date
      character*10 time
      character*5 zone
      integer values(8)
      call date_and_time(date,time,zone,values)
      str=' '
      str(21:24)= date(1:4)
      str(6:7)=date(5:6)
      str(8:8)='/'
      str(9:10)=date(7:8)
      str(12:13)=time(1:2)
      str(14:14)=':'
      str(15:16)=time(3:4)
      str(17:17)=':'
      str(18:19)=time(5:6)
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text: Rumsey is an unincorporated community 18 miles northwest of Esparto in Yolo County.  Should you want to take a road trip in your automobiles car, youll be taking Highway 113 North to County 27W, then County 93N, and finally Hwy 16, first west, then north.  Youll pass through Esparto, Brooks, and Guinda, seeing farms farmland where some of the Farmers Market produce is grown!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! I think this is using F03 syntax for procedure pointers.
! So this might really be a FFortran 2003 test code
! (syntax fails for gfortran, but is part of the gfortran test suite).

module myMod34

  CONTAINS

  real function proc3( arg1 )
     integer :: arg1
     proc3 = arg1+7
  end function proc3

  subroutine proc4( arg1 )
     procedure(real), pointer :: arg1
     if (arg1(0)/=7) call abort()
  end subroutine proc4

end module myMod34

program myProg
  use myMod34
  PROCEDURE (real), POINTER :: p => NULL()
  p => proc3
  call proc4( p )
end program myProg

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C Copyright (c) 2007-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C=======================================================================
      SUBROUTINE WRNPS (NTXT, NUMNPS, LNPSNL, LNPSDF, IDNPS, NNNPS,
     &                  NDNPS, IXNNPS, IXDNPS, LTNNPS, FACNPS)
C=======================================================================

C   --*** WRNPS *** (EXOTXT) Write database nodal point sets
C   --   Written by Amy Gilkey - revised 02/27/86
C   --   Modified for ExodusIIv2 database format 10/12/95
C   --
C   --WRNPS writes the nodal point set information from the database.
C   --
C   --Parameters:
C   --   LNPSNL - IN - the length of the nodal point sets node list
C   --   IDNPS - IN - the nodal point set ID for each set
C   --   NNNPS - IN - the number of nodes for each set
C   --   IXNNPS - IN - the index of the first node for each set
C   --   LSTNPS - IN - the nodes for all sets
C   --   FACNPS - IN - the distribution factors for all sets

C   --   NTXT   - IN - the text file
C   --   NUMNPS - IN - the number of nodal point sets
C   --   LNPSNL - IN - length of node set node list
C   --   LNPSDF - IN - length of node set distribution factor list
C   --   IDNPS  - IN - array containing the node set ID's for each node set
C   --   NNNPS  - IN - array containing the number of nodes for each node set
C   --   NDNPS  - IN - array containing number of dist. fact for each node set
C   --   IXNNPS - IN - array containing indices into the LTNNPS array which
C                      are the location of the 1st nodes for each set
C   --   IXDNPS - IN - array containing indices into the FACNPS array which
C                      are the location of the 1st dist factor for each set
C   --   LTNNPS - IN - array containing the nodes for all node sets
C                      Internal node IDs
C   --   FACNPS - IN - array containing the distribtion factors
C                      for all node sets

      INTEGER NTXT, NUMNPS, LNPSNL, LNPSDF
      INTEGER IDNPS(*)
      INTEGER NNNPS(*)
      INTEGER NDNPS(*)
      INTEGER IXNNPS(*)
      INTEGER IXDNPS(*)
      INTEGER LTNNPS(*)
      REAL    FACNPS(*)

      WRITE (NTXT, 10005) '! Node sets                  ', NUMNPS
      WRITE (NTXT, 10005) '! Len node set node list     ', LNPSNL
      WRITE (NTXT, 10005) '! Len node set dist fact list', LNPSDF
      DO 100 INPS = 1, NUMNPS
         WRITE (NTXT, '(A, I10)') '! Nodal point set', INPS
         WRITE (NTXT, 10000) IDNPS(INPS), NNNPS(INPS), NDNPS(INPS),
     &      '! ID, nodes, dist factors'

         INS = IXNNPS(INPS)
         INE = IXNNPS(INPS) + NNNPS(INPS) - 1
         if (ins .le. ine)
     &      WRITE (NTXT, 10020) (LTNNPS(NL), FACNPS(NL), NL=INS,INE)

  100 CONTINUE

      RETURN
10000 FORMAT (3I10, 6X, A)
10005 FORMAT (A, I10)
10020 FORMAT (I10, 2X, 1pE16.7)
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
module type_chidg_vector_recv_comm
#include <messenger.h>
    use mod_kinds,          only: ik
    use mod_constants,      only: INTERIOR, CHIMERA, BOUNDARY
    use type_mesh,          only: mesh_t
    use type_ivector,       only: ivector_t
    use type_domain_vector, only: domain_vector_t
    use mod_chidg_mpi,      only: ChiDG_COMM
    use mpi_f08,            only: MPI_Recv, MPI_INTEGER4, MPI_STATUS_IGNORE, MPI_ANY_TAG
    implicit none

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   6/30/2016
    !!
    !!
    !-------------------------------------------------------------------------------
    type, public :: chidg_vector_recv_comm_t

        integer(ik)                         :: proc
        type(domain_vector_t),  allocatable :: dom(:)

    contains

        procedure,  public  :: init
        procedure,  public  :: clear
        procedure,  public  :: ndomains
        procedure,  public  :: restrict
        procedure,  public  :: prolong

    end type chidg_vector_recv_comm_t
    !********************************************************************************

contains

    !>  Receive information from 'proc' about what it is sending, so we can prepare
    !!  here about what to receive.
    !!
    !!  Communication here is initiated from 'type_chidg_vector_send_comm'.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   6/30/2016
    !!
    !!
    !---------------------------------------------------------------------------------
    subroutine init(self,mesh,proc,comm)
        class(chidg_vector_recv_comm_t), intent(inout)   :: self
        type(mesh_t),                    intent(inout)   :: mesh
        integer(ik),                     intent(in)      :: proc
        integer(ik),                     intent(in)      :: comm

        character(:),   allocatable :: user_msg
        integer(ik)                 :: idom, idom_recv, ndom_recv, ierr, ielem, iface,  &
                                       ChiID, donor_domain_g, donor_element_g,          &
                                       idomain_g, ielement_g, dom_store, idom_loop,     &
                                       ielem_loop, ielem_recv, nelem_recv,              &
                                       neighbor_domain_g, neighbor_element_g,           &
                                       bc_domain_g, bc_element_g,                       &
                                       recv_element, recv_domain, idonor,               &
                                       group_ID, patch_ID, face_ID, elem_ID
        integer(ik),    allocatable :: comm_procs_dom(:)
        logical                     :: proc_has_domain, domain_found, element_found,    &
                                       is_interior, is_chimera, is_boundary,            &
                                       comm_neighbor, comm_donor, comm_elem,            &
                                       donor_recv_found, bc_recv_found

        !
        ! Set processor being received from
        !
        self%proc = proc

        !
        ! Get number of domains being received from proc. This is sent from chidg_vector_send_comm
        !
        call MPI_Recv(ndom_recv,1,MPI_INTEGER4,self%proc, MPI_ANY_TAG, ChiDG_COMM, MPI_STATUS_IGNORE, ierr)

        !
        ! Allocate number of recv domains to recv information from proc
        !
        if (allocated(self%dom)) deallocate(self%dom)
        allocate(self%dom(ndom_recv), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Call initialization for each domain to be received
        !
        do idom_recv = 1,ndom_recv
            call self%dom(idom_recv)%init(proc)
        end do ! idom

        !
        ! Loop through mesh and initialize recv indices
        !
        do idom = 1,mesh%ndomains()
            do ielem = 1,mesh%domain(idom)%nelem
                do iface = 1,size(mesh%domain(idom)%faces,2)

                    is_interior = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )
                    is_chimera  = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA  )
                    is_boundary = ( mesh%domain(idom)%faces(ielem,iface)%ftype == BOUNDARY )

                    !
                    ! Initialize recv for comm INTERIOR neighbors
                    !
                    if (is_interior) then

                        comm_neighbor = (proc == mesh%domain(idom)%faces(ielem,iface)%ineighbor_proc)

                        ! If neighbor is being communicated, find it in the recv domains
                        if (comm_neighbor) then
                            neighbor_domain_g  = mesh%domain(idom)%faces(ielem,iface)%ineighbor_domain_g
                            neighbor_element_g = mesh%domain(idom)%faces(ielem,iface)%ineighbor_element_g

                            ! Loop through domains being received to find the right domain
                            do idom_recv = 1,ndom_recv
                                recv_domain = self%dom(idom_recv)%vecs(1)%dparent_g()
                                if (recv_domain == neighbor_domain_g) then

                                    ! Loop through the elements in the recv domain to find the right neighbor element
                                    do ielem_recv = 1,size(self%dom(idom_recv)%vecs)

                                        recv_element = self%dom(idom_recv)%vecs(ielem_recv)%eparent_g()

                                        ! Set the location where a face can find its off-processor neighbor
                                        if (recv_element == neighbor_element_g) then
                                            mesh%domain(idom)%faces(ielem,iface)%recv_comm    = comm
                                            mesh%domain(idom)%faces(ielem,iface)%recv_domain  = idom_recv
                                            mesh%domain(idom)%faces(ielem,iface)%recv_element = ielem_recv
                                        end if

                                    end do !ielem_recv
                                end if
                            end do !idom_recv

                        end if

                    !
                    ! Initialize recv for comm CHIMERA receivers
                    !
                    else if (is_chimera) then

                        ChiID = mesh%domain(idom)%faces(ielem,iface)%ChiID
                        do idonor = 1,mesh%domain(idom)%chimera%recv(ChiID)%ndonors()

                            !comm_donor = (proc == mesh%domain(idom)%chimera%recv(ChiID)%donor_proc%at(idonor) )
                            comm_donor = (proc == mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%iproc )

                            donor_recv_found = .false.
                            if (comm_donor) then
                                !donor_domain_g  = mesh%domain(idom)%chimera%recv(ChiID)%donor_domain_g%at(idonor)
                                !donor_element_g = mesh%domain(idom)%chimera%recv(ChiID)%donor_element_g%at(idonor)
                                donor_domain_g  = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%idomain_g
                                donor_element_g = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%ielement_g

                                ! Loop through domains being received to find the right domain
                                do idom_recv = 1,ndom_recv
                                    recv_domain = self%dom(idom_recv)%vecs(1)%dparent_g()
                                    if (recv_domain == donor_domain_g) then

                                        ! Loop through the elements in the recv domain to find the right neighbor element
                                        do ielem_recv = 1,size(self%dom(idom_recv)%vecs)
                                            recv_element = self%dom(idom_recv)%vecs(ielem_recv)%eparent_g()

                                            ! Set the location where a face can find its off-processor neighbor
                                            if (recv_element == donor_element_g) then
                                                !mesh%domain(idom)%chimera%recv(ChiID)%donor_recv_comm%data_(idonor)    = comm
                                                !mesh%domain(idom)%chimera%recv(ChiID)%donor_recv_domain%data_(idonor)  = idom_recv
                                                !mesh%domain(idom)%chimera%recv(ChiID)%donor_recv_element%data_(idonor) = ielem_recv
                                                mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_comm    = comm
                                                mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_domain  = idom_recv
                                                mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_element = ielem_recv
                                                donor_recv_found = .true.
                                                exit
                                            end if

                                        end do !ielem_recv

                                    end if

                                    if (donor_recv_found) exit
                                end do !idom_recv

                                user_msg = "chidg_vector_recv_comm%init: chimera receiver did not &
                                            find parallel donor."
                                if (.not. donor_recv_found) call chidg_signal_three(FATAL,user_msg,IRANK,donor_domain_g,donor_element_g)
                            end if !comm_donor

                        end do ! idonor

                    !
                    ! Initialize recv for comm BOUNDARY coupling
                    !
                    else if (is_boundary) then

                        !
                        ! Get boundary group/patch/face identifier
                        !
                        group_ID = mesh%domain(idom)%faces(ielem,iface)%group_ID
                        patch_ID = mesh%domain(idom)%faces(ielem,iface)%patch_ID
                        face_ID  = mesh%domain(idom)%faces(ielem,iface)%face_ID

                        !
                        ! Loop through the coupling for this face:
                        !   - detect off-processor coupled elements
                        !   - if off-processor, set recv indices
                        !
                        do elem_ID = 1,mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ncoupled_elements()

                            comm_elem = (proc == mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%proc(elem_ID) )

                            bc_recv_found = .false.
                            if (comm_elem) then
                                ! Get global indices for parallel coupled element
                                bc_domain_g  = mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%idomain_g(elem_ID)
                                bc_element_g = mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ielement_g(elem_ID)

                                ! Loop through domains being received to find the right domain
                                do idom_recv = 1,ndom_recv
                                    recv_domain = self%dom(idom_recv)%vecs(1)%dparent_g()
                                    if (recv_domain == bc_domain_g) then

                                        ! Loop through the elements in the recv domain to find the right neighbor element
                                        do ielem_recv = 1,size(self%dom(idom_recv)%vecs)
                                            recv_element = self%dom(idom_recv)%vecs(ielem_recv)%eparent_g()

                                            ! Set the location where a face can find its off-processor neighbor
                                            if (recv_element == bc_element_g) then
                                                call mesh%bc_patch_group(group_ID)%patch(patch_ID)%set_coupled_element_recv(face_ID,      &
                                                                                bc_domain_g,  &
                                                                                bc_element_g, &
                                                                                comm,         &
                                                                                idom_recv,    &
                                                                                ielem_recv )
                                                bc_recv_found = .true.
                                                exit
                                            end if

                                        end do !ielem_recv

                                    end if

                                    if (bc_recv_found) exit
                                end do !idom_recv

                                user_msg = "chidg_vector_recv_comm%init: boundary face did not &
                                            find parallel element."
                                if (.not. bc_recv_found) call chidg_signal_three(FATAL,user_msg,IRANK,bc_domain_g,bc_element_g)

                            end if !comm_elem
                        end do !elem_ID

                    end if ! is_interior, is_chimera, is_boundary

                end do ! iface
            end do ! ielem
        end do ! idom

    end subroutine init
    !*********************************************************************************

    !>  Zero all entries.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   7/7/2016
    !!
    !!
    !----------------------------------------------------------------------------------
    subroutine clear(self)
        class(chidg_vector_recv_comm_t), intent(inout)   :: self

        integer(ik) :: idom

        do idom = 1,size(self%dom)
            call self%dom(idom)%clear()
        end do

    end subroutine clear
    !**********************************************************************************

    !>  Return number of domains being received from a given processor.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   7/21/2017
    !!
    !----------------------------------------------------------------------------------
    function ndomains(self) result(ndomains_)
        class(chidg_vector_recv_comm_t),    intent(in)  :: self

        integer(ik) :: ndomains_

        if (allocated(self%dom)) then
            ndomains_ = size(self%dom)
        else
            ndomains_ = 0
        end if

    end function ndomains
    !***********************************************************************************

    !>  Return a restricted copy of the chidg_vector_recv_comm_t object.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   7/21/2017
    !!
    !----------------------------------------------------------------------------------
    function restrict(self,nterms_r) result(restricted)
        class(chidg_vector_recv_comm_t),    intent(in)  :: self
        integer(ik),                        intent(in)  :: nterms_r

        integer(ik)                     :: ierr, idom
        type(chidg_vector_recv_comm_t)  :: restricted

        restricted%proc = self%proc

        ! Allocate storage for each recv domain
        allocate(restricted%dom(self%ndomains()), stat=ierr)
        if (ierr /= 0) call AllocationError

        ! For each block return a restricted version
        do idom = 1,self%ndomains()
            restricted%dom(idom) = self%dom(idom)%restrict(nterms_r)
        end do

    end function restrict
    !**********************************************************************************

    !>  Return a prolonged copy of the chidg_vector_recv_comm_t object.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   7/24/2017
    !!
    !----------------------------------------------------------------------------------
    function prolong(self,nterms_p) result(prolonged)
        class(chidg_vector_recv_comm_t),    intent(in)  :: self
        integer(ik),                        intent(in)  :: nterms_p

        integer(ik)                     :: ierr, idom
        type(chidg_vector_recv_comm_t)  :: prolonged

        prolonged%proc = self%proc

        ! Allocate storage for each recv domain
        allocate(prolonged%dom(self%ndomains()), stat=ierr)
        if (ierr /= 0) call AllocationError

        ! For each block return a restricted version
        do idom = 1,self%ndomains()
            prolonged%dom(idom) = self%dom(idom)%prolong(nterms_p)
        end do

    end function prolong
    !**********************************************************************************

end module type_chidg_vector_recv_comm

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
	SUBROUTINE GDPTYAX  ( ptype, yaxis, ymin, ymax, iyaxis,
     +			      ratio, ystrt, ystop, ylbl, nylbl, rmargn,
     +	 		      iylbsf, iyglsf, iytmsf, iret )
C************************************************************************
C* GDPYAX								*
C*									*
C* This subroutine returns the parameters to use for the y axis 	*
C* in GDPTPDF								*
C*									*
C* GDPTYAX ( PTYPE, YAXIS, IVCORD, YMIN, YMAX, IYAXIS, RATIO, YSTRT,  	*
C*          YSTOP, YLBL, NYLBL, RMARGN, IYLBSF, IYGLSF, IYTMSF, IRET )	*
C*									*
C* Input parameters:							*
C*	PTYPE		CHAR*		Y axis type			*
C*	YAXIS		CHAR*		Ymin / ymax / yinc		*
C*	YMIN		REAL		Y min				*
C*	YMAX		REAL		Y MAX				*
C*									*
C* Output parameters:							*
C*	IYAXIS		INTEGER		Y axis integer type		*
C*	RATIO		REAL		Height to width ratio		*
C*	YSTRT		REAL		Bottom y value			*
C*	YSTOP		REAL		Top y value			*
C*	YLBL  (NYLBL)	REAL		Y axis label values		*
C*	NYLBL		INTEGER		Number of y axis labels		*
C*	RMARGN (4)	REAL		Margins				*
C*	IYLBSF		INTEGER		Label frequency			*
C*	IYGLSF		INTEGER		Grid line frequency		*
C*	IYTMSF          INTEGER		Tick mark frequency		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -7 = invalid axis type		*
C*                                      -11 = invalid axis request	*
C**									*
C* Log:									*
C* M. Li/SAIC		08/07						*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	ptype, yaxis
	REAL		ylbl (*), rmargn (*)
C*
	CHARACTER	parm*32
	LOGICAL         skewt
C*
C------------------------------------------------------------------------
	iret = 0
	skewt = .false.
	parm = ' '
C
C*	Get the values from PTYPE.
C
	CALL IN_PTYP ( ptype, iyaxis, ratio, rmargn, ier )

	IF  ( ier .ne. 0 ) THEN
	    iret = -7
	    RETURN
        END IF
C
C*	Get the range of values to use for the y-axis.
C
	ivc    = 0
	ilfdef = 1
        igfdef = 0
        itfdef = 1
	CALL IN_AXIS  ( yaxis, ivc, skewt, parm, ymin, ymax,
     +                  ilfdef, igfdef, itfdef, ystrt, ystop, ylbl,
     +                  nylbl, iylbsf, iyglsf, iytmsf, ier )

C*
	IF ( ier .ne. 0 ) THEN
	   iret = -11
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
French Cuff Consignment is a small chain of thrift stores consignment stores dealing in womens clothing stores clothing. French Cuff specializes in quality fashion items from top designers from around the world. Owned and operated by Darcy McNie and Corey DeRoo, a motherdaughter team. The wiki:sacramento:French Cuff Consignment original location is in Sacramento, and the newest shop is in Davis.

20100817 13:05:02 nbsp French Cuff is now in Davis, too!  I visited Davis for the first time a month ago and the store owner said they had just opened that location.  Their clothes, shoes, purses.... they are fantastic!  Be prepared to spend money because everything there is a find! Users/SteffanyBarrus

20101123 14:18:31 nbsp The Davis location of French Cuff puts all other Davis consignment stores to shame. The selection is amazing  the first time I went in I stayed for over an hour  and clothes are grouped by color/type, not size, so a bit of digging is required to yield treasures. But there are indeed treasures to be found, and some of them are real bargains.

I wish Id known, when I went on my first shopping spree, that as part of their customer recognition your receipt will have a coupon attached for your next purchase. A good strategy: buy what you need the first time and come back with that coupon for more (and there will be more  new stuff arrives all the time).

Most of the time Corey is there; shes warm, friendly and helpful. She really knows her stuff, but is also very chatty. In fact, it can take a while to extricate yourself from a conversation with her  but at least Corey provides ample personal contact with her customers... which is sorely lacking in some of our local stores (Im talkin bout you specially, Nestware)! Users/ZoraDillon

20120110 14:14:53 nbsp Great boutique! Helpful staff, great bargains and gorgeous clothes. Users/Katiehess

20120219 10:35:54 nbsp Love this boutique! The girl who worked there was super friendly and helpful and there was a great variety of cute clothes. The coupon they give you at purchase doesnt hurt either! Users/KelseyCondell

20120430 13:24:22 nbsp While Ill admit, I dont have any experience shopping here, consigning here was a complete waste of time. I brought in 30+ excellent pieces (some designer, some professionally made and some mainstream), clean, pressed and excellent condition. The woman setting up the account hardly spoke to me and started sorting the clothes by stacking them on the counter in a big pile. Once I realized what she was doing, I asked, Are those the ones youre not going to take? She said, Yes, we dont take Gap. To be honest only one or two items were Gap the rest were Ann Taylor, Chicos, etc... I quickly realized that she was only going to take about 4 items so I strolled around to see what other labels they did take for next time. What was the first tag I saw? Gap. The next? Gap. I wanted to walk out, but thought Id give it a chance. After 3 months none of the items she selected sold (two vintage items and two Laundry by Shelli Segal). When I went to pick my items up after my consignment period, she handed all 4 of them to me wrinkled and shoved into a Target bag. Thanks. Users/DebWestergaard


text:
*DK quadxy
      subroutine quadxy(imsgin,xmsgin,cmsgin,msgtype,nwds,ierr2)
C
C#######################################################################
C
C      PURPOSE -
C
C         THIS SUBROUTINE DEFINES AN ARBITRARY, LOCICAL QUAD OF POINTS
C            IN 2D(XY) SPACE.
C
C
C     FORMAT: QUADXY/# points in x/# points in y/                      &
C                    x1,y1,z1/x2,y2,z2/x3,y3,z3/x4,y4,z4/
C                             (counter clockwise)
C                    ixz,iyz/ (cell centered=0, end points =1)
C                    ixratio,iyratio/(regular spacing=1, ratio spacing=0)
C                    xrz,yrz (ratio spacing - xrz is ratio of second
C                              interval to first interval along
C                              line from point 1 to point2
C                              yrz -- ditto along line from point2
C                                     to point3)
C          Note:  the direction from point 1 to point 2 is the 'x' direction
C          Note:  the direction from point 2 to point 3 is the 'y' direction
C      INPUT ARGUMENTS -
C
C        imsgin()  - Integer array of command input tokens
C        xmsgin()  - Real array of command input tokens
C        cmsgin()  - Character array of command input tokens
C        msgtype() - Integer array of command input token types
C        nwds      - Number of command input tokens
C
C
C      OUTPUT ARGUMENTS -
C
C         ierr2 - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C
C      CHANGE HISTORY -
C
C         $Log: quadxy.f,v $
C         Revision 2.00  2007/11/09 20:03:58  spchu
C         Import to CVS
C
CPVCS
CPVCS       Rev 1.6   05 Mar 2002 13:06:22   gable
CPVCS    Changed log file output format from i6 to i9
CPVCS
CPVCS       Rev 1.5   Thu Apr 06 09:28:08 2000   dcg
CPVCS    replace get_info_i and set_info_i calls
CPVCS
CPVCS       Rev 1.4   Mon Apr 14 16:57:04 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.3   Tue Nov 05 09:38:58 1996   jan
CPVCS    modified line 171 - added error flag to arguments
CPVCS
CPVCS       Rev 1.2   Thu Oct 31 14:21:36 1996   dcg
CPVCS    add missing end statement
CPVCS
CPVCS       Rev 1.1   Wed Oct 09 10:25:44 1996   dcg
CPVCS    add ratio spacing to quad routines
CPVCS
CPVCS       Rev 1.0   Mon Jul 29 15:17:30 1996   dcg
CPVCS    Initial revision.
C
C#######################################################################
C
      implicit none
C
      character*132 logmess
C
C#######################################################################
C
      integer nwds, imsgin(nwds), msgtype(nwds),ierr2
      REAL*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
C
      integer ierr,ixratio,iyratio,ixz,iyz,nx,ny,iprint
      integer npoints,ilen,ityp,istart,ipointi,ipointj
C
C#######################################################################
C
      character*32 cmo
      character*8 isubname
C
      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      real*8 xic(1000000), yic(1000000), zic(1000000)
      real*8 x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,xrz,yrz,d
C
C#######################################################################
C
C
      ierr2=0
C
      call cmo_get_name(cmo,ierr2)
      call cmo_get_info('nnodes',cmo,npoints,ilen,ityp,ierr2)
C
      call cmo_get_info('ipointi',cmo,
     *                ipointi,ilen,ityp,ierr)
      if (ierr .ne. 0) call x3d_error(isubname,'cmo_get_info')
      call cmo_get_info('ipointj',cmo,
     *                ipointj,ilen,ityp,ierr)
      if (ierr .ne. 0) call x3d_error(isubname,'cmo_get_info')
C
      istart=2
      nx=imsgin(istart)
         istart=istart+1
      ny=imsgin(istart)
         istart=istart+1
      x1=xmsgin(istart)
         istart=istart+1
      y1=xmsgin(istart)
         istart=istart+1
      z1=xmsgin(istart)
         istart=istart+1
      x2=xmsgin(istart)
         istart=istart+1
      y2=xmsgin(istart)
         istart=istart+1
      z2=xmsgin(istart)
         istart=istart+1
      x3=xmsgin(istart)
         istart=istart+1
      y3=xmsgin(istart)
         istart=istart+1
      z3=xmsgin(istart)
         istart=istart+1
      x4=xmsgin(istart)
         istart=istart+1
      y4=xmsgin(istart)
         istart=istart+1
      z4=xmsgin(istart)
         istart=istart+1
      if(nwds.le.15) then
         ixz=1
         iyz=1
      else
         ixz=imsgin(istart)
            istart=istart+1
         iyz=imsgin(istart)
            istart=istart+1
      endif
      if(nwds.le.17) then
        ixratio=0
        iyratio=0
        xrz=1.
        yrz=1.
      else
         ixratio=imsgin(istart)
            istart=istart+1
         iyratio=imsgin(istart)
            istart=istart+1
         xrz=xmsgin(istart)
            istart=istart+1
         yrz=xmsgin(istart)
            istart=istart+1
      endif

C
C     ******************************************************************
C
C     FILTER THE ZONING FLAGS AND MULTIPLIES TO ALLOW ONLY THE
C        PROPER VALUES.
C
      if(ixz.ne.0.and.ixz.ne.1) ixz=1
      if(iyz.ne.0.and.iyz.ne.1) iyz=1
      if(xrz.le.0.0) xrz=1.0
      if(yrz.le.0.0) yrz=1.0
      iprint=0
      if(ixz.eq.0.and.ixratio.ne.0) then
         iprint=1
         ixratio=0
         xrz=1.
      endif
         if(iyz.eq.0.and.iyratio.ne.0) then
         iprint=1
         ixratio=0
         yrz=1.
      endif
      if(iprint.eq.1) then
         logmess='Ratio spacing not allowed with cell centered points'
         call writloga('default',0,logmess,0,ierr)
      endif
C
C     Adjust end points for cell centered data
C
      if(ixz.eq.0) then
         d=0.5*(x2-x1)/nx
         x1=x1+d
         x2=x2-d
         d=0.5*(y2-y1)/nx
         y1=y1+d
         y2=y2-d
         d=0.5*(z2-z1)/nx
         z1=z1+d
         z2=z2-d
         d=0.5*(x3-x4)/nx
         x4=x4+d
         x3=x3-d
         d=0.5*(y3-y4)/nx
         y4=y4+d
         y3=y3-d
         d=0.5*(z3-z4)/nx
         z4=z4+d
         z3=z3-d
      endif
      if(iyz.eq.0) then
         d=0.5*(x3-x2)/ny
         x2=x2+d
         x3=x3-d
         d=0.5*(y3-y2)/ny
         y2=y2+d
         y3=y3-d
         d=0.5*(z3-z2)/ny
         z2=z2+d
         z3=z3-d
         d=0.5*(x4-x1)/ny
         x4=x4-d
         x1=x1+d
         d=0.5*(y4-y1)/ny
         y4=y4-d
         y1=y1+d
         d=0.5*(z4-z1)/ny
         z4=z4-d
         z1=z1+d
      endif
C
C     Set nnodes in mesh object and get pointers
C
      call cmo_set_info('nnodes',cmo,npoints+nx*ny,1,1,ierr2)
      call cmo_newlen(cmo,ierr2)
      call cmo_get_info('xic',cmo,ipxic,ilen,ityp,ierr2)
      call cmo_get_info('yic',cmo,ipyic,ilen,ityp,ierr2)
      call cmo_get_info('zic',cmo,ipzic,ilen,ityp,ierr2)
C
C     CREATE the points.
C
      call quadpts(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,nx,ny,xrz,yrz,
     x             xic(npoints+1),yic(npoints+1),zic(npoints+1))

C
C     ******************************************************************
C     PRINT OUT THE POINT NUMBERS GENERATED
C
      write(logmess,6000) npoints+1,npoints+nx*ny
 6000 format('  QUADXY GENERATED POINTS ',i9,' TO ',i9)
      call writloga('default',0,logmess,0,ierr)
C
      ipointi=npoints+1
      ipointj=npoints+nx*ny
      call cmo_get_name(cmo,ierr2)
C
      call set_info_i('ipointi',cmo,
     *                ipointi,1,1,ierr)
      if (ierr .ne. 0) call x3d_error(isubname,'cmo_get_info')
      call set_info_i('ipointj',cmo,
     *                ipointj,1,1,ierr)
      if (ierr .ne. 0) call x3d_error(isubname,'ccmo_get_info')
C
      goto 9999
 9999 continue
      return
      end
      subroutine getdelta (xi,yi,zi,xf,yf,zf,n,r,deltax,
     x                     deltay,deltaz)
C   ******************************************************
C
C   PURPOSE:
C     return proportional spacing delta
C
C  INPUT:
C     xi,yi,zi  coordinates of first point
C     xf,yf,zf  coordinates of last point
C     n         number of interval
C     r         ratio of i+1 interval distance to i interval
C
C  OUTPUT:
C     deltax,deltay,deltaz incremental distances to add to
C                          first point to get next point
C
C   *******************************************************
      implicit none
      real*8 xi,yi,zi,xf,yf,zf,r,deltax,deltay,deltaz
      integer n,i
      real*8 sum,xd,yd,zd
C
      xd=xf-xi
      yd=yf-yi
      zd=zf-zi
      if (n.le.1) then
         deltax=xd
         deltay=yd
         deltaz=zd
      else
         sum=0.0
         do i=1,n-1
            sum=sum+r**i
         enddo
         sum=sum+1.0
         deltax=xd/sum
         deltay=yd/sum
         deltaz=zd/sum
      endif
      return
      end
      subroutine quadpts (x1,y1,z1,x2,y2,z2,
     x                    x3,y3,z3,x4,y4,z4,
     x                    nx,ny,rx,ry,x,y,z)
C   ******************************************************
C
C   PURPOSE:
C     create nx by ny set of points using ratio spacing.
C
C  INPUT:
C     x1,y1,z1  coordinates of first point
C     x2,y2,z2  coordinates of second point (counterclockwise)
C     x3,y3,z3  coordinates of third point
C     x4,y4,z4  coordinates of last point
C     nx        number of points in x direction
C     ny        number of points in y direction
C     rx        ratio of i+1 interval distance to i interval
C                x-direction
C     rx        ratio of i+1 interval distance to i interval
C                y-direction
C
C  OUTPUT:
C     x,y,z	    coordinates of created points
C
C   ******************************************************
      implicit none
      real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,rx,ry,
     x       x(*),y(*),z(*),
     x       deltax,deltay,deltaz,factor
      integer nx,ny,i,j,icount,icount1,icount2
C
C   get delta spacing for first nx points
C
      x(1)=x1
      y(1)=y1
      z(1)=z1
      icount=1
      if(nx.ge.2) then
         call getdelta(x1,y1,z1,x2,y2,z2,nx-1,rx,
     *              deltax,deltay,deltaz)
         factor=1.0
         do i=2,nx-1
            if(i.ge.3) factor=factor*rx
            icount=icount+1
            x(icount)=x(icount-1)+deltax*factor
            y(icount)=y(icount-1)+deltay*factor
            z(icount)=z(icount-1)+deltaz*factor
         enddo
         icount=icount+1
         x(icount)=x2
         y(icount)=y2
         z(icount)=z2
      endif
C
C   get delta spacing for last nx points
C
      if(ny.le.1) go to 9999
      if(nx.ge.2) then
         call getdelta(x4,y4,z4,x3,y3,z3,nx-1,rx,
     *              deltax,deltay,deltaz)
      else
         deltax=0.
         deltay=0.
         deltaz=0.
      endif
      icount=nx*(ny-1)+1
      x(icount)=x4
      y(icount)=y4
      z(icount)=z4
      if(nx.ge.2) then
         factor=1.0
         do i=2,nx-1
            if(i.ge.3) factor=factor*rx
            icount=icount+1
            x(icount)=x(icount-1)+deltax*factor
            y(icount)=y(icount-1)+deltay*factor
            z(icount)=z(icount-1)+deltaz*factor
         enddo
         icount=icount+1
         x(icount)=x3
         y(icount)=y3
         z(icount)=z3
      endif
C
C   fill in middle row by using getdelta with y ratio to get
C   end points of rows -- then use getdelta with x ratio
C   for x spacing.
C
      if(ny.gt.1) then
         call getdelta(x1,y1,z1,x4,y4,z4,ny-1,ry,
     *              deltax,deltay,deltaz)
         factor=1.0
         do i=2,ny-1
            if(i.ge.3) factor=factor*ry
            icount=nx*(i-1)+1
            x(icount)=x(icount-nx)+deltax*factor
            y(icount)=y(icount-nx)+deltay*factor
            z(icount)=z(icount-nx)+deltaz*factor
         enddo
         call getdelta(x2,y2,z2,x3,y3,z3,ny-1,ry,
     *              deltax,deltay,deltaz)
         factor=1.0
         do i=2,ny-1
            if(i.ge.3) factor=factor*ry
            icount=nx*i
            x(icount)=x(icount-nx)+deltax*factor
            y(icount)=y(icount-nx)+deltay*factor
            z(icount)=z(icount-nx)+deltaz*factor
         enddo
         do i=2,ny-1
            icount1=nx*(i-1)+1
            icount2=nx*i
            call getdelta(x(icount1),y(icount1),z(icount1),
     *              x(icount2),y(icount2),z(icount2),
     *              nx-1,rx,deltax,deltay,deltaz)
            factor=1.0
            do j=2,nx-1
               if(j.ge.3) factor=factor*rx
               icount1=icount1+1
               x(icount1)=x(icount1-1)+deltax*factor
               y(icount1)=y(icount1-1)+deltay*factor
               z(icount1)=z(icount1-1)+deltaz*factor
            enddo
         enddo
      endif
 9999 return
      end