text:
C     LAST UPDATE 16/03/89
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE PLOT3D
      IMPLICIT NONE
C
C Purpose: Multiple plotting of frames in 3D
C
C AUTHOR: Z.KAM and M.H.J.KOCH, EMBL - HAMBURG
C
      INCLUDE 'COMMON.FOR'
C
C Updates:
C 14/03/84 GRM Modified for ghost graphics
C 19/03/84 GRM Added hard copy facility
C
C Calls  11: FILL   , GETCHN , GETHDR , GETVAL , GRMODE , MINMAX
C            OPNFIL , OT3PLO , SAVPLO , TRMODE , PLOTON
C Ghost   7: ERASE  , FILON  , MAP    , MAPYL  , PLOTON , PICNOW
C            PSPACE
C Called by: OTOKO
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local variables:
C
      REAL    VALUE(10),XMIN,XMAX,FMIN,FMAX,RANGE,DELY,XTOT,YTOT
      INTEGER IRC,IFRAME,IHFMAX,IFRMAX,NCHAN,NFRAME,LENGTH,I,J,K
      INTEGER KREC,ISPEC,LSPEC,MEM,IFFR,ILFR,INCR,IFINC,NPTS,NM
      INTEGER ISAVE1,ISAVE2,ICH1,ICH2,IAUTO,IOPT,ILIN,IMAX,IMIN,NVAL
      CHARACTER*80 TITLE
      CHARACTER*13 HFNAM
      CHARACTER EOFCHAR
      LOGICAL IGETS
C
C IRC    : Return code
C IFRAME : Total nos. of frames
C IFRMAX : Nos. of frames in sequence
C IHFMAX : Nos. of header file in sequence
C KREC   : Output file record
C ISPEC  : First header file of sequence
C LSPEC  : Last header file in sequence
C MEM    : Positional or calibration data indicator
C IFFR   : First frame in sequence
C ILFR   : Last frame in sequence
C INCR   : Header file increment
C IFINC  : Frame increment
C HFNAM  : Header filename
C NCHAN  : Nos. of data points in spectrum
C NFRAME : Nos. of time frames
C ISAVE1 : Temporary storage
C ISAVE2 : Temporary storage
C ICH1   : First channel of interest
C ICH2   : Last channel of interest
C IAUTO  : Automatic scaling required
C IOPT   : Origin position on the screen
C ILIN   : Linear or logrithmic display
C IMIN   : Position of minimum
C IMAX   : Position of maximum
C XMIN   : Minimum value
C XMAX   : Maximum value
C FMIN   : Overall mimimum value of frames
C FMAX   : Overall maximum value of frames
C VALUE  : Numeric values entered at terminal
C NVAL   : Nos. of values entered at terminal
C RANGE  : Range of abscissa values
C DELY   : Ordinate increment
C NM     : ?
C XTOT   : Abscissa span
C YTOT   : Ordinate span
C NPTS   : Nos. of points in frame to plot
C TITLE  : Title of plot
C
C-----------------------------------------------------------------------
10    KREC=0
C
      CALL GETHDR (ITERM,IPRINT,IUNIT,HFNAM,ISPEC,LSPEC,INCR,MEM,IFFR,
     1             ILFR,IFINC,IHFMAX,IFRMAX,NCHAN,IRC)
      IF (IRC.NE.0) GOTO 999
C
      ISAVE1=ISPEC
      ISAVE2=IFFR
      IFRAME=IHFMAX+IFRMAX-1
C
      CALL GETCHN (ITERM,IPRINT,NCHAN,ICH1,ICH2,IRC)
      IF (IRC.NE.0) GOTO 10
      CALL FILL (SP3,NCHAN,0.0)
      IAUTO=0
20    IOPT=0
      WRITE (IPRINT,1000)
      CALL GETVAL (ITERM,VALUE,NVAL,IRC)
      IF (IRC.EQ.1) GOTO 10
      IF (IRC.EQ.2) GOTO 20
      IF (NVAL.GE.1) IOPT=INT(VALUE(1))
      IF (IOPT.NE.0.AND.IOPT.NE.1) GOTO 20
30    ILIN=0
      WRITE (IPRINT,1010)
      CALL GETVAL (ITERM,VALUE,NVAL,IRC)
      IF (IRC.EQ.1) GOTO 10
      IF (IRC.EQ.2) GOTO 30
      IF (NVAL.GE.1) ILIN=INT(VALUE(1))
      IF (ILIN.LT.0.OR.ILIN.GT.2) GOTO 30
      TITLE=' '
      WRITE (IPRINT,1015)
      IF (.NOT.IGETS (ITERM,TITLE)) GOTO 999
C
      DO 50 I=1,IHFMAX
C
         CALL OPNFIL (JUNIT,HFNAM,ISPEC,MEM,IFFR,ILFR,NCHAN,NFRAME,IRC)
         IF (IRC.EQ.0) THEN
C
            ISPEC=ISPEC+INCR
            DO 40 J=1,IFRMAX
               READ (JUNIT,REC=IFFR) (SP1(K),K=1,NCHAN)
               IFFR=IFFR+IFINC
               KREC=KREC+1
C
               IF (ILIN.EQ.2) THEN
                  DO 42 K=1,NCHAN
                     SP1(K)=SP1(K)*REAL(K)*REAL(K)
 42               CONTINUE
               ENDIF
               CALL MINMAX (SP1,NCHAN,ICH1,ICH2,IMIN,IMAX,XMIN,XMAX)
               IF (KREC.EQ.1) THEN
                  FMIN=XMIN
                  FMAX=XMAX
               ENDIF
               IF (XMIN.LT.FMIN) FMIN=XMIN
               IF (XMAX.GT.FMAX) FMAX=XMAX
40          CONTINUE
            CLOSE (UNIT=JUNIT)
         ENDIF
50    CONTINUE
C
      WRITE (IPRINT,1020) FMIN,FMAX
60    WRITE (IPRINT,1030) EOFCHAR ()
      CALL GETVAL (ITERM,VALUE,NVAL,IRC)
      IF (IRC.EQ.1) GOTO 70
      IF (IRC.EQ.2) GOTO 60
      IF (NVAL.GE.1) FMIN=VALUE(1)
      IF (NVAL.GE.2) FMAX=VALUE(2)
      IAUTO=1
70    NPTS=ICH2-ICH1+1
      IF (ILIN.EQ.1) THEN
         IF (FMIN.LT.0.001) FMIN=0.001
         FMIN=ALOG10(FMIN)*1000.0
         FMAX=ALOG10(FMAX)*1000.0
      ENDIF
      RANGE=FMAX-FMIN
      IF (IFRAME.EQ.1) THEN
         DELY=0.
         NM=2
      ELSE
         NM=(NPTS-1)/(IFRAME-1)
         IF (NM.LT.2) NM=2
         DELY=RANGE/(2.*(IFRAME-1))
      ENDIF
      XTOT=NPTS-1+(IFRAME-1)*(NM-1)
      YTOT=DELY*(IFRAME-1)+RANGE
C
      CALL GRMODE
      CALL PLOTON
      CALL PSPACE (0.25,1.1,0.175,0.9)
      CALL ERASE
C
C========PUT AXIS ANNOTATION & TITLE TO PLOT
C
      CALL CTRMAG (20)
      CALL MAP (0.0,1.0,0.0,1.0)
      CALL PCSCEN (0.5,-0.1,TITLE(1:LENGTH(TITLE)))
      CALL MAP   (0.0,XTOT,0.0,YTOT)
C
      KREC=0
      ISPEC=ISAVE1
      IFFR=ISAVE2
      CALL FILL (SPW,NCHAN,0.0)
      DO 110 I=1,IHFMAX
C
         CALL OPNFIL (JUNIT,HFNAM,ISPEC,MEM,IFFR,ILFR,NCHAN,NFRAME,IRC)
         IF (IRC.EQ.0) THEN
C
            ISPEC=ISPEC+INCR
            DO 100 J=1,IFRMAX
               READ (JUNIT,REC=IFFR) (SP1(K),K=1,NCHAN)
               IFFR=IFFR+IFINC
               KREC=KREC+1
C
               IF (ILIN.EQ.1) THEN
                  DO 80 K=ICH1,ICH2
                     IF (SP1(K).LT.0.001) SP1(K)=0.001
                     SP1(K)=ALOG10(SP1(K))*1000.0
80                CONTINUE
               ELSEIF (ILIN.EQ.2) THEN
                  DO 81 K=ICH1,ICH2
                     SP1(K)=SP1(K)*REAL(K)*REAL(K)
 81               CONTINUE
               ENDIF
               IF (IAUTO.NE.0) THEN
                  DO 90 K=ICH1,ICH2
                     IF (SP1(K).GT.FMAX) SP1(K)=FMAX
                     IF (SP1(K).LT.FMIN) SP1(K)=FMIN
90                CONTINUE
               ENDIF
               CALL OT3PLO (SP1(ICH1),SPW,NPTS,KREC,IOPT,ILIN,NM,
     1                      IFRAME,DELY,FMIN)
C
100         CONTINUE
            CLOSE (UNIT=JUNIT)
         ENDIF
110   CONTINUE
C
      CALL PICNOW
      CALL TRMODE
      CALL SAVPLO (ITERM,IPRINT)
      GOTO 10
999   RETURN
C
1000  FORMAT (' Origin on [0] RHS or [1] LHS of screen [0]: ',$)
1010  FORMAT (' [0] X,Y [1] X,Log(Y) [2] X,X^2*Y  [0]: ',$)
1015  FORMAT (' Enter plot title: ',$)
1020  FORMAT (' Minimum value is ',G13.5,' Maximum value is ',G13.5)
1030  FORMAT (' Enter new values or <ctrl-',a1,'>: ',$)
      END


SUBROUTINE D5SPLT (ALM,BLM,clm,MULT,UENORM,num,coord,jatom)
!
!     DECOMPOSITION OF D CHARGE IN DZ2,DX2-Y2,DXY,DXZ,DYZ
!     new (and reversed) definitions for dxz, dyz (6.2.89., Blaha)
!
  use param
  USE chard
  use lo
  IMPLICIT REAL*8 (A-H,O-Z)
  COMPLEX*16       ALM,BLM,clm,CSUMa,csumb,csumc
  DIMENSION      ALM((lmax2+1)*(lmax2+1)),BLM((lmax2+1)*(lmax2+1))
  DIMENSION        cLM((lomax+1)*(lomax+1))
  CHARACTER*5 COORD
!---------------------------------------------------------------------
!
  ipip=max(ilo(2),1)
  SQRT2=SQRT(2.D0)
  SQRT3=SQRT(3.D0)
  IF (COORD.EQ.'OCTAH')  THEN
     !
     ! non standard dsplit
     ! classical octahedral coordinates
     !
     IF (JATOM.EQ.1) THEN
        !     ===========================
        !     1. METALLATOM
        !
        !
        CSUMa=( SQRT3*(ALM(9)+ALM(5))/SQRT2 - ALM(7) )*0.5
        CSUMb=( SQRT3*(bLM(9)+bLM(5))/SQRT2 - bLM(7) )*0.5
        CSUMc=( SQRT3*(cLM(9)+cLM(5))/SQRT2 - cLM(7) )*0.5
        ADz2(num)=ADz2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDz2(num)=BDz2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdz2(num)=cdz2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadz2(num)=cadz2(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdz2(num)=acdz2(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdz2(num)=cbdz2(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdz2(num)=bcdz2(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
!
        CSUMa=(-1)*(ALM(8)+ALM(6))/SQRT2
        CSUMb=(-1)*(BLM(8)+BLM(6))/SQRT2
        CSUMc=(-1)*(cLM(8)+cLM(6))/SQRT2
        ADX2Y2(num)=ADX2Y2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDX2Y2(num)=BDX2Y2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdx2y2(num)=cdx2y2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadx2y2(num)=cadx2y2(num)+ &
             CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdx2y2(num)=acdx2y2(num)+ &
             CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdx2y2(num)=cbdx2y2(num)+ &
             CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdx2y2(num)=bcdx2y2(num)+ &
             CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
!
        CSUMa=( (ALM(9)+ALM(5))/SQRT2 + SQRT3*ALM(7) )*0.5
        CSUMb=( (BLM(9)+BLM(5))/SQRT2 + SQRT3*BLM(7) )*0.5
        CSUMc=( (cLM(9)+cLM(5))/SQRT2 + SQRT3*cLM(7) )*0.5
        ADXY(num)=ADXY(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDXY(num)=BDXY(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdxy(num)=cdxy(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadxy(num)=cadxy(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdxy(num)=acdxy(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdxy(num)=cbdxy(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdxy(num)=bcdxy(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
!
        CSUMa=((ALM(9)-ALM(5))/SQRT2-(ALM(8)-ALM(6))/SQRT2)/SQRT2
        CSUMb=((BLM(9)-BLM(5))/SQRT2-(BLM(8)-BLM(6))/SQRT2)/SQRT2
        CSUMc=((cLM(9)-cLM(5))/SQRT2-(cLM(8)-cLM(6))/SQRT2)/SQRT2
        ADXz(num)=ADXz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDXz(num)=BDXz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdxz(num)=cdxz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadxz(num)=cadxz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdxz(num)=acdxz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdxz(num)=cbdxz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdxz(num)=bcdxz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
        CSUMa=((-1.)*(ALM(9)-ALM(5))/SQRT2-(ALM(8)-ALM(6))/SQRT2)/SQRT2
        CSUMb=((-1.)*(BLM(9)-BLM(5))/SQRT2-(BLM(8)-BLM(6))/SQRT2)/SQRT2
        CSUMc=((-1.)*(cLM(9)-cLM(5))/SQRT2-(cLM(8)-cLM(6))/SQRT2)/SQRT2
        ADYz(num)=ADYz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDYz(num)=BDYz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdyz(num)=cdyz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadyz(num)=cadyz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdyz(num)=acdyz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdyz(num)=cbdyz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdyz(num)=bcdyz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !

     ELSE
        !     ===========================
        !     2.METALLATOM
        !
        CSUMa=(-SQRT3*(ALM(9)+ALM(5))/SQRT2-ALM(7) )*0.5
        CSUMb=(-SQRT3*(BLM(9)+BLM(5))/SQRT2-BLM(7) )*0.5
        CSUMc=(-SQRT3*(cLM(9)+cLM(5))/SQRT2-cLM(7) )*0.5
        ADz2(num)=ADz2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDz2(num)=BDz2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdz2(num)=cdz2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadz2(num)=cadz2(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdz2(num)=acdz2(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdz2(num)=cbdz2(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdz2(num)=bcdz2(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
        CSUMa=(-1.)*(ALM(8)-ALM(6))/SQRT2
        CSUMb=(-1.)*(BLM(8)-BLM(6))/SQRT2
        CSUMc=(-1.)*(cLM(8)-cLM(6))/SQRT2
        ADX2Y2(num)=ADX2Y2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDX2Y2(num)=BDX2Y2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdx2y2(num)=cdx2y2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadx2y2(num)=cadx2y2(num)+ &
             CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdx2y2(num)=acdx2y2(num)+ &
             CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdx2y2(num)=cbdx2y2(num)+ &
             CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdx2y2(num)=bcdx2y2(num)+ &
             CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
        CSUMa=( (-1.)*(ALM(9)+ALM(5))/SQRT2 + SQRT3*ALM(7) )*0.5
        CSUMb=( (-1.)*(BLM(9)+BLM(5))/SQRT2 + SQRT3*BLM(7) )*0.5
        CSUMc=( (-1.)*(cLM(9)+cLM(5))/SQRT2 + SQRT3*cLM(7) )*0.5
        ADXY(num)=ADXY(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDXY(num)=BDXY(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdxy(num)=cdxy(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadxy(num)=cadxy(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdxy(num)=acdxy(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdxy(num)=cbdxy(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdxy(num)=bcdxy(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
        CSUMa=((ALM(9)-ALM(5))/SQRT2-(ALM(8)+ALM(6))/SQRT2 )/SQRT2
        CSUMb=((BLM(9)-BLM(5))/SQRT2-(BLM(8)+BLM(6))/SQRT2 )/SQRT2
        CSUMc=((cLM(9)-cLM(5))/SQRT2-(cLM(8)+cLM(6))/SQRT2 )/SQRT2
        ADXz(num)=ADXz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDXz(num)=BDXz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdxz(num)=cdxz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadxz(num)=cadxz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdxz(num)=acdxz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdxz(num)=cbdxz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdxz(num)=bcdxz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
        CSUMa=((-1.)*(ALM(9)-ALM(5))/SQRT2-(ALM(8)+ALM(6))/SQRT2)/SQRT2
        CSUMb=((-1.)*(BLM(9)-BLM(5))/SQRT2-(BLM(8)+BLM(6))/SQRT2)/SQRT2
        CSUMc=((-1.)*(cLM(9)-cLM(5))/SQRT2-(cLM(8)+cLM(6))/SQRT2)/SQRT2
        ADYz(num)=ADYz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
        BDYz(num)=BDYz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
        cdyz(num)=cdyz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
        cadyz(num)=cadyz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
        acdyz(num)=acdyz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
        cbdyz(num)=cbdyz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
        bcdyz(num)=bcdyz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
        !
     END IF
     !
  ELSE if (COORD.EQ.'TRIGO')  THEN
     !
     CSUMa=(1/SQRT3)*((ALM(8)-ALM(6))+(ALM(9)+ALM(5))/SQRT2)
     CSUMb=(1/SQRT3)*((bLM(8)-bLM(6))+(bLM(9)+bLM(5))/SQRT2)
     CSUMc=(1/SQRT3)*((cLM(8)-cLM(6))+(cLM(9)+cLM(5))/SQRT2)
     ADz2(num)=ADz2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDz2(num)=BDz2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdz2(num)=cdz2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadz2(num)=cadz2(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdz2(num)=acdz2(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdz2(num)=cbdz2(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdz2(num)=bcdz2(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(1/SQRT3)*((aLM(8)+aLM(6))-(aLM(9)-aLM(5))/SQRT2)
     CSUMb=(1/SQRT3)*((bLM(8)+bLM(6))-(bLM(9)-bLM(5))/SQRT2)
     CSUMc=(1/SQRT3)*((cLM(8)+cLM(6))-(cLM(9)-cLM(5))/SQRT2)
     ADX2Y2(num)=ADX2Y2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDX2Y2(num)=BDX2Y2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdx2y2(num)=cdx2y2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadx2y2(num)=cadx2y2(num)+ &
          CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdx2y2(num)=acdx2y2(num)+ &
          CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdx2y2(num)=cbdx2y2(num)+ &
          CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdx2y2(num)=bcdx2y2(num)+ &
          CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=ALM(7)
     CSUMb=BLM(7)
     CSUMc=cLM(7)
     ADXY(num)=ADXY(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDXY(num)=BDXY(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdxy(num)=cdxy(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadxy(num)=cadxy(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdxy(num)=acdxy(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdxy(num)=cbdxy(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdxy(num)=bcdxy(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(1/SQRT3)*((ALM(8)+ALM(6))/SQRT2+(ALM(9)-ALM(5)))
     CSUMb=(1/SQRT3)*((bLM(8)+bLM(6))/SQRT2+(bLM(9)-bLM(5)))
     CSUMc=(1/SQRT3)*((cLM(8)+cLM(6))/SQRT2+(cLM(9)-cLM(5)))
     ADXz(num)=ADXz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDXz(num)=BDXz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdxz(num)=cdxz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadxz(num)=cadxz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdxz(num)=acdxz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdxz(num)=cbdxz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdxz(num)=bcdxz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(1/SQRT3)*((ALM(8)-ALM(6))/SQRT2-(ALM(9)+ALM(5)))
     CSUMb=(1/SQRT3)*((bLM(8)-bLM(6))/SQRT2-(bLM(9)+bLM(5)))
     CSUMc=(1/SQRT3)*((cLM(8)-cLM(6))/SQRT2-(cLM(9)+cLM(5)))
     ADYz(num)=ADYz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDYz(num)=BDYz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdyz(num)=cdyz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadyz(num)=cadyz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdyz(num)=acdyz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdyz(num)=cbdyz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdyz(num)=bcdyz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT

  ELSE

     !
     ! standard-split
     !
     ADZ2(num)=ADZ2(num)+ALM(7)*CONJG(ALM(7))*100.D0/MULT
     BDZ2(num)=BDZ2(num)+BLM(7)*CONJG(BLM(7))*UENORM*100.D0/MULT
     cdz2(num)=cdz2(num)+cLM(7)*CONJG(cLM(7))*100.D0/MULT
     cadz2(num)=cadz2(num)+Clm(7)*CONJG(alm(7))*pi12lo(ipip,2)*100.D0/MULT
     acdz2(num)=acdz2(num)+alm(7)*CONJG(Clm(7))*pi12lo(ipip,2)*100.D0/MULT
     cbdz2(num)=cbdz2(num)+Clm(7)*CONJG(blm(7))*pe12lo(ipip,2)*100.D0/MULT
     bcdz2(num)=bcdz2(num)+blm(7)*CONJG(Clm(7))*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(ALM(9)+ALM(5))/SQRT2
     CSUMb=(BLM(9)+BLM(5))/SQRT2
     CSUMc=(cLM(9)+cLM(5))/SQRT2
     ADX2Y2(num)=ADX2Y2(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDX2Y2(num)=BDX2Y2(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdx2y2(num)=cdx2y2(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadx2y2(num)=cadx2y2(num)+ &
          CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdx2y2(num)=acdx2y2(num)+ &
          CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdx2y2(num)=cbdx2y2(num)+ &
          CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdx2y2(num)=bcdx2y2(num)+ &
          CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(ALM(9)-ALM(5))/SQRT2
     CSUMb=(BLM(9)-BLM(5))/SQRT2
     CSUMc=(cLM(9)-cLM(5))/SQRT2
     ADXY(num)=ADXY(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDXY(num)=BDXY(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdxy(num)=cdxy(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadxy(num)=cadxy(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdxy(num)=acdxy(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdxy(num)=cbdxy(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdxy(num)=bcdxy(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(ALM(8)-ALM(6))/SQRT2
     CSUMb=(BLM(8)-BLM(6))/SQRT2
     CSUMc=(cLM(8)-cLM(6))/SQRT2
     ADXz(num)=ADXz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDXz(num)=BDXz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdxz(num)=cdxz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadxz(num)=cadxz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdxz(num)=acdxz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdxz(num)=cbdxz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdxz(num)=bcdxz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
     !
     CSUMa=(ALM(8)+ALM(6))/SQRT2
     CSUMb=(BLM(8)+BLM(6))/SQRT2
     CSUMc=(cLM(8)+cLM(6))/SQRT2
     ADYz(num)=ADYz(num)+CSUMa*CONJG(CSUMa)*100.D0/MULT
     BDYz(num)=BDYz(num)+CSUMb*CONJG(CSUMb)*UENORM*100.D0/MULT
     cdyz(num)=cdyz(num)+CSUMc*CONJG(CSUMc)*100.D0/MULT
     cadyz(num)=cadyz(num)+CSUMc*CONJG(CSUMa)*pi12lo(ipip,2)*100.D0/MULT
     acdyz(num)=acdyz(num)+CSUMa*CONJG(CSUMc)*pi12lo(ipip,2)*100.D0/MULT
     cbdyz(num)=cbdyz(num)+CSUMc*CONJG(CSUMb)*pe12lo(ipip,2)*100.D0/MULT
     bcdyz(num)=bcdyz(num)+CSUMb*CONJG(CSUMc)*pe12lo(ipip,2)*100.D0/MULT
  endif
  RETURN
END SUBROUTINE D5SPLT


text:
module mod_gridgen_blocks
#include <messenger.h>
    use mod_kinds,              only: rk, ik
    use mod_constants,          only: ZERO, ONE, TWO, THREE, FOUR, SIX, PI
    use mod_string,             only: string_t
    use mod_bc,                 only: create_bc
    use mod_plot3d_utilities,   only: get_block_points_plot3d, &
                                      get_block_elements_plot3d, &
                                      get_block_boundary_faces_plot3d

    use type_bc_state_group,    only: bc_state_group_t
    use type_bc_state,          only: bc_state_t
    use type_bc_state_wrapper,  only: bc_state_wrapper_t

    use mod_hdf_utilities
    use hdf5
    implicit none

contains

    !-------------------------------------------------------------------------------------
    !!
    !!
    !!  Create File: Grid + BC's
    !!  -----------------------------
    !!  create_mesh_file__singleblock
    !!  create_mesh_file__multiblock
    !!  create_mesh_file__D2E8M1
    !!
    !!
    !!  Generate grid: point arrays
    !!  --------------------------
    !!  meshgen_1x1x1_linear
    !!  meshgen_1x1x1_unit_linear
    !!  meshgen_2x2x2_linear
    !!  meshgen_2x2x1_linear
    !!  meshgen_3x3x3_linear
    !!  meshgen_3x3x3_unit_linear
    !!  meshgen_3x3x1_linear
    !!  meshgen_4x1x1_linear
    !!  meshgen_4x2x2_linear
    !!  meshgen_3x1x1_linear
    !!  meshgen_2x1x1_linear
    !!  meshgen_40x15x1_linear
    !!  meshgen_15x15x1_linear
    !!  meshgen_15x15x2_linear
    !!  meshgen_15x15x3_linear
    !!
    !!
    !!
    !**************************************************************************************

    !>  Write a ChiDG-formatted grid file consisting of:
    !!
    !!      - one block domain, D1
    !!      - Linear element mapping, M1
    !!      - boundary conditions initialized to Scalar Extrapolate.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   10/17/2016
    !!
    !!
    !!
    !---------------------------------------------------------------------------------------
    subroutine create_mesh_file__singleblock(filename,equation_sets,group_names,bc_state_groups,nelem_xi,nelem_eta,nelem_zeta,clusterx,x_max_in,x_min_in)
        character(*),               intent(in)              :: filename
        type(string_t),             intent(in), optional    :: equation_sets(:)
        type(string_t),             intent(in), optional    :: group_names(:,:)
        type(bc_state_group_t),     intent(in), optional    :: bc_state_groups(:)
        integer(ik),                intent(in)              :: nelem_xi
        integer(ik),                intent(in)              :: nelem_eta
        integer(ik),                intent(in)              :: nelem_zeta
        integer(ik),                intent(in), optional    :: clusterx
        real(rk),                   intent(in), optional    :: x_max_in
        real(rk),                   intent(in), optional    :: x_min_in

        character(:),                   allocatable :: user_msg
        class(bc_state_t),              allocatable :: bc_state
        character(len=10)                           :: patch_names(6)
        integer(HID_T)                              :: file_id, dom_id, patch_id, bcgroup_id
        integer(ik)                                 :: mapping, bcface, ierr, igroup, istate
        real(rk),                       allocatable :: nodes(:,:)
        integer(ik),                    allocatable :: elements(:,:)
        integer(ik),                    allocatable :: faces(:,:)
        real(rk),   dimension(:,:,:),   allocatable :: xcoords, ycoords, zcoords
        character(len=8)                            :: bc_face_strings(6)
        character(:),   allocatable                 :: bc_face_string

        ! Create/initialize file
        call initialize_file_hdf(filename)
        file_id = open_file_hdf(filename)

        ! Generate coordinates for first block
        call meshgen_NxNxN_linear(nelem_xi,nelem_eta,nelem_zeta,xcoords,ycoords,zcoords,clusterx,x_max_in,x_min_in)

        !
        ! Get nodes/elements
        !
        mapping = 1
        nodes    = get_block_points_plot3d(xcoords,ycoords,zcoords)
        elements = get_block_elements_plot3d(xcoords,ycoords,zcoords,mapping,idomain=1)

        !
        ! Add domains
        !
        if ( present(equation_sets) ) then
            call add_domain_hdf(file_id,'01',nodes,elements,'Cartesian',equation_sets(1)%get())
        else
            call add_domain_hdf(file_id,'01',nodes,elements,'Cartesian','Scalar Advection')
        end if

        !
        ! Set boundary conditions patch connectivities
        !
        dom_id = open_domain_hdf(file_id,'01')

        bc_face_strings = ["XI_MIN  ","XI_MAX  ","ETA_MIN ","ETA_MAX ","ZETA_MIN","ZETA_MAX"]
        do bcface = 1,size(patch_names)

            ! Get face node indices for boundary 'bcface'
            faces = get_block_boundary_faces_plot3d(xcoords,ycoords,zcoords,mapping,bcface)

            ! Set bc patch face indices
            bc_face_string  = trim(bc_face_strings(bcface))
            patch_id = create_patch_hdf(dom_id,bc_face_string)
            call set_patch_hdf(patch_id,faces)
            call close_patch_hdf(patch_id)

        end do !bcface

        !
        ! Create bc_state, 'Scalar Extrapolate'
        !
        call create_bc('Scalar Extrapolate', bc_state)

        !
        ! Add bc_group's
        !
        if (present(bc_state_groups)) then
            do igroup = 1,size(bc_state_groups)
                call create_bc_state_group_hdf(file_id,bc_state_groups(igroup)%name)

                bcgroup_id = open_bc_state_group_hdf(file_id,bc_state_groups(igroup)%name)

                do istate = 1,bc_state_groups(igroup)%nbc_states()
                    call add_bc_state_hdf(bcgroup_id, bc_state_groups(igroup)%bc_state(istate)%state)
                end do
                call close_bc_state_group_hdf(bcgroup_id)
            end do
        else
            call create_bc_state_group_hdf(file_id,'Default')

            bcgroup_id = open_bc_state_group_hdf(file_id,'Default')
            call add_bc_state_hdf(bcgroup_id,bc_state)
            call close_bc_state_group_hdf(bcgroup_id)

        end if

        !
        ! Set boundary condition groups for each patch
        !
        patch_names = ['XI_MIN  ','XI_MAX  ', 'ETA_MIN ', 'ETA_MAX ', 'ZETA_MIN', 'ZETA_MAX']
        do bcface = 1,size(patch_names)

            patch_id = open_patch_hdf(dom_id,trim(patch_names(bcface)))

            ! Set bc_group
            if (present(group_names)) then
                call set_patch_group_hdf(patch_id,group_names(1,bcface)%get())
            else
                call set_patch_group_hdf(patch_id,'Default')
            end if

            call close_patch_hdf(patch_id)

        end do

        ! Set 'Contains Grid'
        call set_contains_grid_hdf(file_id,'True')

        ! Close file
        call close_domain_hdf(dom_id)
        call close_file_hdf(file_id)
        call close_hdf()

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

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   10/19/2016
    !!
    !!
    !!
    !-------------------------------------------------------------------------------------
    subroutine create_mesh_file__multiblock(filename,equation_sets,group_names,bc_state_groups,     &
                                                          nelem_xi,  nelem_eta,  nelem_zeta,  &
                                                          xmax,ymax,zmax,clusterx)
        character(*),           intent(in)              :: filename
        type(string_t),         intent(in), optional    :: equation_sets(:)
        type(string_t),         intent(in), optional    :: group_names(:,:)
        type(bc_state_group_t), intent(in), optional    :: bc_state_groups(:)
        integer(ik),            intent(in)              :: nelem_xi,  nelem_eta,  nelem_zeta
        real(rk),               intent(in), optional    :: xmax, ymax, zmax
        integer(ik),            intent(in), optional    :: clusterx

        class(bc_state_t),  allocatable                 :: bc_state
        character(8)                                    :: faces(6)
        integer(HID_T)                                  :: file_id, dom1_id, dom2_id, bcface1_id, bcface2_id, &
                                                           bcgroup_id, patch1_id, patch2_id
        integer(ik)                                     :: mapping, bcface, ierr, igroup, istate, &
                                                           nxi_max, neta_max, nzeta_max,xi_mid
        real(rk),       allocatable                     :: nodes1(:,:), nodes2(:,:)
        integer(ik),    allocatable                     :: elements1(:,:), elements2(:,:)
        integer(ik),    allocatable                     :: faces1(:,:), faces2(:,:)
        real(rk),       allocatable, dimension(:,:,:)   :: xcoords1, ycoords1, zcoords1, &
                                                           xcoords2, ycoords2, zcoords2
        real(rk)                                        :: xmax_block1,xmin_block2, xmax_current,       &
                                                           xmax_new, ymax_new, zmax_new, ymax_current,  &
                                                           zmax_current
        character(len=8)                                :: bc_face_strings(6)
        character(:),   allocatable                     :: bc_face_string

        !
        ! Create/initialize file
        !
        call initialize_file_hdf(filename)
        file_id = open_file_hdf(filename)

        !
        ! Generate coordinates
        !
        call meshgen_NxNxN_linear(nelem_xi, nelem_eta, nelem_zeta, xcoords1,ycoords1,zcoords1,clusterx)

        !
        ! Split the block:
        !   - Take the second half of the block, store to block 2
        !   - Take the first half of the block, reset to block 1
        !
        nxi_max   = size(xcoords1,1)
        neta_max  = size(xcoords1,2)
        nzeta_max = size(xcoords1,3)
        xi_mid    = int(nint(real(nxi_max)/2.))

        xcoords2 = xcoords1(xi_mid:nxi_max,1:neta_max,1:nzeta_max)
        xcoords1 = xcoords1(1:xi_mid,      1:neta_max,1:nzeta_max)

        ycoords2 = ycoords1(xi_mid:nxi_max,1:neta_max,1:nzeta_max)
        ycoords1 = ycoords1(1:xi_mid,      1:neta_max,1:nzeta_max)

        zcoords2 = zcoords1(xi_mid:nxi_max,1:neta_max,1:nzeta_max)
        zcoords1 = zcoords1(1:xi_mid,      1:neta_max,1:nzeta_max)

!        !
!        ! Translate block2 to end of block1, in x-direction
!        !
!        xmax_block1 = maxval(xcoords1)
!        xmin_block2 = minval(xcoords2)
!        xcoords2 = xcoords2 + (xmax_block1-xmin_block2)

        !
        ! Scale block coordinates to have xmax,ymax,zmax
        !   - normalize current maximum to 1
        !   - then multiply by desired extent
        !   - default extents are (1,1,1)
        !
        xmax_new = 1._rk
        ymax_new = 1._rk
        zmax_new = 1._rk
        if (present(xmax)) xmax_new = xmax
        if (present(ymax)) ymax_new = ymax
        if (present(zmax)) zmax_new = zmax

        xmax_current = maxval(xcoords2)
        xcoords1 = xcoords1/xmax_current
        xcoords2 = xcoords2/xmax_current

        xcoords1 = xmax_new*xcoords1
        xcoords2 = xmax_new*xcoords2

        ymax_current = maxval(ycoords2)
        ycoords1 = ycoords1/ymax_current
        ycoords2 = ycoords2/ymax_current

        ycoords1 = ymax_new*ycoords1
        ycoords2 = ymax_new*ycoords2

        zmax_current = maxval(zcoords2)
        zcoords1 = zcoords1/zmax_current
        zcoords2 = zcoords2/zmax_current

        zcoords1 = zmax_new*zcoords1
        zcoords2 = zmax_new*zcoords2

        !
        ! Get nodes/elements
        !
        mapping = 1
        nodes1    = get_block_points_plot3d(xcoords1,ycoords1,zcoords1)
        nodes2    = get_block_points_plot3d(xcoords2,ycoords2,zcoords2)
        elements1 = get_block_elements_plot3d(xcoords1,ycoords1,zcoords1,mapping,idomain=1)
        elements2 = get_block_elements_plot3d(xcoords2,ycoords2,zcoords2,mapping,idomain=2)

        !
        ! Add domains
        !
        if ( present(equation_sets) ) then
            call add_domain_hdf(file_id,'01',nodes1,elements1,'Cartesian',equation_sets(1)%get())
            call add_domain_hdf(file_id,'02',nodes2,elements2,'Cartesian',equation_sets(2)%get())
        else
            call add_domain_hdf(file_id,'01',nodes1,elements1,'Cartesian','Scalar Advection')
            call add_domain_hdf(file_id,'02',nodes2,elements2,'Cartesian','Scalar Advection')
        end if

        !
        ! Set boundary conditions patch connectivities
        !
        dom1_id = open_domain_hdf(file_id,'01')
        dom2_id = open_domain_hdf(file_id,'02')

        bc_face_strings = ["XI_MIN  ","XI_MAX  ","ETA_MIN ","ETA_MAX ","ZETA_MIN","ZETA_MAX"]
        do bcface = 1,6
            ! Get face node indices for boundary 'bcface'
            faces1 = get_block_boundary_faces_plot3d(xcoords1,ycoords1,zcoords1,mapping,bcface)
            faces2 = get_block_boundary_faces_plot3d(xcoords2,ycoords2,zcoords2,mapping,bcface)

            ! Set bc patch face indices
            bc_face_string  = trim(bc_face_strings(bcface))
            patch1_id = create_patch_hdf(dom1_id,bc_face_string)
            patch2_id = create_patch_hdf(dom2_id,bc_face_string)

            call set_patch_hdf(patch1_id,faces1)
            call set_patch_hdf(patch2_id,faces2)

            call close_patch_hdf(patch1_id)
            call close_patch_hdf(patch2_id)

        end do !bcface

        !
        ! Create bc_state, 'Scalar Extrapolate'
        !
        call create_bc('Scalar Extrapolate', bc_state)

        !
        ! Add bc_group's
        !
        if (present(bc_state_groups)) then
            do igroup = 1,size(bc_state_groups)
                call create_bc_state_group_hdf(file_id,bc_state_groups(igroup)%name)

                bcgroup_id = open_bc_state_group_hdf(file_id,bc_state_groups(igroup)%name)

                do istate = 1,bc_state_groups(igroup)%nbc_states()
                    call add_bc_state_hdf(bcgroup_id, bc_state_groups(igroup)%bc_state(istate)%state)
                end do
                call close_bc_state_group_hdf(bcgroup_id)
            end do
        else
            call create_bc_state_group_hdf(file_id,'Default')

            bcgroup_id = open_bc_state_group_hdf(file_id,'Default')
            call add_bc_state_hdf(bcgroup_id,bc_state)
            call close_bc_state_group_hdf(bcgroup_id)

        end if

        !
        ! Assign groups to boundary condition patches
        !
        faces = ['XI_MIN  ','XI_MAX  ', 'ETA_MIN ', 'ETA_MAX ', 'ZETA_MIN', 'ZETA_MAX']
        do bcface = 1,size(faces)
            patch1_id = open_patch_hdf(dom1_id,trim(adjustl(faces(bcface))))
            patch2_id = open_patch_hdf(dom2_id,trim(adjustl(faces(bcface))))

            if (present(group_names)) then
                call set_patch_group_hdf(patch1_id,group_names(1,bcface)%get())
                call set_patch_group_hdf(patch2_id,group_names(2,bcface)%get())
            else
                call set_patch_group_hdf(patch1_id,'Default')
                call set_patch_group_hdf(patch2_id,'Default')
            end if

            call close_patch_hdf(patch1_id)
            call close_patch_hdf(patch2_id)
        end do

        ! Set 'Contains Grid'
        call set_contains_grid_hdf(file_id,'True')

        ! Close file
        call close_domain_hdf(dom1_id)
        call close_domain_hdf(dom2_id)
        call close_file_hdf(file_id)
        call close_hdf()

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

    !>  Write a ChiDG-formatted grid file consisting of:
    !!
    !!      - two block domains
    !!      - overlapping slightly
    !!      - boundary conditions initialized to Scalar Extrapolate. Interior boundaries
    !!        not set so they are detected as Chimera.
    !!
    !!  Incoming Parameter, 'matching' specifies if the elements should overlap with
    !!  a single, or potentially multiple elements
    !!
    !!     Block 1           Block 2 : matching=True        Block 2 : matching=False
    !!  .-----.-----.             .-----.-----.                  .-----.-----.
    !!  |     |     |             |     |     |                  |     |     |
    !!  |     |     |             |     |     |                  |     |     |
    !!  .-----.-----.             .-----.-----.                  |     |     |
    !!  |     |     |             |     |     |                  .-----.-----.
    !!  |     |     |             |     |     |                  |     |     |
    !!  .-----.-----.             .-----.-----.                  .-----.-----.
    !!
    !!  Abutting/Matching combinations:
    !!  -------------------------------
    !!
    !!       abutting = .true.        abutting = .false.
    !!       matching = .true.        matching = .true.
    !!          ----..----               ----.-.----
    !!              ||                       | |
    !!              ||                       | |
    !!          ----..----               ----.-.----
    !!              ||                       | |
    !!              ||                       | |
    !!          ----..----               ----.-.----
    !!
    !!
    !!       abutting = .true.         abutting = .false.
    !!       matching = .false.        matching = .false.
    !!
    !!         -----..-----              ----.-.----
    !!              ||                       | |
    !!              ||                       | |
    !!         -----.|                   ----|-.
    !!              |.-----                  .-|----
    !!              ||                       | |
    !!         -----..-----              ----.-.----
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   10/17/2016
    !!
    !!
    !!
    !---------------------------------------------------------------------------------------
    subroutine create_mesh_file__D2E8M1(filename,abutting,matching,equation_sets,group_names,bc_state_groups)
        character(*),           intent(in)              :: filename
        logical,                intent(in)              :: abutting
        logical,                intent(in)              :: matching
        type(string_t),         intent(in), optional    :: equation_sets(:)
        type(string_t),         intent(in), optional    :: group_names(:,:)
        type(bc_state_group_t), intent(in), optional    :: bc_state_groups(:)

        class(bc_state_t),  allocatable                 :: bc_state
        character(8)                                    :: faces(5)
        character(:),   allocatable                     :: user_msg
        integer(HID_T)                                  :: file_id, dom1_id, dom2_id, bcgroup_id, &
                                                           patch1_id, patch2_id
        integer(ik)                                     :: mapping, bcface, ierr, igroup, istate
        real(rk),       allocatable                     :: nodes1(:,:), nodes2(:,:)
        integer(ik),    allocatable                     :: elements1(:,:), elements2(:,:)
        integer(ik),    allocatable                     :: faces1(:,:), faces2(:,:)
        real(rk),       allocatable, dimension(:,:,:)   :: xcoords1, xcoords2, ycoords1, ycoords2, zcoords1, zcoords2
        real(rk)                                        :: xmax,ymax, xmax_current, ymax_current, zmax_current
        character(len=8)                                :: bc_face_strings(6)
        character(:),   allocatable                     :: bc_face_string

        ! Create/initialize file
        call initialize_file_hdf(filename)
        file_id = open_file_hdf(filename)

        ! Generate coordinates for first block
        call meshgen_2x2x2_linear(xcoords1,ycoords1,zcoords1)

        !
        ! Create second block by copying and translating first block.
        !
        xmax = maxval(xcoords1)
        ymax = maxval(ycoords1)

        !
        ! If abutting=.true., Create block2 by copying block1 and translating
        ! it by xmax of block1.
        !
        ! If abutting=.false., only translate by a fraction of xmax so there is overlap
        !
        if (abutting) then
            xcoords2 = xcoords1 + xmax
        else
            xcoords2 = xcoords1 + 0.90_rk*xmax
        end if
        ycoords2 = ycoords1

        !
        ! If matching=false, shift center plane of points so that the overlapping
        ! faces between blocks to not match exactly, and might have multiple Chimera donors
        !
        if (.not. matching) then
            ycoords2(:,2,:) = ycoords2(:,2,:) - 0.2*ymax
        end if

        !
        ! Set zcoords
        !
        zcoords2 = zcoords1

        !
        ! Scale coordinates to xmax, ymax, zmax
        !
        xmax_current = maxval(xcoords2)
        ymax_current = maxval(ycoords2)
        zmax_current = maxval(zcoords2)

        xcoords1 = (1.0_rk)*xcoords1/xmax_current
        ycoords1 = (1.0_rk)*ycoords1/ymax_current
        zcoords1 = (1.0_rk)*zcoords1/zmax_current

        xcoords2 = (1.0_rk)*xcoords2/xmax_current
        ycoords2 = (1.0_rk)*ycoords2/ymax_current
        zcoords2 = (1.0_rk)*zcoords2/zmax_current

        !
        ! Get nodes/elements
        !
        mapping = 1
        nodes1    = get_block_points_plot3d(xcoords1,ycoords1,zcoords1)
        nodes2    = get_block_points_plot3d(xcoords2,ycoords2,zcoords2)
        elements1 = get_block_elements_plot3d(xcoords1,ycoords1,zcoords1,mapping,idomain=1)
        elements2 = get_block_elements_plot3d(xcoords2,ycoords2,zcoords2,mapping,idomain=2)

        !
        ! Add domains
        !
        call add_domain_hdf(file_id,'01',nodes1,elements1,'Cartesian','Scalar Advection')
        call add_domain_hdf(file_id,'02',nodes2,elements2,'Cartesian','Scalar Advection')

        !
        ! Set boundary conditions patch connectivities
        !
        dom1_id = open_domain_hdf(file_id,'01')
        dom2_id = open_domain_hdf(file_id,'02')

        bc_face_strings = ["XI_MIN  ","XI_MAX  ","ETA_MIN ","ETA_MAX ","ZETA_MIN","ZETA_MAX"]
        do bcface = 1,6
            ! Get face node indices for boundary 'bcface'
            faces1 = get_block_boundary_faces_plot3d(xcoords1,ycoords1,zcoords1,mapping,bcface)
            faces2 = get_block_boundary_faces_plot3d(xcoords2,ycoords2,zcoords2,mapping,bcface)

            ! Set bc patch face indices
            bc_face_string  = trim(bc_face_strings(bcface))
            patch1_id = create_patch_hdf(dom1_id,bc_face_string)
            patch2_id = create_patch_hdf(dom2_id,bc_face_string)

            call set_patch_hdf(patch1_id,faces1)
            call set_patch_hdf(patch2_id,faces2)

            call close_patch_hdf(patch1_id)
            call close_patch_hdf(patch2_id)

        end do !bcface

        !
        ! Create bc_state, 'Scalar Extrapolate'
        !
        call create_bc('Scalar Extrapolate', bc_state)

        !
        ! Add bc_group's
        !
        if (present(bc_state_groups)) then
            user_msg = 'create_mesh_file__D2E8M1: Not quite ready to accept custom bc_group sets.'
            call chidg_signal(FATAL,user_msg)
        else
            call create_bc_state_group_hdf(file_id,'Default')

            bcgroup_id = open_bc_state_group_hdf(file_id,'Default')
            call add_bc_state_hdf(bcgroup_id,bc_state)
            call close_bc_state_group_hdf(bcgroup_id)

        end if

        !
        ! Set boundary condition states for Domain 1: Leave XI_MAX empty
        !
        faces = ['XI_MIN  ', 'ETA_MIN ', 'ETA_MAX ', 'ZETA_MIN', 'ZETA_MAX']
        do bcface = 1,size(faces)
            patch1_id = open_patch_hdf(dom1_id,trim(adjustl(faces(bcface))))
            call set_patch_group_hdf(patch1_id,'Default')
            call close_patch_hdf(patch1_id)
        end do

        !
        ! Set boundary condition states for Domain 2: Leave XI_MIN empty
        !
        faces = ['XI_MAX  ', 'ETA_MIN ', 'ETA_MAX ', 'ZETA_MIN', 'ZETA_MAX']
        do bcface = 1,size(faces)
            patch2_id = open_patch_hdf(dom2_id,trim(adjustl(faces(bcface))))
            call set_patch_group_hdf(patch2_id,'Default')
            call close_patch_hdf(patch2_id)
        end do

        ! Set 'Contains Grid'
        call set_contains_grid_hdf(file_id,'True')

        ! Close file
        call close_domain_hdf(dom1_id)
        call close_domain_hdf(dom2_id)
        call close_file_hdf(file_id)
        call close_hdf()

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

    !>  Generate a set of points defining a:
    !!      - nelem_xi by nelem_eta by nelem_zeta element, single-block, mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @date   10/20/2016
    !!
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !!
    !--------------------------------------------------------------------------------------
    subroutine meshgen_NxNxN_linear(nelem_xi,nelem_eta,nelem_zeta,xcoords,ycoords,zcoords,clusterx,x_max_in,x_min_in)
        integer(ik),    intent(in)                  :: nelem_xi
        integer(ik),    intent(in)                  :: nelem_eta
        integer(ik),    intent(in)                  :: nelem_zeta
        real(rk),       intent(inout),  allocatable :: xcoords(:,:,:)
        real(rk),       intent(inout),  allocatable :: ycoords(:,:,:)
        real(rk),       intent(inout),  allocatable :: zcoords(:,:,:)
        integer(ik),    intent(in),     optional    :: clusterx
        real(rk),       intent(in),     optional    :: x_max_in
        real(rk),       intent(in),     optional    :: x_min_in

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ierr, &
                       npts_xi, npts_eta, npts_zeta
        real(rk)    :: x,y,z, x_max, x_min

        npts_xi   = nelem_xi   + 1
        npts_eta  = nelem_eta  + 1
        npts_zeta = nelem_zeta + 1

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Set max/min X-Coordinate
        !
        if (present(x_max_in)) then
            x_max = x_max_in
        else
            x_max = ONE
        end if

        if (present(x_min_in)) then
            x_min = x_min_in
        else
            x_min = ZERO
        end if

        !
        ! Generate points
        !
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi

                    if (present(clusterx)) then
                        if ( clusterx == -1 ) then
                            !x = ONE - tanh( (PI/TWO)*(ONE - real(ipt_xi-1,rk)/real(npts_xi-1,rk) ) )/tanh(PI/TWO)
                            x = x_max - (x_max-x_min)*tanh( (PI/TWO)*(ONE - real(ipt_xi-1,rk)/real(npts_xi-1,rk) ) )/tanh(PI/TWO)
                        else if ( clusterx == 1 ) then
                            call chidg_signal(FATAL,"meshgen_NxNxN_linear: 'clusterx'=1 not yet implemented.")
                        else
                            call chidg_signal(FATAL,"meshgen_NxNxN_linear: Invalid value for 'clusterx'. -1,1.")
                        end if
                    else
                        !x = real(ipt_xi-1,rk)/real(npts_xi-1,rk)
                        x = x_min + (x_max - x_min)*real(ipt_xi-1,rk)/real(npts_xi-1,rk)
                    end if

                    if (ipt_xi == npts_xi) then
                        !x = ONE
                        x = x_max
                    end if

                    y = real(ipt_eta-1,rk)/real(npts_eta-1,rk)
                    z = real(ipt_zeta-1,rk)/real(npts_zeta-1,rk)

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                end do
            end do
        end do

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

    !> Generate a set of points defining a 1x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   5/24/2016
    !!
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !--------------------------------------------------------------------------------------
    subroutine meshgen_1x1x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ierr, npts_xi, &
                       npts_eta, npts_zeta
        real(rk)    :: x,y,z, dx, dy, dz

        ! elements (1x1x1) - linear

        npts_xi   = 2
        npts_eta  = 2
        npts_zeta = 2

        dx = 1._rk
        dy = 1._rk
        dz = 1._rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = ZERO
        do ipt_zeta = 1,npts_zeta
            y = ZERO
            do ipt_eta = 1,npts_eta
                x = ZERO
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

    !> Generate a set of points defining a 1x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !--------------------------------------------------------------------------------------
    subroutine meshgen_1x1x1_unit_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ierr, npts_xi, &
                       npts_eta, npts_zeta
        real(rk)    :: x,y,z, dx, dy, dz

        ! elements (1x1x1) - linear
        npts_xi   = 2
        npts_eta  = 2
        npts_zeta = 2

        dx = 2._rk
        dy = 2._rk
        dz = 2._rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = -ONE
        do ipt_zeta = 1,npts_zeta
            y = -ONE
            do ipt_eta = 1,npts_eta
                x = -ONE
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

    !> Generate a set of points defining a 2x2x2 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !--------------------------------------------------------------------------------------
    subroutine meshgen_2x2x2_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 27
        real(rk)                    :: x,y,z

        ! elements (2x2x2) - linear
        !
        !          *-------*-------*
        !         /       /       /|
        !        *-------*-------* |
        !       /       /       /| *
        !      *-------*-------* |/|
        !      |       |       | * |
        !      |       |       |/| *
        !      *-------*-------* |/
        !      |       |       | *
        !      |       |       |/
        !      *-------*-------*
        !
        !
        npts_xi   = 3
        npts_eta  = 3
        npts_zeta = 3

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    x = ONE*real(ipt_xi  -1,rk)/real(npts_xi  -1,rk)
                    y = ONE*real(ipt_eta -1,rk)/real(npts_eta -1,rk)
                    z = ONE*real(ipt_zeta-1,rk)/real(npts_zeta-1,rk)

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 2x2x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !--------------------------------------------------------------------------------------
    subroutine meshgen_2x2x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 18
        real(rk), dimension(npt)    :: x,y,z

        ! elements (2x2x1) - linear
        !
        !          *-------*
        !         /       /|
        !        *-------* |
        !       /       /| *
        !      *-------* |/|
        !      |       | * |
        !      |       |/| *
        !      *-------* |/
        !      |       | *
        !      |       |/
        !      *-------*
        !
        !
        x = [ZERO, ONE, TWO, ZERO, ONE, TWO, ZERO, ONE, TWO, &
             ZERO, ONE, TWO, ZERO, ONE, TWO, ZERO, ONE, TWO]

        y = [ZERO, ZERO, ZERO, ONE, ONE, ONE, TWO, TWO, TWO, &
             ZERO, ZERO, ZERO, ONE, ONE, ONE, TWO, TWO, TWO]

        z = [ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE]

        npts_xi   = 3
        npts_eta  = 3
        npts_zeta = 2

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 3x3x3 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_3x3x3_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 64
        !real(rk), dimension(npt)    :: x,y,z
        real(rk)                    :: x,y,z

        ! elements (3x3x3) - linear
        !
        !            *-------*-------*-------*
        !           /       /       /       /|
        !          *-------*-------*-------* |
        !         /       /       /       /| *
        !        *-------*-------*-------* |/|
        !       /       /       /       /| * |
        !      *-------*-------*-------* |/| *
        !      |       |       |       | * |/|
        !      |       |       |       |/| * |
        !      *-------*-------*-------* |/| *
        !      |       |       |       | * |/
        !      |       |       |       |/| *
        !      *-------*-------*-------* |/
        !      |       |       |       | *
        !      |       |       |       |/
        !      *-------*-------*-------*
        !
        !
!        x = [ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
!             ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
!             ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
!             ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
!             ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
!             ZERO, ONE, TWO, THREE]
!
!        y = [ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, &
!             THREE, THREE, THREE, THREE, &
!             ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, &
!             THREE, THREE, THREE, THREE, &
!             ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, &
!             THREE, THREE, THREE, THREE, &
!             ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, &
!             THREE, THREE, THREE, THREE]
!
!        z = [ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, &
!             ZERO, ZERO, ZERO, ZERO, ZERO, &
!             ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, &
!             ONE, ONE, ONE, &
!             TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, &
!             TWO, TWO, TWO, &
!             THREE, THREE, THREE, THREE, THREE, THREE, THREE, THREE, THREE, THREE, &
!             THREE, THREE, THREE, THREE, THREE, THREE]

        npts_xi   = 4
        npts_eta  = 4
        npts_zeta = 4

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    x = (real(ipt_xi  -1,rk))/real(npts_xi  -1,rk)
                    y = (real(ipt_eta -1,rk))/real(npts_eta -1,rk)
                    z = (real(ipt_zeta-1,rk))/real(npts_zeta-1,rk)

                    !xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    !ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    !zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 3x3x3 unit-element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------
    subroutine meshgen_3x3x3_unit_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 64
        real(rk), dimension(npt)    :: x,y,z

        ! elements (3x3x3) - unit linear
        !
        !            *-------*-------*-------*
        !           /       /       /       /|
        !          *-------*-------*-------* |
        !         /       /       /       /| *
        !        *-------*-------*-------* |/|
        !       /       /       /       /| * |
        !      *-------*-------*-------* |/| *
        !      |       |       |       | * |/|
        !      |       |       |       |/| * |
        !      *-------*-------*-------* |/| *
        !      |       |       |       | * |/
        !      |       |       |       |/| *
        !      *-------*-------*-------* |/
        !      |       |       |       | *
        !      |       |       |       |/
        !      *-------*-------*-------*
        !
        !
        x = [ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, &
             ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, &
             ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, &
             ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX, ZERO, TWO, FOUR, SIX]

        y = [ZERO, ZERO, ZERO, ZERO, TWO, TWO, TWO, TWO, FOUR, FOUR, FOUR, FOUR, SIX, SIX, SIX, SIX, &
             ZERO, ZERO, ZERO, ZERO, TWO, TWO, TWO, TWO, FOUR, FOUR, FOUR, FOUR, SIX, SIX, SIX, SIX, &
             ZERO, ZERO, ZERO, ZERO, TWO, TWO, TWO, TWO, FOUR, FOUR, FOUR, FOUR, SIX, SIX, SIX, SIX, &
             ZERO, ZERO, ZERO, ZERO, TWO, TWO, TWO, TWO, FOUR, FOUR, FOUR, FOUR, SIX, SIX, SIX, SIX]

        z = [ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, &
             TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, TWO, &
             FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, FOUR, &
             SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX, SIX]

        npts_xi   = 4
        npts_eta  = 4
        npts_zeta = 4

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 3x3x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_3x3x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 32
        real(rk), dimension(npt)    :: x,y,z

        ! elements (3x3x1) - linear
        !
        !            *-------*
        !           /       /|
        !          *-------* |
        !         /       /| *
        !        *-------* |/|
        !       /       /| * |
        !      *-------* |/| *
        !      |       | * |/|
        !      |       |/| * |
        !      *-------* |/| *
        !      |       | * |/
        !      |       |/| *
        !      *-------* |/
        !      |       | *
        !      |       |/
        !      *-------*
        !
        !
        x = [ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, &
             ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE, ZERO, ONE, TWO, THREE]

        y = [ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, THREE, THREE, THREE, THREE, &
             ZERO, ZERO, ZERO, ZERO, ONE, ONE, ONE, ONE, TWO, TWO, TWO, TWO, THREE, THREE, THREE, THREE]

        z = [ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE, ONE]

        npts_xi   = 4
        npts_eta  = 4
        npts_zeta = 2

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 4x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_4x1x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 20
        real(rk), dimension(npt)    :: x,y,z

        ! elements (4x1x1) - linear
        !
        !      *------*------*------*------*
        !      |      |      |      |      |
        !      |      |      |      |      |
        !      *------*------*------*------*
        !

        x = [ZERO, ONE, TWO, THREE, FOUR, &
             ZERO, ONE, TWO, THREE, FOUR, &
             ZERO, ONE, TWO, THREE, FOUR, &
             ZERO, ONE, TWO, THREE, FOUR]

        y = [ZERO, ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, ONE, &
             ZERO, ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, ONE]

        z = [ZERO, ZERO, ZERO, ZERO, ZERO, &
             ZERO, ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, ONE, &
             ONE, ONE, ONE, ONE, ONE]

        npts_xi   = 5
        npts_eta  = 2
        npts_zeta = 2

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 4x2x2 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_4x2x2_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        real(rk)    :: x,y,z
        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ierr, &
                       npts_xi, npts_eta, npts_zeta

        ! elements (4x2x2) - linear
        !  *---*---*---*---*
        !  |\   \   \   \   \
        !  * *---*---*---*---*
        !  |\|\   \   \   \   \
        !  * * *---*---*---*---*
        !   \|\|   |   |   |   |
        !    * *---*---*---*---*
        !     \|   |   |   |   |
        !      *---*---*---*---*
        !

        npts_xi   = 5
        npts_eta  = 3
        npts_zeta = 3

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    x = TWO*real(ipt_xi  -1,rk)/real(npts_xi  -1,rk)
                    y = ONE*real(ipt_eta -1,rk)/real(npts_eta -1,rk)
                    z = ONE*real(ipt_zeta-1,rk)/real(npts_zeta-1,rk)

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z
                end do
            end do
        end do

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

    !> Generate a set of points defining a 2x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_2x1x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 12
        real(rk), dimension(npt)    :: x,y,z

        ! elements (3x1x1) - linear
        !
        !      *------*------*
        !      |      |      |
        !      |      |      |
        !      *------*------*
        !

        x = [ZERO, ONE, TWO, &
             ZERO, ONE, TWO, &
             ZERO, ONE, TWO, &
             ZERO, ONE, TWO]

        y = [ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, &
             ZERO, ZERO, ZERO, &
             ONE, ONE, ONE]

        z = [ZERO, ZERO, ZERO, &
             ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, &
             ONE, ONE, ONE]

        npts_xi   = 3
        npts_eta  = 2
        npts_zeta = 2

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 3x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_3x1x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        integer(ik), parameter      :: npt = 16
        real(rk), dimension(npt)    :: x,y,z

        ! elements (3x1x1) - linear
        !
        !      *------*------*------*
        !      |      |      |      |
        !      |      |      |      |
        !      *------*------*------*
        !

        x = [ZERO, ONE, TWO, THREE, &
             ZERO, ONE, TWO, THREE, &
             ZERO, ONE, TWO, THREE, &
             ZERO, ONE, TWO, THREE]

        y = [ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, &
             ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE]

        z = [ZERO, ZERO, ZERO, ZERO, &
             ZERO, ZERO, ZERO, ZERO, &
             ONE, ONE, ONE, ONE, &
             ONE, ONE, ONE, ONE]

        npts_xi   = 4
        npts_eta  = 2
        npts_zeta = 2

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        ipt = 1
        do ipt_zeta = 1,npts_zeta
            do ipt_eta = 1,npts_eta
                do ipt_xi = 1,npts_xi
                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x(ipt)
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y(ipt)
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z(ipt)
                    ipt = ipt + 1
                end do
            end do
        end do

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

    !> Generate a set of points defining a 4x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !--------------------------------------------------------------------------------------
    subroutine meshgen_40x15x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        real(rk)    :: x,y,z, dx, dy, dz

        ! elements (40x15x1) - linear

        npts_xi   = 41
        npts_eta  = 16
        npts_zeta = 2

        dx = 0.5_rk
        dy = 0.5_rk
        dz = 1._rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = ZERO
        do ipt_zeta = 1,npts_zeta
            y = ZERO
            do ipt_eta = 1,npts_eta
                x = ZERO
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

    !> Generate a set of points defining a 4x1x1 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_15x15x1_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        real(rk)    :: x,y,z,dx,dy,dz

        ! elements (15x15x1) - linear

        npts_xi   = 16
        npts_eta  = 16
        npts_zeta = 2

        dx = 0.5_rk
        dy = 0.5_rk
        dz = 1.0_rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = ZERO
        do ipt_zeta = 1,npts_zeta
            y = ZERO
            do ipt_eta = 1,npts_eta
                x = ZERO
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z
                    ipt = ipt + 1

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

    !> Generate a set of points defining a 15x15x2 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_15x15x2_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi,npts_eta,npts_zeta
        real(rk)    :: x,y,z,dx,dy,dz

        ! elements (15x15x2) - linear

        npts_xi   = 16
        npts_eta  = 16
        npts_zeta = 3

        dx = 0.5_rk
        dy = 0.5_rk
        dz = 0.5_rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = ZERO
        do ipt_zeta = 1,npts_zeta
            y = ZERO
            do ipt_eta = 1,npts_eta
                x = ZERO
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

    !> Generate a set of points defining a 15x15x3 element mesh
    !!
    !!  @author Nathan A. Wukie
    !!  @param[inout]   pts     points_t array of rank-3 that gets allocated,
    !!                          filled, and returned
    !---------------------------------------------------------------------------------------
    subroutine meshgen_15x15x3_linear(xcoords,ycoords,zcoords)
        real(rk),   allocatable :: xcoords(:,:,:)
        real(rk),   allocatable :: ycoords(:,:,:)
        real(rk),   allocatable :: zcoords(:,:,:)

        integer(ik) :: ipt_xi, ipt_eta, ipt_zeta, ipt, ierr, &
                       npts_xi, npts_eta, npts_zeta

        real(rk)    :: x,y,z, dx, dy, dz

        ! elements (15x15x3) - linear

        npts_xi   = 16
        npts_eta  = 16
        npts_zeta = 4

        dx = 0.5_rk
        dy = 0.5_rk
        dz = 0.5_rk

        allocate(xcoords(npts_xi,npts_eta,npts_zeta), ycoords(npts_xi,npts_eta,npts_zeta), &
                 zcoords(npts_xi,npts_eta,npts_zeta), stat=ierr)
        if (ierr /= 0) call AllocationError

        z = ZERO
        do ipt_zeta = 1,npts_zeta
            y = ZERO
            do ipt_eta = 1,npts_eta
                x = ZERO
                do ipt_xi = 1,npts_xi

                    xcoords(ipt_xi,ipt_eta,ipt_zeta) = x
                    ycoords(ipt_xi,ipt_eta,ipt_zeta) = y
                    zcoords(ipt_xi,ipt_eta,ipt_zeta) = z

                    x = x + dx
                end do
                y = y + dy
            end do
            z = z + dz
        end do

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

end module mod_gridgen_blocks

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998 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 interpoltet(x,y,z,xo,yo,zo,nx,ny,nz,planfa,ifatet,
     &     netet,kontet,cotet,iparent,co,nkfa,nkfb,konl,ratio,ikf)
!
      implicit none
!
!     10 nearest nodes: neighbor(10)
!
      integer ifatet(4,*),neighbor(10),near,nx(*),ny(*),nz(*),
     &     ifs,iface,i,konl(4,*),ielmax,netet,m,ikf(*),j,kontet(*),
     &     iparent(*),nterms,indexe,nelem,ii,nkfa,nkfb,loopa,
     &     jj,inside,konl_opt(4),k
!
      real*8 cotet(3,*),planfa(4,*),dface,dfacemax,tolerance,dist,
     &     ratio(4,*),pneigh(3,4),pnode(3),xi,et,ze,dist_opt,
     &     xp,yp,zp,x(*),y(*),z(*),xo(*),yo(*),zo(*),co(3,*),
     &     ratio_opt(4)
!
      tolerance=1.d-6
      nterms=4
      loopa=8
c      write(*,*) 'interpoltet'
!
      do m=nkfa,nkfb
!
        jj=ikf(m)
        xp=co(1,jj)
        yp=co(2,jj)
        zp=co(3,jj)
!
        do ii=1,2
!
          if(ii.eq.1) then
            near=1
          else
            near=10
          endif
          call near3d(xo,yo,zo,x,y,z,nx,ny,nz,xp,yp,zp,netet,neighbor,
     &near)
!
          inside=0
!
!         a node is inside the tetrahedron if the substitution
!         of its coordinates in the equation of the tet faces
!         (corrected by ifs/abs(ifs)) is for all faces positive
!
          do i=1,near
            nelem=iparent(neighbor(i))
            dface=0.d0
            do j=1,4
              ifs=ifatet(j,nelem)
              iface=abs(ifs)
              dist=planfa(1,iface)*xp+planfa(2,iface)*yp
     &             +planfa(3,iface)*zp+planfa(4,iface)
              if(dist*ifs.lt.-1.d-10*iface) then
                dface=dface+dist*ifs/iface
              endif
            enddo
            if(dface.gt.-1.d-10) then
              inside=1
c              write(*,*) 'interpoltet inside'
              exit
            endif
            if(i.eq.1) then
              dfacemax=dface
              ielmax=nelem
            else
              if(dface.gt.dfacemax) then
                ielmax=nelem
                dfacemax=dface
              endif
            endif
          enddo
!
!     if no element was found, the element with the smallest
!     dfacemax (summed distance) is taken
!
          if(inside.eq.0) then
            nelem=ielmax
          endif
!
          indexe=4*(nelem-1)
!
          do i=1,4
            konl(i,m)=kontet(indexe+i)
          enddo
!
!     nodes of master element
!
          do i=1,4
            do j=1,3
              pneigh(j,i)=cotet(j,konl(i,m))
            enddo
          enddo
!
!     slave node
!
          pnode(1)=xp
          pnode(2)=yp
          pnode(3)=zp
!
!     attaching slave node to master element
!
          call attach_3d(pneigh,pnode,nterms,ratio(1,m),dist,xi,et,ze,
     &         loopa)
c          write(*,*) nelem,dist,tolerance
!
!     checking the parent elements of the "best" tetrahedra
!     in case the distance between slave node and location of
!     interpolation exceeds "tolerance"
!
          if(dist.gt.tolerance) then
            if(ii.eq.1) cycle
            do i=1,4
              konl_opt(i)=konl(i,m)
              ratio_opt(i)=ratio(i,m)
            enddo
            dist_opt=dist
!
            do k=1,near
!
!     slave node
!
              pnode(1)=xp
              pnode(2)=yp
              pnode(3)=zp
!
              nelem=iparent(neighbor(k))
              indexe=4*(nelem-1)
!
              do i=1,4
                konl(i,m)=kontet(indexe+i)
              enddo
!
!     nodes of master element
!
              do i=1,4
                do j=1,3
                  pneigh(j,i)=cotet(j,konl(i,m))
                enddo
              enddo
!
!     attaching slave node to master element
!
              call attach_3d(pneigh,pnode,nterms,ratio(1,m),dist,
     &             xi,et,ze,loopa)
!
!     check whether the present element yields better results
!
              if(dist.lt.dist_opt) then
                do i=1,4
                  konl_opt(i)=konl(i,m)
                  ratio_opt(i)=ratio(i,m)
                enddo
                dist_opt=dist
              endif
              if(dist.lt.tolerance) exit
            enddo
!
!     storing the optimal configuration
!
            do i=1,4
              konl(i,m)=konl_opt(i)
              ratio(i,m)=ratio_opt(i)
            enddo
          endif
!
          if((ii.eq.2).or.(dist.lt.tolerance)) exit
!
c          if(dist.lt.tolerance) exit
!
        enddo
c        write(*,*) jj,dist,co(1,jj),ratio(1,m)*pneigh(1,1)+
c     &       ratio(2,m)*pneigh(1,2)+ratio(3,m)*pneigh(1,3)+
c     &       ratio(4,m)*pneigh(1,4)
!
      enddo
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
  ! -------
  function nf90_inq_libvers()
    character(len = 80) :: nf90_inq_libvers

    nf90_inq_libvers = nf_inq_libvers()
  end function nf90_inq_libvers
  ! -------
  function nf90_strerror(ncerr)
    integer, intent( in) :: ncerr
    character(len = 80)  :: nf90_strerror

    nf90_strerror = nf_strerror(ncerr)
  end function nf90_strerror
  ! -------
  !
  ! File level control routines:
  !
  function nf90_inq_base_pe(ncid, pe)
    integer, intent( in) :: ncid
    integer, intent(out) :: pe
    integer              :: nf90_inq_base_pe

    nf90_inq_base_pe = nf_inq_base_pe(ncid, pe)
  end function nf90_inq_base_pe
  ! -------
  function nf90_set_base_pe(ncid, pe)
    integer, intent( in) :: ncid, pe
    integer              :: nf90_set_base_pe

    nf90_set_base_pe = nf_set_base_pe(ncid, pe)
  end function nf90_set_base_pe
  ! -------
  function nf90_create(path, cmode, ncid, initialsize, chunksize)
    character (len = *), intent(in   ) :: path
    integer,             intent(in   ) :: cmode
    integer,             intent(  out) :: ncid
    integer, optional,   intent(in   ) :: initialsize
    integer, optional,   intent(inout) :: chunksize
    integer                            :: nf90_create

    integer :: fileSize, chunk

    if(.not. (present(initialsize) .or. present(chunksize)) ) then
      nf90_create = nf_create(path, cmode, ncid)
    else
      ! Default values per man page
      filesize = 0; chunk = nf90_sizehint_default
      if(present(initialsize)) filesize = initialsize
      if(present(chunksize  )) chunk    = chunksize
      nf90_create = nf__create(path, cmode, filesize, chunk, ncid)
      ! Pass back the value actually used
      if(present(chunksize  )) chunksize = chunk
    end if
  end function nf90_create
  ! -------
  function nf90_create_mp(path, cmode, initalsz, basepe, chunksizehint, ncid)
    character (len = *), intent( in) :: path
    integer,             intent( in) :: cmode, initalsz, basepe, chunksizehint
    integer,             intent(out) :: ncid
    integer                          :: nf90_create_mp

    nf90_create_mp = nf__create_mp(path, cmode, initalsz, basepe, chunksizehint, ncid)
  end function nf90_create_mp
  ! -------
  function nf90_open(path, mode, ncid, chunksize)
    character (len = *), intent(in   ) :: path
    integer,             intent(in   ) :: mode
    integer,             intent(  out) :: ncid
    integer, optional,   intent(inout) :: chunksize
    integer                            :: nf90_open

    if(present(chunksize)) then
      nf90_open = nf__open(path, mode, chunksize, ncid)
    else
      nf90_open = nf_open(path, mode, ncid)
    end if
  end function nf90_open
  ! -------
  function nf90_open_mp(path, mode, basepe, chunksizeint, ncid)
    character (len = *), intent( in) :: path
    integer,             intent( in) :: mode, basepe, chunksizeint
    integer,             intent(out) :: ncid
    integer                          :: nf90_open_mp

    nf90_open_mp = nf__open_mp(path, mode, basepe, chunksizeint, ncid)
  end function nf90_open_mp
  ! -------
  function nf90_set_fill(ncid, fillmode, old_mode)
    integer, intent( in) :: ncid, fillmode
    integer, intent(out) :: old_mode
    integer              :: nf90_set_fill

    nf90_set_fill = nf_set_fill(ncid, fillmode, old_mode)
  end function nf90_set_fill
  ! -------
  function nf90_redef(ncid)
    integer, intent( in) :: ncid
    integer              :: nf90_redef

    nf90_redef = nf_redef(ncid)
  end function nf90_redef
   ! -------
  function nf90_enddef(ncid, h_minfree, v_align, v_minfree, r_align)
    integer,           intent( in) :: ncid
    integer, optional, intent( in) :: h_minfree, v_align, v_minfree, r_align
    integer                        :: nf90_enddef

    integer :: hMinFree, VAlign, VMinFree, RAlign

    if(.not. any( (/ present(h_minfree), present(v_align), &
                     present(v_minfree), present(r_align) /) ) )then
      nf90_enddef = nf_enddef(ncid)
    else
      ! Default values per the man page
      hMinFree = 0; VMinFree = 0
      VAlign   = 4; RAlign   = 4
      if(present(h_minfree)) HMinFree = h_minfree
      if(present(v_align  )) VAlign   = v_align
      if(present(v_minfree)) VMinFree = v_minfree
      if(present(r_align  )) RAlign   = r_align
      nf90_enddef = nf__enddef(ncid, hMinFree, VAlign, VMinFree, RAlign)
    end if
  end function nf90_enddef
  ! -------
  function nf90_sync(ncid)
    integer, intent( in) :: ncid
    integer              :: nf90_sync

    nf90_sync = nf_sync(ncid)
  end function nf90_sync
  ! -------
  function nf90_abort(ncid)
    integer, intent( in) :: ncid
    integer              :: nf90_abort

    nf90_abort = nf_abort(ncid)
  end function nf90_abort
  ! -------
  function nf90_close(ncid)
    integer, intent( in) :: ncid
    integer              :: nf90_close

    nf90_close = nf_close(ncid)
  end function nf90_close
  ! -------
  function nf90_delete(name)
    character(len = *), intent( in) :: name
    integer                         :: nf90_delete

    nf90_delete = nf_delete(name)
  end function nf90_delete

  !
  ! A single file level inquiry routine
  !
  function nf90_inquire(ncid, nDimensions, nVariables, nAttributes, unlimitedDimId, formatNum)
    integer,           intent( in) :: ncid
    integer, optional, intent(out) :: nDimensions, nVariables, nAttributes, unlimitedDimId, formatNum
    integer                        :: nf90_inquire

    integer :: nDims, nVars, nGAtts, unlimDimId, frmt

    nf90_inquire = nf_inq(ncid, nDims, nVars, nGAtts, unlimDimId)
    if(present(nDimensions))    nDimensions    = nDims
    if(present(nVariables))     nVariables     = nVars
    if(present(nAttributes))    nAttributes    = nGAtts
    if(present(unlimitedDimId)) unlimitedDimId = unlimDimId
    if(present(formatNum)) then
       nf90_inquire = nf_inq_format(ncid, frmt)
       formatNum = frmt
    endif
  end function nf90_inquire

max_stars_repo_name: liujiamingustc/phd


text:
      module MOD_PRIINT
      contains
C**********  MODULNAME: PRIINT    ******* 24/03/87  21.29.40.******    34 KARTEN
      SUBROUTINE PRIINT (XJC,XJCARR,EDDI,EDDARR,R,ND,XLAMBDA,NF,LSTEP,
     $                   JOBNUM,MODHEAD)
      use MOD_TRADFUN
      use MOD_extrxjc
      IMPLICIT REAL*8(A-H,O-Z)

      DIMENSION XJC(ND),XLAMBDA(NF),EDDI(3,ND),R(ND)
	DIMENSION XJCARR(ND,NF),EDDARR(3,ND,NF)
      CHARACTER MODHEAD*104

      PRINT 2,MODHEAD,JOBNUM
    2 FORMAT (10H1$$$$$$$$$,/,1X,  A  ,20X,'JOB NO.',I5,
     $            //,10X,'MEAN INTENSITY',/,10X,14('-'),/,
     $ ' FREQUENCY     DEPTH       J-NUE          T-RAD   ',
     $ '         EDDINGTON     SPHERICITY      EDDINGTON',/,
     $ '   INDEX       INDEX     (ERG/CM+2)       (KELVIN)',
     $ '         FACTOR F      FACTOR  R.R.Q    FACTOR H/J',/)

      DO 4 K=1,NF
c      ENCODE (8,FMT_KEY,NAME) 'XJC ',K
cc    3 FORMAT (A4,I4)
c      CALL READMS (3,XJC,ND,NAME)
c      ENCODE (7,3,NAME) 'EDDI',K
c      CALL READMS (3,EDDI,3*ND,NAME)
      call extrxjc (XJCARR,XJC,EDDARR,EDDI,nd,nf,K)
      DO 1 L=1,ND
      IF (((L-1)/LSTEP)*LSTEP .NE. (L-1) .AND. L.NE.ND) GOTO 1
      TRAD=TRADFUN (XLAMBDA(K),XJC(L))
      IF (L.EQ.ND .OR. L.EQ.1) THEN
      PRINT 5,K,L,XJC(L),TRAD,EDDI(1,L),EDDI(2,L)*R(L)*R(L),EDDI(3,L)
    5 FORMAT (2I10,E15.3,F15.0,F15.3,F15.3,F15.3)
      ELSE
      PRINT 5,K,L,XJC(L),TRAD,EDDI(1,L),EDDI(2,L)*R(L)*R(L)
      ENDIF
    1 CONTINUE
    4 PRINT 6
    6 FORMAT (1X)

      RETURN
      END subroutine
      end module

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


text:
! { dg-do run }
! Checks the LOGICAL version of dot_product
!
! Contributed by Paul Thomas  <pault@gcc.gnu.org>
!
   logical :: l1(4) = (/.TRUE.,.FALSE.,.TRUE.,.FALSE./)
   logical :: l2(4) = (/.FALSE.,.TRUE.,.FALSE.,.TRUE./)
   if (dot_product (l1, l2)) call abort ()
   l2 = .TRUE.
   if (.not.dot_product (l1, l2)) call abort ()
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
	SUBROUTINE SFESHP  ( lunedt, isffln, nparm, iploc, iret )
C************************************************************************
C* SFESHP								*
C*									*
C* This subroutine reads data from an edit text file and writes it	*
C* to a non-standard surface file.					*
C*									*
C* SFESHP  ( LUNEDT, ISFFLN, NPARM, IPLOC, IRET )			*
C*									*
C* Input parameters:							*
C*	LUNEDT		INTEGER		LUN for edit file		*
C*	ISFFLN		INTEGER		Surface file number		*
C*	NPARM		INTEGER		Number of parameters		*
C*	IPLOC(*)	INTEGER		Location of data in output file	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal			*
C** Log:								*
C* S. Jacobs/NMC	10/94		Copied from SFEDTA		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	INTEGER		iploc (*)
C*
	CHARACTER	record*132, tim*20, stn*8, stid*8
	CHARACTER	curtim*20, cdata (MMPARM)*20
	REAL		rdata ( MMPARM ), data ( MMPARM )
	LOGICAL		first
C------------------------------------------------------------------------
	iret = 0
C
C*	Read all the records in the file.
C
	ier    = 0
	first  = .true.
	knt    = 0
	curtim = ' '
C*
	DO WHILE  ( ier .eq. 0 )
C
C*	    Read the next record.
C
	    READ   ( lunedt, 10, IOSTAT = ier )  record
10	    FORMAT ( A )
C
C*	    Skip blank records.
C
	    iskp = INDEX ( record, 'YYMMDD' )
	    IF  ( iskp .ne. 0 )  record = ' '
	    IF  ( first )  THEN
		iskp = INDEX ( record, '/' )
		IF  ( iskp .eq. 0 )  record = ' '
	    END IF
	    IF  ( ( ier .eq. 0 ) .and. ( record .ne. ' ' ) )  THEN
C
C*		Break record into character strings.
C
		CALL ST_C2C  ( record, MMPARM, cdata, ndat, ier )
C
C*		Get station and time from first data record.
C
		IF  ( first )  THEN
		    stn   = cdata (1)
		    tim   = cdata (2) (1:6) // '/' // cdata (3) (1:4)
		    first = .false.
		    i     = 4
		  ELSE
		    i     = 1
		END IF
C
C*		Convert rest of variables into real numbers.
C
		DO WHILE  ( ( i .le. ndat ) .and. ( knt .lt. nparm ) )
		    knt = knt + 1
		    CALL ST_CRNM  ( cdata (i), rdata (knt), ier )
		    i = i + 1
		END DO
C
C*		If all the variables are found, write to file.
C
		IF  ( knt .ge. nparm )  THEN
C
C*		    Reset counters.
C
		    first = .true.
		    knt   = 0
C
C*		    Set the lat/lon from the edit file data array.
C
		    rlat = rdata ( iploc (1) )
		    rlon = rdata ( iploc (2) )
C
C*		    Fill the rest of the data array.
C
		    DO  i = 3, nparm
			data ( iploc (i) ) = rdata (i)
		    END DO
C
C*		    Write the data to the file.
C
		    CALL ST_NUMB ( stn, isnum, ierr )
		    IF  ( ierr .eq. 0 )  THEN
			stid  = ' '
		    ELSE
			stid  = stn
			isnum = RMISSD
		    END IF
		    CALL SF_WSDD  ( isffln, tim, stn, isnum,
     +				    rlat, rlon, 0, ' ', ' ',
     +				    ihhmm, data, ier )
		END IF
	    END IF
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak