text:
program test_ieee_arithmetic
use, intrinsic :: ieee_arithmetic
implicit none
character (len=20) :: fmt = "(a16,*(l9))"
integer, parameter :: wp = kind(1.0) ! same output for wp = kind(1.0d0)
real(kind=wp)      :: z,vec(6),NaN,t
z = 0.0_wp
t = tiny(z)
vec = [z/z,z/1.0_wp,-z/1.0_wp,1.0_wp/z,-1.0_wp/z,t/10]
print*,ieee_value(0.0,ieee_positive_inf),1.0/z ! 2 ways of geting +Inf
NaN = ieee_value(0.0,ieee_quiet_nan) ! get NaN
print*,"NaN == NaN?",NaN == NaN ! demonstrate that NaN /= NaN
print*
print "(19x,*(a9))","0.0/0.0","0.0/1.0","-0.0/1.0","1.0/0.0","-1.0/0.0","tiny/10"
print fmt,"ieee_is_nan"     ,ieee_is_nan(vec)
print fmt,"ieee_is_negative",ieee_is_negative(vec)
print fmt,"ieee_is_finite"  ,ieee_is_finite(vec)
print fmt,"ieee_is_normal"  ,ieee_is_normal(vec)
end program test_ieee_arithmetic
! gfortran output:
!         Infinity         Infinity
! NaN == NaN? F
!
!                      0.0/0.0  0.0/1.0 -0.0/1.0  1.0/0.0 -1.0/0.0  tiny/10
!      ieee_is_nan        T        F        F        F        F        F
! ieee_is_negative        F        F        T        F        T        F
!   ieee_is_finite        F        T        T        F        F        T
!   ieee_is_normal        F        T        T        F        F        F

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
         IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      parameter(nn=340,nch=1005,ntp=200,pi=3.141592654,
     $xlambd=0.0005543,ro= 1.987d-03,x1=1.2d0,x2=3.d0,acc=1.d-5)
	dimension fe(10,ntp),dzxe(10,ntp),tinv(ntp)
	dimension c(10),af(ntp),zita(ntp)
	dimension t(ntp),f(0:ntp),zx(ntp),dtim(ntp),e(10),d0(10)
	character*9 tab1
	double precision an, bessj0, bessj1, rtsafe
      external an, funcd, bessj0, bessj1, rtsafe
	common /integ/ imod

	tab1=char(09)
	open(unit=7,file='temstep.in',status='old')
	open(unit=15,file='arr-me.in',status='old')
	open(unit=16,file='arr-me.out',status='unknown')
	open(unit=8,file='each-me.gr',status='unknown')
	open(unit=10,file='arr-me.dat',status='unknown')
	open(unit=1,file='each-me.dat',status='unknown')
	open(unit=11,file='logr-me.dat',status='unknown')
c	open(unit=20,file='test.dat',status='unknown')
	imp = 1
	print *, 'slabs= 1, spheres= 2 , cilinders= 3'
	read *, imod
	print *, 'model =', imod
      b = 10.d0 - 2.d0 * imod
      if(imod.eq.1)then
         imp = 2
         acut = 0.5
      endif
	if(imod.eq.2)then
	   acut = 0.85
	elseif(imod.eq.3)then
	   acut = 0.6
	endif

	read(15,120)nsamp
	do 11 j=1, nsamp
	  zita(j)=0.
   	  read(15,140)e(j),d0(j),c(j)
	  tc = tc + c(j)
	  d0(j) = 10.**d0(j)
	  if (imod.eq.1) d0(j)=d0(j)/4.
11    continue
	if(abs(tc-1.).gt.1e-8) then
	print *,'Error  -  your concentration c(j) sum ',tc
	stop
	endif
	read(15,*)slop,ord
  120	format(I20)
	read(7,120)nstep
	do 10 nt=1,nstep
	f(nt)=0.
	read(7,*)t(nt),dtim(nt)
	dtim(nt)=60.*dtim(nt)
	do 20 j=1,nsamp
	fe(j,nt)=1.
	zita(j)=d0(j)*dtim(nt)*dexp(-e(j)/ro/t(nt))+zita(j)
	do 30 m=1,40000
	xu=an(m)*zita(j)
	if (xu.gt.99.)goto 20
	fe(j,nt)=fe(j,nt)-b*dexp(-xu)/an(m)
   30	f(nt)=f(nt)-b*c(j)*dexp(-xu)/an(m)
   20	continue
	nnt=nt
	f(nt)=1.+ f(nt)
	if (f(nt).gt..9999)then
	nnt=nt-1
	goto 40
	endif
   10	continue
   40	continue

c	INVERSION OF THE DATA

 	do 50 k=1,nnt
        if (f(k).le.acut) then
	    if(imod.eq.1)zx(k+1)=pi*(f(k)/4.)**2
	    if(imod.eq.2)zx(k+1)=(2.-pi/3.*f(k)-2.*
     $    dsqrt(1.-pi/3.*f(k)))/pi
          if(imod.eq.3)then
            sb = b**2/pi
            zx(k+1)=((sb-2.d0*f(k))-dsqrt(sb*(sb-4.d0*f(k))))/2.d0
          endif
        else
          zx(k+1)=-dlog(an(1)/b*(1.-f(k)))/an(1)
        endif
          if(imod.eq.3)then
            sb = b**2/pi
            z1=((sb-2.d0*f(k))-dsqrt(sb*(sb-4.d0*f(k))))/2.d0
            z2=-dlog(an(1)/b*(1.-f(k)))/an(1)
c           write(20,*)f(k),z1,z2
          endif

   50   continue
        zx(1)=0.
        do 80 k=1,nnt
        af(k)=(f(k)+f(k-1))/2.*100.
        dzx = dlog10((zx(k+1)-zx(k))/dtim(k)*imp**2)
        tinv(k)=1./t(k)*10000.
        xlogr=(ord-slop*tinv(k)-dzx)/2.
	write(10,110)tinv(k),tab1,dzx
	write(11,110)100.*f(k-1),tab1,xlogr
	write(11,110)100.*(f(k)-0.000001),tab1,xlogr
  80   write(16,110)tinv(k),tab1,dzx,tab1,f(k)*100.,tab1,xlogr,tab1
     $  ,af(k)

	do 200 j=1,nsamp
	do 150 k=1,nnt
	  if (fe(j,k).le.acut) then
	    if(imod.eq.1)zx(k+1)=pi*(fe(j,k)/4.)**2
	    if(imod.eq.2)zx(k+1)=(2.-pi/3.*fe(j,k)-2.*
     $    dsqrt(1.-pi/3.*fe(j,k)))/pi
          if(imod.eq.3)then
            sb = b**2/pi
            zx(k+1)=((sb-2.d0*fe(j,k))-dsqrt(sb*(sb-4.d0*fe(j,k))))/2.d0
          endif
        else
          if (fe(j,k).lt.1.)then
             zx(k+1)=-dlog(an(1)/b*(1.-fe(j,k)))/an(1)
	    else
	      zx(k+1) = 0.
	    endif
        endif

150   continue
      zx(1)=0.
      do 180 k=1,nnt
	if (zx(k+1).gt.0.) then
        dzxe(j,k) = dlog10((zx(k+1)-zx(k))*imp**2/dtim(k))
	endif
  180	continue
  200	continue
	do 195 j=1,nsamp
	do 190 k=1,nnt
	write(8,110)tinv(k),tab1,fe(j,k)
	if (dzxe(j,k).eq.0.)goto 190
        write(1,110)tinv(k),tab1,dzxe(j,k)
  190	continue
	write(8,130)
	write(1,130)
  195	continue

  110	format(1x,5(f12.8,a1))
  130   format(1x,'$')
  140	format(g20.8)
 	end

       double precision function an(m)
       implicit double precision (a-h,o-z)
       parameter(pi=3.14159265,x1=1.2d0,x2=3.d0,acc=1.d-5)
       external funcd
       common /integ/ imod
       if(imod.eq.1)an = ((2*m-1)*pi)**2
       if(imod.eq.2)an = (m*pi)**2
       if(imod.eq.3.and.m.gt.1)then
         den=pi*(4*m-1)
         zero= den/4.+1./(2.*den)-31./(6*den**3)+3779./(15*den**5)
         an = zero**2
       elseif(imod.eq.3)then
         an=rtsafe(funcd,x1,x2,acc)**2
       endif
       return
       end

       subroutine funcd(x1,f,df)
	implicit double precision (a-h,o-z)
           f=bessj0(x1)
           df= - bessj1(x1)
       return
       end

      double precision function bessj0(x)
       implicit double precision (a-h,o-z)
c      real*8 y,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
C     *    s1,s2,s3,s4,s5,s6
      data p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     *    -.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-
     *1,
     *    .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      data r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,651619640.7d
     *0,
     *    -11214424.18d0,77392.33017d0,-184.9052456d0/,
     *    s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,
     *    9494680.718d0,59272.64853d0,267.8532712d0,1.d0/
      if(dabs(x).lt.8.)then
        y=x**2
        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
     *      /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
      else
        ax=dabs(x)
        z=8./ax
        y=z**2
        xx=ax-.785398164
        bessj0=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+y*(p4+y
     *      *p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      return
      end

      double precision function bessj1(x)
	implicit double precision (a-h,o-z)
C      real*8 y,p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
c     *    s1,s2,s3,s4,s5,s6
      data r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,242396853.1d0
     *,
     *    -2972611.439d0,15704.48260d0,-30.16036606d0/,
     *    s1,s2,s3,s4,s5,s6/144725228442.d0,2300535178.d0,
     *    18583304.74d0,99447.43394d0,376.9991397d0,1.d0/
      data p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,.2457520174d-5
     *,
     *    -.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,-.2002690873d-3
     *,
     *    .8449199096d-5,-.88228987d-6,.105787412d-6/
      if(dabs(x).lt.8.)then
        y=x**2
        bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
     *      /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
      else
        ax=dabs(x)
        z=8./ax
        y=z**2
        xx=ax-2.356194491
        bessj1=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+y*(p4+y
     *      *p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
     *      *dsign(1.d0,x)
      endif
      return
      end

      double precision function rtsafe(funcd,x1,x2,xacc)
        implicit double precision (a-h,o-z)
      parameter (maxit=100)
      call funcd(x1,fl,df)
      call funcd(x2,fh,df)
      if(fl*fh.ge.0.) pause 'root must be bracketed'
      if(fl.lt.0.)then
        xl=x1
        xh=x2
      else
        xh=x1
        xl=x2
        swap=fl
        fl=fh
        fh=swap
      endif
      rtsafe=.5*(x1+x2)
      dxold=dabs(x2-x1)
      dx=dxold
      call funcd(rtsafe,f,df)
      do 11 j=1,maxit
        if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0.
     *      .or. dabs(2.*f).gt.dabs(dxold*df) ) then
          dxold=dx
          dx=0.5*(xh-xl)
          rtsafe=xl+dx
          if(xl.eq.rtsafe)return
        else
          dxold=dx
          dx=f/df
          temp=rtsafe
          rtsafe=rtsafe-dx
          if(temp.eq.rtsafe)return
        endif
        if(dabs(dx).lt.xacc) return
        call funcd(rtsafe,f,df)
        if(f.lt.0.) then
          xl=rtsafe
          fl=f
        else
          xh=rtsafe
          fh=f
        endif
11    continue
      pause 'rtsafe exceeding maximum iterations'
      return
      end

ext: f

max_stars_repo_name: ASUPychron/pychron

max_issues_repo_name: ASUPychron/pychron

max_forks_repo_name: UIllinoisHALPychron/pychron


text:
!+++++++++++++++++++++++++++++++++++++!
!                                     !
!     Global variable definitions     !
!       for the mesh generator        !
!                                     !
!+++++++++++++++++++++++++++++++++++++!
!..RCS/CVS ident
! $Id: gen_mod.h90,v 1.8 2000/03/22 21:10:47 bojan Exp $
! $Source: /home/muhamed/.CVSROOT/T-Rex/Library/gen_mod.h90,v $
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>!
MODULE gen_mod

  USE allp_mod

  IMPLICIT NONE

  !INTEGER,PARAMETER :: MAXP   = 200

  REAL,ALLOCATABLE,PUBLIC    :: x(:),  y(:),  z(:)   ! node coordinates
  REAL,ALLOCATABLE,PUBLIC    :: walln(:)             ! node distance from the wall
  INTEGER,ALLOCATABLE,PUBLIC :: SideN(:,:)           ! numb, n1, n2, n3, n4
  INTEGER,ALLOCATABLE,PUBLIC :: SideCc(:,:)

  INTEGER,ALLOCATABLE,PUBLIC :: CellC(:,:)           ! cell's neighbours
  INTEGER,ALLOCATABLE,PUBLIC :: CellN(:,:)           ! cell nodes

  INTEGER,ALLOCATABLE,PUBLIC :: TwinN(:,:)

  INTEGER,ALLOCATABLE,PUBLIC :: NewN(:)    ! new number for the nodes and cells
  INTEGER,ALLOCATABLE,PUBLIC :: NewC(:)    ! new number for cells
  INTEGER,ALLOCATABLE,PUBLIC :: NewS(:)    ! new number for sides
  INTEGER,ALLOCATABLE,PUBLIC :: CelMar(:)  ! cell marker

  INTEGER,ALLOCATABLE,PUBLIC :: NodeN2(:,:)
  INTEGER,ALLOCATABLE,PUBLIC :: NodeN4(:,:)
  INTEGER,ALLOCATABLE,PUBLIC :: NodeN8(:,:)

  INTEGER,ALLOCATABLE,PUBLIC :: level(:)   ! refinement level

  INTEGER,PUBLIC :: MAXN, MAXB, MAXS

  INTEGER,PUBLIC :: NR(MAXP)               ! refin. levels, refin. regions

  REAL,PUBLIC    :: xp(MAXP), yp(MAXP), zp(MAXP)       ! point coordinates
  REAL,PUBLIC    :: xl(MAXP,MAXL),yl(MAXP,MAXL),zl(MAXP,MAXL),LinWgt(MAXP)
  REAL,PUBLIC    :: BlkWgt(MAXL,3), BlFaWt(MAXL,3)     ! leave this
  REAL,PUBLIC    :: FRegio(MAXP,MAXP,0:6)              ! levels, regions

  REAL,PUBLIC    :: SRegio(MAXP,0:6), Srelax(MAXP)  ! levels, regions
  LOGICAL,PUBLIC :: SdirX(MAXP), SdirY(MAXP), SdirZ(MAXP)
  INTEGER,PUBLIC :: Siter(MAXP)

  INTEGER,PUBLIC :: BlkPnt(MAXP,0:8),  & ! 0 for orientation
	     BlkRes(MAXP,6),    & ! NI,NJ,NK,NI*NJ*NK,NNo,NVo
	     BlkFac(MAXP,6,4),  &
	     BlFaLa(MAXP),      &
	     Bound(MAXP,8),     &
	     Period(MAXP,8),    &
	     Copy(MAXP,0:8)

  INTEGER,PUBLIC :: LinPnt(MAXL,2), LinRes(MAXL)
  INTEGER,PUBLIC :: Nbloc, NP, Nline, Nsurf, Nboun, Nperi
  INTEGER,PUBLIC :: NN, NN2, NN4, NN8
  INTEGER,PUBLIC:: NSR                   ! smoothing regions
  INTEGER,PUBLIC :: NSsh                  ! number of shadow faces

  INTEGER,PUBLIC :: WallFacFst, WallFacLst

  INTEGER,PUBLIC :: ELIPSO, RECTAN, PLANE,YES,NO

  CHARACTER*4,PUBLIC :: BndFac(MAXP)

END MODULE

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      LOGICAL FUNCTION PLTGTG(INDX,BUFF)
      CHARACTER*16 IERROR
      REAL DEVCAP(23)
      REAL DEFOUT(7)
      COMMON /STATUS/DEVCAP,DEFOUT
      REAL DEVP(5)
      COMMON /DEVICE/DEVP
      REAL COLP(3)
      REAL PALETT(3,16)
      COMMON /COLOR/COLP,PALETT
      REAL TEXTP(40)
      COMMON /TEXT/TEXTP
      REAL VECTP(5)
      REAL XCUR
      REAL YCUR
      COMMON /VECTRC/VECTP,XCUR,YCUR
      INTEGER IDEX(200,2)
      INTEGER NVECT(200,2)
      REAL XSIZE(200,2)
      REAL YSIZE(200,2)
      REAL X0(2300,2)
      REAL Y0(2300,2)
      REAL X1(2300,2)
      REAL Y1(2300,2)
      COMMON /FONT/IDEX,NVECT,XSIZE,YSIZE,X0,Y0,X1,Y1
      REAL GRAPHP(100)
      COMMON /GRAPH/GRAPHP
      COMMON /MAPPAR/MAPP(11)
      REAL MAPP
      COMMON /STORAG/MEMORY(1000)
      REAL BUFF(*)

      PLTGTG = .TRUE.
      IF (INDX.EQ.1) THEN
         BUFF(1) = GRAPHP(1)

      ELSE IF (INDX.EQ.2) THEN
         BUFF(1) = GRAPHP(2)

      ELSE IF (INDX.EQ.3) THEN
         BUFF(1) = GRAPHP(3)

      ELSE IF (INDX.EQ.4) THEN
         BUFF(1) = GRAPHP(4)

      ELSE IF (INDX.EQ.5) THEN
         BUFF(1) = GRAPHP(5)

      ELSE IF (INDX.EQ.6) THEN
         BUFF(1) = GRAPHP(38)

      ELSE IF (INDX.EQ.7) THEN
         IF (GRAPHP(6).EQ.1.) THEN
            BUFF(1) = GRAPHP(47) - 4.

         ELSE
            BUFF(1) = 0.
         END IF

      ELSE IF (INDX.EQ.8) THEN
         BUFF(1) = GRAPHP(23)

      ELSE IF (INDX.EQ.9) THEN
         BUFF(1) = GRAPHP(21)

      ELSE IF (INDX.EQ.10) THEN
         BUFF(1) = GRAPHP(37)

      ELSE IF (INDX.EQ.11) THEN
         IF (BUFF(1).EQ.3.) THEN
            BUFF(2) = GRAPHP(24)
            BUFF(3) = GRAPHP(25)
            BUFF(4) = GRAPHP(26)
            BUFF(5) = GRAPHP(27)
            BUFF(6) = GRAPHP(28)
            BUFF(7) = GRAPHP(29)
            BUFF(8) = GRAPHP(30)
            BUFF(9) = GRAPHP(31)
         END IF

         IF (BUFF(1).EQ.4.) THEN
            BUFF(2) = GRAPHP(78)
            BUFF(3) = GRAPHP(79)
            BUFF(4) = GRAPHP(80)
            BUFF(5) = GRAPHP(81)
            BUFF(6) = GRAPHP(82)
            BUFF(7) = GRAPHP(83)
            BUFF(8) = GRAPHP(84)
            BUFF(9) = GRAPHP(85)
            BUFF(10) = GRAPHP(86)
            BUFF(11) = GRAPHP(87)
         END IF

         BUFF(1) = GRAPHP(22)

      ELSE IF (INDX.EQ.12) THEN
         BUFF(1) = GRAPHP(32)

      ELSE IF (INDX.EQ.13) THEN
         BUFF(1) = GRAPHP(91)

      ELSE IF (INDX.EQ.14) THEN
         BUFF(1) = GRAPHP(90)

      ELSE IF (INDX.EQ.15) THEN
         BUFF(1) = GRAPHP(35)

      ELSE IF (INDX.EQ.16) THEN
         BUFF(1) = GRAPHP(36)

      ELSE IF (INDX.EQ.17) THEN
         BUFF(1) = GRAPHP(39)

      ELSE IF (INDX.EQ.18) THEN
         BUFF(1) = GRAPHP(40)

      ELSE IF (INDX.EQ.19) THEN
         BUFF(1) = GRAPHP(41)

      ELSE IF (INDX.EQ.20) THEN
         BUFF(1) = GRAPHP(42)

      ELSE IF (INDX.EQ.21) THEN
         BUFF(1) = GRAPHP(92)

      ELSE IF (INDX.EQ.22) THEN
         BUFF(1) = GRAPHP(44)

      ELSE IF (INDX.EQ.23) THEN
         BUFF(1) = GRAPHP(45)

      ELSE IF (INDX.EQ.47) THEN
         BUFF(1) = GRAPHP(88)

      ELSE IF (INDX.EQ.48) THEN
         BUFF(1) = GRAPHP(89)

      ELSE IF (INDX.EQ.24) THEN
         BUFF(1) = GRAPHP(46)

      ELSE IF (INDX.EQ.27) THEN
         DO 2280 I = 0,13
            BUFF(I+1) = GRAPHP(7+I)
 2280    CONTINUE

      ELSE IF (INDX.EQ.28) THEN
         BUFF(1) = GRAPHP(62)

      ELSE IF (INDX.EQ.29) THEN
         BUFF(1) = GRAPHP(63)

      ELSE IF (INDX.EQ.30) THEN
         BUFF(1) = GRAPHP(64)

      ELSE IF (INDX.EQ.31) THEN
         BUFF(1) = GRAPHP(65)

      ELSE IF (INDX.EQ.32) THEN
         BUFF(1) = GRAPHP(66)

      ELSE IF (INDX.EQ.33) THEN
         BUFF(1) = GRAPHP(67)

      ELSE IF (INDX.EQ.34) THEN
         BUFF(1) = GRAPHP(68)

      ELSE IF (INDX.EQ.35) THEN
         BUFF(1) = GRAPHP(69)

      ELSE IF (INDX.EQ.36) THEN
         BUFF(1) = GRAPHP(70)

      ELSE IF (INDX.EQ.37) THEN
         BUFF(1) = GRAPHP(71)

      ELSE IF (INDX.EQ.38) THEN
         BUFF(1) = GRAPHP(72)

      ELSE IF (INDX.EQ.39) THEN
         BUFF(1) = GRAPHP(73)

      ELSE IF (INDX.EQ.43) THEN
         BUFF(1) = GRAPHP(74)

      ELSE IF (INDX.EQ.44) THEN
         BUFF(1) = GRAPHP(75)

      ELSE IF (INDX.EQ.45) THEN
         BUFF(1) = GRAPHP(76)

      ELSE IF (INDX.EQ.46) THEN
         BUFF(1) = GRAPHP(77)

      ELSE
         CALL CHRIC(INDX,IERROR,L)
         CALL PLTFLU
         CALL SIORPT('PLTGTG','Illegal index '//IERROR(1:L)//'.',2)
         PLTGTG = .FALSE.
         RETURN

      END IF

      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!##############################################################################
!# ****************************************************************************
!# <name> AllenCahn_partridiscr </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module contains routines to read the parametrisation, create the
!# triangulation and set up the discretisation for the heat conduction problem.
!# The following routines can be found here:
!#
!# 0.) AC5_initSolution
!# 1.) AC5_initParamTriang
!#     -> Read .PRM/.TRI files. Generate meshes on all levels.
!#
!# 2.) AC5_doneParamTriang
!#     Clean up parametrisation/triangulation, release memory.
!#
!# 3.) AC5_initDiscretisation
!#     -> Initialise the spatial discretisation.
!#
!# 4.) AC5_doneDiscretisation
!#     -> Clean up the discretisation, release memory.
!# </purpose>
!##############################################################################

MODULE AllenCahn_partridiscr

  USE fsystem
  USE storage
  USE linearsolver
  USE boundary
  USE bilinearformevaluation
  USE linearformevaluation
  USE cubature
  USE matrixfilters
  USE vectorfilters
  USE bcassembly
  USE triangulation
  USE spatialdiscretisation
  USE sortstrategy
  USE coarsegridcorrection
  USE ucd
  USE timestepping
  USE genoutput

  use mprimitives

  USE collection
  USE paramlist
  use convection
  use vectorio

  USE AllenCahn_callback

  USE AllenCahn_basic

  IMPLICIT NONE

CONTAINS

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

!<subroutine>
  subroutine AC5_initSolution(rproblem, rvector)
!<description>
  ! Initialises the initial solution vector into rvector. Depending on the settings
  ! in the DAT file this is either zero or read from a file.
  !
  ! The routine assumes that basic mass matrices have already been constructed.
!</description>

!<input>
  ! A problem structure saving problem-dependent information.
  type(t_ACproblem), intent(INOUT) :: rproblem
!</input>

!<inputoutput>
  ! The solution vector to be initialised. Must be set up according to the
  ! maximum level NLMAX in rproblem!
  type(t_vectorBlock), intent(INOUT) :: rvector
!</inputoutput>

!</subroutine>

    ! local variables
    integer(I32) :: istart,ctypeInitialSolution
    type(t_vectorBlock) :: rvector1,rvector2
    type(t_vectorScalar) :: rvectorTemp
    character(LEN=SYS_STRLEN) :: sarray,sfile,sfileString
    integer :: ilev,ierror
    integer(PREC_VECIDX) :: NEQ
    type(t_interlevelProjectionBlock) :: rprojection
    type(t_linearForm) :: rlinform
    type(t_blockDiscretisation), pointer :: p_rdiscretisation
    type(t_linsolNode), pointer :: p_rsolverNode,p_rpreconditioner
    type(t_matrixBlock), dimension(1) :: Rmatrices
!    type(t_matrixBlock) :: Rmatrices

    type(t_vectorBlock) :: rsingleRHS,rsingleSol
    type(t_vectorscalar) :: rtemp

!     logical :: mcai
    ! Get the parameter what to do with rvector

    ctypeInitialSolution=3

	!istart
    !sfileString

    ! What is the type of the initial solution?
    select case (ctypeInitialSolution)
    case (0)
      ! Init with zero
      call lsysbl_clearVector (rvector)

    case (1,2)
      ! We have to read a file -- formatted or unformatted.
      !
      ! Create a temp vector at level NLMAX-istart+1.
!      if (istart .gt. 0) then
!        ilev = istart
!      else
!        ilev = rproblem%NLMAX-abs(istart)+1
!      end if

!      if (ilev .lt. rproblem%NLMIN) then
!        call output_line (&
!            'Level of start vector is < NLMIN! Initialising with zero!', &
!            OU_CLASS_WARNING,OU_MODE_STD,'cc_initInitialSolution')
!        istart = 0
!      end if

!      if (ilev .gt. rproblem%NLMAX) then
!        call output_line (&
!            'Level of start vector is > NLMAX! Initialising with zero!', &
!            OU_CLASS_WARNING,OU_MODE_STD,'cc_initInitialSolution')
!        istart = 0
!      end if

      ! Remove possible ''-characters
!      read(sfileString,*) sfile

!      call lsysbl_createVectorBlock (&
!          rproblem%RlevelInfo(ilev)%rdiscretisation,rvector1,.false.)

      ! Read in the vector
!      call vecio_readBlockVectorHR (&
!          rvector1, sarray, .true., 0, sfile, istart .gt. 0)

      ! If the vector is on level < NLMAX, we have to bring it to level NLMAX
!      do while (ilev .lt. rproblem%NLMAX)

        ! Initialise a vector for the higher level and a prolongation structure.
!        call lsysbl_createVectorBlock (&
!            rproblem%RlevelInfo(ilev+1)%rdiscretisation,rvector2,.false.)

!        call mlprj_initProjectionVec (rprojection,rvector2)

        ! Prolongate to the next higher level.

!        NEQ = mlprj_getTempMemoryVec (rprojection,rvector1,rvector2)
!        if (NEQ .ne. 0) call lsyssc_createVector (rvectorTemp,NEQ,.false.)
!        call mlprj_performProlongation (rprojection,rvector1, &
!                                        rvector2,rvectorTemp)
!        if (NEQ .ne. 0) call lsyssc_releaseVector (rvectorTemp)

        ! Swap rvector1 and rvector2. Release the coarse grid vector.
!        call lsysbl_swapVectors (rvector1,rvector2)
!        call lsysbl_releaseVector (rvector2)

!        call mlprj_doneProjection (rprojection)

        ! rvector1 is now on level ilev+1
!        ilev = ilev+1

!      end do

      ! Copy the resulting vector rvector1 to the output.
!      call lsysbl_copyVector (rvector1,rvector)

      ! Release the temp vector
!      call lsysbl_releaseVector (rvector1)

    case (3)

      ! We have to create the solution by analytical callback functions.
      ! To do this for an arbitrary finite element, we have to do an
      ! L2 projection for X,Y and P. That means we have to solve:
      !     <u,phi> = <u_analytic,phi>
      ! which means in the FE context:
      !     Mu = f
      ! with f a RHS created by analytical point values.

!      if (rproblem%MSHOW_Initialisation .ge. 1) &
      call output_line('AC problem, Preparing L2 projection of analytical initial solution...')

      ! Get a pointer to the RHS on the finest level as well as to the
      ! block discretisation structure:
      p_rdiscretisation => rvector%p_rblockDiscr

      ! The vector structure is already prepared, but the entries are missing.
      !
      ! At first set up the corresponding linear form (u_analytical,Phi_j):
      rlinform%itermCount = 1
      rlinform%Idescriptors(1) = DER_FUNC

      ! Create a temp vector as RHS.

	  !call lsysbl_createVectorBlock (&
      !    rproblem%RlevelInfo(rproblem%NLMAX)%rdiscretisation,rvector1,.false.)
	  ! if rdiscretisation is not a field of t_ACproblem_lvl, we use
      call lsysbl_createVecBlockIndirect (rvector,rvector1,.TRUE., ST_double)

      ! Assemble the RHS.
      !
      ! Initialise the collection for the assembly process with callback routines.
      ! Basically, this stores the simulation time in the collection if the
      ! simulation is nonstationary.
      call AC5_initCollectForAssembly (rproblem,rproblem%rcollection)

      ! Discretise the X-velocity part:
      call linf_buildVectorScalar (&
                p_rdiscretisation%RspatialDiscr(1),rlinform,.true.,&
                rvector1%RvectorBlock(1),coeff_AnalyticSolution_AC,&
                rproblem%rcollection)

      ! Clean up the collection (as we are done with the assembly.
      call AC5_doneCollectForAssembly (rproblem,rproblem%rcollection)

      ! Now prepare a linear solver for solving with the mass matrix.
      ! This is just a simple 1-level Jacobi solver as the mass matrix
      ! is usually well conditioned.
      call linsol_initJacobi (p_rpreconditioner, 0.8_DP)
      call linsol_initDefCorr (p_rsolverNode, p_rpreconditioner)

      p_rsolverNode%depsRel = SYS_EPSREAL_DP * 100.0_DP
      p_rsolverNode%nmaxiterations = 1000
      p_rsolverNode%ioutputLevel = 1

      call output_line('AC problem, Solving L2 projection for \phi...')

      ! We solve separately for the three components.
      ! Prepare the solver for the velocity.
      ! Prepare the solver for the X-velocity.

! We should make sure that rmatrixMass has been assembled.

     call lsysbl_createMatFromScalar(&
        rproblem%RlevelInfo(rproblem%NLMAX)%rmatrixMass%RmatrixBlock(1,1),Rmatrices(1))

      call lsyssc_copyvector(rvector%RvectorBlock(1), rtemp)
      call lsyssc_scalevector(rtemp, 0.0_DP)
!      call lsyssc_synchroniseSortMatVec (Rmatrices(1)%RmatrixBlock(1,1), &
!	                                     rvector%RvectorBlock(1), rtemp)

      call linsol_setMatrices (p_rsolverNode,Rmatrices)
      call linsol_initStructure (p_rsolverNode,ierror)
      call linsol_initData (p_rsolverNode,ierror)

      ! -----
      ! Solve for the X-velocity.
      call lsysbl_createVecFromScalar(rvector1%RvectorBlock(1),rsingleRHS)
      call lsysbl_createVecFromScalar(rvector%RvectorBlock(1),rsingleSol)
      call lsysbl_duplicateVector (rsingleSol,rvector2,&
          LSYSSC_DUP_COPY,LSYSSC_DUP_EMPTY)

      call linsol_solveAdaptively (p_rsolverNode,rsingleSol,rsingleRHS,rvector2)
!~~~~~We should copy back rsingleSol to rvecto%RvectorBlock(1)

      if (p_rsolverNode%iresult .ne. 0) then
        ! Cancel, there is something wrong.
        call output_line('Cannot compute L2 projection, solver broke down. Using zero!')
        call lsysbl_clearVector (rvector)

!      else
!        call lsysbl_copyVector(rsingleSol, rvector)
      end if
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      ! -----
      ! That's it, cleanup.
      call lsysbl_releaseVector (rsingleRHS)
      call lsysbl_releaseVector (rsingleSol)
      call lsysbl_releaseVector (rvector2)
      call lsysbl_releaseMatrix (Rmatrices(1))

	  call lsyssc_releaseVector (rtemp)

      call linsol_doneData (p_rsolverNode,ierror)
      call linsol_doneStructure (p_rsolverNode,ierror)
      call linsol_releaseSolver (p_rsolverNode)

    end select

  end subroutine

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

!<subroutine>

  subroutine AC5_initCollectForAssembly (rproblem,rcollection)

!<description>
  ! This subroutine is an auxiliary subroutine called by the CC2D framework
  ! and has usually not to be changed by the user.
  !
  ! The subroutine prepares the collection rcollection to be passed to callback
  ! routines for assembling boundary conditions or RHS vectors. It is
  ! called directly prior to the assembly to store problem-specific information
  ! in the quick-access arrays of the collection.
  ! Basically speaking, the routine stores information about whether thw problem
  ! is stationary, nonstationary or about the current simulation time.
!</description>

!<input>
  ! Problem structure with all problem relevant data.
  type(t_ACproblem), intent(IN) :: rproblem
!</input>

!<inputoutput>
  ! Collection structure to be initialised.
  type(t_collection), intent(INOUT) :: rcollection
!</inputoutput>

!</subroutine>

    ! In a nonstationary simulation, save the simulation time as well as the
    ! minimum and maximum time to the quick-access array of the collection,
    ! so it can be accessed in the callback routines!
!    rcollection%Iquickaccess(1) = rproblem%itimedependence
    rcollection%Dquickaccess(1) = rproblem%rtimedependence%dtime

  end subroutine

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

!<subroutine>

  subroutine AC5_doneCollectForAssembly (rproblem,rcollection)

!<description>
  ! This subroutine is an auxiliary subroutine called by the CC2D framework
  ! and has usually not to be changed by the user.
  !
  ! After the assembly process, this subroutine is called to release temporary
  ! information from the collection which was stored there by
  ! cc_initCollectForAssembly.
!</description>

!<input>
  ! Problem structure with all problem relevant data.
  type(t_ACproblem), intent(IN) :: rproblem
!</input>

!<inputoutput>
  ! Collection structure to be cleaned up.
  type(t_collection), intent(INOUT) :: rcollection
!</inputoutput>

!</subroutine>

    ! Currently, this subroutine is empty as all information stored in
    ! the collection in cc_initCollectForAssembly is put to the quick-access
    ! arrays -- which do not have to be cleaned up.
    ! This might change in future...

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

!<subroutine>

  SUBROUTINE AC5_initParamTriang (NLMIN,NLMAX,rACproblem)

!<description>
  ! This routine initialises the parametrisation and triangulation of the
  ! domain. The corresponding .prm/.tri files are read from disc and
  ! the triangulation is refined as described by the parameter ilv.
!</description>

!<input>
  ! Minimum refinement level of the mesh; = coarse grid = level 1
  INTEGER, INTENT(IN) :: NLMIN

  ! Maximum refinement level
  INTEGER, INTENT(IN) :: NLMAX
!</input>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  TYPE(t_ACproblem), INTENT(INOUT) :: rACproblem
!</inputoutput>

!</subroutine>

  ! local variables
  INTEGER :: i

    ! Initialise the level in the problem structure

    rACproblem%NLMIN = NLMIN
    rACproblem%NLMAX = NLMAX

    ! At first, read in the parametrisation of the boundary and save
    ! it to rboundary.
    CALL boundary_read_prm(rACproblem%rboundary, './pre/QUAD.prm')

    ! Now read in the basic triangulation.
    CALL tria_readTriFile2D (rACproblem%RlevelInfo(rACproblem%NLMIN)%rtriangulation, &
        './pre/QUAD.tri', rACproblem%rboundary)

    ! Refine the mesh up to the minimum level
    CALL tria_quickRefine2LevelOrdering(rACproblem%NLMIN-1,&
        rACproblem%RlevelInfo(rACproblem%NLMIN)%rtriangulation,rACproblem%rboundary)

    ! Create information about adjacencies and everything one needs from
    ! a triangulation. Afterwards, we have the coarse mesh.
    CALL tria_initStandardMeshFromRaw (&
        rACproblem%RlevelInfo(rACproblem%NLMIN)%rtriangulation,rACproblem%rboundary)

    ! Now, refine to level up to nlmax.
    DO i=rACproblem%NLMIN+1,rACproblem%NLMAX
      CALL tria_refine2LevelOrdering (rACproblem%RlevelInfo(i-1)%rtriangulation,&
          rACproblem%RlevelInfo(i)%rtriangulation, rACproblem%rboundary)
      CALL tria_initStandardMeshFromRaw (rACproblem%RlevelInfo(i)%rtriangulation,&
          rACproblem%rboundary)
    END DO

  END SUBROUTINE

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

!<subroutine>

  SUBROUTINE AC5_initDiscretisation (rACproblem)

!<description>
  ! This routine initialises the discretisation structure of the underlying
  ! problem and saves it to the problem structure.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  TYPE(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem
!</inputoutput>

!</subroutine>

  ! local variables
  INTEGER :: I

    ! An object for saving the domain:
    TYPE(t_boundary), POINTER :: p_rboundary

    ! An object for saving the triangulation on the domain
    TYPE(t_triangulation), POINTER :: p_rtriangulation

    ! An object for the block discretisation on one level
    TYPE(t_blockDiscretisation), POINTER :: p_rdiscretisation

    DO i=rACproblem%NLMIN,rACproblem%NLMAX
      ! Ask the problem structure to give us the boundary and triangulation.
      ! We need it for the discretisation.
      p_rboundary => rACproblem%rboundary
      p_rtriangulation => rACproblem%RlevelInfo(i)%rtriangulation

      ! Now we can start to initialise the discretisation. At first, set up
      ! a block discretisation structure that specifies the blocks in the
      ! solution vector. In this simple problem, we only have one block.
      ALLOCATE(p_rdiscretisation)
      CALL spdiscr_initBlockDiscr (p_rdiscretisation,1,&
                                    p_rtriangulation, p_rboundary)

      ! Save the discretisation structure to our local LevelInfo structure
      ! for later use.
      rACproblem%RlevelInfo(i)%p_rdiscretisation => p_rdiscretisation

      ! p_rdiscretisation%Rdiscretisations is a list of scalar
      ! discretisation structures for every component of the solution vector.
      ! Initialise the first element of the list to specify the element
      ! and cubature rule for this solution component:

!       CALL spdiscr_initDiscr_simple ( &
!                   p_rdiscretisation%RspatialDiscr(1), &
!                   EL_E011,CUB_G4X4, &
!                   p_rtriangulation, p_rboundary)

      CALL spdiscr_initDiscr_simple ( &
                  p_rdiscretisation%RspatialDiscr(1), &
                  EL_Q1,CUB_G4X4, &
                  p_rtriangulation, p_rboundary)

!     CALL spdiscr_initDiscr_simple ( &
!                  p_rdiscretisation%RspatialDiscr(1), &
!                  EL_EM30,CUB_G4X4, &
!                  p_rtriangulation, p_rboundary)

    END DO

  END SUBROUTINE

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

!<subroutine>

  SUBROUTINE AC5_doneDiscretisation (rACproblem)

!<description>
  ! Releases the discretisation from the heap.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  TYPE(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem
!</inputoutput>

!</subroutine>

  ! local variables
  INTEGER :: i

    DO i=rACproblem%NLMAX,rACproblem%NLMIN,-1
      ! Delete the block discretisation together with the associated
      ! scalar spatial discretisations....
      CALL spdiscr_releaseBlockDiscr(rACproblem%RlevelInfo(i)%p_rdiscretisation)

      ! and remove the allocated block discretisation structure from the heap.
      DEALLOCATE(rACproblem%RlevelInfo(i)%p_rdiscretisation)
    END DO

  END SUBROUTINE

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

!<subroutine>

  SUBROUTINE AC5_doneParamTriang (rACproblem)

!<description>
  ! Releases the triangulation and parametrisation from the heap.
!</description>

!<inputoutput>
  ! A problem structure saving problem-dependent information.
  TYPE(t_ACproblem), INTENT(INOUT), TARGET :: rACproblem
!</inputoutput>

!</subroutine>

  ! local variables
  INTEGER :: i

    DO i=rACproblem%NLMAX,rACproblem%NLMIN,-1
      ! Release the triangulation
      CALL tria_done (rACproblem%RlevelInfo(i)%rtriangulation)
    END DO

    ! Finally release the domain.
    CALL boundary_release (rACproblem%rboundary)

  END SUBROUTINE

END MODULE

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
implicit none
real :: x, xold, f, fprime
integer :: m

m = 5

stop
end

max_stars_repo_path: 4670 Numerical Analysis/Class Examples/modified.f90


text:
module quadrature_module

      use polynomial_module
      implicit none

      contains

!===================================================================
subroutine trapzpointsandweights(x,w,np)

	real*8,allocatable,intent(out) :: w(:)
	real*8,allocatable,intent(out) :: x(:)

	integer,intent(out) :: np
	real*8 :: dx
	integer :: i

	np=10 !for now

	dx=2.d0/(np-1)

	allocate(w(np))
	allocate(x(np))

	w(1)=0.5*dx
	x(1)=-1.0

	w(np)=0.5*dx
	x(np)=1.0

	do i=2,np-1
		w(i)=dx
		x(i)=-1.d0+(i-1)*dx
	enddo

end subroutine trapzpointsandweights
!===================================================================
subroutine gausspointsandweights(x,w,npoints)

		integer,intent(in) :: npoints
		real*8, allocatable, intent(out) :: w(:)
		real*8, allocatable, intent(out) :: x(:)

		allocate(w(npoints))
		allocate(x(npoints))

		if(npoints .eq. 1) then
			w(1)=2.0
			x(1)=0.0
		endif

		if(npoints .eq. 2) then
			w(1)=1.0
			w(2)=1.0

			x(1)=-1.0/sqrt(3.0)
			x(2)=1.0/sqrt(3.0)
		endif

		if(npoints .eq. 3) then
			w(1)=5.0/9.0
			w(2)=8.0/9.0
			w(3)=5.0/9.0

			x(1)=-sqrt(3.0/5.0)
			x(2)=0.0
			x(3)=sqrt(3.0/5.0)
		endif

		if(npoints .eq. 4) then
			w(1) = (18.0-sqrt(30.0))/36.0
			w(2) = (18.0+sqrt(30.0))/36.0
			w(3) = (18.0+sqrt(30.0))/36.0
			w(4) = (18.0-sqrt(30.0))/36.0

			x(1) = -sqrt((3.0+2.0*sqrt(6.0/5.0))/7.0)
			x(2) = -sqrt((3.0-2.0*sqrt(6.0/5.0))/7.0)
			x(3) =  sqrt((3.0-2.0*sqrt(6.0/5.0))/7.0)
			x(4) =  sqrt((3.0+2.0*sqrt(6.0/5.0))/7.0)
		endif

end subroutine gausspointsandweights
!===================================================================
subroutine gll_gen(x,w,N)

      integer, intent(in) :: N
      real*8,allocatable, intent(out) :: x(:),w(:)

      real*8 :: dleg(N+1)
      real*8 :: pi,tol,x_it,xold
      integer :: maxit,N1,i,j,k

      allocate(x(N+1))
      allocate(w(N+1))

      pi = acos(-1.)

      tol = 1d-15

      N1 = N+1

      maxit = 1000   ! max iterations for newton-raphson

      x(1) = -1.d0
      x(N1) = 1.d0

      do i = 1, N+1

        x_it = -cos(pi * float(i-1) / N) ! initial guess - chebyshev points

        do j = 1, maxit

    	  xold = x_it
          dleg(1) = 1.d0
          dleg(2) = x_it

	  do k = 2,N

            dleg(k+1) = (  (2.d0*float(k) - 1.d0) * dleg(k) * x_it &
                    - (float(k)-1.d0)*dleg(k-1) ) / float(k)

          enddo

          x_it = x_it - ( x_it * dleg(N1) - dleg(N) ) / &
                (float(N1) * dleg(N1) )

          if (abs(x_it - xold) .lt. tol) goto 10

        enddo
        write(*,*) 'max itertions reached'
10      continue

        x(i) = x_it
        w(i) = 2.d0 / (float(N * N1) * dleg(N1)**2 )

      enddo

      return

end subroutine gll_gen
!===================================================================
 subroutine gl_gen(x,w,n)
      integer n
      double precision x1,x2
      double precision eps
      parameter (eps=3.d-14)
      integer i,j,m
      double precision p1,p2,p3,pp,xl,xm,z,z1
      real*8,allocatable :: x(:),w(:)

      allocate(x(n))
      allocate(w(n))

      m=(n+1)/2

      x1 = -1.d0
      x2 = +1.d0

      xm=0.5d0*(x2+x1)
      xl=0.5d0*(x2-x1)
      do 12 i=1,m
        z=cos(3.141592653589793d0*(i-.25d0)/(n+.5d0))
1       continue
          p1=1.d0
          p2=0.d0
          do 11 j=1,n
            p3=p2
            p2=p1
            p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
11        continue
          pp=n*(z*p1-p2)/(z*z-1.d0)
          z1=z
          z=z1-p1/pp
        if(abs(z-z1).gt.eps)goto 1
        x(i)=xm-xl*z
        x(n+1-i)=xm+xl*z
        w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
        w(n+1-i)=w(i)
12    continue
      return
      end
!===================================================================
subroutine globattopointsandweights(x,w,npoints)

		real*8,allocatable,intent(out) :: w(:),x(:)
		integer, intent(in) :: npoints

		allocate(w(npoints))
		allocate(x(npoints))

		if(npoints .eq. 2) then
			w(1)=1.0
			w(2)=1.0

			x(1)=-1.0
			x(2)=1.0
		endif

		if(npoints .eq. 3) then

			w(1)=1.0/3.0
			w(2)=4.0/3.0
			w(3)=1.0/3.0

			x(1)=-1.0
			x(2)=0.0
			x(3)=1.0
		endif

		if(npoints .eq. 4) then

			w(1)=1.0/6.0
			w(2)=5.0/6.0
			w(3)=5.0/6.0
			w(4)=1.0/6.0

			x(1)=-1.0
			x(2)=-sqrt(1.0/5.0)
			x(3)=sqrt(1.0/5.0)
			x(4)=1.0
		endif

		if(npoints .eq. 5) then

			w(1)=1.0/10.0
			w(2)=49.0/90.0
			w(3)=32.0/45.0
			w(4)=49.0/90.0
			w(5)=1.0/10.0

			x(1)=-1.0
			x(2)=-sqrt(3.0/7.0)
			x(3)=0.0
			x(4)=sqrt(3.0/7.0)
			x(5)=1.0
		endif

end subroutine globattopointsandweights
!===================================================================
subroutine getuniformlagrangepoints(n,x)

		integer :: n
		real*8 :: x(:)

		integer :: i
		real*8 :: delx
		real*8 :: xmin,xmax

		xmin = -1.d0
		xmax =  1.d0
		delx= (xmax-xmin)/n

		do i=1,n+1
			x(i)=xmin+(i-1)*delx
		enddo

end subroutine getuniformlagrangepoints
!===================================================================
subroutine convolvetrapz(p1,p2,a,b,npoints,answer)

	      type(polynomial), intent(in) :: p1,p2
	      real*8, intent(in) :: a,b
	      integer, intent(in) :: npoints
	      real*8 , intent(out) :: answer

	      integer :: i
	      real*8 :: delx,integval
	      real*8 :: h1,h2

	      delx=(b-a)/(npoints-1)
	      integval=0.0

	      do i=1,npoints-1

	      	h1 = polyfindval(p1,a+(i-1)*delx)*polyfindval(p2,a+(i-1)*delx)
		h2 = polyfindval(p1,a+i*delx)*polyfindval(p2,a+i*delx)

		integval = integval + 0.5*delx*(h1+h2)

	      enddo

	      answer=integval

end subroutine convolvetrapz
!===================================================================

end module quadrature_module

max_stars_repo_name: NREL/VirtualEngineering

max_issues_repo_name: NREL/VirtualEngineering

max_forks_repo_name: NREL/VirtualEngineering


text:
C$Procedure ZZSECPRT ( Calculate dot terms for DPSPCE perturbation )

      SUBROUTINE ZZSECPRT ( ISYNFL, DG   , DEL   , XNI  ,
     .                      OMEGAO, ATIME, OMGDOT, XLI  ,
     .                      XFACT , XLDOT, XNDOT , XNDDT )

C$ Abstract
C
C    Routine to calculate the dot terms for the secular perturbation
C    of a vehicle.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     SECULAR PERTURBATION
C
C$ Declarations

      IMPLICIT NONE

      INTEGER               ISYNFL

      DOUBLE PRECISION      OMEGAO
      DOUBLE PRECISION      ATIME
      DOUBLE PRECISION      OMGDOT
      DOUBLE PRECISION      XLI
      DOUBLE PRECISION      DG     ( 10 )
      DOUBLE PRECISION      XFACT
      DOUBLE PRECISION      XLDOT
      DOUBLE PRECISION      XNI
      DOUBLE PRECISION      DEL    (  3 )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     ISYNFL     I   Resonance flag
C     DG         I   Parameter array
C     DEL        I   Parameter array of delta values
C     XNI        I   An intermediate linear term passed from the
C                    calling routine
C     OMEGAO     I   Original argument of perigee
C     ATIME      I   An intermediate time term passed from the calling
C                    routine
C     OMGDOT     I   Time rate of change of argument of perigee
C     XLI        I   An intermediate angular term passed from the
C                    calling routine
C     XFACT      I   The value BFACT - XNQ
C     XLDOT      O   Time rate of change of XL
C     XNDOT      O   Time rate of change of XN
C     XNDDT      O   Time rate of change of XNDOT
C
C$ Detailed_Input
C
C     ISYNFL      is the flag used to indicate the need for resonance
C                 calculations.
C
C     DG          is the parameter array replacing the Dxxxx values.
C
C     DEL         is the parameter array replacing DEL1, DEL2 and DEL3.
C
C     XNI         is an intermediate linear term passed from the main
C                 term for the calculation of XLDOT = XNI  +  XFACT
C
C     OMEGAO      is the original value for the argument of perigee.
C
C     ATIME       is an intermediate time term passed from the main
C                 routine used to calculate the time dependent
C                 argument of perigee term XOMI
C
C     OMGDOT      is the time derivative of the argument of the perigee.
C
C     XLI         is an intermediate angular term
C
C     XFACT       is the value BFACT - XNQ calculated in ZZDPINIT
C
C$ Detailed_Output
C
C     XLDOT       time derivative of the XL term.
C
C     XNDOT       time derivative of the XN term.
C
C     XNDDT       second time derivative of XN, time derivative of the
C                 time derivative.
C
C$ Parameters
C
C     None.
C
C$ Files
C
C     None.
C
C$ Exceptions
C
C     None.
C
C$ Particulars
C
C     This subroutine was constructed from a section of code in ZZDPSEC
C     in the original Spacetrack 3 report.  The code block was called
C     using a set of conditional GO TO's.  The block has been written as
C     this subroutine to improve clarity and maintainability and to
C     conform to the NAIF style standard.
C
C$ Examples
C
C     None needed.
C
C$ Restrictions
C
C     1)  This routine should be called only by ZZDPSEC as part of the
C     DPSPCE subroutine package.  It has no other use.
C
C$ Author_and_Institution
C
C     E.D. Wright      (JPL)
C
C$ Literature_References
C
C     Spacetrack 3 report.
C
C$ Version
C
C-   SPICELIB Version 1.1.0, 24-MAR-1999 (EDW)
C
C       Correction made to format of Version descriptions.
C       Name of variable OMGDT changed to OMGDOT to be consistent
C       with name useage in other deep space two line elements
C       routines.
C
C-   SPICELIB Version 1.0.0, MAY-19-1997 (EDW)
C
C
C-&

C$ Index_Entries
C
C     perturbed dot terms
C
C-&

C
C     Local variables.
C
      DOUBLE PRECISION      X2LI
      DOUBLE PRECISION      X2OMI
      DOUBLE PRECISION      XNDDT
      DOUBLE PRECISION      XNDOT
      DOUBLE PRECISION      XOMI

      DOUBLE PRECISION      G22
      PARAMETER           ( G22    = 5.7686396D0  )

      DOUBLE PRECISION      G32
      PARAMETER           ( G32    = 0.95240898D0 )

      DOUBLE PRECISION      G44
      PARAMETER           ( G44    = 1.8014998D0  )

      DOUBLE PRECISION      G52
      PARAMETER           ( G52    = 1.0508330D0  )

      DOUBLE PRECISION      G54
      PARAMETER           ( G54    = 4.4108898D0  )

      DOUBLE PRECISION      FASX2
      PARAMETER           ( FASX2  = 0.13130908D0 )

      DOUBLE PRECISION      FASX4
      PARAMETER           ( FASX4  = 2.8843198D0  )

      DOUBLE PRECISION      FASX6
      PARAMETER           ( FASX6  = 0.37448087D0 )

C
C     Calculate the dot terms with respect to the state of the
C     resonance flag.
C

      IF ( ISYNFL .EQ. 0 ) THEN

C
C        Resonance flag set.
C

         XOMI  = OMEGAO + OMGDOT * ATIME
         X2OMI = XOMI   + XOMI
         X2LI  = XLI    + XLI

         XNDOT =  DG( 1  )*DSIN( X2OMI + XLI  - G22 )
     .          + DG( 2  )*DSIN( XLI          - G22 )
     .          + DG( 3  )*DSIN( XOMI  + XLI  - G32 )
     .          + DG( 4  )*DSIN(-XOMI  + XLI  - G32 )
     .          + DG( 5  )*DSIN( X2OMI + X2LI - G44 )
     .          + DG( 6  )*DSIN( X2LI         - G44 )
     .          + DG( 7  )*DSIN( XOMI  + XLI  - G52 )
     .          + DG( 8  )*DSIN(-XOMI  + XLI  - G52 )
     .          + DG( 9  )*DSIN( XOMI  + X2LI - G54 )
     .          + DG( 10 )*DSIN(-XOMI  + X2LI - G54 )

         XNDDT =  DG( 1  )*DCOS( X2OMI + XLI  - G22 )
     .          + DG( 2  )*DCOS( XLI          - G22 )
     .          + DG( 3  )*DCOS( XOMI  + XLI  - G32 )
     .          + DG( 4  )*DCOS(-XOMI  + XLI  - G32 )
     .          + DG( 7  )*DCOS( XOMI  + XLI  - G52 )
     .          + DG( 8  )*DCOS(-XOMI  + XLI  - G52 )
     .   +2.D0*(  DG( 5  )*DCOS( X2OMI + X2LI - G44 )
     .          + DG( 6  )*DCOS( X2LI         - G44 )
     .          + DG( 9  )*DCOS( XOMI  + X2LI - G54 )
     .          + DG( 10 )*DCOS( XOMI  + X2LI - G54 ))

      ELSE

C
C        Resonance flag not set
C

         XNDOT =  DEL(1)*DSIN(        XLI - FASX2   )
     .          + DEL(2)*DSIN( 2.D0*( XLI - FASX4 ) )
     .          + DEL(3)*DSIN( 3.D0*( XLI - FASX6 ) )

         XNDDT =      DEL(1)*DCOS(        XLI - FASX2   )
     .         + 2.D0*DEL(2)*DCOS( 2.D0*( XLI - FASX4 ) )
     .         + 3.D0*DEL(3)*DCOS( 3.D0*( XLI - FASX6 ) )

      END IF

      XLDOT = XNI  +  XFACT
      XNDDT = XNDDT * XLDOT

C
C     Hi!  What are you doing way down here?  Did you bring pizza?
C

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
	SUBROUTINE G2T_SQUAD ( iret )
C************************************************************************
C* G2T_SQUAD								*
C*									*
C* This subroutine sets quadrant info to the common.  The quad index	*
C* is from the G2T_ZONE.TBL.  For most four quad subzones, the first 	*
C* is NE quad, and the 2nd/3rd/4th quads are SE/SW/NW zones.  The quad	*
C* info is output in the EXCEPTION part of text.			*
C*									*
C* G2T_SQUAD ( IRET )							*
C*									*
C* Input parameters:							*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -6 = cannot find the subzone	*
C**									*
C* Log:									*
C* T. Lee/SAIC		 1/07	Created					*
C* T. Lee/SAIC		 7/07	Updated for multiple zones		*
C* T. Lee/SAIC		11/07	Updated for TPC zones			*
C* T. Lee/SAIC		12/07	Updated for '_0_' PORTION flag (porflg)	*
C************************************************************************
	INCLUDE		'goftxt.cmn'
	LOGICAL		found
C------------------------------------------------------------------------
	iret = 0
C
	DO ii = 1, numzon
	  found = .false.
	  DO jj = 1, nsubzn ( ii )
C
C*	    Blank for entire zone.
C
	    IF  ( .not. found .and.
     +	        ( INDEX ( bndtag (ii,jj), '_ENTIRE' ) ) .ne. 0 )  THEN
		found = .true.
		squad ( ii, jj ) = ' '
C
C*	      First quadrant (except of S_OF_NE to avoid NE error).
C

	      ELSE IF ( INDEX ( bndtag (ii,jj), '_100_70W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), 'MAINE_N' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), 'MAINE_S' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_EOF120' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_NCNTRL' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_CNTRL' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_68_70W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_31N_E' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_31N_W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_E_68W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_E_70W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_W_70W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_W_1000' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_NESA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_SESA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_SWSA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_CSA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_ESA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_NSA' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_NE' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_SE' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_SW' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_NW' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_E' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_S' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_W' ) .ne. 0 .or.
     +			INDEX ( bndtag (ii,jj), '_N' ) .ne. 0 )  THEN
C
C*		Set 'OVER' flags.
C
		IF ( INDEX ( bndtag (ii,jj), '_0_' ) .ne. 0 )  THEN
		    porflg ( ii, jj ) = .false.
		  ELSE
		    porflg ( ii, jj ) = .true.
		END IF
C
C*		Set ID names.
C
		IF ( INDEX ( bndtag (ii,jj), '_100_70W' ) .ne. 0 )  THEN
		    squad ( ii, jj ) = 'FROM 100NM OF THE COAST TO 70W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), 'MAINE_N' ) .ne. 0 )
     +			    THEN
		    squad ( ii, jj ) = 'NORTH'
		  ELSE IF ( INDEX ( bndtag (ii,jj), 'MAINE_S' ) .ne. 0 )
     +			    THEN
		    squad ( ii, jj ) = 'SOUTH'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_EOF120' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST OF 120W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_NCNTRL' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NORTH CENTRAL'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_CNTRL' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'CENTRAL'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_68_70W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'BETWEEN 68W AND 70W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_31N_E' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_31N_W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'WEST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_E_68W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST OF 68W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_E_70W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST OF 70W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_W_70W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'WEST OF 70W'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_W_1000' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'WEST OF 100NM OF THE COAST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_NESA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NE'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_SESA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'SE'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_SWSA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'SW'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_CSA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'CENTRAL'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_ESA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_NSA' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NORTH'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_NE' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NE'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_SE' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'SE'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_SW' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'SW'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_NW' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NW'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_E' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'EAST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_W' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'WEST'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_S' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'SOUTH'
		  ELSE IF ( INDEX ( bndtag (ii,jj), '_N' )
     +				.ne. 0 )  THEN
		    squad ( ii, jj ) = 'NORTH'
		  ELSE
		    iret = -6
		END IF
	      ELSE
		iret = -6
	    END IF
	  END DO
	END DO
C
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine interp (matl,dnt,temp,opacity)
c***********************************************************************
c     calculates material opacities as a function of temperature and
c     density by interpolating in the llnl library.
c     arguments:
c     matl: material number
c     temp: temperature (units?),
c     dnt :  density (units?)
c     returns: 13 group opacity vector into array opacity.
c***********************************************************************
      implicit double precision (a-h,o-z)

      common /opactab/ x(14,32,13,4), t(32,2), d(14)
      dimension opacity(13)
c     scatter, shared opacity

      ityp = 1
      if (matl .eq. 2 .or. matl .eq. 4) ityp = 2

c--------------------find bounding temperatures
      do 20 it = 1,32
         if (temp .lt. t(it,ityp)) goto 30
 20   continue
 30   tmax = t(it,ityp)
      it1 = it-1
      tmin = t(it1,ityp)
      trat = (temp - tmin)/(tmax - tmin)

c--------------------find bounding densities
      do 40 id = 1,14
c         if (dnt .lt. d(id)) goto 50   Replaced with following due to
c addressing problems
         if (dnt .le. d(id)) goto 50
 40   continue
 50   dmax = d(id)
      id1 = id-1
      dmin = d(id1)
      drat = (dnt - dmin)/(dmax - dmin)

c--------------------interpolate
      do 60 ie = 1,13
         xt = x(id1,it,ie,matl) + (x(id,it,ie,matl)-x(id1,it,ie,matl))
     |        * drat
         xt1 = x(id1,it1,ie,matl)
     |        + (x(id,it1,ie,matl)-x(id1,it1,ie,matl))*drat
         opacity (ie) = exp(xt1 + (xt-xt1)*trat)
 60   continue

      return
      end

max_stars_repo_name: gmlueck/llvm-test-suite

max_issues_repo_name: gmlueck/llvm-test-suite

max_forks_repo_name: gmlueck/llvm-test-suite


! { dg-do compile }
program main
  type t
    integer :: x
  end type t
  type(t), dimension(2) :: a, b
  a(1)%x = 1
  a(2)%x = 2
  b = eoshift(a,1) ! { dg-error "Missing 'boundary' argument to 'eoshift' intrinsic" }
  print *,b%x
end program main


text:
! RUN: %flang_fc1 -fsyntax-only -fdebug-dump-symbols %s 2>&1 | FileCheck %s
! Verify that the closure of EQUIVALENCE'd symbols with any DATA
! initialization produces a combined initializer, with explicit
! initialization overriding any default component initialization.
! CHECK: .F18.0, SAVE (CompilerCreated) size=8 offset=0: ObjectEntity type: INTEGER(4) shape: 1_8:2_8 init:[INTEGER(4)::456_4,234_4]
! CHECK: ja (InDataStmt) size=8 offset=0: ObjectEntity type: INTEGER(4) shape: 1_8:2_8
! CHECK-NOT: x0, SAVE size=8 offset=8: ObjectEntity type: TYPE(t1) init:t1(m=123_4,n=234_4)
! CHECK: x1 size=8 offset=16: ObjectEntity type: TYPE(t1) init:t1(m=345_4,n=234_4)
! CHECK: x2 size=8 offset=0: ObjectEntity type: TYPE(t1)
! CHECK-NOT: x3a, SAVE size=8 offset=24: ObjectEntity type: TYPE(t3) init:t3(t2=t2(k=567_4),j=0_4)
! CHECK: x3b size=8 offset=32: ObjectEntity type: TYPE(t3) init:t3(k=567_4,j=678_4)
! CHECK: Equivalence Sets: (x2,ja(1)) (.F18.0,x2)
type :: t1
  sequence
  integer :: m = 123
  integer :: n = 234
end type
type :: t2
  integer :: k = 567
end type
type, extends(t2) :: t3
  integer :: j ! uninitialized
end type
type(t1), save :: x0 ! not enabled
type(t1) :: x1 = t1(m=345)
type(t1) :: x2
type(t3), save :: x3a ! not enabled
type(t3) :: x3b = t3(j=678)
integer :: ja(2)
equivalence(x2, ja)
data ja(1)/456/
end

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
! this function is no longer used, and will be removed soon
      function num_parthds_stochy()
      integer:: number_of_openMP_threads
      character(2) :: omp_threads
      integer :: stat
      call get_environment_variable("OMP_NUM_THREADS",omp_threads)
      read(omp_threads,*,iostat=stat)number_of_openMP_threads
      num_parthds_stochy = number_of_openMP_threads
      return
      end

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


program sobol_index_test

  use mo_kind,        only: dp,i4, i8
  use mo_ansi_colors, only: color, c_red, c_green
  use mo_sobol,       only: sobol
  use mo_model,       only: getrange, model
  use mo_sobol_index, only: sobol_index

  implicit none

  integer(i4), parameter                 :: npara = 3                       ! number of parameter
  integer(i4), parameter                 :: nsets = 1000                    ! number of parametersets
  integer(i4), parameter                 :: nx = 100                        ! number of variable values
  integer(i8)                            :: skip = 30000_i8                 ! used for generating parameter sets
  !                                                                         ! via sobol sequences

  ! Samples and Model Output
  real(dp),    dimension(npara,2)        :: ParaRange                       ! array: Parameter ranges (Min, Max)
  real(dp),    dimension(nx)             :: x                               ! array: variable values

  ! Sobol index
  real(dp),    dimension(2*npara)        :: sample_sobol                    ! sobol sequence
  real(dp),    dimension(npara)          :: sample_a                        ! parameter sample A
  real(dp),    dimension(npara)          :: sample_b                        ! parameter sample B
  real(dp),    dimension(npara)          :: sample_ci                       ! parameter sample Ci
  real(dp),    dimension(nsets,nx)       :: ya, yb                          ! model output A and B,
  !                                                                         !    i.e. time series with nx time points
  real(dp),    dimension(nsets,npara,nx) :: yc                              ! model output C(i), i=1,npara,
  !                                                                         !    i.e. time series with nx time points
  real(dp),    dimension(npara,nx)       :: si                              ! Sobol index (main effect)  per time point and param.
  !                                                                         !    based on ya, yb, and yc
  real(dp),    dimension(npara,nx)       :: sti                             ! Sobol index (total effect) per time point and param.
  !                                                                         !    based on ya, yb, and yc
  real(dp),    dimension(npara,2)        :: smean                           ! Mean SI and STI
  real(dp),    dimension(npara,2)        :: wmean                           ! Variance weighted mean SI and STI

  integer(i8) :: seed

  ! Dummy variables
  integer(i4)                            :: set, para, i
  logical                                :: isgood = .true.

  ParaRange = GetRange()

  ! ------------------------------------------------------------------
  ! variable values, e.g. time points
  ! ------------------------------------------------------------------
  do i=1,nx
     x(i) = real(i-1,dp) * 1.0_dp + 50.0_dp
  end do

  ! ------------------------------------------------------------------
  ! Generating model outputs for parameter sets A, B, and C(i)
  ! ------------------------------------------------------------------
  seed = 0_i8
  call sobol(int(2*npara,i8), seed, sample_sobol)
  do set = 1,nsets
     ! sobol sequence
     call sobol(int(2*npara,i8), skip, sample_sobol)
     sample_a = sample_sobol(1:npara)
     sample_b = sample_sobol(npara+1:2*npara)
     ! scaling
     sample_a = ParaRange(:,1) + sample_a * (ParaRange(:,2) - ParaRange(:,1))
     sample_b = ParaRange(:,1) + sample_b * (ParaRange(:,2) - ParaRange(:,1))

     ya(set,:) = model(sample_a,x)
     yb(set,:) = model(sample_b,x)

     do para=1,npara
        sample_ci       = sample_b
        sample_ci(para) = sample_a(para)
        yc(set,para,:)  = model(sample_ci,x)
     end do
  end do

  ! ------------------------------------------------------------------
  ! Calculate Sobol index: Main Effect (Si) and Total Effect (STi)
  ! ------------------------------------------------------------------
  call sobol_index(ya, yb, yc, si, sti, method=6, smean=smean, wmean=wmean)

  ! ------------------------------------------------------------------
  ! Printing and checking
  ! ------------------------------------------------------------------
  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A24,F5.1,A8)') '     Sobol index   (x = ',x(3),')       '
  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A40)')         'para  SI         STI                    '
  do para=1,npara
     write(*,'(I4,A2,F7.4,A4,F7.4)') para, '   ', si(para,3),'  ', sti(para,3)
  end do
  write(*,'(A40)') '----------------------------------------'
  write(*,*) ' '

  if (nint(si(1,3)*10000._dp)      .ne.  3226_i4)  isgood = .false.
  if (nint(si(2,3)*10000._dp)      .ne.  1404_i4)  isgood = .false.
  if (nint(si(3,3)*10000._dp)      .ne. -  58_i4)  isgood = .false.
  if (nint(sti(1,3)*10000._dp)     .ne.  7826_i4)  isgood = .false.
  if (nint(sti(2,3)*10000._dp)     .ne.  1490_i4)  isgood = .false.
  if (nint(sti(3,3)*10000._dp)     .ne.  5348_i4)  isgood = .false.

  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A40)')         '     Mean Sobol index                   '
  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A40)')         'para  SI         STI                    '
  do para=1,npara
     write(*,'(I4,A2,F7.4,A4,F7.4)') para, '   ', smean(para,1),'  ', smean(para,2)
  end do
  write(*,'(A40)') '----------------------------------------'
  write(*,*) ' '

  if (nint(smean(1,1)*10000._dp)     .ne.  2417_i4)  isgood = .false.
  if (nint(smean(2,1)*10000._dp)     .ne.  3689_i4)  isgood = .false.
  if (nint(smean(3,1)*10000._dp)     .ne.    13_i4)  isgood = .false.
  if (nint(smean(1,2)*10000._dp)     .ne.  5784_i4)  isgood = .false.
  if (nint(smean(2,2)*10000._dp)     .ne.  3820_i4)  isgood = .false.
  if (nint(smean(3,2)*10000._dp)     .ne.  3980_i4)  isgood = .false.

  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A40)')         '     Var. weighted mean Sobol index     '
  write(*,'(A40)')         '----------------------------------------'
  write(*,'(A40)')         'para  SI         STI                    '
  do para=1,npara
     write(*,'(I4,A2,F7.4,A4,F7.4)') para, '   ', wmean(para,1),'  ', wmean(para,2)
  end do
  write(*,'(A40)') '----------------------------------------'
  write(*,*) ' '

  if (nint(wmean(1,1)*10000._dp)     .ne.  2319_i4)  isgood = .false.
  if (nint(wmean(2,1)*10000._dp)     .ne.  3963_i4)  isgood = .false.
  if (nint(wmean(3,1)*10000._dp)     .ne.    20_i4)  isgood = .false.
  if (nint(wmean(1,2)*10000._dp)     .ne.  5535_i4)  isgood = .false.
  if (nint(wmean(2,2)*10000._dp)     .ne.  4102_i4)  isgood = .false.
  if (nint(wmean(3,2)*10000._dp)     .ne.  3813_i4)  isgood = .false.

  if (isgood) then
     write(*,*) 'mo_sobol_index ', color('o.k.', c_green)
  else
     write(*,*) 'mo_sobol_index ', color('failed', c_red)
  end if

end program sobol_index_test


text:
      SUBROUTINE FINDGB
C
C      SUBROUTINE FINDGB CALCULATES THE GRAVITY GRADIENT MATRIX,
C             SUNLINE VECTOR FOR INTERNAL ORBIT GENERATOR, AND
C             THE SOLAR ILLUMINATION FACTOR
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 J2,J3,J4,J22
C
C
      COMMON/CNBODY/ J2,J3,J4,J22,ZJ20,ZMU,WWO,FLAT,AEARTH
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/CSTVAL/ TSTART,ZL0(10),ZL1(10),ZL2(10),ZLA(10)
C
      COMMON/IMAIN1/ IDATE,LSAVE,INOPT,IPLOT,NUMEQS,IPLTPE,IORB,ITAPE
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/ISHADE/ IPLANS,ISATSH,IWRTTF
C
      COMMON/OUTTWO/ SOLILL,EPSERR
C
      COMMON/PRESUR/ DTOO,POO
C
      COMMON/RPOOL1/ RHOK(10),TIME,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
      COMMON/RPOOL2/ PO,SD(3),DAN(3,10),DBN(3,10),CFMT(3,3),DIY1(3),
     .               SD1(3),DT1,P1,AERO,DTO1,YIZK(3),PO1
C
      COMMON /RSUNCL/ WE,TVER,ECLPTC
C
      COMMON/VECTRS/ XSAT(3),XSATDT(3),AD(3)
C
      COMMON/ORBNEW/ROD(3),VOD(3),TST
C
      DIMENSION XSAT1(3),XSAT2(3)
C
C
      DATA CONV/7.905355D0/
      IF(IORB.GT.0) GO TO 5
      IF(TIME-TSTART.GT.0) GO TO 4
C
      DO 3 I=1,3
      ROD(I)=XSAT(I)
    3 VOD(I)=XSATDT(I)
      TST=TSTART
    4 CALL XFIND(ROD,VOD,ZMU,TST   ,TIME,XSAT,XSATDT)
      GO TO 559
 5    CONTINUE
      TORB=TIME
      CALL TCNVRT(DDATE,MM,SEC,TLAST,TORB,2)
      CALL DTAPRE(ITAPE,IERR,DDATE,SEC,RX,RY,RZ,VX,VY,VZ)
      XSAT(1)=RX*AEARTH
      XSAT(2)=RY*AEARTH
      XSAT(3)=RZ*AEARTH
      XSATDT(1)=VX*CONV
      XSATDT(2)=VY*CONV
      XSATDT(3)=VZ*CONV
      IF(IERR.LT.3) GO TO 2000
      IF(IERR.LT.6) GO TO 3000
 559  CONTINUE
      RADSQR=0
      DO 10 I=1,3
   10 RADSQR=RADSQR + XSAT(I)*XSAT(I)
      RADIUS=DSQRT(RADSQR)
C
      DO 20 I=1,3
      XSAT1(I)=XSAT(I)/RADIUS
   20 XSAT2(I)=XSAT1(I)*XSAT1(I)
C
      DO 30 I=1,3
      DO 30 J=1,3
   30 GBAR(I,J)=0.D0
C
      IF(IGRAV.EQ.0) GO TO 60
C
      RAD3RD=RADIUS*RADSQR
      RAD4TH=RADIUS*RAD3RD
      RAD5TH=RADIUS*RAD4TH
C
      ZJ2OMU=-(ZJ20*ZMU*AEARTH*AEARTH)/(2.D0*RAD5TH)
      ZMUR3=-ZMU/RAD3RD
C
      DO 40 I=1,2
C
      GBAR(I,I)=  ZMUR3*(1.D0 - 3.D0*XSAT2(I))
     .          + ZJ2OMU*(3.D0 - 15.D0*(XSAT2(I) + XSAT2(3))
     .                  + 105.D0*XSAT2(I)*XSAT2(3))
C
   40 GBAR(3,I)=- 3.D0*ZMUR3*XSAT1(I)*XSAT1(3)
     .          + ZJ2OMU*(105.D0*XSAT1(I)*XSAT1(3)*XSAT2(3)
     .                  - 45.D0*XSAT1(I)*XSAT1(3))
C
C
C
      X21=XSAT1(1)*XSAT1(2)
      GBAR(2,1)=- 3.D0*ZMUR3*X21 + ZJ2OMU*(105.D0*X21*XSAT2(3)
     .          - 15.D0*X21)
C
      X34=XSAT2(3)*XSAT2(3)
C
      GBAR(3,3)=  ZMUR3*(1.D0 - 3.D0*XSAT2(3))
     .          + ZJ2OMU*(9.D0 - 90.D0*XSAT2(3) + 105.D0*X34)
C
C
      DO 50 I=1,3
      DO 50 J=I,3
   50 GBAR(I,J)=GBAR(J,I)
C
C                       SOLAR ILLUMINATION
C
   60 COSC=0.D0
C
      EPS=ECLPTC*RADIAN
      CEPS=DCOS(EPS)
      SEPS=DSIN(EPS)
      WS=WE*(TIME-TVER)*RADIAN
      SINWS=DSIN(WS)
      COSWS=DCOS(WS)
      SD(1)=-COSWS
      SD(2)=-SINWS*CEPS
      SD(3)=-SINWS*SEPS
C
   65 DO 70 I=1,3
   70 COSC=COSC + SD(I)*XSAT1(I)
C
      IF(COSC.LE.0.D0) GO TO 80
C
      SINC=DSQRT(1.D0 - COSC*COSC)
      RADCOS=RADIUS*COSC
      RADSIN=RADIUS*SINC - AEARTH
C
      H1=RADSIN + 0.00461D0*RADCOS
C
      IF(H1.LE.0.D0) GO TO 90
C
      H2=RADSIN - 0.004695D0*RADCOS
C
      IF(H2.GE.0.D0) GO TO 80
C
      H3=0.004653D0*RADCOS
      H4=1.D0 + H2/H3
C
      SOLILL=1.D0 - DARCOS(H4)/PI + H4/PI*DSQRT(1.D0 - H4*H4)
      GO TO 100
   80 SOLILL=1.D0
      GO TO 100
   90 SOLILL=0.D0
  100 DTO=SOLILL*DTOO
      PO=SOLILL*POO
      DT1=0.D0
C
      IF(IPLANS.EQ.0.AND.ISATSH.EQ.0) RETURN
      CALL SHADEP(SOLILL)
C
      RETURN
 2000 WRITE(6,2001)
 2001 FORMAT(' I/O ERROR ON EPHERMIS FILE')
 3000 WRITE(6,3001)
 3001 FORMAT(' INFORMATION NOT AVAILABLE FOR REQUESTED TIME')
      RETURN
      END

ext: for


text:
C
C $Id: nasort.f,v 1.4 2008-07-27 03:11:53 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 NASORT(XA,IP,N)
C
C  Given an array of REAL numbers of dimension N in XA,
C  this routine returns a permutation vector IP such that
C
C      XA(IP(I)) .LE. XA(IP(J))
C           for all I,J such that  1. LE. I .LE. J .LE. N .
C
C  This is an implementation of the shellsort (the details of this
C  algorithm can be found in the book "Algorithms" by Robert Sedgewick).
C
      INTEGER  H,I,N,IP(N),IP1,IP2,IPH,J,M
      REAL     XA(N)
C
C  Check on the input argument N.
C
      IF (N .LT. 0) THEN
        CALL SETER(' ASCSRT -- N TOO SMALL',1,1)
      ELSE IF (N .LE. 1) THEN
        IP(1) = 1
        RETURN
      ENDIF
C
C  Calculate the largest H less than N using the recursive definition
C               H(L+1) = 3*H(L)+1
C
      H = 1
   10 CONTINUE
      H = 3*H+1
      IF (H .LT. N) GO TO 10
      H = (H-1)/3
C
C  Initialize permutation vector.
C
      DO 20 I=1,N
      IP(I) = I
   20 CONTINUE
C
C  Main loop.
C
   30 CONTINUE
      M = N-H
      DO 50 J=1,M
      I = J
   40 CONTINUE
      IPH = I+H
      IP1 = IP(I)
      IP2 = IP(IPH)
      IF (XA(IP1) .LE. XA(IP2)) GO TO 50
C
C  Interchange.
C
      IP(I)   = IP2
      IP(IPH) = IP1
C
C  Move interchanged elements up the line.
C
      I = I-H
      IF (I .GE. 1) GO TO 40
   50 CONTINUE
C
C  Decrease H, continue sort as long as H is greater than 0.
C
      H = (H-1)/3
      IF (H .LT. 1) RETURN
      GO TO 30
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! Update one interger `i` at the first rank and send it to the last rank in a
! standard send-receive pair.
!
! Compile:
!
!   $ mpifort -g -O0 -Wall -Wextra -Wpedantic -fcheck=all -fbacktrace \
!     send_recv.f90
!
! Run with 2 processes (-np 2), then with 3 or more:
!
!   $ mpirun -np 2 --oversubscribe ./a.out
program send_recv
    use mpi_f08
    implicit none
    integer :: my_rank, n_ranks, first, last
    integer :: i = 0
    type(MPI_Status) :: status

    call MPI_Init()

    call MPI_Comm_size(MPI_COMM_WORLD, n_ranks)
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank)
    first = 0
    last  = n_ranks - 1

    print "(2(a, i0))", "Before: rank=", my_rank, ", i=", i

    if (my_rank == first) then
        i = 999
        call MPI_Send(i, 1, MPI_INTEGER, last, 0, MPI_COMM_WORLD)
    else if (my_rank == last) then
        call MPI_Recv(i, 1, MPI_INTEGER, first, 0, MPI_COMM_WORLD, status)
    end if

    print "(2(a, i0))", "After:  rank=", my_rank, ", i=", i

    call MPI_Finalize()
end program send_recv


text:
A new owner bought the restaurant and they changed their name to Jade Cafe in early 2008.  They introduced a new menu that appeared to have more choices and planned on expanding their lunch buffet to the weekends.  The restaurant closed abruptly after this change in ownership and by May 2008 the phone had been disconnected and the mail slot taped shut.  In June 2008, a banner appeared over the door that said Coming Soon  Red Orchid Restaurant.  THE RED ORCHID IS OPEN AS OF 6/13/08.

ShangriLa was a Restaurants restaurant featuring a fairly typical Chinese food menu with reasonable prices and large portions. Like many Chinese restaurants, they offer three family dinners (Peking Dinner $6.95, Szechuan Dinner $7.25, Mandarin Dinner $8.55), as well as individual entrees. They also have a Buffets Lunch Buffet ($5.99+tax) MonFri 11:30AM3:00PM. The staff is friendly and welcoming and better yet, the food tastes good.

Menu

20050409 11:58:29 nbsp I just thought Id take this opportunity to vent about the shopping center on 8th street.  WOW, that is the most hideous looking shopping center.  Ive seen such colors used well in other places, but that place looks like a complete monstrosity. There is just no attempt at balance.  Absolutely none. Users/JaimeRaba

20050517 21:42:16 nbsp They are very friendly, but largely due to the location, they dont seem to get a lot of business. As a result, the buffet gets kind of nasty by 1:30PM.  Users/CraigBrozinsky

20050705 19:05:25 nbsp This place has great, authentic Chinese food, perhaps the best Davis.  Their fish and fresh vegetable dishes are excellent, but ask the friendly owner/operator what is fresh and it will be a sure delight.  The prices are very reasonable but the ambiance could use some work. Users/TeddyWong

20051016 11:03:57 nbsp The food is some of the best ive had in davis, and the service was top. The place is not that busy, so service is really good. They have a hard time doing business because theyre not in a high profile location, so go enjoy a good meal and help out a local  Users/JimmyPan

20051208 18:57:53 nbsp Good food, large portions, much better than other local chinese restaurants, not very busy Users/JonathanKaplan

20060216 13:34:00 nbsp There are three things I like about the lunch buffet: (1) for $7.53 I filled a large togo box to the brim, (2) the waitress was supernice, and (3) the beef (huge pieces) and broccoli and the shrimp were yummy.  Thats it!  Some of the string beans were cold, the chefs couldnt seem to decide whether to peel and devein the shrimp or not, and the brightorange chicken (not to be confused with the darkorange chicken) had too much batter, which was made soggy by the sauce.  Even the fortune cookie (Super K brand) tasted awful and the chopsticks had dried orange residue on the eating ends.  Ill go back, but next time Ill order from the menu. Users/JennaChan

20060620 13:58:49 nbsp I had the curry chicken last Friday, and it was incredibly bland, so Im having trouble making much sense out of the glowing reviews above. Curry dishes especially should never be bland. Those people who own and operate the place are surely sweethearts, though. Users/RickEle

20060811 18:47:26 nbsp Had the lunch buffet today, 5.99 per adult.  Had the usual chow mein, white rice, battered chicken in sauce, corn on the cob, various fruits, shrimp, green beans, stir fry chicken with baby corn.  Sorry to say that it wasnt fantastic food and I probably wouldnt get the lunch buffet again.  It was just too salty for me; however, I might try the menu another time their service was great though!  Staff was friendly and accomodating. Users/KatYu

20070131 08:58:14 nbsp Their buffet vs Great Wall  which one is better? Users/FrankYoung

20070313 10:50:31 nbsp i thought that the food (not during the buffet was mixed).  The service was the best in davis (baring upscale places) complimentary cream cheese wontons (pretty good) and lemon water was a nice touch.  Very attentative server, like 3 rice refills and i cant count the water refills, and when we were done and got a to go box re added extra rice in the empty compartments.  As for the food the flavors are quite decent, but one dish was really bad.  The first bite was great but as you ate it, it had a major aftertaste that sort of was a sour tingling, like floride at the dentist, lemony water only compounded this, this dish ended up being almost unettible eventhough the first bite was great.  (maybe it was slightly going bad chilis or something, since it was the dish we ordered spicy).  I have to say it was quite odd.  i think im going to give them a second chance just to see if that was a fluke. Users/MattHh

20070326 20:10:27 nbsp I love shangrila. If you are surprised that the buffet doesnt fulfill your culinary fantasies, try reminding yourself that a: its a buffet and b: its Davis. To say that food is uninteresting here in Davis is an insult to uninteresting food.  Anyways, on to the positives. These folks make great vegetarian dishes, which is a major plus. I love meat and eat it daily but I dont need it 3X a day. Dont get the buffet food unless you like buffet food, plain as that. And as everyone else noted, the people who run the place are extra nice and friendly. Overall, Aminus.  Users/YesItsMe

20070426 12:07:55 nbsp Had their lunch buffet. I dont know if it was typical for them or not but it was the worst Chinese food Ive ever had. The pot stickers, won tons, chow mein, and orange chicken were all cold and Im thinking theyre purposely doing that to keep it from getting dried out. The wontons barely had any cream cheese in them at all and the shell was cold, hard and stale. The orange chicken and sweet and sour pork were mostly batter with very little meat and the pork was stale and tough. The brocoli beef was the only decent dish. Users/MattMuller

20070429 19:10:57 nbsp Sigh. What can you do when people dont read the comments? News flash, folks: BUFFETS SUCK! What do you expect? THEYRE BUFFETS! Shangrila is good for its regular dishes. I wouldnt go for the buffet here or anywhere else on this planet. Im on my way out the door this minute to pick up the General Chicken I just ordered and am certain it will be just what I want it to be: a good main dish, not BUFFET FOOD! LOL
I mean all this in a goodnatured way. Why complain about something that sucks when you know its going to suck in the first place?? Users/YesItsMe

20070802 22:06:29 nbsp With these mixed reviews I tried this place to see how it really is and I have to say it isnt as great as Great Wall. So sad, I was hoping for a good Chinese restaurant (. I too experienced a weird taste with their food like Matt. I got the salt & pepper pork ribs which first was not even the name for the dish...but also after the first bite the rest tasted so bad I had to throw it away...thats how bad it tasted...the ingredients didnt taste very fresh either, price is higher than Great Wall, but service was friendly. Waited for our food and the lady gave us lemon water. Overall I wouldnt come here again. Users/BrandonToo

20071017 21:35:21 nbsp Ive eaten at this place a few times, and I have to say that it is AUTHENTIC chinese food. I know the buffet isnt great, because they dont have a lot of customers so the food doesnt get changed too frequently. BUT, i ordered the Szechuan beef there, and it is AWESOME. Im chinese, I would know. I spoke to the owner/chef, hes FROM Szechuan, China. Completely authentic, great food. Try the Mapo Tofu, Twice Cooked pork. Those are authentic chinese dishes. Not the orange chicken crap most ppl get. If you go there, try the spicy foods. Dont get the usual Americanized fried rice and chow mein stuff... Users/Jeanne

20071219 19:43:44 nbsp My family and I went there not long ago. The food was adequate and the service was attentive, but my dad got food poisoning from his meal and we havent been back since. Users/thelonepiper555
  Are you sure it was food poisoning?  Did it appear several hours to a day later?  Are you sure it wasnt an earlier meal?

It was about 12 hours later, we went there after my dad had a blood test that morning (fasting was necessary), and we are without a doubt that it was food poisoning. Thats what his doctor told him the next day. I wouldnt be surprised if that was just a oneoff situation that was of no fault of their own, but I dont blame my dad for feeling weary of the restaurant. Users/thelonepiper555
  Yeah 12 hours is like textbook timing for food poisoning. One dish i had once at this restaurant had a really bizare flouride aftertaste that was highly suspect, so i believed you, i just post that comment after food poisoning comments on all pages because many people claim food poisoning when they experience symptons imediately (food poisoning symptoms generally manifest several hours after infection). I think that places that do very little business are often more likely to have food poisoning incodents because the ingredients are not as fresh or in the case of a buffet the food is set out for a long period of time. I just got food poisoning from a place in boston (dimsum at 3:00pm on a weekday in an empty restaurant is a bad idea) Users/MattHh

20071227 13:54:47 nbsp My girlfriend and I love their vegetarian dishes.  The eggplant and tofu chicken are especially good.  Dont order the tofu sweet & sour pork, though, which just consists of cubes of plain tofu with sweet & sour sauce on it.  My only complaint is that they are often closed at random times when their scheduled hours says they will be open, but I suppose this is probably because of the fact that they do very little business, and probably just cant afford to stay open all the time.  Its kind of really sad actually, because I bet that if they were downtown they would have a thriving business.  The food is really good, they just happen to have a really bad location. Users/hokusman

20080205 01:06:02 nbsp This place definitely not open until 11pm. Ive tried to phone in an order twice. I called around 10:30 one night and was told that they were closing, then I called another night at 9:45 but got no answer. I drove by around 10pm another night and saw no lights. This is very disappointing. Ive been craving Chinese food and I really wanted to check this place out because these comments made it sound worth trying. Users/HopeBravo

20080212 13:23:52 nbsp Appears to be closed? Has been closed for several days now. No menu in the window. Users/MupDup

20080519 21:25:52 nbsp The phone number is disconnected Users/TiffanyWarrick

20080831 20:19:08 nbsp As of 20080609, the Red Orchid is open. See their page Red Orchid Restaurant for more info. Users/TomHinds

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


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

C $Id: inbody.f,v 1.1 1990/11/30 11:09:17 gdsjaar Exp $
C $Log: inbody.f,v $
C Revision 1.1  1990/11/30 11:09:17  gdsjaar
C Initial revision
C
C
CC* FILE: [.MAIN]INBODY.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE INBODY (MR, N9, IIN, IFOUND, IRPB, ADDOLD, NOROOM)
C***********************************************************************
C
C  SUBROUTINE INBODY = INPUTS A BODY LIST INTO THE DATABASE
C
C***********************************************************************
C
      DIMENSION IIN (IFOUND), IRPB (MR)
C
      LOGICAL NOROOM, ADDOLD
C
      NOROOM = .TRUE.
      IF (.NOT.ADDOLD)N9 = 0
      DO 120 I = 1, IFOUND
         JJ = IIN (I)
         IF (JJ .EQ. 0)GOTO 130
         IF (N9 + 1 .GT. MR)RETURN
C
C  SEE IF THE REGION IS ALREADY IN THE BODY LIST
C
         DO 100 J = 1, N9
            IF (IRPB (J) .EQ. JJ)GOTO 110
  100    CONTINUE
C
         N9 = N9 + 1
         IRPB (N9) = JJ
  110    CONTINUE
  120 CONTINUE
  130 CONTINUE
C
      NOROOM = .FALSE.
      RETURN
C
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT