text:
      SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
     $                   BETA, C, LDC )
*     .. SCALAR ARGUMENTS ..
      CHARACTER*1        TRANSA, TRANSB
      INTEGER            M, N, K, LDA, LDB, LDC
      DOUBLE PRECISION   ALPHA, BETA
*     .. ARRAY ARGUMENTS ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
*     ..
*
*  PURPOSE
*  =======
*
*  DGEMM  PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS
*
*     C := ALPHA*OP( A )*OP( B ) + BETA*C,
*
*  WHERE  OP( X ) IS ONE OF
*
*     OP( X ) = X   OR   OP( X ) = X',
*
*  ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A )
*  AN M BY K MATRIX,  OP( B )  A  K BY N MATRIX AND  C AN M BY N MATRIX.
*
*  PARAMETERS
*  ==========
*
*  TRANSA - CHARACTER*1.
*           ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSA = 'N' OR 'N',  OP( A ) = A.
*
*              TRANSA = 'T' OR 'T',  OP( A ) = A'.
*
*              TRANSA = 'C' OR 'C',  OP( A ) = A'.
*
*           UNCHANGED ON EXIT.
*
*  TRANSB - CHARACTER*1.
*           ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSB = 'N' OR 'N',  OP( B ) = B.
*
*              TRANSB = 'T' OR 'T',  OP( B ) = B'.
*
*              TRANSB = 'C' OR 'C',  OP( B ) = B'.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY,  M  SPECIFIES  THE NUMBER  OF ROWS  OF THE  MATRIX
*           OP( A )  AND OF THE  MATRIX  C.  M  MUST  BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY,  N  SPECIFIES THE NUMBER  OF COLUMNS OF THE MATRIX
*           OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE
*           AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  K      - INTEGER.
*           ON ENTRY,  K  SPECIFIES  THE NUMBER OF COLUMNS OF THE MATRIX
*           OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K MUST
*           BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - DOUBLE PRECISION.
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS
*           K  WHEN  TRANSA = 'N' OR 'N',  AND IS  M  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSA = 'N' OR 'N',  THE LEADING  M BY K
*           PART OF THE ARRAY  A  MUST CONTAIN THE MATRIX  A,  OTHERWISE
*           THE LEADING  K BY M  PART OF THE ARRAY  A  MUST CONTAIN  THE
*           MATRIX A.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSA = 'N' OR 'N' THEN
*           LDA MUST BE AT LEAST  MAX( 1, M ), OTHERWISE  LDA MUST BE AT
*           LEAST  MAX( 1, K ).
*           UNCHANGED ON EXIT.
*
*  B      - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS
*           N  WHEN  TRANSB = 'N' OR 'N',  AND IS  K  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSB = 'N' OR 'N',  THE LEADING  K BY N
*           PART OF THE ARRAY  B  MUST CONTAIN THE MATRIX  B,  OTHERWISE
*           THE LEADING  N BY K  PART OF THE ARRAY  B  MUST CONTAIN  THE
*           MATRIX B.
*           UNCHANGED ON EXIT.
*
*  LDB    - INTEGER.
*           ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSB = 'N' OR 'N' THEN
*           LDB MUST BE AT LEAST  MAX( 1, K ), OTHERWISE  LDB MUST BE AT
*           LEAST  MAX( 1, N ).
*           UNCHANGED ON EXIT.
*
*  BETA   - DOUBLE PRECISION.
*           ON ENTRY,  BETA  SPECIFIES THE SCALAR  BETA.  WHEN  BETA  IS
*           SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT.
*           UNCHANGED ON EXIT.
*
*  C      - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ).
*           BEFORE ENTRY, THE LEADING  M BY N  PART OF THE ARRAY  C MUST
*           CONTAIN THE MATRIX  C,  EXCEPT WHEN  BETA  IS ZERO, IN WHICH
*           CASE C NEED NOT BE SET ON ENTRY.
*           ON EXIT, THE ARRAY  C  IS OVERWRITTEN BY THE  M BY N  MATRIX
*           ( ALPHA*OP( A )*OP( B ) + BETA*C ).
*
*  LDC    - INTEGER.
*           ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   LDC  MUST  BE  AT  LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 3 BLAS ROUTINE.
*
*  -- WRITTEN ON 8-FEBRUARY-1989.
*     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
*     IAIN DUFF, AERE HARWELL.
*     JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD.
*     SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD.
*
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     .. LOCAL SCALARS ..
      LOGICAL            NOTA, NOTB
      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
      DOUBLE PRECISION   TEMP
*     .. PARAMETERS ..
      DOUBLE PRECISION   ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     SET  NOTA  AND  NOTB  AS  TRUE IF  A  AND  B  RESPECTIVELY ARE NOT
*     TRANSPOSED AND SET  NROWA, NCOLA AND  NROWB  AS THE NUMBER OF ROWS
*     AND  COLUMNS OF  A  AND THE  NUMBER OF  ROWS  OF  B  RESPECTIVELY.
*
      NOTA  = LSAME( TRANSA, 'N' )
      NOTB  = LSAME( TRANSB, 'N' )
      IF( NOTA )THEN
         NROWA = M
         NCOLA = K
      ELSE
         NROWA = K
         NCOLA = M
      END IF
      IF( NOTB )THEN
         NROWB = K
      ELSE
         NROWB = N
      END IF
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF(      ( .NOT.NOTA                 ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.NOTB                 ).AND.
     $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
         INFO = 2
      ELSE IF( M  .LT.0               )THEN
         INFO = 3
      ELSE IF( N  .LT.0               )THEN
         INFO = 4
      ELSE IF( K  .LT.0               )THEN
         INFO = 5
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 8
      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
         INFO = 10
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
         INFO = 13
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'DGEMM ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     AND IF  ALPHA.EQ.ZERO.
*
      IF( ALPHA.EQ.ZERO )THEN
         IF( BETA.EQ.ZERO )THEN
            DO 20, J = 1, N
               DO 10, I = 1, M
                  C( I, J ) = ZERO
   10          CONTINUE
   20       CONTINUE
         ELSE
            DO 40, J = 1, N
               DO 30, I = 1, M
                  C( I, J ) = BETA*C( I, J )
   30          CONTINUE
   40       CONTINUE
         END IF
         RETURN
      END IF
*
*     START THE OPERATIONS.
*
      IF( NOTB )THEN
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B + BETA*C.
*
            DO 90, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 50, I = 1, M
                     C( I, J ) = ZERO
   50             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 60, I = 1, M
                     C( I, J ) = BETA*C( I, J )
   60             CONTINUE
               END IF
               DO 80, L = 1, K
                  IF( B( L, J ).NE.ZERO )THEN
                     TEMP = ALPHA*B( L, J )
                     DO 70, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
   70                CONTINUE
                  END IF
   80          CONTINUE
   90       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B + BETA*C
*
            DO 120, J = 1, N
               DO 110, I = 1, M
                  TEMP = ZERO
                  DO 100, L = 1, K
                     TEMP = TEMP + A( L, I )*B( L, J )
  100             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  110          CONTINUE
  120       CONTINUE
         END IF
      ELSE
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B' + BETA*C
*
            DO 170, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 130, I = 1, M
                     C( I, J ) = ZERO
  130             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 140, I = 1, M
                     C( I, J ) = BETA*C( I, J )
  140             CONTINUE
               END IF
               DO 160, L = 1, K
                  IF( B( J, L ).NE.ZERO )THEN
                     TEMP = ALPHA*B( J, L )
                     DO 150, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
  150                CONTINUE
                  END IF
  160          CONTINUE
  170       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B' + BETA*C
*
            DO 200, J = 1, N
               DO 190, I = 1, M
                  TEMP = ZERO
                  DO 180, L = 1, K
                     TEMP = TEMP + A( L, I )*B( J, L )
  180             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  190          CONTINUE
  200       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF DGEMM .
*
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
C     Scalar replacement in the presence of conditional control flow
C
C     Steve Carr and Ken Kennedy
C
C     Software - Practive and Experience, Vol. 24, No. 1, pp. 51-77, 1994
C
C     Second example, which PIPS decides to be unprofitable

      subroutine scalarization31(a, b, c, d, e, m, n)
      real a(n), b(n), c(n), d(n), e(n), m(n)

      do i = 1,n
         if(m(i).lt.0.) a(i) = b(i) + c(i)
         d(i) = a(i) + e(i)
      enddo

      do i = 1,n
         if(m(i).lt.0.) a(i) = b(i) + c(i)
         d(i) = a(i) + e(i)
      enddo

      print *, d

      do i = 1,n
         if(m(i).lt.0.) a(i) = b(i) + c(i)
         d(i) = a(i) + e(i)
         d(i) = d(i) + a(i) * e(i)
      enddo

      end

ext: f


text:
!!
!! The compiler is stumbling over the type-is statement in the example below.
!! For some reason it does not see int32 (an intrinsic module named constant)
!! as a compile time constant (it is).
!!
!! If the subroutine is made an external procedure rather than a module
!! procedure, there is no compilation error.
!!
!! If the intrinsic module is used by the module instead of the procedure,
!! the compilation error also goes away.
!!
!! % ifort --version
!! ifort (IFORT) 13.1.2 20130514
!!
!! % ifort intel-bug-20130802a.f90
!! intel-bug-20130802a.f90(34): error #6683: A kind type parameter must be a compile-time constant.   [INT32]
!!     type is (integer(int32))
!! ---------------------^
!! compilation aborted for intel-bug-20130802a.f90 (code 1)
!!

module fubar_module

  ! put it here and everything is okay
  !use,intrinsic :: iso_fortran_env, only: int32

contains

  subroutine fubar (cs)
    ! put it here and the compiler stumbles
    use,intrinsic :: iso_fortran_env, only: int32
    class(*) :: cs
    select type (cs)
    type is (integer(int32))
    end select
  end subroutine

end module

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
use statistical

implicit none

integer,parameter :: n=8
integer :: i
real(8) :: a(n),mu,mu1,mu2

do i=1,n
   a(i) = i
enddo

call circ_mean_bootstrap_smallN(a,mu,mu1,mu2)

write(*,*) mu,mu1,mu2

end

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
!-----------------------------------------------------------------------
! Fortran interface to CUDA FFT Library
! written by Viktor K. Decyk, UCLA
      module gpupfft2_h
      implicit none
!
      interface
         subroutine gpufft2rrcuinit(nx,ny,ndim)
         implicit none
         integer :: nx, ny, ndim
         end subroutine
      end interface
!
      interface
         subroutine gpufft2cuinit(nx,ny,ndim)
         implicit none
         integer :: nx, ny, ndim
         end subroutine
      end interface
!
      interface
         subroutine gpufft2rrcudel()
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine gpufft2cudel()
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine gpupfft2rrcux(gp_f,gp_bsm,isign,indx,indy,kstrt,nvp,&
     &kxp1,kyp,nxh1d,kypd)
         implicit none
         integer :: isign, indx, indy, kstrt, nvp, kxp1, kyp
         integer :: nxh1d, kypd
         integer, dimension(2) :: gp_f, gp_bsm
         end subroutine
      end interface
!
      interface
         subroutine gpupfft2rrcuy(gp_g,gp_brm,isign,indx,indy,kstrt,nvp,&
     &kxp1,kyp,nyd)
         implicit none
         integer :: isign, indx, indy, kstrt, nvp, kxp1, kyp, nyd
         integer, dimension(2) :: gp_g, gp_brm
         end subroutine
      end interface
!
      interface
         subroutine gpupfft2rrcuxn(gp_fn,gp_bsm,isign,indx,indy,ndim,   &
     &kstrt,nvp,kxp1,kyp,nxh1d,kypd)
         implicit none
         integer :: isign, indx, indy, ndim, kstrt, nvp, kxp1, kyp
         integer :: nxh1d, kypd
         integer, dimension(2) :: gp_fn, gp_bsm
         end subroutine
      end interface
!
      interface
         subroutine gpupfft2rrcuyn(gp_gn,gp_brm,isign,indx,indy,ndim,   &
     &kstrt,nvp,kxp1,kyp,nyd)
         implicit none
         integer :: isign, indx, indy, ndim, kstrt, nvp, kxp1, kyp, nyd
         integer, dimension(2) :: gp_gn, gp_brm
         end subroutine
      end interface
!
      interface
         subroutine cgpuppsltpose(gp_f,gp_g,ani,nx,ny,kxp,kyp,kstrt,nxv,&
     &nyv)
         implicit none
         integer :: nx, ny, kxp, kyp, kstrt, nxv, nyv
         real :: ani
         integer, dimension(2) :: gp_f, gp_g
         end subroutine
      end interface
!
      interface
         subroutine cgpuppsltposen(gp_fn,gp_gn,ani,nx,ny,kxp,kyp,kstrt, &
     &ndim,nxv,nyv)
         implicit none
         integer :: nx, ny, kxp, kyp, kstrt, ndim, nxv, nyv
         real :: ani
         integer, dimension(2) :: gp_fn, gp_gn
         end subroutine
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


text:
C
C
      SUBROUTINE ANN1AS(IT)
C
C***********************************************************************
C
C     SETUP FOR TYPE -1 (Low temperature solid)
C
C     Inputs:
C
C     IT  =  pointer to ack array for this material
C
C     Parameters for pre 6/89 models and additional input are defined
C     before entry to this routine.
C
C************************************** 11/89 version slt **************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.D0)
      PARAMETER (ONE=1.D0)
      PARAMETER (TWO=2.D0)
      PARAMETER (THREE=3.D0)
      PARAMETER (FOUR=4.D0)
      PARAMETER (HALF=0.5D0)
      PARAMETER (P1OV10=-0.1D0)
      PARAMETER (P2OV3=0.66667D0)
      PARAMETER (PMATCH=5.D12)
C     number of intervals in Th table
      PARAMETER (NH=20)
      COMMON /FILEOS/ KLST, KINP
C     maximum number of aneos eos - must be 10 for csq
      PARAMETER (MAXMAT=10)
      COMMON /ANES/  ACK(99*MAXMAT),ZZS(30*MAXMAT),COT(30*MAXMAT)
     1 ,FNI(30*MAXMAT),RCT(MAXMAT+1),TCT(MAXMAT+1),RSOL(100*MAXMAT)
     2 ,RVAP(100*MAXMAT),TTWO(100*MAXMAT),SAVER(92),BOLTS,EIP(4370)
     3 ,LOCSV(MAXMAT+1),LOCKP(MAXMAT+1),LOCKPL(MAXMAT+1),NOTRAD
      SAVE   /FILEOS/,/ANES/
C     check user input
C
      IF(ACK(IT+33).LE.ZERO) THEN
        KTYPE=1
        IF(ACK(IT+44).EQ.ZERO) THEN
          CALL ANWARN(1)
          CALL ANMARK('Constant K1 cannot be zero')
          ACK(IT+44)=1.E-6
        END IF
      ELSE
        KTYPE=0
      END IF
      RR=ACK(IT+37)
      XX=ACK(IT+48)
      IF((ACK(IT+48).LE.ZERO).OR.(ACK(IT+48).GE.ACK(IT+11))) THEN
        ACK(IT+42)=P1OV10
        ACK(IT+48)=ACK(IT+11)/(ONE-ACK(IT+42))
      ELSE
        ACK(IT+42)=ONE-ACK(IT+11)/ACK(IT+48)
      END IF
      IF(ACK(IT+37).LE.ACK(IT+11)) THEN
        IF(KTYPE.EQ.0) THEN
          IF(ACK(IT+33).GT.ZERO) THEN
            R4R4=ONE+ACK(IT+21)/(TWO*ACK(IT+33)*PMATCH)
            R4R5=R4R4**2-ONE/ACK(IT+33)
            IF(R4R5.GT.ZERO) THEN
              ACK(IT+43)=R4R4-SQRT(R4R5)
            ELSE
              ACK(IT+43)=PMATCH/ACK(IT+21)
            END IF
          ELSE
            ACK(IT+43)=PMATCH/ACK(IT+21)
          END IF
        ELSEIF(ACK(IT+33).GT.ONE) THEN
          ACK(IT+43)=P2OV3/ACK(IT+33)
        ELSE
          ACK(IT+43)=P2OV3
        END IF
        ACK(IT+37)=ACK(IT+11)/(ONE-ACK(IT+43))
      ELSE
        ACK(IT+43)=ONE-ACK(IT+11)/ACK(IT+37)
      END IF
      IF((XX.NE.ZERO).AND.(XX.NE.ACK(IT+48))) THEN
        CALL ANWARN(0)
        WRITE(KLST,7) XX,ACK(IT+48)
    7   FORMAT(' Input IN(13) ignored',1PE13.6,' changed to',E13.6)
      END IF
      IF((RR.NE.ZERO).AND.(RR.NE.ACK(IT+37))) THEN
        CALL ANWARN(0)
        WRITE(KLST,6) RR,ACK(IT+37)
    6   FORMAT(' Input IN(14) ignored',1PE13.6,' changed to',E13.6)
      END IF
      IF(ACK(IT+18).LE.ZERO) ACK(IT+18)=123.456789
      IF(ACK(IT+15).LE.ZERO) ACK(IT+15)=ONE
      ACK(IT+34)=REAL(NH)
      ETA=ACK(IT+42)
        DO 10 I=1,NH+1
        IF(I.EQ.NH+1) ETA=ACK(IT+43)
        CALL ANN2AS(ETA,KTYPE,ACK(IT+33),ACK(IT+44),ACK(IT+45),
     &       ACK(IT+53),ACK(IT+54),ACK(IT+15),ACK(IT+I+74))
   10   ETA=ETA+(ACK(IT+43)-ACK(IT+42))/ACK(IT+34)
      XX=HALF*ACK(IT+32)**2/ACK(IT+4)
      IF(KTYPE.EQ.0) THEN
C       Us-Up low density form
        ETA=ACK(IT+42)
        RR=ACK(IT+11)/(ONE-ETA)
        TH=EXP(ACK(IT+15)*ETA)*(ACK(IT+12)
     &    +ACK(IT+75)*XX*(ETA/(ONE-ACK(IT+33)*ETA))**3)
        PC=ACK(IT+21)*ETA/(ONE-ACK(IT+33)*ETA)**2
        ACK(IT+51)=ETA*PC/(TWO*ACK(IT+11))+ACK(IT+3)-ACK(IT+4)*TH
        DPHDR=(ACK(IT+21)*(ONE+ACK(IT+33)*ETA)*(ONE-HALF*ACK(IT+15)*ETA)
     &  /(ONE-ACK(IT+33)*ETA)**3-ACK(IT+15)*(ACK(IT+5)*TH-HALF*PC))
     &  *ACK(IT+11)/RR**2
C        write(klst,3434) eta,rr,th,pc,dphdr
C3434     format(' eta,rho,th,ph,dphdr',/,1x,1p5e15.7)
        IF(DPHDR.LE.ZERO) THEN
          CALL ANWARN(1)
          CALL ANMARK('dPh/dRHO negative at lower density limit')
          WRITE(KLST,70) RR,PC,DPHDR
   70     FORMAT(' match point  rho=',1PE13.6,' Ph=',E13.6,/,
     &    ' dPh/dRHO=',E13.6,' To correct use larger density')
        END IF
        PC=PC-ACK(IT+5)*TH
        ACK(IT+49)=PC+DPHDR*RR
        ACK(IT+50)=-RR*RR*DPHDR
C       Us-Up high density form
        ETA=ACK(IT+43)
        RR=ACK(IT+11)/(ONE-ETA)
        TH=EXP(ACK(IT+15)*ETA)*(ACK(IT+12)
     &    +ACK(IT+75+NH)*XX*(ETA/(ONE-ACK(IT+33)*ETA))**3)
        PC=ACK(IT+21)*ETA/(ONE-ACK(IT+33)*ETA)**2
        ACK(IT+40)=ETA*PC/(TWO*ACK(IT+11))+ACK(IT+3)-ACK(IT+4)*TH
        DPHDR=(ACK(IT+21)*(ONE+ACK(IT+33)*ETA)*(ONE-HALF*ACK(IT+15)*ETA)
     &  /(ONE-ACK(IT+33)*ETA)**3-ACK(IT+15)*(ACK(IT+5)*TH-HALF*PC))
     &  *ACK(IT+11)/RR**2
C        write(klst,3434) eta,rr,th,pc,dphdr
        PC=PC-ACK(IT+5)*TH
        ACK(IT+38)=PC-HALF*DPHDR*RR
        ACK(IT+39)=DPHDR/(TWO*RR)
      ELSE
C       power series low density form
        ETA=ACK(IT+42)
        RR=ACK(IT+11)/(ONE-ETA)
        TH=EXP(ACK(IT+15)*ETA)*(ACK(IT+12)
     &    +ACK(IT+75)*XX*ACK(IT+44)*ETA**3/THREE)
        PC=ACK(IT+21)*ETA*(ONE+ACK(IT+44)*ETA+ACK(IT+45)*ETA**2
     &     +ACK(IT+53)*ETA**3+ACK(IT+54)*ETA**4)
        ACK(IT+51)=ETA*PC/(TWO*ACK(IT+11))+ACK(IT+3)-ACK(IT+4)*TH
        DPHDR=ACK(IT+21)*(ACK(IT+44)+TWO*ACK(IT+45)*ETA
     &     +THREE*ACK(IT+53)*ETA**2+FOUR*ACK(IT+54)*ETA**3)
        DPHDR=(PC/ETA-ACK(IT+15)*ACK(IT+5)*TH+(ONE-HALF*ACK(IT+15))
     &        *ETA*DPHDR)*ACK(IT+11)/RR**2
C        write(klst,3434) eta,rr,th,pc,dphdr
        IF(DPHDR.LE.ZERO) THEN
          CALL ANWARN(1)
          CALL ANMARK('dPh/dRHO negative at lower density limit')
          WRITE(KLST,70) RR,PC,DPHDR
        END IF
        PC=PC-ACK(IT+5)*TH
        ACK(IT+49)=PC+DPHDR*RR
        ACK(IT+50)=-RR*RR*DPHDR
C       power series high density form
        ETA=ACK(IT+43)
        RR=ACK(IT+11)/(ONE-ETA)
        TH=EXP(ACK(IT+15)*ETA)*(ACK(IT+12)
     &    +ACK(IT+75+NH)*XX*ACK(IT+44)*ETA**3/THREE)
        PC=ACK(IT+21)*ETA*(ONE+ACK(IT+44)*ETA+ACK(IT+45)*ETA**2
     &     +ACK(IT+53)*ETA**3+ACK(IT+54)*ETA**4)
        ACK(IT+40)=ETA*PC/(TWO*ACK(IT+11))+ACK(IT+3)-ACK(IT+4)*TH
        DPHDR=ACK(IT+21)*(ACK(IT+44)+TWO*ACK(IT+45)*ETA
     &     +THREE*ACK(IT+53)*ETA**2+FOUR*ACK(IT+54)*ETA**3)
        DPHDR=(PC/ETA-ACK(IT+15)*ACK(IT+5)*TH+(ONE-HALF*ACK(IT+15))
     &        *ETA*DPHDR)*ACK(IT+11)/RR**2
C        write(klst,3434) eta,rr,th,pc,dphdr
        IF(DPHDR.LE.ZERO) THEN
          CALL ANWARN(1)
          CALL ANMARK('dPh/dRHO negative at high density limit')
          WRITE(KLST,74) RR,PC,DPHDR
   74     FORMAT(' match point  rho=',1PE13.6,' Ph=',E13.6,/,
     &    ' dPh/dRHO=',E13.6,' To correct use smaller density')
        END IF
        PC=PC-ACK(IT+5)*TH
        ACK(IT+38)=PC-HALF*DPHDR*RR
        ACK(IT+39)=DPHDR/(TWO*RR)
      END IF
C
C     fix melt point if necessary to get it below low density join point
C
      PC=ACK(IT+18)
      ACK(IT+18)=MAX(ACK(IT+18),
     &    -(ACK(IT+49)+ACK(IT+50)/ACK(IT+48))/ACK(IT+5))
      IF(PC.LT.ACK(IT+18)) THEN
        CALL ANWARN(0)
        WRITE(KLST,68) PC,ACK(IT+18),ACK(IT+48)
   68   FORMAT(' Melt temperature increased from',1PE13.6,' to',E13.6,/,
     &  ' to yield RHOmelt < RHOlow = IN(13) =',E13.6)
      END IF
      ACK(IT+47)=-ACK(IT+50)/(ACK(IT+49)+ACK(IT+5)*ACK(IT+18))
      ACK(IT+23)=MIN(ACK(IT+23),ACK(IT+47))
      ACK(IT+96)=ACK(IT+5)/ACK(IT+47)
     &      +ACK(IT+4)*(LOG(ACK(IT+18)/ACK(IT+12))+ONE)
C
      RETURN
      END

max_stars_repo_name: chreinhardt/M-ANEOS

max_issues_repo_name: chreinhardt/M-ANEOS

max_forks_repo_name: chreinhardt/M-ANEOS


text:
subroutine rk4(SNPP)
    ! integrate in time using RK4 scheme
    use global, only : tt,Npts,iswitch,xyz,dt,Nx,Ny,Nz,uvwp,dt_file,t_amend,useKh,saveTSG

    implicit none
    real*8, dimension(3) :: x0,x1,uvw
    integer*8 :: t0,t1,ip,IPP
    integer*8, intent(in) :: SNPP

    t0=abs(iswitch-1)
    t1=iswitch
    do IPP=1,SNPP
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(x0,x1,uvw,ip) SHARED(IPP,SNPP,Npts,xyz,t_amend,t0,t1,dt,useKh)
        do ip=1,Npts
            x0=xyz(ip,:,IPP)
            x1=xyz(ip,:,IPP)

            call find_particle_uvw(0.0d0,ip,IPP,t0,t1,uvw)

            x1=x1+dt*uvw/6d0
            xyz(ip,:,IPP)=x0+dt*uvw/2d0

            call find_particle_uvw(t_amend,ip,IPP,t0,t1,uvw)
            x1=x1+dt*uvw/3d0
            xyz(ip,:,IPP)=x0+dt*uvw/2d0

            call find_particle_uvw(t_amend,ip,IPP,t0,t1,uvw)
            x1=x1+dt*uvw/3d0
            xyz(ip,:,IPP)=x0+dt*uvw

            call find_particle_uvw(t_amend*2d0,ip,IPP,t0,t1,uvw)
            xyz(ip,:,IPP)=x1+dt*uvw/6d0
        enddo
!$OMP END PARALLEL DO
    if (useKh) then
       call apply_mixing(IPP)
    endif
    call set_boundary(IPP)
    enddo

end subroutine rk4

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
      subroutine init_ug(AL,iwg_in,workr_n,kpt,iranm)
******************************************
cc     Written by Lin-Wang Wang, March 30, 2001.
*************************************************************************
**  copyright (c) 2003, The Regents of the University of California,
**  through Lawrence Berkeley National Laboratory (subject to receipt of any
**  required approvals from the U.S. Dept. of Energy).  All rights reserved.
*************************************************************************

******************************************

****************************************
cccccc this is just for one kpt.
****************************************
****  It stores the wavefunction in G space.
******************************************
      use fft_data
      use load_data
      use data
      implicit double precision (a-h,o-z)

      include 'mpif.h'
      include 'param.escan_real'

      integer status(MPI_STATUS_SIZE)
******************************************
       real*8 AL(3,3),ALt(3,3)
***********************************************
       complex*16 workr_n(mr_n)

*************************************************
****  input ug from file  and check the consistency
*************************************************

       if(iwg_in.eq.1) then
       call read_ug()
       endif
************************************************
*** generate the initial wavefunction ug from random if iwg_in.eq.0
************************************************
      if(iwg_in.eq.0) then
****  this is the one should be used later
	do m=1,mx
	do i=1,nr_n
        x1=ran1(iranm)
        x2=ran1(iranm)
        workr_n(i)=dcmplx(x1-0.5d0,x2-0.5d0)
cc        workr_n(i)=dcmplx(x1-0.5d0,0.d0)
	enddo
************************************************

        call mpi_barrier(MPI_COMM_K,ierr)

        call d3fft_comp(ug_n(1,m),workr_n,1,kpt)

      enddo
      endif
*************************************************
**** end generate the initial wavefunction from random
*************************************************
ccccccc it is not necessary to have the orthogonalization here
c        do m=1,mx
c        call orth_comp(ug_n(1,m),ug_n,m-1,1,kpt,
c     &   sug_n,swg,1,0,fnorm)     ! ipsp_all always be 1 here, since sug_n is not known yet
c	enddo
*************************************************
      return
      contains
*************************************************

       subroutine read_ug()

       complex*16,allocatable,dimension(:) :: ugtemp

        call mpi_barrier(MPI_COMM_K,ierr)
cccccc  file (11) should have been opened before the do loop for kpt=1,nkpt

ccccccccccccc check the header.
        if(inode==nnodes.and.kpt.eq.1) then

         read(11)n1t,n2t,n3t,mxt
         read(11)Ecutt
         read(11)ALt
         read(11)nnodes_o

         if(n1t.ne.n1.or.n2t.ne.n2.or.n3t.ne.n3) then
         write(6,*) "n1t,n2t,n3t changed, stop", n1t,n2t,n3t
         call mpi_abort(MPI_COMM_WORLD,ierr)
         endif

         if(Ecut.ne.Ecutt.or.mxt.ne.mx) then
         write(6,*) "Ecutt,mxt changed, stop", Ecutt,mxt
         call mpi_abort(MPI_COMM_WORLD,ierr)
         endif

         diff=0.d0
         do i=1,3
            do j=1,3
               diff=diff+dabs(AL(i,j)-ALt(i,j))
             enddo
         enddo

         if(diff.gt.0.001d0) then
            write(6,*) "the AL.ne.ALt, stop"
         call mpi_abort(MPI_COMM_WORLD,ierr)
         endif

         if(nnodes_o.ne.nnodes) then
           write(6,*) "nnodes_o changed, stop", nnodes_o,nnodes
         call mpi_abort(MPI_COMM_WORLD,ierr)
         endif

      endif     ! check the header, make sure it is the same
ccccccccccccccccccccccccccccccccccccccccc
         call mpi_barrier(MPI_COMM_K,ierr)

         if(inode==nnodes) then
            allocate(ugtemp(mg_nx))

            do i=0,nnodes-2
               do iwavefun=1,mx
                  read(11)ugtemp    ! assuming the mg_nx is the same as before
                  ug_n(:,iwavefun)=ugtemp
               end do

          call mpi_send(ug_n,mg_nx*mx,MPI_DOUBLE_COMPLEX,i,
     &            102,MPI_COMM_K,ierr)
            end do
c     Now read in the data for node nnodes
            do iwavefun=1,mx
               read(11)ugtemp
               ug_n(:,iwavefun)=ugtemp
            end do
            deallocate(ugtemp)

         else   ! for other nodes

         call mpi_recv(ug_n,mg_nx*mx,MPI_DOUBLE_COMPLEX,nnodes-1,
     &       102,MPI_COMM_K,status,ierr)

         end if

        return
        end subroutine read_ug

        end

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
        subroutine nl_to_RGB(rad,glwref,contrast,cntref,iverbose,rc,gc,bc)

        include 'wa.inc'
        include 'wac.inc'

        real luminance ! candela / m**2
        real luma_of_counts
        integer init /0/
        save init,fasun

        parameter (day_int = day_int0)

!       Exponential clipping of gamma function
!       expgamma(x) = 1. - exp(-x**(1./gamma))
!       expgaminv(x) = (-(log(1.-x)))**gamma ! inverse of expgamma
!       radhi(glwref,cntref) = (10.**glwref) / (expgaminv(cntref/255.))
!       rad_to_counts2(rad1) = expgamma(rad1/radhi(glwref,cntref)) * 255.

!       Straight gamma function
        fgamma(x) = x**(1./gamma)
        fgaminv(x) = x**gamma ! inverse of expgamma
        radhi(glwref,cntref) = (10.**glwref) / (fgaminv(cntref/255.))
        rad_to_counts2(rad1) = fgamma(rad1/radhi(glwref,cntref)) * 255.

!       Convert nl to rintensity

!       Radiance units are watts/(m**2 sr nm)

!       Day_int0 is luminance of sunlight from 1AU distance uniformly
!       scattered throughout a 360 degree sphere. This will be defined as
!       the spectral radiance matching that of the solar spectrum.
!       Note solar constant is 1361.5 watts/m**2
        real rad(nc)
        real rel_solar(nc), cdm2(nc)
        real fasun(nct), rel_solar_extrap(nct), farad(nct)
        real nl_int

        if(init .eq. 0)then
            write(6,*)' initializing a call to get_fluxsun ',init
            call get_fluxsun(wa_tri,nct,1,fasun)
            init = 1
        endif

        if(maxval(rad(:)) .eq. 0.)then
            rc = 0.; gc = 0.; bc = 0.
            return
        endif

!       nl_int = rgb_to_y(rad(1),rad(2),rad(3))

        rel_solar(:) = rad(:) / day_int
        rel_solar_luminance = rgb_to_y(rel_solar(1),rel_solar(2) &
                                      ,rel_solar(3))
        if(iverbose .eq. 1)then
          write(6,*)'rel_solar_luminance = ',rel_solar_luminance
        endif

!       1 lambert = 3183.0988 candela/m**2        (luminance)
!       1 lux = 1 lumen/m**2                      (illuminance)
!       1 lux = 1 lumen/m**2                      (luminous exitance)
!       1 candela = 1/683.002 watts/sr (at 550nm) (luminous intensity)
!       1 watt = 683 lumens (555nm source)        (radiant flux)
!       1 watt =  93 lumens (solar spectrum)      (radiant flux)
!       1 watt/m**2                               (irradiance)
!       1 lumen = 1 cd*sr                         (luminous flux)
!       1 candela = 1 lumen per steradian         (luminous intensity)
        cdm2(:) = (rad(:)/1e9) * 3183.0988 ! nl to candela/m**2

!       Extrapolate color
        call extrap_color(rel_solar,rel_solar_extrap) ! dimensionless
        farad(:) = fasun(:) * rel_solar_extrap(:)     ! spect irradiance
        if(iverbose .eq. 1)write(6,*)' farad = ',farad

!       Convert sprad to xyz
        call get_tricolor(farad,iverbose,xx,yy,zz,x,y,z,luminance)
        if(iverbose .eq. 1)write(6,*)'xyzfarad = ',x,y,z
        if(iverbose .eq. 1)write(6,*)'luminance farad = ',luminance

!       Convert xyz to (linear) rgb
        ct = 5780.
        if(iverbose .eq. 1)write(6,*)'ct = ',ct
        call xyztosrgb(x,y,z,ct,r,g,b)

!       Convert rintensity to rgb (gamma correction)
        call linearrgb_to_counts(r,g,b,rc,gc,bc)

        if(iverbose .eq. 1)write(6,*)'glwref/contrast = ',glwref,contrast

!       if(iverbose .eq. 1)write(6,*)'desired rad = ',desired_rad

        luma_of_counts = RGB2luma(rc,gc,bc)
        if(iverbose .eq. 1)write(6,*)'luma_of_counts = ',luma_of_counts

!       solar_counts = 240.
        solar_counts = rad_to_counts2(day_int)
        if(iverbose .eq. 1)then
            write(6,*)'day_int = ',day_int
            write(6,*)'glwref = ',glwref
            write(6,*)'cntref = ',cntref
            write(6,*)'solar_counts = ',solar_counts
        endif

!       desired_luma = 240. * rel_solar_luminance**(1./gamma)
!       desired_luma = solar_counts + log10(rel_solar_luminance) * contrast
        desired_luma = rad_to_counts2(day_int*rel_solar_luminance)
        if(iverbose .eq. 1)write(6,*)'desired luma = ',desired_luma
        scale_luma = desired_luma / luma_of_counts
        if(iverbose .eq. 1)write(6,*)'scale_luma = ',scale_luma
        rc = rc * scale_luma
        gc = gc * scale_luma
        bc = bc * scale_luma

        return
        end

        subroutine extrap_color(a_nc,a_nct)

        include 'wa.inc'
        include 'wac.inc'

        real a_nc(nc), a_nct(nct)

!       Assuming nc = 3, find polynomial fit
        X1 = wa(1);   X2 = wa(2);   X3 = wa(3)
        Y1 = a_nc(1); Y2 = a_nc(2); Y3 = a_nc(3)

        a = ((Y2-Y1)*(X1-X3) + (Y3-Y1)*(X2-X1)) / &
            ((X1-X3)*(X2**2-X1**2) + (X2-X1)*(X3**2-X1**2))
        b = ((Y2 - Y1) - A*(X2**2 - X1**2)) / (X2-X1)
        c = Y1 - A*X1**2 - B*X1

        a_nct(:) = a*wa_tri(:)**2 + b*wa_tri(:) + c

        return
        end

        subroutine rad_to_xyz(rad_in,x,y,z)

        real rad_in(3),rad(3)

        rad(1) = rad_in(1) + 0.2 * rad_in(3)
        rad(2) = rad_in(2)
        rad(3) = rad_in(3)

        x = rad(1) / sum(rad)
        y = rad(2) / sum(rad)
        z = rad(3) / sum(rad)

        return
        end

        subroutine linearrgb_to_counts(r,g,b,rc,gc,bc)

        parameter (gamma = 2.2)

        rr = r ** (1./gamma)
        gg = g ** (1./gamma)
        bb = b ** (1./gamma)

        rc = rr*255.
        gc = gg*255.
        bc = bb*255.

        return
        end

        subroutine xyztorgb(x,y,z,r,g,b)

        xx = x; yy = y; zz = z

        R=  (0.41847   *XX)-(0.15866  *YY)-(0.082835*ZZ)
        G= -(0.091169  *XX)+(0.25243  *YY)+(0.015708*ZZ)
        B=  (0.00092090*XX)-(0.0025498*YY)+(0.17860 *ZZ)

        return
        end

        subroutine rgbtoxyz(r,g,b,x,y,z)

        x = (.49    * r + .31    * g + .20    * b) / 0.17697
        y = (.17697 * r + .81240 * g + .01063 * b) / 0.17697
        z = (             .01    * g + .99    * b) / 0.17697

        return
        end

        subroutine xyztosrgb(x,y,z,ct,r,g,b)

        colfrac = (ct - 5000.) / (6500. - 5000.)

        r65 =  3.2406*x -1.5372*y -0.4986*z
        g65 = -0.9689*x +1.8758*y +0.0415*z
        b65 =  0.0557*x -0.2040*y +1.0570*z

        r50 =  3.1339*x -1.6169*y -0.4906*z
        g50 = -0.9785*x +1.9160*y +0.0333*z
        b50 =  0.0720*x -0.2290*y +1.4057*z

        r = r65 * colfrac + r50 * (1.-colfrac)
        g = g65 * colfrac + g50 * (1.-colfrac)
        b = b65 * colfrac + b50 * (1.-colfrac)

!       Clip at 0. when outside the gamut
        r = max(r,0.)
        g = max(g,0.)
        b = max(b,0.)

        return
        end

        subroutine srgbtoxyz(r,g,b,ct,x,y,z)

        colfrac = (ct - 5000.) / (6500. - 5000.)

        x65 =  0.4124*r + 0.3576*g + 0.1805*b
        y65 =  0.2126*r + 0.7152*g + 0.0722*b
        z65 =  0.0193*r + 0.1192*g + 0.9505*b

        x50 =  0.4361*r + 0.3851*g + 0.1430*b
        y50 =  0.2225*r + 0.7169*g + 0.0606*b
        z50 =  0.0139*r + 0.0971*g + 0.7141*b

        x = x65 * colfrac + x50 * (1.-colfrac)
        y = y65 * colfrac + y50 * (1.-colfrac)
        z = z65 * colfrac + z50 * (1.-colfrac)

        return
        end

        subroutine get_fluxsun(wa,nc,iverbose,fa)

!       include 'wa.inc'
!       include 'wac.inc'

!       Returns solar spectral irradiance for selected wavelengths

        real fa(nc) ! watts/(m**2 nm)
        real wa(nc) ! wavelengths in microns

        if(iverbose .eq. 1)write(6,*)'  subroutine get_fluxsun',nc,wa

        TS = 5900. ! 5780.
        DS = 1.0 * 1.496e13 ! au to cm
        RS = 1.0 * 6.96e10  ! solar radii to cm
        xc = 0.; yc = 0.; zc = 0.

        do ic = 1,nc
          W = wa(ic) * 1e-4 ! microns to cm
          w_ang = wa(ic) * 10000.
          BB = (.0000374/(W**5.)) / (EXP(1.43/(W*TS))-1.)
          if(iverbose .eq. 1)write(6,*)'ic/w/bb',ic,w,bb
          FA(IC)=((RS/DS)**2)*BB*1E-8 ! erg/cm2/s/A
          FA(IC)=FA(IC) * .01         ! convert to W/m2/nm
        enddo ! ic

        if(iverbose .eq. 1)write(6,*)'fa sun is ',fa

        return
        end

        subroutine get_tricolor(fa,iverbose,xc,yc,zc,x,y,z,luminance)

        include 'wa.inc'
        include 'wac.inc'

!       Integrate spectral irradiance with the color matching functions

        real fa(nct)   ! spectral irradiance array vs wavelength (Input)
        real luminance ! candela / m**2                    (Output = yc)
        character*255 static_dir

!       3500 to 8000 Angstroms in 10 steps (color matching functions)
        real x1(nct)!/0, .014, .336, .005, .433, 1.062, .283, .011, 0, 0/
        real y1(nct)!/0, .0004, .038, .323, .995, .631, .107, .004, .0001, 0/
        real z1(nct)!/0, .068, 1.773, .272, .009, 0,0,0,0,0/
        integer init/0/
        save x1,y1,z1,init

        if(init .eq. 0)then
          call get_directory('static',static_dir,len_dir)
          open(11,file=trim(static_dir)//'/cie2.txt',status='old')
!         open(11,file='cie2.txt',status='old')
          do ict = 1,nct
            read(11,*)inm,x1(ict),y1(ict),z1(ict)
          enddo ! ict
          init = 1
        endif

        xc = 0.; yc = 0.; zc = 0.

!       Integrate with trapezoidal rule
        do ic = 2,nct
          icm = ic-1
          W = wa_tri(ic) * 1e-4 ! microns to cm
          w_ang = wa_tri(ic) * 10000.

          xs = 0.5 * (x1(icm)*fa(icm) + x1(ic)*fa(ic))
          ys = 0.5 * (y1(icm)*fa(icm) + y1(ic)*fa(ic))
          zs = 0.5 * (z1(icm)*fa(icm) + z1(ic)*fa(ic))

          xc = xc + xs
          yc = yc + ys
          zc = zc + zs

          if(iverbose .eq. 1)then
            write(6,1)ic,w_ang,fa(ic),xs,ys,zs,xc,yc,zc
1           format(' ic,wa,fa',i4,f8.0,f10.5,' xyzs',3f10.5,' xyzc',3f10.5)
          endif

        enddo ! ic

        x = xc / (xc+yc+zc)
        y = yc / (xc+yc+zc)
        z = zc / (xc+yc+zc)

        luminance = yc

        if(iverbose .eq. 1)write(6,11)x,y,z,luminance
11      format(' xyzl tricolor is ',3f10.6,f12.0)

        return
        end

      subroutine hsl_to_rgb(hue,sat,rintens,red,grn,blu)

      real luma

!     Hue is 0:R, 1:B, 2:G, 3:R

      red1 = max(1.0 - abs(hue - 0.0),0.0)
      red2 = max(1.0 - abs(hue - 3.0),0.0)
      red = max(red1,red2)
      grn = max(1.0 - abs(hue  - 2.0),0.0)
      blu = max(1.0 - abs(hue  - 1.0),0.0)

!     Normalize to the max intensity
      colmax = max(red,grn,blu)
      if(colmax .gt. 0.)then
          red = red/colmax
          grn = grn/colmax
          blu = blu/colmax
      endif

      red = (red*sat) + 1.0*(1.0-sat)
      grn = (grn*sat) + 1.0*(1.0-sat)
      blu = (blu*sat) + 1.0*(1.0-sat)

      red = red * rintens
      grn = grn * rintens
      blu = blu * rintens

!     Scale the RGB intensity according to the luma calculation

      luma = .30 * red + .59 * grn + .11 * blu

      if(luma .gt. 0.)then
          ratio_corr = rintens / luma
      else
          ratio_corr = 1.0
      endif

      red = red * ratio_corr
      grn = grn * ratio_corr
      blu = blu * ratio_corr

      return
      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      PROGRAM runEEMD

      IMPLICIT NONE
      INCLUDE 'eemd.common'
      REAL,    PARAMETER :: STDEV=0.2
      INTEGER, PARAMETER :: EMEMBERS=100
      INTEGER, PARAMETER :: IMFS=12

      INTEGER :: i,j,TimeSeriesLength
      REAL, DIMENSION(MAXSIZE) :: indata
      INTEGER:: idum
      REAL :: eemd_result(MAXSIZE,IMFS+2)

      idum=3571
      open (4, file='tempTS.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=1)
      OPEN (15,FILE='tempEEMD.dat',FORM='UNFORMATTED',ACCESS='DIRECT', RECL=1)
!      OPEN(UNIT = 15, FILE = 'newfort.txt', STATUS = 'UNKNOWN')
      open (19, file='tempEEMD.ctl')

      do i = 1, MAXSIZE
         read(4, rec=i, ERR=10) indata(i)
      enddo
  10  continue
      TimeSeriesLength = i-1

      CALL eemd(TimeSeriesLength,indata,STDEV,EMEMBERS,IMFS,idum,eemd_result)
      DO i=1,TimeSeriesLength
        DO j=1,IMFS+2
          WRITE(15,REC=i+(j-1)*TimeSeriesLength) (eemd_result(i,j))
        ENDDO
      ENDDO

 500  FORMAT(14F12.8)

      WRITE (19,*) 'DSET      ^tempEEMD.dat'
      WRITE (19,*) 'TITLE     Sample Model Data'
      WRITE (19,*) 'UNDEF    -9.99E+30 '
      WRITE (19,*) 'XDEF     1 linear   0 1'
      WRITE (19,*) 'YDEF     1 linear   0 1'
      WRITE (19,*) 'EDEF       16	NAMES	TS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15'
      WRITE (19,*) 'ZDEF       1  linear 1 1'
      WRITE (19,*) 'TDEF       ',TimeSeriesLength,' LINEAR 1JAN2000 1MO'
      WRITE (19,*) 'VARS      1 '
      WRITE (19,*) '   imf  0  99  100 Data Points'
      WRITE (19,*) 'ENDVARS'
      WRITE (19,*) ' '

      END PROGRAM runEEMD

max_stars_repo_name: mathnathan/EEMD

max_issues_repo_name: mathnathan/EEMD

max_forks_repo_name: mathnathan/EEMD


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program dart_to_clm

!----------------------------------------------------------------------
! purpose: interface between DART and the CLM model
!
! method: Read DART state vector and overwrite values in a CLM restart file.
!
!----------------------------------------------------------------------

use        types_mod, only : r8
use    utilities_mod, only : initialize_utilities, finalize_utilities, &
                             find_namelist_in_file, check_namelist_read, &
                             logfileunit, open_file, close_file, &
                             error_handler, E_MSG
use time_manager_mod, only : time_type, print_time, print_date, operator(-), get_time
use        model_mod, only : static_init_model, fill_missing_r8_with_orig, &
                             get_model_size, get_clm_restart_filename, &
                             read_model_time

implicit none

! version controlled file description for error handling, do not edit
character(len=*), parameter :: source   = &
   "$URL$"
character(len=*), parameter :: revision = "$Revision$"
character(len=*), parameter :: revdate  = "$Date$"

!------------------------------------------------------------------
! The namelist variables
!------------------------------------------------------------------

character(len=256) :: dart_to_clm_input_file = 'model_restart.nc'
logical            :: advance_time_present   = .false.

namelist /dart_to_clm_nml/ dart_to_clm_input_file, &
                           advance_time_present

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

character(len=256) :: clm_restart_filename
integer            :: iunit, io
type(time_type)    :: model_time, adv_to_time

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

call initialize_utilities(progname='dart_to_clm')

!----------------------------------------------------------------------
! Call model_mod:static_init_model() which reads the clm namelists
! to set grid sizes, etc.
!----------------------------------------------------------------------

call static_init_model()

! Read the namelist to get the input filename.

call find_namelist_in_file("input.nml", "dart_to_clm_nml", iunit)
read(iunit, nml = dart_to_clm_nml, iostat = io)
call check_namelist_read(iunit, io, "dart_to_clm_nml")

call get_clm_restart_filename( clm_restart_filename )

write(*,*)
write(*,'(''dart_to_clm:converting DART file '',A, &
      &'' to clm restart file '',A)') &
     trim(dart_to_clm_input_file), trim(clm_restart_filename)

!----------------------------------------------------------------------
! Reads the valid time, the state, and the target time.
!----------------------------------------------------------------------

model_time = read_model_time(clm_restart_filename)

!----------------------------------------------------------------------
! update the current clm state vector
!----------------------------------------------------------------------

call fill_missing_r8_with_orig(clm_restart_filename, dart_to_clm_input_file, model_time)

!----------------------------------------------------------------------
! Log what we think we're doing, and exit.
!----------------------------------------------------------------------

call print_date( model_time,'dart_to_clm:clm  model date')
call print_time( model_time,'dart_to_clm:DART model time')
call print_date( model_time,'dart_to_clm:clm  model date',logfileunit)
call print_time( model_time,'dart_to_clm:DART model time',logfileunit)

if ( advance_time_present ) then
   call error_handler(E_MSG,'dart_to_clm','warning: DART not configured to advance CLM', &
            source, revision, revdate)
   call print_time(adv_to_time,'dart_to_clm:advance_to time')
   call print_date(adv_to_time,'dart_to_clm:advance_to date')
   call print_time(adv_to_time,'dart_to_clm:advance_to time',logfileunit)
   call print_date(adv_to_time,'dart_to_clm:advance_to date',logfileunit)
endif

call finalize_utilities('dart_to_clm')

end program dart_to_clm

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
SUBROUTINE beam_mm(mm,fs,ell)
!
! This subroutine forms the consistent mass matrix of a beam element.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::fs,ell
 REAL(iwp),INTENT(OUT)::mm(:,:)
 REAL(iwp)::fac
 fac=(fs*ell)/420.0_iwp
 mm(1,1)=156.0_iwp*fac
 mm(3,3)=mm(1,1)
 mm(1,2)=22.0_iwp*ell*fac
 mm(2,1)=mm(1,2)
 mm(3,4)=-mm(1,2)
 mm(4,3)=mm(3,4)
 mm(1,3)=54.0_iwp*fac
 mm(3,1)=mm(1,3)
 mm(1,4)=-13.0_iwp*ell*fac
 mm(4,1)=mm(1,4)
 mm(2,3)=-mm(1,4)
 mm(3,2)=mm(2,3)
 mm(2,2)=4.0_iwp*(ell**2)*fac
 mm(4,4)=mm(2,2)
 mm(2,4)=-3.0_iwp*(ell**2)*fac
 mm(4,2)=mm(2,4)
RETURN
END SUBROUTINE beam_mm

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
== Mathematical constants
PI  << 3.141 592 653 59
E   << 2.718 281 828 46
PHI << 1.618 033 988 75
SQ2 << 1.414 213 562 37

π << PI
φ << PHI

== Physical constants
PLANCK        CONSTANT << 6.626 070 04e-34
GRAVITATIONAL CONSTANT << 6.674 08e-10


submodule (pathlib) pathlib_intel

implicit none (type, external)

contains

module procedure directory_exists

inquire(directory=expanduser(path), exist=exists)

end procedure directory_exists

end submodule pathlib_intel


text:
        SUBROUTINE SPCOPN(ROOT,ISTAT)
*
*  Module number:
*
*  Module name: SPCOPN
*
*  Keyphrase:
*  ----------
*       open input data files for polarimetry reduction
*
*  Description:
*  ------------
*       This routine opens the input data files (except the .c1h file
*	which is already open).
*
*  FORTRAN name: spcopn.for
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       ymsput,yfname
*  SDAS:
*       uimopn, uhdgsr, uimgid
*
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1         Oct 92      D. Bazell       Modified version of YCLOPN
*     2		Sep 93	    H. Bushouse	    Modified to read all input files
*     2.1       Feb 98      M. De La Pena   Added KYDPLY to CONFG1.
*-------------------------------------------------------------------------------
* inputs:
*       ROOT - root name of the file
*
* outputs:
*       ISTAT - error status
*
        CHARACTER*64 ROOT
        INTEGER ISTAT
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
      INTEGER RDONLY
      PARAMETER (RDONLY = 1)
C
C Common block containing confiquration parameters
C
        CHARACTER*5 DET
        CHARACTER*3 FGWAID,APERID,YTYPE(3)
        CHARACTER*1 POLID
        INTEGER FCHNL,NCHNLS,OVERSN,NXSTEP,YBASE,YRANGE,YSTEPS,
     *          INTS,SLICES,NPAT,NREAD,NCLEAR,LVTIME
        LOGICAL HEADER,TRAILR,DEFDDT,KYDPLY
        COMMON /CONFG1/KYDPLY,DET,FGWAID,APERID,POLID,YTYPE
        COMMON /CONFG2/FCHNL,NCHNLS,OVERSN,NXSTEP,YBASE,YRANGE,YSTEPS,
     *          INTS,SLICES,NPAT,NREAD,NCLEAR,LVTIME,HEADER,TRAILR,
     *          DEFDDT
C
C Common block containing input/output file descriptions
C
	INTEGER IDS(30),NAXIS(30),NAXIS1(30),NAXIS2(30),GCOUNT(30)
	REAL FILL(30)
	COMMON /IOCOM/IDS,GCOUNT,NAXIS,NAXIS1,NAXIS2,FILL
C
C Common block containing input file name
C
        CHARACTER*64 INFILE, INEXT
        COMMON /CINFILE/INFILE, INEXT
C
C Local variables
C
        CHARACTER*3 QUAL(4)
        INTEGER DIMEN(8),DTYPE,I
        CHARACTER*64 NAME
        CHARACTER*80 CONTXT
	DATA QUAL/'c0h','c1h','c2h','cqh'/
C---------------------------------------------------------------------------
C
C loop on input files (skip c1h, it's already open)
C
        DO 100 I=1,4
	   IF (I .EQ. 2) GO TO 100
C
C Open the file
C
           CALL YFNAME(ROOT,QUAL(I),1,0,NAME)
c          CALL UIMOPN(NAME,RDONLY,INID(I),ISTAT)
           CALL UIMOPN(NAME,RDONLY,IDS(I+10),ISTAT)
           IF(ISTAT.NE.0)THEN
              CONTXT='ERROR opening data file '//NAME
              GO TO 999
           ENDIF
C
C get dimensions of the file
C
c          CALL UIMGID(INID(I),DTYPE,INNAX(I),DIMEN,ISTAT)
           CALL UIMGID(IDS(I+10),DTYPE,NAXIS(I+10),DIMEN,ISTAT)
           IF(ISTAT.NE.0)THEN
              CONTXT='ERROR reading file '//NAME
              GO TO 999
           ENDIF
c          IF(INNAX(I).GT.2)THEN
           IF(NAXIS(I+10).GT.2)THEN
              CONTXT='INVALID naxis for file '//NAME
              GOTO 999
           ENDIF
c          INNAX1(I)=DIMEN(1)
c          INNAX2(I)=DIMEN(2)
           NAXIS1(I+10)=DIMEN(1)
           NAXIS2(I+10)=DIMEN(2)
C
C get number of groups in the file
C
c          CALL UHDGSI(INID(I),'GCOUNT',INGCNT(I),ISTAT)
           CALL UHDGSI(IDS(I+10),'GCOUNT',GCOUNT(I+10),ISTAT)
           IF(ISTAT.NE.0)THEN
              CONTXT='ERROR getting GCOUNT from file '//NAME
              GO TO 999
           ENDIF
 100    CONTINUE
        ISTAT=0
        GO TO 1000
999     CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine allocate_damage              *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 12/5/2020  rhd             *
c     *                                                              *
c     *     allocates information for the damage routines as needed  *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine allocate_damage ( dowhat )
c
      use global_data, only: out, mxedof, nelblk, nonode, noelem,
     &                       mxconn, props, incflg
      use elem_extinct_data
      use node_release_data
      use damage_data
      use dam_param_code, only : dam_param
c
      implicit none
c
      integer :: dowhat
c
      integer :: dum, elem, elem_ptr, i
      double precision ::
     &      dumd1, dumd2, dumd3, dumd4, dumd5, dumd6, dumd7,
     &     dumd8, dumd9, dumd10, porosity, plast_strain,
     &     values(20)
      double precision, parameter :: zero = 0.d0
      logical :: debug, duml
      real :: dumr
c
      debug = .true.
c
c                               allocate dam_state, dam_ifv,
c                               dam_blk_killed
c
      select case( dowhat )
c
      case( 1 ) ! dowhat
c
         if( allocated(dam_state) ) then
            call errmsg(215,dum,'dam_state',dumr,dumd1)
            return
         end if
         allocate( dam_state(num_kill_elem) )
         dam_state(1:num_kill_elem) = 0
c
         if( allocated(dam_ifv) ) then
            call errmsg(215,dum,'dam_ifv',dumr,dumd1)
            return
         end if
         allocate( dam_ifv(mxedof,num_kill_elem) )
         dam_ifv(1:mxedof,1:num_kill_elem) = zero
c
         if( allocated(dam_blk_killed) ) then
            call errmsg(215,dum,'dam_blk_killed',dumr,dumd1)
            return
         end if
         allocate( dam_blk_killed(nelblk) )
         dam_blk_killed(1:nelblk) = .false.
c
c                               allocate dam_print_list
c                               also allocate old_mises and
c                               old_mean to calculate change in
c                               stress over step.
c
      case( 2 ) ! dowhat
         if( allocated(dam_print_list) ) deallocate( dam_print_list )
         allocate( dam_print_list(num_print_list) )
c
         if( allocated(old_mises) ) deallocate( old_mises )
         allocate( old_mises( num_print_list ) )
         old_mises(1:num_print_list) = zero
c
         if( allocated( old_mean ) ) deallocate( old_mean )
         allocate( old_mean( num_print_list ) )
         old_mean(1:num_print_list) = zero
c
c
c                                allocate kill_order_list
c
      case( 3 ) ! dowhat
         if( allocated(kill_order_list) ) deallocate(kill_order_list)
         allocate(kill_order_list(num_kill_order_list))
c
c                               allocate dam_node_elecnt
c
      case( 4 ) ! dowhat
         if( allocated(dam_node_elecnt) ) then
            call errmsg(215,dum,'dam_node_elecnt',dumr,dumd1)
            return
         end if
         allocate(dam_node_elecnt(nonode))
         dam_node_elecnt(1:nonode) = 0
c
c
c                               allocate dam_face_nodes and
c                               dam_dbar_elems. these store the
c                               nodes for each element to compute
c                               average extension normal to crack plane,
c                               the average extension when the
c                               element reaches critical porosity,
c                               and the average extension on the
c                               previous load step.
c
      case( 5 ) ! dowhat
         if( allocated(dam_face_nodes) ) then
            call errmsg(215,dum,'dam_face_nodes_elecnt',dumr,dumd1)
            return
         end if
         if( allocated(dam_dbar_elems) ) then
            call errmsg(215,dum,'dam_dbar_elems',dumr,dumd1)
            return
         end if
         allocate( dam_face_nodes(4,num_kill_elem) )
         allocate( dam_dbar_elems(2,num_kill_elem) )
c
c                                allocate: (for node_release crack growth)
c                                  crack_plane_nodes -- nodes on crack plane
c                                  inv_crkpln_nodes -- inverse of
c                                      crack_plane_nodes
c                                  num_neighbors -- # of neighboring nodes
c                                      for each crack plane node
c                                  neighbor_nodes -- the neighboring nodes
c                                  crack_front_nodes -- linked list of
c                                      the nodes on the crack front
c                                  crkpln_nodes_react -- an array to hold
c                                      the reaction force at the released
c                                      node at release time
c                                  crkpln_nodes_state -- the current release
c                                      state for the killed nodes
c
      case( 6 ) ! dowhat
         if( allocated(crack_plane_nodes) )
     &       deallocate( crack_plane_nodes )
         allocate( crack_plane_nodes(num_crack_plane_nodes) )
         if( allocated(inv_crkpln_nodes) )
     &       deallocate( inv_crkpln_nodes )
         allocate( inv_crkpln_nodes(nonode) )
c
         if( allocated(num_neighbors) ) deallocate( num_neighbors )
         allocate( num_neighbors(num_crack_plane_nodes) )
c
         if( allocated(neighbor_nodes) ) deallocate(neighbor_nodes)
         allocate( neighbor_nodes(mxconn,num_crack_plane_nodes) )
c
         if( allocated(crack_front_nodes) )
     &       deallocate( crack_front_nodes )
         allocate( crack_front_nodes(num_crack_plane_nodes,2) )
c
         if( allocated(crkpln_nodes_react) )
     &         deallocate( crkpln_nodes_react )
         allocate( crkpln_nodes_react(num_crack_plane_nodes) )
         if( allocated(crkpln_nodes_state) )
     &         deallocate( crkpln_nodes_state )
         allocate( crkpln_nodes_state(num_crack_plane_nodes) )
c
c                            allocate: (for node_release crack growth
c                                       and traction separation law)
c                                node_release_frac(crack_plane_nodes) --
c                                     stores the last fraction of release
c                                     for a released node under traction
c                                     separation.  We initialize it to
c                                     zero in this routine
c
      case( 7 ) ! dowhat
         if( allocated(node_release_frac) )
     &       deallocate( node_release_frac )
         allocate( node_release_frac(num_crack_plane_nodes) )
         node_release_frac(1:num_crack_plane_nodes) = zero
c
c                            allocate: (for overshoot control algorithm)
c                                  old_angles_at_front -- Stores the old
c                                       angle at each crack front node for
c                                       each of its neighbors.
c                                  overshoot_allocated -- logical variable
c                                       indicating if overshoot data
c                                       structures have been allocated
c
      case( 8 ) ! dowhat
         if( allocated(old_angles_at_front) )
     &       deallocate( old_angles_at_front )
         allocate( old_angles_at_front(num_crack_plane_nodes,
     &             mxconn) )
         old_angles_at_front(1:num_crack_plane_nodes,1:mxconn) = zero
         overshoot_allocated = .true.
c
c
c                            allocate: (for step size control)
c                                g_stp_cntrl_allocated: indicates that
c                                     the structures for step size control
c                                     for crack growth or have been
c                                     allocated.
c                                - for gurson crack growth:
c                                     old_porosity(num_kill_elem) --
c                                       stores the porosity at the previous
c                                       step for comparison to current
c                                       porosity
c                                - for smcs crack growth:
c                                     old_plast_strain(num_kill_elem) --
c                                       stores the plastic strain at the
c                                       previous step for comparison to
c                                       current plast strain
c                                - for cohesive crack growth:
c                                     old_deff(num_kill_elem) --
c                                       stores the effective (relative)
c                                       displacement at the previous step
c                                       for comparison to the current value
c
      case( 9 ) ! dowhat
         select case( crack_growth_type )
           case( 1 )   !   crack_growth_type
c
c            gurson growth, initialize storage of old porosities
c
             if ( allocated( old_porosity ) ) deallocate( old_porosity )
             allocate( old_porosity( num_kill_elem ) )
c
c             Now initialize the old_porosity with the current
c             porosity values...  If we are reallocating due to a
c             restart, then don't bother to fill the array;
c             we will be refilling from the stored values. If
c             we are here before the computation of the first
c             load step, then fill array with the initial porosity.
c
             if ( g_stp_cntrl_allocated ) return
             if ( debug ) write(out,*) '>>> initializing porosity'
             if ( .not. incflg ) then
                do elem  = 1, noelem
                   elem_ptr = dam_ptr(elem)
                   if( elem_ptr .eq. 0 ) cycle
                   old_porosity(elem_ptr) = props(26,elem)
                end do
             else
                do elem  = 1, noelem
                   elem_ptr = dam_ptr(elem)
                   if ( elem_ptr .eq. 0 ) cycle
                   call dam_param( elem, duml, debug, porosity, dumd1,
     &                     dumd2, dumd3, dumd4, duml, dumd5,
     &                     dumd6, dumd7, dumd8, dumd9)
                   old_porosity( elem_ptr ) = porosity
                   if( debug )
     &                  write(out,'("   Poros. for:",i5,"=",e13.6)')
     &                  elem, porosity
                end do   ! elem
             end if ! incflg
c
           case( 3 )  !   crack_growth_type
c
c            crack growth is smcs, initialize storage of
c            old plastic strains
c
              if( allocated( old_plast_strain ) )
     &            deallocate( old_plast_strain )
              allocate( old_plast_strain( num_kill_elem ) )
c
c             Now initialize old_plast_strain with the current
c             plastic strain values...  if we are reallocating
c             due to a restart, then don't bother to fill the array;
c             we will be refilling from the stored values. If
c             we are here before the computation of the first
c             load step, then fill array with zeroes.
c
c
              if( g_stp_cntrl_allocated ) return
              if( .not. incflg ) then
                do elem_ptr = 1, num_kill_elem
                   old_plast_strain(elem_ptr) = zero
                end do
              else
                if( debug ) write(out,*) '>>> init plast_strain'
                do elem = 1, noelem
                   elem_ptr = dam_ptr( elem )
                   if( elem_ptr .eq. 0 ) cycle
                   call dam_param( elem, duml, debug, dumd5,
     &                     plast_strain, dumd2, dumd3, dumd4, duml,
     &                     dumd5, dumd6, dumd7, dumd8, dumd9 )
                   old_plast_strain( elem_ptr ) = plast_strain
                   if( debug )
     &                  write(out,'("   Pl. e for:",i5,"is",e13.6)')
     &                       elem, plast_strain
                end do
              end if  ! incflg
c
           case( 4 ) ! crack_growth_type
c
c             crack growth using cohesive elements
c
              if( allocated(old_deff) ) deallocate(old_deff)
              allocate( old_deff(num_kill_elem) )
c
c             Now initialize the old_deff with the current
c             deff values...  If we are reallocating due to a
c             restart, then don't bother to fill the array;
c             we will be refilling from the stored values. If
c             we are here before the computation of the first
c             load step, then fill array with zero
c
              if( g_stp_cntrl_allocated ) return
              if( debug ) write(out,*) '>>> initializing deff'
              if( .not. incflg ) then
                 do elem  = 1, noelem
                   elem_ptr = dam_ptr(elem)
                   if( elem_ptr .eq. 0 ) cycle
                   old_deff(elem_ptr) = zero
                 end do
              else
                 do elem  = 1, noelem
                    elem_ptr = dam_ptr(elem)
                    if( elem_ptr .eq. 0 ) cycle
                    call dam_param_cohes( elem, duml, debug, values, 1 )
                    old_deff(elem_ptr) = values(6)
                    if( debug )
     &                 write(out,'("   deff. for:",i5,"=",e13.6)')
     &                       elem, values(6)
                 end do
              end if    ! incflg
c
           case default  !  crack_growth_type
             write(out,*) '>>> invalid case in allocate_damage, type 9'
             call die_abort
c
         end select   !   crack_growth_type
c
      g_stp_cntrl_allocated = .true.

c
c
c                            allocate: (constant front growth)
c                                  master_nodes -- holds the list of the
c                                       master nodes, one per crack front,
c                                       which will control the growth of
c                                       its front
c
      case( 10 ) !  dowhat
         if( allocated(master_nodes) ) deallocate( master_nodes )
         allocate( master_nodes(num_crack_fronts) )
         master_nodes(1:num_crack_fronts) = zero
c
c
c                            allocate: (constant front growth)
c                                  crack_front_list -- holds the list of the
c                                       nodes on each of the crack fronts.
c                                  master_lines -- holds a list of
c                                       num_nodes_back nodes of the nodes
c                                       on the master line behind each
c                                       master node.
c
      case( 11 )  ! dowhat
c
c           calculate number of nodes in between growth increments.
c
         num_nodes_grwinc = max( init_ctoa_dist,ctoa_dist ) /
     &                      char_length + 1
         num_nodes_grwinc = max( num_nodes_grwinc,100 )
c
         if( allocated(crack_front_list) )
     &       deallocate( crack_front_list )
         allocate( crack_front_list(num_crack_fronts * num_nodes_grwinc,
     &             num_nodes_thick) )
         crack_front_list(1:num_crack_fronts*num_nodes_grwinc,
     &                    1:num_nodes_thick) = 0
c
         if( allocated(master_lines) ) deallocate( master_lines )
         allocate( master_lines(num_crack_fronts, num_nodes_back + 1) )
         master_lines(1:num_crack_fronts, 1:num_nodes_back + 1) = 0
c
c                            allocate: ( dam_ptr )
c
       case( 12 )  ! dowhat
          if( allocated(dam_ptr) )  deallocate( dam_ptr )
          allocate( dam_ptr(noelem) )
          dam_ptr(1:noelem) = 0
c
c                            SMCS data
c
       case( 13 )  ! dowhat
          if( allocated( smcs_weighted_T ) )
     &        deallocate( smcs_weighted_T )
          if( allocated( smcs_old_epsplas ) )
     &       deallocate( smcs_old_epsplas )
          if( allocated( smcs_weighted_zeta ) )
     &        deallocate( smcs_weighted_zeta  )
          if( allocated( smcs_weighted_bar_theta ) )
     &        deallocate( smcs_weighted_bar_theta  )
          allocate( smcs_weighted_T(num_kill_elem),
     &              smcs_old_epsplas(num_kill_elem),
     &              smcs_weighted_zeta(num_kill_elem),
     &              smcs_weighted_bar_theta(num_kill_elem) )
          do i = 1, num_kill_elem
            smcs_weighted_T(i)          = zero
            smcs_old_epsplas(i)         = zero
            smcs_weighted_zeta(i)       = zero
            smcs_weighted_bar_theta(i)  = zero
          end do
c
      end select  !  dowhat
c
      return
      end


text:
This account appears to be an organizational account for a Woodlandbased web design company.

20120221 12:01:32 nbsp Welcome to the Wiki/Business Owner Welcome to the Wiki!  Please take a moment to read over that welcome page  it has useful information for business owners/representatives.  Some of the changes you made on the Executive Cleaners page got altered or undone.  As that welcome page explains, the wiki cant be used promotionally due to its nonprofit arrangement with the IRS.  There isnt really a set of hardandfast rules, but there are some wiki norms we use to try to keep things running smoothly.  One of those is having links to category pages.  There were also some concerns about some of the images you uploaded and possible problems with copyrights. Users/TomGarberson

20120221 12:09:41 nbsp It looks like youre taking money to create a wiki page for someone. The wiki is not like Facebook. This is a community resource, not a phone directory. You are totally tone deaf to this community if you think what youre doing fits in here. Users/WilliamLewis

20120221 12:25:49 nbsp William points out that it looks like youre associated with the Woodlandbased web design company of the same name.  Please see the notice above regarding organizational accounts. Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module type_solverdata
#include <messenger.h>
    use mod_kinds,                      only: rk,ik
    use mod_constants,                  only: NFACES
    use mod_string,                     only: string_t
    use type_chidg_vector,               only: chidg_vector_t
    use type_chidg_matrix,               only: chidg_matrix_t
    use type_mesh,                  only: mesh_t
    use type_function_status,           only: function_status_t
    use type_equationset_function_data, only: equationset_function_data_t
    use type_element_info,              only: element_info_t
    use type_face_info,                 only: face_info_t
    use type_function_info,             only: function_info_t
    implicit none

    !> Container for solver data.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   3/15/2016
    !!
    !!
    !------------------------------------------------------------------------------------------
    type, public  :: solverdata_t

        !
        ! Base solver data
        !
        type(chidg_vector_t)             :: q              ! Solution vector
        type(chidg_vector_t)             :: dq             ! Change in solution vector
        type(chidg_vector_t)             :: rhs            ! Residual of the spatial scheme
        type(chidg_matrix_t)             :: lhs            ! Linearization of the spatial scheme

        !
        ! Container for reading data
        !
        type(chidg_vector_t)            :: q_in

        !
        ! Container for storing data for post processing
        !
        type(chidg_vector_t)            :: q_out

        !
        ! Auxiliary fields
        !
        type(string_t),         allocatable :: auxiliary_field_name(:)
        type(chidg_vector_t),   allocatable :: auxiliary_field(:)

        !
        ! Time information
        !
        !real(rk)                        :: t               ! Global time
        real(rk),   allocatable         :: dt(:,:)         ! Element-local time-step, (ndomains,maxelems)

        !
        ! Function registration
        !
        type(function_status_t)         :: function_status ! Status of function residuals and linearizations

        logical                         :: solverInitialized = .false.

    contains

        generic, public     :: init => init_base
        procedure, private  :: init_base

        procedure           :: add_auxiliary_field

        procedure           :: nauxiliary_fields
        procedure           :: get_auxiliary_field_index
        procedure           :: get_auxiliary_field_name

        procedure           :: release

    end type solverdata_t
    !*******************************************************************************************

contains

    !>  Initialize solver base data structures
    !!      - allocate and initialize q, dq, rhs, and linearization.
    !!      - Should be called by specialized 'init' procedure for derived solvers.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!  @param[in]  mesh            Array of mesh_t instances which define storage requirements.
    !!  @param[in]  bcset_coupling  Array of bcset_coupling instances which describe the
    !!                              coupling of elements in bcs.
    !!  @param[in]  function_data   Array of containers that hold information on number of
    !!                              each function in eqnset.
    !!
    !-------------------------------------------------------------------------------------------
    subroutine init_base(self,mesh,function_data)
        class(solverdata_t),                intent(inout)   :: self
        type(mesh_t),                       intent(inout)   :: mesh
        type(equationset_function_data_t),  intent(in)      :: function_data(:)

        integer(ik) :: ierr, ndom, maxelems, idom
        logical     :: increase_maxelems = .false.

        ! Initialize vectors
        call self%q%init(    mesh,mesh%ntime_)
        call self%dq%init(   mesh,mesh%ntime_)
        call self%rhs%init(  mesh,mesh%ntime_)
        call self%q_in%init( mesh,mesh%ntime_)
        call self%q_out%init(mesh,mesh%ntime_)

        ! Initialize matrix and parallel recv data
        call self%lhs%init(mesh,'full')
        call self%lhs%init_recv(self%rhs)

        !
        ! Find maximum number of elements in any domain
        !
        ndom = mesh%ndomains()
        maxelems = 0
        do idom = 1,ndom

            increase_maxelems = ( mesh%domain(idom)%nelem > maxelems )

            if (increase_maxelems) then
                maxelems = mesh%domain(idom)%nelem
            end if

        end do

        !
        ! Allocate timestep storage
        !
        if (allocated(self%dt)) deallocate(self%dt)
        allocate(self%dt(ndom,maxelems),stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Initialize storage on flux and linearization registration
        !
        call self%function_status%init( mesh, function_data)

        !
        ! Confirm solver initialization
        !
        self%solverInitialized = .true.

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

    !>  Add chidg_vector for storing an auxiliary field for the problem.
    !!
    !!  One could call this as:
    !!      call solverdata%add_auxiliary_field('my field')
    !!
    !!  which would just add an empty chidg_vector for storing the auxiliary field
    !!
    !!  One could also call this as:
    !!      call solverdata%add_auxiliary_field('my field', my_vector)
    !!
    !!  which would create space for a new auxiliary vector and assign the incoming
    !!  chidg_vector, my_vector, to the field.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/1/2016
    !!
    !!
    !!
    !------------------------------------------------------------------------------------------
    subroutine add_auxiliary_field(self,fieldname,auxiliary_vector)
        class(solverdata_t),    intent(inout)           :: self
        character(*),           intent(in)              :: fieldname
        type(chidg_vector_t),   intent(in), optional    :: auxiliary_vector

        integer(ik) :: naux_vectors, ierr

        type(string_t),         allocatable :: temp_names(:)
        type(chidg_vector_t),   allocatable :: temp_vectors(:)

        !
        ! Resize array storage
        !
        if (allocated(self%auxiliary_field)) then
            naux_vectors = size(self%auxiliary_field) + 1
        else
            naux_vectors = 1
        end if

        allocate(temp_names(naux_vectors), &
                 temp_vectors(naux_vectors), stat=ierr)
        if (ierr /= 0) call AllocationError

        !
        ! Copy previously added auxiliary fields to new array
        !
        if (naux_vectors > 1) then
            temp_names(1:size(self%auxiliary_field_name)) = self%auxiliary_field_name(1:size(self%auxiliary_field_name))
            temp_vectors(1:size(self%auxiliary_field))    = self%auxiliary_field(1:size(self%auxiliary_field))
        end if

        !
        ! Set field name
        !
        temp_names(naux_vectors) = string_t(trim(fieldname))

        !
        ! Move resized temp allocation back to solverdata_t container.
        !
        call move_alloc(temp_names,self%auxiliary_field_name)
        call move_alloc(temp_vectors,self%auxiliary_field)

        !
        ! Store incoming vector if present
        !
        if (present(auxiliary_vector)) then
            self%auxiliary_field(naux_vectors) = auxiliary_vector
        end if

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

    !>  Given the name of an auxiliary field, return the index of the vector containing
    !!  the field in self%auxiliary_field.
    !!
    !!  The field represented as a chidg_vector would then be accessed as:
    !!      solverdata%auxiliary_field(index)
    !!
    !!  If not found, returns 0.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   11/1/2016
    !!
    !-----------------------------------------------------------------------------------------
    function get_auxiliary_field_index(self,fieldname) result(field_index)
        class(solverdata_t),    intent(in)  :: self
        character(*),           intent(in)  :: fieldname

        integer(ik)                 :: field_index, ifield

        !
        ! Loop through names to try and find field
        !
        if (allocated(self%auxiliary_field_name)) then

            field_index = 0
            do ifield = 1,size(self%auxiliary_field_name)
                if (self%auxiliary_field_name(ifield)%get() == trim(fieldname)) then
                    field_index = ifield
                    exit
                end if
            end do

        else

            field_index = 0

        end if

    end function get_auxiliary_field_index
    !******************************************************************************************

    !>  Given the index of an auxiliary field, return the name of the field the auxiliary
    !!  vector is representing.
    !!
    !!  The field represented as a chidg_vector would then be accessed as:
    !!      solverdata%auxiliary_field(index)
    !!
    !!  @author Nathan A. Wukie
    !!  @date   12/21/2016
    !!
    !-----------------------------------------------------------------------------------------
    function get_auxiliary_field_name(self,field_index) result(field_name)
        class(solverdata_t),    intent(in)  :: self
        integer(ik),            intent(in)  :: field_index

        character(:),   allocatable :: field_name, user_msg

        ! Check bounds
        user_msg = "solverdata%get_auxiliary_field_name: Index is out of bounds."
        if (field_index > self%nauxiliary_fields()) call chidg_signal(FATAL,user_msg)

        ! Return the appropriate field name
        field_name = trim(self%auxiliary_field_name(field_index)%get())

    end function get_auxiliary_field_name
    !******************************************************************************************

    !>  Return the number of auxiliary fields that have been added to the solverdata_t
    !!  container.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   12/21/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    function nauxiliary_fields(self) result(nfields)
        class(solverdata_t),    intent(in)  :: self

        integer(ik) :: nfields

        if (allocated(self%auxiliary_field)) then
            nfields = size(self%auxiliary_field)
        else
            nfields = 0
        end if

    end function nauxiliary_fields
    !*****************************************************************************************

    !>  Release allocated data.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   3/3/2017
    !!
    !!
    !----------------------------------------------------------------------------------------
    subroutine release(self)
       class(solverdata_t), intent(inout)   :: self

        ! Release chidg_vector data
        call self%q%release()
        call self%dq%release()
        call self%q_in%release()
        call self%rhs%release()

        ! Release chidg_matrix data
        call self%lhs%release()

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

end module type_solverdata

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
program main

 character(len=:),allocatable :: str

 str = 'hello'
 str = str//' world'

 write(*,*) str

end program main

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


! Chemfiles, an efficient IO library for chemistry file formats
! Copyright (C) 2015-2019 Guillaume Fraux -- BSD licence
module chemfiles_selection
    use iso_c_binding
    use chemfiles_ffi
    use chemfiles_strings
    use chemfiles_frame
    implicit none

    private
    public :: chfl_selection

    type, extends(chfl_ptr) :: chfl_selection
    contains
        procedure :: init
        procedure :: copy
        procedure :: size => get_size
        procedure :: string
        procedure :: evaluate
        procedure :: matches
    end type

contains
    subroutine init(this, selection, status)
        implicit none
        class(chfl_selection), intent(inout) :: this
        character(len=*), intent(in) :: selection
        integer(chfl_status), intent(out), optional :: status

        call this%unsafe_set_ptr(c_chfl_selection(f_to_c_str(selection)), status)
    end subroutine

    subroutine copy(this, selection, status)
        implicit none
        class(chfl_selection), intent(inout) :: this
        class(chfl_selection), intent(in) :: selection
        integer(chfl_status), intent(out), optional :: status

        call this%unsafe_set_ptr(c_chfl_selection_copy(selection%unsafe_const_ptr()), status)
    end subroutine

    integer(kind=c_int64_t) function get_size(this, status)
        implicit none
        class(chfl_selection), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_selection_size(this%unsafe_const_ptr(), get_size)

        if (present(status)) then
            status = status_tmp
        end if
    end function

    character(len=CHFL_STRING_LENGTH) function string(this, status)
        implicit none
        class(chfl_selection), intent(in) :: this
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_selection_string(                                 &
            this%unsafe_const_ptr(), string, int(CHFL_STRING_LENGTH, c_int64_t)   &
        )
        string = rm_null_in_str(string)
        if (present(status)) then
            status = status_tmp
        end if
    end function

    subroutine evaluate(this, frame, count, status)
        implicit none
        class(chfl_selection), intent(inout) :: this
        class(chfl_frame), intent(in) :: frame
        integer(kind=c_int64_t), intent(out) :: count
        integer(chfl_status), intent(out), optional :: status
        integer(chfl_status) :: status_tmp

        status_tmp = c_chfl_selection_evaluate(this%unsafe_ptr(), frame%unsafe_const_ptr(), count)

        if (present(status)) then
            status = status_tmp
        end if
    end subroutine

    subroutine matches(this, data, status)
        implicit none
        class(chfl_selection), intent(in) :: this
        type(chfl_match), intent(out), dimension(:), target :: data
        integer(chfl_status), intent(out), optional :: status
        integer(kind=c_int64_t) :: count
        integer(chfl_status) :: status_tmp

        count = size(data, 1, c_int64_t)
        status_tmp = c_chfl_selection_matches(this%unsafe_const_ptr(), c_loc(data), count)

        if (present(status)) then
            status = status_tmp
        end if
    end subroutine
end module


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

C   --*** PRELB *** (BLOT) Display database element blocks
C   --
C   --PRELB displays the element blocks.
C   --
C   --Parameters:
C   --   OPTION - IN - '*' to print all, else print options:
C   --      ' ' to print block summary
C   --      'N' to print element block name
C   --      'V' to print element variable truth table
C   --      'C' to print connectivity
C   --      'A' to print attributes
C   --   NOUT - IN - the output file, <=0 for standard
C   --   NELBLK - IN - the number of element blocks
C   --   NLISEL - IN - the number of selected elements by block
C   --   LISEL - IN - the indices of the selected elements by block
C   --   IDELB - IN - the element block ID for each block
C   --   LENE - IN - the cumulative element counts by element block
C   --   NUMLNK - IN - the number of nodes per element for each block
C   --   NUMATR - IN - the number of attributes for each block
C   --   LINK - IN - the connectivity array for all blocks
C   --   ATRIB - IN - the attribute array for all blocks
C   --   EBTYPE - IN - the names of the element block types
C   --   NVAREL - IN - the number of element variables
C   --   NAMEEV - IN - the names of the element variables
C   --   ISEVOK - IN - the element block variable truth table;
C   --      variable i of block j exists iff ISEVOK(i,j) is NOT 0
C   --   LISEV - SCRATCH - size = NVAREL (if 'V' in OPTION)

      include 'exodusII.inc'
      CHARACTER*(*) OPTION
      INTEGER NLISEL(0:NELBLK)
      INTEGER LISEL(0:*)
      INTEGER IDELB(*)
      INTEGER LENE(0:NELBLK)
      INTEGER NUMLNK(*)
      INTEGER NUMATR(*)
      INTEGER LINK(*)
      REAL ATRIB(*)
      CHARACTER*(MXSTLN) EBTYPE(*)
      CHARACTER*(*) EBNAME(*)
      CHARACTER*(*) NAMEEV(*)
      CHARACTER*(*) ATNAME(*)
      INTEGER ISEVOK(NVAREL,*)
      INTEGER LISEV(*)
      INTEGER MAPN(*), MAPE(*)
      LOGICAL DOMAPN, DOMAPE
      LOGICAL ISABRT
      LOGICAL DONAM, DOVTBL, DOCONN, DOATR
      LOGICAL BLK1
      CHARACTER*40 STRA, STRB

      DONAM  = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'N') .GT. 0))
      DOVTBL = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'V') .GT. 0))
      DOCONN = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'C') .GT. 0))
      DOATR  = ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'A') .GT. 0))

      DONAM = .TRUE.
      BLK1 = .TRUE.
      IF (NOUT .GT. 0) THEN
        IF (DOCONN .AND. DOATR) THEN
          WRITE (NOUT, 10020) 'CONNECTIVITY and ATTRIBUTES'
        ELSE IF (DOCONN) THEN
          WRITE (NOUT, 10020) 'CONNECTIVITY'
        ELSE IF (DOATR) THEN
          WRITE (NOUT, 10020) 'ATTRIBUTES'
        ELSE
          WRITE (NOUT, 10020)
        END IF
      END IF

      if (nout .gt. 0) then
        if (domape) write (nout, 10060) 'Element'
        if (domapn) write (nout, 10060) 'Nodal'
      else
        if (domape) write (*, 10060) 'Element'
        if (domapn) write (*, 10060) 'Nodal'
      end if

      WRITE (STRA, 10000, IOSTAT=IDUM) NELBLK
10000 FORMAT ('(#', I10, ')')
      CALL PCKSTR (1, STRA)
      LSTRA = LENSTR (STRA)
      WRITE (STRB, 10010, IOSTAT=IDUM) LENE(NELBLK), LENE(NELBLK)
10010 FORMAT ('(', I10, '..', I10, ')')
      CALL PCKSTR (1, STRB)
      LSTRB = LENSTR (STRB)

      ISLNK = 1
      ISATR = 1
      ISATN = 1

      DO 110 IELB = 1, NELBLK
        NUME = LENE(IELB) - LENE(IELB-1)
        if (ebtype(ielb) .eq. 'nsided' .or.
     *    ebtype(ielb) .eq. 'NSIDED') THEN
          numnod = numlnk(ielb)
        else
          numnod = numlnk(ielb)*nume
        end if
        IF (ISABRT ()) RETURN
        IF (NLISEL(IELB) .GT. 0) THEN
          IEL = LENE(IELB-1)+1
          LEL = LENE(IELB)

          IF (BLK1 .OR. DOCONN .OR. DOATR) THEN
            BLK1 = .FALSE.
            IF (NOUT .GT. 0) THEN
              WRITE (NOUT, *)
            ELSE
              WRITE (*, *)
            END IF
          END IF

          WRITE (STRA, 10000, IOSTAT=IDUM) IELB
          CALL PCKSTR (1, STRA)
          WRITE (STRB, 10010, IOSTAT=IDUM) IEL, LEL
          CALL PCKSTR (1, STRB)
          IF (NOUT .GT. 0) THEN
            WRITE (NOUT, 10030, IOSTAT=IDUM)
     &        IDELB(IELB), STRA(:LSTRA),
     &        NUME, STRB(:LSTRB), NUMLNK(IELB), NUMATR(IELB)
          ELSE
            WRITE (*, 10030, IOSTAT=IDUM)
     &        IDELB(IELB), STRA(:LSTRA),
     &        NUME, STRB(:LSTRB), NUMLNK(IELB), NUMATR(IELB)
          END IF

          IF (DONAM) THEN
            LNAM = LENSTR(EBTYPE(IELB))
            LNM =  LENSTR(EBNAME(IELB))
            IF (NOUT .GT. 0) THEN
              WRITE (NOUT, 10040) EBNAME(IELB)(:LNM),
     $          EBTYPE(IELB)(:LNAM)
            ELSE
              WRITE (*, 10040) EBNAME(IELB)(:LNM), EBTYPE(IELB)(:LNAM)
            END IF
            if (numatr(ielb) .gt. 0) then
              mxnam = 0
              DO I = 1, numatr(ielb)
                LNAM = lenstr(atname(isatn+i-1))
                mxnam = max(lnam, mxnam)
              END DO
              if (mxnam .gt. 1) then
                IF (NOUT .GT. 0) THEN
                  WRITE (NOUT, 10090, IOSTAT=IDUM)
     &              (ATNAME(I)(:mxnam),
     *              I=ISATN,ISATN+NUMATR(IELB)-1)
                ELSE
                  WRITE (*, 10090, IOSTAT=IDUM)
     &              (ATNAME(I)(:mxnam),
     *              I=ISATN,ISATN+NUMATR(IELB)-1)
                END IF
              end if
            end if
          END IF

          IF (DOVTBL) THEN
            NSEL = 0
            mxnam = 0
            DO 100 I = 1, NVAREL
              IF (ISEVOK(I,IELB) .NE. 0) THEN
                NSEL = NSEL + 1
                LISEV(NSEL) = I
                LNAM = lenstr(nameev(i))
                mxnam = max(lnam, mxnam)
              END IF
 100        CONTINUE

            IF (NOUT .GT. 0) THEN
              if (mxnam .gt. 24) then
                WRITE (NOUT, 10070, IOSTAT=IDUM)
     &            (NAMEEV(LISEV(I))(:mxnam), I=1,NSEL)
              else
                WRITE (NOUT, 10050, IOSTAT=IDUM)
     &            (NAMEEV(LISEV(I))(:mxnam), I=1,NSEL)
              endif
              write (nout, 10080)
            ELSE
              if (mxnam .gt. 24) then
                WRITE (*, 10070, IOSTAT=IDUM)
     &            (NAMEEV(LISEV(I))(:mxnam), I=1,NSEL)
              else
                WRITE (*, 10050, IOSTAT=IDUM)
     &            (NAMEEV(LISEV(I))(:mxnam), I=1,NSEL)
              endif
              write (*, 10080)
            END IF
          END IF

          IF (DOCONN .OR. DOATR) THEN
            if (ebtype(ielb) .eq. 'nsided' .or.
     *        ebtype(ielb) .eq. 'NSIDED') THEN
              CALL PREBN (OPTION, NOUT, IELB,
     *          nume, IDELB(IELB),
     *          IEL-1, NLISEL(IELB), LISEL(IEL),
     *          NUMATR(IELB),
     &          LINK(ISLNK), ATRIB(ISATR), max(1,numatr(ielb)),
     *          MAPN, DOMAPN, MAPE, DOMAPE)
            else
              CALL PREB1 (OPTION, NOUT, IEL-1,
     &          NLISEL(IELB), LISEL(IEL), NUMLNK(IELB), NUMATR(IELB),
     &          LINK(ISLNK), ATRIB(ISATR), max(1,numatr(ielb)),
     *          MAPN, DOMAPN, MAPE, DOMAPE)
            end if
            IF (ISABRT ()) RETURN
          END IF
        END IF

        IF (DOCONN) ISLNK = ISLNK + numnod
        IF (DOATR)  ISATR = ISATR + NUME * NUMATR(IELB)
        ISATN = ISATN + NUMATR(IELB)
 110  CONTINUE

      RETURN

10020 FORMAT (/, 1X, 'ELEMENT BLOCKS', :, ' - ', A)
C ... NOTE: Normal elements hava a mximum of <100 nodes, but superelements
C           Or other "strange" elements may have lots of nodes; keep format high
10030 FORMAT (1X, 'Block', I10, 1X, A, ':',
     &  I10, ' elements', 1X, A,
     &  I6, '-node', I4, ' attributes')
10040 FORMAT (2X, 'Element block name = "',A
     $  ,'", Element type = "', A, '"')
10050 FORMAT ((2x,4(2X, A)))
10060 FORMAT (1x, A, ' IDs are global')
10070 FORMAT ((2x,3(2X, A)))
10080 FORMAT (1X)
10090 FORMAT (2X, 'Attributes: ', 10(2X, A))
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT