text:
        WRITE (7,350)
  999	format(5(f16.8,x))
c
c
c
c	..attempt to calculate density fluctuations from the continuity
c		equation	11/23/92
c
c
c	..fixed up on 6/18/93
c
c
	emag = sqrt(e(1)**2. + e(2)**2. + e(3)**2. + e(4)**2. + e(5)**2.
	1		+ e(6)**2.)
c
c	print*,'emag = ',emag		!this should be 1
c
	del_nr = 0
	del_ni = 0
c
c********************
c	ALL, N = 1,NC
c
c	..first calculate the complex matrix product
c
c
	DO N = 1,NC
	   do jj=1,3
		del_nr(N) = del_nr(N) + 			!real part
     1		XKX*(dieli(N,1,jj)*E(jj) + dielr(N,1,jj)*E(jj+3))
     1		+ XKZ*(dieli(N,3,jj)*E(jj)+ dielr(N,3,jj)*E(jj+3))
c
		del_ni(N) = del_ni(N) +				!imaginary part
     1		XKX*(dielr(N,1,jj)*E(jj) - dieli(N,1,jj)*E(jj+3))
     1		+ XKZ*(dielr(N,3,jj)*E(jj) - dielr(N,3,jj)*E(jj+3))
c
	   end do
c
c	..then include the real scalars
c
c
	del_nr(N) = del_nr(N)*ratiom(N)/(dens(N)*abs(ratiom(N))
     1		*(WR**2.))
	del_ni(N) = del_ni(N)*ratiom(N)/(dens(N)*abs(ratiom(N))
     1		*( WR**2.))
	ENDDO
c
c****************

ext: for


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

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

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

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

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

! End MIT license text.

      SUBROUTINE BCHECK_2D ( B, BTYPE, ID, NROWB, NCOLB, NUM_GRIDS, XB, XL, BW )

! Checks strain-displacement matrices for rigid body motion and constant strain for 2-D shell elements
! (6 DOF per grid point and up to 4 grid points)

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  BUG, F04, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  BCHECK_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO, TWO
      USE MODEL_STUF, ONLY            :  ELDOF, NELGP, TE
      USE MODEL_STUF, ONLY            :  AGRID, ELGP

      USE BCHECK_2D_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'BCHECK_2D'
      CHARACTER(LEN=*), INTENT(IN)    :: BTYPE             ! Type of B matrix ('M' for membrane, 'B' for bending, 'S' for shear)
      CHARACTER(45*BYTE)              :: MESSAG(14)        ! Output messages for the 14 modes of deformation (6 RB + 8 const strain)

      INTEGER(LONG), INTENT(IN)       :: NROWB             ! Number of rows in the input B matrix
      INTEGER(LONG), INTENT(IN)       :: NCOLB             ! Number of cols in the input B matrix
      INTEGER(LONG), INTENT(IN)       :: NUM_GRIDS         ! Number of grids for the input B matrix
      INTEGER(LONG), INTENT(IN)       :: ID(NCOLB)         ! List of elem DOF's for each of the elem grids (e.g 3,4,5 for each of
!                                                            4 grids for a 4 node plate bending elem
      INTEGER(LONG)                   :: I,J,K             ! DO loop indices
      INTEGER(LONG)                   :: KK                ! A computed index into array W
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = BCHECK_BEGEND

      REAL(DOUBLE) , INTENT(IN)       :: B(NROWB,NCOLB)    ! Strain-displ matrix
      REAL(DOUBLE) , INTENT(IN)       :: XB(NUM_GRIDS,3)   ! Basic coords of elem grids (diff than XEB for TPLT2's in a MIN4T QUAD4)
      REAL(DOUBLE) , INTENT(IN)       :: XL(NUM_GRIDS,3)   ! Local coords of elem grids (diff than XEL for TPLT2's in a MIN4T QUAD4)
      REAL(DOUBLE) , INTENT(OUT)      :: BW(NROWB,14)      ! Output from subr BCHECK_2D (matrix of NROWB elem strains for various
!                                                            elem rigid body motions/constant strain distortions)
      REAL(DOUBLE)                    :: GRD_COORDS(3)     ! 3 coords from XEB for one node of the element
      REAL(DOUBLE)                    :: REF_COORDS(3)     ! 3 coords from XEB for node 1
      REAL(DOUBLE)                    :: RB_DISP(6,6)      ! 6 x 6 RB matrix for one grid for this element
      REAL(DOUBLE)                    :: W(24,14)          ! Displs for the 14 modes of elem deformation (6 RB + 8 constant strain)

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

! **********************************************************************************************************************************
! Initialize outputs

      DO I=1,NROWB
         DO J=1,14
            BW(I,J) = ZERO
         ENDDO
      ENDDO

! Init W

      DO I=1,24
         DO J=1,14
            W(I,J) = ZERO
         ENDDO
      ENDDO

      MESSAG( 1) = '      Rigid Body Displacement,      T1  = 1.0'
      MESSAG( 2) = '      Rigid Body Displacement,      T2  = 1.0'
      MESSAG( 3) = '      Rigid Body Displacement,      T3  = 1.0'
      MESSAG( 4) = '      Rigid Body Rotation,          R1  = 1.0'
      MESSAG( 5) = '      Rigid Body Rotation,          R2  = 1.0'
      MESSAG( 6) = '      Rigid Body Rotation,          R3  = 1.0'
      MESSAG( 7) = '      Constant in-plane strain,     Exx = 1.0'
      MESSAG( 8) = '      Constant in-plane strain,     Eyy = 1.0'
      MESSAG( 9) = '      Constant in-plane strain,     Exy = 1.0'
      MESSAG(10) = '      Constant curvature,           Cxx = 1.0'
      MESSAG(11) = '      Constant curvature,           Cyy = 1.0'
      MESSAG(12) = '      Constant curvature,           Cxy = 1.0'
      MESSAG(13) = '      Constant transv shear strain, Gxz = 1.0'
      MESSAG(14) = '      Constant transv shear strain, Gyz = 1.0'

! Calc RB modes (cols 1-6 of W)

      DO I=1,NUM_GRIDS
         DO J=1,3
            GRD_COORDS(J) = XB(I,J)
            REF_COORDS(J) = XB(1,J)
         ENDDO
         CALL RIGID_BODY_DISP_MAT ( GRD_COORDS, REF_COORDS, RB_DISP )
         DO J=1,6
            DO K=1,6
               W(6*(I-1)+J,K) = RB_DISP(J,K)
            ENDDO
         ENDDO
      ENDDO

! Calc constant strain modes (cols 7-14 of W)

      DO K=1,NUM_GRIDS

         KK = 6*(K-1)

         W(KK+1, 7) =  XL(K,1)                            ! This gives constant direct x strain

         W(KK+2, 8) =  XL(K,2)                            ! This gives constant direct y strain

         W(KK+1, 9) =  XL(K,2)/TWO                        ! The next 2 give constant shear xy strain
         W(KK+2, 9) =  XL(K,1)/TWO

         W(KK+3,10) = -XL(K,1)*XL(K,1)/TWO                ! This gives constant curvature, Cxx
         W(KK+5,10) =  XL(K,1)

         W(KK+3,11) = -XL(K,2)*XL(K,2)/TWO                ! This gives constant curvature, Cyy
         W(KK+4,11) = -XL(K,2)

         W(KK+3,12) = -XL(K,1)*XL(K,2)/TWO                ! This gives constant curvature, Cxy
         W(KK+4,12) = -XL(K,1)/TWO
         W(KK+5,12) =  XL(K,2)/TWO

         W(KK+3,13) =  XL(K,1)                            ! This gives constant transverse shear strain

         W(KK+3,14) =  XL(K,2)                            ! This gives constant transverse shear strain

      ENDDO

! Calc BW = B*W (but B has fewer cols since element has no stiffness for some DOF's)

      DO I=1,NROWB
         DO J=1,14
            BW(I,J) = ZERO
            DO K=1,NCOLB
               BW(I,J) = BW(I,J) + B(I,K)*W(ID(K),J)
            ENDDO
         ENDDO
      ENDDO

! Write results

      IF ((BTYPE == 'B') .OR. (BTYPE == 'M')) THEN
         DO J=1,6
            WRITE(BUG,9101) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ELSE IF (BTYPE == 'S') THEN
         DO J=1,6
            WRITE(BUG,9201) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ENDIF
      WRITE(BUG,*)

      IF (BTYPE == 'B') THEN
         DO J=7,9
            WRITE(BUG,9101) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ELSE IF (BTYPE == 'M') THEN
         WRITE(BUG,9102) MESSAG(7),(BW(I,7),I=1,NROWB)
         WRITE(BUG,9103) MESSAG(8),(BW(I,8),I=1,NROWB)
         WRITE(BUG,9104) MESSAG(9),(BW(I,9),I=1,NROWB)
      ELSE IF (BTYPE == 'S') THEN
         DO J=7,9
            WRITE(BUG,9201) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ENDIF
      WRITE(BUG,*)

      IF (BTYPE == 'B') THEN
         WRITE(BUG,9102) MESSAG(10),(BW(I,10),I=1,NROWB)
         WRITE(BUG,9103) MESSAG(11),(BW(I,11),I=1,NROWB)
         WRITE(BUG,9104) MESSAG(12),(BW(I,12),I=1,NROWB)
      ELSE IF (BTYPE == 'M') THEN
         DO J=10,12
            WRITE(BUG,9101) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ELSE IF (BTYPE == 'S') THEN
         DO J=10,12
            WRITE(BUG,9201) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ENDIF
      WRITE(BUG,*)

      IF ((BTYPE == 'B') .OR. (BTYPE == 'M')) THEN
         DO J=13,14
           WRITE(BUG,9101) MESSAG(J),(BW(I,J),I=1,NROWB)
         ENDDO
      ELSE IF (BTYPE == 'S') THEN
         WRITE(BUG,9202) MESSAG(13),(BW(I,13),I=1,NROWB)
         WRITE(BUG,9203) MESSAG(14),(BW(I,14),I=1,NROWB)
      ENDIF
      WRITE(BUG,*)

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

      RETURN

! **********************************************************************************************************************************
 9101 FORMAT(1X,A,3(1ES15.6),'  (should be 0, 0, 0)')

 9102 FORMAT(1X,A,3(1ES15.6),'  (should be 1, 0, 0)')

 9103 FORMAT(1X,A,3(1ES15.6),'  (should be 0, 1, 0)')

 9104 FORMAT(1X,A,3(1ES15.6),'  (should be 0, 0, 1)')

 9201 FORMAT(1X,A,2(1ES15.6),'  (should be 0, 0)')

 9202 FORMAT(1X,A,2(1ES15.6),'  (should be 1, 0)')

 9203 FORMAT(1X,A,2(1ES15.6),'  (should be 0, 1)')

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

      END SUBROUTINE BCHECK_2D


text:
!
! The Laboratory of Algorithms
!
! The MIT License
!
! Copyright 2011-2016 Andrey Pudov.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the 'Software'), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!

submodule (Foundation) Array
contains
    !
    ! Initializes a newly allocated array.
    !
    module subroutine array_init(self)
        class(Array), intent(in out) :: self
    end subroutine

    !
    ! Releases allocated array.
    !
    module subroutine array_destroy(self)
        class(Array), intent(in out) :: self

        integer index

        if (self%selfAllocated) then
            do index = 1, size(self%list)
                deallocate(self%list(index)%link)
            end do
        end if

        deallocate(self%list)
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_character(self, list)
        class(Array), intent(in out)        :: self
        character, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_complex(self, list)
        class(Array), intent(in out)      :: self
        complex, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_double(self, list)
        class(Array), intent(in out)               :: self
        double precision, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_integer(self, list)
        class(Array), intent(in out)      :: self
        integer, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_logical(self, list)
        class(Array), intent(in out)      :: self
        logical, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! Initializes a newly allocated array by placing in it the objects contained in a given Fortran array.
    !
    module subroutine array_initWithFArray_real(self, list)
        class(Array), intent(in out)   :: self
        real, dimension(:), intent(in) :: list

        class(Number), pointer :: buffer
        integer index

        allocate(self%list(size(list)))
        self%selfAllocated = .true.

        do index = 1, size(list)
            allocate(buffer)

            buffer = list(index)
            self%list(index)%link => buffer
        end do
    end subroutine

    !
    ! The number of objects in the array.
    !
    module function array_count(self) result(value)
        class(Array), intent(in) :: self
        integer                  :: value

        value = size(self%list)
    end function

    !
    ! Returns the object located at the specified index.
    !
    module function array_objectAtIndex(self, index) result(value)
        class(Array), intent(in) :: self
        integer, intent(in)      :: index
        class(Object), pointer   :: value

        if ((index < 1) .or. (index > size(self%list))) then
            ! TODO
        end if

        value => self%list(index)%link
    end function

    !
    ! The new array contains references to the receiving array’s elements, not copies of them.
    !
    module function array_sortedArrayUsingFunction(self, comparator) result(value)
        class(Array), intent(in) :: self
        interface
            function comparator(value1, value2) result(value)
                import Object

                class(Object), pointer, intent(in) :: value1
                class(Object), pointer, intent(in) :: value2
                integer                   :: value
            end function
        end interface
        class(Array), pointer :: value

        type(ObjectLink), dimension((size(self%list) + 1) / 2) :: buffer

        allocate(value)
        allocate(value%list, source = self%list)

        call mergeSort(value%list, buffer, size(value%list), comparator)
    end function

    recursive subroutine mergeSort(list, buffer, length, comparator)
        integer, intent(in)                                           :: length
        type(ObjectLink), dimension(length), intent(in out)           :: list
        type(ObjectLink), dimension((length + 1) / 2), intent(in out) :: buffer
        interface
            function comparator(value1, value2) result(order)
                import Object

                class(Object), pointer, intent(in) :: value1
                class(Object), pointer, intent(in) :: value2
                integer                            :: order
            end function
        end interface

        type(ObjectLink) link
        integer :: NA, NB

        select case (length)
        case (:1)
            return
        case (2)
            if (comparator(list(1)%link, list(2)%link) == ORDERED_DESCENDING) then
                link = list(1)
                list(1) = list(2)
                list(2) = link
            end if

            return
        case default
            NA = (length + 1) / 2
            NB = length - NA

            call mergeSort(list, buffer, NA, comparator)
            call mergeSort(list(NA + 1), buffer, NB, comparator)

            if (comparator(list(NA)%link, list(NA + 1)%link) == ORDERED_DESCENDING) then
                buffer(1:NA) = list(1:NA)
                call merge(buffer, NA, list(NA + 1), NB, list, length, comparator)
            endif
            return
        end select
    end subroutine

    subroutine merge(A, NA, B, NB, C, NC, comparator)
        integer, intent(in) :: NA,NB,NC           ! Normal usage: NA+NB = NC
        type(ObjectLink), intent(in out) :: A(NA) ! B overlays C(NA+1:NC)
        type(ObjectLink), intent(in)     :: B(NB)
        type(ObjectLink), intent(in out) :: C(NC)
        interface
            function comparator(value1, value2) result(order)
                import Object

                class(Object), pointer, intent(in) :: value1
                class(Object), pointer, intent(in) :: value2
                integer                            :: order
            end function
        end interface

        integer :: I,J,K
        integer order

        I = 1; J = 1; K = 1;
        do while(I <= NA .and. J <= NB)
            order = comparator(A(I)%link, B(J)%link)
            if ((order == ORDERED_ASCENDING) .or. (order == ORDERED_SAME)) then
                C(K) = A(I)
                I = I+1
            else
                C(K) = B(J)
                J = J+1
            endif

            K = K + 1
        enddo

        do while (I <= NA)
            C(K) = A(I)
            I = I + 1
            K = K + 1
        enddo
    end subroutine
end submodule

ext: f

max_stars_repo_name: andreypudov/Algorithms

max_issues_repo_name: jlokimlin/Algorithms

max_forks_repo_name: jlokimlin/Algorithms


text:
C
C	$Id: g01msg.f,v 1.4 2008-07-27 00:21:05 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE G01MSG(IACT)
C
C  Put out MESSAGE elements.
C
      include 'g01prm.h'
      include 'gksin.h'
      include 'g01ins.h'
C
C
C  IACT = 0  for no action, or
C       = 1  for action
C
      INTEGER IACT,G01PBL
C
C  Find the first non-blank, non-null and set the length.
C
      ILEN = 1
      DO 10 I=80,1,-1
        IF (STR(I:I).NE.' ' .AND. STR(I:I).NE.CHAR(0)) THEN
          ILEN = I
          GO TO 20
        ENDIF
   10 CONTINUE
   20 CONTINUE
C
C  Put out opcode (class and id) and length.
C
      NBYTES = G01PBL (ILEN, 1+(MEFW-1)/8 )
      CALL GPUTNI (7, 1, NBYTES, IERR)
C
C  Put out action/no-action flag.
C
      CALL GPUTPR (IACT, MEFW, 1, RERR)
C
C  Put out the message string.
C
      CALL GPUTPS (STR(1:ILEN), ILEN, ILEN, 0, IERR)
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


!=============================================================================!
!
!  These routines generate the factored LHS for the Navier-Stokes Equations
!  using 4th order differencing.
!
!=============================================================================!
        subroutine lhs1f( mat, Ah, Dh, Vh11, dtl, v )
!
!  Form the LHS for the \xi direction
!
!=============================================================================!
        use global
        use buffer
        use stencil
        implicit none

        real :: mat(5,ndof,ndof,nx,ny), v(ndof,nx,ny)
        real :: Ah(ndof,ndof,nx,ny), Dh(ndof,ndof,nx,ny), Vh11(6,nx,ny)
        real :: dtl(nx,ny)

!.... first derivative operator

        real :: a1, a2, a3, a4, a5
        real :: b1, b2, b3, b4, b5

        integer :: i, j, idof, jdof
!=============================================================================!

!.... fourth-order stencil

        a1 = alfa * ga1  / dxi
        a2 = alfa * ga2  / dxi
        a3 = alfa * zero / dxi
        a4 = alfa * ga3  / dxi
        a5 = alfa * ga4  / dxi

        b1 =  alfa * da1 / dxi**2
        b2 =  alfa * da2 / dxi**2
        b3 =  alfa * da3 / dxi**2
        b4 =  alfa * da4 / dxi**2
        b5 =  alfa * da5 / dxi**2

        !$omp parallel do private (i,idof,jdof)
        loop_j: do j = 1, ny
        loop_i: do i = 1, nx

!.... \hat{A} term

        do idof = 1, ndof
          do jdof = 1, ndof
            mat(1,idof,jdof,i,j) = a1 * dtl(i,j) * Ah(idof,jdof,i,j)
            mat(2,idof,jdof,i,j) = a2 * dtl(i,j) * Ah(idof,jdof,i,j)
            mat(3,idof,jdof,i,j) = a3 * dtl(i,j) * Ah(idof,jdof,i,j)
            mat(4,idof,jdof,i,j) = a4 * dtl(i,j) * Ah(idof,jdof,i,j)
            mat(5,idof,jdof,i,j) = a5 * dtl(i,j) * Ah(idof,jdof,i,j)
          end do
        end do

!.... \hat{D} term

        do idof = 1, ndof
          do jdof = 1, ndof
            mat(3,idof,jdof,i,j) = mat(3,idof,jdof,i,j) + &
                                   alfa * dtl(i,j) * Dh(idof,jdof,i,j)
          end do
        end do

!.... \hat{V}_{\xi\xi} term

        mat(1,2,2,i,j) = mat(1,2,2,i,j) - b1 * dtl(i,j) * Vh11(1,i,j)
        mat(1,2,3,i,j) = mat(1,2,3,i,j) - b1 * dtl(i,j) * Vh11(5,i,j)
        mat(1,3,2,i,j) = mat(1,3,2,i,j) - b1 * dtl(i,j) * Vh11(6,i,j)
        mat(1,3,3,i,j) = mat(1,3,3,i,j) - b1 * dtl(i,j) * Vh11(2,i,j)
        mat(1,4,4,i,j) = mat(1,4,4,i,j) - b1 * dtl(i,j) * Vh11(3,i,j)
        mat(1,5,5,i,j) = mat(1,5,5,i,j) - b1 * dtl(i,j) * Vh11(4,i,j)

        mat(2,2,2,i,j) = mat(2,2,2,i,j) - b2 * dtl(i,j) * Vh11(1,i,j)
        mat(2,2,3,i,j) = mat(2,2,3,i,j) - b2 * dtl(i,j) * Vh11(5,i,j)
        mat(2,3,2,i,j) = mat(2,3,2,i,j) - b2 * dtl(i,j) * Vh11(6,i,j)
        mat(2,3,3,i,j) = mat(2,3,3,i,j) - b2 * dtl(i,j) * Vh11(2,i,j)
        mat(2,4,4,i,j) = mat(2,4,4,i,j) - b2 * dtl(i,j) * Vh11(3,i,j)
        mat(2,5,5,i,j) = mat(2,5,5,i,j) - b2 * dtl(i,j) * Vh11(4,i,j)

        mat(3,2,2,i,j) = mat(3,2,2,i,j) - b3 * dtl(i,j) * Vh11(1,i,j)
        mat(3,2,3,i,j) = mat(3,2,3,i,j) - b3 * dtl(i,j) * Vh11(5,i,j)
        mat(3,3,2,i,j) = mat(3,3,2,i,j) - b3 * dtl(i,j) * Vh11(6,i,j)
        mat(3,3,3,i,j) = mat(3,3,3,i,j) - b3 * dtl(i,j) * Vh11(2,i,j)
        mat(3,4,4,i,j) = mat(3,4,4,i,j) - b3 * dtl(i,j) * Vh11(3,i,j)
        mat(3,5,5,i,j) = mat(3,5,5,i,j) - b3 * dtl(i,j) * Vh11(4,i,j)

        mat(4,2,2,i,j) = mat(4,2,2,i,j) - b4 * dtl(i,j) * Vh11(1,i,j)
        mat(4,2,3,i,j) = mat(4,2,3,i,j) - b4 * dtl(i,j) * Vh11(5,i,j)
        mat(4,3,2,i,j) = mat(4,3,2,i,j) - b4 * dtl(i,j) * Vh11(6,i,j)
        mat(4,3,3,i,j) = mat(4,3,3,i,j) - b4 * dtl(i,j) * Vh11(2,i,j)
        mat(4,4,4,i,j) = mat(4,4,4,i,j) - b4 * dtl(i,j) * Vh11(3,i,j)
        mat(4,5,5,i,j) = mat(4,5,5,i,j) - b4 * dtl(i,j) * Vh11(4,i,j)

        mat(5,2,2,i,j) = mat(5,2,2,i,j) - b5 * dtl(i,j) * Vh11(1,i,j)
        mat(5,2,3,i,j) = mat(5,2,3,i,j) - b5 * dtl(i,j) * Vh11(5,i,j)
        mat(5,3,2,i,j) = mat(5,3,2,i,j) - b5 * dtl(i,j) * Vh11(6,i,j)
        mat(5,3,3,i,j) = mat(5,3,3,i,j) - b5 * dtl(i,j) * Vh11(2,i,j)
        mat(5,4,4,i,j) = mat(5,4,4,i,j) - b5 * dtl(i,j) * Vh11(3,i,j)
        mat(5,5,5,i,j) = mat(5,5,5,i,j) - b5 * dtl(i,j) * Vh11(4,i,j)

!.... I term

        mat(3,1,1,i,j) = mat(3,1,1,i,j) + one
        mat(3,2,2,i,j) = mat(3,2,2,i,j) + one
        mat(3,3,3,i,j) = mat(3,3,3,i,j) + one
        mat(3,4,4,i,j) = mat(3,4,4,i,j) + one
        mat(3,5,5,i,j) = mat(3,5,5,i,j) + one

!.... implicit damping term updated to fourth order

        if (eps_e .ne. zero) then

        mat(1,1,1,i,j) = mat(1,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,2,2,i,j) = mat(1,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,3,3,i,j) = mat(1,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,4,4,i,j) = mat(1,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,5,5,i,j) = mat(1,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1

        mat(2,1,1,i,j) = mat(2,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,2,2,i,j) = mat(2,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,3,3,i,j) = mat(2,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,4,4,i,j) = mat(2,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,5,5,i,j) = mat(2,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2

        mat(3,1,1,i,j) = mat(3,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,2,2,i,j) = mat(3,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,3,3,i,j) = mat(3,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,4,4,i,j) = mat(3,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,5,5,i,j) = mat(3,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3

        mat(4,1,1,i,j) = mat(4,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,2,2,i,j) = mat(4,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,3,3,i,j) = mat(4,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,4,4,i,j) = mat(4,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,5,5,i,j) = mat(4,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4

        mat(5,1,1,i,j) = mat(5,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,2,2,i,j) = mat(5,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,3,3,i,j) = mat(5,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,4,4,i,j) = mat(5,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,5,5,i,j) = mat(5,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5

        end if

        end do loop_i
        end do loop_j

!.... apply boundary treatment to the LHS

        call lhsbt1f( mat, Ah, Dh, Vh11, dtl )

!.... apply boundary conditions to the LHS

        call lhsbc1f( mat, v, vm )

        return
        end

!=============================================================================!
        subroutine lhs2f( mat, Bh, Dh, Vh22, dtl, v)
!
!  Form the LHS for the \eta direction
!
!=============================================================================!
        use global
        use buffer
        use stencil
        implicit none

        real :: mat(5,ndof,ndof,nx,ny), v(ndof,nx,ny)
        real :: Bh(ndof,ndof,nx,ny), Dh(ndof,ndof,nx,ny), Vh22(6,nx,ny)
        real :: dtl(nx,ny)

!.... first derivative operator

        real :: a1, a2, a3, a4, a5
        real :: b1, b2, b3, b4, b5

        integer :: i,j, idof, jdof
!=============================================================================!

!.... fourth-order stencil

        a1 = alfa * ga1 / deta
        a2 = alfa * ga2 / deta
        a3 = alfa * zero / deta
        a4 = alfa * ga3 / deta
        a5 = alfa * ga4 / deta

        b1 = alfa * da1 / deta**2
        b2 = alfa * da2 / deta**2
        b3 = alfa * da3 / deta**2
        b4 = alfa * da4 / deta**2
        b5 = alfa * da5 / deta**2

!$omp parallel do private (i,idof,jdof)
loop_j: do j = 1, ny
loop_i: do i = 1, nx

!.... \hat{B} term

        do idof = 1, ndof
          do jdof = 1, ndof
            mat(1,idof,jdof,i,j) = a1 * dtl(i,j) * Bh(idof,jdof,i,j)
            mat(2,idof,jdof,i,j) = a2 * dtl(i,j) * Bh(idof,jdof,i,j)
            mat(3,idof,jdof,i,j) = a3 * dtl(i,j) * Bh(idof,jdof,i,j)
            mat(4,idof,jdof,i,j) = a4 * dtl(i,j) * Bh(idof,jdof,i,j)
            mat(5,idof,jdof,i,j) = a5 * dtl(i,j) * Bh(idof,jdof,i,j)
          end do
        end do

!.... \hat{V}_{\eta\eta} term

        mat(1,2,2,i,j) = mat(1,2,2,i,j) - b1 * dtl(i,j) * Vh22(1,i,j)
        mat(1,2,3,i,j) = mat(1,2,3,i,j) - b1 * dtl(i,j) * Vh22(5,i,j)
        mat(1,3,2,i,j) = mat(1,3,2,i,j) - b1 * dtl(i,j) * Vh22(6,i,j)
        mat(1,3,3,i,j) = mat(1,3,3,i,j) - b1 * dtl(i,j) * Vh22(2,i,j)
        mat(1,4,4,i,j) = mat(1,4,4,i,j) - b1 * dtl(i,j) * Vh22(3,i,j)
        mat(1,5,5,i,j) = mat(1,5,5,i,j) - b1 * dtl(i,j) * Vh22(4,i,j)

        mat(2,2,2,i,j) = mat(2,2,2,i,j) - b2 * dtl(i,j) * Vh22(1,i,j)
        mat(2,2,3,i,j) = mat(2,2,3,i,j) - b2 * dtl(i,j) * Vh22(5,i,j)
        mat(2,3,2,i,j) = mat(2,3,2,i,j) - b2 * dtl(i,j) * Vh22(6,i,j)
        mat(2,3,3,i,j) = mat(2,3,3,i,j) - b2 * dtl(i,j) * Vh22(2,i,j)
        mat(2,4,4,i,j) = mat(2,4,4,i,j) - b2 * dtl(i,j) * Vh22(3,i,j)
        mat(2,5,5,i,j) = mat(2,5,5,i,j) - b2 * dtl(i,j) * Vh22(4,i,j)

        mat(3,2,2,i,j) = mat(3,2,2,i,j) - b3 * dtl(i,j) * Vh22(1,i,j)
        mat(3,2,3,i,j) = mat(3,2,3,i,j) - b3 * dtl(i,j) * Vh22(5,i,j)
        mat(3,3,2,i,j) = mat(3,3,2,i,j) - b3 * dtl(i,j) * Vh22(6,i,j)
        mat(3,3,3,i,j) = mat(3,3,3,i,j) - b3 * dtl(i,j) * Vh22(2,i,j)
        mat(3,4,4,i,j) = mat(3,4,4,i,j) - b3 * dtl(i,j) * Vh22(3,i,j)
        mat(3,5,5,i,j) = mat(3,5,5,i,j) - b3 * dtl(i,j) * Vh22(4,i,j)

        mat(4,2,2,i,j) = mat(4,2,2,i,j) - b4 * dtl(i,j) * Vh22(1,i,j)
        mat(4,2,3,i,j) = mat(4,2,3,i,j) - b4 * dtl(i,j) * Vh22(5,i,j)
        mat(4,3,2,i,j) = mat(4,3,2,i,j) - b4 * dtl(i,j) * Vh22(6,i,j)
        mat(4,3,3,i,j) = mat(4,3,3,i,j) - b4 * dtl(i,j) * Vh22(2,i,j)
        mat(4,4,4,i,j) = mat(4,4,4,i,j) - b4 * dtl(i,j) * Vh22(3,i,j)
        mat(4,5,5,i,j) = mat(4,5,5,i,j) - b4 * dtl(i,j) * Vh22(4,i,j)

        mat(5,2,2,i,j) = mat(5,2,2,i,j) - b5 * dtl(i,j) * Vh22(1,i,j)
        mat(5,2,3,i,j) = mat(5,2,3,i,j) - b5 * dtl(i,j) * Vh22(5,i,j)
        mat(5,3,2,i,j) = mat(5,3,2,i,j) - b5 * dtl(i,j) * Vh22(6,i,j)
        mat(5,3,3,i,j) = mat(5,3,3,i,j) - b5 * dtl(i,j) * Vh22(2,i,j)
        mat(5,4,4,i,j) = mat(5,4,4,i,j) - b5 * dtl(i,j) * Vh22(3,i,j)
        mat(5,5,5,i,j) = mat(5,5,5,i,j) - b5 * dtl(i,j) * Vh22(4,i,j)

!.... I term

        mat(3,1,1,i,j) = mat(3,1,1,i,j) + one
        mat(3,2,2,i,j) = mat(3,2,2,i,j) + one
        mat(3,3,3,i,j) = mat(3,3,3,i,j) + one
        mat(3,4,4,i,j) = mat(3,4,4,i,j) + one
        mat(3,5,5,i,j) = mat(3,5,5,i,j) + one

!.... implicit damping term updated to fourth-order

        if (eps_e .ne. zero) then

        mat(1,1,1,i,j) = mat(1,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,2,2,i,j) = mat(1,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,3,3,i,j) = mat(1,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,4,4,i,j) = mat(1,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1
        mat(1,5,5,i,j) = mat(1,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb1

        mat(2,1,1,i,j) = mat(2,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,2,2,i,j) = mat(2,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,3,3,i,j) = mat(2,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,4,4,i,j) = mat(2,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2
        mat(2,5,5,i,j) = mat(2,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb2

        mat(3,1,1,i,j) = mat(3,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,2,2,i,j) = mat(3,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,3,3,i,j) = mat(3,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,4,4,i,j) = mat(3,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3
        mat(3,5,5,i,j) = mat(3,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb3

        mat(4,1,1,i,j) = mat(4,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,2,2,i,j) = mat(4,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,3,3,i,j) = mat(4,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,4,4,i,j) = mat(4,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4
        mat(4,5,5,i,j) = mat(4,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb4

        mat(5,1,1,i,j) = mat(5,1,1,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,2,2,i,j) = mat(5,2,2,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,3,3,i,j) = mat(5,3,3,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,4,4,i,j) = mat(5,4,4,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5
        mat(5,5,5,i,j) = mat(5,5,5,i,j) + alfa * dtl(i,j) * eps_e * buff(i,j) * fb5

        end if

!.... sponge

        if (ispg .eq. 1) then

          mat(3,1,1,i,j) = mat(3,1,1,i,j) + alfa * dtl(i,j) * spg(i,j)
          mat(3,2,2,i,j) = mat(3,2,2,i,j) + alfa * dtl(i,j) * spg(i,j)
          mat(3,3,3,i,j) = mat(3,3,3,i,j) + alfa * dtl(i,j) * spg(i,j)
          mat(3,4,4,i,j) = mat(3,4,4,i,j) + alfa * dtl(i,j) * spg(i,j)
          mat(3,5,5,i,j) = mat(3,5,5,i,j) + alfa * dtl(i,j) * spg(i,j)

        else if (ispg .ge. 2) then

          mat(3,1,1,i,j) = mat(3,1,1,i,j) + alfa * dtl(i,j) * ( spg(i,j) + spg2(i,j) )
          mat(3,2,2,i,j) = mat(3,2,2,i,j) + alfa * dtl(i,j) * ( spg(i,j) + spg2(i,j) )
          mat(3,3,3,i,j) = mat(3,3,3,i,j) + alfa * dtl(i,j) * ( spg(i,j) + spg2(i,j) )
          mat(3,4,4,i,j) = mat(3,4,4,i,j) + alfa * dtl(i,j) * ( spg(i,j) + spg2(i,j) )
          mat(3,5,5,i,j) = mat(3,5,5,i,j) + alfa * dtl(i,j) * ( spg(i,j) + spg2(i,j) )

        end if

        end do loop_i
        end do loop_j

!.... apply boundary treatment to the LHS

        call lhsbt2f( mat, Bh, Dh, Vh22, spg, spg2, dtl )

!.... apply boundary conditions to the LHS

        call lhsbc2f( mat, v )

        return
        end


text: #redirect Users/JoSealwood

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details
      SUBROUTINE MXLIST (UNIT, OFFSET,
     *   DICT, DPOINT, LDICT, NNAMES,
     *   VOID, LVOID, NVOIDS, CHRCOL)

      IMPLICIT INTEGER (A-Z)

C     This routine lists the contents of the tables of the
C     memory manager.

C***********************************************************************

C     UNIT     Output unit number
C     OFFSET   Offset to internal reference vector
C     DICT     Dictionary table
C     DPOINT   Pointer table
C     LDICT    Dimension of dictionary
C     NNAMES   Number of names in the dictionary
               CHARACTER*8 DICT(LDICT,CHRCOL)
               DIMENSION DPOINT(LDICT,CHRCOL,3), NNAMES(2)
C     VOID     Void table
C     LVOID    Dimension of void table
C     NVOIDS   Number of voids
               DIMENSION VOID(LVOID,CHRCOL,2), NVOIDS(2)
C     CHRCOL   Column number for character tables

C***********************************************************************

      CHARACTER*8 TNAME1, TNAME2

      TOFF = OFFSET
      DO 290 ICOL = 1, CHRCOL
         IF (ICOL .EQ. 2) THEN
            TOFF = 0
            WRITE (UNIT, 10090)
         END IF
C        DICTIONARY.

         MNDICT = 0
         MCDICT = 0
         DO 100 I = 1, NNAMES(ICOL)
            MNDICT = MNDICT + MAX(0,DPOINT(I,ICOL,2))
            MCDICT = MCDICT + MAX(0,DPOINT(I,ICOL,3))
  100    CONTINUE
         WRITE (UNIT,10000)
         DO 110 I = 1, NNAMES(ICOL)
            IF (DPOINT(I,ICOL,2) .GE. 0) THEN
               WRITE (UNIT,10010) I, DICT(I,ICOL),
     *            DPOINT(I,ICOL,1)+TOFF,
     *            DPOINT(I,ICOL,2), DPOINT(I,ICOL,3)
            ELSE
               WRITE (UNIT,10020) I, DICT(I,ICOL),
     *            DPOINT(I,ICOL,2), DPOINT(I,ICOL,3)
            END IF
  110    CONTINUE
         WRITE (UNIT,10040) MNDICT, MCDICT

C        VOID TABLE.

         MNVOID = 0
         DO 120 I = 1, NVOIDS(ICOL)
            MNVOID = MNVOID + VOID(I,ICOL,2)
  120    CONTINUE
         WRITE (UNIT,10070)
         WRITE (UNIT,10060) (I,VOID(I,ICOL,1)+TOFF,
     *      VOID(I,ICOL,2),I=1,NVOIDS(ICOL))
         WRITE (UNIT,10050) MNVOID

C        OUTPUT ORDERED LIST OF TABLES.

C        First sort dictionary into location order with unresolved
C        allocations first.

         JSTRT = 2
         DO 150 I = 1, NNAMES(ICOL)-1
            IF (DPOINT(I,ICOL,2) .GE. 0) THEN
               DO 130 J = MAX(JSTRT,I+1), NNAMES(ICOL)
                  IF (DPOINT(J,ICOL,2) .LT. 0) THEN
                     JSTRT = J+1
                     TNAME1 = DICT(I,ICOL)
                     TNAME2 = DICT(J,ICOL)
                     DICT(I,ICOL) = TNAME2
                     DICT(J,ICOL) = TNAME1
                     TEMP = DPOINT(I,ICOL,1)
                     DPOINT(I,ICOL,1) = DPOINT(J,ICOL,1)
                     DPOINT(J,ICOL,1) = TEMP
                     TEMP = DPOINT(I,ICOL,2)
                     DPOINT(I,ICOL,2) = DPOINT(J,ICOL,2)
                     DPOINT(J,ICOL,2) = TEMP
                     TEMP = DPOINT(I,ICOL,3)
                     DPOINT(I,ICOL,3) = DPOINT(J,ICOL,3)
                     DPOINT(J,ICOL,3) = TEMP
                     GO TO 140
                  END IF
  130          CONTINUE
               GO TO 160
            END IF
  140       CONTINUE
  150    CONTINUE
  160    CONTINUE
         DO 170 NDEFER = 1, NNAMES(ICOL)
            IF (DPOINT(NDEFER,ICOL,2) .GE. 0) GO TO 180
  170    CONTINUE
  180    NDEFER = NDEFER - 1
         DO 200 I = NDEFER+1, NNAMES(ICOL)-1
            DO 190 J = I+1, NNAMES(ICOL)
               IF (DPOINT(J,ICOL,1) .LT. DPOINT(I,ICOL,1)) THEN
                  TNAME1 = DICT(I,ICOL)
                  TNAME2 = DICT(J,ICOL)
                  DICT(I,ICOL) = TNAME2
                  DICT(J,ICOL) = TNAME1
                  TEMP = DPOINT(I,ICOL,1)
                  DPOINT(I,ICOL,1) = DPOINT(J,ICOL,1)
                  DPOINT(J,ICOL,1) = TEMP
                  TEMP = DPOINT(I,ICOL,2)
                  DPOINT(I,ICOL,2) = DPOINT(J,ICOL,2)
                  DPOINT(J,ICOL,2) = TEMP
                  TEMP = DPOINT(I,ICOL,3)
                  DPOINT(I,ICOL,3) = DPOINT(J,ICOL,3)
                  DPOINT(J,ICOL,3) = TEMP
               END IF
  190       CONTINUE
  200    CONTINUE

C        STARTING STUFF FOR LOOP

         DO 210 IDICT = 1, NNAMES(ICOL)
            IF (DPOINT(IDICT,ICOL,2) .GE. 0) GO TO 220
  210    CONTINUE
  220    CONTINUE
         IF (IDICT .LE. NNAMES(ICOL) .AND. NVOIDS(ICOL) .GT. 0) THEN
            NXTLOC = MIN(DPOINT(IDICT,ICOL,1), VOID(1,ICOL,1))
         ELSE IF (IDICT .LE. NNAMES(ICOL)) THEN
            NXTLOC = DPOINT(IDICT,ICOL,1)
         ELSE IF (NVOIDS(ICOL) .GT. 0) THEN
            NXTLOC = VOID(1,ICOL,1)
         ELSE
            NXTLOC = 0
         END IF
         IVOID = 1

         WRITE (UNIT, 10080)
         ILIST = 0

C        Deferred space names first.

         DTOT = 0
         DCTOT = 0
         DO 230 IDICT = 1, NNAMES(ICOL)
            IF (DPOINT(IDICT,ICOL,2) .LT. 0) THEN
               ILIST = ILIST + 1
               WRITE (UNIT,10020) ILIST, DICT(IDICT,ICOL),
     *            DPOINT(IDICT,ICOL,2),
     *            DPOINT(IDICT,ICOL,3)
               DTOT = DTOT - DPOINT(IDICT,ICOL,2)
               DCTOT = DCTOT + MAX (0, DPOINT(IDICT,ICOL,3))
            ELSE
               GO TO 240
            END IF
  230    CONTINUE
  240    CONTINUE
         IF (IDICT .GT. 1) WRITE (UNIT, 10030)
     *      'DEFERRED TOTAL', DTOT, DCTOT

C        LOOP

         TOTAL = 0
         SUBTOT = 0
         CSTOT = 0
         CTOT = 0
  250    CONTINUE
            IF (IDICT .LE. NNAMES(ICOL)
     *         .AND. IVOID .LE. NVOIDS(ICOL)) THEN
               IF (NXTLOC .EQ. DPOINT(IDICT,ICOL,1)) THEN
                  ILIST = ILIST + 1
                  WRITE (UNIT,10010) ILIST, DICT(IDICT,ICOL),
     *               DPOINT(IDICT,ICOL,1)+TOFF, DPOINT(IDICT,ICOL,2),
     *               DPOINT(IDICT,ICOL,3)
                  NXTLOC = NXTLOC + DPOINT(IDICT,ICOL,2)
                  SUBTOT = SUBTOT + DPOINT(IDICT,ICOL,2)
                  TOTAL = TOTAL + DPOINT(IDICT,ICOL,2)
                  CSTOT = CSTOT + MAX(0,DPOINT(IDICT,ICOL,3))
                  CTOT = CTOT + MAX(0,DPOINT(IDICT,ICOL,3))
                  IDICT = IDICT + 1
               ELSE IF (NXTLOC .EQ. VOID(IVOID,ICOL,1)) THEN
                  ILIST = ILIST + 1
                  WRITE (UNIT,10010) ILIST, ' ',
     *               VOID(IVOID,ICOL,1)+TOFF,
     *               VOID(IVOID,ICOL,2)
                  NXTLOC = NXTLOC + VOID(IVOID,ICOL,2)
                  SUBTOT = SUBTOT + VOID(IVOID,ICOL,2)
                  TOTAL = TOTAL + VOID(IVOID,ICOL,2)
                  IVOID = IVOID + 1
               ELSE
                  NXTLOC = MIN( DPOINT(IDICT,ICOL,1),
     *               VOID(IVOID,ICOL,1) )
                  WRITE (UNIT, 10030) 'BLOCK SIZE', SUBTOT, CSTOT
                  SUBTOT = 0
               END IF
            ELSE IF (IDICT .LE. NNAMES(ICOL)) THEN
               IF (NXTLOC .EQ. DPOINT(IDICT,ICOL,1)) THEN
                  ILIST = ILIST + 1
                  WRITE (UNIT,10010) ILIST, DICT(IDICT,ICOL),
     *               DPOINT(IDICT,ICOL,1)+TOFF, DPOINT(IDICT,ICOL,2),
     *               DPOINT(IDICT,ICOL,3)
                  NXTLOC = NXTLOC + DPOINT(IDICT,ICOL,2)
                  SUBTOT = SUBTOT + DPOINT(IDICT,ICOL,2)
                  TOTAL = TOTAL + DPOINT(IDICT,ICOL,2)
                  CSTOT = CSTOT + MAX(0,DPOINT(IDICT,ICOL,3))
                  CTOT = CTOT + MAX(0,DPOINT(IDICT,ICOL,3))
                  IDICT = IDICT + 1
               ELSE
                  NXTLOC = DPOINT(IDICT,ICOL,1)
                  WRITE (UNIT, 10030) 'BLOCK SIZE', SUBTOT, CSTOT
                  SUBTOT = 0
               END IF
            ELSE IF (IVOID .LE. NVOIDS(ICOL)) THEN
               IF (NXTLOC .EQ. VOID(IVOID,ICOL,1)) THEN
                  ILIST = ILIST + 1
                  WRITE (UNIT,10010) ILIST, ' ',
     *               VOID(IVOID,ICOL,1)+TOFF,
     *               VOID(IVOID,ICOL,2)
                  NXTLOC = NXTLOC + VOID(IVOID,ICOL,2)
                  SUBTOT = SUBTOT + VOID(IVOID,ICOL,2)
                  TOTAL = TOTAL + VOID(IVOID,ICOL,2)
                  IVOID = IVOID + 1
               ELSE
                  NXTLOC = VOID(IVOID,ICOL,1)
                  WRITE (UNIT, 10030) 'BLOCK SIZE', SUBTOT, CSTOT
                  SUBTOT = 0
               END IF
            ELSE
               GO TO 260
            END IF
         GO TO 250
  260    CONTINUE
         WRITE (UNIT, 10030) 'BLOCK SIZE', SUBTOT, CSTOT
         WRITE (UNIT, 10030) 'ALLOCATED TOTAL', TOTAL, CTOT
         WRITE (UNIT, 10030) '    GRAND TOTAL', TOTAL+DTOT, CTOT+DCTOT

C        SORT DICTIONARY BACK INTO NAME ORDER

         DO 280 I = 1, NNAMES(ICOL)-1
            DO 270 J = I+1, NNAMES(ICOL)
               IF (DICT(J,ICOL) .LT. DICT(I,ICOL)) THEN
                  TNAME1 = DICT(I,ICOL)
                  TNAME2 = DICT(J,ICOL)
                  DICT(I,ICOL) = TNAME2
                  DICT(J,ICOL) = TNAME1
                  TEMP = DPOINT(I,ICOL,1)
                  DPOINT(I,ICOL,1) = DPOINT(J,ICOL,1)
                  DPOINT(J,ICOL,1) = TEMP
                  TEMP = DPOINT(I,ICOL,2)
                  DPOINT(I,ICOL,2) = DPOINT(J,ICOL,2)
                  DPOINT(J,ICOL,2) = TEMP
                  TEMP = DPOINT(I,ICOL,3)
                  DPOINT(I,ICOL,3) = DPOINT(J,ICOL,3)
                  DPOINT(J,ICOL,3) = TEMP
               END IF
  270       CONTINUE
  280    CONTINUE
  290 CONTINUE
      RETURN
10000 FORMAT(//1X,50('*')/1X,7('* '),
     *   ' D I C T I O N A R Y  ',7('* ')//
     *   T45, 'NUMERIC', T61, 'CHARACTER'/
     *   T10,'NAME      LOCATION',t45, 'LENGTH', T61,'LENGTH'/)
10010 FORMAT((1X,I4,2X,A8,3(2X,I16)))
10020 FORMAT((1X,I4,2X,A8,12X,2(2X,I16)))
10030 FORMAT (1X, A, T36, I16, 2X, I16/)
10040 FORMAT(/,T10,'TOTAL',T36,I16, 2X, I16)
10050 FORMAT(/,T10,'TOTAL',T24,I16)
10060 FORMAT((1X,I6,2(3X,I16)))
10070 FORMAT(///' * * * V O I D   T A B L E * * *'//
     *   '           LOCATION       LENGTH'/)
10080 FORMAT(//1X,50('*')/1X,6('* '),
     *   ' O R D E R E D   L I S T  ',6('* ')//
     *   ,T45, 'NUMERIC', T61, 'CHARACTER'/
     *   T10,'NAME      LOCATION',t45, 'LENGTH', T61,'LENGTH'/)
10090 FORMAT(//1X,79('*')//,' SPLIT CHARACTER STORAGE'//1X,79('*')//)
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!>  ChiDG Interpolation Utilities
!!
!!  Purpose:
!!  ---------------------------------------
!!  Facilitate interpolation from modal representations
!!  to node discrete representations of variables.
!!
!!
!!  Public routines:
!!  ---------------------------------------
!!      interpolate_element_autodiff
!!      interpolate_face_autodiff
!!
!!      interpolate_element_standard
!!      interpolate_face_standard
!!
!!
!!  Module utility routines:
!!  ---------------------------------------
!!      get_face_interpolation_info
!!      get_face_interpolation_interpolator
!!      get_face_interpolation_mask
!!      get_face_interpolation_comm
!!      get_face_interpolation_ndonors
!!      get_face_interpolation_style
!!      get_interpolation_nderiv
!!
!!
!---------------------------------------------------------------------------------------------
module mod_interpolate
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: CHIMERA, INTERIOR, BOUNDARY, &
                                      ME, NEIGHBOR, ONE, ZERO

    use mod_chidg_mpi,          only: IRANK
    use mod_DNAD_tools,         only: compute_neighbor_face
    use DNAD_D

    use type_mesh,              only: mesh_t
    use type_element_info,      only: element_info_t
    use type_face_info,         only: face_info_t
    use type_function_info,     only: function_info_t
    use type_recv,              only: recv_t
    use type_chidg_vector,      only: chidg_vector_t
    implicit none

contains

    !>  Interpolate polynomial expansion to volume quadrature node set, initializing the
    !!  automatic differentiation process if needed.
    !!
    !!  This returns an array of automatic differentiation(AD) values at the quadrature nodes,
    !!  that also have their derivatives initialized.
    !!
    !!  Some interpolation parameters to note that are used inside here:
    !!      - interpolation_type:   'value', 'grad1', 'grad2', 'grad3'
    !!
    !!  Partial mode expansions can be used to perform interpolations. To do this, pass the
    !!  Pmin and Pmax optional arguments. Pmin is the minimum 1D polynomial order to be used
    !!  in the interpolation. Pmax is the maximum 1D polynomial order to be used in the
    !!  interpolation. One might pass Pmin=0, Pmax=1. This would return an interpolation that
    !!  used only the first 8 modes in the polynomial expansion. The constant mode and all
    !!  the linear modes, but no modes of higher order.
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   11/5/2016
    !!
    !!  @param[in]  itime - Index for the time step in solution
    !!
    !-----------------------------------------------------------------------------------------
    function interpolate_element_autodiff(mesh,vector,elem_info,fcn_info,ifield,itime,interpolation_type,mode_min,mode_max) result(var_gq)
        type(mesh_t),           intent(in)              :: mesh
        type(chidg_vector_t),   intent(in)              :: vector
        type(element_info_t),   intent(in)              :: elem_info
        type(function_info_t),  intent(in)              :: fcn_info
        integer(ik),            intent(in)              :: ifield
        integer(ik),            intent(in)              :: itime
        character(*),           intent(in)              :: interpolation_type
        integer(ik),            intent(in), optional    :: mode_min
        integer(ik),            intent(in), optional    :: mode_max

        character(:),   allocatable :: user_msg

        type(AD_D), allocatable :: qdiff(:), var_gq(:)
        real(rk),   allocatable :: qtmp(:)

        integer(ik) :: nderiv, set_deriv, iterm, nterms, ierr
        logical     :: differentiate_me

        associate( idom => elem_info%idomain_l, ielem => elem_info%ielement_l )

        !
        ! Get number of derivatives to initialize for automatic differentiation
        !
        nderiv = get_interpolation_nderiv(mesh,fcn_info)

        !
        ! Determine nterms from interpolator matrix
        !
        nterms = mesh%domain(idom)%elems(ielem)%basis_s%nterms_i()

        ! Allocate qdiff buffer and initialize derivative allocations + zero
        allocate(qdiff(nterms), stat=ierr)
        if (ierr /= 0) call AllocationError
        do iterm = 1,nterms
            qdiff(iterm) = AD_D(nderiv)
        end do
        qdiff = ZERO

        !
        ! Retrieve modal coefficients representing ifield in vector to 'qtmp'
        !
        qtmp = vector%dom(idom)%vecs(ielem)%getvar(ifield,itime)

        !
        ! 1: If mode_min,mode_max are present, use a subset of the expansion modes.
        ! 2: If not present, use the entire expansion.
        !
        if (present(mode_min) .and. present(mode_max)) then
            qdiff(mode_min:mode_max) = qtmp(mode_min:mode_max)
        else
            qdiff(1:nterms) = qtmp(1:nterms)
        end if

        !
        ! If the current element is being differentiated (ielem == ielem_seed)
        ! then copy the solution modes to local AD variable and seed derivatives
        !
        differentiate_me = ( (elem_info%idomain_g  == fcn_info%seed%idomain_g ) .and. &
                             (elem_info%ielement_g == fcn_info%seed%ielement_g) )

        if (differentiate_me) then
            ! Differentiating element, initialize appropriate derivatives to ONE
            do iterm = 1,size(qdiff)
                set_deriv = (ifield - 1)*mesh%domain(idom)%elems(ielem)%nterms_s + iterm
                qdiff(iterm)%xp_ad_(set_deriv) = ONE
            end do
        else
            ! Not differentiating element, set all derivatives to ZERO
            do iterm = 1,size(qdiff)
                qdiff(iterm)%xp_ad_ = ZERO
            end do
        end if

        !
        ! Select interpolation type
        !
        select case (interpolation_type)
            case('value')
                var_gq = matmul(mesh%domain(idom)%elems(ielem)%basis_s%interpolator('Value'),qdiff)

            case('grad1')
                var_gq = matmul(mesh%domain(idom)%elems(ielem)%grad1,qdiff)

            case('grad2')
                var_gq = matmul(mesh%domain(idom)%elems(ielem)%grad2,qdiff)

            case('grad3')
                var_gq = matmul(mesh%domain(idom)%elems(ielem)%grad3,qdiff)

            case default
                user_msg = "interpolate_element_autodiff: The 'interpolation_type' incoming&
                            parameter was not a valid string. Valid strings for &
                            'interpolation_type' include 'value', 'grad1', 'grad2', 'grad3'."
                call chidg_signal_one(FATAL,user_msg,interpolation_type)
        end select

        end associate

    end function interpolate_element_autodiff
    !*****************************************************************************************

    !>  Interpolate variable from polynomial expansion to explicit values at quadrature
    !!  nodes. The automatic differentiation process really starts here, when the polynomial
    !!  expansion is evaluated.
    !!
    !!  The interpolation process occurs through a matrix-vector multiplication. That is an
    !!  interpolation matrix multiplied by a vector of modes from the polynomial expansion.
    !!  To start the automatic differentiation, the derivative arrays of the values for the
    !!  polynomial modes must be initialized before any computation. So, before the
    !!  matrix-vector multiplication.
    !!
    !!  Some interpolation parameters to note that a user might select:
    !!      - interpolation_type:   'value', 'grad1', 'grad2', 'grad3'
    !!      - interpolation_source: ME, NEIGHBOR
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]      mesh                    Array of mesh instances.
    !!  @param[in]      face                    Face indices for locating the face in mesh.
    !!  @param[in]      vector                  A chidg_vector containing a modal representation of fields on mesh
    !!  @param[in]      ifield                  Index of field being interpolated
    !!  @param[inout]   var_gq                  Autodiff values of field evaluated at gq points
    !!  @param[in]      interpolation_type      Interpolate 'value', 'grad1', 'grad2', 'grad3'
    !!  @param[in]      interpolation_source    ME/NEIGHBOR indicating element to interpolate from
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   11/5/2016
    !!
    !!  @param[in]      itime                   Index for time step in solution
    !!
    !------------------------------------------------------------------------------------------
    function interpolate_face_autodiff(mesh,vector,face_info,fcn_info, ifield, itime, interpolation_type, interpolation_source) result(var_gq)
        type(mesh_t),           intent(in)              :: mesh
        type(chidg_vector_t),   intent(in)              :: vector
        type(face_info_t),      intent(in)              :: face_info
        type(function_info_t),  intent(in)              :: fcn_info
        integer(ik),            intent(in)              :: ifield
        integer(ik),            intent(in)              :: itime
        character(*),           intent(in)              :: interpolation_type
        integer(ik),            intent(in)              :: interpolation_source

        type(face_info_t)   :: iface_info
        type(recv_t)        :: recv_info

        type(AD_D),         allocatable :: qdiff(:), var_gq(:)
        real(rk),           allocatable :: qtmp(:)
        real(rk),           allocatable :: interpolator(:,:)
        character(:),       allocatable :: interpolation_style

        integer(ik) :: nderiv, set_deriv, iterm, igq, nterms_s, ierr, nnodes
        logical     :: differentiate_me, conforming_interpolation, chimera_interpolation, parallel_interpolation

        ! Chimera data
        integer(ik)                 :: ndonors, idonor
        logical,    allocatable     :: mask(:)          ! node mask for distributing Chimera quadrature points
        type(AD_D), allocatable     :: var_gq_chimera(:)

        !
        ! Allocate output array
        !
        nnodes   = mesh%domain(face_info%idomain_l)%elems(face_info%ielement_l)%basis_s%nnodes_if()
        nterms_s = mesh%domain(face_info%idomain_l)%elems(face_info%ielement_l)%basis_s%nterms_i()
        allocate(var_gq(nnodes), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Get number of donors for the interpolation
        !
        ndonors = get_face_interpolation_ndonors(mesh,face_info,interpolation_source)

        !
        ! Get interpolation style. Conforming or Chimera
        !
        interpolation_style = get_face_interpolation_style(mesh,face_info,interpolation_source)
        conforming_interpolation = (interpolation_style == 'conforming')
        chimera_interpolation    = (interpolation_style == 'chimera')

        !
        ! Get number of derivatives to initialize for automatic differentiation
        !
        nderiv = get_interpolation_nderiv(mesh,fcn_info)

        !
        ! For each donor element to the interpolation.
        ! (ndonors could be > 1 for Chimera interpolations)
        !
        do idonor = 1,ndonors

            !
            ! Get face info for face being interpolated to(ME, NEIGHBOR),
            ! interpolation matrix, and recv data for parallel access
            !
            iface_info   = get_face_interpolation_info(        mesh,face_info,interpolation_source,idonor)
            mask         = get_face_interpolation_mask(        mesh,face_info,interpolation_source,idonor)
            recv_info    = get_face_interpolation_comm(        mesh,face_info,interpolation_source,idonor)
            interpolator = get_face_interpolation_interpolator(mesh,face_info,interpolation_source,idonor,interpolation_type,iface_info)

            parallel_interpolation = (recv_info%comm /= 0)

            !
            ! Deduce nterms from interpolator matrix
            !
            nterms_s = size(interpolator,2)

            !
            ! Allocate solution and derivative arrays for temporary solution variable
            !
            if ( allocated(qdiff) ) deallocate(qdiff)
            allocate(qdiff(nterms_s), stat=ierr)
            if (ierr /= 0) call AllocationError

            do iterm = 1,nterms_s
                qdiff(iterm) = AD_D(nderiv)
            end do

            !
            ! Retrieve modal coefficients for ifield from vector
            !
            if (parallel_interpolation) then
                qtmp = vector%recv%comm(recv_info%comm)%dom(recv_info%domain)%vecs(recv_info%element)%getvar(ifield,itime)
            else
                qtmp = vector%dom(iface_info%idomain_l)%vecs(iface_info%ielement_l)%getvar(ifield,itime)
            end if

            !
            ! Copy correct number of modes. We use this because there is the possibility
            ! that 'q' could be an auxiliary vector with nterms > nterms_s. For example,
            ! if wall distance was computed for P1, and we are running Navier Stokes P0.
            ! So we take only up to the modes that we need.
            !
            qdiff(1:nterms_s) = qtmp(1:nterms_s)

            !
            ! If the current element is being differentiated (ielem == ielem_seed)
            ! then copy the solution modes to local AD variable and seed derivatives
            !
            differentiate_me = ( (iface_info%idomain_g  == fcn_info%seed%idomain_g ) .and. &
                                 (iface_info%ielement_g == fcn_info%seed%ielement_g) )

            if ( differentiate_me ) then
                ! Loop through the terms in qdiff, seed appropriate derivatives to ONE
                do iterm = 1,size(qdiff)
                    ! For the given term, seed its appropriate derivative
                    set_deriv = (ifield - 1)*nterms_s + iterm
                    qdiff(iterm)%xp_ad_(set_deriv) = ONE
                end do

            else
                ! Loop through the terms in qdiff. Set all derivatives to ZERO
                do iterm = 1,size(qdiff)
                    qdiff(iterm)%xp_ad_ = ZERO
                end do

            end if

            !
            ! Interpolate solution to GQ nodes via matrix-vector multiplication
            !
            if ( conforming_interpolation ) then
                var_gq = matmul(interpolator,  qdiff)

            elseif ( chimera_interpolation ) then
                ! Perform interpolation
                var_gq_chimera = matmul(interpolator,  qdiff)

                ! Scatter chimera nodes to appropriate location in var_gq
                var_gq = unpack(var_gq_chimera,mask,var_gq)

            else
                call chidg_signal(FATAL,"interpolate_face: face interpolation type error")
            end if

        end do ! idonor

    end function interpolate_face_autodiff
    !*****************************************************************************************

    !> Compute variable at quadrature nodes.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   11/5/2016
    !!
    !----------------------------------------------------------------------------------------
    function interpolate_element_standard(mesh,q,idomain_l,ielement_l,ifield,itime,interpolation_type) result(var_gq)
        type(mesh_t),           intent(in)      :: mesh
        type(chidg_vector_t),   intent(in)      :: q
        integer(ik),            intent(in)      :: idomain_l
        integer(ik),            intent(in)      :: ielement_l
        integer(ik),            intent(in)      :: ifield
        integer(ik),            intent(in)      :: itime
        character(*),           intent(in)      :: interpolation_type

        real(rk),   allocatable :: var_gq(:)

        !
        ! Use quadrature instance to compute variable at quadrature nodes.
        ! This takes the form of a matrix multiplication of the quadrature matrix
        ! with the array of modes for the given variable.
        !
        select case (interpolation_type)
            case('value')
                var_gq = matmul(mesh%domain(idomain_l)%elems(ielement_l)%basis_s%interpolator('Value'), q%dom(idomain_l)%vecs(ielement_l)%getvar(ifield,itime))
            case('grad1')
                var_gq = matmul(mesh%domain(idomain_l)%elems(ielement_l)%grad1,      q%dom(idomain_l)%vecs(ielement_l)%getvar(ifield,itime))
            case('grad2')
                var_gq = matmul(mesh%domain(idomain_l)%elems(ielement_l)%grad2,      q%dom(idomain_l)%vecs(ielement_l)%getvar(ifield,itime))
            case('grad3')
                var_gq = matmul(mesh%domain(idomain_l)%elems(ielement_l)%grad3,      q%dom(idomain_l)%vecs(ielement_l)%getvar(ifield,itime))
            case default
                call chidg_signal(FATAL,"interpolate_element_standard: invalid interpolation_type. Options are 'value', 'grad1', 'grad2', 'grad3'.")
        end select

    end function interpolate_element_standard
    !*****************************************************************************************

    !> Compute variable at quadrature nodes.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   11/5/2016
    !!
    !-----------------------------------------------------------------------------------------
    function interpolate_face_standard(mesh,q,idomain_l,ielement_l,iface,ifield,itime) result(var_gq)
        type(mesh_t),           intent(in)      :: mesh
        type(chidg_vector_t),   intent(in)      :: q
        integer(ik),            intent(in)      :: idomain_l, ielement_l, iface, ifield
        integer(ik),            intent(in)      :: itime

        real(rk),   allocatable :: var_gq(:)

        !
        ! Use quadrature instance to compute variable at quadrature nodes.
        ! This takes the form of a matrix multiplication of the face quadrature matrix
        ! with the array of modes for the given variable
        !
        var_gq = matmul(mesh%domain(idomain_l)%faces(ielement_l,iface)%basis_s%interpolator('Value',iface), q%dom(idomain_l)%vecs(ielement_l)%getvar(ifield,itime))

    end function interpolate_face_standard
    !*****************************************************************************************

!
!    !> Interpolate variable from polynomial expansion to explicit values at quadrature nodes. The automatic
!    !! differentiation process really starts here, when the polynomial expansion is evaluated.
!    !!
!    !! The interpolation process occurs through a matrix-vector multiplication. That is an interpolation matrix
!    !! multiplied by a vector of modes from the polynomial expansion. To start the automatic differentiation,
!    !! the derivative arrays of the values for the polynomial modes must be initialized before any computation.
!    !! So, before the matrix-vector multiplication.
!    !!
!    !!  @author Nathan A. Wukie
!    !!  @date   2/1/2016
!    !!
!    !!  @param[in]      mesh        Array of mesh instances.
!    !!  @param[in]      face        Face info, such as indices for locating the face in the mesh.
!    !!  @param[in]      q           Solution vector
!    !!  @param[in]      ieqn        Index of the equation variable being interpolated
!    !!  @param[inout]   var_gq      Array of auto-diff values of the equation evaluated at gq points that is passed back
!    !!  @param[in]      source      ME/NEIGHBOR indicating which element to interpolate from
!    !!
!    !-----------------------------------------------------------------------------------------------------------
!    subroutine interpolate_boundary_autodiff(mesh,face,q,ieqn,points,var)
!        type(mesh_t),           intent(in)              :: mesh(:)
!        type(face_info_t),      intent(in)              :: face
!        type(chidg_vector_t),    intent(in)              :: q
!        integer(ik),            intent(in)              :: ieqn
!        type(point_t),          intent(in)              :: points(:)
!        type(AD_D),             intent(inout)           :: var(:)
!
!        integer(ik)     :: idomain_l, ielement_l, iface
!        type(seed_t)    :: seed
!
!        type(AD_D), allocatable  :: qdiff(:)
!        type(AD_D), allocatable  :: tmp(:)
!        real(rk),   allocatable  :: interpolator(:,:)
!
!        real(rk)        :: xi, eta, zeta
!        integer(ik)     :: nderiv, set_deriv, iterm, nterms_s, ierr, neqns_seed, nterms_s_seed, ipnt, nterms
!        integer(ik)     :: idomain_l_seed, ielement_l_seed, idonor, idomain_l_interp, ielement_l_interp, igq
!        integer(ik)     :: idom_d, ielem_d
!        type(point_t)   :: point_comp, node
!        type(ivector_t) :: donor_domains, donor_elements
!        type(rvector_t) :: donor_xi, donor_eta, donor_zeta
!        logical         :: linearize_me             = .false.
!
!
!
!
!
!        idomain_l  = face%idomain_l
!        ielement_l = face%ielement_l
!        iface      = face%iface
!        seed       = face%seed
!
!
!        !
!        ! Get domain/element index that is being differentiated
!        !
!        idomain_l_seed  = seed%idomain_l
!        ielement_l_seed = seed%ielement_l
!
!
!
!        !
!        ! Get the number of degrees of freedom for the seed element
!        ! and set this as the number of partial derivatives to track
!        !
!        if (ielement_l_seed == 0) then
!            !
!            ! If ielem_seed == 0 then we aren't interested in tracking derivatives
!            !
!            nderiv = 1
!
!        else
!            !
!            ! Get number of equations and terms in solution expansions
!            !
!            neqns_seed    = mesh(idomain_l_seed)%elems(ielement_l_seed)%neqns
!            nterms_s_seed = mesh(idomain_l_seed)%elems(ielement_l_seed)%nterms_s
!
!            !
!            ! Compute number of unknowns in the seed element, which is the number of partial derivatives we are tracking
!            !
!            nderiv = neqns_seed  *  nterms_s_seed
!        end if
!
!
!
!        !
!        ! Allocate the derivative array for each autodiff variable
!        ! MIGHT NOT NEED THIS IF IT GETS AUTOMATICALLY ALLOCATED ON ASSIGNMENT -- TEST
!        !
!        do igq = 1,size(var)
!            allocate(var(igq)%xp_ad_(nderiv))
!        end do
!
!
!
!
!        !
!        ! Find elements associated with boundary nodes
!        !
!        do ipnt = 1,size(points)
!
!
!            call compute_element_donor(mesh, points(ipnt), idom_d, ielem_d, point_comp)
!
!
!            call donor_domains%push_back(  idom_d  )
!            call donor_elements%push_back( ielem_d )
!            call donor_xi%push_back(   point_comp%c1_ )
!            call donor_eta%push_back(  point_comp%c2_ )
!            call donor_zeta%push_back( point_comp%c3_ )
!
!        end do ! ipnt
!
!
!
!
!        !
!        ! For each donor element to the interpolation.
!        !
!        do idonor = 1,donor_elements%size()
!
!
!            !
!            ! Get donor element and coordinate
!            !
!            idom_interp  = donor_domains%at(idonor)
!            ielem_interp = donor_elements%at(idonor)
!
!            xi   = donor_xi%at(  idonor)
!            eta  = donor_eta%at( idonor)
!            zeta = donor_zeta%at(idonor)
!
!            call node%set(xi,eta,zeta)
!
!
!            !
!            ! Get interpolation matrix from quadrature instance
!            !
!            nterms = mesh(idom_interp)%elems(ielem_interp)%nterms_s
!            if (allocated(interpolator)) deallocate(interpolator)
!            allocate(interpolator(1,nterms), stat=ierr)
!            if (ierr /= 0) call AllocationError
!
!            do iterm = 1,nterms
!                interpolator(1,iterm) = polynomial_val(3,nterms,iterm,node)
!            end do ! imode
!
!
!
!
!
!            !
!            ! If the current element is being differentiated (ielem == ielem_seed)
!            ! then copy the solution modes to local AD variable and seed derivatives
!            !
!
!            linearize_me = ( (idom_interp == idom_seed) .and. (ielem_interp == ielem_seed) )
!
!            if ( linearize_me ) then
!
!
!                !
!                ! Allocate AD array to store a copy of the solution which starts the differentiation
!                !
!                nterms_s = mesh(idom_interp)%elems(ielem_interp)%nterms_s
!                if ( allocated(qdiff) ) deallocate(qdiff)
!                allocate(qdiff(nterms_s), stat=ierr)
!                if (ierr /= 0) call AllocationError
!
!
!                !
!                ! Allocate derivative arrays for temporary solution variable
!                !
!                do iterm = 1,nterms_s
!                    allocate(qdiff(iterm)%xp_ad_(nderiv))
!                end do
!
!
!                !
!                ! Copy the solution variables from 'q' to 'qdiff'
!                !
!                qdiff = q%dom(idom_interp)%vecs(ielem_interp)%getvar(ieqn)
!
!
!                !
!                ! Loop through the terms in qdiff and seed derivatives
!                !
!                do iterm = 1,size(qdiff)
!                    ! For the given term, seed its appropriate derivative
!                    set_deriv = (ieqn - 1)*nterms_s + iterm
!                    qdiff(iterm)%xp_ad_(set_deriv) = 1.0_rk
!                end do
!
!
!
!
!                !
!                ! Interpolate solution to GQ nodes via matrix-vector multiplication
!                !
!                !var(idonor) = matmul(interpolator(1,:),  qdiff)
!                tmp = matmul(interpolator,  qdiff)
!                var(idonor) = tmp(1)
!
!            else
!
!                !
!                ! If the solution variable derivatives dont need initialized
!                ! then just use the q(ielem) values and derivatives get
!                ! initialized to zero
!                !
!
!                if ( allocated(tmp) ) deallocate(tmp)
!                allocate(tmp(1))
!                allocate(tmp(1)%xp_ad_(nderiv))
!
!
!                !var(idonor) = matmul(interpolator(1,:),  q%dom(idom_interp)%vecs(ielem_interp)%getvar(ieqn))
!                tmp = matmul(interpolator,  q%dom(idom_interp)%vecs(ielem_interp)%getvar(ieqn))
!                var(idonor) = tmp(1)
!
!            end if
!
!
!
!
!        end do ! idonor
!
!
!    end subroutine interpolate_boundary_autodiff
!    !*************************************************************************************************************

    !>  This routine returns face_info data for the face being interpolated to.
    !!
    !!  Given a face, interpolation type, and interpolation source, this routine returns
    !!      - a face_info_t of the face to be interpolated to
    !!      - an interpolation matrix
    !!      - receiver info to be used if the interpolation is remote
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/16/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    function get_face_interpolation_info(mesh,face_info,interpolation_source,idonor) result(iface_info)
        type(mesh_t),       intent(in)                  :: mesh
        type(face_info_t),  intent(in)                  :: face_info
        integer(ik),        intent(in)                  :: interpolation_source
        integer(ik),        intent(in)                  :: idonor

        type(face_info_t)   :: iface_info
        integer(ik)         :: ChiID
        logical             :: conforming_interpolation, chimera_interpolation, parallel_interpolation

        associate( idom => face_info%idomain_l, ielem => face_info%ielement_l, iface => face_info%iface )

        !
        ! Compute neighbor access indices
        !
        if ( interpolation_source == ME ) then

            ! Interpolate from ME element
            iface_info%idomain_l  = face_info%idomain_l
            iface_info%ielement_l = face_info%ielement_l
            iface_info%idomain_g  = face_info%idomain_g
            iface_info%ielement_g = face_info%ielement_g
            iface_info%iface      = face_info%iface

        elseif ( interpolation_source == NEIGHBOR ) then

            chimera_interpolation    = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA )
            conforming_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )

            ! Interpolate from conforming NEIGHBOR element
            if ( conforming_interpolation ) then
                iface_info%idomain_g  = mesh%domain(idom)%faces(ielem,iface)%ineighbor_domain_g
                iface_info%idomain_l  = mesh%domain(idom)%faces(ielem,iface)%ineighbor_domain_l
                iface_info%ielement_g = mesh%domain(idom)%faces(ielem,iface)%ineighbor_element_g
                iface_info%ielement_l = mesh%domain(idom)%faces(ielem,iface)%ineighbor_element_l
                iface_info%iface      = compute_neighbor_face(mesh,idom,ielem,iface,idonor)         ! THIS PROBABLY NEEDS IMPROVED

            ! Interpolate from CHIMERA donor element
            elseif ( chimera_interpolation ) then
                ChiID = mesh%domain(idom)%faces(ielem,iface)%ChiID
                iface_info%idomain_g  = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%idomain_g
                iface_info%idomain_l  = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%idomain_l
                iface_info%ielement_g = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%ielement_g
                iface_info%ielement_l = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%ielement_l
            else
                call chidg_signal(FATAL,"get_face_interpolation_info: neighbor conforming_interpolation nor chimera_interpolation were detected")
            end if

        else
            call chidg_signal(FATAL,"get_face_interpolation_info: invalid source. ME or NEIGHBOR.")
        end if

        end associate

    end function get_face_interpolation_info
    !*****************************************************************************************

    !>  This returns an interpolation matrix that is used to actually perform the interpolation
    !!  from a modal expansion to a set of quadrature nodes.
    !!
    !!  The interpolation from a modal expansion to a set of quadrature nodes takes the form
    !!  of a matrix-vector multiplication, where the vector entries are modal coefficients of
    !!  the polynomial expansion, and the matrix contains entries that evaluate the modes of
    !!  the polynomial expansion at the quadrature nodes. This routine returns the
    !!  interpolation matrix. Additionally, an interpolation could be computing the actual
    !!  value of the expansion at the nodes('value'), or it could be computing derivatives
    !!  ('grad1', 'grad2', 'grad3'). The interpolation_type specifies what kind of interpolation
    !!  to perform.
    !!
    !!  Given a face, interpolation type, and interpolation source, this routine returns
    !!      - a face_info_t of the face to be interpolated to
    !!      - an interpolation matrix
    !!      - receiver info to be used if the interpolation is remote
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/16/2016
    !!
    !!
    !----------------------------------------------------------------------------------------
    function get_face_interpolation_interpolator(mesh,source_face,interpolation_source,idonor,interpolation_type,donor_face) result(interpolator)
        type(mesh_t),       intent(in)  :: mesh
        type(face_info_t),  intent(in)  :: source_face
        integer(ik),        intent(in)  :: interpolation_source
        integer(ik),        intent(in)  :: idonor
        type(face_info_t),  intent(in)  :: donor_face
        character(*),       intent(in)  :: interpolation_type

        real(rk),       allocatable :: interpolator(:,:)
        integer(ik),    allocatable :: gq_node_indices(:)
        integer(ik)                 :: inode, ChiID, donor_proc
        logical                     :: conforming_interpolation, chimera_interpolation, parallel_interpolation

        associate( idom => source_face%idomain_l, ielem => source_face%ielement_l, iface => source_face%iface )

        !
        ! Compute neighbor access indices
        !
        if ( interpolation_source == ME ) then
            select case(interpolation_type)
                case('value')
                    interpolator = mesh%domain(idom)%faces(ielem,iface)%basis_s%interpolator('Value',iface)
                case('grad1')
                    interpolator = mesh%domain(idom)%faces(ielem,iface)%grad1
                case('grad2')
                    interpolator = mesh%domain(idom)%faces(ielem,iface)%grad2
                case('grad3')
                    interpolator = mesh%domain(idom)%faces(ielem,iface)%grad3
                case default
                    call chidg_signal(FATAL,"get_face_interpolation_interpolator: Invalid interpolation_type. Options are 'value', 'grad1', 'grad2', 'grad3'.")
            end select

        elseif ( interpolation_source == NEIGHBOR ) then

            chimera_interpolation    = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA )
            conforming_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )

            ! Interpolate from conforming NEIGHBOR element
            if ( conforming_interpolation ) then

                parallel_interpolation   = ( IRANK /= mesh%domain(idom)%faces(ielem,iface)%ineighbor_proc )

                ! If parallel use iface_interp to get an interpolation for an opposite face. Assumes same order and eqns.
                if (parallel_interpolation) then
                    select case(interpolation_type)
                        case('value')
                            interpolator = mesh%domain(idom)%faces(ielem,iface)%basis_s%interpolator('Value',donor_face%iface)    ! THIS PROBABLY NEEDS IMPROVED
                        case('grad1')
                            interpolator = mesh%domain(idom)%faces(ielem,iface)%neighbor_grad1
                        case('grad2')
                            interpolator = mesh%domain(idom)%faces(ielem,iface)%neighbor_grad2
                        case('grad3')
                            interpolator = mesh%domain(idom)%faces(ielem,iface)%neighbor_grad3
                        case default
                            call chidg_signal(FATAL,"get_face_interpolation_interpolator: Invalid interpolation_type. Options are 'value', 'grad1', 'grad2', 'grad3'.")
                    end select
                else
                    select case(interpolation_type)
                        case('value')
                            interpolator = mesh%domain(donor_face%idomain_l)%faces(donor_face%ielement_l,donor_face%iface)%basis_s%interpolator('Value',donor_face%iface)
                        case('grad1')
                            interpolator = mesh%domain(donor_face%idomain_l)%faces(donor_face%ielement_l,donor_face%iface)%grad1
                        case('grad2')
                            interpolator = mesh%domain(donor_face%idomain_l)%faces(donor_face%ielement_l,donor_face%iface)%grad2
                        case('grad3')
                            interpolator = mesh%domain(donor_face%idomain_l)%faces(donor_face%ielement_l,donor_face%iface)%grad3
                        case default
                            call chidg_signal(FATAL,"get_face_interpolation_interpolator: Invalid interpolation_type. Options are 'value', 'grad1', 'grad2', 'grad3'.")
                    end select
                end if

            ! Interpolate from CHIMERA donor element
            elseif ( chimera_interpolation ) then
                ChiID = mesh%domain(idom)%faces(ielem,iface)%ChiID
                    select case(interpolation_type)
                        case('value')
                            interpolator = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%value
                        case('grad1')
                            interpolator = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%grad1
                        case('grad2')
                            interpolator = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%grad2
                        case('grad3')
                            interpolator = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%grad3
                        case default
                            call chidg_signal(FATAL,"get_face_interpolation_interpolator: Invalid interpolation_type. Options are 'value', 'grad1', 'grad2', 'grad3'.")
                    end select

            else
                call chidg_signal(FATAL,"get_face_interpolation_interpolator: neighbor conforming_interpolation nor chimera_interpolation were detected")
            end if

        else
            call chidg_signal(FATAL,"get_face_interpolation_interpolator: invalid source. ME or NEIGHBOR.")
        end if

        end associate

    end function get_face_interpolation_interpolator
    !*****************************************************************************************

    !>  This routine returns a logical mask that indicates which nodes in the full
    !!  quadrature node set are to be filled by a chimera donor. This defines the scatter
    !!  pattern for the nodes returned by a chimera donor to the nodes in the full node set.
    !!
    !!  If the interpolation is not a CHIMERA interpolation, then the mask is simply allocated
    !!  to length one and has no meaning.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/16/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    function get_face_interpolation_mask(mesh,face_info,interpolation_source,idonor) result(mask)
        type(mesh_t),       intent(in)                  :: mesh
        type(face_info_t),  intent(in)                  :: face_info
        integer(ik),        intent(in)                  :: interpolation_source
        integer(ik),        intent(in)                  :: idonor

        logical, allocatable :: mask(:) !< This gets returned if CHIMERA interpolation

        integer(ik), allocatable    :: gq_node_indices(:)
        integer(ik)                 :: inode, ChiID, nnodes, ierr
        logical                     :: chimera_interpolation

        associate( idom => face_info%idomain_l, ielem => face_info%ielement_l, iface => face_info%iface )

        if ( interpolation_source == NEIGHBOR ) then

            nnodes = mesh%domain(idom)%faces(ielem,iface)%basis_s%nnodes_if()
            allocate(mask(nnodes), stat=ierr)
            mask = .false.
            if (ierr /= 0) call AllocationError

            chimera_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA )

            !
            ! Interpolate from CHIMERA donor element
            !
            if ( chimera_interpolation ) then
                ChiID = mesh%domain(idom)%faces(ielem,iface)%ChiID
                gq_node_indices = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%node_index(:)

                ! Create mask over full GQ vector of only those nodes that are filled by the current element
                do inode = 1,size(gq_node_indices)
                    mask(gq_node_indices(inode)) = .true.
                end do

            end if

        else
            ! Otherwise, just allocate something so it doens't segfault when it tries to return.
            allocate(mask(1))

        end if

        end associate

    end function get_face_interpolation_mask
    !*****************************************************************************************

    !>  This routine returns a recv_t communication structure indicating whether the
    !!  interpolation donor is LOCAL or REMOTE.
    !!
    !!  The recv_t contains recv_comm, recv_domain, and recv_element components. If these are
    !!  set, then the interpolation is remote and these indices specify where in the 'recv'
    !!  container to find the solution modes for the interpolation. If these indices are not
    !!  set, then the interpolation is LOCAL and the main interpolation routine can use the
    !!  local element indices to locate the solution modes.
    !!
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/16/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    function get_face_interpolation_comm(mesh,face_info,interpolation_source,idonor) result(recv_info)
        type(mesh_t),       intent(in)                  :: mesh
        type(face_info_t),  intent(in)                  :: face_info
        integer(ik),        intent(in)                  :: interpolation_source
        integer(ik),        intent(in)                  :: idonor

        type(recv_t)    :: recv_info                !< This gets set if REMOTE interpolation
        integer(ik)     :: ChiID, donor_proc
        logical         :: conforming_interpolation, chimera_interpolation, parallel_interpolation

        associate( idom => face_info%idomain_l, ielem => face_info%ielement_l, iface => face_info%iface )

        !
        ! Initialize recv_info container to null, indicating LOCAL interpolation.
        ! Always the case, if interpolation_source=ME.
        !
        recv_info = recv_t(0,0,0)

        if ( interpolation_source == NEIGHBOR ) then

            chimera_interpolation    = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA )
            conforming_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )

            ! Interpolate from conforming NEIGHBOR element
            if ( conforming_interpolation ) then

                parallel_interpolation   = ( IRANK /= mesh%domain(idom)%faces(ielem,iface)%ineighbor_proc )

                if (parallel_interpolation) then
                    recv_info%comm    = mesh%domain(idom)%faces(ielem,iface)%recv_comm
                    recv_info%domain  = mesh%domain(idom)%faces(ielem,iface)%recv_domain
                    recv_info%element = mesh%domain(idom)%faces(ielem,iface)%recv_element
                end if

            ! Interpolate from CHIMERA donor element
            elseif ( chimera_interpolation ) then
                ChiID = mesh%domain(idom)%faces(ielem,iface)%ChiID
                donor_proc = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%iproc

                parallel_interpolation = (IRANK /= donor_proc)
                if (parallel_interpolation) then
                     recv_info%comm    = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_comm
                     recv_info%domain  = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_domain
                     recv_info%element = mesh%domain(idom)%chimera%recv(ChiID)%donor(idonor)%recv_element
                end if

            else
                call chidg_signal(FATAL,"get_face_interpolation_comm: neighbor conforming_interpolation nor chimera_interpolation were detected")
            end if

        end if

        end associate

    end function get_face_interpolation_comm
    !*****************************************************************************************

    !>  Return the number of elements contributing to the interpolation.
    !!
    !!  For standard conforming interpolations this will be ==1.
    !!  For Chimera interpolations, this could be >=1.
    !!
    !!  @author Nathan A. Wukie(AFRL)
    !!  @date   8/17/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    function get_face_interpolation_ndonors(mesh,face_info,interpolation_source) result(ndonors)
        type(mesh_t),       intent(in)  :: mesh
        type(face_info_t),  intent(in)  :: face_info
        integer(ik),        intent(in)  :: interpolation_source

        integer(ik) :: ndonors, ChiID
        logical     :: chimera_interpolation, conforming_interpolation

        associate( idom => face_info%idomain_l, ielem => face_info%ielement_l, iface => face_info%iface )

        !
        ! Test if interpolating from local element
        !
        if ( interpolation_source == ME ) then

            ndonors = 1

        !
        ! Test if interpolating from neighbor element(s)
        !
        elseif (interpolation_source == NEIGHBOR ) then

            chimera_interpolation    = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA )
            conforming_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )

            ! Test for standard conforming interpolation from neighbor
            if ( conforming_interpolation ) then
                ndonors = 1

            ! Test for chimera interpolation from neighbor
            elseif ( chimera_interpolation ) then
                ChiID   = mesh%domain(idom)%faces(ielem,iface)%ChiID
                ndonors = mesh%domain(idom)%chimera%recv(ChiID)%ndonors()

            else
                call chidg_signal(FATAL,"get_face_interpolation_ndonors: invalid value for 'face%ftype'")
            end if

        else
            call chidg_signal(FATAL,"get_face_interpolation_ndonors: invalid value for incoming parameter 'source'")
        end if

        end associate

    end function get_face_interpolation_ndonors
    !*****************************************************************************************

    !>  Return the style of the interpolation.
    !!
    !!  Interpolation styles include:
    !!      - 'conforming'    conforming interpolation
    !!      - 'chimera'       chimera interpolation
    !!
    !!  @author Nathan A. Wukie(AFRL)
    !!  @date   8/17/2016
    !!
    !!
    !----------------------------------------------------------------------------------------
    function get_face_interpolation_style(mesh,face_info,interpolation_source) result(interpolation_style)
        type(mesh_t),       intent(in)  :: mesh
        type(face_info_t),  intent(in)  :: face_info
        integer(ik),        intent(in)  :: interpolation_source

        character(len=:),   allocatable :: interpolation_style
        logical                         :: conforming_interpolation, chimera_interpolation

        associate( idom => face_info%idomain_l, ielem => face_info%ielement_l, iface => face_info%iface )

        if ( interpolation_source == ME ) then
            conforming_interpolation = ( (mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR) .or. &
                                         (mesh%domain(idom)%faces(ielem,iface)%ftype == BOUNDARY) .or. &
                                         (mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA) )         ! including chimera here because in the ME case, it doesn't matter
        elseif (interpolation_source == NEIGHBOR ) then
            chimera_interpolation    = ( mesh%domain(idom)%faces(ielem,iface)%ftype == CHIMERA  )
            conforming_interpolation = ( mesh%domain(idom)%faces(ielem,iface)%ftype == INTERIOR )
        else
            call chidg_signal(FATAL,"get_face_interpolation_style: Invalid interpolation_source. ME or NEIGHBOR")
        end if

        if (conforming_interpolation) then
            interpolation_style = "conforming"
        else if (chimera_interpolation) then
            interpolation_style = "chimera"
        else
            call chidg_signal(FATAL,"get_face_interpolation_style: Error in selecting interpolation style")
        end if

        end associate

    end function get_face_interpolation_style
    !******************************************************************************************

    !>  Determine the number of derivatives being computed for the automatic differentiation
    !!  process.
    !!
    !!  This is the number of degrees of freedom in the element being differentiated with
    !!  respect to.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/17/2016
    !!
    !-----------------------------------------------------------------------------------------
    function get_interpolation_nderiv(mesh,function_info) result(nderiv)
        type(mesh_t),           intent(in)  :: mesh
        type(function_info_t),  intent(in)  :: function_info

        integer(ik) :: nderiv, neqns_seed, nterms_s_seed
        logical     :: parallel_seed

        !
        ! If ielem_seed == 0 then we aren't interested in tracking derivatives.
        !   !So set it to lowest number possible while still having something
        !   !allocated in the AD type so the operations are valid.
        !
        !   Actually, allocating with size 0 is okay in fortran. Just need to be
        !   careful not to try and access anything as var%xp_ad_(1)
        !
        if (function_info%seed%ielement_l == 0) then
            !nderiv = 1
            nderiv = 0

        else

            !
            ! Compute number of unknowns in the seed element, which is the number of
            ! partial derivatives we are tracking.
            !
            neqns_seed    = function_info%seed%neqns
            nterms_s_seed = function_info%seed%nterms_s

            nderiv = neqns_seed  *  nterms_s_seed

        end if

    end function get_interpolation_nderiv
    !*****************************************************************************************

end module mod_interpolate

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


module Tables_frac

   ! define table variables for storing melt liquid water content
   ! and melt rate fractions based upon Liang Liao and Bob Meneghinis
   ! melting model calculations

   ! Bill Olson  July, 2013

      integer :: d0_nx, height_ny, mu_nz  ! number of d0, height, and mu values represented in tables

      real :: mu_init, mu_interval ! first mu [] and mu interval [] values
      real :: d0_init, d0_interval ! first d0 [mm] and d0 interval [mm] values
      real :: height_init, height_interval ! first height [m] and height interval [m] values

      integer, allocatable :: max_j(:, :) !level of maximum reflectivity
      real, allocatable :: mlwc_fraction(:, :, :), mrate_fraction(:, :, :) ! mass fraction of melt lwc and rate

end module Tables_frac

      subroutine read_melt_percentages

   ! read tables of melt liquid water content
   ! and melt rate fractions based upon Liang Liao and Bob Meneghinis
   ! melting model calculations

   ! Bill Olson  July, 2013

      use Tables_frac

      implicit none

      logical :: printflag
      integer :: i, j, k

      printflag = .true.

!begin  WSO 12/12/13 updated melting fraction file to Dm indexing
      open(11, file = 'MeltFrac/melt_lwc_rate_percentages_Dm.rho100.dat', &
       form = 'formatted', status='old', action='read')
!end    WSO 12/12/13
      read(11, '(i10, f7.2, f10.4)') mu_nz, mu_init, mu_interval
      read(11, '(i10, f7.4, f12.8)') d0_nx, d0_init, d0_interval
      read(11, '(i10, f7.1, f10.6)') height_ny, height_init, height_interval

      if(printflag) write(6, '("mu_nz:  ", i5, "  mu_init: ", f7.2, "  mu_interval: ", f10.4)'), mu_nz, mu_init, mu_interval
      if(printflag) write(6, '("d0_nx:  ", i5, "  d0_init: ", f7.4, "  d0_interval: ", f12.8)'), d0_nx, d0_init, d0_interval
      if(printflag) write(6, '("height_ny:  ", i5, "  height_init: ", f7.1, "  height_interval: ", f10.6)'), &
       height_ny, height_init, height_interval

      allocate(max_j(d0_nx, mu_nz), mlwc_fraction(d0_nx, height_ny, mu_nz), mrate_fraction(d0_nx, height_ny, mu_nz))

      do k = 1, mu_nz

    ! read lwc melt fractions
        do i = 1, d0_nx
          read(11, '(i5, 2x, 51f5.2)') max_j(i, k), (mlwc_fraction(i, j, k), j = 1, height_ny)
        end do

    ! read rate melt fractions
        do i = 1, d0_nx
          read(11, '(i5, 2x, 51f5.2)') max_j(i, k), (mrate_fraction(i, j, k), j = 1, height_ny)
        end do

      end do

      if(printflag) write(6, '("read the percentages in")')

      return
      end subroutine read_melt_percentages