text:
program ex4p1
    use mt_mod
    implicit none
    integer, parameter:: ik=selected_int_kind(11), rk=selected_real_kind(11,100)
    integer(kind=ik) :: i, j, h, l, Np, runs=500
    real(kind=rk), parameter :: pi = 4.0*atan(1.0), k=8.6173324E-5, T=600
    real(kind=rk) :: x, x0, x1, trial, state, u, u1, a=1/real(2,rk)*k*T, mean, E_theor=3/2*k*T, sumrun=0, x_max, b
    real(kind=rk), allocatable :: E_MC(:), means(:), rms(:)

    b=sqrt(2/pi)*1/(a**3)

    call init_genrand64(int(123465, ik))

    open(1, action="write", file="data_max.txt", status="replace", position="append")
    open(2, action="write", file="data_points.txt", status="replace", position="append")

    allocate(means(runs))
    allocate(rms(1000))

    !Data for (E_max, RMS) plot with N=10^3
    Np = 10**3
    allocate(E_MC(Np))
    !Vary E_max from 0.001 to 0.1
    do l=1, 100
        !Markov Chain
        x_max=0.001*l       !0b
        do h=1, runs
            !Inital config.
            x = 3/real(2,rk)*k*T    !0a
            do i=1, Np
                state = rho(x)      !0c
                x0 = x              !1
                u = genrand64_real1()*2-1   !2
                x1 = x + u*x_max    !3
                trial = rho(x1)     !4

                if (trial/state >= 1) then
                    x = x1
                    state = trial   !5
                else if (trial/state >= genrand64_real1()) then
                    x = x1
                    state = trial   !5
                else
                    x = x0
                end if
                E_MC(i) = abs(x)    !6
            end do

            means(h)=sum(E_MC)/real(Np, rk) !Mean energy of distribution
        end do

        !Sum inside RMS
        do i=1, runs
            sumrun=sumrun+(means(i)-E_theor)**2
        end do

        !RMS for x_max=0.001*l
        rms(l)=sqrt(1/real(runs,rk)*sumrun)
        sumrun=0

        print *, l*0.001,"done."
    end do

    !Save data
    do l=1, 1000
        write(1, *) rms(h), 0.001*l
    end do
    deallocate(E_MC)
    deallocate(rms)

    print*, "E_max done"

    allocate(rms(6))

    !Set x_max according to previous results
    x_max=0.05

    !Markov Chain as before, but now with looping over N_points
    do j=2, 7
        Np = 10**j
        allocate(E_MC(Np))
        do h=1, runs
            x = 3/real(2,rk)*k*T
            do i=1, Np
                x0 = x
                E_MC(i) = abs(x)
                state = rho(x)
                u = genrand64_real1()*2-1
                x1 = x + u*x_max
                trial = rho(x1)
                if (trial/state >= 1) then
                    x = x1
                    state = trial
                else if (trial/state >= genrand64_real1()) then
                    x = x1
                    state = trial
                else
                    x = x0
                end if
            end do
            means(h)=sum(E_MC)/real(Np, rk)
        end do
        do i=1, runs
            sumrun=sumrun+(means(i)-E_theor)**2
        end do
        rms(j-1)=sqrt(1/real(runs,rk)*sumrun)
        sumrun=0
        write(2, *) rms(j-1), 10**j
        deallocate(E_MC)
        print*, j, "-power done"
    end do

    close(1)
    close(2)

contains

real(kind=rk) function rho(E)
    implicit none
    real(kind=rk), intent(in) :: E
    rho=b*E**2*exp(-1/real(2,rk)*(E/a)**2)
end function rho
end program ex4p1

max_stars_repo_name: toicca/basics-of-mc

max_issues_repo_name: toicca/basics-of-mc

max_forks_repo_name: toicca/basics-of-mc


text:
*DECK SGER
      SUBROUTINE SGER (M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
C***BEGIN PROLOGUE  SGER
C***PURPOSE  Perform rank 1 update of a real general matrix.
C***LIBRARY   SLATEC (BLAS)
C***CATEGORY  D1B4
C***TYPE      SINGLE PRECISION (SGER-S)
C***KEYWORDS  LEVEL 2 BLAS, LINEAR ALGEBRA
C***AUTHOR  Dongarra, J. J., (ANL)
C           Du Croz, J., (NAG)
C           Hammarling, S., (NAG)
C           Hanson, R. J., (SNLA)
C***DESCRIPTION
C
C  SGER   performs the rank 1 operation
C
C     A := alpha*x*y' + A,
C
C  where alpha is a scalar, x is an m element vector, y is an n element
C  vector and A is an m by n matrix.
C
C  Parameters
C  ==========
C
C  M      - INTEGER.
C           On entry, M specifies the number of rows of the matrix A.
C           M must be at least zero.
C           Unchanged on exit.
C
C  N      - INTEGER.
C           On entry, N specifies the number of columns of the matrix A.
C           N must be at least zero.
C           Unchanged on exit.
C
C  ALPHA  - REAL            .
C           On entry, ALPHA specifies the scalar alpha.
C           Unchanged on exit.
C
C  X      - REAL             array of dimension at least
C           ( 1 + ( m - 1)*abs( INCX)).
C           Before entry, the incremented array X must contain the m
C           element vector x.
C           Unchanged on exit.
C
C  INCX   - INTEGER.
C           On entry, INCX specifies the increment for the elements of
C           X. INCX must not be zero.
C           Unchanged on exit.
C
C  Y      - REAL             array of dimension at least
C           ( 1 + ( n - 1 )*abs( INCY ) ).
C           Before entry, the incremented array Y must contain the n
C           element vector y.
C           Unchanged on exit.
C
C  INCY   - INTEGER.
C           On entry, INCY specifies the increment for the elements of
C           Y. INCY must not be zero.
C           Unchanged on exit.
C
C  A      - REAL             array of DIMENSION ( LDA, n ).
C           Before entry, the leading m by n part of the array A must
C           contain the matrix of coefficients. On exit, A is
C           overwritten by the updated matrix.
C
C  LDA    - INTEGER.
C           On entry, LDA specifies the first dimension of A as declared
C           in the calling (sub) program. LDA must be at least
C           max( 1, m ).
C           Unchanged on exit.
C
C***REFERENCES  Dongarra, J. J., Du Croz, J., Hammarling, S., and
C                 Hanson, R. J.  An extended set of Fortran basic linear
C                 algebra subprograms.  ACM TOMS, Vol. 14, No. 1,
C                 pp. 1-17, March 1988.
C***ROUTINES CALLED  XERBLA
C***REVISION HISTORY  (YYMMDD)
C   861022  DATE WRITTEN
C   910605  Modified to meet SLATEC prologue standards.  Only comment
C           lines were modified.  (BKS)
C***END PROLOGUE  SGER
C     .. Scalar Arguments ..
      REAL               ALPHA
      INTEGER            INCX, INCY, LDA, M, N
C     .. Array Arguments ..
      REAL               A( LDA, * ), X( * ), Y( * )
C     .. Parameters ..
      REAL               ZERO
      PARAMETER        ( ZERO = 0.0E+0 )
C     .. Local Scalars ..
      REAL               TEMP
      INTEGER            I, INFO, IX, J, JY, KX
C     .. External Subroutines ..
      EXTERNAL           XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C***FIRST EXECUTABLE STATEMENT  SGER
C
C     Test the input parameters.
C
      INFO = 0
      IF     ( M.LT.0 )THEN
         INFO = 1
      ELSE IF( N.LT.0 )THEN
         INFO = 2
      ELSE IF( INCX.EQ.0 )THEN
         INFO = 5
      ELSE IF( INCY.EQ.0 )THEN
         INFO = 7
      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
         INFO = 9
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'SGER  ', INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
     $   RETURN
C
C     Start the operations. In this version the elements of A are
C     accessed sequentially with one pass through A.
C
      IF( INCY.GT.0 )THEN
         JY = 1
      ELSE
         JY = 1 - ( N - 1 )*INCY
      END IF
      IF( INCX.EQ.1 )THEN
         DO 20, J = 1, N
            IF( Y( JY ).NE.ZERO )THEN
               TEMP = ALPHA*Y( JY )
               DO 10, I = 1, M
                  A( I, J ) = A( I, J ) + X( I )*TEMP
   10          CONTINUE
            END IF
            JY = JY + INCY
   20    CONTINUE
      ELSE
         IF( INCX.GT.0 )THEN
            KX = 1
         ELSE
            KX = 1 - ( M - 1 )*INCX
         END IF
         DO 40, J = 1, N
            IF( Y( JY ).NE.ZERO )THEN
               TEMP = ALPHA*Y( JY )
               IX   = KX
               DO 30, I = 1, M
                  A( I, J ) = A( I, J ) + X( IX )*TEMP
                  IX        = IX        + INCX
   30          CONTINUE
            END IF
            JY = JY + INCY
   40    CONTINUE
      END IF
C
      RETURN
C
C     End of SGER  .
C
      END

ext: f


text:
implicit none
integer :: i,n
real :: y=0, term

read *, n
do i=1,n
	term=1./(i*(i+1))
	y=y+term
	if(abs(term)<=1.e-5) exit
end do
print *, 'sum',y,' term=',i

end

max_stars_repo_name: quchunguang/test


text:
module linalg_mod
    use, intrinsic :: iso_fortran_env, only : output_unit, error_unit, DP => REAL64
    implicit none

    private

    public :: read_array, write_array, generate_random_array

contains

    subroutine read_array(file_name, array, dimensions)
        implicit none
        character(len=*), intent(in) :: file_name
        real(kind=DP), dimension(..), allocatable, intent(out) :: array
        integer, dimension(:), intent(in) :: dimensions

        if (rank(array) /= size(dimensions)) then
            write (unit=error_unit, fmt='(A)') '# error: array rank does not match dimensions'
            stop 3
        end if
        select rank(array)
            rank(1)
                call read_vector(file_name, dimensions(1), array)
            rank(2)
                call read_matrix(file_name, dimensions(1), dimensions(2), array)
            rank default
                write (unit=error_unit, fmt='(A, I0, A)') &
                    '# error: rank ', rank(array), ' not supported'
                stop 4
        end select
    end subroutine read_array

    subroutine write_array(array, file_name)
        implicit none
        real(kind=DP), dimension(..), intent(in) :: array
        character(len=*), intent(in), optional :: file_name

        select rank(array)
            rank(1)
                call write_vector(array, file_name)
            rank(2)
                call write_matrix(array, file_name)
            rank default
                write (unit=error_unit, fmt='(A, I0, A)') &
                    '# error: rank ', rank(array), ' not supported'
                stop 4
        end select
    end subroutine write_array

    subroutine read_vector(file_name, n, vector)
        implicit none
        character(len=*), intent(in) :: file_name
        integer, value :: n
        real(kind=DP), dimension(:), allocatable, intent(out) :: vector
        integer :: i, status, unit_nr
        character(len=1024) :: msg

        open (newunit=unit_nr, file=trim(file_name), access='sequential', &
            action='read', status='old', iostat=status, iomsg=msg)
        if (status /= 0) then
            write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
            stop 1
        end if

        if (.not. allocated(vector)) then
            allocate (vector(n), stat=status)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 3
            end if
        end if

        do i = 1, n
            read (unit=unit_nr, fmt=*, iostat= status, iomsg=msg) vector(i)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 2
            end if
        end do

        close (unit=unit_nr)

    end subroutine read_vector

    subroutine read_matrix(file_name, nr_rows, nr_cols, matrix)
        implicit none
        character(len=*), intent(in) :: file_name
        integer, value :: nr_rows, nr_cols
        real(kind=DP), dimension(:, :), allocatable, intent(out) :: matrix
        integer :: i, status, unit_nr
        character(len=1024) :: msg

        open (newunit=unit_nr, file=trim(file_name), access='sequential', &
            action='read', status='old', iostat=status, iomsg=msg)
        if (status /= 0) then
            write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
            stop 1
        end if

        if (.not. allocated(matrix)) then
            allocate (matrix(nr_rows, nr_cols), stat=status)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 3
            end if
        end if

        do i = 1, nr_rows
            read (unit=unit_nr, fmt=*, iostat= status, iomsg=msg) matrix(i, :)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 2
            end if
        end do

        close (unit=unit_nr)

    end subroutine read_matrix

    subroutine write_vector(vector, file_name)
        implicit none
        real(kind=DP), dimension(:), intent(in) :: vector
        character(len=*), intent(in), optional :: file_name
        integer :: unit_nr, i, status
        character(len=1024) :: msg

        if (present(file_name)) then
            open (newunit=unit_nr, file=trim(file_name), access='sequential', &
                action='write', status='replace', iostat=status, iomsg=msg)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 1
            end if
        else
            unit_nr = output_unit
        end if

        do i = 1, size(vector)
            write (unit=unit_nr, fmt='(E27.15)') vector(i)
        end do

        if (present(file_name)) close (unit_nr)
    end subroutine write_vector

    subroutine write_matrix(matrix, file_name)
        implicit none
        real(kind=DP), dimension(:, :), intent(in) :: matrix
        character(len=*), intent(in), optional :: file_name
        integer :: unit_nr, i, status
        character(len=1024) :: msg

        if (present(file_name)) then
            open (newunit=unit_nr, file=trim(file_name), access='sequential', &
                action='write', status='replace', iostat=status, iomsg=msg)
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') '# error: ', trim(msg)
                stop 1
            end if
        else
            unit_nr = output_unit
        end if

        do i = 1, size(matrix, 1)
            write (unit=unit_nr, fmt='(*(E27.15))') matrix(i, :)
        end do

        if (present(file_name)) close (unit_nr)
    end subroutine write_matrix

    subroutine generate_random_array(array, dimensions)
        implicit none
        real(kind=DP), dimension(..), allocatable, intent(out) :: array
        integer, dimension(:), intent(in) :: dimensions
        integer :: status

        if (rank(array) /= size(dimensions)) then
            write (unit=error_unit, fmt='(A)') '# error: array rank does not match dimensions'
            stop 3
        end if
        select rank(array)
            rank(1)
                if (.not. allocated(array)) then
                    allocate(array(dimensions(1)), stat=status)
                    if (status /= 0) then
                        write (unit=error_unit, fmt='(A, I0)') &
                            '# error: can not allocate vector of size ', dimensions(1)
                        stop 6
                    end if
                end if
                call random_number(array)
            rank(2)
                if (.not. allocated(array)) then
                    allocate(array(dimensions(1), dimensions(2)), stat=status)
                    if (status /= 0) then
                        write (unit=error_unit, fmt='(A,I0, A, I0)') &
                            '# error: can not allocate matrix of size ', &
                            dimensions(1), ' by ', dimensions(2)
                        stop 6
                    end if
                end if
                call random_number(array)
            rank default
                write (unit=error_unit, fmt='(A, I0, A)') &
                    '# error: rank ', rank(array), ' not supported'
                stop 4
        end select
    end subroutine generate_random_array

end module linalg_mod

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
C
C $Id: pcffhl.f,v 1.6 2008-07-27 00:17:19 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE PCFFHL (NCI,JNDX,XORIG,YORIG,NCO,JNDXO,RDGU)
C
      include 'pcffdx.h'
      include 'pcffsp.h'
C
C Process holes.
C
      DIMENSION RDGU(*)
C
      JNDXO = JNDX
      NDX = NCI
C
   10 CONTINUE
      NDX = NDX+1
C
      IF (SFLGS(NDX) .EQ. 3) THEN
C
C  End hole, connect the last point back to the initial point if not
C  an outline font.
C
        IF (OUTLIN .EQ. 0) THEN
          JNDXO = JNDXO+1
          RDGU(JNDXO) = XORIG
          JNDXO = JNDXO+1
          RDGU(JNDXO) = YORIG
        ENDIF
C
C  Terminate the line if an outline font has been requested and
C  set the next point in the new line back to the original
C  point upon entry to this subroutine.
C
        IF (OUTLIN .EQ. 1) THEN
          JNDXO = JNDXO+1
          RDGU(JNDXO) = -2048.
          JNDXO = JNDXO+1
          RDGU(JNDXO) = 0.
          JNDXO = JNDXO+1
          RDGU(JNDXO) = XORIG
          JNDXO = JNDXO+1
          RDGU(JNDXO) = YORIG
        ENDIF
        NCO = NDX
        RETURN
      ELSE IF (SFLGS(NDX) .EQ. 0) THEN
C
C  Add coordinate to the current path.
C
        JNDXO = JNDXO+1
        RDGU(JNDXO) = XC(NDX)
        JNDXO = JNDXO+1
        RDGU(JNDXO) = YC(NDX)
        GO TO 10
      ELSE IF (SFLGS(NDX) .EQ. 4) THEN
C
C  Process the Bezier curve.
C
          BCNTLX(1) = RDGU(JNDXO-1)
          BCNTLY(1) = RDGU(JNDXO)
          DO 30 IB=2,4
            BCNTLX(IB) = XC(NDX+IB-2)
            BCNTLY(IB) = YC(NDX+IB-2)
   30     CONTINUE
C
          CALL BCCURV(BCNTLX,BCNTLY,IBZL,BZXC,BZYC,NBP)
          DO 60 K=2,NBP
            JNDXO = JNDXO+1
            RDGU(JNDXO) = BZXC(K)
            JNDXO = JNDXO+1
            RDGU(JNDXO) = BZYC(K)
   60     CONTINUE
          NDX = NDX+2
          GO TO 10
      ELSE
          PRINT * , 'PCFFHL -- Invalid fontcap encoding encountered.'
          STOP
      ENDIF
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      FUNCTION irbit2(iseed)
      INTEGER irbit2,iseed,IB1,IB2,IB5,IB18,MASK
      PARAMETER (IB1=1,IB2=2,IB5=16,IB18=131072,MASK=IB1+IB2+IB5)
      if(iand(iseed,IB18).ne.0)then
        iseed=ior(ishft(ieor(iseed,MASK),1),IB1)
        irbit2=1
      else
        iseed=iand(ishft(iseed,1),not(IB1))
        irbit2=0
      endif
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/irbit2.for


text:
      subroutine UMESS(TEXT, MACT, IVAR)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-03-25 UMESS  Krogh  Included full set of parameters.
c>> 1992-01-29 UMESS  Krogh  Initial code.
c
c Dummy routine called by the error message routine MESS.  This gives
c the default action when an error message is requested, which is to do
c nothing here and take the default actions in MESS.  One can change the
c actions taken by changing the values in MACT or IVAR.  Comments here
c refer to things defined in the listing for MESS.
c Variables in the calling sequence are defined as follows:
c
c TEXT  Same as TEXT passed into MESS.  The characters up to the first
c       "$" sign are the name of the subroutine that made the call to
c       MESS.
c MACT  Contents of MACT connected with the error message.  The first
c       location is the value of K52 = 10*s + p, where s is the stop
c       level specified for the error message, and p the print level.
c       The next is the value of L52, which is an error index associated
c       with an error message, and the next is the value of M52.  These
c       values or a subset of them together with TEXT should serve to
c       uniquely identify the error message.  The values of MACT(1)
c       (=K52), or  MACT(2) (=L52) can be changed and the new values
c       will be used if a return is made to MESS.  (For example MACT(1)
c       = MACT(1) - mod(MACT(1), 10), would turn off printing of a
c       message.)  Note that if MACT(1) > 90, and it is replaced with a
c       value which results in MESS returning rather than stopping,
c       unpredictable results are to be expected.
c IVAR  An array containing the parameters defining default values for
c       MESS.  The parameter values below can be used to access data in
c       this array.  For example, IVAR(MESTOP) gives the current stop
c       value, so one can tell (using s from MACT(1) above) if an error
c       is going to stop, and perhaps do some clean up before the stop
c       occurs.  Values in IVAR can be changed, but this should be done
c       with great care, since if MESS doesn't stop, the next error or
c       diagnostic message will use the new value, and diagnostic
c       messages do not call this routine.
c
      integer   MESUNI, MEHEAD, MEDDIG, MEMLIN, MEELIN, MEMUNI, MEEUNI,
     1  MESCRN, MEDIAG, MEMAXE, MESTOP, MEPRNT, METDIG, MENTXT, MEIDAT,
     2  MEFDAT, MEMDAT, MEMDA1, MEMDA2, MEMDA3, MEMDA4, MEMDA5, METABS,
     3  MEERRS, MECONT, MERET , MEEMES, METEXT, METABL, MERES3, MEIVCI,
     4  MEIVEC, MEIMAT, MEJVCI, MEJVEC, MEJMAT, MEFVCI, MEFVEC, MEFMAT,
     5  MEGVCI, MEGVEC, MEGMAT, MEMAXI, MEGBAS, MEVBAS, MEVLAS, MEFSPV
c Parameters for changing the environment.
      parameter (MESUNI=10,MEHEAD=11,MEDDIG=12,MEMLIN=13,MEELIN=14,
     1 MEMUNI=15,MEEUNI=16,MESCRN=17,MEDIAG=18,MEMAXE=19,MESTOP=20,
     2 MEPRNT=21,METDIG=22,MENTXT=23,MEIDAT=24,MEFDAT=25,MEMDAT=26,
     3 MEMDA1=27,MEMDA2=28,MEMDA3=29,MEMDA4=30,MEMDA5=31,METABS=32,
     4 MEERRS=33)
c Parameters for actions.
      parameter (MECONT=50, MERET=51,MEEMES=52,METEXT=53,MEFSPV=54,
     1 METABL=55,MERES3=56,MEIVEC=57,MEIMAT=58,MEJVEC=59,MEJMAT=60,
     2 MEFVEC=61,MEFMAT=62,MEGVEC=63,MEGMAT=64,MEIVCI=65,MEJVCI=66,
     2 MEFVCI=67,MEGVCI=68)
c Parameter derived from those above.
      parameter (MEMAXI=68,MEGBAS=49,MEVBAS=10,MEVLAS=33)
c
      character*(*) TEXT(*)
      integer   MACT(*), IVAR(MEVBAS:MEVLAS)
c
c Code for your actions goes here.
c
      return
      end


text:
! When "len" appears in the parameter list of an include file the the parser fails.
! Also the name of the line number appears without a correct file name.
  function nf90_def_dim(ncid, name, len, dimid)
    integer,             intent( in) :: ncid
    character (len = *), intent( in) :: name
    integer,             intent( in) :: len
    integer,             intent(out) :: dimid
    integer                          :: nf90_def_dim

    nf90_def_dim = nf_def_dim(ncid, name, len, dimid)
  end function nf90_def_dim

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
The Art in Public Places program has sponsored many pieces of town art in Davis.  According to the http://www.city.davis.ca.us City of Davis website, Since 1973, the City has set aside 1% of all funds allocated for Capital Improvement Projects to purchase art that enhances the quality of the Davis environment. The artists selected through this program represent a diversity of style and media. Sculpture, murals, paintings, tile work, and functional pieces are located around the city for your enjoyment.  The city maintains a walking guide to the art located around the city, only available in print (not online).

The aesthetics of many of these pieces are often critiqued by a diverse range of people.

Sponsored Works

This is an incomplete list

  Arroyo Park UFO, 2004.
  Alluvium, 2003
  Clepsydra, 2003
  Nurturing the Seeds of Tomorrow, 2002
  Cnawan Stone, 1995
  Ancient Shadows, 1995
  Davis Family Mural, 1991
  Watchman ＃8 Watchman #8, 1988
  The Loneliness of Integrity, 1986
  Zenobias Throne

For art in public places that was not sponsored by this program, see Urban Art.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
20080905 19:18:30 nbsp Welcome to the Wiki. Thanks for adding the information about the River City Rickshaws rickshaws. Users/JasonAller

20100314 18:57:49 nbsp You ... are everywhere, haha! Users/ChristyMarsden


text:
c     tmpl matching to calculate sigma_TOA, for rms calu both assuming re-chi^2=1 and from baseline in frequency domain

      subroutine tmplmatch(profile,template,nbins,shift,eshift,
     +                     eshift2,scale)

      implicit none

      integer nbins
      real*8 profile(nbins),template(nbins)
      real*8 eshift, shift, snr, esnr,eshift2,scale
      real*8 amp(nbins/2), pha(nbins/2)

      call makeFtempl(amp,pha,template,nbins)

c-----Do the frequency domain template fit, get shift in rad
      call fftfit(profile,amp,pha,nbins,shift,eshift,snr,esnr,
     +            eshift2,scale)

c     if(snr<5.0) write(*,*)"Warning: individual channel SNR low"

      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine makeFtempl(amp,pha,da,nbins)

C     Converts a standard profile into normalized harmonic-squared form
c
c     based on a subroutine by AW
c

      implicit none

      integer maxsam,i,nh,ipos, nbins
      real*8 fixoff, twopi,da,amp,pha,pmax

      parameter (MAXSAM=4096,twopi=3.141592653589793*2,fixoff=0.0)

      dimension da(nbins),amp(MAXSAM/2),pha(MAXSAM/2)

      complex cstd(0:MAXSAM/2)

      nh=nbins/2

c-------Computefrequency-domaintemplate
      call cprof_old(da,nbins,nh,cstd,amp,pha)

      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine cprof_old(y,nmax,nh,c,amp,pha)

C  Compute FFT of profile in array y(nmax), and return amplitude and phase
C  in arrays amp(nh) and pha(nh).  Note that nh=nmax/2, and that the DC term
C  is returned in c(0), fundamental in c(1), ..., Nyquist freq in c(nh).

	parameter(MAXSAM=4096)
	real*8 y(nmax),amp(nh),pha(nh)
	complex c(0:nh),temp(MAXSAM)

	do 10 i=1,nh
10	 temp(i)=cmplx(y(2*i-1),y(2*i))
	call ffft_cdft(temp,nmax,1,1)
	c(0)=temp(1)
	do 20 i=1,nh
	 c(i)=temp(i+1)
	 amp(i)=cabs(c(i))
	 pha(i)=0.
20	if(amp(i).gt.0.) pha(i)=aimag(clog(c(i)))
	return
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine fftfit(prof,s,phi,nmax,shift,eshift,snr,esnr,eshift2,scale)

C  Fourier transform domain routine for determining pulse TOAs.
C  Input data:
C	prof(nmax)	profile
C	s(nh)		standard profile amplitude
C	phi(nh)	standard profile phase
C	nmax		length of prof

C  Outputs:
C	shift		shift required to align std with prof, in bins
C	eshift	uncertainty in shift
C	snr		signal to noise ratio of prof
C	esnr		uncertainty in snr

C  Method:
C  It is assumed that prof(j)=a + b*std(j-tau) + noise(j), where the
C  (j-tau) subscript is supposed to mean that the time-shift between the
C  observed and standard profiles need not be an integer number of bins.
C  The algorithm is a straightforward Chi-squared minimization with respect
C  to a, b, and tau.  (The result for a, being of no physical or instrumental
C  interest, is not actually evaluated -- though it easily could be.)
C  First and second partial derivatives of Chisqr with respect to b and tau
C  are computed, and used to evaluate s, snr, and their "standard errors,"
C  or one-sigma uncertainties.  The only special trick is that the expression
C  for the best-fit value of tau defines it implicitly, and cannot be solved
C  analytically.  It is solved numerically instead, finding the minimum of
C  Chisqr near a best guess from a CCF at 32 lags done in the Fourier domain.

C  Also note that it may
C  be desirable to return the scale factor b relating prof to std, instead
C  of snr.  In that case you could also return the noise estimate rms.
c
c  Changed a lower bound to iterate the transcendental equation from 5 to 4,
c  in accord with changing nprof from 64 to 32 in fccf32.f.
c  AW, March 1992.

	parameter (twopi=3.141592653589793*2,MAXSAM=4096)
	real*8 prof(MAXSAM),p(MAXSAM/2),theta(MAXSAM/2)
	real*8 s(MAXSAM/2),phi(MAXSAM/2),r(MAXSAM/2),tmp(MAXSAM/2)
        real*8 shift,eshift,snr,esnr,eshift2,scale
	complex cp(0:MAXSAM/2)
	logical low,high

c       if (phi(1).ne.0) then
c	  write (0,*) ' Phase of fundamental not zero, check .hm file'
c	  stop
c	end if

	nh=nmax/2
	call cprof_old(prof,nmax,nh,cp,p,theta)

        rmsbegin =nh/8*5
c       calculate rms from the high frequency harmonics
        readrms=0.0
        do i=rmsbegin,nh
           readrms=readrms+(p(i))**2
        end do
        readrms=sqrt(readrms/(nh-rmsbegin+1))

c-------DefineTOAsuch that phase of fundamental = 0
c       do i=1,nh
c        theta(i)=mod(theta(i)-float(i)*theta(1),twopi)
c       enddo

	do 10 k=1,nh
	 tmp(k)=p(k)*s(k)
10	 r(k)=theta(k)-phi(k)

        fac=nmax/twopi

	call fccf_cdft(tmp,r,shift)

C  The "DO 60" loop solves the transcendental equation yielding the best-fit
C  value of tau.  Here the number starts at 16 (number used in CCF)
C  and increases by factors of 2, at each step finding the function
C  zero closest to the starting point from the previous iteration.

 	tau=shift

	do 60 isum=4,99
	 nsum=2.0**isum
	 if(nsum.gt.nh) go to 70
	 dtau=twopi/(nsum*5)
	 edtau=1./(2.*nsum+1.)
         if (nsum.gt.(nh/2.+.5)) edtau=1.e-5

	 ntries=0
	 low=.false.
	 high=.false.
50	 ftau=dchisqr(tau,tmp,r,nsum)
	 ntries=ntries+1
	 if(ftau.lt.0.0) then
	   a=tau
	   fa=ftau
	   tau=tau+dtau
	   low=.true.
	 else
	   b=tau
	   fb=ftau
	   tau=tau-dtau
	   high=.true.
	 end if
	 if (ntries.gt.10) then
	   shift=0.
	   eshift=999.
	   snr=0.
	   esnr=0.
	   return
	 end if
	 if (low.neqv.high) go to 50
	 tau=zbrent(a,b,fa,fb,edtau,tmp,r,nsum)
60	continue

70	s1=0.
	s2=0.
	s3=0.
	do 80 k=1,nh
	 cosfac=cos(-r(k)+k*tau)
	 s1=s1 + tmp(k)*cosfac
	 s2=s2 + s(k)**2
80	 s3=s3 + k**2 *tmp(k)*cosfac
	b=s1/s2
	s1=0.
	do 90 k=1,nh
	 sq=p(k)**2-2.*b*p(k)*s(k)*cos(r(k)-k*tau)+(b*s(k))**2
90	 if(s(k).ne.0.) s1=s1+sq
	rms=sqrt(s1/nh)
	errb=rms/sqrt(2.0*s2)
	errtau=rms/sqrt(2.0*b*s3)
	snr=2.0*sqrt(2.0*nh)*b/rms

	shift=tau
	eshift=errtau
	esnr=snr*errb/b
        scale=b

c       new eshift from readrms
	errtau=readrms/sqrt(2.0*b*s3)
        eshift2=errtau

        rchisq = (s1/readrms**2)/(nh-3)
c       write(*,*)"rms (re chi^2=1):",rms," real:",readrms,
c    +  " re chi^2:",rchisq

	return
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine fccf_cdft(amp,pha,shift)

C  Calculates CCF in Fourier domain using 8 harmonics of amp and pha arrays
C  amp=p*s,pha=theta-phi.  Finds maximum of CCF at 32 lags over the pulsar
C  period, and returns value of shift in radians.
c
c  Changed nprof from 64 to 32 (and all the appropriate hardwired constants)
c  to deal with 32-point profiles.  AW -- March 1992

	parameter (nprof=32,nprof1=nprof-1)
	parameter (MAXSAM=4096,twopi=3.141592653589793*2)
	real*8 amp(MAXSAM/2),pha(MAXSAM/2),shift
	complex ccf(0:nprof1)

	nh=nprof/2
	ccf(0)=(0.,0.)
	do 10 i=1,nh/2
	 ccf(i)=cmplx(amp(i)*cos(pha(i)),amp(i)*sin(pha(i)))
10	 ccf(nprof-i)=conjg(ccf(i))
	do 20 i=nh/2+1,nh
	 ccf(i)=(0.,0.)
20	 ccf(nprof-i)=(0.,0.)
	call ffft_cdft(ccf,nprof,-1,0)
	cmax=-1.e30
	do 30 i=0,nprof1
	 rc=real(ccf(i))
	 if (rc.gt.cmax) then
	   cmax=rc
	   imax=i
	 end if
30	continue

	fb=cmax
	ia=imax-1
	if(ia.eq.-1) ia=nprof-1
	fa=real(ccf(ia))
	ic=imax+1
	if(ic.eq.nprof) ic=0
	fc=real(ccf(ic))
	if ((2*fb-fc-fa).ne.0) then
	  shift=imax+0.5*(fa-fc)/(2*fb-fc-fa)
	else
	  shift=imax
	end if
c	if(shift.gt.nh) shift=shift-nprof
	shift=shift*twopi/nprof

	return
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	function zbrent(x1,x2,f1,f2,tol,tmp,pha,nsum)

C Brent's method root finding, calls dchisqr(x,tmp,r,nsum) function for fftfit
C Fit refined till output accuracy is tol

	parameter (itmax=100,eps=6.e-8,MAXSAM=4096)
	real*8 tmp(MAXSAM/2),pha(MAXSAM/2)

	a=x1
	b=x2
	fa=f1
	fb=f2
	fc=fb
	do 11 iter=1,itmax
	if(fb*fc.gt.0.) then
	  c=a
	  fc=fa
	  d=b-a
	  e=d
	end if
	if(abs(fc).lt.abs(fb)) then
	  a=b
	  b=c
	  c=a
	  fa=fb
	  fb=fc
	  fc=fa
	end if
	tol1=2.*eps*abs(b)+0.5*tol
	xm=.5*(c-b)
	if(abs(xm).le.tol1 .or. fb.eq.0.) then
	  zbrent=b
	  return
	end if
	if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
	  s=fb/fa
	  if(a.eq.c) then
	    p=2.*xm*s
	    q=1.-s
	  else
	    q=fa/fc
	    r=fb/fc
	    p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
	    q=(q-1.)*(r-1.)*(s-1.)
	  end if
	  if(p.gt.0.) q=-q
	  p=abs(p)
	  if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
	    e=d
	    d=p/q
	  else
	    d=xm
	    e=d
	  end if
	else
	  d=xm
	  e=d
	end if
	a=b
	fa=fb
	if(abs(d) .gt. tol1) then
	  b=b+d
	else
	  b=b+sign(tol1,xm)
	end if
	fb=dchisqr(b,tmp,pha,nsum)
11	continue
	zbrent=b
	return
	end

	function dchisqr(tau,tmp,r,nsum)

	parameter (MAXSAM=8192)
	real*8 tmp(MAXSAM/2),r(MAXSAM/2)

	s=0.
	do 40 k=1,nsum
40	 s=s+k*tmp(k)*sin(-r(k)+k*tau)
	dchisqr=s
	return
	end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


text:
An Elections Election was held on Tuesday, Nov. 2, 2010.

 Races

 Federal
  United States Senator
  United States Representative 1st District

 State
  Governor
  Lieutenant Governor
  Secretary of State
  Controller
  Treasurer
  Attorney General
  Insurance Commissioner
  Superintendent of Public Instruction
  Member of the State Assembly, 2nd District
  Member of the State Assembly, 8th District
  Member of the State Board of Equalization, District 1

 Propositions
  Proposition 19
  Proposition 20
  Proposition 21
  Proposition 22
  Proposition 23
  Proposition 24
  Proposition 25
  Proposition 26
  Proposition 27

 Judicial


text:
program test
  character(len=80)     :: filename
  integer               :: my_id

  my_id  = 0
  write(filename,'(A,I1,A)') 'OBJ_', my_id, '.dat'

  print *, filename

end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
Noahs Bagels, located in the Marketplace, is a neighborhood deli that serves award winning bagels, artisan breads, homestyle soups, signature sandwiches and tasty sweets. On Fridays, they also have freshly baked Challah loaves. The employees are friendly and typically comprised of college students.

The bagels here are baked fresh daily (the bagels are thrown out at endofday) Noahs will, at times, donate their dayoldbagels to charities if requested of the manager. Bagels are also madetoorder on the spot. It takes an hour to make a bagel, so allow leeway for that if requesting just out of the oven bagels for an order. On weekend mornings Noahs is usually very busy and the line extends to the door, but it moves so quickly that you are in and out quickly. One can even call in and preorder their dozen or a meal such as a sandwich or breakfast bagel so they can skip the long order line and just walk in and pay for and pick up their food at the register.

There is ample seating at Noahs. Inside is counter seating and well as the large patio outside. Many patrons sit outside to enjoy the day, read the newspaper, like the The Sacramento Bee Sac Bee, or just relax and enjoy the good food and friends. Noahs is dogs dog friendly as well. They offer bowls of water to the their furry friends and they also carry Doggie bagels, made especially with your canine in mind. The posts outside the store have metal latches to secure dog leashes for owner who wish to come inside. Everyone deserves a treat when they visit Noahs. Like every Noahs Bagels, photographs from turn of the century New York line the walls.

Original Bagel (Whole Wheat or Plain): $1.05
Bagel Thins (Whole Wheat, Plain or Everything): $1.05
Signature Bagels: $1.15
Small Challah Roll or Gourmet Bagel: $1.80
Bakers Dozens (13 bagels per dozen): 1/2 dozen$6.00 , regular dozen$10.00, dozen with two cream cheeses(shmears)$14.50. An individual container of shmear are $3.29 (plain), $3.49 (flavoured), or $3.99 (lox), however if you purchase bulk bagels (6 or more) the price drops to $2.50 per tub, no matter what flavour.
(prices current as of 5/2/12)

They have an extensive menu of toppings which includes: smoked salmon (1 oz piece of lox), turkey, ham, roast beef, pastrami, corned beef, capers, lettuce, red onion, tomatoes, sprouts, cheese, cucumbers, pickles, and many more.

They feature sweets baked daily: Blueberry muffins, Strawberry white chocolate muffins, LTO muffins, cinnamon walnut strudels, cinnamon twists, Heavenly Chocolate Chip Cookies, snickerdoodle cookies, mudslide cookies, and three flavors of coffee cakes (blueberry, applecinnamon, and chocolate chip).

Egg mits (turkey sausage, spinach and swiss, baconcheddar, plain egg, lox, or with just egg) are enormous breakfast sandwiches on a bagel of your choice. They are very filling and served all day. the egg also can be substituted for egg whites for the more healthconscientious among us. Their breakfast wrap, the Sante Fe, weighs over a pound!
They also offer soups, sandwiches, salads and pizza bagels which are incredible. They are made to order and also come with a pickle and a side of either fruit salad, chips, potato salad, or cole slaw talk about a complete meal (and very filling).

Noahs also has Catering Nosh Boxes. These are perfect for breakfast or lunch. You can get a variety of bagels and shmears, sandwiches, salads and pastries. They also deliver for a nominal fee.

They also sell a selection of coffees. (Though why one would choose to drink this coffee over Peets is beyond me.)  Users/JenKao |

My favorite is the group picture of a NY womens sport team. They just have this amused look on their faces, as if the photographer was having trouble getting them to stop laughing.  Users/JenKao )

I just went to Noahs and tried their Reuben ... yummy! SS

20060306 16:20:42 nbsp I really think that Noah’s has gone downhill in the last couple years. In the last 3 or 4 visits, I can’t think of one that has been completely accurate. The wrong creamcheese, forgetting the lox, not doing any toasting
 The staff is always friendly, but they don’t seem to receive nearly enough training, and there seems to be a high rate of turnover. I can’t count the number of different prices they’ve tried to charge us for a bagel with lox and cream cheese. They’ve recently introduced an expanded sandwich line, which seems to push the staff to their limits and create even longer lines. They have tasty bagels, but I think my fiancé and I are going to try other places for awhile, since we’re tired of being disappointed. Users/MattCzarnowski

20060306 16:30:33 nbsp Yeah Ive been disappointed by Noahs ever since they got rid of their four cheese grilled cheese bagel. Was one of the greatest sandwiches in the history of the sandwich. Users/ZacMorris

20060408 19:46:14 nbsp I love bagels!  Tony was so very helpful!   yay! Users/VeronicaVillegas.

20060420 02:10:54 nbsp This is a revelation, I didnt know that there was a Noahs in Davis. Excellent! Users/TusharRawat

20060616 12:05:16 nbsp I second the going downhill comment. When they did all their expanded sandwiches they took the smoked whitefish salad off the Menu. Le Boo. Bagels are decent, I like mine with lox creamcheese, tomato and onion. Avoid their Egg Mits They microwave the egg till its flavoless and has an awfull texture. Users/CharlesGeorge

20060928 17:23:10 nbsp Gone downhill? Wait, you hear that? Whaaaaaaaaa!!!! Oh that was the whambulence, haters. This place has tasted the same until high school, ITS BAGELS. Oh, the jalapeno salsa schmear, actually pretty yummy. Users/SteveBlais

20060928 17:24:35 nbsp Since high school, sorry Im tired stayed up late to study for a midterm. Users/SteveBlais

20061108 20:49:22 nbsp bagels are good. but if your have ever been to 52nd street bagel in so cal, its much better because theyre more flavorful Users/AnnThiNguyen
    You posted this two years ago, I realize, but still, who cares what the bagels in so cal are like seeing as they arent an option for anyone who reads this wiki.

    Hes using it as a reference just in case anyone who reads this wiki is from Southern California, such as myself. And yes, I do care, so shut your mouth, stop being such a smartass.
    Their goal wasnt to describe the quality of noahs bagels, it was to complain about bagels in Davis not being as wonderful as in So Cal. Thats not the intent of this particular wiki page. I wasnt being a smartass, i was questioning the relevance of their post. Also, shut my mouth..? wow, you really take bagel defense personally.
    Uhm, another question, I am from soCal which fifty second street do you mean? ~DP
    Wegmans Marco Polo bagels are also great.  Upstate New York and Central PA.  Wonderful.  Theres also a couple really good Jewish delis in South Florida that make really tasty ones.  Its a bitch of a morning drive from Davis, though.

20071215 17:07:17 nbsp The staff there is very friendly.  The service is quick; even when theres a line out the door, I dont wait for more than 10 minutes.
The spinach, egg and bacon panini is delicious Users/drpottsiv

20080117 21:52:35 nbsp yumm!! even the plain bagel with just plain cream cheese is awesome...never to sweet or too heavy like others... Users/SimplyIncognito

20081109 14:16:08 nbsp Noahs doesnt sell bagels, they sell bagel shaped bread.  There was a great SNR article a year ago, a blind taste test which Noahs almost unanimously lost (against local Sacramento area bagel shops).  Bagels are supposed to be boiled and then baked  dense and chewy.  Noahs differs greatly in method (not boiled! either steamed or steam treated or use steamovens or some such thing) to getem lighter and fluffier, and I think Ive read about using some sort of oil to help keepem from being dry.  Quite a difference! Im surprised no one commented about that over the years.  I ate here reasonably frequently when I lived in North and West Davis, but please people: be sure to try real bagels at other places! Im sure JW knows what Im talking about.  Users/EdWins

20090303 13:51:30 nbsp SOOOOOOOOOOOOOOOOOOOO I recently had a BAD customer service experience at Noahs.  READ THIS STORY AND BOYCOTT NOAHS FOR SHITTY CUSTOMER SERVICE!!!!!!!!!!!!!!!!

One of the employees handed me the wrong bagel, and i was holding it in my hand as I paid.  (Which is fine, ok, people make mistakes)...

When I got back to my office I noticed that it was the wrong bagel with smear i was allergic to.  So I had to go BACK to Noahs to get the RIGHT bagel, in the rain, got soaking wet, and it took about a half hour to drive back, find parking, wait for the bagel that they said would be ready when I got there but it wasnt... all the meanwhile Im on a deadline at work which has now passed... (Which is also fine, shit happens).... But all this....

...TO BE COMPLETELY INSULTED BY THE EMPLOYEE!

First, he accused me of grabbing the wrong bagel even though he HANDED it to me and saw me HOLDING the bagel while I was paying for it.  He basically called me stupid for not noticing it was not mine, even though it was in the bag already at that time (It SAID Vikki on it, not your name... not my fault!). If that wasnt RUDE ENOUGH!!!!!! he noticed my angst as I was waiting for the bagel and missing my work deadline and said dude, dont you work like, down the street? What an ass.  I NEVER make rude comments on wiki, but when it comes to treating people, especially CUSTOMERS with respect, that irks me.

I work in customer service, and you should always treat that person with at LEAST respect, if not like royalty.  GET A CLUE, NOAH EMPLOYEE! Never going there again.   Users/Squeekyclean
   Im not a big fan of Noahs either, but I think your anger is misdirected here. Boycott Noahs for shitty customer service? More like, Theres one Noahs employee who is an asshole. I know you probably wrote your comment so you could vent somewhere, but you should probably complain to the management and describe the employee or get his name. Now, if management does nothing about the problem, then you do have a valid grudge against that particular Noahs. Users/KellyM

20090304 19:11:34 nbsp No, Im upset at Noahs for not training their staff to be professional.  So, I still hold my comment to be valid.  Im very passionate about great customer service so bad service is a peeve. If Noahs does not train their staff on the importance of customer service, then they suck and people should not go there!   Users/Squeekyclean

20090606 06:58:46 nbsp I just got a bowl of chicken noodle soup and a bagel and the soup was BAD! I called them and they said I can get a refund or they would give me a fresh bowl (it wasnt even fresh, i saw them microwave it). I went back to get my refund, they didnt even apologize. NOTHING.  Users/AngelaJ

20090730 09:19:10 nbsp Ive now had two equally horrendous experiences at this joint and my wife and I both decided were never going back. Both times Ive gone, Ive had to wait 20 minutesplus for a bagel people, not a bagel with anything special  just a smear or whatever they call it. Also, both times, the person taking the order recognized someone in line behind me and made a Bline to take their order before mine. Finally, this last time, after waiting like I said for more than 20 minutes, I see a guy place an order on the counter, but not call anyones name. I assumed maybe the order wasnt finished. He was going to put more up there, or whatever, anyway, they call your name when your orders ready. This time, he placed the order on the counter, looked around sheepishly, and just walked back to do his slowpaced work. After it sat up there for five minutes getting cold, I walked up and saw that it was my order. To add insult to injury, the coffee tastes like slightly darkened water. No wonder they give you bottomless coffee. Youd have to drink an entire batch to get the buzz you get from a cup from anywhere else. I tried two different blends, thinking the first was just a fluke, and it was still horrible.

I think the only reason so many people go there is its the only breakfast joint on that end of Davis. BTW: I grew up in Philly and New York, and these things, I dont know what you call them, but bagel is a charitable word at best. More like chewybread with a hole in the middle. The only positive I give is the strawberry smear, which is delicious no doubt. I rarely leave such a negative comment but Im just so disappointed in a place that tries to come off as some New York delicatessen when they know and I know, this place has a long way to go to meet those standards.

In NYC, they know that smiling and making small talk is no substitute for good food and timely service. Users/condemned2bfree

20091009 10:00:55 nbsp I really dont get why every one is so upset with Noahs.  Its a freaking bagel chain not a three star michelin.  We dont complain about the quality of the meat at McDonalds because we know its crap and its supposed to be.  We dont cry about how bad the service is at Carls Junior because we know that they get paid nothing to do demeaning work.  I mean seriously... Im from New York and blah blah blah... They made we miss my work deadline blah blah blah... Bagels are supposed to be boiled not steamed blah... you get the point.  I guess I just want to clarify to those who are reading this and arent familiar with Noahs; ITS A CHAIN dont expect anything more than what you would get at any other corporate eatery.  If you want quality and kindness support small business.   Users/swiftling

20091119 23:38:29 nbsp Im addicted to their bagels with LOX and cream Cheese. OMG their lox is so amazing.  Users/BryceH

20091124 20:42:50 nbsp THE BEST COFFEE IN DAVIS!  I wish they would sell their Chelsea coffee beans to the public.   OHH and in regards to people saying they have bad service? ARE YOU KIDDING ME?  GO TO a Chinese Restaurant and see how they treat people there.Users/JSlice

20100129 16:15:14 nbsp Taste is taste, so if youve eaten here and liked it more power to you, however if you are wondering if they make true bagels by boiling then baking them the answer is no.  You can get a nice roll with a hole at Noahs, (it tastes just like the one youd get at Safeway, BTW) but if youve had actual bagels, youll be disappointed as I was. Users/rfrazier

20100326 11:05:55 nbsp Please bring back the sundried tomatoes bagel! Apparently those are discontinued. Its the best. What happened?   Users/LeeY

20100430 16:37:43 nbsp I wouldnt know where else to get my bagels. Filling enough to make it to lunch, plain or shmeared. Users/sjpacheco

20100430 17:51:47 nbsp I agree with LeeY, I always got sundried tomato bagels there but now they are discontinued. The sundried tomato cream cheese just doesnt do it for me (if they even have that anymore). Users/JenniferCook

20100502 09:42:57 nbsp I love Noahs Veg Out sandwich but they consistently mess up my order in one form or another. Just be sure to double check your order. Users/EmmyMelton

20100509 14:28:09 nbsp Visiting Salt Lake City, I went with my sister to Einstein Bagel and was surprised to see that they have glutenfree bagels! As I am a Celiac, I decided to try one. Very, very tasty! And since Einsteins and Noahs are the same, I am wondering if anyone has spotted glutenfree bagels at this location.  Users/nataliesadler

20101108 09:45:02 nbsp I was hungry,counting on a tasty bagles...but a very bad customer service...i ask for toasted bagle and....they didnt toast my bagle and also gave me a wrong bagle...fresh tomato was miss...looool..i have paid almost $8.00 for breakfast...very young team with no experience...i am not happy....  Users/amirhero

20101117 21:19:34 nbsp I have always loved Noahs, BUT was severely disappointed to discover last week (via Nutrition Info on Noahs website) that almost all, if not 100% of their bagels, have HIGH FRUCTOSE CORN SYRUP... buyers, beware!!  Im writing the company now to let them know; hoping enough people will do so and theyll offer bagels without it.   Users/lalamusic

20101230 05:07:35 nbsp Poor management. Staff is not provided the proper resources to do their jobs well, to do their jobs as is instructed by the training videos. Im under the impression that this has to do with minimizing costs. Employees are not paid well and cant really hope for raises, even after years of working there. Its a pretty terrible place to work and there are better places to spend your money in Davis... but if youre desperate for work, there really is a high turn over. Walk in with a big smile and a modest resume. If you cant speak with the manager, come back when shes there. And go back frequently. Users/Janny

20110112 10:20:53 nbsp My better half and I are regulars at Noahs. The food is usually very tasty and prepared rapidly and we do enjoy it. Today, however, we came across the first slight snag. We ordered sweet bagels (blueberry and cinnamon raisin) and, once we bit into them, we were greeted with a mouthful of onion and garlic vapor. The bagels certainly looked like blueberry and cinnamon raisin. We exchanged and the replacements (untreated  just plain) suffered the same problem. It could be that they are blindtesting a new cinnamongarliconion bagel and if so Im not a fan. I just hope they havent decided to bake sweet and savory together. Users/Mazzle

20110205 21:25:52 nbsp pumpernickle bagels with a smattering of cream cheese or jam paired with a cup of your favorite coffee will put you a bit closer toward eternal contentment.  Users/blastoff

20110514 17:12:20 nbsp It tastes good but in my opinion its too expensive. I pay like $5 for a bagel sandwich and eat it in about 2 minutes and im still hungry after. Users/rachealgalante

20110727 09:06:14 nbsp Like another poster already noted, this places uses HIGH FRUCTOSE CORN SYRUP in their bagels...pretty much defeating the purpose of being a healthy alternative to traditional fast foods restaurants! Users/JohnVanderMolen

20110728 16:02:28 nbsp Maybe I dont buy bagels enough to understand the process, but I feel this is a terrible inefficient food service business. I came in to pick up a single bagel for my wife, with nothing done to it, and I had to wait in a line (expected) and then was given my bagel. I then had to wait in another line, at the other end of the store, by the same person, who just moved to the other register... While I was paying they were taking time to make the 3 people in front of me their bigger orders.... I dont really understand the double line thing? Cant I just pay at the same time I give my order? Isnt that the logical way? All I know is I had to wait over 10 minutes to get a single bagel with nothing done to it, however, I got to hold the bagel (which I didnt own yet) within the first 2 minutes?!?! From now on, Ill just have my wife get her own bagels... Hear that Noahs? Youre destroying chivalry.  Users/elliottPollard
  I havent been to Noahs in awhile, but they used to have a line for those who are only getting bagels.  Admittedly, I think that most of those people were buying a half dozen or a dozen bagels, but Im sure the same rule would apply for just one bagel.  So, either they took that line away, or you just didnt see it.   Users/CovertProfessor
   Yeah, one of the other customers in line kinda glared at me because I wasnt behind the rest of the line. So moved to the back and figured they know best, because they go there more .... Users/elliottPollard
    Nah, the sandwich people dont really understand the bagel people.  But going to a bagel store just to buy bagels is an old and proud tradition that supersedes the bagelassandwich movement.  Next time, Id just ignore the glare.  cp

20110914 06:14:12 nbsp Noahs wasnt in Davis when I went to school there but I discovered it (and Peets) when I went to grad school at Cal. I instantly loved the bagels far more than the tougher traditional bagel. After finishing up I move to SoCal and it took some time but they finally have a shop in La Jolla (actually called Einsteins Bagels, same company same product) and I love it. I remember when there was only one Noahs and one Peets coffee and people would happily stand in line at both places.  Now, as a matter of fact I serve Noahs/Einsteins bagels and Peets coffee at our weekly lab meetings and everyone is now addicted! Users/PaulWilliamDiaz

20111001 01:27:03 nbsp I think that some of the feedback posted is rather inaccurate. First, Noahs and Einsteins are not the same company. They were however started by the same person who sold Noahs and ten started Einsteins. You would think that creating one bagel business that doesnt actually make traditional bagels would be enough for anyone. Real bagels are boiled and then baked. They are super soft on the inside and a little bit chewy on the outside. They probably dont maintain there texture as long as the fake kind. If you think it doesnt matter than try microwaving your bread rather than baking it just to see how big a difference process makes. Might as well roll up wonder bread and roll it into a circle and call that a bagel or cut a round pier of cake cut a hole into it and call it a donut. Lies are lies. But if you like it then you are no different then the many who think that high fructose corn sweetener tastes. And there are upsides to having unrefined tastebuds and lower bills! Users/RudiJ

20120403 11:43:53 nbsp This place is tasty but I cant figure out their prices. Love their garlicherb spread, but why does it cost an extra $1.50 to get spread on a bagel, and why is that not noted on the menu? Also, their potato salad is delicious (I normally dont even like potato salad), but Im always charged a different amount when I get it, from $1 to $2.50. I would also recommend they get more vegetarian sandwich options. Users/MeggoWaffle

20120505 08:38:39 nbsp You know, the food doesnt taste terrible. Not worth a long line, but decent. However, my wife likes toasted bagels with cream cheese, lox, and capers, and they cant do it right. They charge for the capers, and dont put them on the bagel. This has happened 4 out of 5 times over the past yearandahalf (I learned my lesson, and babysit the employee while they make the bagel, and check it before I leave). Users/EricWecker

20121218 08:38:01 nbsp The best breakfast bagel sandwiches in Davis! I see to it that I come here whenever I come back home!! Users/EmThornhill

20140803 12:22:00 nbsp The last three times I have been to Noahs I have waited over 30 minutes for small orders. I worked at a similar bagel store for years and know what efficiency looks like. This place needs to figure out a better system. The employees are nice but much to slow and reliant on the order displays. Also, if you are going to have an order here sign then follow it. Nobody knew where to go.  Anonymous


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Computes norms of error in test mode case.
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0
!> @param[in] nnod_loc number of local nodes
!> @param[in] u1 current displacement
!> @param[in] v1 current velocity
!> @param[in] time discrete instant time
!> @param[in] ne_loc number of local elements
!> @param[in] cs_loc spectral local connectivity
!> @param[in] cs_nnz_loc length of cs_loc
!> @param[in] nm number of materials
!> @param[in] prop_mat material properties (rho, lambda, mu, gamma)
!> @param[in] sdeg_mat polynomial degree vector
!> @param[in] tag_mat label for material properties
!> @param[in] loc_n_num  local node numbering (local to global)
!> @param[in] local_el_num local element numbering (local to global)
!> @param[in] nelem_dg_local number of DG elements
!> @param[in] alfa11 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa12 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa13 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa21 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa22 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa23 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa31 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa32 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] alfa33 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta11 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta12 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta13 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta21 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta22 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta23 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta31 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta32 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] beta33 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] gamma1 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] gamma2 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] gamma3 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] delta1 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] delta2 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] delta3 costant value for the bilinear map from (-1,1)^3 to the current element
!> @param[in] xs_loc x-coordinate spectral node
!> @param[in] ys_loc y-coordinate spectral node
!> @param[in] zs_loc z-coordinate spectral node
!> @param[in] mpi_id MPI id for process
!> @param[in] el_new struct for DG elements
!> @param[in] nelem_dg_global number of total DG elemnents
!> @param[in] nsd_jump number of node to send for jumps
!> @param[in] node_sd_jump list of node to send for jums
!> @param[in] nrv_jump number of node to receive for jumps
!> @param[in] node_rv_jump list of node to receive for jums
!> @param[in] send_length_jump number of node to send for jumps to each MPI process
!> @param[in] recv_length_jump number of node to receive for jumps to each MPI process
!> @param[in] mpi_np  number of MPI process
!> @param[in] mpi_comm MPI common world
!> @param[in] cs_nnz_dg length of cs_dg
!> @param[in] cs_dg spectral connectivity for DG elements
!> @param[out] ... Write EN.ERR file containing the computed error norms.

!> @warning Only for serial computation. Need to be tested for parallel computation

     subroutine COMPUTE_ENERGY_ERROR(nnod_loc, u1, v1, time, ne_loc, cs_loc, cs_nnz_loc,&
                              nm,prop_mat, sdeg_mat,tag_mat, &
                              loc_n_num, local_el_num, nelem_dg_local, &
                              alfa11,alfa12,alfa13,&
                              alfa21,alfa22,alfa23,&
                              alfa31,alfa32,alfa33,&
                              beta11,beta12,beta13,&
                              beta21,beta22,beta23,&
                              beta31,beta32,beta33,&
                              gamma1,gamma2,gamma3,&
                              delta1,delta2,delta3,&
                              xs_loc,ys_loc,zs_loc,&
                              mpi_id,el_new, &
                              nelem_dg_global, &
                              nsd_jump,node_sd_jump, &
                              nrv_jump,node_rv_jump, &
                              send_length_jump, recv_length_jump, &
                              mpi_np, mpi_comm, &
                              cs_nnz_dg, cs_dg)

     use max_var
     use DGJUMP
     use speed_exit_codes

     implicit none

     include 'SPEED.MPI'

     type(el4loop), dimension(ne_loc), intent(in) :: el_new

     character*100000 :: input_line
     character*4 :: keyword
     character*70 :: file_mpi

     integer*4 :: nnod_loc, nfunc, mpi_id, cs_nnz_loc, nelem_dg_global, nsd_jump, nrv_jump
     integer*4 :: nm, unit_mpi, ne_loc, status, ic, id, ie_curr, jb,kb,ib,isb
     integer*4 :: mpi_np, mpi_comm, mpi_ierr, trofa, posit, cs_nnz_dg
     integer*4 :: ielem, iene, iene_curr, imne, jshift,mm, nelem_dg_local
     integer*4 :: ie, im, nn, is,in,iaz,i,j,k,ileft,iright,nval
     integer*4, dimension(0:cs_nnz_loc) :: cs_loc
     integer*4, dimension(0:cs_nnz_dg) ::  cs_dg
     integer*4, dimension(nnod_loc) :: loc_n_num
     integer*4, dimension(ne_loc) :: local_el_num
     integer*4, dimension(nsd_jump):: node_sd_jump
     integer*4, dimension(nrv_jump):: node_rv_jump
     integer*4, dimension(nm) :: sdeg_mat
     integer*4, dimension(nm) :: tag_mat
     integer*4, dimension(SPEED_STATUS_SIZE) :: mpi_stat
     integer*4, dimension(mpi_np) :: send_length_jump, recv_length_jump

     real*8 :: time, term1,term2,term3, L2_err, EN_err, term_ener,term_l2
     real*8 :: Linf_err, EN_err_tot, L2_err_tot, Linf_err_tot, xp, yp, zp, phi

     real*8 :: L2_err_vel, H1_err, term_H1, H1_err_tot, term_l2_vel
     real*8 :: Linf_err_vel, L2_err_vel_tot, Linf_err_vel_tot
     real*8 :: e11_ex, e12_ex, e13_ex, e22_ex, e23_ex, e33_ex
     real*8 :: sxx_ex, syy_ex, szz_ex, sxy_ex, szx_ex, syz_ex

     real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3,w1, rho,lambda,mu, cost1,cost2,cost3,pi,x,y,z, u1_ctm, u2_ctm, u3_ctm
     real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j, u1_ex,u2_ex,u3_ex, v1_ex, v2_ex, v3_ex

     real*8, dimension(4) :: Linf_array, Linf_array_vel
     real*8, dimension(nm,4) :: prop_mat
     real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
     real*8, dimension(ne_loc) :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
     real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3,delta1,delta2,delta3
     real*8, dimension(3*nnod_loc) :: u1, v1
     real*8, dimension(nnod_loc) :: xs_loc, ys_loc, zs_loc
     real*8, dimension(3*nsd_jump) :: send_buffer_error
     real*8, dimension(3*nrv_jump) :: recv_buffer_error, jump_minus_error

     real*8, dimension(:), allocatable :: ct, ww, ctm, wwm, u_p, u_m, jump_all
     real*8, dimension(:,:), allocatable :: dd, ddm
     real*8, dimension(:), allocatable :: val_data, basis_functions
     real*8, dimension(:,:,:), allocatable :: ux_el,uy_el,uz_el
     real*8, dimension(:,:,:), allocatable :: duxdx_el,duydx_el,duzdx_el
     real*8, dimension(:,:,:), allocatable :: duxdy_el,duydy_el,duzdy_el
     real*8, dimension(:,:,:), allocatable :: duxdz_el,duydz_el,duzdz_el
     real*8, dimension(:,:,:), allocatable :: sxx_el,syy_el,szz_el
     real*8, dimension(:,:,:), allocatable :: syz_el,szx_el,sxy_el
     real*8, dimension(:,:,:), allocatable :: lambda_el,mu_el

     pi = 4.d0*datan(1.d0)

     L2_err = 0.d0; EN_err = 0.d0; Linf_err = 0.d0;
     L2_err_tot = 0.d0; EN_err_tot = 0.d0; Linf_err_tot = 0.d0;

     L2_err_vel = 0.d0;  H1_err = 0.d0;  Linf_err_vel = 0.d0;
     L2_err_vel_tot = 0.d0;  H1_err_tot = 0.d0; Linf_err_vel_tot = 0.d0;

     do ie = 1,ne_loc

         im = cs_loc(cs_loc(ie -1))
         nn = sdeg_mat(im) +1;
         mm = nn ! 8 ! gll nodes for computing norms
         allocate(ct(nn),ww(nn),dd(nn,nn));  allocate(ctm(mm),wwm(mm),ddm(mm,mm));

         allocate(ux_el(nn,nn,nn),uy_el(nn,nn,nn),uz_el(nn,nn,nn))
     	 allocate(duxdx_el(nn,nn,nn),duydx_el(nn,nn,nn),duzdx_el(nn,nn,nn))
    	 allocate(duxdy_el(nn,nn,nn),duydy_el(nn,nn,nn),duzdy_el(nn,nn,nn))
         allocate(duxdz_el(nn,nn,nn),duydz_el(nn,nn,nn),duzdz_el(nn,nn,nn))
     	 allocate(sxx_el(nn,nn,nn),syy_el(nn,nn,nn),szz_el(nn,nn,nn))
      	 allocate(syz_el(nn,nn,nn),szx_el(nn,nn,nn),sxy_el(nn,nn,nn))
      	 allocate(lambda_el(nn,nn,nn),mu_el(nn,nn,nn))

         call MAKE_LGL_NW(nn,ct,ww,dd)
         call MAKE_LGL_NW(mm,ctm,wwm,ddm)

         rho = prop_mat(im,1); lambda = prop_mat(im,2); mu = prop_mat(im,3)

         do k = 1,nn
            do j = 1,nn
               do i = 1,nn
                  is = nn*nn*(k -1) +nn*(j -1) +i
                  in = cs_loc(cs_loc(ie -1) + is)

                  iaz = 3*(in -1) +1; ux_el(i,j,k) = u1(iaz)
                  iaz = 3*(in -1) +2; uy_el(i,j,k) = u1(iaz)
                  iaz = 3*(in -1) +3; uz_el(i,j,k) = u1(iaz)

               enddo
            enddo
         enddo

        call MAKE_STRAIN_TENSOR(nn,ct,ww,dd,&
        		  alfa11(ie),alfa12(ie),alfa13(ie),&
                          alfa21(ie),alfa22(ie),alfa23(ie),&
                     	  alfa31(ie),alfa32(ie),alfa33(ie),&
                     	  beta11(ie),beta12(ie),beta13(ie),&
                     	  beta21(ie),beta22(ie),beta23(ie),&
                     	  beta31(ie),beta32(ie),beta33(ie),&
                     	  gamma1(ie),gamma2(ie),gamma3(ie),&
                     	  delta1(ie),delta2(ie),delta3(ie),&
                     	  ux_el,uy_el,uz_el,&
                     	  duxdx_el,duydx_el,duzdx_el,&
		     	  duxdy_el,duydy_el,duzdy_el,&
		     	  duxdz_el,duydz_el,duzdz_el)

        lambda_el = lambda; mu_el = mu;

        call MAKE_STRESS_TENSOR(nn,lambda_el,mu_el,&
                                duxdx_el,duydx_el,duzdx_el,&
	                        duxdy_el,duydy_el,duzdy_el,&
				duxdz_el,duydz_el,duzdz_el,&
                                sxx_el,syy_el,szz_el,&
				syz_el,szx_el,sxy_el)

         do k = 1,mm
            do j = 1,mm
               do i = 1,mm
                      dxdx = alfa11(ie) +beta12(ie)*ctm(k) &
                             + beta13(ie)*ctm(j) +gamma1(ie)*ctm(j)*ctm(k)
                      dydx = alfa21(ie) +beta22(ie)*ctm(k) &
                             + beta23(ie)*ctm(j) +gamma2(ie)*ctm(j)*ctm(k)
                      dzdx = alfa31(ie) +beta32(ie)*ctm(k) &
                             + beta33(ie)*ctm(j) +gamma3(ie)*ctm(j)*ctm(k)

                      dxdy = alfa12(ie) +beta11(ie)*ctm(k) &
                             + beta13(ie)*ctm(i) +gamma1(ie)*ctm(k)*ctm(i)
                      dydy = alfa22(ie) +beta21(ie)*ctm(k) &
                             + beta23(ie)*ctm(i) +gamma2(ie)*ctm(k)*ctm(i)
                      dzdy = alfa32(ie) +beta31(ie)*ctm(k) &
                             + beta33(ie)*ctm(i) +gamma3(ie)*ctm(k)*ctm(i)

                      dxdz = alfa13(ie) +beta11(ie)*ctm(j) &
                           + beta12(ie)*ctm(i) +gamma1(ie)*ctm(i)*ctm(j)
                      dydz = alfa23(ie) +beta21(ie)*ctm(j) &
                           + beta22(ie)*ctm(i) +gamma2(ie)*ctm(i)*ctm(j)
                      dzdz = alfa33(ie) +beta31(ie)*ctm(j) &
                           + beta32(ie)*ctm(i) +gamma3(ie)*ctm(i)*ctm(j)

                      det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
                            - dydz * (dxdx*dzdy - dzdx*dxdy) &
                            + dzdz * (dxdx*dydy - dydx*dxdy)

                      !Local to global node transformation
                      xp = alfa11(ie)*ctm(i) + alfa12(ie)*ctm(j) &
                           + alfa13(ie)*ctm(k) + beta11(ie)*ctm(j)*ctm(k) &
                           + beta12(ie)*ctm(i)*ctm(k) + beta13(ie)*ctm(i)*ctm(j) &
                           + gamma1(ie)*ctm(i)*ctm(j)*ctm(k) + delta1(ie)

                      yp = alfa21(ie)*ctm(i) + alfa22(ie)*ctm(j) &
                           + alfa23(ie)*ctm(k) + beta21(ie)*ctm(j)*ctm(k) &
                           + beta22(ie)*ctm(i)*ctm(k) + beta23(ie)*ctm(i)*ctm(j) &
                           + gamma2(ie)*ctm(i)*ctm(j)*ctm(k) + delta2(ie)

                      zp = alfa31(ie)*ctm(i) + alfa32(ie)*ctm(j) &
                           + alfa33(ie)*ctm(k) + beta31(ie)*ctm(j)*ctm(k) &
                           + beta32(ie)*ctm(i)*ctm(k) + beta33(ie)*ctm(i)*ctm(j) &
                           + gamma3(ie)*ctm(i)*ctm(j)*ctm(k) + delta3(ie)

                      !STATIONARY POLYNOMIAL SOLUTIONS
                      !u1_ex = dsin(3.d0*pi*time)*((xp*(xp-1.d0))*(yp*(yp-1.d0))*(zp*(zp-1.d0)));
                      !u2_ex = 0.d0;
                      !u3_ex = 0.d0;

                      is = nn*nn*(k -1) +nn*(j -1) +i
                      in = cs_loc(cs_loc(ie -1) +is)

      !POPLYNOMIAL SOLUTION
      !u1_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
      !u2_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
      !u3_ex = dsin(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
      !
      !v1_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
      !v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);
      !v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*((xs_loc(in)*(xs_loc(in)-1.d0))**2.d0 &
      !                                      *(ys_loc(in)*(ys_loc(in)-1.d0))**2.d0 &
      !                                      *(zs_loc(in)*(zs_loc(in)-1.d0))**2.d0);

      !PAPER WITH BLANCA
      u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xp))**2.d0*dsin(2.d0*pi*yp)*dsin(2.d0*pi*zp);
      u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*yp))**2.d0*dsin(2.d0*pi*xp)*dsin(2.d0*pi*zp);
      u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zp))**2.d0*dsin(2.d0*pi*xp)*dsin(2.d0*pi*yp);

      v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));

      x = xs_loc(in); y = ys_loc(in); z = zs_loc(in);
      e11_ex = -2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z)
      e22_ex =  2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z)

      e33_ex =  2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)

      e12_ex =  pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(pi*y)**2.d0*dsin(2.d0*pi*z) &
                              - pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*z)

      e13_ex = pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(2.d0*pi*y)*dsin(pi*z)**2.d0 &
                             - pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*y)

      e23_ex = pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)**2.d0 &
                             + pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*z)**2.d0

      sxx_ex = lambda*(2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
                               - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
                               + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
                               - 4.d0*pi*mu*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z)

      syy_ex = lambda*(2.d0*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
                               - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
                               + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
                               + 4.d0*pi*mu*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z)

      szz_ex = lambda*(2*pi*dcos(pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)*dsin(2.d0*pi*z) &
                               - 2.d0*pi*dcos(pi*x)*dsin(3.d0*pi*time)*dsin(pi*x)*dsin(2.d0*pi*y)*dsin(2.d0*pi*z) &
                               + 2.d0*pi*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)) &
                               + 4.d0*pi*mu*dcos(pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(2.d0*pi*y)*dsin(pi*z)

      sxy_ex = mu*(2.d0*pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(pi*y)**2.d0*dsin(2.d0*pi*z) &
                               - 2.d0*pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*z))

      szx_ex = mu*(2.d0*pi*dcos(2.d0*pi*x)*dsin(3.d0*pi*time)*dsin(2.d0*pi*y)*dsin(pi*z)**2.d0 &
                               - 2.d0*pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(pi*x)**2.d0*dsin(2.d0*pi*y))

      syz_ex = mu*(2.d0*pi*dcos(2.d0*pi*z)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*y)**2.d0 &
                               + 2.d0*pi*dcos(2.d0*pi*y)*dsin(3.d0*pi*time)*dsin(2.d0*pi*x)*dsin(pi*z)**2.d0)

                      u1_ctm = 0.d0; u2_ctm = 0.d0; u3_ctm = 0.d0

                      allocate(basis_functions(mm**3)) ! basis function evaluation
                      do kb = 1, nn
                         do jb = 1, nn
                            do ib =1, nn
                               isb = nn*nn*(kb -1) + nn*(jb -1) +ib
                               in = cs_loc(cs_loc(ie -1) +isb)
                               iaz = 3*(in -1) + 1

                               basis_functions(isb) = phi(nn-1, ct(ib), ctm(i)) * &
                                                      phi(nn-1, ct(jb), ctm(j)) * &
                                                      phi(nn-1, ct(kb), ctm(k))

                               u1_ctm =  u1_ctm +  basis_functions(isb)*u1(iaz);
                               u2_ctm =  u2_ctm +  basis_functions(isb)*u1(iaz+1);
                               u3_ctm =  u3_ctm +  basis_functions(isb)*u1(iaz+2);

                            enddo
                         enddo
                      enddo
                      deallocate(basis_functions)

                      !L2 NORM OLD
                      is = nn*nn*(k -1) +nn*(j -1) +i
                      in = cs_loc(cs_loc(ie -1) +is)

                      iaz = 3*(in -1) + 1

                      term_l2_vel =  det_j * ww(i) * ww(j) * ww(k) * &
                        ( (v1_ex - v1(iaz))**2.d0 + (v2_ex - v1(iaz+1))**2.d0 +(v3_ex - v1(iaz+2))**2.d0)

                      L2_err_vel = L2_err_vel + term_l2_vel

                      !LINF NORM - DISPL & VEL
                      Linf_array = (/ Linf_err,abs(u1_ex - u1_ctm),abs(u2_ex - u2_ctm),abs(u3_ex - u3_ctm) /)
                      Linf_err = maxval(Linf_array)

                      !L2 NORM NEW
                      term_l2 =  det_j * wwm(i) * wwm(j) * wwm(k) * &
                        ( (u1_ex - u1_ctm)**2.d0 + (u2_ex - u2_ctm)**2.d0 +(u3_ex - u3_ctm)**2.d0)
                      L2_err = L2_err + term_l2

                      !ENERGY NORM
                      iaz = 3*(in -1) + 1
                      term_ener =  det_j * ww(i) * ww(j) * ww(k) * ( &
                                   rho * (v1_ex - v1(iaz+0))**2.d0 &
                                 + rho * (v2_ex - v1(iaz+1))**2.d0 &
                                 + rho * (v3_ex - v1(iaz+2))**2.d0 &
                                 + (sxx_el(i,j,k) - sxx_ex) * (duxdx_el(i,j,k) - e11_ex) &
                                 + (syy_el(i,j,k) - syy_ex) * (duydy_el(i,j,k) - e22_ex) &
                                 + (szz_el(i,j,k) - szz_ex) * (duzdz_el(i,j,k) - e33_ex) &
                                 + 2.d0*(sxy_el(i,j,k) - sxy_ex) * (0.5d0*duxdy_el(i,j,k) + 0.5d0*duydx_el(i,j,k) - e12_ex) &
                                 + 2.d0*(szx_el(i,j,k) - szx_ex) * (0.5d0*duxdz_el(i,j,k) + 0.5d0*duzdx_el(i,j,k) - e13_ex) &
                                 + 2.d0*(syz_el(i,j,k) - syz_ex) * (0.5d0*duydz_el(i,j,k) + 0.5d0*duzdy_el(i,j,k) - e23_ex) &
                                )

                      EN_err = EN_err + term_ener

                      term_H1 =  det_j * ww(i) * ww(j) * ww(k) * ( &
                                 + (sxx_el(i,j,k) - sxx_ex) * (duxdx_el(i,j,k) - e11_ex) &
                                 + (syy_el(i,j,k) - syy_ex) * (duydy_el(i,j,k) - e22_ex) &
                                 + (szz_el(i,j,k) - szz_ex) * (duzdz_el(i,j,k) - e33_ex) &
                                 + 2.d0*(sxy_el(i,j,k) - sxy_ex) * (0.5d0*duxdy_el(i,j,k) + 0.5d0*duydx_el(i,j,k) - e12_ex) &
                                 + 2.d0*(szx_el(i,j,k) - szx_ex) * (0.5d0*duxdz_el(i,j,k) + 0.5d0*duzdx_el(i,j,k) - e13_ex) &
                                 + 2.d0*(syz_el(i,j,k) - syz_ex) * (0.5d0*duydz_el(i,j,k) + 0.5d0*duzdy_el(i,j,k) - e23_ex) &
                                )

                      H1_err = H1_err + term_H1

                      enddo
                   enddo
                enddo

         deallocate(ct,ww,dd, ctm, wwm, ddm)
         deallocate(ux_el,uy_el,uz_el)
     	 deallocate(duxdx_el,duydx_el,duzdx_el,duxdy_el,duydy_el,duzdy_el,duxdz_el,duydz_el,duzdz_el)
     	 deallocate(sxx_el,syy_el,szz_el,syz_el,szx_el,sxy_el)
         deallocate(lambda_el,mu_el)

     enddo

!----------------------------------------------------------------------------------------------------------
!    NOW COMPUTE NON CONFORMING TERM S_e int_e [u-uh][u-uh]
!----------------------------------------------------------------------------------------------------------
     if(nelem_dg_global .gt. 0) then

          jump_minus_error = 0.d0

          do i = 1,nsd_jump
             in = node_sd_jump(i)
                  !PAPER WITH BLANCA
                  u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
                  u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
                  u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));

                  send_buffer_error(3*(i -1) +1) = abs(u1(3*(in -1) +1) - u1_ex)
                  send_buffer_error(3*(i -1) +2) = abs(u1(3*(in -1) +2) - u2_ex)
                  send_buffer_error(3*(i -1) +3) = abs(u1(3*(in -1) +3) - u3_ex)

          enddo

         call EXCHANGE_DOUBLE(3*nsd_jump, send_buffer_error, 3*nrv_jump, recv_buffer_error,&
                                mpi_np, send_length_jump, recv_length_jump,&
                                mpi_comm, mpi_stat, mpi_ierr,mpi_id)

          do i = 1,nrv_jump
             in = node_rv_jump(i)
             jump_minus_error(3*(i -1) +1) = recv_buffer_error(3*(in -1) +1)
             jump_minus_error(3*(i -1) +2) = recv_buffer_error(3*(in -1) +2)
             jump_minus_error(3*(i -1) +3) = recv_buffer_error(3*(in -1) +3)
          enddo

      endif

      do ie = 1, nelem_dg_local

        ielem = el_new(ie)%ind; nn = el_new(ie)%deg
        call GET_INDLOC_FROM_INDGLO(local_el_num, ne_loc, ielem, ie_curr)

        allocate(u_p(3*nn**3)); u_p = 0.d0;

               do k = 1,nn
                  do j = 1,nn
                     do i = 1,nn
                        is = nn*nn*(k -1) +nn*(j -1) +i
                        in = cs_loc(cs_loc(ie_curr -1) + is)

                  !PAPER WITH BLANCA
                  u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
                  u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
                  u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));

                        iaz = 3*(in -1) +1; u_p(is) = abs(u1(iaz) - u1_ex)
                        iaz = 3*(in -1) +2; u_p(nn**3+is) = abs(u1(iaz) - u2_ex)
                        iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = abs(u1(iaz) - u3_ex)

                     enddo
                  enddo
                enddo

               allocate(u_m(el_new(ie)%nnz_col_only_uv)); u_m = 0.d0; jshift = 0

               do ic = 1, el_new(ie)%num_of_ne

                   iene = el_new(ie)%el_conf(ic,1)
                   imne = el_new(ie)%el_conf(ic,0)

                   mm = 2
                   do i = 1, nm
                      if(tag_mat(i) .eq. imne ) then
                         mm = sdeg_mat(i) +1
                      endif
                   enddo

                   trofa  = 0
                   if(mpi_np .gt. 1) call CHECK_MPI(cs_nnz_dg, cs_dg, iene, trofa, posit)

                   if(trofa .eq. 1) then
                     do k = 1,mm
                        do j = 1,mm
                           do i = 1,mm
                              is =  mm*mm*(k -1) +mm*(j -1) +i
                              id = cs_dg(cs_dg(posit-1)+is)
                              iaz = 3*(id -1) +1; u_m(is + jshift) = jump_minus_error(iaz)
                              iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = jump_minus_error(iaz)
                              iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = jump_minus_error(iaz)
                            enddo
                         enddo
                      enddo
                   else

                     call GET_INDLOC_FROM_INDGLO(local_el_num, ne_loc, iene, iene_curr)

                     do k = 1,mm
                        do j = 1,mm
                           do i = 1,mm
                              is =  mm*mm*(k -1) +mm*(j -1) +i
                              id = cs_loc(cs_loc(iene_curr-1)+is)

                  !PAPER WITH BLANCA
                  u1_ex = -dsin(3.d0*pi*time)*(dsin(pi*xs_loc(id)))**2.d0*dsin(2.d0*pi*ys_loc(id))*dsin(2.d0*pi*zs_loc(id));
                  u2_ex = +dsin(3.d0*pi*time)*(dsin(pi*ys_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*zs_loc(id));
                  u3_ex = +dsin(3.d0*pi*time)*(dsin(pi*zs_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*ys_loc(id));

                              iaz = 3*(id -1) +1; u_m(is + jshift) = abs(u1(iaz) - u1_ex)
                              iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = abs(u1(iaz) - u2_ex)
                              iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = abs(u1(iaz) - u3_ex)

                           enddo
                        enddo
                      enddo
                    endif

                      !  write(*,*) '--------------- fine +- ------------------------------'
                      !  read(*,*)

                   jshift = jshift + 3*mm**3

                enddo

                allocate(jump_all(3*nn**3))
                jump_all = 0.d0

                call MATMUL_SPARSE(el_new(ie)%matMin_only_uv, el_new(ie)%nnz_minus_only_uv, el_new(ie)%JMin_only_uv, &
                                  el_new(ie)%IMin_only_uv, jump_all, 3*nn**3,  &
                                  u_m, el_new(ie)%nnz_col_only_uv,1)

                EN_err = EN_err + dot_product(jump_all,u_p)
                H1_err = H1_err + dot_product(jump_all,u_p)

                deallocate(u_m, jump_all)
                allocate(jump_all(3*nn**3)); jump_all = 0.d0

                call MATMUL_SPARSE(el_new(ie)%matPlus_only_uv, el_new(ie)%nnz_plus_only_uv,el_new(ie)%JPlus_only_uv, &
                                  el_new(ie)%IPlus_only_uv, jump_all, 3*nn**3, &
                                  u_p, 3*nn**3,1)

                EN_err = EN_err + dot_product(jump_all,u_p)
                H1_err = H1_err + dot_product(jump_all,u_p)

                 deallocate(u_p, jump_all)

      enddo

!----------------------------------------------------------------------------------------------------------
!    NOW COMPUTE NON CONFORMING TERM S_e int_e [v-vh][v-vh] v=du/dt
!----------------------------------------------------------------------------------------------------------
     if(nelem_dg_global .gt. 0) then

          jump_minus_error = 0.d0

          do i = 1,nsd_jump
             in = node_sd_jump(i)
                  !PAPER WITH BLANCA
      v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));

                  send_buffer_error(3*(i -1) +1) = abs(v1(3*(in -1) +1) - v1_ex)
                  send_buffer_error(3*(i -1) +2) = abs(v1(3*(in -1) +2) - v2_ex)
                  send_buffer_error(3*(i -1) +3) = abs(v1(3*(in -1) +3) - v3_ex)

          enddo

         call EXCHANGE_DOUBLE(3*nsd_jump, send_buffer_error, 3*nrv_jump, recv_buffer_error,&
                                mpi_np, send_length_jump, recv_length_jump,&
                                mpi_comm, mpi_stat, mpi_ierr,mpi_id)

          do i = 1,nrv_jump
             in = node_rv_jump(i)
             jump_minus_error(3*(i -1) +1) = recv_buffer_error(3*(in -1) +1)
             jump_minus_error(3*(i -1) +2) = recv_buffer_error(3*(in -1) +2)
             jump_minus_error(3*(i -1) +3) = recv_buffer_error(3*(in -1) +3)
          enddo

      endif

      do ie = 1, nelem_dg_local

        ielem = el_new(ie)%ind; nn = el_new(ie)%deg
        call GET_INDLOC_FROM_INDGLO(local_el_num, ne_loc, ielem, ie_curr)

        allocate(u_p(3*nn**3)); u_p = 0.d0;

               do k = 1,nn
                  do j = 1,nn
                     do i = 1,nn
                        is = nn*nn*(k -1) +nn*(j -1) +i
                        in = cs_loc(cs_loc(ie_curr -1) + is)

                  !PAPER WITH BLANCA
      v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(in)))**2.d0*dsin(2.d0*pi*ys_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*zs_loc(in));
      v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(in)))**2.d0*dsin(2.d0*pi*xs_loc(in))*dsin(2.d0*pi*ys_loc(in));

                        iaz = 3*(in -1) +1; u_p(is) = abs(v1(iaz) - v1_ex)
                        iaz = 3*(in -1) +2; u_p(nn**3+is) = abs(v1(iaz) - v2_ex)
                        iaz = 3*(in -1) +3; u_p(2*nn**3 + is) = abs(v1(iaz) - v3_ex)

                     enddo
                  enddo
                enddo

               allocate(u_m(el_new(ie)%nnz_col_only_uv)); u_m = 0.d0; jshift = 0

               do ic = 1, el_new(ie)%num_of_ne

                   iene = el_new(ie)%el_conf(ic,1)
                   imne = el_new(ie)%el_conf(ic,0)

                   mm = 2
                   do i = 1, nm
                      if(tag_mat(i) .eq. imne ) then
                         mm = sdeg_mat(i) +1
                      endif
                   enddo

                   trofa  = 0
                   if(mpi_np .gt. 1) call CHECK_MPI(cs_nnz_dg, cs_dg, iene, trofa, posit)

                   if(trofa .eq. 1) then
                     do k = 1,mm
                        do j = 1,mm
                           do i = 1,mm
                              is =  mm*mm*(k -1) +mm*(j -1) +i
                              id = cs_dg(cs_dg(posit-1)+is)
                              iaz = 3*(id -1) +1; u_m(is + jshift) = jump_minus_error(iaz)
                              iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = jump_minus_error(iaz)
                              iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = jump_minus_error(iaz)
                            enddo
                         enddo
                      enddo
                   else

                     call GET_INDLOC_FROM_INDGLO(local_el_num, ne_loc, iene, iene_curr)

                     do k = 1,mm
                        do j = 1,mm
                           do i = 1,mm
                              is =  mm*mm*(k -1) +mm*(j -1) +i
                              id = cs_loc(cs_loc(iene_curr-1)+is)

                  !PAPER WITH BLANCA
      v1_ex = - 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*xs_loc(id)))**2.d0*dsin(2.d0*pi*ys_loc(id))*dsin(2.d0*pi*zs_loc(id));
      v2_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*ys_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*zs_loc(id));
      v3_ex = + 3.d0*pi*dcos(3.d0*pi*time)*(dsin(pi*zs_loc(id)))**2.d0*dsin(2.d0*pi*xs_loc(id))*dsin(2.d0*pi*ys_loc(id));

                              iaz = 3*(id -1) +1; u_m(is + jshift) = abs(v1(iaz) - v1_ex)
                              iaz = 3*(id -1) +2; u_m(mm**3 + is + jshift) = abs(v1(iaz) - v2_ex)
                              iaz = 3*(id -1) +3; u_m(2*(mm**3) + is + jshift) = abs(v1(iaz) - v3_ex)

                           enddo
                        enddo
                      enddo
                    endif

                      !  write(*,*) '--------------- fine +- ------------------------------'
                      !  read(*,*)

                   jshift = jshift + 3*mm**3

                enddo

                allocate(jump_all(3*nn**3))
                jump_all = 0.d0

                call MATMUL_SPARSE(el_new(ie)%matMin_only_uv, el_new(ie)%nnz_minus_only_uv, el_new(ie)%JMin_only_uv, &
                                  el_new(ie)%IMin_only_uv, jump_all, 3*nn**3,  &
                                  u_m, el_new(ie)%nnz_col_only_uv,1)

                EN_err = EN_err + dot_product(jump_all,u_p)
                !write(*,*) dot_product(jump_all,u_p)
                !H1_err = H1_err + dot_product(jump_all,u_p)

                deallocate(u_m, jump_all)
                allocate(jump_all(3*nn**3)); jump_all = 0.d0

                call MATMUL_SPARSE(el_new(ie)%matPlus_only_uv, el_new(ie)%nnz_plus_only_uv,el_new(ie)%JPlus_only_uv, &
                                  el_new(ie)%IPlus_only_uv, jump_all, 3*nn**3, &
                                  u_p, 3*nn**3,1)

                EN_err = EN_err + dot_product(jump_all,u_p)
                !write(*,*) dot_product(jump_all,u_p)
                !read(*,*)
                !H1_err = H1_err + dot_product(jump_all,u_p)

                 deallocate(u_p, jump_all)

      enddo

     Linf_err_tot = dabs(Linf_err)
     !Linf_err_vel_tot = dabs(Linf_err_vel)

     L2_err_tot = dabs(L2_err)
     L2_err_vel_tot = dabs(L2_err_vel)
     H1_err_tot = dabs(H1_err)
     EN_err_tot = dabs(En_err)

     if (EN_err_tot .ge. 100.d0) call EXIT(EXIT_ENERGY_ERROR)

     if(mpi_id .eq. 0) then
        file_mpi = 'EN.ERR'; unit_mpi = 40;
        open(unit_mpi,file=file_mpi, position='APPEND')
        write(unit_mpi,*)  sqrt(H1_err_tot), sqrt(EN_err_tot)!, sqrt(L2_err_tot),sqrt(L2_err_vel_tot)

        close(unit_mpi)
     endif

     end subroutine COMPUTE_ENERGY_ERROR


text:
      double complex function Faxsl(st,j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      integer j1,j2,j3,j4,j5,j6
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'masses.f'
      character*9 st
      double complex L1
      double precision t,mtsq
      mtsq=mt**2
      if(st.eq.'q+qb-g+g-') then
      Faxsl=
     .(2d0*(1/(24d0*mtsq)-L1(-t(j1,j2,j4),-s(j5,j6))/(2d0*s(j5,j6)))*zb(
     .j1,j3)*
     .(za(j2,j5)*zb(j1,j2)+za(j4,j5)*zb(j1,j4))*zb(j3,j6))/
     .(s(j5,j6)*zb(j1,j4)*zb(j2,j4))+
     .(2d0*(1/(24d0*mtsq)-L1(-t(j1,j2,j3),-s(j5,j6))/(2d0*s(j5,j6)))*za(
     .j2,j4)*za(j4,j5)*
     .(-(za(j1,j2)*zb(j1,j6))+za(j2,j3)*zb(j3,j6)))/
     .(s(j5,j6)*za(j1,j3)*za(j2,j3))
      elseif(st.eq.'q+qb-g+g+') then
      Faxsl=
     .(2d0*(1/(24d0*mtsq)-L1(-t(j1,j2,j4),-s(j5,j6))/(2d0*s(j5,j6)))*za(
     .j2,j5)*
     .(-(za(j1,j2)*zb(j1,j3))-za(j2,j4)*zb(j3,j4))*zb(j3,j6))/
     .(s(j5,j6)*za(j1,j4)*za(j2,j4))+
     .(2d0*(1/(24d0*mtsq)-L1(-t(j1,j2,j3),-s(j5,j6))/(2d0*s(j5,j6)))*za(
     .j2,j5)*
     .(-(za(j1,j2)*zb(j1,j4))+za(j2,j3)*zb(j3,j4))*zb(j4,j6))/
     .(s(j5,j6)*za(j1,j3)*za(j2,j3))
      endif
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen