text:
c *****getKofEMatrix******************
c    fills rate matricies

      subroutine getKofEMatrix(E,dE,rKofEI,rKofEP)
      implicit none

c    get header files
      include 'cdparams.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdwell2.fh'
      include 'cdisprop.fh'
      include 'cdcontrl.fh'

c   local variables
      integer iWell,ioWell,iProd
      real*8 E,dE,rhoDenom
      real*8 rKofEI(mxWells,mxWells),rKofEP(mxWells,mxProds)
      real*8 calRho, xNOS
c    end declaratons

c    note we kludge it if the ratio goes awry

      do 100 iWell = 1,nWells

c       make sure bottom of well is above E
         if (E.gt.EWell(iWell)) then

c dmm 20000404 replace calRho with xNOS
c            rhoDenom = calRho(E-EWell(iWell),dE,iWell,0)
c            rhoDenom = xNOS(E-Ewell(iWell), dE, iWell)

c tal 20000621 switch
             if (oldrho) then
                rhoDenom = calRho(E-EWell(iWell),dE,iWell,0)
             else
                rhoDenom = xNOS(E-Ewell(iWell), dE, iWell)
             end if
c tal 20000621

c          first, do Product channels
            do 60 iProd = 1,nProds(iWell)

c tal 20000621
             if (oldrho) then
               rKofEP(iWell,iProd) =
     2            dmin1(calRho(E-EWell(iWell)-EProd(iWell,iProd),dE,
     3               iWell,nWells+iProd) / rhoDenom, 1.d0)
             else
               rKofEP(iWell,iProd) =
     2            dmin1(xNOS(E-EWell(iWell)-EProd(iWell,iProd),dE,
     3               iWell) / rhoDenom, 1.d0)
             end if
c tal 20000621

60           continue

c          second, do isomerization channels
            do 80 ioWell = 1, nWells
               if (Aisom(iWell,ioWell).gt.1.d-175) then

c tal 20000621
                 if (oldrho) then
                  rKofEI(iWell,ioWell) =
     2               dmin1(calRho(E-EWell(iWell)-EIsom(iWell,ioWell),dE,
     3                  iWell,ioWell) / rhoDenom, 1.d0)
                 else
                  rKofEI(iWell,ioWell) =
     2               dmin1(xNOS(E-EWell(iWell)-EIsom(iWell,ioWell),dE,
     3                  iWell) / rhoDenom, 1.d0)
                 end if
               else
                  rKofEI(iWell,ioWell) = 0.d0
               endif
c temporary k(E) printing line (commented out by pey, 25/3/04
c               write(99,81) iWell, ioWell, E, EWell(iWell), Eisom(iWell
c     $              ,ioWell), rKofEI(iWell,ioWell)
80          continue
 81         format(I2,1x,I2,1x,F7.2,1x, F7.2, 1x, F7.2, 1x, E13.4)
c     dmm 20000404 end replacement of calRho with xNOS

c      zero out for well whose depths are above E
         else
            do 90 iProd = 1,nProds(iwell)
               rkofEP(iWell,iProd) = 0.d0
90          continue

            do 95 ioWell = 1,nWells
               rkofEI(iWell,ioWell) = 0.d0
95          continue
         endif

100   continue

      return
      end

c ******FILLAR1******************
c    fills rate vectors

      subroutine fillAr1(it,rKofEI,RKofEP,RateI,RateP,ACDepI,ACDepP)
      implicit none

c    get header files
      include 'cdparams.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdwell2.fh'
c     include 'cdisprop.fh'
      include 'cdrange.fh'

c   local variables
      integer iWell,ioWell,iProd,it
      real*8 rateP(mxWells,0:mxProds),rateI(mxWells,mxWells)
      real*8 rKofEI(mxWells,mxWells),rKofEP(mxWells,mxProds),
     2   ACDepI(mxWells),ACDepP(mxWells)

c    end declaratons

      do 100 iWell = 1,nWells
         ACdepI(iWell) = 0.d0
         ACdepP(iWell) = 0.d0

c       first, do Product channels
         do 60 iProd = 1,nProds(iWell)
            RateP(iWell,iProd) = Aprod(iWell,iProd)*
     2         temp(it)**rNProd(iWell,iProd)*
     3         dexp(-alphaProd(iWell,iProd)*temp(it))*
     4         rKofEP(iWell,iProd)

            ACdepP(iWell) = ACdepP(iWell) +
     2        RateP(iWell,iProd)
60        continue

c       second, do isomerization channels
         do 80 ioWell = 1, nWells
            if (Aisom(iWell,ioWell).gt.1.d-175) then
               RateI(iWell,ioWell) = AIsom(iWell,ioWell)*
     2            temp(it)**rNIsom(iWell,ioWell)*
     3            dexp(-alphaIsom(iWell,ioWell)*temp(it))*
     4            rKofEI(iWell,ioWell)

               ACdepI(iWell) = ACdepI(iWell) +
     2            RateI(iWell,ioWell)
            else
               RateI(iWell,ioWell) = 0.d0
            endif
80       continue
100   continue

      return
      end

c ******FILLAR2******************
c    fills ACarray

      subroutine fillAr2(it,ip,ACDepP,ACDepI,stab,ACarray,
     2   RateI,RateP,noAcces,noExit,Rleak)

      implicit none

c    get header files
      include 'cdparams.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdwell2.fh'
      include 'cdisprop.fh'
      include 'cdrange.fh'

c   local variables
      logical noAcces(mxWells),noexit
      integer iWell,ioWell,it,ip
      real*8 ACdepI(mxWells),ACdepP(mxWells),rateP(mxWells,0:mxProds),
     2   rateI(mxWells,mxWells),ACarray(mxWells,mxWells),
     3   stab(mxWells,mxTpts),Rleak
c    end declaratons

c    we provide some leakage current of M for the low P limit
c    hopefully, this will help avoid singularities; compare this
c    to typical M's:  P(atm) /(80*T)
      Rleak = 1.d-70
      noexit = .false.

c    first case (ip=0) is low p limit
      if (ip.eq.0) then
         do 130 iWell = 1,nWells
c          first write input channels
            do 120 ioWell = 1,nWells
               ACarray(iWell,ioWell) = RateI(ioWell,iWell)
120         continue

c          then write depletion channel but
c          don't count stabilization rate for low p limit
            RateP(iWell,0) = stab(iWell,it)
            ACarray(iWell,iWell) =  -ACdepP(iWell) -ACdepI(iWell)
     2          -Rleak*stab(iWell,it)
130      continue

c    second case is for general p
      else if (ip.le.npres) then
         do 160 iWell = 1,nWells
c          first write input channels
            do 150 ioWell = 1,nWells
                ACarray(iWell,ioWell) = RateI(ioWell,iWell)
150         continue

c          then write depletion channel
            RateP(iWell,0) = stab(iWell,it) * pres(it,ip) /
     2         (82.1d0*temp(it))
            ACarray(iWell,iWell) =  -ACdepP(iWell) -ACdepI(iWell)
     2         - RateP(iWell,0)
160      continue

c       third case (ip=npres+1) is h p limit, here only count
c       input channels from Wells of one less depth
         else
            do 190 iWell = 1,nWells
               do 180 ioWell = 1,nWells
                   if (idpWell(ioWell).eq.idpWell(iWell)-1) then
                      ACarray(iWell,ioWell) = RateI(ioWell,iWell)
                   else
                      ACarray(iWell,ioWell) = 0.d0
                   endif
180             continue

c           in high P limit deactivation channel domimates
             RateP(iWell,0) = stab(iWell,it)
             ACarray(iWell,iWell) = -RateP(iWell,0)
190       continue
       endif

c     NEW!!! since we can make arbitrary linear combination of
c     our equations, we change the equation corresponding to inpWell
c     to be the sum of all equations - the right hand side does not
c     change because all other entries are zero; on the matrix side,
c     all isomerization rates must cancel, leaving only the stabilization
c     and Product channels.  we put these in explicitly.
c     this should make things better behaved at low pressures

c IMPORTANT:  we put an amazing kludge to make the low pressure rate
c             come out in the low P case for dissoc - if there are
c             no exit channels the stabilization rates are the only
c             things across the inpWell line - in this case we can't
c             put in zero but we put in a small number times stab
c             and flag it;
c             above, we must cancel it

       if (ip.eq.0) then
          noexit = .true.
          do 200 iWell = 1,nWells
             if ((ACdepP(iWell).gt.1.d-175).and.(.not.noAcces(iWell)))
     2          noexit = .false.
200       continue

          if (noexit) then
             do 210 iWell = 1,nWells
                ACarray(inpWell,iWell) = -Rleak*stab(iWell,it)
210          continue
          else
             do 220 iWell = 1,nWells
                ACarray(inpWell,iWell) =  -ACdepP(iWell)
     2             -Rleak*stab(iWell,it)
220          continue
          endif

       else if (ip.le.npres) then
          do 250 iWell = 1,nWells
             ACarray(inpWell,iWell) =  -ACdepP(iWell) -RateP(iWell,0)
250       continue
      endif

c    we leave hp limit alone - this matrix should be okay by itself,
c    plus it doesn't satisfy the same summation condition anyway

      return
      end

c***SOLVE*******************
      subroutine solve(N,A,B,id,Nd,X,elim,ierr,lecho, E)
c    id corresponds to idpWell; we need this for elimination procedure
c    Nd is largest value of id
c    IMPORTANT: parameter mxWells must be identical to that of main

      implicit none
      include 'cdparams.fh'

      real rNegThresh
      parameter (rNegThresh = 1.d-10)

c    local variables
      integer i,j
      integer N,Nsubmit,index,jndex,id(mxWells),Nd,ierr,lecho,lout
      logical elim(mxWells),neg
      real*8 A(mxWells,mxWells),Asubmit(mxWells,mxWells),
     2   B(mxWells),X(mxWells),Bsubmit(mxWells),Xout(mxWells), E
c    end local variables

      neg = .false.

c temp debug line
c      write(*,*) 'mxWells in solve', mxWells

c    we use input vector elim to get rid of inactive channels
c    note this could still leave some singularities near threshold
c    in the low p limit where there may not be any exit channels
c    to an isomer

c    now eliminate rows and columns corresponding to inactive channels
      index = 0
      do 200 i = 1,N
         if (.not.elim(i)) then
            index = index+1
             Bsubmit(index) = B(i)

             jndex = 0
             do 190 j = 1,N
               if (.not.elim(j)) then
                  jndex = jndex+1
                  Asubmit(index,jndex) = A(i,j)
               endif
190         continue

         endif
200   continue
      Nsubmit = index

c    now get solution to Asubmit * Xout = Bsubmit

      if (Nsubmit.gt.0)
     2   call getB(Nsubmit,Asubmit,Bsubmit,Xout,ierr,lout, E)

c    reconstruct solution, note X's have already been zeroed for default
c    NOTE:  now - we try something - if any of the X's come out really
c    negative (more than rNegThresh: this is unphysical and we zero them
c    all out; else we just zero out that channel

      index = 0
      do 300 i=1,N
         if (.not.elim(i)) then
            index = index+1
	    if (Xout(index).lt.-rNegThresh) neg = .true.
            X(i) = Xout(index)
            if (X(i).lt.0.d0) X(i) = 0.d0
         endif
300   continue

      if (neg) then
         do 350 i = 1,N
            X(i) = 0.d0
350      continue
      endif

      return
      end

c****GETB1**********************************************
c    here, we attempt to link to gauss elim routines
c    from Numerical recipes, by press et al.
      subroutine getB(N,A,B,X,ierr,lout, E)

      implicit none
      include 'cdparams.fh'
      include 'cdwell2.fh'
      include 'cdwell0.fh'

c    local variables
      integer i
      integer N,ierr,lout
      real*8 A(mxWells,mxWells),B(mxWells),X(mxWells), E
      character*8 name
c    end local variables

c***********************************************************************
c     20000809 tlada code for printing out A arrays ... commented out
c     for now ...
c      name = char(int((E - Exmin(inpwell))/10)+48) //
c     2       char(int(mod((E-Exmin(inpwell)),10)) + 48)
c      open (78, file=name, status='unknown')

c      write(78,10) E
c 10   format('E = ', D10)
c      if (N .eq. 1) then
c         write(78,25) A(1,1)
c      else if (N .eq. 2) then
c         write(78,26) A(1,1),A(1,2),A(2,1),A(2,2)
c      else if (N.eq. 3) then
c         write (78,27) A(1,1),A(1,2),A(1,3),A(2,1),A(2,2),A(2,3),
c     2        A(3,1),A(3,2),A(3,3)
c      else
c         write (78,28) A(1,1),A(1,2),A(1,3),A(1,4),A(2,1),A(2,2),
c     2        A(2,3),A(2,4),A(3,1),A(3,2),A(3,3),A(3,4),A(4,1),A(4,2),
c     3        A(4,3),A(4,4)
c      end if
c 25   format('A = [', D10, ']')
c 26   format('A = [', D10, D10, ' ; ', D10, D10, ']')
c 27   format('A = [',D10,D10,D10,';',D10,D10,D10,';',D10,D10,D10,']')
c 28   format('A = [',D10,d10,d10,d10,';',d10,d10,d10,d10,';',
c     2     d10,d10,d10,d10,';',d10,d10,d10,d10,']')
c***********************************************************************

      ierr = 0

c    get gaussian elimin solver
      call gaussj(A,N,mxWells,B,1,1)

c    return X's
      do 300 i=1,N
         X(i) = B(i)
300   continue

c***********************************************************************
c     20000809 tlada code for printing matrices ...
c      if (N.eq.1) then
c         write(78,305) B(1)
c      else if (N.eq.2) then
c         write(78,306) B(1),B(2)
c         else if (N.eq.3) then
c            write(78,307) B(1),B(2),B(3)
c            else
c               write(78,308) B(1),B(2),B(3), B(4)
c      end if
c 305  format('B = [',D10,']')
c 306  format('B = [',D10,';',D10,']')
c 307  format('B = [',D10,';',D10,';',D10,']')
c 308  format('B= [',D10,';',D10,';',D10,';',D10,']')
c      close(78)
c***********************************************************************

      return
      end

c****GETB2**********************************************
c    here we attempt link to ldu routines
c    from NUMERICAL RECIPES by Press et al.
      subroutine getB2(N,A,B,X,ierr,lout)

      implicit none
      include 'cdparams.fh'

c    local variables
      integer i,j
      integer N,ierr,lout,Ndim,Idum(mxWells)
      real*8 A(mxWells,mxWells),Adum(mxWells,mxWells),Bdum(mxWells),
     2   B(mxWells),X(mxWells),Vdum(mxWells),D
c    end local variables

c    we use dummy arrays to leave originals uncorrupted
c    we initialize D and Vdum just to surpress warnings

      Ndim = mxWells
      D = 0.d0
      ierr = 0

c    initialize dummy arrays for input matrices
      do 260 i = 1,N
	 Vdum(i) = 0.d0
         Bdum(i) = B(i)
         do 260 j = 1,N
            Adum(i,j) = A(i,j)
260   continue

c    get LU decomposed Adum
      call ludcmp(Adum,N,Ndim,Idum,D,Vdum,ierr,lout)

c    proceed if no errors encountered
      if(ierr.eq.0) then
c       get solution for X
         call lubksb(Adum,N,Ndim,Idum,Bdum)

c      reconstruct solution
         do 300 i = 1,N
            X(i) = Bdum(i)
300      continue
      endif

      return
      end

c***REMOVED GETB3 for successful linking on SGI <ayc 1/94>

c***CLEANUP**********************
      subroutine cleanup(option)
      implicit none

c    get header files
      include 'cdparams.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'
      include 'cdwell2.fh'
      include 'cdrates.fh'

c    local variables
      character*8 option
      integer iWell,iProd
c    end declarations

c       get integer powers for lindemann forms
c       (depends on depth of Well and whether
c       stabilization or Product channel)

      if (option.eq.'Chemact') then

         do 500 iWell = 1,nWells
            do 500 iProd = 0,nProds(iWell)

               nPr(iWell,iProd) = idpWell(iWell)

               if (iProd.eq.0) then
                  nRKinf(iWell,iProd) = 1 - idpWell(iWell)
               else
                  nRKinf(iWell,iProd) = -idpWell(iWell)
               endif
500      continue

      else
c       all dissoc channels

         do 600 iWell = 1,nWells
            do 600 iProd = 0,nProds(iWell)

               nPr(iWell,iProd) = idpWell(iWell)

               if (iProd.eq.0) then
                  nRKinf(iWell,iProd) = 2 - idpWell(iWell)
               else
                  nRKinf(iWell,iProd) = 1 - idpWell(iWell)
               endif
600      continue

      endif

      return
      end

c***********************************************************************
c***********************************************************************
c dmm 20000908
c     initCDcalc
c
c     After debugging we have found that certain arrays need to be 0'd
c     out at the beginning of each CHEMDIS run.  This function does that
c     .  Called from cddiss.f, cdchem.f
c
c***********************************************************************
c***********************************************************************

      subroutine initCDcalc(ACpop, RateP)

      implicit none
      include 'cdparams.fh'
      include 'cdrates.fh'

      real*8 ACpop(mxWells), RateP(mxWells, 0:mxProds)
      integer iWell, iProd, it, ip

      do 10 iWell = 1, mxWells
         ACpop(iWell) = 0.0D0
         do 20 iProd = 0, mxProds
            RateP(iWell, iProd) = 0.0D0
c$$$            do 30 it = 1, mxTPts
c$$$               do 40 ip = -1, mxPPts+1
c$$$                  RK(iWell, iProd, it, ip) = 0.0D0
c$$$ 40            continue
c$$$ 30         continue
 20      continue
 10   continue

      return
      end

c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** START                                     c
c----------------------------------------------------------------------c
c***********************************************************************
c***********************************************************************
c     dmm 20011022
c
c     kofERRKM
c
c     This function uses the sum of states stored in xrhoofn as 'well'
c     iTS and the sum of states for iWell to compute a simple version
c     of the RRKM k(E) for printing to the Barker file.
c
c     Function evalues the density of states for the denominator as the
c     average of DOS over Evr + Ebar to Evr + Ebar + Egrain.  Sum of
c     states evaluated for the arithmetic middle of the grain ...
c***********************************************************************
c***********************************************************************

      subroutine kofERRKM(Evr, Ebar, Egrain, iWell, iTS, xkofE, xkrev)

      implicit none

      include 'cdparams.fh'
      include 'straightBS.fh'
      include 'cdnos.fh'
      include 'cdwell0.fh'
      include 'cdwell1.fh'

      real*8 xNOScm, xkrev, rhode2, Ebar2
      real*8 Evr, Ebar, Egrain, xkofE, hplanc, pdgn, sumTS, rhodenom
      integer iWell, iTS, nETS

c     Degeneracy is needed here, since this would normally be embedded
c     in the A-factor for QRRK but is never available explicitly

      pfflag = .true.
      pdgn = dfloat(idgRRKM)
      hplanc = (6.6260755D-34/(1000.0*4.18))*349.755*6.02D+23
      rhodenom = xNOScm( (Evr + Ebar),Egrain,i1RRKM)/Egrain
c     dmm 20011102  added +1 offset to nETS
      nETS = nint((Evr + 0.5*Egrain)/bsgs) + 1
      sumTS = xrhoofn(nETS, nWells+1)
      if(sumTS.le.0.0D0) then
         sumTS = 1.0D0
      endif
      xkofE = pdgn*sumTS/(hplanc*rhodenom)
      pfflag = .false.

c***********************************************************************
c detailed balance k(E) calculation
c***********************************************************************

      Ebar2 = Eisom(i2RRKM,i1RRKM)*349.755
      rhode2 = xNOScm( (Evr + Ebar2), Egrain, i2RRKM)/Egrain

      xkrev = xkofE*rhodenom/rhode2
c      write(91,100) Evr, xkrev
 100   format(1x, F9.2, 1x, E16.9)
      return
      end
c----------------------------------------------------------------------c
c *** RRKM 1 k(E) Thread *** END                                       c
c----------------------------------------------------------------------c

c $Id$
c $Author$
c $Date$
c $Source$
c $Revision$
c $Log$
c Revision 1.1  2007-02-20 23:10:23  sandeeps
c Initial revision
c
c Revision 1.8  2001/11/21 03:16:45  dmmatheu
c after RRKM changes (successful?); before adding reverse rate constant feature
c
c Revision 1.7  2001/10/22 17:55:44  dmmatheu
c before adding RRKM k(E) function
c
c Revision 1.6  2001/09/20 18:03:19  dmmatheu
c before modifying initCDcalc to see whether we can do without
c time-consuming zeroing of the very large RK matrix
c
c Revision 1.5  2000/08/09 16:57:29  dmmatheu
c tlada changes -- matrix printing code (for matlab-based solver
c debugging) is incomplete and here commented out (dmm)
c
c Revision 1.4  2000/06/23 15:16:51  dmmatheu
c test of VC after merging ...
c
c Revision 1.3  2000/06/23 14:59:39  dmmatheu
c before merging with Tom's code ... working debugged version
c
c Revision 1.2  2000/04/04 17:25:02  dmmatheu
c before replacing calRho with xNOS
c
c Revision 1.1  2000/04/04 17:24:43  dmmatheu
c Initial revision
c

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


MODULE test_module
!  Purpose:
!    To illustrate the use of allocatable arguments.
!
CONTAINS

   SUBROUTINE test_alloc(array)
   IMPLICIT NONE
   REAL,DIMENSION(:),ALLOCATABLE,INTENT(INOUT) :: array
                                     ! Test array

   ! Local variables
   INTEGER :: i             ! Loop index
   INTEGER :: istat         ! Allocate status

   ! Get the status of this array
   IF ( ALLOCATED(array) ) THEN
      WRITE (*,'(A)') 'Sub: the array is allocated'
      WRITE (*,'(A,6F4.1)') 'Sub: Array on entry = ', array
   ELSE
      WRITE (*,*) 'In sub: the array is not allocated'
   END IF

   ! Deallocate the array
   IF ( ALLOCATED(array) ) THEN
      DEALLOCATE( array, STAT=istat )
   END IF

   ! Reallocate as a 5 element vector
   ALLOCATE(array(5), STAT=istat )

   ! Save data
   DO i = 1, 5
       array(i) = 6 - i
   END DO

   ! Display contents of array a on exit
   WRITE (*,'(A,6F4.1)') 'Sub: Array on exit = ', array

   ! Return to caller
   END SUBROUTINE test_alloc

END MODULE test_module

PROGRAM allocatable_arguments
!
!  Purpose:
!    To illustrate the use of allocatable arguments.
!
USE test_module
IMPLICIT NONE

! Declare local variables
REAL,ALLOCATABLE,DIMENSION(:) :: a
INTEGER :: istat         ! Allocate status

! Allocate array
ALLOCATE( a(6), STAT=istat )

! Initialize array
a = (/ 1., 2., 3., 4., 5., 6. /)

! Display a before call
WRITE (*,'(A,6F4.1)') 'Main: Array a before call = ', a

! Call subroutine
CALL test_alloc(a)

! Display a after call
WRITE (*,'(A,6F4.1)') 'Main: Array a after call = ', a

END PROGRAM allocatable_arguments


text:
c     This is the API for foo subroutine. We take in some inputs, apply some cool things and return other things.
      subroutine foo
        integer*2 a ! This is my variable 'a' ... it stores numbers both larger and small
        CHARACTER*2 some_string(15)

        some_string = "some_string"

        if (foo) then
            IF (some_string(0).eq.'s') foo=9                              Text after col 72 are comments and should remain
        endif
      end

ext: f

max_stars_repo_name: envp/fortran-src

max_issues_repo_name: envp/fortran-src

max_forks_repo_name: envp/fortran-src


text:
c***INITX************************************
c    sub to return initial-guess solution and scaling parameters
      subroutine initX(Npar,X0,Xscale,Fscale,IDstr,option)
      implicit none

c    local variables
      integer Npar
      real*8 X0(*),Xscale(*),Fscale
      character IDstr*(*),option*8
c    end declarations

c    IMPORTANT: these options are set in getparams

      if (option.eq.'Troe') then
         call initTroe(Npar,X0,Xscale,Fscale,IDstr)

      else if (option.eq.'SRI') then
         call initSRI(Npar,X0,Xscale,Fscale,IDstr)

      else
         call initDEF(Npar,X0,Xscale,Fscale,IDstr)

      endif

      return
      end

c***GETFFIT************************************
c    routine to return fitted F
      function getFfit(Npar,X,T,P,option)
      implicit none

c    local variables
      integer Npar
      real*8 X(*),T,P
      real*8 getFfit,getTroe,getSRI,getDef
      character option*8
c    end declarations

      if (option.eq.'Troe') then
         getFfit = getTroe(Npar,X,T,P)

      else if (option.eq.'SRI') then
         getFfit = getSRI(Npar,X,T,P)

      else
         getFfit = getDef(Npar,X,T,P)

      endif

      return
      end

c***INITDEF************************************
c    companion routine to getDef for initializing variables
      subroutine initDef(Npar,X0,Xscale,Fscale,IDstr)
      implicit none

c    local variables
      integer iter
      integer Npar,iamp,ioff,iwidth,iasym
      real*8 X0(*),Xscale(*),Fscale
      character IDstr*(*)
c    end declarations

c    make consistent with getDef
      data iamp,iwidth,iasym,ioff/1,6,9,12/
      Npar = 12

c    no better idea so set Xscale to 1 - see IMSL description;
c    also zero out all X0's before setting non-zero values
c    ok - let's try not quite zeroing out (we want init vector
c    all the same order of magnitude)

c    enter FIT function ID string here
      IDstr = 'G (12): 5 amp, 3 width, 3 asym, 1 shift, exp = 1/2'

      do 10 iter = 1,Npar
         Xscale(iter) = 1.d0
         X0(iter) = 1.d-1
10    continue
      Fscale = 1.d0

c    Now we reset the amplitude and width
      X0(iamp) = 1.d0
      X0(iwidth) = 2.d0

      return
      end

c***GETDEF************************************
c    function to return defaulf F fitting function
      function getDef(Npar,X,T,P)
      implicit none

c    local variables
      integer Npar,iamp,ioff,iwidth,iasym
      real*8 X(*),T,P,FcentLn,FLShape,Plgoff,Width,asym,Plg,
     2   arg,Tnm
      real*8 getDef
c    end declarations

      data iamp,iwidth,iasym,ioff/1,6,9,12/

c    IMPORTANT: data statement must be consistent with initDef
c               IDstr set in initDef

c    we invent new scaled T
      Tnm = (1.d-3*T)**0.5

c    use base 10 log for P
      Plg = dlog10(P)

c    a: we use first 5 parameters for amplitude function
c    and we try fitting to cubic in 1000./T
      FcentLn = X(iamp) + X(iamp+1)/(Tnm) + X(iamp+2)/(Tnm*Tnm)
     2   + X(iamp+3)/(Tnm*Tnm*Tnm) + X(iamp+4)/(Tnm*Tnm*Tnm*Tnm)

c    the rest are lineshape parameters:
c    b: we use 3 parameters for Half-width
      Width = X(iwidth) + X(iwidth+1)/(Tnm) + X(iwidth+2)/(Tnm*Tnm)

c    c: we use 3 parameters for asymmetry parameter
      asym = X(iasym) + X(iasym+1)/(Tnm) + X(iasym+2)/(Tnm*Tnm)

c    d: we use 1 parameter for offset from 0
      Plgoff =  X(ioff)

c    now the lineshape function

c    Lorentzian: we don't use it but we keep it just in case
c     FLShape = 1.d0 /
c    2   ( 1.d0 + ((Plg + Plgoff)/(Width + asym*Plg))**2 )

c    Gaussian:
      FLShape = dexp( -((Plg + Plgoff)/(Width + asym*Plg))**2 )

c    finally return the solution; but constrain argument
c    to exponential to avoid overflows

      arg = -FLShape*FcentLn
      if (arg.gt.100.d0) arg = 100.d0
      if (arg.lt.-100.d0) arg = -100.d0

      getDef = dexp(arg)

      return
      end

c***INITTROE************************************
c    companion routine to getTroe for initializing variables
      subroutine initTroe(Npar,X0,Xscale,Fscale,IDstr)
      implicit none

c    local variables
      integer iter
      integer Npar
      real*8 X0(*),Xscale(*),Fscale
      character IDstr*(*)
c    end declarations

      Npar = 4

c    no better idea so set Xscale to 1 - see IMSL description;
c    set all x's to order 1 (we define T's as x's * 1000)

c    enter FIT function ID string here
      IDstr = '4-Parameter Troe Fit in T/1000 (w/atan)'

      do 10 iter = 1,Npar
         Xscale(iter) = 1.d0
         X0(iter) = 5.d-1
10    continue
      Fscale = 1.d0

c    Now we reset the a factor
      X0(1) = 4.d-1

      return
      end

c***GETTROE************************************
c    function to return defaulf F fitting function
      function getTroe(Npar,X,T,P)
      implicit none

c    local variables
      integer Npar,icoef
      real*8 X(*),T,P,TStar(3),a,Fcent,
     2   FLShape,Rn,Rc,Rd,Plg
      real*8 getTroe
c    end declarations

      data Rd /0.14/

c    we try atan function to constrain a between 0 and 1
      a = datan(X(1))/3.14159265359d0 + .5d0

c    we wish to constrain T's to be positive
      do 10 icoef = 1,3
         Tstar(icoef) = 1.d3*dabs(X(icoef+1))
10    continue

c    use base 10 log for P
      Plg = dlog10(P)

c    this is right out of the Chemkin II manual, pg 23
      Fcent = (1.d0 - a)*dexp(-T/TStar(3)) + a*dexp(-T/TStar(1))
     2   + dexp(-TStar(2)/T)

      Rc = -0.4d0 - 0.67d0*dlog10(Fcent)
      Rn = 0.75 - 1.27*dlog10(Fcent)

      FLShape = 1.d0 /
     2   ( 1.d0 + ( (Plg + Rc)/ (Rn - Rd*(Plg + Rc)) )**2 )

      getTroe = Fcent**(FLShape)

      return
      end

c***INITSRI************************************
c    companion routine to getSRI for initializing variables
      subroutine initSRI(Npar,X0,Xscale,Fscale,IDstr)
      implicit none

c    local variables
      integer iter
      integer Npar
      real*8 X0(*),Xscale(*),Fscale
      character IDstr*(*)
c    end declarations

      Npar = 3

c    no better idea so set Xscale to 1 - see IMSL description;

c    enter FIT function ID string here
      IDstr = '3-Parameter SRI Fit in T/1000'

      do 10 iter = 1,Npar
         Xscale(iter) = 1.d0
         X0(iter) = 1.d0
10    continue
      Fscale = 1.d0

c    Now we reset the a factor
      X0(1) = 4.d-1

      return
      end

c***GETSRI************************************
c    function to return defaulf F fitting function
      function getSRI(Npar,X,T,P)
      implicit none

c    local variables
      integer Npar
      real*8 X(Npar),T,P,a,b,c,Fcent,FLShape,Plg
      real*8 getSRI
c    end declarations

c    (tried an equiv statement here, but compiler didn't
c    like it - probably because of passed variable)

      a = dabs(X(1))
      b = 1.d3*dabs(X(2))
      c = 1.d3*dabs(X(3))

c    use base 10 log for P
      Plg = dlog10(P)

c    this is right out of the Chemkin II manual, pg 24 -
c    but we have not included the optional d and e factors,
c    because we don't see that kind of limiting behavior

      Fcent = a*dexp(-b/T) + dexp(-T/c)

      FLShape = 1.d0 / ( 1.d0 + Plg**2 )

      getSRI = Fcent**(FLShape)

      return
      end

c  we tried these subs but we don't use them anymore
c**BOUNCR************************************
c    func to reflect fitting parameter to right of of boundary
      function bouncR(xval,xlim,xzone)
      implicit none

      real*8 xval,xlim,xzone
      real*8 bouncR

c    bouncR is the reflection of xval to the other side
c    of xlim but within the region defined by xzone

c    xzone should be small compared to expected magnitude of
c    xval (because legitimate xval's could be stuck in the
c    zone); still xzone should be large enough that things
c    put there don't simple vanish into the noise

c    the test condition (xval < xlim) is assumed

      bouncR = xlim + xzone/(1.d0 + (xlim - xval))

      return
      end

c**BOUNCL************************************
c    func to reflect fitting parameter to left of of boundary
      function bouncL(xval,xlim,xzone)
      implicit none

      real*8 xval,xlim,xzone
      real*8 bouncL

c    bouncL is the reflection of xval to the other side
c    of xlim but within the region defined by xzone

c    xzone should be small compared to expected magnitude of
c    xval (because legitimate xval's could be stuck in the
c    zone); still xzone should be large enough that things
c    put there don't simple vanish into the noise

c    the test condition (xval > xlim) is assumed

      bouncL = xlim - xzone/(1.d0 + (xval - xlim))

      return
      end


text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module pao_array_utility
! ------------------------------------------------------------------------------
! Code area 11: basis functions
! ------------------------------------------------------------------------------

!!****h* Conquest/pao_array_utility
!!  NAME
!!   pao_array_utility
!!  PURPOSE
!!   contains routines to reproject an array of given grid spacing onto
!!   a different grid spacing
!!  AUTHOR
!!   R Choudhury
!!  CREATION DATE
!!   24/07/03
!!  MODIFICATION HISTORY
!!   2008/02/06 08:30 dave
!!    Changed for output to file not stdout
!!   2019/08/16 12:12 dave
!!    Now uses splines, also changed name
!!  SOURCE
!!
module pao_array_utility

  use global_module, ONLY: io_lun

  implicit none

!!***
contains

!!****f* pao_array_utility/reproject_array *
!!
!!  NAME
!!   reproject_array
!!  USAGE
!!   reproject_array(yin,del_old,n_old,yout,del_new,n_new)
!!  PURPOSE
!!   Reprojects an array from a given grid spacing onto new grid spacing
!!  INPUTS
!!   yin - array to be reprojected
!!   del_old - old grid spacing
!!   n_old - old no of grid points
!!   yout - reprojected array values
!!   del_new - new array grid spacing
!!   n_new - new no of grid points
!!  USES
!!   datatypes
!!  AUTHOR
!!   R Choudhury
!!  CREATION DATE
!!   24/07/03
!!  MODIFICATION HISTORY
!!
!!  SOURCE
!!
  subroutine reproject_array(yin,del_old,n_old,yout,del_new,n_new)

    use datatypes
    use numbers
    use splines, ONLY: spline_nonU

    implicit none

    integer, intent(in) :: n_old,n_new
    real(double), intent(in) :: del_old,del_new
    real(double), intent(in), dimension(n_old) :: yin
    real(double), intent(out), dimension(n_new) :: yout

    integer :: i,j,l,m
    real(double), allocatable, dimension(:) :: y2,xin
    real(double) :: yp1,ypn, x,y, xx
    real(double) :: a, b, c, d, r1, r2, r3, r4

    allocate(xin(n_old))
    allocate(y2(n_old))
    !setting up x grid for spline_new subroutine
    do i=1,n_old
       xin(i) = (i-1)*del_old
    enddo

    !setting first derivatives of endpoints for spline
    !nb this DOES introduce an inaccuracy (depending on grid fineness)
    yp1 = (yin(2)-yin(1))/del_old
    ypn = (yin(n_old)-yin(n_old-1))/del_old

    call spline_nonU(n_old,xin,yin,yp1,ypn,y2)

    do i=1,n_new
       x=real(i-1,double)*del_new
       j = floor(x/del_old) + 1
       if(j+1<=n_old) then
          xx = real(j,double)*del_old
          a = (xx - x)/del_old
          b = one - a
          c = a * ( a * a - one ) * del_old * del_old / six
          d = b * ( b * b - one ) * del_old * del_old / six
          r1 = yin(j)
          r2 = yin(j+1)
          r3 = y2(j)
          r4 = y2(j+1)
          yout(i) = a*r1 + b*r2 + c*r3 + d*r4
       end if
    enddo
    deallocate(xin)
    deallocate(y2)

  end subroutine reproject_array
!!***

!!****f* pao_array_utility/matcharrays *
!!
!!  NAME
!!   matcharrays
!!  USAGE
!!   matcharrays(y1,n1,d1,y2,n2,d2,d12,n1_new,n2_new,n12)
!!
!!  PURPOSE
!!   Takes two arrays and projects them onto (same) specified grid spacing
!!  INPUTS
!!   y1, y2 - (values of) the two arrays
!!   n1, n2 - no of points in each array
!!   d1, d2 - spacing of original arrays
!!   d12 - new grid spacing
!!   n1_new, n2_new  - no of points in reprojected arrays
!!   n12 - size of reprojected arrays - power of two for FFT algorithms
!!  USES
!!   datatypes
!!  AUTHOR
!!   R Choudhury
!!  CREATION DATE
!!   24/07/03
!!  MODIFICATION HISTORY
!!
!!  SOURCE
!!
  subroutine matcharrays(y1,n1,d1,y2,n2,d2,d12,n1_new,n2_new,n12)

    use datatypes

    implicit none

    integer, intent(in) :: n1,n2,n12
    real(double), intent(inout), dimension(n12) :: y1,y2
    real(double), allocatable, dimension(:) :: ytemp
    real(double), intent(in) :: d1,d2,d12
    integer,intent(inout) :: n1_new,n2_new
    integer i

    n1_new = ((n1-1)*d1/d12)+1
    allocate(ytemp(n1_new))
    call reproject_array(y1(1:n1),d1,n1,ytemp,d12,n1_new)
    y1(1:n1_new) = ytemp(1:n1_new)
    y1(n1_new+1:n12) = 0.0_double
    deallocate(ytemp)

    n2_new = ((n2-1)*d2/d12)+1
    allocate(ytemp(n2_new))
    call reproject_array(y2(1:n2),d2,n2,ytemp,d12,n2_new)
    y2(1:n2_new) = ytemp(1:n2_new)
    y2(n2_new+1:n12) = 0.0_double
    deallocate(ytemp)

  end subroutine matcharrays
!!***

end module pao_array_utility

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
      module matdat_def_m
c
        implicit none
c
        integer, parameter :: MAXMAT  = 6
     &,                       MAXPROP = 10
     &,                       MAXTS   = 100
        real*8  :: datmat(MAXMAT,MAXPROP,MAXTS)
        integer :: nummat
     &,            matflg(MAXMAT,MAXTS)
     &,            mat_tag(MAXMAT,MAXTS)
     &,            mat_eos(MAXMAT,MAXTS)
        real*8  :: mat_prop(MAXMAT,MAXPROP,MAXTS)
        logical :: mexist
        real*8  :: datelas(1,2)
        common /matdat/ datmat,   matflg, mat_tag, mat_eos,
     &                  mat_prop, nummat, mexist,  datelas
        integer, parameter :: ieos_ideal_gas = 1
     &,                       ieos_ideal_gas_2 = 2
     &,                       ieos_liquid_1  = 3
     &,                       ieos_solid_1   = 4
c
        integer, parameter :: iprop_ideal_gas_mw    = 1  ! molecular weight
     &,                       iprop_ideal_gas_gamma = 2
     &,                       iprop_gas_visc  = 3
     &,                       iprop_liquid_1_rho_ref = 1
     &,                       iprop_liquid_1_p_ref   = 2
     &,                       iprop_liquid_1_T_ref   = 3
     &,                       iprop_liquid_1_cv      = 4
     &,                       iprop_liquid_1_alphaP  = 5
     &,                       iprop_liquid_1_betaT   = 6
c
c NOTE: This should only be an integer array BUT
c       since vtk reader does not support int yet,
c       so, we cast the integer numbers into this
c       real array!!!
c
        real*8, pointer :: mattype_interior(:)
c
      end module matdat_def_m


text:
!
! barrier program
!
! Example of a barrier directive
!
! 2008
! Creation: A. Muller, 2008
!
! deadlock if N not a multiple of nb processes

      program barrier
      implicit none
      integer N
      parameter (N=12)
      integer i,a(N)

!$omp parallel do
      do 10 i=1,N
         a(i)=i
         print *,i
!$omp barrier
 10   continue

      print *, a
      end


text:
!=======================================================================
!
!  Calculate the column density (cm^-2) of each species along each
!  HEALPix ray to the PDR surface (or simulation boundary) for all
!  particles. Calculations can be performed using either the Euler
!  or Trapezium integration scheme (specified by a compiler flag).
!
!-----------------------------------------------------------------------
SUBROUTINE CALCULATE_COLUMN_DENSITIES(NPART,NRAYS,NSPEC,PARTICLE)

   USE HEALPIX_TYPES
   USE PARTICLE_MODULE

   IMPLICIT NONE

   INTEGER(KIND=I4B),   INTENT(IN)    :: NPART,NRAYS,NSPEC
   TYPE(PARTICLE_TYPE), INTENT(INOUT) :: PARTICLE(*)

   INTEGER(KIND=I4B) :: I,J,K,L,M,P

   REAL(KIND=DP) :: STEP_SIZE

!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(DYNAMIC) &
!$OMP    SHARED(NPART,NRAYS,NSPEC,PARTICLE) &
!$OMP    PRIVATE(I,J,K,L,M,STEP_SIZE)
   DO P=1,NPART ! Loop over particles

!     Calculate the total hydrogen column density (N_H) and
!     column density of each species along each HEALPix ray
      PARTICLE(P)%TOTAL_COLUMN = 0.0D0
      PARTICLE(P)%COLUMN_DENSITY = 0.0D0
      DO J=0,NRAYS-1 ! Loop over rays
         DO I=2,PARTICLE(P)%NPOINT(J) ! Loop over evaluation points
            L = PARTICLE(P)%EVALUATION_POINT(J,I-1)
            M = PARTICLE(P)%EVALUATION_POINT(J,I)
            STEP_SIZE = SQRT((PARTICLE(M)%COORDINATES(1)-PARTICLE(L)%COORDINATES(1))**2 &
                         & + (PARTICLE(M)%COORDINATES(2)-PARTICLE(L)%COORDINATES(2))**2 &
                         & + (PARTICLE(M)%COORDINATES(3)-PARTICLE(L)%COORDINATES(3))**2)
#ifdef EULER
!           Total hydrogen column density using the Euler integration method
            PARTICLE(P)%TOTAL_COLUMN(J) = PARTICLE(P)%TOTAL_COLUMN(J) &
                                      & + PARTICLE(M)%GAS_DENSITY &
                                      &  *STEP_SIZE
#else
!           Total hydrogen column density using the Trapezium integration method (default)
            PARTICLE(P)%TOTAL_COLUMN(J) = PARTICLE(P)%TOTAL_COLUMN(J) &
                                       & + (PARTICLE(L)%GAS_DENSITY+PARTICLE(M)%GAS_DENSITY)/2.0 &
                                       &  *STEP_SIZE
#endif
            DO K=1,NSPEC ! Loop over species
#ifdef EULER
!              Column density of each species using the Euler integration method
               PARTICLE(P)%COLUMN_DENSITY(J,K) = PARTICLE(P)%COLUMN_DENSITY(J,K) &
                                             & + PARTICLE(M)%ABUNDANCE(K) &
                                             &  *PARTICLE(M)%GAS_DENSITY &
                                             &  *STEP_SIZE
#else
!              Column density of each species using the Trapezium integration method (default)
               PARTICLE(P)%COLUMN_DENSITY(J,K) = PARTICLE(P)%COLUMN_DENSITY(J,K) &
                                             & + (PARTICLE(L)%ABUNDANCE(K)*PARTICLE(L)%GAS_DENSITY &
                                             &   +PARTICLE(M)%ABUNDANCE(K)*PARTICLE(M)%GAS_DENSITY)/2.0 &
                                             &  *STEP_SIZE
#endif
            END DO ! End of loop over species
         END DO ! End of loop over evaluation points
      END DO ! End of loop over rays
   END DO ! End of loop over particles
!$OMP END PARALLEL DO

   RETURN
END SUBROUTINE CALCULATE_COLUMN_DENSITIES
!=======================================================================

max_stars_repo_name: uclchem/uclpdr

max_issues_repo_name: uclchem/uclpdr

max_forks_repo_name: uclchem/uclpdr


text:
subroutine lmn_cvgunk(a,itask)
   !-----------------------------------------------------------------------
   ! DESCRIPTION
   !    This routine performs several convergence checks for the
   !    low Mach number equations
   !-----------------------------------------------------------------------
   use MPI
   use typre
   use def_parame
   use Mod_LowMach
   implicit none
   class(LowMachProblem)   :: a
   integer(ip), intent(in) :: itask
   integer(ip), save       :: ipass=0
   real(rp)                :: rplmn,rvlmn,rtlmn,ripre,rivel,ritem,rilmn
   integer(ip)             :: npoinLocal,ndime,ierr

   select case(itask)

   !Check convergence of the inner iterations, always in the L2 norm:
   !|| u(n,i,j) - u(n,i,j-1)|| / ||u(n,i,j)||
   case(1)

      call a%InnerResiduals(rvlmn,rplmn,rtlmn)
      rilmn = max(rvlmn,rtlmn,rplmn)
      a%rilmn = rilmn
      call MPI_ALLREDUCE(MPI_IN_PLACE,a%rilmn,1,MPI_REAL8,MPI_MAX,a%MPIcomm,ierr)

      if(((a%rilmn<a%cotol) .or. (a%itera>=a%maxit))) a%kfl_goite = 0
      if (a%MPIrank == a%MPIroot) then
         if(ipass==0) then
            ipass=1
            write(a%lun_conve,100)
            write(a%lun_outph,200)
            write(a%lun_outph,201)
            write(a%lun_nolin,300)
         end if

         write(a%lun_conve,101) a%istep,a%itera,a%ctime,rvlmn,rtlmn,rplmn

         write(a%lun_outph,202) a%istep,a%itera,a%ctime,&
            a%tamin,a%tamax,a%tamea,&
            a%remin,a%remax,a%remea
         ! Write in the log file.
         if (a%kfl_flush == 1) then
            call flush(a%lun_outph)
            call flush(a%lun_conve)
         endif
      end if

   !Check convergence of the outer iterations in the norm selected by the user:
   !|| u(n,i,*) - u(n,i-1,*)|| / ||u(n,i,*)||
   case(2)
      call a%Mesh%GetNpoinLocal(npoinLocal)
      call a%Mesh%GetNdime(ndime)
      call vecresMPI(a%kfl_normc,ndime*npoinLocal,a%veloc(1,1,1),a%veloc(1,1,2),a%resid,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      call vecresMPI(a%kfl_normc,npoinLocal,      a%press(1,1),  a%press(1,2),  a%resip,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      call vecresMPI(a%kfl_normc,npoinLocal,      a%tempe(1,1),  a%tempe(1,2),  a%resit,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      a%resid = max(a%resiv,a%resit,a%resip)

   !Check residual of the time evolution, always in the L2 norm:
   !|| u(n,*,*) - u(n-1,*,*)|| / ||u(n,*,*)||
   case(3)
      call a%Mesh%GetNpoinLocal(npoinLocal)
      call a%Mesh%GetNdime(ndime)
      call vecresMPI(two,ndime*npoinLocal,a%veloc(1,1,1),a%veloc(1,1,3),rivel,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      call vecresMPI(two,npoinLocal,      a%press(1,1),  a%press(1,3),  ripre,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      call vecresMPI(two,npoinLocal,      a%tempe(1,1),  a%tempe(1,3),  ritem,zelmn,a%kfl_MPIComType,a%MPIcomm,a%MPIroot,a%MPIrank,a%MPIsize)
      rilmn = max(rivel,ritem,ripre)

      if (a%MPIrank == a%MPIroot) then
         write(a%lun_conve,101) a%istep,0,a%ctime,rivel,ritem,ripre
         if(rilmn.le.a%sstol) then
            a%kfl_stead = 1
            write(a%lun_conve,102) a%istep
         end if
      endif
      CALL MPI_BCAST(a%kfl_stead, 1, MPI_INTEGER4, a%MPIroot, a%MPIcomm, ierr)

   end select

   !Formats
100 format('$ ','  Time','      Inner',&
       &      '       Current','      Velocity','    Temperature','     Pressure',/,&
       & '$ ','  step','  iteration',&
       &      '          time','      residual','      residual','      residual')
103 format('$ ','  Time','     Global','      Inner',' Non-Linear',&
       &      '       Current','      Velocity','    Temperature','     Pressure',/,&
       & '$ ','  step','  iteration','  iteration',&
       &      '          time','      residual','      residual','      residual')

101 format(i9,1x,i9,5(2x,e12.6),1(i9))
104 format(2x,4(2x,i9),3(2x,e12.6),4x,a)
102 format('$ >>>  VELOCITY IS STATIONARY AT TIME STEP ',i5)

200 format(///,&
       5x, '>>>  EVOLUTION OF NUMERICAL PARAMETERS:',///,&
       10x,a)
201 format('$ ','       Time','      Inner',&
       &      '       Current','    Tau_min       Tau_max',&
       &      '       Tau_mean       Re_min        Re_max',&
       &      '        Re_mean       yp_min        yp_max',&
       &      '        yp_mean',/,&
       & '$ ','       step','  iteration',&
       &      '          time')
202 format(4x,i9,2x,i9,7(2x,e12.6))
300 format('$ ','   Step','   Iteration','     SGS iterations',&
       &        '    Non Linear Subscales Maximum Residual',/,&
       &        45x,  ' Velocity','     Temperature')

end subroutine lmn_cvgunk

subroutine lmn_InnerResiduals(a,rvlmn,rplmn,rtlmn)
   use typre
   use Mod_LowMach
   implicit none
   class(LowMachProblem) :: a
   real(rp) :: rvlmn,rplmn,rtlmn

   integer(ip) :: ndime, npoinLocal

   call a%Mesh%GetNpoinLocal(npoinLocal)
   call a%Mesh%GetNdime(ndime)
   call vecresMPIHeterogeneous(2_ip,ndime+2,ndime,npoinLocal,1,1,ndime,a%unkno,a%veloc(1,1,1),rvlmn,zelmn,a%MPIcomm,a%MPIroot)
   call vecresMPIHeterogeneous(2_ip,ndime+2,1,npoinLocal,ndime+2,1,1,a%unkno,a%press(1,1),rplmn,zelmn,a%MPIcomm,a%MPIroot)
   call vecresMPIHeterogeneous(2_ip,ndime+2,1,npoinLocal,ndime+1,1,1,a%unkno,a%tempe(1,1),rtlmn,zelmn,a%MPIcomm,a%MPIroot)

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


module seismogram

  use common, only : wp
  use datatypes, only : block_grid_t, seismogram_type
  implicit none

contains

  subroutine init_seismogram(input,S,name,G)

    implicit none

    integer,intent(in) :: input
    type(block_grid_t), intent(in) :: G
    type(seismogram_type),intent(inout) :: S
    character(*),intent(in) :: name

    integer :: mx,my,mz,px,py,pz
    logical :: output_exact_moment
    logical :: output_seismograms,output_fault_topo,output_fields_block1,output_fields_block2
    integer :: stride_fields,n,stat
    character(256) :: temp,filename
    character(3) :: field

    !real(kind = wp) :: xmin, xmax, ymin, ymax, zmin, zmax

    namelist /output_list/ output_exact_moment, output_seismograms, output_fault_topo,   &
                           output_fields_block1,output_fields_block2,stride_fields

    mx = G%C%mq
    my = G%C%mr
    mz = G%C%ms
    px = G%C%pq
    py = G%C%pr
    pz = G%C%ps

    ! defaults

    output_exact_moment = .false.
    output_seismograms = .false.
    output_fields_block1 = .false.
    output_fields_block2 = .false.
    stride_fields = 1

    rewind(input)
    read(input,nml=output_list,iostat=stat)
    if (stat>0) stop 'error reading namelist output_list'

    S%output_exact_moment = output_exact_moment
    S%output_seismograms = output_seismograms
    S%output_fields_block1 = output_fields_block1
    S%output_fields_block2 = output_fields_block2
    S%stride_fields = stride_fields

    ! initialization for seismogram output

    if (S%output_seismograms) then

       ! move to start of station list

       if (S%block_num == 1) then
          rewind(input)
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---begin:station_listU---') exit
          end do

          ! determine number of stations in list

          S%nstations = 0
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---end:station_listU---') exit
             S%nstations = S%nstations+1
          end do
       end if

       if (S%block_num == 2) then
          rewind(input)
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---begin:station_listV---') exit
          end do

          ! determine number of stations in list

          S%nstations = 0
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---end:station_listV---') exit
             S%nstations = S%nstations+1
          end do
       end if

       ! allocate station indices array and output file unit array

       if (S%output_exact_moment) then
         allocate(S%file_unit(2*S%nstations))
       else
         allocate(S%file_unit(S%nstations))
       end if

         allocate(S%i(S%nstations),S%j(S%nstations),S%k(S%nstations))
         allocate(S%i_phys(S%nstations),S%j_phys(S%nstations), &
              S%k_phys(S%nstations))

       ! read station list again, this time storing station indices

       rewind(input)

       if(S%block_num == 1) then
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---begin:station_listU---') exit
          end do
       end if

       if(S%block_num == 2) then
          do
             read(input,'(a)',end=100) temp
             if (temp=='!---begin:station_listV---') exit
          end do
       end if

      if(S%nstations > 0) then
       do n = 1,S%nstations
          read(input,*) S%i_phys(n),S%j_phys(n),S%k_phys(n)
          !print *, S%i_phys(n),S%j_phys(n),S%k_phys(n)
       end do

       !xmin = minval(G%X(mx:px,my:py,mz:pz,1))
       !xmax = maxval(G%X(mx:px,my:py,mz:pz,1))
       !ymin = minval(G%X(mx:px,my:py,mz:pz,2))
       !ymax = maxval(G%X(mx:px,my:py,mz:pz,2))
       !zmin = minval(G%X(mx:px,my:py,mz:pz,3))
       !zmax = maxval(G%X(mx:px,my:py,mz:pz,3))
       S%i(:) = -10000
       S%j(:) = -10000
       S%k(:) = -10000

       call Find_Coordinates(G%X, S%i_phys, S%j_phys, S%k_phys, &
            S%i, S%j, S%k, S%nstations, mx, my, mz, px, py, pz)

       ! open file units for output

       do n = 1,S%nstations
          if (S%i(n) > 0 .and. S%j(n) > 0 .and. S%k(n) > 0) then
            write(filename,'(a,i0,a,i0,a,i0,a,i0,a)') trim(adjustl(name)) // '_', &
                S%i(n),'_',S%j(n),'_',S%k(n),'_block',S%block_num,'.dat'
            open(newunit=S%file_unit(n),file=filename)
          end if
       end do

       if( S%output_exact_moment) then
       do n = 1,S%nstations
          if (S%i(n) > 0 .and. S%j(n) > 0 .and. S%k(n) > 0) then
            write(filename,'(a,i0,a,i0,a,i0,a,i0,a)') trim(adjustl(name)) // '_exact_', &
                S%i(n),'_',S%j(n),'_',S%k(n),'_block',S%block_num,'.dat'
            open(newunit=S%file_unit(S%nstations+n),file=filename)
          end if
       end do
       end if
     end if

    end if

    ! initialization for body fields output

    if (S%output_fields_block1.or.S%output_fields_block2) then

       ! open file units for output

       do n = 1,9

          select case(n)
          case(1)
             field = 'vx'
          case(2)
             field = 'vy'
          case(3)
             field = 'vz'
          case(4)
             field = 'sxx'
          case(5)
             field = 'sxy'
          case(6)
             field = 'sxz'
          case(7)
             field = 'syy'
          case(8)
             field = 'syz'
          case(9)
             field = 'szz'
          end select

          if (S%output_fields_block1 .and. S%block_num == 1) then
             write(filename,'(a,a,a)') trim(adjustl(name)) // '_block1_',trim(adjustl(field)),'.dat'
             open(newunit=S%file_unit_block1(n),file=filename)
          end if

          if (S%output_fields_block2 .and. S%block_num == 2) then
             write(filename,'(a,a,a)') trim(adjustl(name)) // '_block2_',trim(adjustl(field)),'.dat'
             open(newunit=S%file_unit_block2(n),file=filename)
          end if

       end do

    end if

    return

100 stop 'error reading seismogram station list'

  end subroutine init_seismogram

  subroutine write_seismogram(S, t, F)

    use mpi3dbasic, only : rank
    implicit none

    type(seismogram_type),intent(in) :: S
    real(kind = wp),intent(in) :: t
    real(kind = wp), dimension(:,:,:,:), allocatable, intent(in) :: F
    integer :: n

    if (S%output_seismograms .and. S%nstations > 0) then

       do n = 1,S%nstations
          if (S%i(n) > 0 .and. S%j(n) > 0 .and. S%k(n) > 0) then
            write(S%file_unit(n),'(4f15.10)') t, F(S%i(n),S%j(n),S%k(n),1:3)
          end if
       end do

    end if

    if (S%output_fields_block1 .and. S%block_num == 1) then
       do n = 1,9
          write(S%file_unit_block1(n),*) F(:,:,:,n)
          print *, S%file_unit_block1(n)
          flush(S%file_unit_block1(n))
       end do

    end if

  end subroutine write_seismogram

  subroutine destroy_seismogram(S)

    implicit none

    type(seismogram_type),intent(inout) :: S

    integer :: n

    if (S%output_seismograms) then

      if (S%nstations > 0) then
         do n = 1,S%nstations
            if (S%i(n) > 0 .and. S%j(n) > 0 .and. S%k(n) > 0) close(S%file_unit(n))
         end do
       end if

      if (allocated(S%i)) deallocate(S%i,S%j,S%k,S%file_unit)

    end if

    if (S%output_fields_block1) then

       do n = 1,9
          close(S%file_unit_block1(n))
       end do

    end if

    if (S%output_fields_block2) then

       do n = 1,9
          close(S%file_unit_block2(n))
       end do

    end if

  end subroutine destroy_seismogram

  subroutine Find_Coordinates(XX, x1, y1, z1, x_i, y_j, z_k, nstations, mx, my, mz, px, py, pz)

    ! Given the physical positions of the receivers x1, y1, z1
    ! Find the corresponding indices in x_i, y_j, z_k in the mesh XX

    implicit none

    integer, intent(in) :: mx, my, mz, px, py, pz               ! size of the grid-block
    integer, intent(in) :: nstations
    real(kind = wp), dimension(:), intent(in) :: x1, y1, z1                ! receiver  positions
    integer, dimension(:), intent(out) :: x_i, y_j, z_k         ! spatial indices of receiver  positions to be found
    real(kind = wp), dimension(:,:,:,:), allocatable, intent(in) :: XX                  ! grid
    integer :: i, j, k, c,  i0, j0, k0
    real(kind = wp) :: vec(3), dist, mindist,xmin, xmax, ymin, ymax, zmin, zmax
    real(kind = wp) :: hx,hy,hz, dist0

    xmin = minval(XX(mx:px,my:py,mz:pz,1))
    xmax = maxval(XX(mx:px,my:py,mz:pz,1))
    ymin = minval(XX(mx:px,my:py,mz:pz,2))
    ymax = maxval(XX(mx:px,my:py,mz:pz,2))
    zmin = minval(XX(mx:px,my:py,mz:pz,3))
    zmax = maxval(XX(mx:px,my:py,mz:pz,3))

    hx = 1d0
    hy = 1d0
    hz = 1d0

    !print *,  xmin, xmax, ymin, ymax, zmin, zmax
    !print *, nstations
    do c = 1, nstations

       i0 = -9999
       j0 = -9999
       k0 = -9999

       mindist = 1.0e8_wp

       !if ((xmin <= x1(c) .and. x1(c) <= xmax) .and. &
       !     (ymin <= y1(c) .and. y1(c) <= ymax) .and. &
       !     (zmin <= z1(c) .and. z1(c) <= zmax)) then

          k_loop: do k = mz, pz
             do j = my, py
                do i = mx, px

                   hx = XX(px, j, k, 1)-XX(px-1, j, k, 1)
                   hy = XX(i, py, k, 2)-XX(i, py-1, k, 2)
                   hz = XX(i, j, pz, 3)-XX(i, j, pz-1, 3)

                   if (i<px) hx = XX(i+1, j, k, 1)-XX(i, j, k, 1)
                   if (j<py) hy = XX(i, j+1, k, 2)-XX(i, j, k, 2)
                   if (k<pz) hz = XX(i, j, k+1, 3)-XX(i, j, k, 3)

                   vec = [XX(i, j, k, 1)-x1(c), XX(i, j, k, 2)-y1(c), XX(i, j, k, 3)-z1(c)]
                   dist = sqrt(dot_product(vec, vec))
                   dist0 = sqrt(hx**2 + hy**2 + hz**2)

                   !if (((abs(XX(i, j, k, 1)-x1(c))<=0.8_wp*hx) .and. &
                   !     (abs(XX(i, j, k, 2)-y1(c))<=0.8_wp*hy) .and. &
                   !     (abs(XX(i, j, k, 3)-z1(c))<=0.8_wp*hz)) .and. &
                   !      (dist <= mindist)) then

                   if ((dist <= mindist) .and. (dist <= 0.5*dist0))  then
                   !if ((dist <= mindist))  then

                      i0 = i
                      j0 = j
                      k0 = k

                      x_i(c) = i0
                      y_j(c) = j0
                      z_k(c) = k0

                      mindist = dist

                   end if

                end do
             end do

          end do k_loop
          !print*, mindist, c, i0, j0, k0, XX(i0, j0, k0, 1), XX(i0, j0, k0, 2), XX(i0, j0, k0, 3), x1(c), y1(c), z1(c)
       !end if
       if (i0 > 0 .and. j0>0 .and. k0 > 0) then
          print*, mindist, c, i0, j0, k0, XX(i0, j0, k0, 1), XX(i0, j0, k0, 2), XX(i0, j0, k0, 3), x1(c), y1(c), z1(c)
       end if

        !x_i(c) = i0
        !y_j(c) = j0
        !z_k(c) = k0
     end do

  end subroutine Find_Coordinates

    subroutine Find_Coordinates_moment(XX, x1, y1, z1, x_i, y_j, z_k, nstations,nq,nr,ns, mx, my, mz, px, py, pz)

    ! Given the physical positions of the receivers x1, y1, z1
    ! Find the corresponding indices in x_i, y_j, z_k in the mesh XX

    implicit none

    integer, intent(in) :: mx, my, mz, px, py, pz, nq, nr, ns               ! size of the grid-block
    integer, intent(in) :: nstations
    real(kind = wp), dimension(:), intent(in) :: x1, y1, z1                ! receiver  positions
    integer, dimension(:), intent(out) :: x_i, y_j, z_k         ! spatial indices of receiver  positions to be found
    real(kind = wp), dimension(:,:,:,:), allocatable, intent(in) :: XX                  ! grid
    integer :: i, j, k, c
    real(kind = wp) :: vec(3), dist, mindist,xmin, xmax, ymin, ymax, zmin, zmax
    real(kind = wp) :: hx,hy,hz

    xmin = minval(XX(mx:px,my:py,mz:pz,1))
    xmax = maxval(XX(mx:px,my:py,mz:pz,1))
    ymin = minval(XX(mx:px,my:py,mz:pz,2))
    ymax = maxval(XX(mx:px,my:py,mz:pz,2))
    zmin = minval(XX(mx:px,my:py,mz:pz,3))
    zmax = maxval(XX(mx:px,my:py,mz:pz,3))

    hx = 1d0
    hy = 1d0
    hz = 1d0

    do c = 1, nstations
       x_i(c) = -9999
       y_j(c) = -9999
       z_k(c) = -9999

       if ((xmin <= x1(c) .and. x1(c) <= xmax) .and. &
            (ymin <= y1(c) .and. y1(c) <= ymax) .and. &
            (zmin <= z1(c) .and. z1(c) <= zmax)) then

          mindist = 1.0e8_wp
          !        hx = XX(px, j, k, 1)-XX(px-1, j, k, 1)
          !        hy = XX(i, py, k, 1)-XX(i, py-1, k, 2)
          !        hz = XX(i, j, pz, 1)-XX(i, j, pz-1, 3)

          k_loop: do k = mz, pz
             do j = my, py
                do i = mx, px

                   !print*, i,j,k,XX(i,j,k,1)

                   if (i<nq) hx = XX(i+1, j, k, 1)-XX(i, j, k, 1)
                   if (j<nr) hy = XX(i, j+1, k, 2)-XX(i, j, k, 2)
                   if (k<ns) hz = XX(i, j, k+1, 3)-XX(i, j, k, 3)

                   vec = [XX(i, j, k, 1)-x1(c), XX(i, j, k, 2)-y1(c), XX(i, j, k, 3)-z1(c)]
                   dist = sqrt(dot_product(vec, vec))

                   !                    if (((abs(XX(i, j, k, 1)-x1(c))<=0.9d0*hx) .and. &
                   !                         (abs(XX(i, j, k, 2)-y1(c))<=0.9d0*hy) .and. &
                   !                         (abs(XX(i, j, k, 3)-z1(c))<=0.9d0*hz))) then

                   if (mx == 1 .and. px == 25 .and. my == 1 .and. py == 25) then
                      !print*,XX(i,j,k,1)
                      ! print*,'dist_min = :',i,j,k,dist
                   end if
                   if (dist <= mindist) then

                      x_i(c) = i
                      y_j(c) = j
                      z_k(c) = k

                      mindist = dist

                   end if
                end do
             end do
          end do k_loop
       end if
    end do

  end subroutine Find_Coordinates_moment

end module seismogram


text:
      SUBROUTINE CHKSN4(MBDCND,NBDCND,ALPHA,BETA,COFX,SINGLR)
C
C     THIS SUBROUTINE CHECKS IF THE PDE   SEPX4
C     MUST SOLVE IS A SINGULAR OPERATOR
C
      COMMON /SPL4/   KSWX       ,KSWY       ,K          ,L          ,
     1                AIT        ,BIT        ,CIT        ,DIT        ,
     2                MIT        ,NIT        ,IS         ,MS         ,
     3                JS         ,NS         ,DLX        ,DLY        ,
     4                TDLX3      ,TDLY3      ,DLX4       ,DLY4
      LOGICAL         SINGLR
      EXTERNAL COFX
      SINGLR = .FALSE.
C
C     CHECK IF THE BOUNDARY CONDITIONS ARE
C     ENTIRELY PERIODIC AND/OR MIXED
C
      IF ((MBDCND.NE.0 .AND. MBDCND.NE.3) .OR.
     1    (NBDCND.NE.0 .AND. NBDCND.NE.3)) RETURN
C
C     CHECK THAT MIXED CONDITIONS ARE PURE NEUMAN
C
      IF (MBDCND .NE. 3) GO TO  10
      IF (ALPHA.NE.0.0 .OR. BETA.NE.0.0) RETURN
   10 CONTINUE
C
C     CHECK THAT NON-DERIVATIVE COEFFICIENT FUNCTIONS
C     ARE ZERO
C
      DO  30 I=IS,MS
         XI = AIT+FLOAT(I-1)*DLX
         CALL COFX (XI,AI,BI,CI)
         IF (CI .NE. 0.0) RETURN
   30 CONTINUE
C
C     THE OPERATOR MUST BE SINGULAR IF THIS POINT IS REACHED
C
      SINGLR = .TRUE.
      RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
Lieutenant Thomas W. Waltz, badge number 51, is the investigations lieutenant for the Davis Police Department, where he has worked since at least 1994. He was also the Vice President of the Davis Police Officers Association ( http://www.city.davis.ca.us/meetings/councilpackets/20040802/03K_Middle_Income_Housing.pdf source), which gives an indication of the respect his fellow officers have for him.  He used to be the watch commander.

Hes a decorated police officer, and is known to many for his kind, helpful, and responsible service. He is also 411 tall. Waltz was mentioned in an http://www.maxim.com/articles/index.aspx?a_id3551 article in wiki:wikipedia::Maxim Magazine on the Davis Police Department Trading Cards.

He has provided the training on behalf of the Police Department for the Community Emergency Response Team.  He also worked at the 2005 Picnic Day.

20061024 17:01:07 nbsp So yeah, the officer Waltz chased one of my friends down 45 years ago and caught him in a foot pursuit Users/StevenDaubert

20061024 17:01:43 nbsp oh, hes deep into some martial arts, so he will handle you with the proverbial quickness if you try to step to him ^____^ Users/StevenDaubert

20070301 19:20:40 nbsp I saw him today. It was magical. Users/MollyeKt

20070505 23:42:10 nbsp Was waltz at the protest with the single round tear gas launcher? Users/StevenDaubert

20110203 07:19:06 nbsp Saw the LT on tv online because Linda Wolf was one if my wifes bridesmaides...very composed and in control. I am willing to donate time money or both from my agency to Davis PD departments...someone should ping me jeremy.b.belk@gmail.com Users/JeremyBrandonBelk

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
great_ne(jj1,bb1).
great_ne(z1,p1).
great_ne(bb1,h1).
great_ne(ee1,m1).
great_ne(cc1,b1).
great_ne(dd1,a1).
great_ne(y1,b1).
great_ne(k1,b1).
great_ne(y1,h1).
great_ne(q1,kk1).
great_ne(ee1,hh1).
great_ne(bb1,m1).
great_ne(cc1,dd1).
great_ne(z1,ee1).
great_ne(k1,ii1).
great_ne(kk1,y1).
great_ne(cc1,hh1).
great_ne(kk1,hh1).
great_ne(y1,w1).
great_ne(ff1,a1).
great_ne(cc1,w1).
great_ne(jj1,k1).
great_ne(x1,v1).
great_ne(u1,v1).
great_ne(cc1,bb1).
great_ne(m1,y1).
great_ne(jj1,c1).
great_ne(q1,p1).
great_ne(cc1,a1).
great_ne(y1,ee1).
great_ne(ii1,u1).
great_ne(a1,f1).
great_ne(v1,z1).
great_ne(ff1,u1).
great_ne(h1,kk1).
great_ne(hh1,t1).
great_ne(aa1,p1).
great_ne(x1,cc1).
great_ne(jj1,cc1).
great_ne(h1,x1).
great_ne(y1,u1).
great_ne(m1,q1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE WRITE_INTEGER_VEC ( ARRAY_DESCR, INT_VEC, NROWS )

! Writes an integer vector to F06 in a format that has 10 terms across the page (repeated until vector is completely written)

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  WRITE_MATRIX_BY_COLS_BEGEND

      USE WRITE_INTEGER_VEC_USE_IFs                        ! Corrected 2019/07/14 (was WRITE_VECTOR_USE_IFs)

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'WRITE_INTEGER_VEC'
      CHARACTER(LEN=*), INTENT(IN)    :: ARRAY_DESCR       ! Character descriptor of the integer array to be printed
      CHARACTER(131*BYTE)             :: HEADER            ! MAT_DESCRIPTOR centered in a line of output

      INTEGER(LONG), INTENT(IN)       :: NROWS             ! Number of rows in matrix MATOUT
      INTEGER(LONG), INTENT(IN)       :: INT_VEC(NROWS)    ! Integer vector to write out
      INTEGER(LONG)                   :: I,J               ! DO loop indices
      INTEGER(LONG)                   :: INT_VEC_LINE(10)  ! 10 members of INT_VEC
      INTEGER(LONG)                   :: NUM_LEFT          ! Count of the number of rows of INT_VEC left to write out
      INTEGER(LONG)                   :: PAD               ! Number of spaces to pad in HEADER to center MAT_DESCRIPTOR
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = WRITE_MATRIX_BY_COLS_BEGEND

      INTRINSIC                       :: LEN

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      PAD = (132 - LEN(ARRAY_DESCR))/2
      HEADER(1:) = ' '
      HEADER(PAD+1:) = ARRAY_DESCR
      WRITE(F06,101) HEADER

      NUM_LEFT = NROWS
      DO I=1,NROWS,10
         IF (NUM_LEFT >= 10) THEN
            DO J=1,10
               INT_VEC_LINE(J) = INT_VEC(I+J-1)
            ENDDO
            WRITE(F06,103) (INT_VEC_LINE(J),J=1,10)
         ELSE
            DO J=1,NUM_LEFT
               INT_VEC_LINE(J) = INT_VEC(I+J-1)
            ENDDO
            WRITE(F06,103) (INT_VEC_LINE(J),J=1,NUM_LEFT)
         ENDIF
         NUM_LEFT = NUM_LEFT - 10
      ENDDO
      WRITE(F06,*)

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  101 FORMAT(1X,/,1X,A,/)

  103 FORMAT(16X,10(I10))

! **********************************************************************************************************************************

      END SUBROUTINE WRITE_INTEGER_VEC


text:
Karen Spiess served as a Trustee for the Davis JUSD School Board from 1987 to 1995.  She was initially appointed by the board in January 1987 to fill out the remainder of Helen Thomson Helen Thomsons term when she was elected to the Yolo County Board of Supervisors.  Spiess ran unopposed for her seat the following November and was reelected in 1991.

Source:

McEwing, Shearon; 16 Oct. 1991, Softspoken candidate is outspoken on the issue of education, Davis Enterprise, pg. A1.


text:
PROGRAM test

   IMPLICIT NONE

   integer             :: nkpts
   real                :: kpoints(3,1)
   real                :: amat(3,3)
   real                :: bmat(3,3)
   integer             :: num(3)
   integer             :: num_bands
   integer             :: ntype
   character(len=50)   :: seedname
   integer             :: num_atoms
   character(len=2)    :: atom_symbols(1)
   logical             :: gamma_only
   logical             :: spinors
   integer             :: nntot
   integer             :: nnlist(1,1)
   integer             :: nncell(3,1,1)
   integer             :: num_bands2
   integer             :: num_wann2
   real                :: proj_site(3,1)
   integer             :: proj_l(1)
   integer             :: proj_m(1)
   integer             :: proj_radial(1)
   real                :: proj_z(3,1)
   real                :: proj_x(3,1)
   real                :: proj_zona(1)
   integer             :: exclude_bands(1)
   real                :: pos(3,1)

   ! Taken from wannier90-1.2/src/wannier_lib.F90
   interface
   subroutine wannier_setup(seed__name, mp_grid_loc, num_kpts_loc,&
                            real_lattice_loc, recip_lattice_loc, kpt_latt_loc,&
                            num_bands_tot, num_atoms_loc, atom_symbols_loc,&
                            atoms_cart_loc, gamma_only_loc, spinors_loc, nntot_loc,&
                            nnlist_loc, nncell_loc, num_bands_loc, num_wann_loc,&
                            proj_site_loc, proj_l_loc, proj_m_loc, proj_radial_loc,&
                            proj_z_loc, proj_x_loc, proj_zona_loc, exclude_bands_loc)

      implicit none

      integer, parameter :: dp = selected_real_kind(15,300)
      integer, parameter :: num_nnmax=12
      character(len=*), intent(in) :: seed__name
      integer, dimension(3), intent(in) :: mp_grid_loc
      integer, intent(in) :: num_kpts_loc
      real(kind=dp), dimension(3,3), intent(in) :: real_lattice_loc
      real(kind=dp), dimension(3,3), intent(in) :: recip_lattice_loc
      real(kind=dp), dimension(3,num_kpts_loc), intent(in) :: kpt_latt_loc
      integer, intent(in) :: num_bands_tot
      integer, intent(in) :: num_atoms_loc
      character(len=*), dimension(num_atoms_loc), intent(in) :: atom_symbols_loc
      real(kind=dp), dimension(3,num_atoms_loc), intent(in) :: atoms_cart_loc
      logical, intent(in) :: gamma_only_loc
      logical, intent(in) :: spinors_loc

      integer, intent(out) :: nntot_loc
      integer, dimension(num_kpts_loc,num_nnmax), intent(out) :: nnlist_loc
      integer,dimension(3,num_kpts_loc,num_nnmax), intent(out) :: nncell_loc
      integer, intent(out) :: num_bands_loc
      integer, intent(out) :: num_wann_loc
      real(kind=dp), dimension(3,num_bands_tot), intent(out) :: proj_site_loc
      integer, dimension(num_bands_tot), intent(out) :: proj_l_loc
      integer, dimension(num_bands_tot), intent(out) :: proj_m_loc
      integer, dimension(num_bands_tot), intent(out) :: proj_radial_loc
      real(kind=dp), dimension(3,num_bands_tot), intent(out) :: proj_z_loc
      real(kind=dp), dimension(3,num_bands_tot), intent(out) :: proj_x_loc
      real(kind=dp), dimension(num_bands_tot), intent(out) :: proj_zona_loc
      integer, dimension(num_bands_tot), intent(out) :: exclude_bands_loc
   end subroutine wannier_setup
   end interface

   seedname='blahblahblah'
   num = 0
   nkpts = 1
   amat = 0.0
   bmat = 0.0
   kpoints = 0.0
   num_bands = 0
   num_atoms = 0
   atom_symbols = ':)'
   pos = 0.0
   gamma_only = .FALSE.
   spinors = .FALSE.
   nntot = 0
   nnlist = 0
   nncell = 0
   num_bands2 = 0
   num_wann2 = 0
   proj_site = 0.0
   proj_l = 0
   proj_m = 0
   proj_radial = 0
   proj_z = 0.0
   proj_x = 0.0
   proj_zona = 0.0
   exclude_bands = 0

   CALL wannier_setup(seedname,num,&
                      nkpts,&
                      transpose(amat),bmat,&
                      kpoints,num_bands,&
                      num_atoms,atom_symbols,pos,&
                      gamma_only,spinors,&
                      nntot,nnlist,nncell,num_bands2,&
                      num_wann2,&
                      proj_site,proj_l,proj_m,&
                      proj_radial,proj_z,&
                      proj_x,proj_zona,exclude_bands)

END

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
! { dg-do compile }
! { dg-options "-fcray-pointer" }

subroutine err1
  integer :: in_common1, in_common2, v, w, equiv1, equiv2
  common /in_common1/ in_common1
  pointer (ipt1, in_common1)		! { dg-error "conflicts with COMMON" }
  pointer (ipt2, in_common2)
  common /in_common2/ in_common2	! { dg-error "conflicts with COMMON" }
  equivalence (v, equiv1)
  pointer (ipt3, equiv1)		! { dg-error "conflicts with EQUIVALENCE" }
  pointer (ipt4, equiv2)
  equivalence (w, equiv2)		! { dg-error "conflicts with EQUIVALENCE" }
end subroutine err1

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text: The Mel Gagnon Family Bench is located in the ValleyWise Garden just north west of the Arboretum Gazebo. There are several benches in the garden.