text:
c ________________________________________________________________
c rutina leelineasii
c datos de la linea

c Se supone que el fichero nomlineas tiene una cabecera con DOS
c lineas de menos de 80 caracteres

c       iln es el numero de la linea dentro del fichero lineas
c       atom es el simbolo del atomo,istage el estado de ionizacion,
c       wlengt la l.d.o. en angstroms
c       zeff es la correccion multiplicativa al damping (gamma6)
c       energy potencial de excitacion del nivel inferior en ev.
c       loggf=log10(peso estadistico (g)* fuerza oscilador (f))
c       mult(i) da el spin total (spin=.5*float(mult(i)-1)),
c       design(i) da el momento angular orbital (ej:S=0,P=1,D=2,F=3...),
c       tam(i) da el momento angular total.
c       code1 representa mediante letras minusculas los momentos angulares
c       orbitales semienteros, asi p=1/2, f=3/2, h=5/2,k=7/2,m=9/2,o=11/2
c                              r=13/2,t=15/2, u=17/2,v=19/2,w=21/2
c
c  Basilio 22 Junio 1994
c  Basilio 23 Junio 1994 modificacion oam semienteros
c ...................................................................

        subroutine leelineasii(filelineas,iln,atom,istage,wlengt,zeff,energy,
     &                       loggf,mult,design,tam,alfa,sigma)
        integer mult(2)
        real tam(2),loggf,wlengt,zeff,alfa,sigma
        character design(2)*1,atom*2,linea*100
        character*1 code1(21)
        character*200 nomlineas
        character*200 filelineas
        character*33 mensaje
c       character*20 cadena
        common/ficlineas/nomlineas
        common/numerito/num
        data code1/'p','1','f','2','h','3','k','4','m','5','o'
     & ,'6','r','7','t','8','u','9','v','0','w'/

        ican=1
        mensaje=' containing the atomic parameters'

        call cabecera(ican,filelineas,mensaje,ifail)
        if(ifail.eq.1)goto 999

c contamos las lineas
        num=0
        nxx=0
        do while(nxx.ne.iln)
           read(ican,'(a)',err=997)linea
           call busco('=',linea,1,igu)
           read (linea(1:igu-1),*,err=998) nxx
           num=num+1
        end do
        close(ican)
        call cabecera(ican,filelineas,mensaje,ifail)

c saltamos las num-1 primeras lineas

        do i=1,num-1
           read(ican,'(a)',err=997)linea
        end do

c leemos la linea como una cadena

        read(ican,'(a)',err=998)linea
c	print*,'leelineasii busco',iln,'encuentro',linea
        close(ican)

c        print*,'leelineasii :',linea

c buscamos el '=' y leemos nxx y atom
        call busco('=',linea,1,igu)
        read (linea(1:igu-1),*,err=998) nxx
        atom=linea(igu+1:igu+3)

c buscamos el proximo blanco
        call busco(' ',linea,igu,ibl1)

        istage=nteroleido(ibl1,ibl2,linea)

c	print*,'istage=',istage
        if(istage.gt.2)then
            print*,' '
            print*,'STOP: No ionization stages higher than 2 can be handled'
            write(*,100) 'The ionization stage',istage,'has been found in the file containing the atomic parameters'
            print*,' '
            print*,'_______________________________________________________________________________'
            stop
        endif
        wlengt=realleido(ibl2,ibl3,linea)
c        print*,'wave=',wlengt,ibl2,ibl3
        zeff=realleido(ibl3,ibl30,linea)
        energy=realleido(ibl30,ibl4,linea)
        loggf=realleido(ibl4,ibl5,linea)

c	print*,'wave=',wlengt
c	print*,'zeff=',zeff
c	print*,'energy=',energy
c	print*,'loggf=',loggf

c buscamos el proximo '-'.
        call busco('-',linea,ibl5,iguion)

c buscamos si existe un corchete antes del guion.El entero esta a la izquierda del corchete
        call busco('[',linea,ibl5,iopen1)
        if(iopen1.lt.iguion)then
           read (linea(iopen1-2:iopen1-1),'(i2)',err=998)mult(1)
           call busco('/',linea,iopen1,ibar1)
           read (linea(iopen1+1:ibar1-1),'(i2)',err=998)ntero1
           design(1)=code1(ntero1)
        else             !no tenemos corchete a la izquierda del -
c El entero empieza 7 characteres antes del guion
           call busco('-',linea,ibl5,iguion)
           iint1=iguion-7
           read (linea(iint1:iint1+1),'(i2)',err=998)mult(1)
           ich1=iint1+2
           design(1)=linea(ich1:ich1)
        end if
        read (linea(ich1+1:ich1+4),'(f4.1)',err=998)tam(1)

c buscamos si existe un corchete despues del guion.
c        print*,'____________________________________'
        call busco('[',linea,iguion,iopen2)
c        print*,'iopen2=',iopen2
c        print*,'____________________________________'

        if(iopen2.lt.100)then
           read (linea(iopen2-2:iopen2-1),'(i2)',err=998)mult(2)
           call busco('/',linea,iopen2,ibar2)
           read (linea(iopen2+1:ibar2-1),'(i2)',err=998)ntero2
           design(2)=code1(ntero2)
        else             !no tenemos corchete a la derecha del -
c El entero empieza 7 characteres despues del guion
           iint2=iguion+1
           read (linea(iint2:iint2+1),'(i2)',err=998)mult(2)
           ich2=iint2+2
           design(2)=linea(ich2:ich2)
        end if
        read (linea(ich2+1:ich2+4),'(f4.1)',err=998)tam(2)

c        print*,'en leelineas2',mult(1),design(1),tam(1)
c        print*,'en leelineas2',mult(2),design(2),tam(2)

        alfa=realleido(ich2+5,ial1,linea)
        sigma=realleido(ial1,ial2,linea)

c        print*,'alfa,sigma=',alfa,sigma
c        print*,'leelineasii num:',num
c        print*,'leelineasii iln:',iln

        return

999     print*,' '
        print*,'STOP: The file containing the atomic parameters does NOT exist.'
        print*,' '
        print*,'_______________________________________________________________________________'
        stop

998     print*,' '
        print*,'STOP: The file containing the atomic parameters is not written correctly'
        print*,'_______________________________________________________________________________'
        stop

997     print*,' '
        write(*,110) 'STOP: The line ',iln,'appearing in the profiles and/or wavelength'
        print*,'grid does NOT exist in the file containing the atomic parameters.'
        print*,' '
        print*,'_______________________________________________________________________________'
        stop

100     format(1x,a20,i2,1x,a59)
110     format(1x,a15,i2,1x,a43)

        end

c ___________________________________________________________________
c realleido      -lee un real si esta entre dos blancos.
c
c entradas:
c     ini        - es la posicion de un blanco anterior al real
c     linea      - es la cadena a leer
c salidas:
c     ifi        - es la posicion del primer blanco posterior al real
c     realleido  - es el real buscado
c
c Basilio -22 de Junio 1994
c ...................................................................

        real*4 function realleido(ini,ifi,linea)
        character linea*100
	character str*13

	str='             '
c buscamos el primer 'no blanco' siguiente a ini
        call nobusco(' ',linea,ini,inbl1)
c buscamos el proximo 'blanco'
        call busco(' ',linea,inbl1,ifi)
c leemos el real contenido en linea(inbl1:ifi-1)
   	str(1:ifi-inbl1)=linea(inbl1:ifi-1)
        read (str,'(f13.5)',err=998) realleido
c        print*,'inibl1=',inbl1,ifi-1
c        read (linea(inbl1:ifi-1),'(f13.5)',err=998) realleido
c	Esto lo tengo que hacer, en lugar de la linea anterior,
c	porque el puto compilador del KIS no le da la gana de
c	cargar linea(inbl1:ifi-1) en realleido.

c	open(1,file='gilipollas')
c	write(1,*) linea(inbl1:ifi-1)
c	close(1)

c	open(1,file='gilipollas')
c	read(1,*) realleido
c	close(1,status='DELETE')

        return

 998    print*,' '
        print*,'STOP: The file containing the atomic parameters is not written correctly.'
        print*,' '
        print*,'_______________________________________________________________________________'
        stop
       end
c ___________________________________________________________________
c nteroleido      -lee un entero si esta entre dos blancos.
c
c entradas:
c     ini        - es la posicion de un blanco anterior al entero
c     linea      - es la cadena a leer
c salidas:
c     ifi        - es la posicion del primer blanco posterior al entero
c     nteroleido  - es el entero buscado
c
c Basilio -22 de Junio 1994
c ...................................................................

        integer*4 function nteroleido(ini,ifi,linea)
        character linea*100

c buscamos el primer 'no blanco' siguiente a ini
        call nobusco(' ',linea,ini,inbl1)
c buscamos el proximo 'blanco'
        call busco(' ',linea,inbl1,ifi)
c leemos el real contenido en linea(inbl1:ifi-1)

        read (linea(inbl1:ifi-1),'(i1)') nteroleido

        return

998     print*,' '
        print*,'STOP: The file containing the atomic parameters is not written correctly.'
        print*,' '
        print*,'_______________________________________________________________________________'
        stop
        end

c _______________________________________________________________
c busco          -da la posicion de un character que se encuentre
c                 a partir de una posicion dada en un string
c entradas:
c     ch         - es el character buscado
c     ini        - es la posicion a partir de la cual comienza la
c                  busqueda
c salidas:
c     ifi        - es la posicion buscada
c
c Basilio -22 de Junio 1994
c ................................................................
        subroutine busco(ch,linea,ini,ifi)
        character ch*1,linea*100

        i=ini
        do while(linea(i:i).ne.ch.and.i.lt.100)
           i=i+1
        end do
        ifi=i
c       print*,ifi
        return
        end
c _______________________________________________________________
c nobusco         -da la primera posicion en que no aparece
c                  un character dado, comenzando la busqueda
c                  a partir duna posicion dada de un string
c entradas:
c     ch         - es el character (NO) buscado
c     ini        - es la posicion a partir de la cual comienza la
c                  busqueda
c salidas:
c     ifi        - es la posicion buscada
c
c Basilio -22 de Junio 1994
c ................................................................

        subroutine nobusco(ch,linea,ini,ifi)
        character ch*1,linea*100

        i=ini
        do while(linea(i:i).eq.ch.and.i.lt.100)
           i=i+1
        end do
        ifi=i
        return
        end


text:
!==============================================================================!
  subroutine Save_Cgns_Cells(grid, sub)
!------------------------------------------------------------------------------!
!   Dummy subroutine                                                           !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
  integer         :: sub
!==============================================================================!

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! 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.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License

!
!////////////////////////////////////////////////////////////////////////
!
!      DictionaryTests.f90
!      Created: January 29, 2013 9:37 AM
!      By: David Kopriva
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE FTDictionaryClassTests
         USE FTDictionaryClass
         USE FTValueClass
         USE FTAssertions
         IMPLICIT NONE

         TYPE(FTDictionary)                               :: dict
         CLASS(FTObject)                        , POINTER :: obj
         TYPE (FTValue)                         , POINTER :: v
         TYPE (FTMutableObjectArray)            , POINTER :: storedObjects
         CHARACTER(LEN=FTDICT_KWD_STRING_LENGTH), POINTER :: storedKeys(:)

         CHARACTER(LEN=FTDICT_KWD_STRING_LENGTH), DIMENSION(4) :: keys   = ["first ","second","third ","fourth"]
         CHARACTER(LEN=FTDICT_KWD_STRING_LENGTH), DIMENSION(4) :: values = ["one  ","two  ","three","four "]
         CHARACTER(LEN=FTDICT_KWD_STRING_LENGTH)               :: s, msg, storedKey, sExpected, sActual
         INTEGER                                               :: i

         CALL dict % initWithSize(64)
!
!        ----------------------
!        Check empty dictionary
!        ----------------------
!
         CALL FTAssertEqual(expectedValue = 0,             &
                            actualValue =  dict % COUNT(), &
                            msg = "Empty dictionary count")
         obj => dict % objectForKey(keys(1))
         CALL FTAssert(test = .NOT.ASSOCIATED(obj),msg = "No object for empty dictionary")

!
!        -----------------------------------------
!        Add the keys and values to the dictionary
!        -----------------------------------------
!
         DO i = 1, 4
            ALLOCATE(v)
            CALL v % initWithValue(values(i))
            obj => v
            CALL dict % addObjectForKey(obj,keys(i))
            CALL releaseFTValue(v)
            CALL FTAssertEqual(1,v % refCount(),"Reference Counting: Addition of object and release")
            CALL FTAssertEqual(i,dict % count(),"Adding to dictionary object count")
         END DO
!
!        ------------------
!        Get them back out
!        ------------------
!
         DO i = 1,4
            obj => dict % objectForKey(keys(i))
            v   => valueFromObject(obj)
            IF ( ASSOCIATED(v) )     THEN
               s = v % stringValue(FTDICT_KWD_STRING_LENGTH)
               CALL FTAssertEqual(values(i),s,"Value for key in dictionary class")
            ELSE
               msg = "Value for key "//TRIM(values(i))// " not of correct type"
               CALL FTAssert(.false.,msg)
            END IF
         END DO
!
!        -------------------------
!        Find the keys and objects
!        -------------------------
!
         storedKeys    => dict % AllKeys()
         storedObjects => dict % AllObjects()

         DO i = 1, 4
            storedKey = storedKeys(i)
            obj       => dict % objectForKey(storedKey)
            v         => valueFromObject(obj)
            sExpected = v % stringValue(FTDICT_KWD_STRING_LENGTH)

            obj      => storedObjects % objectAtIndex(indx = i)
            v        => valueFromObject(obj)
            sActual  =  v % stringValue(FTDICT_KWD_STRING_LENGTH)

            CALL FTAssertEqual(sExpected, sActual,"String for stored key")
         END DO
!
!        ----------------------
!        Check inquiry routines
!        ----------------------
!
         CALL FTAssert( dict % containsKey("first")    ,msg = "dictionary contains key")
         CALL FTAssert( .NOT. dict % containsKey("bob"),msg = "dictionary doesn't contain key")
!
!        ---------------
!        Clean up memory
!        ---------------
!
         DEALLOCATE(storedKeys)
         CALL releaseFTMutableObjectArray(storedObjects)

      END SUBROUTINE FTDictionaryClassTests

max_stars_repo_name: HigherOrderMethods/FTObjectLibrary

max_issues_repo_name: HigherOrderMethods/FTObjectLibrary

max_forks_repo_name: HigherOrderMethods/FTObjectLibrary


text:
Our mission is to promote holistic Medical Services health and healing through our practice and programs. Each of these programs and techniques applied singly or combined give you an uniquely individualized wellness program to support your goals for freedom from pain and restoration of health and vitality for life.
Massage Program  Swedish & Esalen, Cross Fiber
Massage and Biofeedback Somatic Touch© , Aura and Energy Work
Postural and Structural Assessment and Realignment
Bodywork and Energy Work Techniques
OrthoBionomy©
Polarity Therapy
Intuitive Energy Work
Reflexology
Biofeedback Somatic Touch©, Relaxation Breath Techniques and Isolated Active Stretch and Movement Release Techniques, Shock and Trauma Injury Release and Recovery
All Services and techniques are appropriately and deeply effective, yet gentle enough to be able to assist all ages and needs: From the MothertoBe, the New Born, babies, teens, and the Elder Generation
Other: Bach Flower Essences, Affirmations, Gems, Crystals, & Stones, Essential Oils and Nutrition, Life & Health Counseling
Life and Health Consultations: Simple Selfcare Solutions, Harmony and Balance in Life, Dietary and Nutrition, Spiritual Coaching
Seminars, Classes and Tutoring: For Health Maintenance: Self Care Concepts, DeStress, Pain Management and Core Stabilization and Centering Exercises

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
        program wtype11
        call foo(x,i)
        end
        subroutine foo(j,y)
        j=1
        end


! ##################################################################################################################################
! 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 WRITE_PCOMP_EQUIV ( PCOMP_TM, PCOMP_IB, PCOMP_TS )

! Write equiv PSHELL and MAT2's for a PCOMP used, if requested, based on user Bulk Data PARAM PCOMPEQ

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE CONSTANTS_1, ONLY           :  TWELVE
      USE IOUNT1, ONLY                :  ERR, F04, F06, WRT_ERR, WRT_LOG, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, MEMATC, MID1_PCOMP_EQ, MID2_PCOMP_EQ, MID3_PCOMP_EQ,                        &
                                         MID4_PCOMP_EQ, MID1_PCOMP_EQ, MID2_PCOMP_EQ, MID3_PCOMP_EQ, MID4_PCOMP_EQ
      USE PARAMS, ONLY                :  EPSIL, PCOMPEQ, SUPINFO
      USE MODEL_STUF, ONLY            :  INTL_PID, PCOMP, RHO, SHELL_ALP, SHELL_A, SHELL_B, SHELL_D, SHELL_T, SHELL_T_MOD,         &
                                         TREF, ZS

      USE TIMDAT, ONLY                :  TSEC

      USE WRITE_PCOMP_EQUIV_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'WRITE_PCOMP_EQUIV'
      CHARACTER( 8*BYTE)              :: C8FLD_GIJ(3,3,4)   ! Char representation of MAT2 Gij entries
      CHARACTER( 8*BYTE)              :: C8FLD_ALP(6,MEMATC)! Char representation of MAT2 Ai (CTE) entries
      CHARACTER( 8*BYTE)              :: C8FLD_RHO(4)       ! Char representation of MAT2 RHO entry
      CHARACTER( 8*BYTE)              :: C8FLD_TREF(4)      ! Char representation of MAT2 TREF entry
      CHARACTER( 8*BYTE)              :: C8FLD_TM           ! Char representation of MAT2 TM entry
      CHARACTER( 8*BYTE)              :: C8FLD_ZS(2)        ! Char representation of MAT2 ZS entry
      CHARACTER(18*BYTE)              :: NAME1(4)           ! Name of a SHELL matrix (A, D, T, or B)
      CHARACTER( 7*BYTE)              :: NAME2(4)           ! Name of a SHELL matrix (A, D, T, or B)
      CHARACTER( 1*BYTE)              :: FINITE_MAT_PROPS(4)! Indicator of whether any SHELL matrix had zero diag terms

      INTEGER(LONG)                   :: I,J                ! DO loop indices
      INTEGER(LONG)                   :: ICONT              ! Continuation mnemonic for PSHELL equivalent B.D. entry
      INTEGER(LONG)                   :: IERR(4)            ! Error indicator if SHELL_ALP was calculated

      REAL(DOUBLE), INTENT(IN)        :: PCOMP_TM           ! Membrane thickness of PCOMP for equivalent PSHELL
      REAL(DOUBLE), INTENT(IN)        :: PCOMP_IB           ! Bending MOI of PCOMP for equivalent PSHELL
      REAL(DOUBLE), INTENT(IN)        :: PCOMP_TS           ! Transverse shear thickness of PCOMP for equivalent PSHELL
      REAL(DOUBLE)                    :: EPS1               ! Small number
      REAL(DOUBLE)                    :: PCOMP_IBP          ! 12*IB/TM^3
      REAL(DOUBLE)                    :: PCOMP_TSTM         ! TM/TS
      REAL(DOUBLE)                    :: PCOMP_BEN_MAT2(3,3)! MAT2 material matrix for bending for equivalent PSHELL
      REAL(DOUBLE)                    :: PCOMP_MBC_MAT2(3,3)! MAT2 material matrix for mem/ben coupling for equivalent PSHELL
      REAL(DOUBLE)                    :: PCOMP_MEM_MAT2(3,3)! MAT2 material matrix for membrane for equivalent PSHELL
      REAL(DOUBLE)                    :: PCOMP_TSH_MAT2(2,2)! MAT2 material matrix for transverse shear for equivalent PSHELL

! **********************************************************************************************************************************
      EPS1 = EPSIL(1)

! Determine if any SHELL_A, D, T, B  matrices are null

      FINITE_MAT_PROPS(1) = 'Y'
      IF ((DABS(SHELL_A(1,1)) < EPS1) .OR. (DABS(SHELL_A(1,1)) < EPS1) .OR. (DABS(SHELL_A(1,1)) < EPS1) .OR.                       &
          (DABS(SHELL_A(1,1)) < EPS1)) THEN
         FINITE_MAT_PROPS(1) = 'N'
      ENDIF

      FINITE_MAT_PROPS(2) = 'Y'
      IF ((DABS(SHELL_D(1,1)) < EPS1) .OR. (DABS(SHELL_D(1,1)) < EPS1) .OR. (DABS(SHELL_D(1,1)) < EPS1) .OR.                       &
          (DABS(SHELL_D(1,1)) < EPS1)) THEN
         FINITE_MAT_PROPS(2) = 'N'
      ENDIF

      FINITE_MAT_PROPS(3) = 'Y'
      IF ((DABS(SHELL_T(1,1)) < EPS1) .OR. (DABS(SHELL_T(1,1)) < EPS1) .OR. (DABS(SHELL_T(1,1)) < EPS1)) THEN
         FINITE_MAT_PROPS(3) = 'N'
      ENDIF

      FINITE_MAT_PROPS(4) = 'Y'
      IF ((DABS(SHELL_B(1,1)) < EPS1) .OR. (DABS(SHELL_B(1,1)) < EPS1) .OR. (DABS(SHELL_B(1,1)) < EPS1) .OR.                       &
          (DABS(SHELL_B(1,1)) < EPS1)) THEN
         FINITE_MAT_PROPS(4) = 'N'
      ENDIF

! Write message if any CTE's were not able to be calculated, but only if that SHELL matrix had props to be output

      NAME1(1) = 'MEMBRANE'  ;   NAME1(2) = 'BENDING'   ;   NAME1(3) = 'TRANSVERSE SHEAR'   ;   NAME1(4) = 'MEMB/BEND COUPLING'
      NAME2(1) = 'SHELL_A'   ;   NAME2(2) = 'SHELL_D'   ;   NAME2(3) = 'SHELL_T'            ;   NAME2(4) = 'SHELL_B'

      CALL SOLVE_SHELL_ALP ( IERR )                         ! Solve for equiv CTE's for this PCOMP

      WRITE(F06,9900)
      DO I=1,4
         IF      (IERR(I) == 1) THEN
            IF (FINITE_MAT_PROPS(I) == 'Y') THEN
               WRITE(ERR,9801) NAME1(I), PCOMP(INTL_PID,1), NAME2(I)
               IF (SUPINFO == 'N') THEN
                  WRITE(F06,9801) NAME1(I), PCOMP(INTL_PID,1), NAME2(I)
               ENDIF
            ENDIF
         ELSE IF (IERR(I) == 2) THEN
            IF (FINITE_MAT_PROPS(I) == 'Y') THEN
               WRITE(ERR,9802) NAME1(I), PCOMP(INTL_PID,1), NAME2(I)
               IF (SUPINFO == 'N') THEN
                  WRITE(F06,9802) NAME1(I), PCOMP(INTL_PID,1), NAME2(I)
               ENDIF
            ENDIF
         ENDIF
      ENDDO

      PCOMP_IBP  = TWELVE*PCOMP_IB/(PCOMP_TM*PCOMP_TM*PCOMP_TM)
      PCOMP_TSTM = PCOMP_TS/PCOMP_TM
      WRITE(F06,9901) PCOMP(INTL_PID,1), PCOMP_TM, PCOMP_IBP, PCOMP_TSTM
      IF (SUPINFO == 'N') THEN
         WRITE(F06,9901) PCOMP(INTL_PID,1), PCOMP_TM, PCOMP_IBP, PCOMP_TSTM
      ENDIF

      DO I=1,3
         DO J =1,3
            PCOMP_MEM_MAT2(I,J) = SHELL_A(I,J)/PCOMP_TM
            PCOMP_MBC_MAT2(I,J) = SHELL_B(I,J)
            PCOMP_BEN_MAT2(I,J) = SHELL_D(I,J)/PCOMP_IB
         ENDDO
      ENDDO

      DO I=1,2
         DO J =1,2
            PCOMP_TSH_MAT2(I,J) = SHELL_T(I,J)/PCOMP_TS
         ENDDO
      ENDDO

      MID1_PCOMP_EQ = MID1_PCOMP_EQ + PCOMP(INTL_PID,1)
      MID2_PCOMP_EQ = MID2_PCOMP_EQ + PCOMP(INTL_PID,1)
      MID3_PCOMP_EQ = MID3_PCOMP_EQ + PCOMP(INTL_PID,1)
      MID4_PCOMP_EQ = MID4_PCOMP_EQ + PCOMP(INTL_PID,1)

      ICONT = MID1_PCOMP_EQ

      IF (PCOMPEQ > 1) THEN

         IF ((IERR(1) == 0) .OR. (IERR(2) == 0) .OR. (IERR(3) == 0) .OR. (IERR(4) == 0)) THEN
            WRITE(F06,9902)
         ENDIF

         IF (FINITE_MAT_PROPS(1) == 'Y') THEN
            IF (IERR(1) == 0) THEN
               WRITE(F06,9903) MID1_PCOMP_EQ, PCOMP_MEM_MAT2(1,1), PCOMP_MEM_MAT2(1,2), PCOMP_MEM_MAT2(1,3),                       &
                                              PCOMP_MEM_MAT2(2,2), PCOMP_MEM_MAT2(2,3), PCOMP_MEM_MAT2(3,3),                       &
                                              SHELL_ALP(1,1), SHELL_ALP(2,1), SHELL_ALP(3,1)
            ELSE
               WRITE(F06,9903) MID1_PCOMP_EQ, PCOMP_MEM_MAT2(1,1), PCOMP_MEM_MAT2(1,2), PCOMP_MEM_MAT2(1,3),                       &
                                              PCOMP_MEM_MAT2(2,2), PCOMP_MEM_MAT2(2,3), PCOMP_MEM_MAT2(3,3)
            ENDIF
         ENDIF

         IF (FINITE_MAT_PROPS(2) == 'Y') THEN
            IF (IERR(2) == 0) THEN
               WRITE(F06,9904) MID2_PCOMP_EQ, PCOMP_BEN_MAT2(1,1), PCOMP_BEN_MAT2(1,2), PCOMP_BEN_MAT2(1,3),                       &
                                              PCOMP_BEN_MAT2(2,2), PCOMP_BEN_MAT2(2,3), PCOMP_BEN_MAT2(3,3),                       &
                                              SHELL_ALP(1,2), SHELL_ALP(2,2), SHELL_ALP(3,2)
            ELSE
               WRITE(F06,9904) MID2_PCOMP_EQ, PCOMP_BEN_MAT2(1,1), PCOMP_BEN_MAT2(1,2), PCOMP_BEN_MAT2(1,3),                       &
                                              PCOMP_BEN_MAT2(2,2), PCOMP_BEN_MAT2(2,3), PCOMP_BEN_MAT2(3,3)
            ENDIF
         ENDIF

         IF (FINITE_MAT_PROPS(3) == 'Y') THEN
            IF (IERR(3) == 0) THEN
               WRITE(F06,9905) MID3_PCOMP_EQ, PCOMP_TSH_MAT2(1,1), PCOMP_TSH_MAT2(1,2), PCOMP_TSH_MAT2(2,2),                       &
                                              SHELL_ALP(5,3), SHELL_ALP(6,3)
            ELSE
               WRITE(F06,9905) MID3_PCOMP_EQ, PCOMP_TSH_MAT2(1,1), PCOMP_TSH_MAT2(1,2), PCOMP_TSH_MAT2(2,2)
            ENDIF
         ENDIF

         IF (FINITE_MAT_PROPS(4) == 'Y') THEN
            IF (IERR(4) == 0) THEN
               WRITE(F06,9906) MID4_PCOMP_EQ, PCOMP_MBC_MAT2(1,1), PCOMP_MBC_MAT2(1,2), PCOMP_MBC_MAT2(1,3),                       &
                                              PCOMP_MBC_MAT2(2,2), PCOMP_MBC_MAT2(2,3), PCOMP_MBC_MAT2(3,3),                       &
                                              SHELL_ALP(1,4), SHELL_ALP(2,4), SHELL_ALP(3,4)
            ELSE
               WRITE(F06,9906) MID4_PCOMP_EQ, PCOMP_MBC_MAT2(1,1), PCOMP_MBC_MAT2(1,2), PCOMP_MBC_MAT2(1,3),                       &
                                              PCOMP_MBC_MAT2(2,2), PCOMP_MBC_MAT2(2,3), PCOMP_MBC_MAT2(3,3)
            ENDIF
         ENDIF

         WRITE(F06,*)

      ENDIF

      CALL GET_CHAR8_OUTPUTS ( C8FLD_GIJ, C8FLD_ALP, C8FLD_RHO, C8FLD_TREF, C8FLD_TM, C8FLD_ZS )

! Write PSHELL

      IF (IERR(4) == 0) THEN
         WRITE(F06,9912) PCOMP(INTL_PID,1), MID1_PCOMP_EQ, C8FLD_TM, MID2_PCOMP_EQ, MID3_PCOMP_EQ,ICONT,                           &
                         ICONT, C8FLD_ZS(2), C8FLD_ZS(2), MID4_PCOMP_EQ

      ELSE
         WRITE(F06,9913) PCOMP(INTL_PID,1), MID1_PCOMP_EQ, C8FLD_TM, MID2_PCOMP_EQ, MID3_PCOMP_EQ,ICONT,                           &
                            ICONT, C8FLD_ZS(1), C8FLD_ZS(2)
      ENDIF

! Write MAT2 for membrane

      IF (FINITE_MAT_PROPS(1) == 'Y') THEN
         IF (IERR(1) == 0) THEN

            WRITE(F06,9922) MID1_PCOMP_EQ, C8FLD_GIJ(1,1,1), C8FLD_GIJ(1,2,1), C8FLD_GIJ(1,3,1),                                   &
                                           C8FLD_GIJ(2,2,1), C8FLD_GIJ(2,3,1), C8FLD_GIJ(3,3,1), C8FLD_RHO(1),ICONT+1,             &
                                           ICONT+1, C8FLD_ALP(1,1), C8FLD_ALP(2,1), C8FLD_ALP(3,1), C8FLD_TREF(1)
         ELSE
            WRITE(F06,9923) MID1_PCOMP_EQ, C8FLD_GIJ(1,1,1), C8FLD_GIJ(1,2,1), C8FLD_GIJ(1,3,1),                                   &
                                           C8FLD_GIJ(2,2,1), C8FLD_GIJ(2,3,1), C8FLD_GIJ(3,3,1), C8FLD_RHO(1),ICONT+1,             &
                                           ICONT+1, C8FLD_TREF(1)
         ENDIF
      ENDIF

! Write MAT2 for bending

      IF (FINITE_MAT_PROPS(2) == 'Y') THEN
         IF (IERR(2) == 0) THEN
            WRITE(F06,9922) MID2_PCOMP_EQ, C8FLD_GIJ(1,1,2), C8FLD_GIJ(1,2,2), C8FLD_GIJ(1,3,2),                                   &
                                           C8FLD_GIJ(2,2,2), C8FLD_GIJ(2,3,2), C8FLD_GIJ(3,3,2), C8FLD_RHO(2),ICONT+2,             &
                                           ICONT+2, C8FLD_ALP(1,2), C8FLD_ALP(2,2), C8FLD_ALP(3,2), C8FLD_TREF(1)
         ELSE

            WRITE(F06,9923) MID2_PCOMP_EQ, C8FLD_GIJ(1,1,2), C8FLD_GIJ(1,2,2), C8FLD_GIJ(1,3,2),                                   &
                                           C8FLD_GIJ(2,2,2), C8FLD_GIJ(2,3,2), C8FLD_GIJ(3,3,2), C8FLD_RHO(2),ICONT+2,             &
                                           ICONT+2, C8FLD_TREF(1)
         ENDIF
      ENDIF

! Write MAT2 for transverse shear

      IF (FINITE_MAT_PROPS(3) == 'Y') THEN
         IF (IERR(3) == 0) THEN
            WRITE(F06,9932) MID3_PCOMP_EQ, C8FLD_GIJ(1,1,3), C8FLD_GIJ(1,2,3), '0.0000+0'      ,                                   &
                                           C8FLD_GIJ(2,2,3), '0.0000+0'      , '0.0000+0', C8FLD_RHO(3), ICONT+3,                  &
                                           ICONT+3, C8FLD_ALP(5,3), C8FLD_ALP(6,3), C8FLD_TREF(1)
         ELSE

            WRITE(F06,9933) MID3_PCOMP_EQ, C8FLD_GIJ(1,1,3), C8FLD_GIJ(1,2,3), '0.0000+0'      ,                                   &
                                           C8FLD_GIJ(2,2,3), '0.0000+0'      , '0.0000+0', C8FLD_RHO(3), ICONT+3,                  &
                                           ICONT+3, C8FLD_TREF(1)
         ENDIF
      ENDIF

! Write MAT2 for membrane/bending coupling

      IF (FINITE_MAT_PROPS(4) == 'Y') THEN
         IF (IERR(4) == 0) THEN
            WRITE(F06,9922) MID4_PCOMP_EQ, C8FLD_GIJ(1,1,4), C8FLD_GIJ(1,2,4), C8FLD_GIJ(1,3,4),                                   &
                                           C8FLD_GIJ(2,2,4), C8FLD_GIJ(2,3,4), C8FLD_GIJ(3,3,4), C8FLD_RHO(4),ICONT+4,             &
                                           ICONT+4, C8FLD_ALP(1,4), C8FLD_ALP(2,4), C8FLD_ALP(3,4), C8FLD_TREF(1)
         ELSE

            WRITE(F06,9923) MID4_PCOMP_EQ, C8FLD_GIJ(1,1,4), C8FLD_GIJ(1,2,4), C8FLD_GIJ(1,3,4),                                   &
                                           C8FLD_GIJ(2,2,4), C8FLD_GIJ(2,3,4), C8FLD_GIJ(3,3,4), C8FLD_RHO(4),ICONT+4,             &
                                           ICONT+4, C8FLD_TREF(1)
         ENDIF
      ENDIF

! Write message if we have changed SHELL_T to reflect approx zero transverse shear flex

      IF (SHELL_T_MOD == 'Y') THEN
         WRITE(ERR,9995) MID3_PCOMP_EQ
         IF (SUPINFO == 'N') THEN
            WRITE(F06,9995) MID3_PCOMP_EQ
         ENDIF
      ENDIF

      WRITE(F06,9900)

      PCOMP(INTL_PID,6) = 1                                 ! Lets future calls to this subr know that PSHELL, MAT2 were written

! **********************************************************************************************************************************
 9801 FORMAT(' *INFORMATION: Cannot calculate equiv CTE''s for ',A,' for PCOMP ',I8,' since the det of matrix ',A,' is zero',/)

 9802 FORMAT(' *INFORMATION: Cannot calculate equiv CTE''s for ',A,' for PCOMP ',I8,' since matrix ',A,' cannot be inverted',/)

 9900 FORMAT('++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++', &
             '+++++++++++++++++++++++++++++++++++++++++++++++++')

 9901 FORMAT(' *INFORMATION: Equivalent PSHELL amd MAT2 entries for PCOMP ',I8,' with:'                                            &
                    ,/,14X,' membrane thickness TM = ',1ES13.6,',   12*IB/(TM^3) = ',1ES13.6,',   TS/TM = ',1ES13.6,':',/)

 9902 FORMAT('        Type         MATL ID       G11            G12            G13            G22            G23            G33',  &
             '             A1             A2             A3',/                                                                     &
             ' ------------------ --------  -------------  -------------  -------------  -------------  -------------',            &
             '  -------------  -------------  -------------  -------------')

 9903 FORMAT('      membrane      ',I8,9(1ES15.6))

 9904 FORMAT('      bending       ',I8,9(1ES15.6))

 9905 FORMAT('  transverse shear  ',I8,2(1ES15.6),15X,1ES15.6,30X,2(1ES15.6))

 9906 FORMAT('  mem/bend coupling ',I8,9(1ES15.6))

 9912 FORMAT('PSHELL   ,',I9,',',I9,',',A9,',',I9,',      1.0,',I9,',','      1.0,         , +',I7,/,' +',I7,2(',',A9),I9,/)

 9913 FORMAT('PSHELL   ,',I9,',',I9,',',A9,',',I9,',      1.0,',I9,',','      1.0,         , +',I7,/,' +',I7,2(',',A9),/)

 9922 FORMAT('MAT2     ,',I9,7(',',A9),', +',I7,/,' +',I7,3(',',A9),',',A9,/)

 9923 FORMAT('MAT2     ,',I9,7(',',A9),', +',I7,/,' +',I7,3(',',8X),',',A9,/)

 9932 FORMAT('MAT2     ,',I9,7(',',A9),', +',I7,/,' +',I7,2(',',A9),',',9X,',',A9,/)

 9933 FORMAT('MAT2     ,',I9,7(',',A9),', +',I7,/,' +',I7,2(',',8X),',',9X,',',A9,/)

 9995 FORMAT(' *INFORMATION: The transverse shear modulii on the above MAT2 ',I8,' entry are the MYSTRAN calculated values to',    &
                           ' replace the zero input values by the user that',/,14X,' were meant to simulate zero transverse shear',&
                           ' flexibility')

! ##################################################################################################################################

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE GET_CHAR8_OUTPUTS ( C8FLD_GIJ, C8FLD_ALP, C8FLD_RHO, C8FLD_TREF, C8FLD_TM, C8FLD_ZS )

      IMPLICIT NONE

      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_GIJ(3,3,4)    ! Char representation of MAT2 Gij entries
      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_ALP(6,MEMATC) ! Char representation of MAT2 CTE entries
      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_RHO(4)        ! Char representation of MAT2 RHO entry
      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_TREF(4)       ! Char representation of MAT2 TREF entry
      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_TM            ! Char representation of MAT2 TM entry
      CHARACTER(8*BYTE), INTENT(OUT)  :: C8FLD_ZS(2)         ! Char representation of MAT2 ZS entry

      INTEGER(LONG)                   :: II,JJ,KK            ! DO loop indices

! **********************************************************************************************************************************
      DO II=1,3
         DO JJ=1,3
            DO KK=1,4
               C8FLD_GIJ(II,JJ,KK)(1:8) = ' '
            ENDDO
         ENDDO
      ENDDO

      DO II=1,6
         DO JJ=1,MEMATC
            C8FLD_ALP(II,JJ)(1:8) = ' '
         ENDDO
      ENDDO

! Put GIJ values in character format

      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(1,1) ,C8FLD_GIJ(1,1,1) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(1,2) ,C8FLD_GIJ(1,2,1) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(1,3) ,C8FLD_GIJ(1,3,1) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(2,2) ,C8FLD_GIJ(2,2,1) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(2,3) ,C8FLD_GIJ(2,3,1) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MEM_MAT2(3,3) ,C8FLD_GIJ(3,3,1) )

      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(1,1) ,C8FLD_GIJ(1,1,2) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(1,2) ,C8FLD_GIJ(1,2,2) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(1,3) ,C8FLD_GIJ(1,3,2) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(2,2) ,C8FLD_GIJ(2,2,2) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(2,3) ,C8FLD_GIJ(2,3,2) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_BEN_MAT2(3,3) ,C8FLD_GIJ(3,3,2) )

      CALL REAL_DATA_TO_C8FLD ( PCOMP_TSH_MAT2(1,1) ,C8FLD_GIJ(1,1,3) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_TSH_MAT2(1,2) ,C8FLD_GIJ(1,2,3) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_TSH_MAT2(2,2) ,C8FLD_GIJ(2,2,3) )

      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(1,1) ,C8FLD_GIJ(1,1,4) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(1,2) ,C8FLD_GIJ(1,2,4) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(1,3) ,C8FLD_GIJ(1,3,4) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(2,2) ,C8FLD_GIJ(2,2,4) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(2,3) ,C8FLD_GIJ(2,3,4) )
      CALL REAL_DATA_TO_C8FLD ( PCOMP_MBC_MAT2(3,3) ,C8FLD_GIJ(3,3,4) )

! Put CTE's in character format

      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(1,1), C8FLD_ALP(1,1) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(2,1), C8FLD_ALP(2,1) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(3,1), C8FLD_ALP(3,1) )

      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(1,2), C8FLD_ALP(1,2) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(2,2), C8FLD_ALP(2,2) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(3,2), C8FLD_ALP(3,2) )

      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(5,3), C8FLD_ALP(5,3) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(6,3), C8FLD_ALP(6,3) )

      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(1,4), C8FLD_ALP(1,4) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(2,4), C8FLD_ALP(2,4) )
      CALL REAL_DATA_TO_C8FLD ( SHELL_ALP(3,4), C8FLD_ALP(3,4) )

! Put RHO's in character format

      CALL REAL_DATA_TO_C8FLD ( RHO(1), C8FLD_RHO(1) )
      CALL REAL_DATA_TO_C8FLD ( RHO(2), C8FLD_RHO(2) )
      CALL REAL_DATA_TO_C8FLD ( RHO(3), C8FLD_RHO(3) )
      CALL REAL_DATA_TO_C8FLD ( RHO(4), C8FLD_RHO(4) )

! Put TREF's in character format

      CALL REAL_DATA_TO_C8FLD ( TREF(1), C8FLD_TREF(1) )
      CALL REAL_DATA_TO_C8FLD ( TREF(2), C8FLD_TREF(2) )
      CALL REAL_DATA_TO_C8FLD ( TREF(3), C8FLD_TREF(3) )
      CALL REAL_DATA_TO_C8FLD ( TREF(4), C8FLD_TREF(4) )

! Put TM in character format

      CALL REAL_DATA_TO_C8FLD ( PCOMP_TM, C8FLD_TM )

! Put ZS's in character format

      CALL REAL_DATA_TO_C8FLD ( ZS(1), C8FLD_ZS(1) )
      CALL REAL_DATA_TO_C8FLD ( ZS(2), C8FLD_ZS(2) )

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

      END SUBROUTINE GET_CHAR8_OUTPUTS

      END SUBROUTINE WRITE_PCOMP_EQUIV


text:
      subroutine gensl1
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  This subroutine sets up for and calls the necessary routines to
CD1  solve the heat and mass transfer equations with full
CD1  jacobian (unsymmetric, 2n by 2n).
CD1
C***********************************************************************
CD2
CD2  REVISION HISTORY
CD2
CD2 $Log:   /pvcs.config/fehm90/src/gensl1.f_a  $
CD2
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:06   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:06:34   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:09:28   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:24:10   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:02:34   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:41:46 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.16   Fri Apr 26 15:20:50 1996   gaz
CD2 coding added for mdnodes
CD2
CD2    Rev 1.15   Wed Jan 17 08:57:30 1996   hend
CD2 Removed nonconstant parameter statement for IBM
CD2
CD2    Rev 1.14   Tue Jan 16 14:31:08 1996   hend
CD2 Added capability for 5,6, and n degrees of freedom
CD2
CD2    Rev 1.13   Fri Jan 12 17:48:52 1996   llt
CD2 changed mmgetblk agruments
CD2
CD2    Rev 1.12   Wed Jan 10 12:28:56 1996   hend
CD2 Corrected ECD Number
CD2
CD2    Rev 1.11   Wed Jan 10 12:01:46 1996   hend
CD2 Added Prolog
CD2
CD2    Rev 1.10   08/02/95 17:30:40   llt
CD2 moved where sto5 was getting released
CD2
CD2    Rev 1.9   06/01/95 16:48:06   gaz
CD2 minor corrections f0=-1
CD2
CD2    Rev 1.8   04/25/95 09:06:34   llt
CD2 retrieved lost log history information
CD2
CD2    Rev 1.7   03/23/95 18:50:46   gaz
CD2 No change.
CD2
CD2    Rev 1.6   03/23/95 18:43:56   gaz
CD2 gaz changed bp to abs(bp in one if block
CD2
CD2    Rev 1.5   03/23/95 18:26:36   gaz
CD2 gaz changes to call to rdof_new and tolerance options
CD2
CD2    Rev 1.4   03/10/95 10:45:48   llt
CD2 removed unneeded pointers - gaz
CD2
CD2    Rev 1.3   05/11/94 16:14:50   llt
CD2 bug fixes - gaz
CD2
CD2    Rev 1.2   03/23/94 14:43:12   robinson
CD2 Additional cleanup of memory management
CD2
CD2    Rev 1.1   03/18/94 15:59:32   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:24:22   pvcs
CD2 original version in process of being certified
CD2
C**********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.3.2 Heat- and mass-transfer equations
CD3  2.5.2 Solve nonlinear equation set at each time step
CD3
C**********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C***********************************************************************

      use davidi
      use comflow
      use comhi
      use comgi
      use comfi
      use comei
      use comdi
      use comci
      use combi
      use comdti
      use comai
      use comriv
      use comcouple
      implicit none

      integer ndex(2)
      integer i, id, neqp1, nsizea, nsizea1
      integer iparchek,kb,jj,i1,i2,i3,i4,j
      real a11
      real*8 dumn(100),term
      real*8, allocatable :: sto5(:,:)
      real*8, allocatable :: dum(:)
      real*8  facr, fdum2, tollr, tolls, tmch_new, tmch_old, bp_max
      save tmch_old
      parameter(iparchek =0)
      neqp1=neq+1
c test of coefficients
       if(iparchek.eq.1) then
       write(ierr,*) 'fe coef. neighbors'
       write(ierr,*)
       do i = 1,neq
        write(ierr,*)' vol node ', i,' = ', sx1(i)
        i1 = nelm(i)+1
        i2 = nelmdg(i)-1
        do j = i1,i2
         kb = nelm(j)
         i3 = nelmdg(kb)+1
         i4 = nelm(kb+1)
         do jj = i3,i4
         if(nelm(jj).eq.i) then
          iw = istrw(jj-neqp1)
          if(iw.gt.0) then
           a11 = sx(iw,isox)+sx(iw,isoy)+sx(iw,isoz)
          else
           a11 = 0.0
          endif
          write(ierr,355)i,kb,iw,a11
          go to 899
         endif
        enddo
899     continue
       enddo
        i1 = nelmdg(i)+1
        i2 = nelm(i+1)
        do j = i1,i2
         kb = nelm(j)
         iw = istrw(j-neqp1)
          if(iw.gt.0) then
           a11 = sx(iw,isox)+sx(iw,isoy)+sx(iw,isoz)
          else
           a11 = 0.0
          endif
         write(ierr,355)i,kb,iw,a11
        enddo
       enddo
355   format(i6,1x,i6,1x,i6,1x,f12.3)
       stop
       endif
c     zero out arrays
c
c  check for possible stopping on variable changes
c
      if(iad.ge.1) then
       if(nr_stop.eq.2) then
        fdum=-1.0
        go to 999
       endif
      endif
c
      do i=1,neq
         bp(i+nrhs(1))=0.0
         bp(i+nrhs(2))=0.0
      enddo
      nsizea1=nelm(neqp1)-neqp1
      nsizea=4*nsizea1
      do i=1,nsizea
         a(i)=0.0
      enddo
      fdum2=0.
      do id=1,neq
        if(iriver.eq.2.and.id.gt.neq_primary) then
         call geneq1_well(id)
        else if(ianpe.ne.0) then
         call geneq1_ani(id)
        else
         call geneq1(id)
        endif
      enddo
      do id=1,neq
         if(ps(id).le.0.0) then
            a(nelmdg(id)-neqp1+nmat(1))=sx1(id)
            bp(id+nrhs(1))=0.0
         endif
      enddo

c     call md_modes to complete equations for multiply defined nodes
c don't need to call geneqmdnode (through md_nodes) gaz
c     if(imdnode.ne.0) call md_nodes(3,0,0)

      call dual(10)
c
c possible correction for rate-limited gdpm model
c
      call gdpm_corr(1)
c
c s kelkar for debugging 3/20/12
c      call extract_adiag
c      call extract_nodal_equations(2627)
c
      if(islord.ne.0) call switch(nmat,nmatb,islord,2,1,nsizea1)
      if(islord.ne.0) call switchb(bp,nrhs,nrhsb,islord,2,1,neq)
c
c normalize equations
c
      call normal_dof(neq,a,bp,nelm,nmat,nrhs,nelmdg
     &     ,ndex,2,dumn(1),dumn(37),dumn(73),0,fdum2)
      if(ndex(1).lt.0) then
         write(ierr,*)
     & '>>> cannot normalize, stopping, i=',abs(ndex(1)),'<<<'
         if(iout.ne.0) then
          write(iout,*)
     & '>>> cannot normalize, stopping, i=',abs(ndex(1)),'<<<'
         endif
         if(iptty.ne.0) then
          write(iptty,*)
     & '>>> cannot normalize, stopping, i=',abs(ndex(1)),'<<<'
         endif
         stop
      endif

      fdum=sqrt(fdum2)
      if(fdum.eq.0.0) go to 999
      if(iad.eq.0) then
         f0=max(fdum*epe,tmch)
      endif
      if(fdum1.lt.0.0.and.iad.ne.0) then
	  bp_max = 0.0
         do i=1,neq
            bp_max= max(abs(bp(i+nrhs(1))),abs(bp(i+nrhs(2))),bp_max)
         enddo
	 if(bp_max.lt.tmch.and.nr_stop.eq.0) then
         fdum=-1.0
         go to 999
       else
         f0=-1.0
	   if(iad.eq.1) tmch_old = bp_max
	   tmch_new = min(bp_max,tmch_old)
	   tmch_old = tmch_new
	 endif
      endif
      if(f0.gt.0.0d00) then
         if(fdum.le.f0.and.iad.ne.0) goto 999
         facr=1.0d00
         tolls=max(facr*f0,min(g1*fdum,g2*fdum2))
         tollr=tolls*g3
         if(tmch*g3.gt.tollr) tollr=tmch*g3
      else if(f0.le.0.0d00) then
         tollr=tmch_new*g3
      endif

c     time=second(ghv)
c
c      times=second(0.0)
      iter=maxsolve
      call storage_derivatives(1,1)

      if(irdof.eq.0) then
c
c       full solution
c

         allocate(dum(neq*8))
         if(gdpm_flag.eq.0) then
            if (igauss .gt. 1) then
               call solve_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,
     &              north,tollr,irb,iirb,npvt,gmres,dum,piv,
     &              h,c,ss,g,y,iter,iback,2,iptty,maxor,accm)
            else
               call solve_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,
     &              north,tollr,irb,iirb,npvt,gmres,dum,piv,
     &              h,c,ss,g,y,iter,iback,2,iptty,maxor,accm)
            end if
         else
            if (igauss .gt. 1) then
               call solve_dual(neq_primary,neq,a,b,bp,nmat,nb,nrhs,
     &              nelm,nelm_primary,nop,north,tollr,irb,iirb,
     &              npvt,gmres,dum,piv,h,c,ss,g,y,iter,iback,2,
     &              iptty,maxor,igdpm,maxgdpmlayers,ngdpm_layers,
     &              nelmdg,accm,mdof_sol)
            else
               call solve_dual(neq_primary,neq,a,b,bp,nmat,nmat,nrhs,
     &              nelm,nelm_primary,nelm_primary,north,tollr,irb,iirb,
     &              npvt,gmres,dum,piv,h,c,ss,g,y,iter,iback,2,iptty,
     &              maxor,igdpm,maxgdpmlayers,ngdpm_layers,nelmdg,
     &              accm,mdof_sol)
            end if
         endif
	   deallocate(dum)
         itert=itert+iter
         itotals=itotals+iter
         minkt=minkt+mink
      else

         allocate(sto5(neq,2))
         allocate(dum(neq*8))
         if (igauss .gt. 1) then
            call rdof_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         else
            call rdof_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,2,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         end if
         deallocate(sto5)
         deallocate(dum)

         itert=itert+iter
         itotals=itotals+iter
         minkt=minkt+mink
      endif

      call storage_derivatives(0,1)
      if(islord.ne.0) call switch(nmat,nmatb,islord,2,2,nsizea1)
      if(islord.ne.0) call switchb(bp,nrhs,nrhsb,islord,2,2,neq)
      call dual(11)

 999  continue

      return
      end
c.............................................................

      subroutine extract_adiag
c s kelkar 3/20/12. extract diagonal elements of 'a'
      use comai, only : neq
      use combi, only : nelmdg
      use comei, only : a
      use davidi, only : nmat
      implicit none

      integer idisk

      real*8 adiagsk(2*neq)

      do idisk=1,neq
         adiagsk(idisk) = a(nelmdg(idisk)-neq-1+nmat(1))
         adiagsk(idisk+neq) = a(nelmdg(idisk)-neq-1+nmat(4))
      enddo

      return
      end
c...........................................................

      subroutine extract_nodal_equations(idisk)
c s kelkar 3/20/12. extract rows of a corrosponding
c to the node is]disk from the 'a' matrix
      use comai, only : neq, ierr, l
      use combi, only : nelmdg, nelm
      use comei, only : a
      use comgi, only : bp
      use davidi, only : nmat, nrhs
      implicit none

      integer idisk,iblock, jdisk, indsk,ii1,ii2, j,k,iq
      integer neqp1
      integer kbsk(20), inda, kbp

      real*8 aeqsk(20,4)

      neqp1=neq+1

      ii1=nelm(idisk)+1
      ii2=nelm(idisk+1)

      iq=0
      do j=ii1,ii2
         iq=iq+1
         kbsk(iq)=nelm(j)
         do k=1,4
            inda = j-neqp1+nmat(k)
            aeqsk(iq,k)=a(inda)
         enddo
      enddo

      write(ierr,*)' Time step is:', l
      write(ierr,*)' Node I: ', idisk
      write(ierr,*)(kbsk(j),j=1,iq)
 9991 format(30(2x,i10))
      do k=1,4
         kbp=1
         if(k.gt.2)kbp=2
         write(ierr,9992)(aeqsk(j,k),j=1,iq),bp(nrhs(kbp)+idisk)
      enddo
 9992 format(30(2x,g12.5))

      return
      end
c...........................................................

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


    Program dggglm_example

!     DGGGLM Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use blas_interfaces, Only: dnrm2
      Use lapack_interfaces, Only: dggglm
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: rnorm
      Integer :: i, info, lda, ldb, lwork, m, n, p
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), b(:, :), d(:), work(:), x(:), &
        y(:)
!     .. Executable Statements ..
      Write (nout, *) 'DGGGLM Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n, p
      lda = m
      ldb = m
      lwork = n + m + nb*(m+p)
      Allocate (a(lda,n), b(ldb,p), d(m), work(lwork), x(n), y(p))

!     Read A, B and D from data file

      Read (nin, *)(a(i,1:n), i=1, m)
      Read (nin, *)(b(i,1:p), i=1, m)
      Read (nin, *) d(1:m)

!     Solve the weighted least squares problem

!     minimize ||inv(B)*(d - A*x)|| (in the 2-norm)

      Call dggglm(m, n, p, a, lda, b, ldb, d, x, y, work, lwork, info)

!     Print least squares solution, x

      Write (nout, *) 'Weighted least squares solution'
      Write (nout, 100) x(1:n)

!     Print residual vector y = inv(B)*(d - A*x)

      Write (nout, *)
      Write (nout, *) 'Residual vector'
      Write (nout, 110) y(1:p)

!     Compute and print the square root of the residual sum of
!     squares

      rnorm = dnrm2(p, y, 1)

      Write (nout, *)
      Write (nout, *) 'Square root of the residual sum of squares'
      Write (nout, 110) rnorm

100   Format (1X, 7F11.4)
110   Format (3X, 1P, 7E11.2)
    End Program


text:
!==============================================================================!
  subroutine ModOut
!------------------------------------------------------------------------------!
!   Modifies the fluxes at outflow boundaries to conserve the mass             !
!   in the whole domain                                                        !
!----------------------------------[Modules]-----------------------------------!
  use allp_mod, only: convect, inflow, outflow, pressure, tiny
  use all_mod
  use pro_mod
  use Moin_Problem_Mod
!------------------------------------------------------------------------------!
  implicit none
!-----------------------------------[Locals]-----------------------------------!
  integer :: m, s, c1, c2
  real    :: fac(256)
  real    :: delta_mass, RhoS
!---------------------------------[Interface]----------------------------------!
  ! these function must be pushed into par_mod, then interface is not needed
  interface
    subroutine glosum(phi)
      implicit none
      include 'mpif.h'
      real    :: phi
      real    :: phinew
      integer :: error

    end subroutine glosum
  end interface
!==============================================================================!

  !--------------------------------------------------------!
  !   Calculate density change in domain = d rho / dt*dV   !
  !--------------------------------------------------------!

  ! - d rho / dt dv
  delta_mass = 0.0
  do m = 1, Nmat
    delta_mass = sum(volume(1: NC)/dt* &
    (1.5*Rho % n(1: NC) - 2.0*Rho % o(1: NC) + 0.5*Rho % oo(1: NC)))

    call glosum(delta_mass)
  end do

  !--------------------------------------!
  !   Calculate the inflow mass fluxes   !
  !--------------------------------------!

  do m = 1, Nmat
    massin(m) = 0.0
    do s = 1, NS
      c1 = SideC(1,s)
      c2 = SideC(2,s)
      if(c2 < 0) then

        if(TypeBC(c2)  ==  INFLOW) then
          if(material(c1) == m) massin(m) = massin(m) + Flux(s)
        endif
        if(TypeBC(c2)  ==  PRESSURE .and. Flux(s) < 0.0) then
          if(material(c1) == m) massin(m) = massin(m) + Flux(s)
        end if
        if(TypeBC(c2)  ==  CONVECT .and. Flux(s) < 0.0) then
          if(material(c1) == m) massin(m) = massin(m) + Flux(s)
        end if
      end if
    end do
    call glosum(massin(m))
  end do

  !---------------------------------------!
  !   Calculate the outflow mass fluxes   !
  !---------------------------------------!

  do m = 1, Nmat
    masout(m) = 0.0
    do s = 1, NS
      c1 = SideC(1,s)
      c2 = SideC(2,s)
      if(c2  < 0) then
        if(TypeBC(c2) == OUTFLOW) then
          if(material(c1) == m) masout(m) = masout(m) + Flux(s)
        endif
        if(TypeBC(c2) == CONVECT .and. Flux(s) > 0.0) then
          if(material(c1) == m) masout(m) = masout(m) + Flux(s)
        endif
        if(TypeBC(c2) == PRESSURE .and. Flux(s) > 0.0) then
          if(material(c1) == m) masout(m) = masout(m) + Flux(s)
        endif

      endif
    end do
    call glosum(masout(m))
  end do

  !---------------------------------------!
  !   Correct it via fac to satisfy the   !
  !   overall mass balance                !
  !---------------------------------------!
  do m = 1, Nmat
    fac(m) = -(massin(m)+delta_mass)/(masout(m)+TINY)
  end do

  do m = 1, Nmat
    do s = 1, NS
      c1 = SideC(1,s)
      c2 = SideC(2,s)
      if(c2  < 0) then
        if(TypeBC(c2) == OUTFLOW .or. TypeBC(c2) == CONVECT .or. TypeBC(c2) == PRESSURE) then
          if(material(c1) == m) then

            U % n(c2) = rho % n(c1) * U % n(c1) / rho % n(c2) * fac(m)
            V % n(c2) = rho % n(c1) * V % n(c1) / rho % n(c2) * fac(m)
            W % n(c2) = rho % n(c1) * W % n(c1) / rho % n(c2) * fac(m)

            RhoS = f(s) * rho % n(c1) + ( 1.0-f(s) ) * rho % n(c2)

            Flux(s) =   f(s) * ( rho % n(c1) * U % n(c1) * Sx(s) +   &
                                 rho % n(c1) * V % n(c1) * Sy(s) +   &
                                 rho % n(c1) * W % n(c1) * Sz(s) ) + &
                ( 1.0-f(s) ) * ( rho % n(c2) * U % n(c2) * Sx(s) +   &
                                 rho % n(c2) * V % n(c2) * Sy(s) +   &
                                 rho % n(c2) * W % n(c2) * Sz(s) )
            Flux_u(s) = Flux(s) / RhoS

          endif
        endif
      endif
    end do
  end do

end subroutine ModOut

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
I am a Davisites Davisite! I know its sad but true ) I grew up in town, went to Valley Oak Elementary School Valley Oak, Holmes Junior High School Holmes & Davis Senior High School Davis High and went off to LA for awhile only to return with my Significant Others hubby Nate in tow. I have a almost 6 year old daughter and am a stay at home mommy although I do work outside the home doing organic skin services at girlgonegreen skincare and was formerly at  Pure Beauty in South Davis. My passions include parenting, skincare, Bikram Yoga Davis Bikram yoga and eating & shopping locally. I am trying to live as green as possible & I enjoy the challenge!

20080816 15:31:00 nbsp Welcome to the Wiki. Have a few links! Users/JasonAller

20110803 09:00:40 nbsp Im confused... first you say that PureBeauty has relocated, and then you say that its closed  which is it?  If its closed, then we still keep the page around for historical reasons and mark it as departed.  If its moved, we can indicate that, too.  But please dont delete the content, which belongs to the Davis community, from the Wiki. Users/CovertProfessor
    She did the ah, a delete button... Ill use it to delete this one bit of text, and accidentally deleted the whole thing.  It happens now and then.  I got a polite but panicy email that I just read. jw

20110803 09:45:32 nbsp Hey, I just got your email.  No worries... people accidentally delete entries now and then.  Its easy to restore them (CovertProfessor already has).  Although sometimes people think youre trying to remove the whole thing (thus the questions). Users/JabberWokky

20111008 21:09:44 nbsp Hi, Maybe you would like to start your own Business page. I would be happy to help if needed. Its a bit confusing in the Pure Beauty page. :) Users/jsbmeb

20111017 16:12:25 nbsp Did you add it to the the salons page? Users/JabberWokky

20111019 06:16:40 nbsp Hi sorry it took me a long time to reply. Added girlgonegreen to Beauty Salons and Barber Shops :) Users/jsbmeb

20120405 10:07:38 nbsp No problem.  You add thumbnail and how big you want it.  You can also add a caption to each of them.  I added an example caption to the door photo, but feel free to add whatever explains what is in the photo.  Just like all thumbnailed photos, you can click on them to view them full size.

The Image() format can be a little confusing.  If you mess anything up, dont worry about it... somebody will notice and come along to help fix it. Users/JabberWokky


text:
MODULE EulerEquationsLimiterUtilitiesModule_DG

  USE KindModule, ONLY: &
    DP
  USE ProgramHeaderModule, ONLY: &
    nX, nNodesX, nDOFX
  USE UtilitiesModule, ONLY: &
    NodeNumberX
  USE ReferenceElementModuleX, ONLY: &
    WeightsX_q
  USE MeshModule, ONLY: &
    MeshX
  USE PolynomialBasisModule_Lagrange, ONLY: &
    L_X1, L_X2, L_X3
  USE PolynomialBasisMappingModule, ONLY: &
    MapNodalToModal_Fluid
  USE FluidFieldsModule, ONLY: &
    uCF, iCF_S1, iCF_S2, iCF_S3, nCF, &
    Shock

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: InitializeDiscontinuityDetector
  PUBLIC :: DetectDiscontinuities
  PUBLIC :: SolveTheta_Bisection

  LOGICAL, PARAMETER :: UseDiscontinuityDetector = .TRUE.
  REAL(DP), DIMENSION(:,:), ALLOCATABLE :: Lagrange_X1_P
  REAL(DP), DIMENSION(:,:), ALLOCATABLE :: Lagrange_X1_M

CONTAINS

  SUBROUTINE InitializeDiscontinuityDetector

    INTEGER  :: iNodeX1, iNodeX2, iNodeX3, iNodeX
    INTEGER  :: jNodeX1, jNodeX2, jNodeX3, jNodeX
    REAL(DP) :: eta_X1, eta_X2, eta_X3

    ALLOCATE( Lagrange_X1_P(nDOFX,nDOFX) )
    ALLOCATE( Lagrange_X1_M(nDOFX,nDOFX) )

    DO jNodeX3 = 1, nNodesX(3)
      DO jNodeX2 = 1, nNodesX(2)
        DO jNodeX1 = 1, nNodesX(1)

          jNodeX = NodeNumberX( jNodeX1, jNodeX2, jNodeX3 )

          eta_X1 = MeshX(1) % Nodes(jNodeX1)
          eta_X2 = MeshX(2) % Nodes(jNodeX2)
          eta_X3 = MeshX(3) % Nodes(jNodeX3)

          DO iNodeX3 = 1, nNodesX(3)
            DO iNodeX2 = 1, nNodesX(2)
              DO iNodeX1 = 1, nNodesX(1)

                iNodeX = NodeNumberX( iNodeX1, iNodeX2, iNodeX3 )

                Lagrange_X1_P(iNodeX,jNodeX) &
                  = L_X1(iNodeX1) % P( eta_X1 + 1.0_DP ) &
                    * L_X2(iNodeX2) % P( eta_X2 ) &
                      * L_X3(iNodeX3) % P( eta_X3 )

                Lagrange_X1_M(iNodeX,jNodeX) &
                  = L_X1(iNodeX1) % P( eta_X1 - 1.0_DP ) &
                    * L_X2(iNodeX2) % P( eta_X2 ) &
                      * L_X3(iNodeX3) % P( eta_X3 )

              END DO
            END DO
          END DO

        END DO
      END DO
    END DO

  END SUBROUTINE InitializeDiscontinuityDetector

  SUBROUTINE DetectDiscontinuities

    INTEGER  :: iX1, iX2, iX3, iCF, k
    REAL(DP) :: F_M, F_P
    REAL(DP), PARAMETER :: alpha = 1.5_DP
    REAL(DP), DIMENSION(nCF)       :: uCF_A
    REAL(DP), DIMENSION(nDOFX,nCF) :: uCF_M
    REAL(DP), DIMENSION(nDOFX,nCF) :: uCF_M_P_X1
    REAL(DP), DIMENSION(nDOFX,nCF) :: uCF_M_N_X1

    Shock(:,:,:) = 0.0_DP

    IF( .NOT. UseDiscontinuityDetector )THEN

      Shock(:,:,:) = 1.0_DP
      RETURN

    END IF

    DO iX3 = 1, nX(3)
      DO iX2 = 1, nX(2)
        !$OMP PARALLEL DO PRIVATE &
        !$OMP&              ( iX1, iCF, k, F_M, F_P, uCF_A, &
        !$OMP&                uCF_M_P_X1, uCF_M, uCF_M_N_X1 )
        DO iX1 = 1, nX(1)

          ! --- Map To Modal Representation ---

          DO iCF = 1, nCF

            CALL MapNodalToModal_Fluid &
                   ( uCF(:,iX1-1,iX2,iX3,iCF), uCF_M_P_X1(:,iCF) )
            CALL MapNodalToModal_Fluid &
                   ( uCF(:,iX1,  iX2,iX3,iCF), uCF_M     (:,iCF) )
            CALL MapNodalToModal_Fluid &
                   ( uCF(:,iX1+1,iX2,iX3,iCF), uCF_M_N_X1(:,iCF) )

          END DO

          ! --- Cell-Averaged Quantities ---

          uCF_A(1:nCF) = uCF_M(1,1:nCF)

          ! --- Detect Discontinuity with Harten's Sub-Cell Method ---

          LOOP: DO iCF = 1, nCF

            IF( ANY( iCF == [ iCF_S1, iCF_S2, iCF_S3 ] ) ) CYCLE LOOP

            F_M = 0.0_DP
            F_P = 0.0_DP
            DO k = 1, nDOFX

              F_M &
                = F_M &
                  + WeightsX_q(k) &
                    * SUM( Lagrange_X1_M(:,k) * uCF(:,iX1+1,iX2,iX3,iCF) )

              F_P &
                = F_P &
                  + WeightsX_q(k) &
                    * SUM( Lagrange_X1_P(:,k) * uCF(:,iX1-1,iX2,iX3,iCF) )

            END DO
            F_M = F_M - uCF_A(iCF)
            F_P = F_P - uCF_A(iCF)

            IF( F_M * F_P <= 0.0_DP )THEN
              DO k = 1, MIN( 2, nDOFX )
                IF( ABS( uCF_M(k,iCF) ) &
                    > alpha * ABS( uCF_M_P_X1(k,iCF) ) .AND. &
                    ABS( uCF_M(k,iCF) ) &
                    > alpha * ABS( uCF_M_N_X1(k,iCF) ) ) &
                THEN
                  Shock(iX1,iX2,iX3) = 1.0_DP
                  EXIT LOOP
                END IF
              END DO
            END IF

          END DO LOOP

        END DO
        !$OMP END PARALLEL DO
      END DO
    END DO

  END SUBROUTINE DetectDiscontinuities

  SUBROUTINE SolveTheta_Bisection &
    ( D_K, S1_K, S2_K, S3_K, E_K, dD, dS1, dS2, dS3, dE, Tol_E, Theta )

    REAL(DP), INTENT(in)  :: D_K, S1_K, S2_K, S3_K, E_K
    REAL(DP), INTENT(in)  :: dD,  dS1,  dS2,  dS3,  dE, Tol_E
    REAL(DP), INTENT(out) :: Theta

    LOGICAL  :: Converged
    INTEGER  :: Iter
    REAL(DP) :: a, b, c, ab, Phi_a, Phi_b, Phi_c, Phi_0
    INTEGER,  PARAMETER :: MaxIter = 128
    REAL(DP), PARAMETER :: Tol_ab  = 1.0d-8
    REAL(DP), PARAMETER :: Tol_Phi = 1.0d-8

    a = 0.0_DP
    Phi_a &
      = ThetaFun &
          ( D_K, S1_K, S2_K, S3_K, E_K, dD, dS1, dS2, dS3, dE, Tol_E, a )

    b = 1.0_DP
    Phi_b &
      = ThetaFun &
          ( D_K, S1_K, S2_K, S3_K, E_K, dD, dS1, dS2, dS3, dE, Tol_E, b )

    Phi_0 = Phi_a
    ab    = b - a

    Converged = .FALSE.
    Iter      = 0
    DO WHILE ( .NOT. Converged )

      Iter = Iter + 1

      ab = 0.5_DP * ab
      c  = a + ab

      Phi_c  &
        = ThetaFun &
            ( D_K, S1_K, S2_K, S3_K, E_K, dD, dS1, dS2, dS3, dE, Tol_E, c )

      IF( Phi_a * Phi_c < 0.0_DP )THEN

        b     = c
        Phi_b = Phi_c

      ELSE

        a     = c
        Phi_a = Phi_c

      END IF

      IF( ab < Tol_ab .AND. ABS( Phi_a ) / Phi_0 < Tol_Phi ) &
        Converged = .TRUE.

      IF( Iter > MaxIter .AND. .NOT. Converged )THEN
        WRITE(*,*)
        WRITE(*,'(A4,A)') &
          '', 'SolveTheta_Bisection (Euler Equations)'
        WRITE(*,'(A6,A21,I4.4,A11)') &
          '', 'No Convergence After ', Iter, ' Iterations'
        WRITE(*,*)
      END IF

    END DO

    Theta = a

  END SUBROUTINE SolveTheta_Bisection

  PURE REAL(DP) FUNCTION ThetaFun &
    ( D_K, S1_K, S2_K, S3_K, E_K, dD, dS1, dS2, dS3, dE, Tol_E, Theta )

    REAL(DP), INTENT(in) :: D_K, S1_K, S2_K, S3_K, E_K
    REAL(DP), INTENT(in) :: dD, dS1, dS2, dS3, dE, Tol_E, Theta

    ThetaFun &
      = E_K+Theta*dE &
        - 0.5_DP*DOT_PRODUCT &
                   ( [S1_K+Theta*dS1,S2_K+Theta*dS2,S3_K+Theta*dS3], &
                     [S1_K+Theta*dS1,S2_K+Theta*dS2,S3_K+Theta*dS3] ) &
          / ( D_K+Theta*dD ) - Tol_E

    RETURN
  END FUNCTION ThetaFun

END MODULE EulerEquationsLimiterUtilitiesModule_DG

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


text:
64 megs buffer: work

: block  #2 rshift work + ;
: blocks  #2 rshift ;
: >slotva  block ;

: reserve  dup constant over constant + ;

0
4 reserve #classes classes       \ 128 x 32bytes
256 reserve #graphics graphics   \ describes images, sprite sheets, and tilesets
4096 reserve #tilemaps tilemaps  \ a tilemap is one block (!6x16 tiles)
256 reserve #sounds sounds       \
4096 reserve #actors actors      \ used for templates as well as gameplay
2048 reserve #dialogs dialogs    \ used for templates as well as gameplay.  point to files, but
                                 \ could also store some state.

cr . .( blocks used)

: graphic  graphics + block ;
: sound  sounds + block ;
: tilemap  tilemaps + block ;

: class  #32 * classes block + ;  \ class 0 is reserved for some variables
: class>name  count ;
: class>details  #28 + @ ;
: class>script  class>name  s" prg/gamester/" s[   +S  s" .f" +S  ]s ;

\ : block-variable  ( block# offset - block# offset+cell )
\     2dup swap block + constant cell+ ;
\
\ classes 0
\     block-variable classSlot
\ 2drop

0 value (class#)

: export  ( class - )
    (class#) class #28 + ! ;

: load ( class# - )
    dup to (class#)
    class
        class>script included ;

: path  s" prg/gamester/test.blk" ;
: write  blocks swap block swap path file! ;
: flush  work 64 megs path file! ;
: cache
    path work 64 megs @file
    classes >slotva @ for
        i 1 +  load
    loop
;

\ '\' as a prefix means it's an "Edit time command"

: \class  ( - <name> )
    1 classes >slotva +!
    bl parse classes >slotva @ dup >r  class place
    r> load
;

: table  create cells /allot does> swap cells + ;

: bclass: struct: lastbody export ;

cache

max_stars_repo_name: RogerLevy/gamester

max_issues_repo_name: RogerLevy/blocks

max_forks_repo_name: RogerLevy/blocks


text:
SUBMODULE (vtk_io) vtk_io_implementation
    !! author: Ian Porter
    !! date: 12/1/2017
    !!
    !! This module implements the ability to read/write VTK formatted files
    !!

    LOGICAL,      SAVE :: file_was_already_open     = .FALSE.
    LOGICAL,      SAVE :: printed_cell_data_header  = .FALSE.
    LOGICAL,      SAVE :: printed_point_data_header = .FALSE.
    INTEGER(i4k), SAVE :: newunit
    CHARACTER(LEN=:), ALLOCATABLE, SAVE :: form

    CONTAINS

        MODULE PROCEDURE vtk_legacy_full_write
        USE Misc,     ONLY : to_uppercase
        USE vtk_vars, ONLY : default_fn, default_title, vtkfilename, vtktitle, version, fcnt, file_extension
        USE XML,      ONLY : convert_format_to_string, file_format_text, file_format, ascii, binary, format_ascii, format_binary
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/1/2017
        !!
        !! This subroutine writes the legacy vtk output file
        !!
        INTEGER(i4k) :: i, inputstat
        CHARACTER(LEN=11) :: fm

        ! Clear out any pre-existing data
        IF (ALLOCATED(vtkfilename))      DEALLOCATE(vtkfilename)
        IF (ALLOCATED(form))             DEALLOCATE(form)
        IF (ALLOCATED(vtktitle))         DEALLOCATE(vtktitle)
        IF (ALLOCATED(file_format_text)) DEALLOCATE(file_format_text)

        IF (PRESENT(format)) THEN
            file_format = format                                !! Calling program provided what file type to use for vtk file
        ELSE
            file_format = ascii                                 !! Default to ascii
        END IF
        IF (PRESENT(filename)) THEN
            ALLOCATE(vtkfilename, source=filename)              !! Calling program provided a filename
        ELSE
            ALLOCATE(vtkfilename, source=default_fn)            !! Calling program did not provide a filename. Use default
        END IF
        IF (PRESENT(multiple_io)) THEN
            IF (multiple_io) THEN
                mio_filename: BLOCK
                    CHARACTER(LEN=8) :: fcnt_char = ''          !! File count character
                    CHARACTER(LEN=:), ALLOCATABLE :: base_fn    !! Base file name
                    WRITE (fcnt_char,FMT='(i8)') fcnt
                    ALLOCATE(base_fn, source=vtkfilename(1:INDEX(to_uppercase(vtkfilename),to_uppercase(file_extension))-1))
                    DEALLOCATE(vtkfilename)
                    ALLOCATE(vtkfilename, source=base_fn // "_" // TRIM(ADJUSTL(fcnt_char)) // file_extension)
                    fcnt = fcnt + 1                             !! Increase timestep file counter by 1
                END BLOCK mio_filename
            END IF
        END IF

        IF (PRESENT(title)) THEN
            ALLOCATE(vtktitle, source=title)                    !! Calling program provided a title
        ELSE
            ALLOCATE(vtktitle, source=default_title)            !! Calling program did not provide a title. Use default
        END IF

        IF (PRESENT(unit)) THEN
            newunit = unit
            INQUIRE(unit=newunit, opened=file_was_already_open) !! Check to see if file is already open
            IF (.NOT. file_was_already_open) THEN               !! File is not yet open. Determine format from file_format
                SELECT CASE (file_format)
                CASE (ascii)
                    ALLOCATE(form, source='formatted')
                CASE (binary)
                    ALLOCATE(form, source='unformatted')
                CASE DEFAULT
                    ERROR STOP 'Error: file_format is incorrectly defined in vtk_legacy_full_write'
                END SELECT
                ALLOCATE(file_format_text, source=convert_format_to_string(file_format))
                OPEN(unit=newunit, file=vtkfilename, iostat=inputstat, status='REPLACE', form=form)
                                                                !! Open the VTK file
            ELSE                                                !! File is already open. Determine format based on file format
                INQUIRE(unit=newunit,form=fm)
                SELECT CASE (TO_UPPERCASE(TRIM(fm)))
                CASE ('FORMATTED')
                    ALLOCATE(file_format_text, source=format_ascii)
                CASE DEFAULT
                    ALLOCATE(file_format_text, source=format_binary)
                END SELECT
            END IF
        ELSE
            !! No unit # provided. Make determination by value set for file_format
            SELECT CASE (file_format)
            CASE (ascii)
                ALLOCATE(form, source='formatted')
            CASE (binary)
                ALLOCATE(form, source='unformatted')
            CASE DEFAULT
                WRITE(*,*) 'Warning: file_format is incorrectly defined. Will default to ASCII'
                ALLOCATE(form, source='formatted')
            END SELECT
            ALLOCATE(file_format_text, source=convert_format_to_string(file_format))
            OPEN(newunit=newunit, file=vtkfilename, iostat=inputstat, status='REPLACE', form=form)
                                                                !! Open the VTK file
        END IF

        WRITE(newunit,100) version                              !! VTK version (currently, 3.0)
        WRITE(newunit,100) vtktitle                             !! VTK title card
        WRITE(newunit,100) file_format_text                     !! VTK file type

        CALL geometry%write(newunit)                            !! Write the geometry information

        IF (PRESENT(celldatasets)) THEN
            WRITE(newunit,101) celldatasets(1)%attribute%nvals
            DO i = 1, SIZE(celldatasets)
                CALL celldatasets(i)%attribute%write(newunit)   !! Write the cell data values
            END DO
        ELSE IF (PRESENT(celldata)) THEN
            WRITE(newunit,101) celldata%nvals
            CALL celldata%write(newunit)                        !! Write the cell data values
        END IF

        IF (PRESENT(pointdatasets)) THEN
            WRITE(newunit,102) pointdatasets(1)%attribute%nvals
            DO I = 1, SIZE(pointdatasets)
                CALL pointdatasets(i)%attribute%write(newunit)
            END DO
        ELSE IF (PRESENT(pointdata)) THEN
            WRITE(newunit,102) pointdata%nvals
            CALL pointdata%write(newunit)                       !! Write the point data values
        END IF

        IF (ANY([ PRESENT(celldatasets), PRESENT(celldata), PRESENT(pointdatasets), PRESENT(pointdata) ])) THEN
            CALL vtk_legacy_finalize (finished=.TRUE.)          !! Full legacy write w/ data. Close file.
        ELSE
            CALL vtk_legacy_finalize (finished=.FALSE.)         !! No data was provided, only geometry info. Do not close file.
        END IF

100     FORMAT(a)
101     FORMAT('CELL_DATA ',i0)
102     FORMAT('POINT_DATA ',i0)

        END PROCEDURE vtk_legacy_full_write

        MODULE PROCEDURE vtk_legacy_append
        USE Misc,     ONLY : to_uppercase
        USE vtk_vars, ONLY : default_fn, default_title, version, file_extension
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/1/2017
        !!
        !! This subroutine writes the legacy vtk output file
        !!
        INTEGER(i4k) :: i, inputstat
        LOGICAL :: file_is_still_open

        INQUIRE(unit=newunit,opened=file_is_still_open)
        IF (.NOT. file_is_still_open) THEN
            !! For some reason, file was closed. Re-open the file to a new unit
            OPEN(newunit=newunit, file=vtkfilename, iostat=inputstat, status='REPLACE', form=form, position='APPEND')
        END IF

        IF (PRESENT(celldatasets)) THEN
            IF (.NOT. printed_cell_data_header) THEN
                WRITE(newunit,101) celldatasets(1)%attribute%nvals
            END IF
            DO i = 1, SIZE(celldatasets)
                CALL celldatasets(i)%attribute%write(newunit)   !! Write the cell data values
            END DO
        ELSE IF (PRESENT(celldata)) THEN
            IF (.NOT. printed_cell_data_header) THEN
                WRITE(newunit,101) celldata%nvals
                printed_cell_data_header = .TRUE.
            END IF
            CALL celldata%write(newunit)                        !! Write the cell data values
        END IF

        IF (PRESENT(pointdatasets)) THEN
            IF (.NOT. printed_point_data_header) THEN
                WRITE(newunit,102) pointdatasets(1)%attribute%nvals
            END IF
            DO I = 1, SIZE(pointdatasets)
                CALL pointdatasets(i)%attribute%write(newunit)
            END DO
        ELSE IF (PRESENT(pointdata)) THEN
            IF (.NOT. printed_point_data_header) THEN
                WRITE(newunit,102) pointdata%nvals
                printed_point_data_header = .TRUE.
            END IF

            CALL pointdata%write(newunit)                       !! Write the point data values
        END IF

101     FORMAT('CELL_DATA ',i0)
102     FORMAT('POINT_DATA ',i0)

        END PROCEDURE vtk_legacy_append

        MODULE PROCEDURE vtk_legacy_finalize
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 06/03/2019
        !!
        !! This subroutine is a finalizer for the vtk file write
        !!

        IF (finished) THEN
            IF (.NOT. file_was_already_open) THEN
                CLOSE(newunit)                                  !! Close the VTK file if file was not open prior to calling vtkmofo
            END IF
            printed_cell_data_header  = .FALSE.
            printed_point_data_header = .FALSE.
        END IF

        END PROCEDURE vtk_legacy_finalize

        MODULE PROCEDURE vtk_legacy_read
        USE Misc,     ONLY : def_len
        USE vtk_vars, ONLY : default_fn, default_title, vtkfilename, vtktitle, version
        USE XML,      ONLY : convert_string_to_format, ascii, binary, file_format, file_format_text
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/20/2017
        !!
        !! This subroutine reads the legacy vtk output file
        !!
        INTEGER(i4k) :: i, inputstat
        LOGICAL      :: file_is_open
        CHARACTER(LEN=:), ALLOCATABLE :: form, vtk_version
        CHARACTER(LEN=def_len) :: line

        INQUIRE(unit = unit, opened = file_is_open)        !! Check to see if file is already open
        IF (.NOT. file_is_open) THEN                       !! File is not yet open. Determine format to open file
            OPEN(unit=unit, file=filename, iostat=inputstat, status='old')
                                                           !! Open the VTK file
        END IF

        READ(unit,100) line                                !! VTK version (currently, 3.0)
        ALLOCATE(vtk_version, source=TRIM(line))
        line = ''
        IF (vtk_version /= version) THEN
                                                           !! If the file is not version 3.0, abort read
            ERROR STOP 'ERROR: vtkmofo only works with version 3.0'
        END IF

        READ(unit,100) title                               !! VTK title card
        READ(unit,100) line                                !! VTK file type
        ALLOCATE(file_format_text, source=TRIM(line))
        line = ''

        CLOSE(unit)                                        !! Close the file to re-open it in the proper format

        file_format = convert_string_to_format(file_format_text)
        SELECT CASE (file_format)
        CASE (ascii)
            ALLOCATE(form, source='formatted')
        CASE (binary)
            ALLOCATE(form, source='unformatted')
        CASE DEFAULT
            ERROR STOP 'Unsupported file type. Must be ASCII or BINARY. Terminated in vtk_legacy_read'
        END SELECT

        OPEN(unit=unit, file=filename, iostat=inputstat, status='old', form=form)
                                                           !! Open the VTK file in the proper format
        READ(unit,*)                                       !! Skip over this line
        READ(unit,*)                                       !! Skip over this line
        READ(unit,*)                                       !! Skip over this line

        CALL geometry%read(unit)                           !! Read the information from the file

        IF (PRESENT(celldatasets)) THEN
            READ(unit,*)
            DO i = 1, SIZE(celldatasets)
                CALL celldatasets(i)%attribute%read(unit)  !! Write the cell data values
            END DO
        ELSE IF (PRESENT(celldata)) THEN
            READ(unit,*)
            CALL celldata%write(unit)                      !! Write the cell data values
        END IF
        IF (PRESENT(pointdatasets)) THEN
            READ(unit,*)
            DO I = 1, SIZE(pointdatasets)
                CALL pointdatasets(i)%attribute%read(unit)
            END DO
        ELSE IF (PRESENT(pointdata)) THEN
            READ(unit,*)
            CALL pointdata%write(unit)                     !! Write the point data values
        END IF

        CLOSE(unit)                                        !! Close the VTK file

        ALLOCATE(vtkfilename, source=filename)             !! Save the filename for future internal use
        ALLOCATE(vtktitle, source=title)                   !! Save the title for future internal use

100     FORMAT(a)

        END PROCEDURE vtk_legacy_read

        MODULE PROCEDURE vtk_serial_full_write
        USE vtk_datasets,    ONLY : struct_pts, struct_grid, rectlnr_grid, polygonal_data, unstruct_grid
        USE VTK_serial_file, ONLY : serial_file
        USE VTK_serial_Grid, ONLY : VTK_serial_RectilinearGrid_dt, VTK_serial_StructuredGrid_dt, &
          &                         VTK_serial_UnstructuredGrid_dt, VTK_serial_ImageData_dt
        USE XML,             ONLY : file_format, file_format_text, convert_format_to_string, ascii, format_ascii
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 5/08/2019
        !!
        !! This subroutine writes the modern serial vtk output file
        !!

        ! Clear out any pre-existing data
        IF (ALLOCATED(vtkfilename))      DEALLOCATE(vtkfilename)
        IF (ALLOCATED(form))             DEALLOCATE(form)
        IF (ALLOCATED(file_format_text)) DEALLOCATE(file_format_text)

        IF (PRESENT(format)) THEN
            ALLOCATE(file_format_text, source=convert_format_to_string(format))
            file_format = format
        ELSE
            ALLOCATE(file_format_text,source=format_ascii)     !! Default to binary
            file_format = ascii
        END IF
        IF (PRESENT(filename)) THEN
            ALLOCATE(vtkfilename, source=filename)              !! Calling program provided a filename
        ELSE
            ALLOCATE(vtkfilename, source=default_fn)            !! Calling program did not provide a filename. Use default
        END IF

        IF (PRESENT(multiple_io)) THEN
            IF (multiple_io) THEN
                mio_filename: BLOCK
                    CHARACTER(LEN=8) :: fcnt_char = ''          !! File count character
                    CHARACTER(LEN=:), ALLOCATABLE :: base_fn    !! Base file name
                    WRITE (fcnt_char,FMT='(i8)') fcnt
                    ALLOCATE(base_fn, source=vtkfilename)
                    DEALLOCATE(vtkfilename)
                    ALLOCATE(vtkfilename, source=base_fn // "_" // TRIM(ADJUSTL(fcnt_char)))
                    fcnt = fcnt + 1                             !! Increase timestep file counter by 1
                END BLOCK mio_filename
            END IF
        END IF

        ALLOCATE(serial_file)

        SELECT TYPE (geometry)
        CLASS IS (struct_pts)
            ALLOCATE(VTK_serial_ImageData_dt::serial_file%vtk_dataset)
        CLASS IS (struct_grid)
            ALLOCATE(VTK_serial_StructuredGrid_dt::serial_file%vtk_dataset)
        CLASS IS (rectlnr_grid)
            ALLOCATE(VTK_serial_RectilinearGrid_dt::serial_file%vtk_dataset)
        CLASS IS (polygonal_data)
            ERROR STOP 'Procedure not yet implemented for: POLYGONAL GRID. Termination in subroutine: vtk_serial_full_write'
        CLASS IS (unstruct_grid)
            ALLOCATE(VTK_serial_UnstructuredGrid_dt::serial_file%vtk_dataset)
        CLASS DEFAULT
            ERROR STOP 'Unsupported geometry type. Termination in subroutine: vtk_serial_full_write'
        END SELECT

        CALL serial_file%vtk_dataset%set_grid(geometry)

        CALL serial_file%setup(filename=vtkfilename // TRIM(serial_file%vtk_dataset%file_extension))
        !! Append data
        CALL vtk_serial_append (celldata, pointdata, celldatasets, pointdatasets)
        !! Finalize the write
        IF (ANY([ PRESENT(celldatasets), PRESENT(celldata), PRESENT(pointdatasets), PRESENT(pointdata) ])) THEN
            CALL vtk_serial_finalize (finished=.TRUE.)          !! Full legacy write w/ data. Close file.
        ELSE
            CALL vtk_serial_finalize (finished=.FALSE.)         !! No data was provided, only geometry info. Do not close file.
        END IF

        END PROCEDURE vtk_serial_full_write

        MODULE PROCEDURE vtk_serial_append
        USE VTK_serial_file, ONLY : serial_file
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 06/24/2019
        !!
        !! This subroutine appends data to the legacy vtk output file
        !!

        IF (.NOT. ALLOCATED(serial_file%vtk_dataset%piece)) ALLOCATE(serial_file%vtk_dataset%piece)
        CALL serial_file%vtk_dataset%piece%add_data(celldata, pointdata, celldatasets, pointdatasets)

        END PROCEDURE vtk_serial_append

        MODULE PROCEDURE vtk_serial_finalize
        USE VTK_serial_file, ONLY : serial_file
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 06/24/2019
        !!
        !! This subroutine is a finalizer for the legacy vtk file write
        !!

        IF (finished) THEN
            CALL serial_file%vtk_dataset%finalize()
                  !! This should write everything inside of the piece
            CALL serial_file%add(serial_file%vtk_dataset)
            CALL serial_file%write()
            CALL serial_file%close_file()                    !! Close the VTK file
            CALL serial_file%me_deallocate()
            IF (ALLOCATED(serial_file)) DEALLOCATE(serial_file)
        END IF

        END PROCEDURE vtk_serial_finalize

END SUBMODULE vtk_io_implementation

max_stars_repo_name: zbeekman/vtkmofo

max_issues_repo_name: zbeekman/vtkmofo

max_forks_repo_name: zbeekman/vtkmofo


SUBROUTINE glob_to_axial(axial,global,coord)
!
! This subroutine transforms the global end reactions
! into an axial force for rod elements (2- or 3-d).
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::global(:),coord(:,:)
 REAL(iwp),INTENT(OUT)::axial
 REAL(iwp)::add,ell,zero=0.0_iwp
 INTEGER::ndim,i
 ndim=UBOUND(coord,2)
 add=zero
 DO i=1,ndim
   add=add+(coord(2,i)-coord(1,i))**2
 END DO
 ell=SQRT(add)
 axial=zero
 DO i=1,ndim
   axial=axial+(coord(2,i)-coord(1,i))/ell*global(ndim+i)
 END DO
RETURN
END SUBROUTINE glob_to_axial


text:
      SUBROUTINE FREQCY(FMATRX,FREQ,CNORML,REDMAS,TRAVEL,EORC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION FMATRX(*), REDMAS(*), FREQ(*), CNORML(*), TRAVEL(*)
      LOGICAL EORC
*********************************************************************
*
*  FRCE CALCULATES THE FORCE CONSTANTS AND VIBRATIONAL FREQUENCIES
*       FOR A MOLECULE.  IT USES THE ISOTOPIC MASSES TO WEIGHT THE
*       FORCE MATRIX
*
* ON INPUT   FMATRX   =  FORCE MATRIX, OF SIZE NUMAT*3*(NUMAT*3+1)/2.
*
*********************************************************************
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /ATMASS/ ATMASS(NUMATM)
      COMMON /SCRACH/ OLDF(MAXHES), DUMMY(MAXPAR**2-MAXHES)
      DIMENSION WTMASS(MAXPAR), SHIFT(6)
      DATA FACT/6.023D23/
C
C    CONVERSION FACTOR FOR SPEED OF LIGHT AND 2 PI.
C
      C2PI=1.D0/(2.998D10*3.141592653598D0*2.D0)
C NOW TO CALCULATE THE VIBRATIONAL FREQUENCIES
C
C   FIND CONVERSION CONSTANTS FOR MASS WEIGHTED SYSTEM
      N3=3*NUMAT
      L=0
      DO 10 I=1,NUMAT
         WEIGHT=1.4142136D0/SQRT(ATMASS(I))
         DO 10 J=1,3
            L=L+1
   10 WTMASS(L)=WEIGHT
C    CONVERT TO MASS WEIGHTED FMATRX
      LINEAR=0
      DO 20 I=1,N3
         DO 20 J=1,I
            LINEAR=LINEAR+1
            OLDF(LINEAR)=  FMATRX(LINEAR)*1.D5
   20 FMATRX(LINEAR)=FMATRX(LINEAR)*WTMASS(I)*WTMASS(J)
C
C    1.D5 IS TO CONVERT FROM MILLIDYNES/ANGSTROM TO DYNES/CM.
C
C    DIAGONALIZE
      CALL FRAME(FMATRX,NUMAT,1, SHIFT)
      CALL RSP(FMATRX,N3,N3,FREQ,CNORML)
      DO 30 I=1,N3
         J=(FREQ(I)+50.D0)*0.01D0
   30 FREQ(I)=FREQ(I)-J*100
      DO 40 I=1,N3
   40 FREQ(I)=FREQ(I)*1.D5
C
C    CALCULATE REDUCED MASSES, STORE IN REDMAS
C
      DO 80 I=1,N3
         II=(I-1)*N3
         SUM=0.D0
         DO 70 J=1,N3
            JII=J+II
            JJ=(J*(J-1))/2
            DO 50 K=1,J
   50       SUM=SUM+CNORML(JII)*OLDF(JJ+K)*CNORML(K+II)
            DO 60 K=J+1,N3
   60       SUM=SUM+CNORML(JII)*OLDF((K*(K-1))/2+J)*CNORML(K+II)
   70    CONTINUE
         SUM1=SUM*2.D0
         IF(ABS(FREQ(I)).GT.ABS(SUM)*1.D-20) THEN
            SUM=1.D0*SUM/FREQ(I)
         ELSE
            SUM=0.D0
         ENDIF
         FREQ(I)=SIGN(SQRT(FACT*ABS(FREQ(I)))*C2PI,FREQ(I))
         IF(ABS(FREQ(I)).LT.ABS(SUM1)*1.D+20) THEN
            SUM1=SQRT(ABS(FREQ(I)/(SUM1*1.D-5)))
         ELSE
            SUM1=0.D0
         ENDIF
         IF(SUM.LT.0.D0.OR.SUM.GT.100)SUM=0.D0
C
C 0.0063024=SQRT(2*A*B*C/N) WHERE
C         A=1.196D8 = CONVERSION OF CM**(-1) TO (ERGS = DYNE.ANGSTROMS)
C         B=1000.0  = MILLIDYNES TO DYNES
C         C=1.D8    = CENTIMETERS TO ANGSTROMS
C         N=6.02205D23 = AVOGADRO'S NUMBER
         TRAVEL(I)=SUM1*0.0063024D0
         IF(TRAVEL(I).GT.1.D0)TRAVEL(I)=0.D0
C#      WRITE(6,*)TRAVEL(I)
   80 REDMAS(I)=SUM
      IF(EORC) THEN
C
C    SWITCH EIGENVALUES TO FREQUENCIES
C
C    CONVERT NORMAL VECTORS TO CARTESIAN COORDINATES
C    AND NORMALIZE SO THAT THE TOTAL MOVEMENT IS 1.0 ANGSTROM.
C
         IJ=0
         DO 120 I=1,N3
            SUM=0.D0
            J=0
            DO 100 JJ=1,NUMAT
               SUM1=0.D0
               DO 90 J1=1,3
                  J=J+1
                  IJ=IJ+1
                  CNORML(IJ)=CNORML(IJ)*WTMASS(J)
   90          SUM1=SUM1+CNORML(IJ)**2
  100       SUM=SUM+SQRT(SUM1)
            SUM=1.D0/SUM
            IJ=IJ-N3
            DO 110 J=1,N3
               IJ=IJ+1
  110       CNORML(IJ)=CNORML(IJ)*SUM
  120    CONTINUE
C
C          RETURN HESSIAN IN MILLIDYNES/ANGSTROM IN FMATRX
C
         DO 130 I=1,LINEAR
  130    FMATRX(I)=OLDF(I)*1.D-5
      ELSE
C
C  RETURN HESSIAN AS MASS-WEIGHTED FMATRIX
         LINEAR=0
C
         DO 140 I=1,N3
            DO 140 J=1,I
               LINEAR=LINEAR+1
  140    FMATRX(LINEAR)=OLDF(LINEAR)*1.D-5*WTMASS(I)*WTMASS(J)
      ENDIF
      RETURN
      END

ext: for

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive