text:
      SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
     *rkqs)
      INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
      REAL eps,h1,hmin,x1,x2,ystart(nvar),TINY
      EXTERNAL derivs,rkqs
      PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=200,TINY=1.e-30)
C	Runge-Kutta driver with adaptive stepsize control.  Integrate the
C	starting values ystart(1:nvar)from x1 to x2 with accuracy eps,
C	storing intermediate results in the common block /path/
C	h1 should be set as a guessed first stepsize.  Hmin as the minimum
C	allowed stesize (can be zero)  On output nok and nbad are the number
C	of good and bad (but retried and fixed) steps taken, and ystart is
C	replaced by values at the end of the integration interval.  Derivs
C	is the user-supplied subrouting for calculating the right hand side
C	derivative, while RKQS is the name of the stepper routine to be
C	used. /path/ contains its own information about how often an
C	intermediate value is to be stored.
      INTEGER i,kmax,kount,nstp
      REAL dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX),
     *yp(NMAX,KMAXX),yscal(NMAX)
      COMMON /path/ kmax,kount,dxsav,xp,yp
      x=x1
      h=sign(h1,x2-x1)
      nok=0
      nbad=0
      kount=0
      do 11 i=1,nvar
        y(i)=ystart(i)
11    continue
      if (kmax.gt.0) xsav=x-2.*dxsav
      do 16 nstp=1,MAXSTP
        call derivs(x,y,dydx)
        do 12 i=1,nvar
          yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
12      continue
        if(kmax.gt.0)then
          if(abs(x-xsav).gt.abs(dxsav)) then
            if(kount.lt.kmax-1)then
              kount=kount+1
              xp(kount)=x
              do 13 i=1,nvar
                yp(i,kount)=y(i)
13            continue
              xsav=x
            endif
          endif
        endif
        if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x
        call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
        if(hdid.eq.h)then
          nok=nok+1
        else
          nbad=nbad+1
        endif
        if((x-x2)*(x2-x1).ge.0.)then
          do 14 i=1,nvar
            ystart(i)=y(i)
14        continue
          if(kmax.ne.0)then
            kount=kount+1
            xp(kount)=x
            do 15 i=1,nvar
              yp(i,kount)=y(i)
15          continue
          endif
          return
        endif
        if(abs(hnext).lt.hmin) THEN
	  PRINT*,'PAUSE AT X=',X,' HNEXT, HMIN',HNEXT,HMIN
	  pause 'stepsize smaller than minimum in odeint'
	ENDIF
        h=hnext
16    continue
      pause 'too many steps in odeint'
      return
      END


text:
      program wave2d
      use module_vernieuw
!*****************************************!
!* 2d shallow-water model                *!
!*                                       *!
!* including:                            *!
!* - horizontal pressure-gradient force  *!
!* - shapiro filter                      *!
!* - flooding algorithm                  *!
!*                                       *!
!* author: j. kaempf, 2008               *!
!*****************************************!
use sub
! local parameters
integer(4), parameter :: nx=500
integer(4), parameter :: ny=500
      real :: dt
      real :: dx
      real :: dy
      real :: eps
      real, dimension(0:ny+1,0:nx+1) :: eta
      real, dimension(0:ny+1,0:nx+1) :: etan
      real :: g
      real, dimension(0:ny+1,0:nx+1) :: h
      real :: hmin
      real, dimension(0:ny+1,0:nx+1) :: hzero
      integer :: j
      integer :: k
      real, dimension(0:ny+1,0:nx+1) :: u
      real, dimension(0:ny+1,0:nx+1) :: un
      real, dimension(0:ny+1,0:nx+1) :: v
      real, dimension(0:ny+1,0:nx+1) :: vn
      integer, dimension(0:ny+1,0:nx+1) :: wet
real :: hmax
real :: time
real :: dtmax
real :: c
real :: lambda
integer :: n
integer :: ntot
integer :: nout
integer :: dummy
! set local parameters
! runtime parameters
ntot = 10000
! output parameter
!nout = 5
! set epsilon for Shapiro filter
eps = 0.05
call init(hmin,dx,dy,dt,g,j,k,hzero,eta,etan,h,wet,u,un,v,vn)

!**********
!debug
open(90,file ='debug.dat',form='formatted')
! output of initial eta distribution
open(10,file ='eta0.dat',form='formatted')
  do j = 0,ny+1
    write(10,'(101F12.6)')(eta(j,k),k=0,nx+1)
  end do
close(10)
! output of initial layer thickness distribution
open(10,file ='h0.dat',form='formatted')
  do j = 0,ny+1
    write(10,'(101F12.6)')(hzero(j,k),k=0,nx+1)
  end do
close(10)
! determine maximum water depth
hmax = 0.
do j = 1,ny
do k = 1,nx
  hmax = max(hmax,h(j,k))
end do
end do
! maximum phase speed
dummy = 0
c = sqrt(2*g*hmax)
!
write(6,*)"c = ",c
! determine stability parameter
lambda = dt*sqrt(g*hmax)/min(dx,dy)
write(6,*)"lambda = ",lambda
if(lambda > 1)then
  write(6,*) "This will not work. Do you know why?"
  stop
end if
! open files for output
open(10,file ='eta.dat',form='formatted')
open(20,file ='h.dat',form='formatted')
open(30,file ='u.dat',form='formatted')
open(40,file ='v.dat',form='formatted')
do j = 26,26
do k = 26,26
eta(j,k) = 1.0
end do
end do
!---------------------------
! simulation loop
!---------------------------
do n = 1,ntot
time = real(n)*dt
! call predictor
  call dyn(j,k,dx,g,eta,dt,dy,un,u,wet,v,vn,h,etan)

! updating including Shapiro filter
  call shapiro(j,k,wet,etan,eps,eta)

      call vernieuw(dt,dx,dy,eps,eta,etan,g,h,hmin,hzero,j,k,u,un,v,vn,wet)

! data output
! Only write once at the end
!if(mod(n,nout)==0)then
!  DO j = 1,ny
!    WRITE(10,'(101F12.6)')(eta(j,k),k=1,nx)
!    WRITE(20,'(101F12.6)')(h(j,k),k=1,nx)
!    WRITE(30,'(101F12.6)')(u(j,k),k=1,nx)
!    WRITE(40,'(101F12.6)')(v(j,k),k=1,nx)
!  END DO
!  WRITE(6,*)"Data output at time = ",time/60.0," min"
!endif
end do
!write results only once at the end....
do j = 0,ny+1
  write(10,'(101F12.6)')(eta(j,k),k=0,nx+1)
  write(20,'(101F12.6)')(h(j,k)  ,k=0,nx+1)
  write(30,'(101F12.6)')(u(j,k)  ,k=0,nx+1)
  write(40,'(101F12.6)')(v(j,k)  ,k=0,nx+1)
end do
write(6,*)"Data output at time = ",time/60.0," min"
end program wave2d

max_stars_repo_name: wimvanderbauwhede/AutoParallel-Fortran

max_issues_repo_name: wimvanderbauwhede/AutoParallel-Fortran

max_forks_repo_name: wimvanderbauwhede/AutoParallel-Fortran


text:
      subroutine runY_0000lli(k,l,i1,Xtwiddle,Gtwiddle,Shat7zz,N0)
      implicit none
C---  Expression for Eq. 5.60b
C---  Calculates D0000lli, requires D0000lll
C---  Small terms of order Xtwiddle(0,k)*D00iiii,Xtwiddle(0,0)*D00iiiii
C---  Denominator Gtwiddle(k,l)
      include 'pvDnames.f'
      include 'pvDv.f'
      include 'Darraydef.f'
      include 'Darrays.f'
      integer ep,N0,k,l,i1,np
      parameter(np=3)
      double precision Xtwiddle(0:np,0:np),Gtwiddle(np,np)
      double complex Shat7zz(np,z4max,-2:0)

      if ((i1 .eq. l) .or. (i1 .eq. 0)) then
      return
      endif

      do ep=-2,0

      Dv(dzzzziii(z3(l,l,i1))+N0,ep)=
     . (-2d0*Gtwiddle(k,i1)*Dv(dzzzziii(z3(l,l,l))+N0,ep)
     . +Gtwiddle(k,1)*Shat7zz(1,z4(l,l,l,i1),ep)
     . +Gtwiddle(k,2)*Shat7zz(2,z4(l,l,l,i1),ep)
     . +Gtwiddle(k,3)*Shat7zz(3,z4(l,l,l,i1),ep)
     . +Xtwiddle(k,0)*Dv(dzziiii(z4(l,l,l,i1))+N0,ep)
     . -Xtwiddle(0,0)*Dv(dzziiiii(z5(k,l,l,l,i1))+N0,ep))
     . /(6d0*Gtwiddle(k,l))

      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C--------------------------------------------------------------------------
        subroutine ftukys(ounit,keywrd,strval,comm,status)

C       update a character string value header record

C       ounit   i  fortran output unit number
C       keywrd  c  keyword name    ( 8 characters, cols.  1- 8)
C       strval  c  keyword value
C       comm    c  keyword comment (47 characters, cols. 34-80)
C       OUTPUT PARAMETERS:
C       status  i  output error status (0 = ok)
C
C       written by Wm Pence, HEASARC/GSFC, Oct 1994

        character*(*) keywrd,strval,comm
        integer ounit,status,tstat

        if (status .gt. 0)return

        tstat=status
C       try modifying the keyword, if it exists
        call ftmkys(ounit,keywrd,strval,comm,status)

        if (status .eq. 202)then
C               keyword doesn't exist, so create it
                status=tstat
C               note that this supports the HEASARC long-string conventions
                call ftpkls(ounit,keywrd,strval,comm,status)
        end if
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!
! EmisCoeff_NC2BIN
!
! Program to convert netCDF format EmisCoeff files to the CRTM Binary
! format.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 19-Jul-2005
!                       paul.vandelst@ssec.wisc.edu
!

PROGRAM EmisCoeff_NC2BIN

  ! ------------------
  ! Environment set up
  ! ------------------
  ! Module usage
  USE File_Utility,        ONLY: File_Exists
  USE Message_Handler,     ONLY: SUCCESS, FAILURE, WARNING, INFORMATION, &
                                 Program_Message, Display_Message
  USE EmisCoeff_Define,    ONLY: EmisCoeff_type, &
                                 Destroy_EmisCoeff, Equal_EmisCoeff
  USE EmisCoeff_Binary_IO, ONLY: Read_EmisCoeff_Binary, Write_EmisCoeff_Binary
  USE EmisCoeff_netCDF_IO, ONLY: Read_EmisCoeff_netCDF, Write_EmisCoeff_netCDF
  USE SignalFile_Utility
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'EmisCoeff_NC2BIN'
  CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = &

  ! ---------
  ! Variables
  ! ---------
  INTEGER :: Error_Status
  CHARACTER(256) :: NC_Filename
  CHARACTER(256) :: BIN_Filename
  TYPE(EmisCoeff_type) :: EmisCoeff
  TYPE(EmisCoeff_type) :: EmisCoeff_Test

  ! Output program header
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to convert netCDF format EmisCoeff files to '//&
                        'their CRTM Binary format.', &
                        '$Revision$' )

  ! Get the input and output filenames
  WRITE( *, FMT     = '( /5x, "Enter the INPUT netCDF EmisCoeff file: " )', &
            ADVANCE = 'NO' )
  READ( *, '( a )' ) NC_Filename
  NC_Filename = ADJUSTL( NC_Filename )
  IF ( .NOT. File_Exists( TRIM( NC_Filename ) ) ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'File '//TRIM( NC_Filename )//' not found.', &
                          FAILURE )
    STOP
  END IF

  WRITE( *, FMT     = '( /5x, "Enter the OUTPUT Binary EmisCoeff file: " )', &
            ADVANCE = 'NO' )
  READ( *, '( a )' ) BIN_Filename
  BIN_Filename = ADJUSTL( BIN_FileNAME )

  ! Check that the netCDF file isn't accidentally overwritten
  IF ( TRIM( NC_Filename ) == TRIM( BIN_Filename ) ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Output filename is the same as the input filename!', &
                          FAILURE )
    STOP
  END IF

  ! Read the input netCDF file
  WRITE( *, '( /5x, "Reading netCDF EmisCoeff data ..." )' )
  Error_Status = Read_EmisCoeff_netCDF( NC_Filename, EmisCoeff )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error reading netCDF EmisCoeff file '//&
                          TRIM( NC_Filename ), &
                          Error_Status )
    STOP
  END IF

  ! Write the binary file
  WRITE( *, '( /5x, "Writing Binary EmisCoeff data ..." )' )
  Error_Status = Write_EmisCoeff_Binary( BIN_Filename, EmisCoeff )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error writing Binary EmisCoeff file '//&
                          TRIM( BIN_Filename ), &
                          Error_Status )
    STOP
  END IF

  ! Test read the binary data file
  WRITE( *, '( /5x, "Test reading the Binary EmisCoeff data file ..." )' )
  Error_Status = Read_EmisCoeff_Binary( BIN_Filename, EmisCoeff_Test )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error reading Binary EmisCoeff file '//&
                          TRIM( BIN_Filename ), &
                          Error_Status )
    STOP
  END IF

  ! Compare the two structures
  WRITE( *, '( /5x, "Comparing the netCDF and Binary EmisCoeff structures ..." )' )
  Error_Status = Equal_EmisCoeff( EmisCoeff_Test, EmisCoeff )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Differences found in netCDF and Binary '//&
                          'file EmisCoeff structure comparison.', &
                          Error_Status )
  ELSE
    CALL Display_Message( PROGRAM_NAME, &
                          'netCDF and Binary file EmisCoeff structures are equal.', &
                          INFORMATION )
  END IF

  ! Create a signal file indicating success
  Error_Status = Create_SignalFile( BIN_Filename )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error creating signal file for '//TRIM(BIN_Filename), &
                          FAILURE )
    STOP
  END IF

  ! Destroy the structures
  Error_Status = Destroy_EmisCoeff( EmisCoeff )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying EmisCoeff structure.', &
                          WARNING )
  END IF

  Error_Status = Destroy_EmisCoeff( EmisCoeff_Test )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying EmisCoeff_Test structure.', &
                          WARNING )
  END IF

END PROGRAM EmisCoeff_NC2BIN

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


! { dg-do run }
! { dg-options "-O -fdump-tree-original" }
program main
  implicit none
  character(len=4) :: c
  integer :: n
  integer :: i
  common /foo/ i

  n = 0
  i = 0
  c = 'abcd'
  if ('a ' // c == 'a' // c) call abort
  if ('a' // c == 'a ' // c) call abort
end program main

! { dg-final { scan-tree-dump-times "gfortran_concat_string" 4 "original" } }
! { dg-final { scan-tree-dump-times "gfortran_compare_string" 2 "original" } }


text:
!****** BASIS_FUNCS.F90 ******************************************************************
module basis_funcs

!===========================================================================
! Parameters
!------------------------------------------------------------
use params

integer, parameter :: nvtk  = 1    ! (was 2)
integer, parameter :: nvtk2 = nvtk*nvtk
integer, parameter :: nvtk3 = nvtk*nvtk*nvtk
!===========================================================================

!===========================================================================
! Definitions
!------------------------------------------------------------
real, dimension(nvtk3,nbastot) :: bfvtk, bfvtk_dx, bfvtk_dy, bfvtk_dz

! TODO: only in init, innerintegral, source_calc
real, dimension(nbastot) :: cbasis

! TODO: only in init + innerintegral
integer ibitri  ! set by set_cbasis_3D using chosen value of iquad

! TODO: only in set_weights_3D
real wgt1d(5)    ! wgt1d: quadrature weights for 1-D integration

! TODO: only in init, set_weights_3D, glflux
real wgt2d(30)   ! wgt2d: quadrature weights for 2-D integration

! TODO: only in init, set_weights_3D, innerintegral
real wgt3d(100)  ! wgt3d: quadrature weights for 3-D integration

! TODO: only in limiter, set_face_vals_3D
real, dimension(nslim,nbastot) :: bf_faces

! TODO: only in:
!   * initialize.f90 (setup)
!   * set_face_vals_3D, set_internal_vals_3D
!   * source_calc, innerintegral
real, dimension(npg,nbastot) :: bfvals_int

! TODO: only in set_internal_vals_3D, set_face_vals_3D
real xquad(20)

real, dimension(npg,nbastot) :: bval_int_wgt  ! used in source_calc
real, dimension(nface,2,nbastot) :: wgtbf_xmp, wgtbf_ymp, wgtbf_zmp  ! used in glflux

! TODO: only in init, flux_cal, prepare_exchange, set_face_vals_3D
real, dimension(nface,nbastot) :: bfvals_zp, bfvals_zm
real, dimension(nface,nbastot) :: bfvals_yp, bfvals_ym
real, dimension(nface,nbastot) :: bfvals_xp, bfvals_xm

! TODO: original types are I4P from LIB_VTK_IO.f90
! integer(I4P), parameter :: nnx=nx*nvtk, nny=ny*nvtk, nnz=nz*nvtk
! integer, parameter :: nnx=nx*nvtk, nny=ny*nvtk, nnz=nz*nvtk

real :: xgrid(20) ! used in set_vtk_vals_3D

! These are used in test_basis_3D
real, dimension(nbastot,nbastot) :: cell_int,xp_int,xm_int,yp_int,ym_int,zp_int,zm_int
real, dimension(nbastot,nbastot) :: cell_int0,xp_int0,xm_int0,yp_int0,ym_int0,zp_int0,zm_int0
real, dimension(nbastot) :: cbas_xp,cbas_xm,cbas_yp,cbas_ym,cbas_zp,cbas_zm
integer, dimension(nbastot) :: ibas_x, ibas_y, ibas_z
!===========================================================================

contains

    subroutine set_bfvals_3D
        ! Defines local basis function values and weights for 1, 2, or 3-point Gaussian quadrature.
        ! Basis functions are evaluated in cell interior and on cell faces.
        implicit none

        call set_vtk_vals_3D()    ! Define basis function values on a 3D grid INTERNAL to the cell.
        call set_internal_vals_3D()     ! Define basis function values at quadrature points INTERNAL to cell.
        call set_face_vals_3D()   ! Define local basis function values at quadrature points on a cell face.
        call set_weights_3D()     ! Define weights for integral approximation using Gaussian quadrature.

        if (test_basis .and. iam == print_mpi) then
    	   print *,'Testing basis...'
    	   call test_basis_3D()
    	   print *,'Done testing basis...'
    	end if

    end subroutine set_bfvals_3D

!----------------------------------------------------------
    subroutine set_vtk_vals_3D
        ! Define basis function values on a rectangular 3D grid INTERNAL to the cell.
        ! For use in output_vtk().

        implicit none
        integer ixyz,i,igrid

        dxvtk = 1./(nvtk*dxi)
        dyvtk = 1./(nvtk*dyi)
        dzvtk = 1./(nvtk*dzi)

        if (nvtk .eq. 1) xgrid(1) = 0.

        if (nvtk .eq. 2) then
            xgrid(1) = -0.5
            xgrid(2) = 0.5
        end if
        if (nvtk .eq. 3) then
            xgrid(1) = -2./3.
            xgrid(2) = 0.
            xgrid(3) = 2./3.
        end if
        if (nvtk .eq. 4) then
            xgrid(1) = -3./4.
            xgrid(2) = -1./4.
            xgrid(3) = 1./4.
            xgrid(4) = 3./4.
        end if

        bfvtk(1,1) = 1.        ! basis function = 1
        bfvtk(1,kx) = 0.        ! basis function = x
        bfvtk(1,ky) = 0.        ! basis function = y
        bfvtk(1,kz) = 0.        ! basis function = z

        bfvtk(1:nvtk3,1) = 1.        ! basis function = 1
        do i=1,nvtk
            bfvtk((i-1)*nvtk+1:i*nvtk,ky) = xgrid(i)        ! basis function = y
            bfvtk((i-1)*nvtk+1:i*nvtk,kz) = xgrid(1:nvtk)        ! basis function = z
        end do
        do i=1,nvtk
            bfvtk((i-1)*nvtk2+1:i*nvtk2,kx) = xgrid(i)        ! basis function = x
            bfvtk((i-1)*nvtk2+1:i*nvtk2,ky) = bfvtk(1:nvtk2,ky)        ! basis function = y
            bfvtk((i-1)*nvtk2+1:i*nvtk2,kz) = bfvtk(1:nvtk2,kz)        ! basis function = z
        end do

        do i=0,2
            bfvtk(1:nvtk3,kxx+i) = 1.5*bfvtk(1:nvtk3,kx+i)**2 - 0.5    ! basis function = P_2(s)
            ! bfvtk(1:nvtk3,kxxx+i) = 2.5*bfvtk(1:nvtk3,kx+i)**3 - 1.5*bfvtk(1:nvtk3,kx+i)   ! basis function = P3(s)
        end do

        bfvtk(1:nvtk3,kyz) = bfvtk(1:nvtk3,ky)*bfvtk(1:nvtk3,kz)        ! basis function = yz
        bfvtk(1:nvtk3,kzx) = bfvtk(1:nvtk3,kz)*bfvtk(1:nvtk3,kx)        ! basis function = zx
        bfvtk(1:nvtk3,kxy) = bfvtk(1:nvtk3,kx)*bfvtk(1:nvtk3,ky)        ! basis function = xy
        bfvtk(1:nvtk3,kxyz) = bfvtk(1:nvtk3,kx)*bfvtk(1:nvtk3,ky)*bfvtk(1:nvtk3,kz)     ! basis function = xyz

        bfvtk(1:nvtk3,kyyz) = bfvtk(1:nvtk3,kyy)*bfvtk(1:nvtk3,kz)     ! basis function = P2(y)z
        bfvtk(1:nvtk3,kyzz) = bfvtk(1:nvtk3,ky)*bfvtk(1:nvtk3,kzz)     ! basis function = P2(z)y
        bfvtk(1:nvtk3,kzzx) = bfvtk(1:nvtk3,kzz)*bfvtk(1:nvtk3,kx)     ! basis function = P2(z)x
        bfvtk(1:nvtk3,kzxx) = bfvtk(1:nvtk3,kz)*bfvtk(1:nvtk3,kxx)     ! basis function = P2(x)z
        bfvtk(1:nvtk3,kxxy) = bfvtk(1:nvtk3,kxx)*bfvtk(1:nvtk3,ky)     ! basis function = P2(x)y
        bfvtk(1:nvtk3,kxyy) = bfvtk(1:nvtk3,kx)*bfvtk(1:nvtk3,kyy)     ! basis function = P2(y)x
        bfvtk(1:nvtk3,kyyzz) = bfvtk(1:nvtk3,kyy)*bfvtk(1:nvtk3,kzz)     ! basis function = P_2(y)P_2(z)
        bfvtk(1:nvtk3,kzzxx) = bfvtk(1:nvtk3,kzz)*bfvtk(1:nvtk3,kxx)     ! basis function = P_2(z)P_2(x)
        bfvtk(1:nvtk3,kxxyy) = bfvtk(1:nvtk3,kxx)*bfvtk(1:nvtk3,kyy)     ! basis function = P_2(x)P_2(y)
        bfvtk(1:nvtk3,kyzxx) = bfvtk(1:nvtk3,kyz)*bfvtk(1:nvtk3,kxx)     ! basis function = yz P_2(x)
        bfvtk(1:nvtk3,kzxyy) = bfvtk(1:nvtk3,kzx)*bfvtk(1:nvtk3,kyy)     ! basis function = zx P_2(y)
        bfvtk(1:nvtk3,kxyzz) = bfvtk(1:nvtk3,kxy)*bfvtk(1:nvtk3,kzz)     ! basis function = xy P_2(z)
        bfvtk(1:nvtk3,kxyyzz) = bfvtk(1:nvtk3,kx)*bfvtk(1:nvtk3,kyy)*bfvtk(1:nvtk3,kzz)   ! basis function = x P_2(y)P_2(z)
        bfvtk(1:nvtk3,kyzzxx) = bfvtk(1:nvtk3,ky)*bfvtk(1:nvtk3,kzz)*bfvtk(1:nvtk3,kxx)   ! basis function = y P_2(z)P_2(x)
        bfvtk(1:nvtk3,kzxxyy) = bfvtk(1:nvtk3,kz)*bfvtk(1:nvtk3,kxx)*bfvtk(1:nvtk3,kyy)   ! basis function = z P_2(x)P_2(y)
        bfvtk(1:nvtk3,kxxyyzz) = bfvtk(1:nvtk3,kx)*bfvtk(1:nvtk3,kyy)*bfvtk(1:nvtk3,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

        do igrid=1,nvtk3
            bfvtk_dx(igrid,1) = 0.
            bfvtk_dx(igrid,kx) = 1.
            bfvtk_dx(igrid,ky:kz) = 0.
            bfvtk_dx(igrid,kyz) = 0.
            bfvtk_dx(igrid,kzx) = bfvtk(igrid,kz)
            bfvtk_dx(igrid,kxy) = bfvtk(igrid,ky)
            bfvtk_dx(igrid,kxyz) = bfvtk(igrid,kyz)

            bfvtk_dx(igrid,kxx) = 3.*bfvtk(igrid,kx)
            bfvtk_dx(igrid,kyy:kzz) = 0.
            bfvtk_dx(igrid,kyyz) = 0.
            bfvtk_dx(igrid,kyzz) = 0.
            bfvtk_dx(igrid,kzzx) = bfvtk(igrid,kzz)
            bfvtk_dx(igrid,kzxx) = bfvtk(igrid,kz)*3.*bfvtk(igrid,kx)
            bfvtk_dx(igrid,kxxy) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,ky)
            bfvtk_dx(igrid,kxyy) = bfvtk(igrid,kyy)
            bfvtk_dx(igrid,kyyzz) = 0.
            bfvtk_dx(igrid,kzzxx) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kzz)
            bfvtk_dx(igrid,kxxyy) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kyy)
            bfvtk_dx(igrid,kyzxx) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kyz)
            bfvtk_dx(igrid,kzxyy) = bfvtk(igrid,kz)*bfvtk(igrid,kyy)
            bfvtk_dx(igrid,kxyzz) = bfvtk(igrid,ky)*bfvtk(igrid,kzz)
            bfvtk_dx(igrid,kxyyzz) = bfvtk(igrid,kyyzz)
            bfvtk_dx(igrid,kyzzxx) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kyzz)
            bfvtk_dx(igrid,kzxxyy) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kyyz)
            bfvtk_dx(igrid,kxxyyzz) = 3.*bfvtk(igrid,kx)*bfvtk(igrid,kyy)*bfvtk(igrid,kzz)

            bfvtk_dy(igrid,1) = 0.
            bfvtk_dy(igrid,ky) = 1.
            bfvtk_dy(igrid,kx) = 0.
            bfvtk_dy(igrid,kz) = 0.
            bfvtk_dy(igrid,kyz) = bfvtk(igrid,kz)
            bfvtk_dy(igrid,kzx) = 0.
            bfvtk_dy(igrid,kxy) = bfvtk(igrid,kx)
            bfvtk_dy(igrid,kxyz) = bfvtk(igrid,kzx)

            bfvtk_dy(igrid,kyy) = 3.*bfvtk(igrid,ky)
            bfvtk_dy(igrid,kxx) = 0.
            bfvtk_dy(igrid,kzz) = 0.
            bfvtk_dy(igrid,kyzz) = bfvtk(igrid,kzz)
            bfvtk_dy(igrid,kyyz) = bfvtk(igrid,kz)*3.*bfvtk(igrid,ky)
            bfvtk_dy(igrid,kzzx) = 0.
            bfvtk_dy(igrid,kzxx) = 0.
            bfvtk_dy(igrid,kxxy) = bfvtk(igrid,kxx)
            bfvtk_dy(igrid,kxyy) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kx)
            bfvtk_dy(igrid,kyyzz) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kzz)
            bfvtk_dy(igrid,kzzxx) = 0.
            bfvtk_dy(igrid,kxxyy) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kxx)
            bfvtk_dy(igrid,kyzxx) = bfvtk(igrid,kz)*bfvtk(igrid,kxx)
            bfvtk_dy(igrid,kzxyy) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kzx)
            bfvtk_dy(igrid,kxyzz) = bfvtk(igrid,kx)*bfvtk(igrid,kzz)
            bfvtk_dy(igrid,kxyyzz) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kzzx)
            bfvtk_dy(igrid,kyzzxx) = bfvtk(igrid,kzzxx)
            bfvtk_dy(igrid,kzxxyy) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kzxx)
            bfvtk_dy(igrid,kxxyyzz) = 3.*bfvtk(igrid,ky)*bfvtk(igrid,kzz)*bfvtk(igrid,kxx)

            bfvtk_dz(igrid,1) = 0.
            bfvtk_dz(igrid,kz) = 1.
            bfvtk_dz(igrid,kx) = 0.
            bfvtk_dz(igrid,ky) = 0.
            bfvtk_dz(igrid,kyz) = bfvtk(igrid,ky)
            bfvtk_dz(igrid,kzx) = bfvtk(igrid,kx)
            bfvtk_dz(igrid,kxy) = 0.
            bfvtk_dz(igrid,kxyz) = bfvtk(igrid,kxy)

            bfvtk_dz(igrid,kzz) = 3.*bfvtk(igrid,kz)
            bfvtk_dz(igrid,kxx) = 0.
            bfvtk_dz(igrid,kyy) = 0.

            bfvtk_dz(igrid,kyzz) = bfvtk(igrid,ky)*3.*bfvtk(igrid,kz)
            bfvtk_dz(igrid,kyyz) = bfvtk(igrid,kyy)
            bfvtk_dz(igrid,kzzx) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kx)
            bfvtk_dz(igrid,kzxx) = bfvtk(igrid,kxx)
            bfvtk_dz(igrid,kxxy) = 0.
            bfvtk_dz(igrid,kxyy) = 0.
            bfvtk_dz(igrid,kyyzz) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kyy)
            bfvtk_dz(igrid,kzzxx) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kxx)
            bfvtk_dz(igrid,kxxyy) = 0.
            bfvtk_dz(igrid,kyzxx) = bfvtk(igrid,ky)*bfvtk(igrid,kxx)
            bfvtk_dz(igrid,kzxyy) = bfvtk(igrid,kx)*bfvtk(igrid,kyy)
            bfvtk_dz(igrid,kxyzz) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kxy)
            bfvtk_dz(igrid,kxyyzz) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kxyy)
            bfvtk_dz(igrid,kyzzxx) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kxxy)
            bfvtk_dz(igrid,kzxxyy) = bfvtk(igrid,kxxyy)
            bfvtk_dz(igrid,kxxyyzz) = 3.*bfvtk(igrid,kz)*bfvtk(igrid,kxx)*bfvtk(igrid,kyy)

            ! bfvtk_dx(igrid,kxxx) = 7.5*bfvtk(igrid,kx)**2 - 1.5
            ! bfvtk_dx(igrid,kyyy:kzzz) = 0.
            ! bfvtk_dy(igrid,kyyy) = 7.5*bfvtk(igrid,ky)**2 - 1.5
            ! bfvtk_dy(igrid,kxxx) = 0.
            ! bfvtk_dy(igrid,kzzz) = 0.
            ! bfvtk_dz(igrid,kzzz) = 7.5*bfvtk(igrid,kz)**2 - 1.5
            ! bfvtk_dz(igrid,kxxx) = 0.
            ! bfvtk_dz(igrid,kyyy) = 0.
        end do

    end subroutine set_vtk_vals_3D

!----------------------------------------------------------

!----------------------------------------------------------
    subroutine set_internal_vals_3D
        ! Define basis function values at quadrature points INTERNAL to cell.

        implicit none
        integer ixyz, i
        real c15d5,c1dsq3,xq4p,xq4m

        c15d5 = sqrt(15.)/5.
        c1dsq3 = 1./sqrt(3.)

        xq4p = sqrt(525. + 70.*sqrt(30.))/35.
        xq4m = sqrt(525. - 70.*sqrt(30.))/35.

        if (iquad .eq. 2) then
            xquad(1) = -c1dsq3
            xquad(2) = c1dsq3
        end if
        if (iquad .eq. 3) then
            xquad(1) = -c15d5
            xquad(2) = 0.
            xquad(3) = c15d5
        end if
        if (iquad .eq. 4) then
            xquad(1) = -xq4p
            xquad(2) = -xq4m
            xquad(3) = xq4m
            xquad(4) = xq4p
        end if

        if (iquad .eq. 1) then            ! 2-point Gaussian quadrature
            bfvals_int(1,1) = 1.        ! basis function = 1
            bfvals_int(1,kx) = 0.        ! basis function = x
            bfvals_int(1,ky) = 0.        ! basis function = y
            bfvals_int(1,kz) = 0.        ! basis function = z
        end if

        if (iquad .gt. 1) then
            bfvals_int(1:npg,1) = 1.        ! basis function = 1
            do i=1,nedge
                bfvals_int((i-1)*nedge+1:i*nedge,ky) = xquad(i)        ! basis function = y
                bfvals_int((i-1)*nedge+1:i*nedge,kz) = xquad(1:nedge)        ! basis function = z
            end do
            do i=1,nedge
                bfvals_int((i-1)*nface+1:i*nface,kx) = xquad(i)        ! basis function = x
                bfvals_int((i-1)*nface+1:i*nface,ky) = bfvals_int(1:nface,ky)        ! basis function = y
                bfvals_int((i-1)*nface+1:i*nface,kz) = bfvals_int(1:nface,kz)        ! basis function = z
            end do
        end if

        do i=0,2
            bfvals_int(1:npg,kxx+i) = 1.5*bfvals_int(1:npg,kx+i)**2 - 0.5    ! basis function = P_2(s)
            ! bfvals_int(1:npg,kxxx+i) = 2.5*bfvals_int(1:npg,kx+i)**3 - 1.5*bfvals_int(1:npg,kx+i)   ! basis function = P_3(s)
        end do

        bfvals_int(1:npg,kyz) = bfvals_int(1:npg,ky)*bfvals_int(1:npg,kz)        ! basis function = yz
        bfvals_int(1:npg,kzx) = bfvals_int(1:npg,kz)*bfvals_int(1:npg,kx)        ! basis function = zx
        bfvals_int(1:npg,kxy) = bfvals_int(1:npg,kx)*bfvals_int(1:npg,ky)        ! basis function = xy
        bfvals_int(1:npg,kxyz) = bfvals_int(1:npg,kx)*bfvals_int(1:npg,ky)*bfvals_int(1:npg,kz)     ! basis function = xyz

        bfvals_int(1:npg,kyyz) = bfvals_int(1:npg,kyy)*bfvals_int(1:npg,kz)     ! basis function = P_2(y)z
        bfvals_int(1:npg,kyzz) = bfvals_int(1:npg,ky)*bfvals_int(1:npg,kzz)     ! basis function = P_2(z)y
        bfvals_int(1:npg,kzzx) = bfvals_int(1:npg,kzz)*bfvals_int(1:npg,kx)     ! basis function = P_2(z)x
        bfvals_int(1:npg,kzxx) = bfvals_int(1:npg,kz)*bfvals_int(1:npg,kxx)     ! basis function = P_2(x)z
        bfvals_int(1:npg,kxxy) = bfvals_int(1:npg,kxx)*bfvals_int(1:npg,ky)     ! basis function = P_2(x)y
        bfvals_int(1:npg,kxyy) = bfvals_int(1:npg,kx)*bfvals_int(1:npg,kyy)     ! basis function = P_2(y)x
        bfvals_int(1:npg,kyyzz) = bfvals_int(1:npg,kyy)*bfvals_int(1:npg,kzz)     ! basis function = P_2(y)P_2(z)
        bfvals_int(1:npg,kzzxx) = bfvals_int(1:npg,kzz)*bfvals_int(1:npg,kxx)     ! basis function = P_2(z)P_2(x)
        bfvals_int(1:npg,kxxyy) = bfvals_int(1:npg,kxx)*bfvals_int(1:npg,kyy)     ! basis function = P_2(x)P_2(y)
        bfvals_int(1:npg,kyzxx) = bfvals_int(1:npg,kyz)*bfvals_int(1:npg,kxx)     ! basis function = yz P_2(x)
        bfvals_int(1:npg,kzxyy) = bfvals_int(1:npg,kzx)*bfvals_int(1:npg,kyy)     ! basis function = zx P_2(y)
        bfvals_int(1:npg,kxyzz) = bfvals_int(1:npg,kxy)*bfvals_int(1:npg,kzz)     ! basis function = xy P_2(z)
        bfvals_int(1:npg,kxyyzz) = bfvals_int(1:npg,kx)*bfvals_int(1:npg,kyy)*bfvals_int(1:npg,kzz)   ! basis function = x P_2(y)P_2(z)
        bfvals_int(1:npg,kyzzxx) = bfvals_int(1:npg,ky)*bfvals_int(1:npg,kzz)*bfvals_int(1:npg,kxx)   ! basis function = y P_2(z)P_2(x)
        bfvals_int(1:npg,kzxxyy) = bfvals_int(1:npg,kz)*bfvals_int(1:npg,kxx)*bfvals_int(1:npg,kyy)   ! basis function = z P_2(x)P_2(y)
        bfvals_int(1:npg,kxxyyzz) = bfvals_int(1:npg,kx)*bfvals_int(1:npg,kyy)*bfvals_int(1:npg,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

    end subroutine set_internal_vals_3D

!----------------------------------------------------------
    subroutine set_face_vals_3D
        ! Define local basis function values at quadrature points on a cell face.
        ! Used in flux_cal() and prepare_exchange() for interpolating from cell center onto cell face.
        ! Also used in glflux().

        implicit none
        integer ixyz,i
        real c15d5,c1dsq3,xq4p,xq4m

        c15d5 = sqrt(15.)/5.
        c1dsq3 = 1./sqrt(3.)

        xq4p = sqrt(525. + 70.*sqrt(30.))/35.
        xq4m = sqrt(525. - 70.*sqrt(30.))/35.

        if (iquad .eq. 2) then
            xquad(1) = -c1dsq3
            xquad(2) = c1dsq3
        end if
        if (iquad .eq. 3) then
            xquad(1) = -c15d5
            xquad(2) = 0.
            xquad(3) = c15d5
        end if
        if (iquad .eq. 4) then
            xquad(1) = -xq4p
            xquad(2) = -xq4m
            xquad(3) = xq4m
            xquad(4) = xq4p
        end if

        ! bfvals_rsp:  positive rs-face
        ! bfvals_rsm:  negative rs-face
        ! bfvals_rsp(,1):  value=1 on positive rs-face
        ! bfvals_rsp(,kx):  x-values on positive rs-face
        ! bfvals_rsp(,ky):  y-values on positive rs-face
        ! bfvals_rsp(,kz):  z-values on positive rs-face
        ! bfvals_rsp(,kxx):  P_2(x)-values on positive rs-face
        ! bfvals_rsp(,kyy):  P_2(y)-values on positive rs-face
        ! bfvals_rsp(,kzz):  P_2(z)-values on positive rs-face
        ! bfvals_rsp(,kyz):  yz-values on positive rs-face
        ! bfvals_rsp(,kzx):  zx-values on positive rs-face
        ! bfvals_rsp(,kxy):  xy-values on positive rs-face

        bfvals_zp(1:nface,1) = 1.
        if (iquad .eq. 1) then
            bfvals_zp(1,kx) = 0.
            bfvals_zp(1,ky) = 0.
        end if

        if (iquad .gt. 1) then
            do i=1,nedge
                bfvals_zp((i-1)*nedge+1:i*nedge,kx) = xquad(1:nedge)
                bfvals_zp((i-1)*nedge+1:i*nedge,ky) = xquad(i)
            end do
        end if

        bfvals_zp(1:nface,kz) = 1.
        bfvals_zp(1:nface,kyz) = bfvals_zp(1:nface,ky)*bfvals_zp(1:nface,kz)
        bfvals_zp(1:nface,kzx) = bfvals_zp(1:nface,kz)*bfvals_zp(1:nface,kx)
        bfvals_zp(1:nface,kxy) = bfvals_zp(1:nface,kx)*bfvals_zp(1:nface,ky)
        bfvals_zp(1:nface,kxyz) = bfvals_zp(1:nface,kx)*bfvals_zp(1:nface,ky)*bfvals_zp(1:nface,kz)     ! basis function = xyz

        bfvals_yp(1:nface,1) = 1.
        bfvals_yp(1:nface,kx) = bfvals_zp(1:nface,ky)
        bfvals_yp(1:nface,ky) = 1.
        bfvals_yp(1:nface,kz) = bfvals_zp(1:nface,kx)
        bfvals_yp(1:nface,kyz) = bfvals_yp(1:nface,ky)*bfvals_yp(1:nface,kz)
        bfvals_yp(1:nface,kzx) = bfvals_yp(1:nface,kz)*bfvals_yp(1:nface,kx)
        bfvals_yp(1:nface,kxy) = bfvals_yp(1:nface,kx)*bfvals_yp(1:nface,ky)
        bfvals_yp(1:nface,kxyz) = bfvals_yp(1:nface,kx)*bfvals_yp(1:nface,ky)*bfvals_yp(1:nface,kz)     ! basis function = xyz

        bfvals_xp(1:nface,1) = 1.
        bfvals_xp(1:nface,kx) = 1.
        bfvals_xp(1:nface,ky) = bfvals_zp(1:nface,kx)
        bfvals_xp(1:nface,kz) = bfvals_zp(1:nface,ky)
        bfvals_xp(1:nface,kyz) = bfvals_xp(1:nface,ky)*bfvals_xp(1:nface,kz)
        bfvals_xp(1:nface,kzx) = bfvals_xp(1:nface,kz)*bfvals_xp(1:nface,kx)
        bfvals_xp(1:nface,kxy) = bfvals_xp(1:nface,kx)*bfvals_xp(1:nface,ky)
        bfvals_xp(1:nface,kxyz) = bfvals_xp(1:nface,kx)*bfvals_xp(1:nface,ky)*bfvals_xp(1:nface,kz)     ! basis function = xyz

        do i=1,3
            bfvals_xp(1:nface,kxx+i-1) = 1.5*bfvals_xp(1:nface,kx+i-1)**2 - 0.5
            bfvals_yp(1:nface,kxx+i-1) = 1.5*bfvals_yp(1:nface,kx+i-1)**2 - 0.5
            bfvals_zp(1:nface,kxx+i-1) = 1.5*bfvals_zp(1:nface,kx+i-1)**2 - 0.5
            ! bfvals_xp(1:nface,kxxx+i-1) = 2.5*bfvals_xp(1:nface,kx+i-1)**3 - 1.5*bfvals_xp(1:nface,kx+i-1)
            ! bfvals_yp(1:nface,kxxx+i-1) = 2.5*bfvals_yp(1:nface,kx+i-1)**3 - 1.5*bfvals_yp(1:nface,kx+i-1)
            ! bfvals_zp(1:nface,kxxx+i-1) = 2.5*bfvals_zp(1:nface,kx+i-1)**3 - 1.5*bfvals_zp(1:nface,kx+i-1)
        end do

        bfvals_xp(1:nface,kyyz) = bfvals_xp(1:nface,kyy)*bfvals_xp(1:nface,kz)     ! basis function = P_2(y)z
        bfvals_xp(1:nface,kyzz) = bfvals_xp(1:nface,ky)*bfvals_xp(1:nface,kzz)     ! basis function = P_2(z)y
        bfvals_xp(1:nface,kzzx) = bfvals_xp(1:nface,kzz)*bfvals_xp(1:nface,kx)     ! basis function = P_2(z)x
        bfvals_xp(1:nface,kzxx) = bfvals_xp(1:nface,kz)*bfvals_xp(1:nface,kxx)     ! basis function = P_2(x)z
        bfvals_xp(1:nface,kxxy) = bfvals_xp(1:nface,kxx)*bfvals_xp(1:nface,ky)     ! basis function = P_2(x)y
        bfvals_xp(1:nface,kxyy) = bfvals_xp(1:nface,kx)*bfvals_xp(1:nface,kyy)     ! basis function = P_2(y)x
        bfvals_xp(1:nface,kyyzz) = bfvals_xp(1:nface,kyy)*bfvals_xp(1:nface,kzz)     ! basis function = P_2(y)P_2(z)
        bfvals_xp(1:nface,kzzxx) = bfvals_xp(1:nface,kzz)*bfvals_xp(1:nface,kxx)     ! basis function = P_2(z)P_2(x)
        bfvals_xp(1:nface,kxxyy) = bfvals_xp(1:nface,kxx)*bfvals_xp(1:nface,kyy)     ! basis function = P_2(x)P_2(y)
        bfvals_xp(1:nface,kyzxx) = bfvals_xp(1:nface,kyz)*bfvals_xp(1:nface,kxx)     ! basis function = yz P_2(x)
        bfvals_xp(1:nface,kzxyy) = bfvals_xp(1:nface,kzx)*bfvals_xp(1:nface,kyy)     ! basis function = zx P_2(y)
        bfvals_xp(1:nface,kxyzz) = bfvals_xp(1:nface,kxy)*bfvals_xp(1:nface,kzz)     ! basis function = xy P_2(z)
        bfvals_xp(1:nface,kxyyzz) = bfvals_xp(1:nface,kx)*bfvals_xp(1:nface,kyy)*bfvals_xp(1:nface,kzz)   ! basis function = x P_2(y)P_2(z)
        bfvals_xp(1:nface,kyzzxx) = bfvals_xp(1:nface,ky)*bfvals_xp(1:nface,kzz)*bfvals_xp(1:nface,kxx)   ! basis function = y P_2(z)P_2(x)
        bfvals_xp(1:nface,kzxxyy) = bfvals_xp(1:nface,kz)*bfvals_xp(1:nface,kxx)*bfvals_xp(1:nface,kyy)   ! basis function = z P_2(x)P_2(y)
        bfvals_xp(1:nface,kxxyyzz) = bfvals_xp(1:nface,kx)*bfvals_xp(1:nface,kyy)*bfvals_xp(1:nface,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

        bfvals_xm = bfvals_xp
        bfvals_xm(1:nface,kx) = -bfvals_xp(1:nface,kx)
        bfvals_xm(1:nface,kzx) = -bfvals_xp(1:nface,kzx)
        bfvals_xm(1:nface,kxy) = -bfvals_xp(1:nface,kxy)
        bfvals_xm(1:nface,kxyz) = -bfvals_xp(1:nface,kxyz)
        bfvals_xm(1:nface,kzzx) = -bfvals_xp(1:nface,kzzx)
        bfvals_xm(1:nface,kxyy) = -bfvals_xp(1:nface,kxyy)
        ! bfvals_xm(1:nface,kxxx) = -bfvals_xp(1:nface,kxxx)
        bfvals_xm(1:nface,kzxyy) = -bfvals_xp(1:nface,kzxyy)
        bfvals_xm(1:nface,kxyzz) = -bfvals_xp(1:nface,kxyzz)
        bfvals_xm(1:nface,kxyyzz) = -bfvals_xp(1:nface,kxyyzz)

        bfvals_yp(1:nface,kyyz) = bfvals_yp(1:nface,kyy)*bfvals_yp(1:nface,kz)     ! basis function = P_2(y)z
        bfvals_yp(1:nface,kyzz) = bfvals_yp(1:nface,ky)*bfvals_yp(1:nface,kzz)     ! basis function = P_2(z)y
        bfvals_yp(1:nface,kzzx) = bfvals_yp(1:nface,kzz)*bfvals_yp(1:nface,kx)     ! basis function = P_2(z)x
        bfvals_yp(1:nface,kzxx) = bfvals_yp(1:nface,kz)*bfvals_yp(1:nface,kxx)     ! basis function = P_2(x)z
        bfvals_yp(1:nface,kxxy) = bfvals_yp(1:nface,kxx)*bfvals_yp(1:nface,ky)     ! basis function = P_2(x)y
        bfvals_yp(1:nface,kxyy) = bfvals_yp(1:nface,kx)*bfvals_yp(1:nface,kyy)     ! basis function = P_2(y)x
        bfvals_yp(1:nface,kyyzz) = bfvals_yp(1:nface,kyy)*bfvals_yp(1:nface,kzz)     ! basis function = P_2(y)P_2(z)
        bfvals_yp(1:nface,kzzxx) = bfvals_yp(1:nface,kzz)*bfvals_yp(1:nface,kxx)     ! basis function = P_2(z)P_2(x)
        bfvals_yp(1:nface,kxxyy) = bfvals_yp(1:nface,kxx)*bfvals_yp(1:nface,kyy)     ! basis function = P_2(x)P_2(y)
        bfvals_yp(1:nface,kyzxx) = bfvals_yp(1:nface,kyz)*bfvals_yp(1:nface,kxx)     ! basis function = yz P_2(x)
        bfvals_yp(1:nface,kzxyy) = bfvals_yp(1:nface,kzx)*bfvals_yp(1:nface,kyy)     ! basis function = zx P_2(y)
        bfvals_yp(1:nface,kxyzz) = bfvals_yp(1:nface,kxy)*bfvals_yp(1:nface,kzz)     ! basis function = xy P_2(z)
        bfvals_yp(1:nface,kxyyzz) = bfvals_yp(1:nface,kx)*bfvals_yp(1:nface,kyy)*bfvals_yp(1:nface,kzz)   ! basis function = x P_2(y)P_2(z)
        bfvals_yp(1:nface,kyzzxx) = bfvals_yp(1:nface,ky)*bfvals_yp(1:nface,kzz)*bfvals_yp(1:nface,kxx)   ! basis function = y P_2(z)P_2(x)
        bfvals_yp(1:nface,kzxxyy) = bfvals_yp(1:nface,kz)*bfvals_yp(1:nface,kxx)*bfvals_yp(1:nface,kyy)   ! basis function = z P_2(x)P_2(y)
        bfvals_yp(1:nface,kxxyyzz) = bfvals_yp(1:nface,kx)*bfvals_yp(1:nface,kyy)*bfvals_yp(1:nface,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

        bfvals_ym = bfvals_yp
        bfvals_ym(1:nface,ky) = -bfvals_yp(1:nface,ky)
        bfvals_ym(1:nface,kyz) = -bfvals_yp(1:nface,kyz)
        bfvals_ym(1:nface,kxy) = -bfvals_yp(1:nface,kxy)
        bfvals_ym(1:nface,kxyz) = -bfvals_yp(1:nface,kxyz)
        bfvals_ym(1:nface,kyzz) = -bfvals_yp(1:nface,kyzz)
        bfvals_ym(1:nface,kxxy) = -bfvals_yp(1:nface,kxxy)
        ! bfvals_ym(1:nface,kyyy) = -bfvals_yp(1:nface,kyyy)
        bfvals_ym(1:nface,kyzxx) = -bfvals_yp(1:nface,kyzxx)
        bfvals_ym(1:nface,kxyzz) = -bfvals_yp(1:nface,kxyzz)
        bfvals_ym(1:nface,kyzzxx) = -bfvals_yp(1:nface,kyzzxx)

        bfvals_zp(1:nface,kyyz) = bfvals_zp(1:nface,kyy)*bfvals_zp(1:nface,kz)     ! basis function = P_2(y)z
        bfvals_zp(1:nface,kyzz) = bfvals_zp(1:nface,ky)*bfvals_zp(1:nface,kzz)     ! basis function = P_2(z)y
        bfvals_zp(1:nface,kzzx) = bfvals_zp(1:nface,kzz)*bfvals_zp(1:nface,kx)     ! basis function = P_2(z)x
        bfvals_zp(1:nface,kzxx) = bfvals_zp(1:nface,kz)*bfvals_zp(1:nface,kxx)     ! basis function = P_2(x)z
        bfvals_zp(1:nface,kxxy) = bfvals_zp(1:nface,kxx)*bfvals_zp(1:nface,ky)     ! basis function = P_2(x)y
        bfvals_zp(1:nface,kxyy) = bfvals_zp(1:nface,kx)*bfvals_zp(1:nface,kyy)     ! basis function = P_2(y)x
        bfvals_zp(1:nface,kyyzz) = bfvals_zp(1:nface,kyy)*bfvals_zp(1:nface,kzz)     ! basis function = P_2(y)P_2(z)
        bfvals_zp(1:nface,kzzxx) = bfvals_zp(1:nface,kzz)*bfvals_zp(1:nface,kxx)     ! basis function = P_2(z)P_2(x)
        bfvals_zp(1:nface,kxxyy) = bfvals_zp(1:nface,kxx)*bfvals_zp(1:nface,kyy)     ! basis function = P_2(x)P_2(y)
        bfvals_zp(1:nface,kyzxx) = bfvals_zp(1:nface,kyz)*bfvals_zp(1:nface,kxx)     ! basis function = yz P_2(x)
        bfvals_zp(1:nface,kzxyy) = bfvals_zp(1:nface,kzx)*bfvals_zp(1:nface,kyy)     ! basis function = zx P_2(y)
        bfvals_zp(1:nface,kxyzz) = bfvals_zp(1:nface,kxy)*bfvals_zp(1:nface,kzz)     ! basis function = xy P_2(z)
        bfvals_zp(1:nface,kxyyzz) = bfvals_zp(1:nface,kx)*bfvals_zp(1:nface,kyy)*bfvals_zp(1:nface,kzz)   ! basis function = x P_2(y)P_2(z)
        bfvals_zp(1:nface,kyzzxx) = bfvals_zp(1:nface,ky)*bfvals_zp(1:nface,kzz)*bfvals_zp(1:nface,kxx)   ! basis function = y P_2(z)P_2(x)
        bfvals_zp(1:nface,kzxxyy) = bfvals_zp(1:nface,kz)*bfvals_zp(1:nface,kxx)*bfvals_zp(1:nface,kyy)   ! basis function = z P_2(x)P_2(y)
        bfvals_zp(1:nface,kxxyyzz) = bfvals_zp(1:nface,kx)*bfvals_zp(1:nface,kyy)*bfvals_zp(1:nface,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

        bfvals_zm = bfvals_zp
        bfvals_zm(1:nface,kz) = -bfvals_zp(1:nface,kz)
        bfvals_zm(1:nface,kyz) = -bfvals_zp(1:nface,kyz)
        bfvals_zm(1:nface,kzx) = -bfvals_zp(1:nface,kzx)
        bfvals_zm(1:nface,kxyz) = -bfvals_zp(1:nface,kxyz)
        bfvals_zm(1:nface,kyyz) = -bfvals_zp(1:nface,kyyz)
        bfvals_zm(1:nface,kzxx) = -bfvals_zp(1:nface,kzxx)
        ! bfvals_zm(1:nface,kzzz) = -bfvals_zp(1:nface,kzzz)
        bfvals_zm(1:nface,kyzxx) = -bfvals_zp(1:nface,kyzxx)
        bfvals_zm(1:nface,kzxyy) = -bfvals_zp(1:nface,kzxyy)
        bfvals_zm(1:nface,kzxxyy) = -bfvals_zp(1:nface,kzxxyy)

        ! Organize local basis values on faces into 1-D vectors.
        ! Used in limiter() and max_lim().

        bf_faces(1:nslim,1) = 1.        ! basis function = 1

        do ixyz=kx,kz
            bf_faces(1:nface,ixyz) = bfvals_xm(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(nface+1:2*nface,ixyz) = bfvals_xp(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(2*nface+1:3*nface,ixyz) = bfvals_ym(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(3*nface+1:4*nface,ixyz) = bfvals_yp(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(4*nface+1:5*nface,ixyz) = bfvals_zm(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(5*nface+1:6*nface,ixyz) = bfvals_zp(1:nface,ixyz)        ! basis function = x,y,z
            bf_faces(6*nface+1:nslim,ixyz) = bfvals_int(1:npg,ixyz)        ! basis function = x,y,z
        end do

        bf_faces(1:nslim,kyz) = bf_faces(1:nslim,ky)*bf_faces(1:nslim,kz)     ! basis function = yz
        bf_faces(1:nslim,kzx) = bf_faces(1:nslim,kz)*bf_faces(1:nslim,kx)     ! basis function = zx
        bf_faces(1:nslim,kxy) = bf_faces(1:nslim,kx)*bf_faces(1:nslim,ky)     ! basis function = xy
        bf_faces(1:nslim,kxyz) = bf_faces(1:nslim,kx)*bf_faces(1:nslim,ky)*bf_faces(1:nslim,kz)     ! basis function = xyz

        do i=0,2
            bf_faces(1:nslim,kxx+i) = 1.5*bf_faces(1:nslim,kx+i)**2 - 0.5    ! basis function = P_2(s)
            ! bf_faces(1:nslim,kxxx+i) = 2.5*bf_faces(1:nslim,kx+i)**3 - 1.5*bf_faces(1:nslim,kx+i)   ! basis function = P_3(s)
        end do

        bf_faces(1:nslim,kyyz) = bf_faces(1:nslim,kyy)*bf_faces(1:nslim,kz)     ! basis function = P_2(y)z
        bf_faces(1:nslim,kyzz) = bf_faces(1:nslim,ky)*bf_faces(1:nslim,kzz)     ! basis function = P_2(z)y
        bf_faces(1:nslim,kzzx) = bf_faces(1:nslim,kzz)*bf_faces(1:nslim,kx)     ! basis function = P_2(z)x
        bf_faces(1:nslim,kzxx) = bf_faces(1:nslim,kz)*bf_faces(1:nslim,kxx)     ! basis function = P_2(x)z
        bf_faces(1:nslim,kxxy) = bf_faces(1:nslim,kxx)*bf_faces(1:nslim,ky)     ! basis function = P_2(x)y
        bf_faces(1:nslim,kxyy) = bf_faces(1:nslim,kx)*bf_faces(1:nslim,kyy)     ! basis function = P_2(y)x
        bf_faces(1:nslim,kyyzz) = bf_faces(1:nslim,kyy)*bf_faces(1:nslim,kzz)     ! basis function = P_2(y)P_2(z)
        bf_faces(1:nslim,kzzxx) = bf_faces(1:nslim,kzz)*bf_faces(1:nslim,kxx)     ! basis function = P_2(z)P_2(x)
        bf_faces(1:nslim,kxxyy) = bf_faces(1:nslim,kxx)*bf_faces(1:nslim,kyy)     ! basis function = P_2(x)P_2(y)
        bf_faces(1:nslim,kyzxx) = bf_faces(1:nslim,kyz)*bf_faces(1:nslim,kxx)     ! basis function = yz P_2(x)
        bf_faces(1:nslim,kzxyy) = bf_faces(1:nslim,kzx)*bf_faces(1:nslim,kyy)     ! basis function = zx P_2(y)
        bf_faces(1:nslim,kxyzz) = bf_faces(1:nslim,kxy)*bf_faces(1:nslim,kzz)     ! basis function = xy P_2(z)
        bf_faces(1:nslim,kxyyzz) = bf_faces(1:nslim,kx)*bf_faces(1:nslim,kyy)*bf_faces(1:nslim,kzz)   ! basis function = x P_2(y)P_2(z)
        bf_faces(1:nslim,kyzzxx) = bf_faces(1:nslim,ky)*bf_faces(1:nslim,kzz)*bf_faces(1:nslim,kxx)   ! basis function = y P_2(z)P_2(x)
        bf_faces(1:nslim,kzxxyy) = bf_faces(1:nslim,kz)*bf_faces(1:nslim,kxx)*bf_faces(1:nslim,kyy)   ! basis function = z P_2(x)P_2(y)
        bf_faces(1:nslim,kxxyyzz) = bf_faces(1:nslim,kx)*bf_faces(1:nslim,kyy)*bf_faces(1:nslim,kzz)  ! basis function = P_2(x)P_2(y)P_2(z)

    end subroutine set_face_vals_3D

!----------------------------------------------------------
    subroutine set_weights_3D

        integer i
        real wq4p,wq4m

        wq4p = (18. + sqrt(30.))/36.
        wq4m = (18. - sqrt(30.))/36.

        ! Define weights for integral approximation using Gaussian quadrature.

        ! Define weights for 1-D integration
        if (iquad .eq. 1) then        ! 1-point quadrature
            wgt1d(1) = 2.
        end if
        if (iquad .eq. 2) then        ! 2-point quadrature
            wgt1d(1:2) = 1.
        end if
        if (iquad .eq. 3) then        ! 3-point quadrature
            wgt1d(1) = 5./9.
            wgt1d(2) = 8./9.
            wgt1d(3) = 5./9.
        end if
        if (iquad .eq. 4) then        ! 4-point quadrature
            wgt1d(1) = wq4m
            wgt1d(2) = wq4p
            wgt1d(3) = wq4p
            wgt1d(4) = wq4m
        end if

        ! Define weights for 2-D integration
        if (iquad .eq. 1) then        ! 1-point quadrature
            wgt2d(1) = 4.
        end if
        if (iquad .eq. 2) then        ! 2-point quadrature
            wgt2d(1:4) = 1.
        end if

        if (iquad .ge. 3) then
            do i= 1,nedge
                wgt2d((i-1)*nedge+1:i*nedge) = wgt1d(1:nedge)*wgt1d(i)
            end do
        end if

        ! Define weights for 3-D integration
        if (iquad .eq. 1) then        ! 1-point quadrature
            wgt3d(1) = 8.
        end if
        if (iquad .eq. 2) then        ! 2-point quadrature
            wgt3d(1:8) = 1.
        end if

        if (iquad .ge. 3) then
            do i= 1,nedge
                wgt3d((i-1)*nface+1:i*nface) = wgt2d(1:nface)*wgt1d(i)
            end do
        end if

    end subroutine set_weights_3D

    !===========================================================================
    ! This routine tests whether the basis functions are orthonormal over the
    ! cell volume and over each cell face.
    !---------------------------------------------------------------------------
    subroutine test_basis_3D
        implicit none
        integer ir,jr,ix1,iy1,iz1,ipass
        real sq3,sq5,sq7,epsm

        sq3 = sqrt(3.)
        sq5 = sqrt(5.)
        sq7 = sqrt(7.)

        ! Indices of basis elements projected onto the x-faces (yz-plane):
        ibas_x(1) = 1
        ibas_x(kx) = 1
        ibas_x(ky) = ky
        ibas_x(kz) = kz
        ibas_x(kxy) = ky
        ibas_x(kyz) = kyz
        ibas_x(kzx) = kz
        ibas_x(kxx) = 1
        ibas_x(kzz) = kzz
        ibas_x(kyy) = kyy
        ibas_x(kxxy) = ky
        ibas_x(kyyz) = kyyz
        ibas_x(kzzx) = kzz
        ibas_x(kxyy) = kyy
        ibas_x(kyzz) = kyzz
        ibas_x(kzxx) = kz
        ibas_x(kxyz) = kyz
        ibas_x(kxxyy) = kyy
        ibas_x(kyyzz) = kyyzz
        ibas_x(kzzxx) = kzz
        ibas_x(kxyzz) = kyzz
        ibas_x(kyzxx) = kyz
        ibas_x(kzxyy) = kyyz
        ibas_x(kxyyzz) = kyyzz
        ibas_x(kyzzxx) = kyzz
        ibas_x(kzxxyy) = kyyz
        ibas_x(kxxyyzz) = kyyzz

        ! Indices of basis elements projected onto the y-faces (zx-plane):
        ibas_y(1) = 1
        ibas_y(kx) = kx
        ibas_y(ky) = 1
        ibas_y(kz) = kz
        ibas_y(kxy) = kx
        ibas_y(kyz) = kz
        ibas_y(kzx) = kzx
        ibas_y(kxx) = kxx
        ibas_y(kzz) = kzz
        ibas_y(kyy) = 1
        ibas_y(kxxy) = kxx
        ibas_y(kyyz) = kz
        ibas_y(kzzx) = kzzx
        ibas_y(kxyy) = kx
        ibas_y(kyzz) = kzz
        ibas_y(kzxx) = kzxx
        ibas_y(kxyz) = kzx
        ibas_y(kxxyy) = kxx
        ibas_y(kyyzz) = kzz
        ibas_y(kzzxx) = kzzxx
        ibas_y(kxyzz) = kzzx
        ibas_y(kyzxx) = kzxx
        ibas_y(kzxyy) = kzx
        ibas_y(kxyyzz) = kzzx
        ibas_y(kyzzxx) = kzzxx
        ibas_y(kzxxyy) = kzxx
        ibas_y(kxxyyzz) = kzzxx

        ! Indices of basis elements projected onto the z-faces (xy-plane):
        ibas_z(1) = 1
        ibas_z(kx) = kx
        ibas_z(ky) = ky
        ibas_z(kz) = 1
        ibas_z(kxy) = kxy
        ibas_z(kyz) = ky
        ibas_z(kzx) = kx
        ibas_z(kxx) = kxx
        ibas_z(kzz) = 1
        ibas_z(kyy) = kyy
        ibas_z(kxxy) = kxxy
        ibas_z(kyyz) = kyy
        ibas_z(kzzx) = kx
        ibas_z(kxyy) = kxyy
        ibas_z(kyzz) = ky
        ibas_z(kzxx) = kxx
        ibas_z(kxyz) = kxy
        ibas_z(kxxyy) = kxxyy
        ibas_z(kyyzz) = kyy
        ibas_z(kzzxx) = kxx
        ibas_z(kxyzz) = kxy
        ibas_z(kyzxx) = kxxy
        ibas_z(kzxyy) = kxyy
        ibas_z(kxyyzz) = kxyy
        ibas_z(kyzzxx) = kxxy
        ibas_z(kzxxyy) = kxxyy
        ibas_z(kxxyyzz) = kxxyy

        ! Normalization of basis elements over the  positive x-face of a cell:
        cbas_xp(1) = 1.
        cbas_xp(kx) = 1.
        cbas_xp(ky) = sq3
        cbas_xp(kz) = sq3
        cbas_xp(kxy) = sq3
        cbas_xp(kyz) = 3.
        cbas_xp(kzx) = sq3
        cbas_xp(kxx) = 1.
        cbas_xp(kzz) = sq5
        cbas_xp(kyy) = sq5
        cbas_xp(kxxy) = sq3
        cbas_xp(kyyz) = sq5*sq3
        cbas_xp(kzzx) = sq5
        cbas_xp(kxyy) = sq5
        cbas_xp(kyzz) = sq3*sq5
        cbas_xp(kzxx) = sq3
        cbas_xp(kxyz) = sq3*sq3
        cbas_xp(kxxyy) = sq5
        cbas_xp(kyyzz) = sq5*sq5
        cbas_xp(kzzxx) = sq5
        cbas_xp(kxyzz) = sq3*sq5
        cbas_xp(kyzxx) = sq3*sq3
        cbas_xp(kzxyy) = sq3*sq5
        cbas_xp(kxyyzz) = sq5*sq5
        cbas_xp(kyzzxx) = sq3*sq5
        cbas_xp(kzxxyy) = sq3*sq5
        cbas_xp(kxxyyzz) = sq5*sq5

        ! Normalization of basis elements over the  positive y-face of a cell:
        cbas_yp(1) = 1.
        cbas_yp(kx) = sq3
        cbas_yp(ky) = 1.
        cbas_yp(kz) = sq3
        cbas_yp(kxy) = sq3
        cbas_yp(kyz) = sq3
        cbas_yp(kzx) = 3.
        cbas_yp(kxx) = sq5
        cbas_yp(kzz) = sq5
        cbas_yp(kyy) = 1.
        cbas_yp(kxxy) = sq5
        cbas_yp(kyyz) = sq3
        cbas_yp(kzzx) = sq5*sq3
        cbas_yp(kxyy) = sq3
        cbas_yp(kyzz) = sq5
        cbas_yp(kzxx) = sq3*sq5
        cbas_yp(kxyz) = sq3*sq3
        cbas_yp(kxxyy) = sq5
        cbas_yp(kyyzz) = sq5
        cbas_yp(kzzxx) = sq5*sq5
        cbas_yp(kxyzz) = sq3*sq5
        cbas_yp(kyzxx) = sq3*sq5
        cbas_yp(kzxyy) = sq3*sq3
        cbas_yp(kxyyzz) = sq3*sq5
        cbas_yp(kyzzxx) = sq5*sq5
        cbas_yp(kzxxyy) = sq3*sq5
        cbas_yp(kxxyyzz) = sq5*sq5

        ! Normalization of basis elements over the  positive z-face of a cell:
        cbas_zp(1) = 1.
        cbas_zp(kx) = sq3
        cbas_zp(ky) = sq3
        cbas_zp(kz) = 1.
        cbas_zp(kxy) = 3.
        cbas_zp(kyz) = sq3
        cbas_zp(kzx) = sq3
        cbas_zp(kxx) = sq5
        cbas_zp(kzz) = 1.
        cbas_zp(kyy) = sq5
        cbas_zp(kxxy) = sq5*sq3
        cbas_zp(kyyz) = sq5
        cbas_zp(kzzx) = sq3
        cbas_zp(kxyy) = sq3*sq5
        cbas_zp(kyzz) = sq3
        cbas_zp(kzxx) = sq5
        cbas_zp(kxyz) = sq3*sq3
        cbas_zp(kxxyy) = sq5*sq5
        cbas_zp(kyyzz) = sq5
        cbas_zp(kzzxx) = sq5
        cbas_zp(kxyzz) = sq3*sq3
        cbas_zp(kyzxx) = sq3*sq5
        cbas_zp(kzxyy) = sq3*sq5
        cbas_zp(kxyyzz) = sq3*sq5
        cbas_zp(kyzzxx) = sq3*sq5
        cbas_zp(kzxxyy) = sq5*sq5
        cbas_zp(kxxyyzz) = sq5*sq5

        ! Normalization of basis elements over the  negative x-face of a cell:
        cbas_xm = cbas_xp
        cbas_xm(kx) = -1.
        cbas_xm(kxy) = -sq3
        cbas_xm(kzx) = -sq3
        cbas_xm(kxyy) = -sq5
        cbas_xm(kzzx) = -sq5
        cbas_xm(kxyz) = -sq3*sq3
        cbas_xm(kxyzz) = -sq3*sq5
        cbas_xm(kzxyy) = -sq3*sq5
        cbas_xm(kxyyzz) = -sq5*sq5

        ! Normalization of basis elements over the  negative y-face of a cell:
        cbas_ym = cbas_yp
        cbas_ym(ky) = -1.
        cbas_ym(kxy) = -sq3
        cbas_ym(kyz) = -sq3
        cbas_ym(kyzz) = -sq5
        cbas_ym(kxxy) = -sq5
        cbas_ym(kxyz) = -sq3*sq3
        cbas_ym(kxyzz) = -sq3*sq5
        cbas_ym(kyzxx) = -sq3*sq5
        cbas_ym(kyzzxx) = -sq5*sq5

        ! Normalization of basis elements over the  negative z-face of a cell:
        cbas_zm = cbas_zp
        cbas_zm(kz) = -1.
        cbas_zm(kyz) = -sq3
        cbas_zm(kzx) = -sq3
        cbas_zm(kzxx) = -sq5
        cbas_zm(kyyz) = -sq5
        cbas_zm(kxyz) = -sq3*sq3
        cbas_zm(kyzxx) = -sq3*sq5
        cbas_zm(kzxyy) = -sq3*sq5
        cbas_zm(kzxxyy) = -sq5*sq5

        ! Below, we specify what the scalar product of any two basis elements should
        ! evaluate to.

        ! cell_int0: scalar product over cell volume
        ! xp_int0: scalar product over positive x-face
        ! xm_int0: scalar product over negative x-face
        ! yp_int0: scalar product over positive y-face
        ! ym_int0: scalar product over negative y-face
        ! zp_int0: scalar product over positive z-face
        ! zm_int0: scalar product over negative z-face

        xp_int0 = 0.
        xm_int0 = 0.
        yp_int0 = 0.
        ym_int0 = 0.
        zp_int0 = 0.
        zm_int0 = 0.

        do jr=1,nbasis
            do ir=1,nbasis
                if (ir .ne. jr) cell_int0(ir,jr) = 0.
                if (ir .eq. jr) cell_int0(ir,jr) = 1.
                ix1 = 0
                iy1 = 0
                iz1 = 0

                if (ibas_x(ir) .eq. 1 .and. ibas_x(jr) .eq. 1) ix1 = 1
                if (ibas_x(ir) .eq. ky .and. ibas_x(jr) .eq. ky) ix1 = 1
                if (ibas_x(ir) .eq. kyy .and. ibas_x(jr) .eq. kyy) ix1 = 1
                if (ibas_x(ir) .eq. kz .and. ibas_x(jr) .eq. kz) ix1 = 1
                if (ibas_x(ir) .eq. kzz .and. ibas_x(jr) .eq. kzz) ix1 = 1
                if (ibas_x(ir) .eq. kyz .and. ibas_x(jr) .eq. kyz) ix1 = 1
                if (ibas_x(ir) .eq. kyyz .and. ibas_x(jr) .eq. kyyz) ix1 = 1
                if (ibas_x(ir) .eq. kyzz .and. ibas_x(jr) .eq. kyzz) ix1 = 1
                if (ibas_x(ir) .eq. kyyzz .and. ibas_x(jr) .eq. kyyzz) ix1 = 1

                if (ibas_y(ir) .eq. 1 .and. ibas_y(jr) .eq. 1) iy1 = 1
                if (ibas_y(ir) .eq. kx .and. ibas_y(jr) .eq. kx) iy1 = 1
                if (ibas_y(ir) .eq. kxx .and. ibas_y(jr) .eq. kxx) iy1 = 1
                if (ibas_y(ir) .eq. kz .and. ibas_y(jr) .eq. kz) iy1 = 1
                if (ibas_y(ir) .eq. kzz .and. ibas_y(jr) .eq. kzz) iy1 = 1
                if (ibas_y(ir) .eq. kzx .and. ibas_y(jr) .eq. kzx) iy1 = 1
                if (ibas_y(ir) .eq. kzxx .and. ibas_y(jr) .eq. kzxx) iy1 = 1
                if (ibas_y(ir) .eq. kzzx .and. ibas_y(jr) .eq. kzzx) iy1 = 1
                if (ibas_y(ir) .eq. kzzxx .and. ibas_y(jr) .eq. kzzxx) iy1 = 1

                if (ibas_z(ir) .eq. 1 .and. ibas_z(jr) .eq. 1) iz1 = 1
                if (ibas_z(ir) .eq. kx .and. ibas_z(jr) .eq. kx) iz1 = 1
                if (ibas_z(ir) .eq. kxx .and. ibas_z(jr) .eq. kxx) iz1 = 1
                if (ibas_z(ir) .eq. ky .and. ibas_z(jr) .eq. ky) iz1 = 1
                if (ibas_z(ir) .eq. kyy .and. ibas_z(jr) .eq. kyy) iz1 = 1
                if (ibas_z(ir) .eq. kxy .and. ibas_z(jr) .eq. kxy) iz1 = 1
                if (ibas_z(ir) .eq. kxxy .and. ibas_z(jr) .eq. kxxy) iz1 = 1
                if (ibas_z(ir) .eq. kxyy .and. ibas_z(jr) .eq. kxyy) iz1 = 1
                if (ibas_z(ir) .eq. kxxyy .and. ibas_z(jr) .eq. kxxyy) iz1 = 1

                if (ix1 .eq. 1) then
                    xp_int0(ir,jr) = 1.  !0.5*3.
                    xm_int0(ir,jr) = 1.  !0.5*3.
                end if
                if (iy1 .eq. 1) then
                    yp_int0(ir,jr) = 1.  !0.5*3.
                    ym_int0(ir,jr) = 1.  !0.5*3.
                end if
                if (iz1 .eq. 1) then
                    zp_int0(ir,jr) = 1.  !0.5*3.
                    zm_int0(ir,jr) = 1.  !0.5*3.
                end if
           end do
        end do

        ! Below, we evaluate the scalar product of all pairs of basis elements and
        ! test whether it's equal to the expected value computed above.

        ! cell_int: scalar product over cell volume
        ! xp_int: scalar product over positive x-face
        ! xm_int: scalar product over negative x-face
        ! yp_int: scalar product over positive y-face
        ! ym_int: scalar product over negative y-face
        ! zp_int: scalar product over positive z-face
        ! zm_int: scalar product over negative z-face

        ipass = 1

        do jr=1,nbasis
            do ir=1,nbasis
                cell_int(ir,jr) = 0.125*cbasis(ir)*sum(wgt3d(1:npg)*bfvals_int(1:npg,ir)*bfvals_int(1:npg,jr))
                xp_int(ir,jr) = 0.25*cbas_xp(ir)*cbas_xp(jr)*sum(wgt2d(1:nface)*bfvals_xp(1:nface,ir)*bfvals_xp(1:nface,jr))
                xm_int(ir,jr) = 0.25*cbas_xm(ir)*cbas_xm(jr)*sum(wgt2d(1:nface)*bfvals_xm(1:nface,ir)*bfvals_xm(1:nface,jr))
                yp_int(ir,jr) = 0.25*cbas_yp(ir)*cbas_yp(jr)*sum(wgt2d(1:nface)*bfvals_yp(1:nface,ir)*bfvals_yp(1:nface,jr))
                ym_int(ir,jr) = 0.25*cbas_ym(ir)*cbas_ym(jr)*sum(wgt2d(1:nface)*bfvals_ym(1:nface,ir)*bfvals_ym(1:nface,jr))
                zp_int(ir,jr) = 0.25*cbas_zp(ir)*cbas_zp(jr)*sum(wgt2d(1:nface)*bfvals_zp(1:nface,ir)*bfvals_zp(1:nface,jr))
                zm_int(ir,jr) = 0.25*cbas_zm(ir)*cbas_zm(jr)*sum(wgt2d(1:nface)*bfvals_zm(1:nface,ir)*bfvals_zm(1:nface,jr))

                epsm = 1.e-5

         !     if (iam .eq. print_mpi) then
            if (abs(cell_int(ir,jr) - cell_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect cell integral: bases',ir,jr
                print *,'computed value: ',cell_int(ir,jr),'  expected value: ',cell_int0(ir,jr)
                ipass = 0
            end if
                if (abs(xp_int(ir,jr) - xp_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect +x-face integral: bases',ir,jr
                print *,'computed value: ',xp_int(ir,jr),'  expected value: ',xp_int0(ir,jr)
                ipass = 0
            end if
                if (abs(xm_int(ir,jr) - xm_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect -x-face integral: bases',ir,jr
                print *,'computed value: ',xm_int(ir,jr),'  expected value: ',xm_int0(ir,jr)
                ipass = 0
            end if
                if (abs(yp_int(ir,jr) - yp_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect +y-face integral: bases',ir,jr
                print *,'computed value: ',yp_int(ir,jr),'  expected value: ',yp_int0(ir,jr)
                ipass = 0
            end if
                if (abs(ym_int(ir,jr) - ym_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect -y-face integral: bases',ir,jr
                print *,'computed value: ',ym_int(ir,jr),'  expected value: ',ym_int0(ir,jr)
                ipass = 0
            end if
                if (abs(zp_int(ir,jr) - zp_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect +z-face integral: bases',ir,jr
                print *,'computed value: ',zp_int(ir,jr),'  expected value: ',zp_int0(ir,jr)
                ipass = 0
            end if
                if (abs(zm_int(ir,jr) - zm_int0(ir,jr)) .gt. epsm) then
                print *,'incorrect -z-face integral: bases',ir,jr
                print *,'computed value: ',zm_int(ir,jr),'  expected value: ',zm_int0(ir,jr)
                ipass = 0
            end if
          !    end if
           end do
        end do

        ! if (ipass .eq. 0) call exit(-1)
    end subroutine test_basis_3D
    !---------------------------------------------------------------------------

end module basis_funcs

max_stars_repo_name: sseyler/HERMESHD

max_issues_repo_name: sseyler/HERMESHD

max_forks_repo_name: sseyler/HERMESHD


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

c      external cgit07
c      external cgimet
      character*10 met1,met2
      integer mf1id, tk1id
      real x(20), y(20)
      integer clr(3)
      integer colors(24)
      character*10 text
      integer lnclr(3)
      integer chhit(10)

      include 'cgidef.f'

c set colors up like vdi:
C           Index  Color  RGB Values
C             0    black   0.,0.,0.
C             1    red     1.,0.,0.
C             2    green   0.,1.,0.
C             3    yellow  1.,1.,0.
C             4    blue    0.,0.,1.
C             5    magenta 1.,0.,1.
C             6    cyan    0.,1.,1.
C             7    white   1.,1.,1.

      data colors / 0, 0, 0,  255, 0, 0, 0, 255, 0,
     *               255, 255, 0, 0, 0, 255,  255, 0, 255,
     *               0, 255, 255,  255, 255, 255 /

      met1 = 'file55'
      met2 = 'meta2'

c      call xcact(cgimet,mf1id)
c      call xcact(cgimet,mf2id)
c      call xcact(cgit07,tk1id)

c      call xcooff(tk1id)
c      call xcooff(mf2id)
c      call cesc(-28372,1,met1)
c      call xcooff(mf1id)
c      call xcoon(mf2id)
c      call cesc(-28372,1,met2)
c      call xcoon(tk1id)
c      call xcoon(mf1id)

c initialize cgi
c      call cesc(-28372,1,met1)
      call ci(1)
      call cvdcx(0., 0., 32767.,32767.)
c      call cv(.5,0.,2.5,1.)
C VDC clip COFF, CON
      call ccl(CON)
C viewsurface clip CDCOFF, CDCREC, CVPORT
      call cdscl(CDCREC)

      do 10 i=1,2
c draw some primitives in default color
c  ..polylines
      x(1) = 100.
      x(2) = 16000.
      y(1) = 16500.
      y(2) = 16500.
      call clnt(1)
      call cpl(2,x,y)
      y(1) = 20500.
      y(2) = 20500.
      call clnt(2)
      call cpl(2,x,y)
      y(1) = 24500.
      y(2) = 24500.
      call clnt(3)
      call cpl(2,x,y)
      y(1) = 28500.
      y(2) = 28500.
      call clnt(4)
      call cpl(2,x,y)

c  ..polygons
      call cis(CHOLLO)
      x(1) = 20000.
      y(1) = 25000.
      x(2) = 29000.
      y(2) = 25000.
      x(3) = 24500.
      y(3) = 32000
      x(4) = 20000.
      y(4) = 25000.
      call cpg(4,x,y)
      x(1) = 20000.
      y(1) = 23500.
      x(2) = 29000.
      y(2) = 23500.
      x(3) = 24500.
      y(3) = 16500.
      x(4) = 20000.
      y(4) = 23500.
      call cis(CSOLID)
      call cpg(4,x,y)

c ...markers
      x(1) = 100.
      y(1) = 100.
      x(2) = 2100.
      y(2) = 2100.
      x(3) = 4100.
      y(3) = 4100.
      x(4) = 6100.
      y(4) = 6100
      x(5) = 8100.
      y(5) = 8100.
      x(6) = 10100.
      y(6) = 10100.
      x(7) = 12100.
      y(7) = 12100.
      x(8) = 14100.
      y(8) = 14100.
      x(9) = 16100.
      y(9) = 16100.
      call cpm(9,x,y)

c ...text
c ...do tektronix first, then metafile
c ...do this because cgtxx is an inquiry
c      do 20 j=1,2
c      if(j.eq.1)then
c        call xcooff( mf1id )
c        call xcooff( mf2id )
c        call xcoon( tk1id )
c        call xcsol( tk1id )
c        call cili( CLOCAT, 1)
c      else
c        call xcooff( tk1id )
c        call xcoon( mf1id )
c        call xcoon( mf2id )
c        call xcsol( mf1id )
        call cili( CLOCAT, 1)
c     endif
      call cqchh(" ",0,1,1,istat,ntot,nlist,chhit)
      text = "TEXT"
      call cchh(327.)
      call ctx( 17000.,100.,1,text(1:4))
      call cgtxx(17000.,100.,text(1:4),ivstat,ivconc,xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      call ctx( xconc,yconc,1,text(1:4))
      call cgtxx(xconc,yconc,text(1:4),ivstat,ivconc, xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      call cchh(3270.)
      call ctx( xconc,yconc,1,text(1:4))

      call cchh(3270.)
      text = "A"
      call ctx( 16500.,10000.,1,text(1:1))
      call cgtxx(16500.,10000.,text(1:1),ivstat,ivconc, xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      text = "B"
      call ctx( xconc,yconc,1,text(1:1))
      call cgtxx(xconc,yconc,text(1:1),ivstat,ivconc, xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      text = "C"
      call ctx( xconc,yconc,1,text(1:1))
      call cgtxx(xconc,yconc,text(1:1),ivstat,ivconc, xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      text = "ABC"
      call ctx( xconc, yconc,1,text(1:3))
      call cgtxx(xconc,yconc,text(1:3),ivstat,ivconc, xconc,yconc,
     *  p1,q1,p2,q2,p3,q3,p4,q4)
      x(1) = p1
      y(1) = q1
      x(2) = p2
      y(2) = q2
      x(3) = p3
      y(3) = q3
      x(4) = p4
      y(4) = q4
      x(5) = p1
      y(5) = q1
      call clnt(1)
      call cpl(5,x,y)

c use request locator as a pause
      call crqlc(1,1,istat,irstat,mvalid,itrig,xx,yy)
      call cpds(0)
20    continue

c turn everything on
c      call xcoon( tk1id )
c      call xcoon( mf1id )
c      call xcoon( mf2id )

c change colors
      if(i.eq.1) then

c does the device support direct color
c ...do it this wierd way cause i'm lazy - if one is only indexed
c     color, then treat them both as indexed color
c        call xcsol( tk1id )
        call cqc( idum, idum, idum, idum, icmode, idum, idum, idum )
        if( icmode .eq. CCLRID ) then
c          call xcsol( mf1id )
          call cqc( idum, idum, idum, idum, icmode, idum, idum, idum )
        endif

c set direct color mode if supported, otherwise use indexed
c only set the color table for indexed color
        if( icmode .eq. CCLRID ) then
          call ccsm( CDRECT )
        else
          call ccsm( CINDEX )
          call cct(0,8,colors)
        endif

c  ...make polylines red
        if( icmode .eq. CCLRID ) then
          clr(1) = 255
          clr(2) = 0
          clr(3) = 0
        else
          clr(1) = 1
        endif
        call clnc(clr)
        call cqlna(istat,lnbi,lntyp,lwmod,lnwid,csmod,lnclr,lcmod)

c  ...make polygons green
        if( icmode .eq. CCLRID ) then
          clr(1) = 0
          clr(2) = 255
          clr(3) = 0
        else
          clr(1) = 2
        endif
        call cflc(clr)

c  ...make polymarkers yellow
        if( icmode .eq. CCLRID ) then
          clr(1) = 255
          clr(2) = 255
          clr(3) = 0
        else
          clr(1) = 3
        endif
        call cmkc(clr)

c  ...make text blue
        if( icmode .eq. CCLRID ) then
          clr(1) = 0
          clr(2) = 0
          clr(3) = 255
        else
          clr(1) = 4
        endif
        call ctxc(clr)
      endif

10    continue
      call ct
c      call xcdact( mf1id )
c      call xcdact( mf2id )
c      call xcdact( tk1id )
      stop
      end

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
Jeff Ito is all about wholesome, family values, peace, and love for all mankind.

Signed up 200702

Jeffs cool and has a hot fiancee Ive been trying to steal. Hi Jeff!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program 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 General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine allocont(ncont,ntie,tieset,nset,set,istartset,
     &  iendset,ialset,lakon,ncone,tietol,ismallsliding,kind1,kind2,
     &  mortar,istep)
!
!     counting the number of triangles needed for the
!     triangulation of the contact master surfaces
!
!     ismallsliding = 0: large sliding
!                   = 1: small sliding
!
      implicit none
!
      logical nodeslavsurf
!
      character*1 kind1,kind2
      character*8 lakon(*)
      character*81 tieset(3,*),mastset,set(*),slavset
!
      integer ncont,ntie,i,j,k,nset,istartset(*),iendset(*),ialset(*),
     &  imast,nelem,jface,ncone,islav,ismallsliding,ipos,mortar,istep,
     &  kflag,idummy,jact
!
      real*8 tietol(3,*)
!
!     number of master triangles
!
      ncont=0
!
!     number of slave entities (nodes for nodal surfaces and
!     faces for facial surfaces)
!
      ncone=0
!
      do i=1,ntie
!
!        check for contact conditions
!
         if((tieset(1,i)(81:81).eq.kind1).or.
     &      (tieset(1,i)(81:81).eq.kind2)) then
            if(tietol(1,i).lt.0.d0) then
               ismallsliding=1
            else
               ismallsliding=0
            endif
            mastset=tieset(3,i)
!
!           determining the master surface
!
            do j=1,nset
               if(set(j).eq.mastset) exit
            enddo
            if(j.gt.nset) then
               ipos=index(mastset,' ')
               write(*,*) '*ERROR in allocont: master surface ',
     &               mastset(1:ipos-2)
               write(*,*) '       does not exist or does not contain'
               write(*,*) '       element faces'
               call exit(201)
            endif
            imast=j
!
!           deleting identical entries in the master facial surface
!           definition (leads otherwise to problems in the triangulation
!           and the creation of imastop)
!
            if(istep.eq.1) then
               kflag=1
               call isortii(ialset(istartset(imast)),idummy,
     &            iendset(imast)-istartset(imast)+1,kflag)
               jact=istartset(imast)
               do j=istartset(imast)+1,iendset(imast)
                  if(ialset(j).eq.ialset(j-1)) cycle
                  jact=jact+1
                  ialset(jact)=ialset(j)
               enddo
               iendset(imast)=jact
            endif
!
            do j=istartset(imast),iendset(imast)
!
               nelem=int(ialset(j)/10.d0)
               jface=ialset(j)-10*nelem
!
               if(lakon(nelem)(4:5).eq.'20') then
                  ncont=ncont+6
               elseif(lakon(nelem)(4:4).eq.'2') then
                  ncont=ncont+8
               elseif(lakon(nelem)(4:4).eq.'8') then
                  ncont=ncont+2
               elseif(lakon(nelem)(4:5).eq.'10') then
                  ncont=ncont+4
               elseif(lakon(nelem)(4:5).eq.'14') then
                  ncont=ncont+6
               elseif(lakon(nelem)(4:4).eq.'4') then
                  ncont=ncont+1
               elseif(lakon(nelem)(4:5).eq.'15') then
                  if(jface.le.2) then
                     ncont=ncont+4
                  else
                     ncont=ncont+6
                  endif
               elseif(lakon(nelem)(4:4).eq.'6') then
                  if(jface.le.2) then
                     ncont=ncont+1
                  else
                     ncont=ncont+2
                  endif
               endif
            enddo
!
!           counting the slave nodes
!
            slavset=tieset(2,i)
            ipos=index(slavset,' ')-1
            if(slavset(ipos:ipos).eq.'T') then
!
!              face-to-face penalty contact (facial slave surface)
!
               mortar=1
               nodeslavsurf=.false.
            elseif(slavset(ipos:ipos).eq.'M') then
!
!              quad-quad Mortar contact (facial slave surface)
!
               mortar=2
               nodeslavsurf=.false.
            elseif(slavset(ipos:ipos).eq.'P') then
!
!              quad-lin Petrov Galerkin Mortar contact (facial slave surface)
!
               mortar=4
               nodeslavsurf=.false.
            elseif(slavset(ipos:ipos).eq.'G') then
!
!              quad-quad Petrov Galerkin Mortar contact (facial slave surface)
!
               mortar=5
               nodeslavsurf=.false.
            elseif(slavset(ipos:ipos).eq.'O') then
!
!              quad-lin Mortar contact (facial slave surface)
!
               mortar=3
               nodeslavsurf=.false.
            else
!
!              node-to-face contact
!              default is a nodal slave surface
!
               nodeslavsurf=.true.
            endif
!
!           determining the slave surface
!
            do j=1,nset
               if(set(j).eq.slavset) exit
            enddo
            if(j.gt.nset) then
               if(mortar.eq.1) then
                  write(*,*)
     &              '*ERROR in allocont: element slave surface ',
     &              slavset(1:ipos-1)
                  write(*,*) '       does not exist'
                  call exit(201)
               endif
               do j=1,nset
                  if((set(j)(1:ipos-1).eq.slavset(1:ipos-1)).and.
     &                 (set(j)(ipos:ipos).eq.'T')) then
                     nodeslavsurf=.false.
                     exit
                  endif
               enddo
               if(j.gt.nset) then
                  write(*,*) '*ERROR in allocont: slave surface ',
     &                 slavset(1:ipos-1)
                  write(*,*) '       does not exist'
                  call exit(201)
               endif
            endif
!
            islav=j
!
!           deleting identical entries in the slave facial surface
!           definition (leads otherwise to problems in the calculation
!           of the are corresponding to the slave nodes)
!
            if((istep.eq.1).and.((mortar.eq.1).or.(.not.nodeslavsurf)))
     &         then
               kflag=1
               call isortii(ialset(istartset(islav)),idummy,
     &            iendset(islav)-istartset(islav)+1,kflag)
               jact=istartset(islav)
               do j=istartset(islav)+1,iendset(islav)
                  if(ialset(j).eq.ialset(j-1)) cycle
                  jact=jact+1
                  ialset(jact)=ialset(j)
               enddo
               iendset(islav)=jact
            endif
!
!           counting the entities (nodes or faces) in the slave
!           surface
!
            do j=istartset(islav),iendset(islav)
               if(ialset(j).gt.0) then
                  ncone=ncone+1
               else
                  k=ialset(j-2)
                  do
                     k=k-ialset(j)
                     if(k.ge.ialset(j-1)) exit
                     ncone=ncone+1
                  enddo
               endif
            enddo
!
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c
c
c =========================================================
      subroutine src1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt)
c =========================================================
      implicit none

      double precision q(1-mbc:maxmx+mbc, meqn)
      double precision aux(1-mbc:maxmx+mbc, *)

c     integrate a source term for friction if ifrictiontype>0.

      integer maxmx,meqn,mbc,mx,maux
      double precision xlower,dx,t,dt

      integer i
      double precision h,hu,hm,p,u,D,tau,rho,theta,tanpsi,prate,peq,dp2
      double precision m,hrate,g

      include "digparamsdec.i"
      include "digparamscommon.i"

      g=grav

      call calcaux(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux)

      do i=0,mx+1

         if (q(i,1).gt.dry_tol) then
            h = q(i,1)
            hu= q(i,2)
            hm= q(i,3)
            p = q(i,4)
            u = hu/h
            D = aux(i,i_D)
            tau = aux(i,i_tau)
            rho = aux(i,i_rho)
            theta = aux(i,i_theta)
            tanpsi = aux(i,i_tanpsi)
            prate = -kappita/(h*mu)*
     &            ((2.d0/(alpha*h))+rho_f*g*cos(theta))
            peq = rho_f*g*h*cos(theta)
            dp2 = -2.d0*u*tanpsi/(alpha*h)

            h  = h + dt*D
            h = max(h,0.d0)
            if (u.gt.0.d0) then
               hu = hu + dt*(g*h*sin(theta) - tau/rho)
            elseif (u.lt.0.d0) then
               hu = hu + dt*(g*h*sin(theta) + tau/rho)
            else
c               hu = hu + dt*max(g*h*sin(theta)-tau/rho,0.d0)
                hu = hu + dt*g*h*sin(theta)
            endif

            p = peq + (p - peq)*exp(prate*dt)
            p = p + dt*dp2

            q(i,1) = h
            q(i,2) = hu
            q(i,4) = max(p,0.d0)
            q(i,4) = min(q(i,4),rho*h*g*cos(theta))
         endif
      enddo

      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
C-------GF 28.08.1980------------------------- P R E S E M S
      INTEGER*2 FSEM,SEMHIB,PNOSEM
     = ,ACCO,ACIN,ACMA,ACTA  ! ACCUMULATOR
     = ,ATCO,ATIN,ATMA,SYAT  ! ATTRIBUTE
     = ,SOCO,REIN,UNCH,SOTA  ! SON, MEMBER
     = ,SYCO,SYIN,SYMA,SEPR  ! NEW SYMBOL
     = ,ACCU,ATTR,EOSY,EOS   ! MISCELLANEOUS
        INTEGER*2 SEMSYM(XXSEMH)
      INTEGER*2 SEM   (XXSEMH)
      COMMON /SEMS/
     =  FSEM      ! 1ST FREE 'SEM'-ENTRY
     = ,SEMHIB  ! HIGH BOUND OF 'SEM'
     = ,PNOSEM  ! -> DUMMY ACTION: UNCH, EOS
     = ,ACCO,ACIN,ACMA,ACTA  ! ACCUMULATOR
     = ,ATCO,ATIN,ATMA,SYAT  ! ATTRIBUTE
     = ,SOCO,REIN,UNCH,SOTA  ! SON, MEMBER
     = ,SYCO,SYIN,SYMA,SEPR  ! NEW SYMBOL
     = ,ACCU,ATTR,EOSY,EOS   ! MISCELLANEOUS
      EQUIVALENCE (PRE(1),SEM(1)),(PRESTA(1),SEMSYM(1))

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
!> @file
!! Include fortran file for gather operations, inner file
!! @author
!!    Copyright (C) 2012-2013 BigDFT group
!!    This file is distributed under the terms of the
!!    GNU General Public License, see ~/COPYING file
!!    or http://www.gnu.org/copyleft/gpl.txt .
!!    For the list of contributors, see ~/AUTHORS
  if (present(root)) then
     iroot=root
  else
     iroot=0
  end if
  if (present(comm)) then
     mpi_comm=comm
  else
     mpi_comm=MPI_COMM_WORLD !or bigdft_mpi%mpi_comm?
  end if

  !verify the size of the receive buffer
  ntasks=mpisize(mpi_comm)
  if (ntotrecv*kind(recvbuf) < ntot*ntasks*kind(sendbuf)) then
     call f_err_throw('Error in mpigather; the size of receive buffer ('//&
          trim(yaml_toa(ntotrecv*kind(recvbuf)))//&
          ') is not enough to contain '//trim(yaml_toa(ntot*kind(sendbuf)))//&
          ' * '//trim(yaml_toa(ntasks))//' elements',err_id=ERR_MPI_WRAPPERS)
     return
  end if
  !then one can proceed with the MPI operation
  ntotrecv=int(int(ntot,kind=8)*kind(sendbuf)/int(kind(recvbuf),kind=8))
  call f_timer_interrupt(TCAT_GATHER)
  call MPI_GATHER(sendbuf,ntot,mpitype(sendbuf),&
       recvbuf,ntotrecv,mpitype(recvbuf),iroot,mpi_comm,ierr)
  call f_timer_resume()
  if (ierr /=0) then
     call f_err_throw('An error in calling to MPI_GATHER occured',&
          err_id=ERR_MPI_WRAPPERS)
     return
  end if

max_stars_repo_name: gimunu/octopus-metric

max_issues_repo_name: gimunu/octopus-metric

max_forks_repo_name: gimunu/octopus-metric


module leq
  implicit none

  private

  public :: leq_init, leqin, gradient, eqitem, bgradient, leq_finish
  public :: invR, Rpos, Zpos, diameter, btori, dbtori,  qfun, pfun
  public :: dpfun, betafun, psi, rcenter, dpdrhofun

  integer :: nr, nt

  real, allocatable, dimension (:)     :: eqpsi, fp, beta, pressure
  real, allocatable, dimension (:,:)   :: R_psi, Z_psi
  real, allocatable, dimension (:,:,:) :: drm, dzm, dbtm, dpm, dtm
  real, allocatable, dimension (:,:,:) :: dpcart, dtcart, dbtcart
  real, allocatable, dimension (:,:,:) :: dpbish, dtbish, dbtbish

  real :: beta_0

  type :: flux_surface
     real :: R_center, R_geo, k, kp, d, dp, r, dr, delp, q, shat, pp, a, ap
     integer :: nt
  end type flux_surface

  type (flux_surface) :: surf

contains
  subroutine leqin(R0, Ra, k, kp, d, dp, r, dr, s, qq, qs, a, ap, nt_used)
    implicit none
    real, intent(in) :: R0, Ra, k, kp, d, dp, r, dr, s, qq, qs, a, ap
    integer, intent(in) :: nt_used

!R0   3.1601153788426086        3.0818707511680938        1.3431053933541011
!      0.10381317967192064        3.3317843820181277E-002   8.0267210436786265E-002
!  0.37777777777777777        1.0000000000000000E-003  -4.6851042901036040E-002
!   1.0642542534450539       0.32875578050587551        0.0000000000000000        0.0000000000000000                4
!f R0   3.1601153788426086        3.0818707511680938        1.3431053933541011
!        2.5126356633383937E-316   3.3317843820181277E-002   4.9406564584124654E-324
!   2.5126356633383937E-316   1.0000000000000000E-003   9.1922790858218116E-317
!   1.0642542534450539        2.3564531823208357E-310   0.0000000000000000        0.0000000000000000                4
!f R0   3.1601153788426086        3.0818707511680938        1.3431053933541011
!        2.2724233178454796E-316   3.3317843820181277E-002   4.9406564584124654E-324
!   2.2724233178454796E-316   1.0000000000000000E-003   9.1922790858218116E-317
!   1.0642542534450539        2.3203788153401778E-310   0.0000000000000000        0.0000000000000000                4

! kp, dp, r, s, qs
    !write (*,*) 'leq inputs', R0,Ra, k, kp, d, dp, r, dr, s, qq, qs, a, ap, nt_used
    surf%R_center = R0
    surf%R_geo = Ra
    surf%delp = s
    surf%k = k
    surf%kp = kp
    surf%q = qq
    surf%shat = qs
    surf%d = d
    surf%dp = dp
    surf%a = a
    surf%ap = ap
    surf%r = r
    surf%dr = dr
    surf%pp = 0.

    beta_0 = 1.

    nr = 3
    nt = nt_used
    if(.not.allocated(beta)) call alloc_arrays(3, nt)
    surf%nt = nt
    call leq_init
  end subroutine leqin

  subroutine alloc_arrays(nr, nt)
    integer, intent(in) :: nr, nt

    allocate(eqpsi(nr), fp(nr), beta(nr), pressure(nr))
    allocate(R_psi(nr, nt), Z_psi(nr, nt))
    allocate(drm(nr, nt, 2), dzm(nr, nt, 2), dbtm(nr, nt, 2), &
         dpm(nr, nt, 2), dtm(nr, nt, 2))
    allocate(dpcart(nr, nt, 2), dtcart(nr, nt, 2), dbtcart(nr, nt, 2))
    allocate(dpbish(nr, nt, 2), dtbish(nr, nt, 2), dbtbish(nr, nt, 2))
  end subroutine alloc_arrays

  subroutine dealloc_arrays
    implicit none
    if(allocated(eqpsi)) deallocate(eqpsi, fp, beta, pressure)
    if(allocated(R_psi)) deallocate(R_psi, Z_psi)
    if(allocated(drm)) deallocate(drm,dzm,dbtm,dpm,dtm)
    if(allocated(dpcart)) deallocate(dpcart,dtcart,dbtcart)
    if(allocated(dpbish)) deallocate(dpbish,dtbish,dbtbish)
  end subroutine dealloc_arrays

  subroutine leq_finish
    implicit none
    call dealloc_arrays
  end subroutine leq_finish

  subroutine leq_init
    use constants, only: pi=>pi
    implicit none
    real, dimension(nr, nt) :: eqpsi1, eqth, eqbtor
    real :: dr(3)
    real :: t, r
    integer :: i, j

    dr(1) = -surf%dr
    dr(2) = 0.
    dr(3) = surf%dr

    do j=1,nt
       do i=1,nr
          r = surf%r + dr(i)
          t = (j-1)*pi/real(nt-1)
          R_psi(i,j) = Rpos(r, t)
          Z_psi(i,j) = Zpos(r, t)
          eqth(i,j) = t
          eqpsi1(i,j) = 1 + dr(i)
          eqbtor(i,j) = surf%r_geo/R_psi(i,j)
       enddo
    enddo

    do i=1,nr
       pressure(i) = -dr(i)
    enddo

    eqpsi(:) = eqpsi1(:,1)

    call derm(eqth,   dtm,  'T')
    call derm(R_psi,  drm,  'E')
    call derm(Z_psi,  dzm,  'O')
    call derm(eqbtor, dbtm, 'E')
    call derm(eqpsi1, dpm,  'E')

! diagnostics
!      do j=1,nt
!         do i=1,nr
!            write(*,*) i,j
!            write(*,100) drm(i,j,1),drm(i,j,2),R_psi(i,j)
!            write(*,101) dzm(i,j,1),dzm(i,j,2),Z_psi(i,j)
!            write(*,102) dtm(i,j,1),dtm(i,j,2),eqth(i,j)
!         enddo
!      enddo
! 100  format('(gr R)1 ',g10.4,' (gr R)2 ',g10.4,' R ',g10.4)
! 101  format('(gr Z)1 ',g10.4,' (gr Z)2 ',g10.4,' Z ',g10.4)
! 102  format('(gr t)1 ',g10.4,' (gr t)2 ',g10.4,' t ',g10.4)
!      write(*,*) nr, nt
!      stop

! below is actually grad(rho) instead of grad(psi),
! and 'cartesian' refers to (R,Z) coordinates -- MAB
! grad(psi) in cartesian form
    call eqdcart(dpm, dpcart)
! grad(psi) in Bishop form
    call eqdbish(dpcart, dpbish)

! grad(BT) in cartesian form
    call eqdcart(dbtm, dbtcart)
! grad(BT) in Bishop form
    call eqdbish(dbtcart, dbtbish)

! grad(theta) in cartesian form
    call eqdcart(dtm, dtcart)
! grad(theta) in Bishop form
    call eqdbish(dtcart, dtbish)

  end subroutine leq_init

  subroutine derm(f, dfm, char)
    use constants, only: pi=>pi
    implicit none
    real, dimension(:,:), intent(in) :: f
    real, dimension(:,:,:), intent(out) :: dfm
    character(1), intent(in) :: char
    integer :: i, j

    i=1
    dfm(i,:,1) = -0.5*(3*f(i,:)-4*f(i+1,:)+f(i+2,:))

    i=nr
    dfm(i,:,1) = 0.5*(3*f(i,:)-4*f(i-1,:)+f(i-2,:))

! assume up-down symmetry for now:

    select case (char)
    case ('E')
       j=1
       dfm(:,j,2) = 0.5*(f(:,j+1)-f(:,j+1))

       j=nt
       dfm(:,j,2) = -0.5*(f(:,j-1)-f(:,j-1))
    case ('O')
       j=1
       dfm(:,j,2) = 0.5*(f(:,j+1)+f(:,j+1))

       j=nt
       dfm(:,j,2) = -0.5*(f(:,j-1)+f(:,j-1))
    case ('T')
       j=1
       dfm(:,j,2) = f(:,j+1)

       j=nt
       dfm(:,j,2) = pi - f(:,j-1)
    end select

    do i=2,nr-1
       dfm(i,:,1)=0.5*(f(i+1,:)-f(i-1,:))
    enddo

    do j=2,nt-1
       dfm(:,j,2)=0.5*(f(:,j+1)-f(:,j-1))
    enddo
  end subroutine derm

  subroutine gradient(theta, grad, char, rp, nth_used, ntm)
    use splines, only: inter_d_cspl
    implicit none
    integer, intent(in) :: nth_used, ntm
    character(1), intent(in) :: char
    real, dimension(-ntm:), intent(in) :: theta
    real, dimension(-ntm:,:), intent(out) :: grad
    real, intent(in) :: rp
    real :: tmp(2), aa(1), daa(1), rpt(1)
    real, dimension(nr,nt,2) :: dcart
    integer :: i

    select case(char)
    case('D')  ! diagnostic
       dcart = dbtcart
    case('P')
       dcart = dpcart
    case('R')
       dcart = dpcart  ! dpcart is correct for 'R'
    case('T')
       dcart = dtcart
    end select

    do i=-nth_used,-1
       call eqitem(theta(i), dcart(:,:,1), tmp(1), 'R')
       call eqitem(theta(i), dcart(:,:,2), tmp(2), 'Z')
       if(char == 'T') then
          grad(i,1)=-tmp(1)
          grad(i,2)=-tmp(2)
       else
          grad(i,1)=tmp(1)
          grad(i,2)=tmp(2)
       endif
    enddo

    do i=0,nth_used
       call eqitem(theta(i), dcart(:,:,1), tmp(1), 'R')
       call eqitem(theta(i), dcart(:,:,2), tmp(2), 'Z')
       grad(i,1)=tmp(1)
       grad(i,2)=tmp(2)
    enddo

!     to get grad(pressure), multiply grad(psi) by dpressure/dpsi

    if(char == 'R') then
       rpt(1) = rp
       call inter_d_cspl(nr, eqpsi, pressure, 1, rpt, aa, daa)
       do i=-nth_used, nth_used
          grad(i,1)=grad(i,1)*daa(1)*0.5*beta_0
          grad(i,2)=grad(i,2)*daa(1)*0.5*beta_0
       enddo
    endif
  end subroutine gradient

  subroutine bgradient(theta, grad, char, rp, nth_used, ntm)
    use splines, only: inter_d_cspl
    implicit none
    integer, intent(in) :: nth_used, ntm
    character(1), intent(in) :: char
    real, dimension(-ntm:), intent(in) :: theta
    real, dimension(-ntm:,:), intent(out) :: grad
    real, intent(in) :: rp
    real :: aa(1), daa(1), rpt(1)
    real, dimension(nr,nt,2) ::  dbish
    integer :: i

    select case(char)
    case('D')  ! diagnostic
       dbish = dbtbish
    case('P')
       dbish = dpbish
    case('R')
       dbish = dpbish  ! dpcart is correct for 'R'
    case('T')
       dbish = dtbish
    end select

    do i=-nth_used,nth_used
       call eqitem(theta(i), dbish(:,:,1), grad(i,1), 'R')
       call eqitem(theta(i), dbish(:,:,2), grad(i,2), 'Z')
    enddo

    if (char == 'T') then
       where (theta(-nth_used:nth_used) < 0.0)
          grad(-nth_used:nth_used,1) = -grad(-nth_used:nth_used,1)
          grad(-nth_used:nth_used,2) = -grad(-nth_used:nth_used,2)
       end where
    end if

!     to get grad(pressure), multiply grad(psi) by dpressure/dpsi

    if(char == 'R') then
       rpt(1) = rp
       call inter_d_cspl(nr, eqpsi, pressure, 1, rpt, aa, daa)
       do i=-nth_used, nth_used
          grad(i,1)=grad(i,1)*daa(1) * 0.5*beta_0
          grad(i,2)=grad(i,2)*daa(1) * 0.5*beta_0
       enddo
    endif
  end subroutine bgradient

  subroutine eqitem(theta_in, f, fstar, char)
    use constants, only: pi=>pi
    implicit none
    real, intent(in) :: theta_in
    real, dimension(:,:), intent(in) :: f
    real, intent(out) :: fstar
    character(1), intent(in) :: char
    integer :: j, istar, jstar
    real :: thet, tp, tps
    real :: st, dt
    real, dimension(size(f,2)) :: mtheta
! find r on psi mesh

    istar = 2

! Now do theta direction

    thet = mod2pi(theta_in)

! assume up-down symmetry

    tp=abs(thet)
    tps=1.
    if(char == 'Z' .and. thet /= 0.) tps=thet/abs(thet)

! get thet on theta mesh

    mtheta = (/ ( real(j-1)*pi/real(nt-1), j=1,nt) /)

! note that theta(1)=0 for local_eq theta

    jstar=-1
    do j=1,nt
       if(tp < mtheta(j)) then
          dt = tp - mtheta(j-1)
          st = mtheta(j) - tp
          jstar=j-1
          exit
       endif
       if(jstar /= -1) write(*,*) 'exit error j'
    enddo

! treat theta = pi separately

    if(jstar == -1) then
       jstar=nt-1
       dt=mtheta(jstar+1)-mtheta(jstar)
       st=0.
    endif

! use opposite area stencil to interpolate

    fstar=f(istar    , jstar    )  * st &
         +f(istar    , jstar + 1)  * dt
    fstar=fstar*tps &
         /(mtheta(jstar+1)-mtheta(jstar))
!     write(*,*) i, dr, dt, sr, st
!     write(*,*) f(istar,jstar+1),f(istar+1,jstar+1)
!     write(*,*) f(istar,jstar),f(istar+1,jstar)
!     write(*,*) eqpsi(istar),eqpsi(istar+1)
!     write(*,*) mtheta(jstar),mtheta(jstar+1)
  end subroutine eqitem

  subroutine eqdcart(dfm, dfcart)
    implicit none
    real, dimension (:,:,:), intent(in)  :: dfm
    real, dimension (:,:,:), intent(out) :: dfcart
    real, dimension (size(dfm,1),size(dfm,2)) :: denom
    integer :: i, j

    denom(:,:) = drm(:,:,1)*dzm(:,:,2) - drm(:,:,2)*dzm(:,:,1)

    dfcart = 0.

    dfcart(:,:,1) =   dfm(:,:,1)*dzm(:,:,2) - dzm(:,:,1)*dfm(:,:,2)
    dfcart(:,:,2) = - dfm(:,:,1)*drm(:,:,2) + drm(:,:,1)*dfm(:,:,2)

    do j=1,nt
       do i=2,nr
          dfcart(i,j,:)=dfcart(i,j,:)/denom(i,j)
       enddo
    enddo
  end subroutine eqdcart

  subroutine eqdbish(dcart, dbish)
    implicit none
    real, dimension(:, :, :), intent (in) :: dcart
    real, dimension(:, :, :), intent(out) :: dbish
    real, dimension(size(dcart,1),size(dcart,2)) :: denom
    integer :: i, j

    denom(:,:) = sqrt(dpcart(:,:,1)**2 + dpcart(:,:,2)**2)

    dbish(:,:,1) = dcart(:,:,1)*dpcart(:,:,1) + dcart(:,:,2)*dpcart(:,:,2)
    dbish(:,:,2) =-dcart(:,:,1)*dpcart(:,:,2) + dcart(:,:,2)*dpcart(:,:,1)

    do j=1,nt
       do i=2,nr
          dbish(i,j,:) = dbish(i,j,:)/denom(i,j)
      enddo
    enddo
  end subroutine eqdbish

  function invR (r, theta)
    real, intent (in) :: r, theta
    real :: invR

    invR=1./Rpos(r, theta)
  end function invR

  function rcenter()
    real :: rcenter

    rcenter = surf%R_center
  end function rcenter

  function Rpos (r, theta)
    use constants, only: pi
    implicit none
    real, intent (in) :: r, theta
    real :: Rpos
    real :: g, gp, dr

    dr = r - surf%r

! For Y Xiao:
!    g = surf%delp/surf%r + surf%d * sin(theta)**2
!    Rpos = surf%R_center*(1.+r*(cos(theta)-g)-g*dr)

    g = cos(theta + surf%d * sin(theta-0.5*pi*surf%a))
    gp = -sin(theta + surf%d * sin(theta-0.5*pi*surf%a))*(surf%dp*sin(theta-0.5*surf%a)-surf%d*0.5*pi*surf%ap*cos(theta-0.5*pi*surf%a))

    Rpos = surf%R_center + surf%delp*dr + g*surf%r + (g+surf%r*gp)*dr
  end function Rpos

  function Zpos (r, theta)
    real, intent (in) :: r, theta
    real :: Zpos, dr

    dr = r - surf%r
    Zpos = surf%k*sin(theta)*surf%r + (surf%r*surf%kp + surf%k)*sin(theta)*dr
  end function Zpos

  function psi (r)
    real, intent (in) :: r
    real :: psi

    psi = r - surf%r
  end function psi

  function mod2pi (theta)
    use constants, only: pi=>pi
    real, intent(in) :: theta
    real :: th, mod2pi
    real, parameter :: theta_tol = 1.e-6
    logical :: out
    if(theta <= pi .and. theta >= -pi) then
       mod2pi = theta
       return
    endif

    if(theta - theta_tol <= pi .and. theta >= -pi) then
       mod2pi = pi
       return
    endif

    if(theta <= pi .and. theta + theta_tol >= -pi) then
       mod2pi = -pi
       return
    endif

    th=theta
    out=.true.
    do while(out)
       if(th > pi) th = th - 2.*pi
       if(th <-pi) th = th + 2.*pi
       if(th <= pi .and. th >= -pi) out=.false.
    enddo
    mod2pi=th
  end function mod2pi

  function diameter (rp)
    real, intent(in) :: rp
    real :: diameter

    diameter = 2.*rp
  end function diameter

  function dbtori ()
    real :: dbtori
    dbtori = 1.
  end function dbtori

  function btori ()
    real :: btori
    btori = surf%r_geo
  end function btori

  function qfun ()
    real :: qfun
    qfun = surf%q
  end function qfun

  function pfun ()
    real :: pfun
    pfun = 0.5*beta_0
  end function pfun

  function dpfun ()
    real :: dpfun
    dpfun = -1.
  end function dpfun

  function dpdrhofun()
    real :: dpdrhofun
    dpdrhofun = surf%pp
  end function dpdrhofun

  function betafun ()
    real :: betafun
    betafun = beta_0
  end function betafun
end module leq


text:
	SUBROUTINE DSCLIP  ( iwndw, iret )
C************************************************************************
C* DSCLIP								*
C*									*
C* This subroutine saves the clipping bounds to be used in /DVWNDW/.	*
C*									*
C* DSCLIP  ( IWNDW, IRET )						*
C*									*
C* Input parameters:							*
C*	IWNDW		INTEGER		Coordinate system number	*
C*					  1 = D/N			*
C*					  2 = V				*
C*					  3 = P 			*
C*									*
C* Output parameters:							*
C*	IRET			INTEGER	Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	5/85	GEMPLT Version 3.1			*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DVWNDW.CMN'
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Set clipping bounds.
C
	IF  ( ( iwndw .ge. 1 ) .and. ( iwndw .le. 3 ) ) THEN
	    icleft = iwleft (iwndw)
	    icbot  = iwbot  (iwndw)
	    icrght = iwrght (iwndw)
	    ictop  = iwtop  (iwndw)
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak