text:
MODULE generic_procedure_module
IMPLICIT NONE
! Declare type vector
TYPE :: vector
  REAL :: x
  REAL :: y
CONTAINS
  GENERIC :: add => vector_plus_vector, vector_plus_scalar
  PROCEDURE, PASS :: vector_plus_scalar
  PROCEDURE, PASS :: vector_plus_vector
END TYPE vector

CONTAINS
  TYPE (vector) FUNCTION vector_plus_vector(this, v2)
  IMPLICIT NONE
  CLASS(vector), INTENT(IN) :: this
  CLASS(vector), INTENT(IN) :: v2
  vector_plus_vector%x = this%x + v2%x
  vector_plus_vector%y = this%y + v2%y
  END FUNCTION vector_plus_vector

  TYPE (vector) FUNCTION vector_plus_scalar(this, scalar)
  IMPLICIT NONE
  CLASS(vector), INTENT(IN) :: this
  REAL, INTENT(IN) :: scalar
  vector_plus_scalar%x = this%x + scalar
  vector_plus_scalar%y = this%y + scalar
  END FUNCTION vector_plus_scalar
END MODULE generic_procedure_module

PROGRAM test_generic_procedures
USE generic_procedure_module
IMPLICIT NONE
TYPE(vector) :: v1
TYPE(vector) :: v2
REAL :: s

WRITE(*, *) 'ENTER FIRST VECTOR (X, Y):'
READ(*, *) v1%x, v1%y
WRITE(*, *) 'ENTER SECOND VECTOR (X, Y):'
READ(*, *) v2%x, v2%y
WRITE(*, *) 'ENTER A SCALAR:'
READ(*, *) s
WRITE(*, 1000) v1%add(v2)
1000 FORMAT ('The sum of the vectors is (',F8.2,',',F8.2,')')
WRITE(*, 1010) v1%add(s)
1010 FORMAT ('The sum of the vector and the scalar is (',F8.2,',',F8.2,')')
END PROGRAM test_generic_procedures

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice


text:
        SUBROUTINE UMD2CO (N, NZ, TRANSA, XX, XSIZE, INFO, ICNTL,
     $          II, ISIZE, W, WP, WHO)
        INTEGER ISIZE, II (ISIZE), N, NZ, W (N), WP (N+1), INFO (40),
     $          ICNTL (20), XSIZE, WHO
        DOUBLE PRECISION
     $          XX (XSIZE)
        LOGICAL TRANSA

C=== UMD2CO ============================================================
C
C  Unsymmetric-pattern MultiFrontal Package (UMFPACK). Version 2.2d
C  Copyright (C) 1997, Timothy A. Davis, University of Florida, USA.
C  ALL RIGHTS RESERVED.
C  Joint work with Iain S. Duff, Rutherford Appleton Laboratory, UK.
C  July 7, 1997. Work supported by the National Science Foundation
C  (DMS-9223088 and DMS-9504974) and the State of Florida; and by CRAY
C  Research Inc. through the allocation of supercomputing resources.

C***********************************************************************
C* NOTICE:  "The UMFPACK Package may be used SOLELY for educational,   *
C* research, and benchmarking purposes by non-profit organizations and *
C* the U.S. government.  Commercial and other organizations may make   *
C* use of UMFPACK SOLELY for benchmarking purposes only.  UMFPACK may  *
C* be modified by or on behalf of the User for such use but at no time *
C* shall UMFPACK or any such modified version of UMFPACK become the    *
C* property of the User.  UMFPACK is provided without warranty of any  *
C* kind, either expressed or implied.  Neither the Authors nor their   *
C* employers shall be liable for any direct or consequential loss or   *
C* damage whatsoever arising out of the use or misuse of UMFPACK by    *
C* the User.  UMFPACK must not be sold.  You may make copies of        *
C* UMFPACK, but this NOTICE and the Copyright notice must appear in    *
C* all copies.  Any other use of UMFPACK requires written permission.  *
C* Your use of UMFPACK is an implicit agreement to these conditions."  *
C*                                                                     *
C* The MA38 Package in Release 12 of the Harwell Subroutine Library    *
C* (HSL) has equivalent functionality (and identical calling interface)*
C* as UMFPACK (the HSL has single and double precision versions only,  *
C* however).  It is available for commercial use.   Technical reports, *
C* information on HSL, and matrices are available via the World Wide   *
C* Web at http://www.cis.rl.ac.uk/struct/ARCD/NUM.html, or by          *
C* anonymous ftp at seamus.cc.rl.ac.uk/pub.  Also contact Dr. Scott    *
C* Roberts, Harwell Subroutine Library, B 552, AEA Technology,         *
C* Harwell, Didcot, Oxon OX11 0RA, England.                            *
C* telephone (44) 1235 434988, fax (44) 1235 434136                    *
C* email Scott.Roberts@aeat.co.uk, who will provide details of price   *
C* and conditions of use.                                              *
C***********************************************************************

C=======================================================================
C  NOT USER-CALLABLE.

C=======================================================================
C  DESCRIPTION:
C=======================================================================
C
C  Convert input matrix (II,XX,n,nz) into from triplet form to column-
C  oriented form, optionally transposing the input matrix (single and
C  double precision versions only).  Remove invalid entries and
C  duplicate entries.

C=======================================================================
C  INPUT:
C=======================================================================
C
C       n:              size of the matrix
C       nz:             number of nonzeros in the input matrix

C       transa:         if true then transpose input matrix

C       XX (1..nz):     values of triplet form
C       xsize:          size of XX, must be >= 2*nz
C       II (1..2*nz):   row and col indices of triplet form
C       isize:          size of II, must be >= max (2*nz,n+1) + nz
C       Icntl:          integer control parameters
C       who:            who called UMD2CO, 1: UMD2FA, 2: UMD2RF
C
C       II must be at least of size (nz + max (2*nz, n+1))
C       XX must be at least of size (nz + max (  nz, n+1))
C
C       input triplet matrix:

C          if (transa) is false:
C               II (p)          i, row index, for p = 1..nz
C               II (nz+p)       j, col index
C               XX (p)          a_ij
C          if (transa) is true:
C               II (p)          j, col index, for p = 1..nz
C               II (nz+p)       i, row index
C               XX (p)          a_ij

C=======================================================================
C  WORKSPACE:
C=======================================================================
C
C       W (1..n)

C=======================================================================
C  OUTPUT:
C=======================================================================
C
C       nz:             number of nonzeros in the output matrix,
C                       after removing invalid entries, and summing up
C                       duplicate entries
C       II (n+2..nz+n+1): row indices in column-form
C       XX (1..nz):     values in column-form.
C       Info (1):       error flag
C       Info (3):       invalid entries
C       Info (2):       duplicate entries
C       Info (5):       remaining valid entries
C       Info (6):       remaining valid entries
C       Info (7):       0
C       Wp (1..n+1)     column pointers for column form
C       II (1..n+1)     column pointers for column form

C=======================================================================
C  SUBROUTINES AND FUNCTIONS CALLED / CALLED BY:
C=======================================================================
C
C       called by subroutines:  UMD2FA, UMD2RF
C       subroutines called:     UMD2ER, UMD2P2
C       functions called:       MAX
        INTRINSIC MAX

C=======================================================================
C  LOCAL SCALARS:
C=======================================================================

        INTEGER ROW, COL, PDEST, P, NZ1, PCOL, IP, XP, IO, PRL, NINVLD,
     $          NDUPL, I
        LOGICAL PR3

C  row:     row index
C  col:     column index
C  pdest:   location of an entry in the column-form, for dupl. removal
C  p:       pointer
C  nz1:     number of entries after removing invalid or duplic. entries
C  pcol:    column col starts here after duplicates removed
C  ip:      column-form copy of matrix placed in II (ip...ip+nz-1)
C  xp:      column-form copy of matrix placed in XX (xp...xp+nz-1)
C  ninvld:  number of invalid entries
C  ndupl:   number of duplicate entries
C  i:       a row index if transa is true, a column index otherwise
C  io:      I/O unit for warning messages (for invalid or dupl. entries)
C  prl:     printing level
C  pr3:     true if printing invalid and duplicate entries

C=======================================================================
C  EXECUTABLE STATEMENTS:
C=======================================================================

C-----------------------------------------------------------------------
C  get arguments and check memory sizes
C-----------------------------------------------------------------------

        IO = ICNTL (2)
        PRL = ICNTL (3)
        PR3 = PRL .GE. 3 .AND. IO .GE. 0

C-----------------------------------------------------------------------
C  count nonzeros in columns and check for invalid entries
C-----------------------------------------------------------------------

        NINVLD = 0
        NDUPL = 0
        DO 10 COL = 1, N
           W (COL) = 0
10      CONTINUE
        NZ1 = NZ
        DO 20 P = NZ, 1, -1
           ROW = II (P)
           COL = II (NZ+P)
           IF (ROW.LT.1.OR.ROW.GT.N.OR.COL.LT.1.OR.COL.GT.N) THEN
C             this entry is invalid - delete it
              IF (PR3) THEN
C                print the offending entry on the diagnostic I/O unit
                 CALL UMD2P2 (WHO, 99, ROW, COL, XX(P), IO)
              ENDIF
              II (P)    = II (NZ1)
              II (NZ+P) = II (NZ+NZ1)
              XX (P)    = XX (NZ1)
              NZ1 = NZ1 - 1
           ELSE

              IF (TRANSA) THEN
C                factorizing A transpose
                 W (ROW) = W (ROW) + 1
              ELSE
C                factorizing A
                 W (COL) = W (COL) + 1
              ENDIF

           ENDIF
20      CONTINUE
        NINVLD = NZ - NZ1
        IF (NINVLD .NE. 0) THEN
C          invalid entries found - set warning flag and continue
           CALL UMD2ER (WHO, ICNTL, INFO, 1, NINVLD)
        ENDIF

C-----------------------------------------------------------------------
C  convert triplet form to column-form
C-----------------------------------------------------------------------

        WP (1) = 1
        DO 30 I = 1, N
           WP (I+1) = WP (I) + W (I)
30      CONTINUE
        DO 40 I = 1, N
           W (I) = WP (I)
40      CONTINUE

C       ----------------------------------------------------------------
C       construct column-form in II (2*nz+1..3*nz) and XX (nz+1..2*nz)
C       ----------------------------------------------------------------

        IP = MAX (2*NZ, N+1)
        XP = NZ

        IF (TRANSA) THEN
           DO 50 P = 1, NZ1
              ROW = II (P)
              COL = II (NZ+P)
              II (IP + W (ROW)) = COL
              XX (XP + W (ROW)) = XX (P)
              W (ROW) = W (ROW) + 1
50         CONTINUE
        ELSE

           DO 60 P = 1, NZ1
              ROW = II (P)
              COL = II (NZ+P)
              II (IP + W (COL)) = ROW
              XX (XP + W (COL)) = XX (P)
              W (COL) = W (COL) + 1
60         CONTINUE

        ENDIF

C       ----------------------------------------------------------------
C       shift the matrix back to II (n+2..nz+n+1) and XX (n+2..nz+n+1)
C       ----------------------------------------------------------------

        NZ = NZ1
CFPP$ NODEPCHK L
        DO 70 P = 1, NZ
           II (N+1+P) = II (IP+P)
           XX (P) = XX (XP+P)
70      CONTINUE

C-----------------------------------------------------------------------
C  remove duplicate entries by adding them up
C-----------------------------------------------------------------------

        DO 80 ROW = 1, N
           W (ROW) = 0
80      CONTINUE
        PDEST = 1
        DO 100 COL = 1, N
           PCOL = PDEST
           DO 90 P = WP (COL), WP (COL+1)-1
              ROW = II (N+1+P)
              IF (W (ROW) .GE. PCOL) THEN
C                this is a duplicate entry
                 XX (W (ROW)) = XX (W (ROW)) + XX (P)
                 IF (PR3) THEN
C                   print the duplicate entry on the diagnostic I/O
C                   unit.  The row and column indices printed reflect
C                   the input matrix.

                    IF (TRANSA) THEN
                       CALL UMD2P2 (WHO, 98, COL, ROW, XX (P), IO)
                    ELSE
                       CALL UMD2P2 (WHO, 98, ROW, COL, XX (P), IO)
                    ENDIF

                 ENDIF
              ELSE
C                this is a new entry, store and record where it is
                 W (ROW) = PDEST
                 IF (PDEST .NE. P) THEN
                    II (N+1+PDEST) = ROW
                    XX (PDEST) = XX (P)
                 ENDIF
                 PDEST = PDEST + 1
              ENDIF
90         CONTINUE
           WP (COL) = PCOL
100     CONTINUE
        WP (N+1) = PDEST
        NZ1 = PDEST - 1
        NDUPL = NZ - NZ1
        IF (NDUPL .NE. 0) THEN
C          duplicate entries found - set warning flag and continue
           CALL UMD2ER (WHO, ICNTL, INFO, 2, NDUPL)
        ENDIF
        NZ = NZ1

C-----------------------------------------------------------------------
C  save column pointers in II (1..n+1)
C-----------------------------------------------------------------------

        DO 110 COL = 1, N+1
           II (COL) = WP (COL)
110     CONTINUE

        INFO (2) = NDUPL
        INFO (3) = NINVLD
        INFO (5) = NZ
        INFO (6) = NZ
        INFO (7) = 0
        IF (NZ .EQ. 0) THEN
C          set error flag if all entries are invalid
           CALL UMD2ER (WHO, ICNTL, INFO, -2, -1)
        ENDIF
        RETURN
        END

ext: f

max_stars_repo_name: kuanshi/ductile-fracture


text:
Name: Russ Bowlus

Age: 38

Occupation: Civil Engineer for the CA DWR, Division of Safety of Dams

Hobbies: Cars and Music

Background: Graduated UCD in 1997 with a BS in Civil Engineering.  Member of the California Aggie Marching BandUh Cal Aggie Marching Banduh! for the coughsixcough years he attended school, playing in the drum section and serving two years on the Officer Council as Assistant Drum Major (94) and Drum Major (95).  Also drove for Unitrans for four years (employee number 390, 2000 driving hours, booyah!) and Fleet Services (driving the UCD charter buses) for three.  Lives in East Davis with his wife (also a UCD graduate), new son, three cars (currently: 98 BMW 328i, 00 Protege, 94 Miata), and no pets.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine shape26h(xi,et,ze,xl,xsj,shp,iflag,konl)
!
!     shape functions and derivatives for a 26-node quadratic
!     isoparametric brick element. -1<=xi,et,ze<=1
!
!     iflag=1: calculate only the value of the shape functions
!     iflag=2: calculate the value of the shape functions and
!              the Jacobian determinant
!     iflag=3: calculate the value of the shape functions, the
!              value of their derivatives w.r.t. the global
!              coordinates and the Jacobian determinant
!
      implicit none
!
      integer i,j,k,iflag,konl(26),jf(3,6),ifaceq(8,6)
!
      real*8 shp(4,26),xs(3,3),xsi(3,3),xl(3,26),shpe(4,26),dd,
     &  dd1,dd2,dd3,xi,et,ze,xsj,omg,omh,omr,opg,oph,opr,
     &  tpgphpr,tmgphpr,tmgmhpr,tpgmhpr,tpgphmr,tmgphmr,tmgmhmr,tpgmhmr,
     &  omgopg,omhoph,omropr,omgmopg,omhmoph,omrmopr,fxi(3),fet(3),
     &  fze(3),dfxi(3),dfet(3),dfze(3)
!
      intent(in) xi,et,ze,xl,iflag,konl
!
      intent(out) shp,xsj
!
      jf=reshape((/2,2,1,2,2,3,2,1,2,3,2,2,2,3,2,1,2,2/),(/3,6/))
!
      ifaceq=reshape(( /4,3,2,1,11,10,9,12,
     &                  5,6,7,8,13,14,15,16,
     &                  1,2,6,5,9,18,13,17,
     &                  2,3,7,6,10,19,14,18,
     &                  3,4,8,7,11,20,15,19,
     &                  4,1,5,8,12,17,16,20/),(/8,6/))
!
!     shape functions in one dimension
!
      fxi(1)=xi*(xi-1.d0)/2.d0
      fxi(2)=(1.d0-xi)*(1.d0+xi)
      fxi(3)=xi*(xi+1.d0)/2.d0
!
      fet(1)=et*(et-1.d0)/2.d0
      fet(2)=(1.d0-et)*(1.d0+et)
      fet(3)=et*(et+1.d0)/2.d0
!
      fze(1)=ze*(ze-1.d0)/2.d0
      fze(2)=(1.d0-ze)*(1.d0+ze)
      fze(3)=ze*(ze+1.d0)/2.d0
!
!     shape functions and their glocal derivatives
!
      omg=1.d0-xi
      omh=1.d0-et
      omr=1.d0-ze
      opg=1.d0+xi
      oph=1.d0+et
      opr=1.d0+ze
      tpgphpr=opg+oph+ze
      tmgphpr=omg+oph+ze
      tmgmhpr=omg+omh+ze
      tpgmhpr=opg+omh+ze
      tpgphmr=opg+oph-ze
      tmgphmr=omg+oph-ze
      tmgmhmr=omg+omh-ze
      tpgmhmr=opg+omh-ze
      omgopg=omg*opg/4.d0
      omhoph=omh*oph/4.d0
      omropr=omr*opr/4.d0
      omgmopg=(omg-opg)/4.d0
      omhmoph=(omh-oph)/4.d0
      omrmopr=(omr-opr)/4.d0
!
!     shape functions
!
      shp(4, 1)=-omg*omh*omr*tpgphpr/8.d0
      shp(4, 2)=-opg*omh*omr*tmgphpr/8.d0
      shp(4, 3)=-opg*oph*omr*tmgmhpr/8.d0
      shp(4, 4)=-omg*oph*omr*tpgmhpr/8.d0
      shp(4, 5)=-omg*omh*opr*tpgphmr/8.d0
      shp(4, 6)=-opg*omh*opr*tmgphmr/8.d0
      shp(4, 7)=-opg*oph*opr*tmgmhmr/8.d0
      shp(4, 8)=-omg*oph*opr*tpgmhmr/8.d0
      shp(4, 9)=omgopg*omh*omr
      shp(4,10)=omhoph*opg*omr
      shp(4,11)=omgopg*oph*omr
      shp(4,12)=omhoph*omg*omr
      shp(4,13)=omgopg*omh*opr
      shp(4,14)=omhoph*opg*opr
      shp(4,15)=omgopg*oph*opr
      shp(4,16)=omhoph*omg*opr
      shp(4,17)=omropr*omg*omh
      shp(4,18)=omropr*opg*omh
      shp(4,19)=omropr*opg*oph
      shp(4,20)=omropr*omg*oph
!
!     correction for the extra nodes in the middle of the faces
!
      do i=1,6
         if(konl(20+i).eq.konl(20)) then
!
!           no extra node in this face
!
            shp(4,20+i)=0.d0
         else
            shp(4,20+i)=fxi(jf(1,i))*fet(jf(2,i))*fze(jf(3,i))
            do j=1,4
               shp(4,ifaceq(j,i))=shp(4,ifaceq(j,i))+shp(4,20+i)/4.d0
            enddo
            do j=5,8
               shp(4,ifaceq(j,i))=shp(4,ifaceq(j,i))-shp(4,20+i)/2.d0
            enddo
         endif
      enddo
!
      if(iflag.eq.1) return
!
!     derivative of the shape functions in one dimension
!
      dfxi(1)=(2.d0*xi-1.d0)/2.d0
      dfxi(2)=-2.d0*xi
      dfxi(3)=(2.d0*xi+1.d0)/2.d0
!
      dfet(1)=(2.d0*et-1.d0)/2.d0
      dfet(2)=-2.d0*et
      dfet(3)=(2.d0*et+1.d0)/2.d0
!
      dfze(1)=(2.d0*ze-1.d0)/2.d0
      dfze(2)=-2.d0*ze
      dfze(3)=(2.d0*ze+1.d0)/2.d0
!
!     local derivatives of the shape functions: xi-derivative
!
      shpe(1, 1)=omh*omr*(tpgphpr-omg)/8.d0
      shpe(1, 2)=(opg-tmgphpr)*omh*omr/8.d0
      shpe(1, 3)=(opg-tmgmhpr)*oph*omr/8.d0
      shpe(1, 4)=oph*omr*(tpgmhpr-omg)/8.d0
      shpe(1, 5)=omh*opr*(tpgphmr-omg)/8.d0
      shpe(1, 6)=(opg-tmgphmr)*omh*opr/8.d0
      shpe(1, 7)=(opg-tmgmhmr)*oph*opr/8.d0
      shpe(1, 8)=oph*opr*(tpgmhmr-omg)/8.d0
      shpe(1, 9)=omgmopg*omh*omr
      shpe(1,10)=omhoph*omr
      shpe(1,11)=omgmopg*oph*omr
      shpe(1,12)=-omhoph*omr
      shpe(1,13)=omgmopg*omh*opr
      shpe(1,14)=omhoph*opr
      shpe(1,15)=omgmopg*oph*opr
      shpe(1,16)=-omhoph*opr
      shpe(1,17)=-omropr*omh
      shpe(1,18)=omropr*omh
      shpe(1,19)=omropr*oph
      shpe(1,20)=-omropr*oph
!
!     correction for the extra nodes in the middle of the faces
!
      do i=1,6
         if(konl(20+i).eq.konl(20)) then
!
!           no extra node in this face
!
            shpe(1,20+i)=0.d0
         else
            shpe(1,20+i)=dfxi(jf(1,i))*fet(jf(2,i))*fze(jf(3,i))
            do j=1,4
               shpe(1,ifaceq(j,i))=shpe(1,ifaceq(j,i))+shpe(1,20+i)/4.d0
            enddo
            do j=5,8
               shpe(1,ifaceq(j,i))=shpe(1,ifaceq(j,i))-shpe(1,20+i)/2.d0
            enddo
         endif
      enddo
!
!     local derivatives of the shape functions: eta-derivative
!
      shpe(2, 1)=omg*omr*(tpgphpr-omh)/8.d0
      shpe(2, 2)=opg*omr*(tmgphpr-omh)/8.d0
      shpe(2, 3)=opg*(oph-tmgmhpr)*omr/8.d0
      shpe(2, 4)=omg*(oph-tpgmhpr)*omr/8.d0
      shpe(2, 5)=omg*opr*(tpgphmr-omh)/8.d0
      shpe(2, 6)=opg*opr*(tmgphmr-omh)/8.d0
      shpe(2, 7)=opg*(oph-tmgmhmr)*opr/8.d0
      shpe(2, 8)=omg*(oph-tpgmhmr)*opr/8.d0
      shpe(2, 9)=-omgopg*omr
      shpe(2,10)=omhmoph*opg*omr
      shpe(2,11)=omgopg*omr
      shpe(2,12)=omhmoph*omg*omr
      shpe(2,13)=-omgopg*opr
      shpe(2,14)=omhmoph*opg*opr
      shpe(2,15)=omgopg*opr
      shpe(2,16)=omhmoph*omg*opr
      shpe(2,17)=-omropr*omg
      shpe(2,18)=-omropr*opg
      shpe(2,19)=omropr*opg
      shpe(2,20)=omropr*omg
!
!     correction for the extra nodes in the middle of the faces
!
      do i=1,6
         if(konl(20+i).eq.konl(20)) then
!
!           no extra node in this face
!
            shpe(2,20+i)=0.d0
         else
            shpe(2,20+i)=fxi(jf(1,i))*dfet(jf(2,i))*fze(jf(3,i))
            do j=1,4
               shpe(2,ifaceq(j,i))=shpe(2,ifaceq(j,i))+shpe(2,20+i)/4.d0
            enddo
            do j=5,8
               shpe(2,ifaceq(j,i))=shpe(2,ifaceq(j,i))-shpe(2,20+i)/2.d0
            enddo
         endif
      enddo
!
!     local derivatives of the shape functions: zeta-derivative
!
      shpe(3, 1)=omg*omh*(tpgphpr-omr)/8.d0
      shpe(3, 2)=opg*omh*(tmgphpr-omr)/8.d0
      shpe(3, 3)=opg*oph*(tmgmhpr-omr)/8.d0
      shpe(3, 4)=omg*oph*(tpgmhpr-omr)/8.d0
      shpe(3, 5)=omg*omh*(opr-tpgphmr)/8.d0
      shpe(3, 6)=opg*omh*(opr-tmgphmr)/8.d0
      shpe(3, 7)=opg*oph*(opr-tmgmhmr)/8.d0
      shpe(3, 8)=omg*oph*(opr-tpgmhmr)/8.d0
      shpe(3, 9)=-omgopg*omh
      shpe(3,10)=-omhoph*opg
      shpe(3,11)=-omgopg*oph
      shpe(3,12)=-omhoph*omg
      shpe(3,13)=omgopg*omh
      shpe(3,14)=omhoph*opg
      shpe(3,15)=omgopg*oph
      shpe(3,16)=omhoph*omg
      shpe(3,17)=omrmopr*omg*omh
      shpe(3,18)=omrmopr*opg*omh
      shpe(3,19)=omrmopr*opg*oph
      shpe(3,20)=omrmopr*omg*oph
!
!     correction for the extra nodes in the middle of the faces
!
      do i=1,6
         if(konl(20+i).eq.konl(20)) then
!
!           no extra node in this face
!
            shpe(3,20+i)=0.d0
         else
            shpe(3,20+i)=fxi(jf(1,i))*fet(jf(2,i))*dfze(jf(3,i))
            do j=1,4
               shpe(3,ifaceq(j,i))=shpe(3,ifaceq(j,i))+shpe(3,20+i)/4.d0
            enddo
            do j=5,8
               shpe(3,ifaceq(j,i))=shpe(3,ifaceq(j,i))-shpe(3,20+i)/2.d0
            enddo
         endif
      enddo
!
!     computation of the local derivative of the global coordinates
!     (xs)
!
      do i=1,3
        do j=1,3
          xs(i,j)=0.d0
          do k=1,26
            xs(i,j)=xs(i,j)+xl(i,k)*shpe(j,k)
          enddo
        enddo
      enddo
!
!     computation of the jacobian determinant
!
      dd1=xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2)
      dd2=xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3)
      dd3=xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1)
      xsj=xs(1,1)*dd1+xs(1,2)*dd2+xs(1,3)*dd3
!
      if(iflag.eq.2) return
!
      dd=1.d0/xsj
!
!     computation of the global derivative of the local coordinates
!     (xsi) (inversion of xs)
!
      xsi(1,1)=dd1*dd
      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd
      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd
      xsi(2,1)=dd2*dd
      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd
      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd
      xsi(3,1)=dd3*dd
      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd
      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd
!
!     computation of the global derivatives of the shape functions
!
      do k=1,26
        do j=1,3
          shp(j,k)=shpe(1,k)*xsi(1,j)+shpe(2,k)*xsi(2,j)
     &          +shpe(3,k)*xsi(3,j)
        enddo
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
You may recognize Amy Zimmerman but you may be confused because she has a twin sister named Users/SharonZimmerman Sharon (who affectionately calls Amy hobag).

Amy graduated from UC Davis in December 2006. She was the copy chief for The California Aggie The Aggie and thinks good writing is sexy.

Amy now lives in San Francisco and is a contributor to wiki:sf SF Wiki

20051018 16:26:13 nbsp miss amy you apparently know who this is already...I do want to say that I did have fun working with someone with whom I share the same love of writing.  Anyway I do want to apologize for the comments about our work.  No its not that bad at all.  In fact I love the people we work with, but some of the comments posted upon the site were merely aimed at our beloved employees.  these people I have no respect for. but until next time...Chris Users/ChrisJernigan

20051022 20:42:08 nbsp Amy,Amy,Amy. Well our time grows short in the lovely little city of Davis...And truth is, my time spent talking and bullshitting with you in the kitchen, must and has come to an end.  I love writing. Point blank. And you can guarantee that Ill do my best to keep in touch with one of the most honest, one of the realest little women that Ive met along this trip.  I will write to you amy...you can bet on it.  And if you ever need an honest(maybe blunt) opinion, Well Ill send you my email.  Until then darlin.  have a good one po Users/ChrisJernigan

20060101 11:21:07 nbsp Just a random note, but you and Kristy have a link to the same DeviantArt account on the Photography page. Not sure if this was intentional. Users/MichaelGiardina

20060317 14:07:47 nbsp good call ill replicate the McDonalds experiment with InNOut coke. Users/CraigBrozinsky

20060321 21:26:59 nbsp All she cares about in life is Sushi  that goes the same for her sister. Users/RobRoy

20060708 17:01:28 nbsp Hi, amy...you dont know me, but Ill be moving to Davis in September to start grad school.  Ive been looking for good resources for rat care, breeders, and enthusiasts.  I see you have a few pretty girls of your own!  Im currently ratsitting for an ex, two lovely himalayans, Rosa and Lenna.  Ill be starting fresh come september.  Any advice? Users/PatrickDragon

20060802 20:45:59 nbsp Hi Amy, could you add Mike Steinman and Brian Ang to the fundraiser donor list. We donated during the fundraiser. Users/BrianAng

20071221 00:32:43 nbsp Long time no see.  Happy Hanukkah (belatedly)  ~ David Users/gdnystrom

20091101 20:55:03 nbsp Business have massive water filters that come thru the line before it hits any part where its used or consumed, Davis water shouldnt affect soda from a fountain, too much... Users/StevenDaubert


!______________________________________________________________________________
      subroutine post_axial(icase,iq1,iq2,jprmi,jprma,jdf)

!..............................................................................
!     write matrix elements for Hill-Wheeler-Griffin solver to file lou
!..............................................................................
      USE VAPHFB_PAR
      implicit real*8 (a-h,o-z)

  401 format (10i5)
  402 format (8e15.8)

      do jproj = jprmi,jprma,jdf
         if(jproj.eq.1) cycle
          write (lou,401) jproj
          write (lou,402) Const%beta2t_mesh(iq1),Const%gamma2t_mesh(iq1),&
     &                    Const%beta2t_mesh(iq2),Const%gamma2t_mesh(iq2)
!    .........................................................
          k1 = 0
          k2 = 0
          write (lou,401) k1,k2
          write (lou,402) Kernel%njkk(jproj,k1,k2), &
     &                    Kernel%hjkk(jproj,k1,k2)
          write (lou,402) Kernel%nn(jproj,k1,k2), &
     &                    Kernel%zz(jproj,k1,k2)
! .................... newly added for J^2
          write (lou,402) Kernel%J2(jproj,k1,k2)
          write (lou,402) Kernel%q0p(jproj,k1,k2),Kernel%q0pc(jproj,k1,k2)
          write (lou,402) Kernel%qpred(jproj,k1,k2,0),Kernel%qcpred(jproj,k1,k2,0)
          write (lou,402) Kernel%qpred(jproj,k1,k2,1),Kernel%qcpred(jproj,k1,k2,1)
          write (lou,402) Kernel%qpred(jproj,k1,k2,2),Kernel%qcpred(jproj,k1,k2,2)

!     ......................................
  99  enddo ! jproj
      return
      end

!______________________________________________________________________________
      subroutine post_general(icase,iq1,iq2,jprmi,jprma,jdf)

!.............................................................................
!     write matrix elements for Hill-Wheeler-Griffin solver to file lou      .
!.............................................................................
      USE VAPHFB_PAR
      implicit real*8 (a-h,o-z)

  401 format (10i5)
  402 format (8e15.8)

      do jproj = jprmi,jprma,jdf
         if(jproj.eq.1) cycle
!       .................................... write matrix elements to tape
        write (lou,401) jproj
        write (lou,402) Const%beta2t_mesh(iq1),Const%gamma2t_mesh(iq1), &
     &                  Const%beta2t_mesh(iq2),Const%gamma2t_mesh(iq2)
!    .........................................................
        if(iv(jproj).gt.zero) then         ! for even spin
           k2max  =  jproj
        endif
        if(iv(jproj).lt.zero) then       ! for odd spin
           k2max  =  abs(jproj-1)
        endif
!    .........................................................
         do k2  =  -k2max, k2max, 1
         do k1  =  -k2max, k2max, 1
               write (lou,401) k1,k2
               write (lou,402) Kernel%njkk(jproj,k1,k2),                    &
     &                         Kernel%hjkk(jproj,k1,k2) !/Kernel%njkk(jproj,k1,k2)
               write (lou,402) Kernel%nn(jproj,k1,k2), &
     &                         Kernel%zz(jproj,k1,k2)
! .................... newly added for J^2
               write (lou,402) Kernel%J2(jproj,k1,k2)
               write (lou,402) Kernel%qpred(jproj,k1,k2,0),Kernel%qcpred(jproj,k1,k2,0)
               write (lou,402) Kernel%q0p(jproj,k1,k2),Kernel%q0pc(jproj,k1,k2)
              ! added for test
             ! write (lou,402) Kernel%qpred(jproj,k1,k2,2),Kernel%qcpred(jproj,k1,k2,2)
              !write (*,'(i4,4f12.8)') jproj,Kernel%qpred(jproj,k1,k2,0),Kernel%qcpred(jproj,k1,k2,0)
          enddo ! k2
         enddo ! k1

         !  <J2,K2|| Q2 || J1,K1>, where J2=J1,J1+1,J1+2
         do k2  =  -k2max, k2max, 1
              do k1  =  -k2max-1, k2max+1, 1
                 write (lou,402) Kernel%qpred(jproj,k1,k2,1),Kernel%qcpred(jproj,k1,k2,1)
              enddo ! k1
              do k1  =  -k2max-2, k2max+2, 1
                 write (lou,402) Kernel%qpred(jproj,k1,k2,2),Kernel%qcpred(jproj,k1,k2,2)
              enddo ! k1
          enddo ! k2

!     ......................................
  99  enddo ! jproj
      return
      end

!.......................................................................

      subroutine ReadCoordinators(lpr)
      use vaphfb_par
      implicit none
      logical lpr
      integer NQ,iq
!    .......... read coordinates
         open(15,file='betgam.dat',status='old')
         read(15,*)
         read(15,*)
         read(15,*)
         read(15,*)
!     ............
         read(15,*) NQ
         if(.NOT. ALLOCATED(Const%beta2t_mesh))   ALLOCATE(Const%beta2t_mesh(1:NQ))
         if(.NOT. ALLOCATED(Const%gamma2t_mesh))  ALLOCATE(Const%gamma2t_mesh(1:NQ))
         if(.NOT. ALLOCATED(Const%P00_mesh))      ALLOCATE(Const%P00_mesh(1:NQ))
         if(.NOT. ALLOCATED(Const%hw_mesh))      ALLOCATE(Const%hw_mesh(1:NQ))
         do iq=1,NQ
            read(15,*) Const%beta2t_mesh(iq),Const%gamma2t_mesh(iq),&
                       Const%P00_mesh(iq), Const%hw_mesh(iq)
            call Filename4wfs(Const%beta2t_mesh(iq),Const%gamma2t_mesh(iq),&
                              Const%P00_mesh(iq), Const%hw_mesh(iq))
            File%cwf(iq) = File%wf    ! name of files for wfs
            if(lpr) print *,File%cwf(iq)
         enddo !
          Input%NGCM=NQ
            if(lpr) print *,'Input%NGCM=',Input%NGCM
         if(Input%nq0f.gt.NQ .or. Input%nq1f.gt.NQ) &
     &   stop ' Error: Input%nq0f or Input%nq1f are larger than the total &
         &number of configurations!'
      end

!________________________________________________________________________


text:
subroutine sld_turnofGeneral(a)
   ! DESCRIPTION
   !    This routine closes the run for the elastic solid prob
   !-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_Solids
   implicit none
   class(SolidsProblem) :: a
   integer(ip)          :: ndime,npoin,ncomp
   integer(ip)          :: ielem,nelem,pgaus,pnode
   integer(ip)          :: nboun,sz

   !Unknowns
   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNboun(nboun)
   call a%Mesh%GetNelem(nelem)
   ncomp = a%ncomp

   sz = (ndime*(ndime+1))/2

   !Write tail for formatted files
   if (a%MPIrank == a%MPIroot ) then
       if (a%kfl_isSlave .eqv. .false.) then
           write(a%lun_outph,'(//,a,/)') '     * * * END OF ELASTIC SOLID RUN * * *'
           write(a%lun_solve,'(//,a,/)') '     * * * END OF ELASTIC SOLID RUN * * *'
       endif
   endif

   call a%Memor%dealloc(ndime,npoin,a%ncomp,a%disp,     'disp','sld_turnof')
   if(a%kfl_NodalStress) call a%Memor%dealloc(sz   ,npoin  ,a%stress,'stress','sld_turnof')
   if(a%kfl_NodalStrain) call a%Memor%dealloc(sz   ,npoin  ,a%strain,'strain','sld_turnof')

   call a%Memor%dealloc(ndime,npoin,a%btraction_nodal,'Bnodaltraction','sld_turnof')

   if(a%kfl_docoupconv) then
       call a%Memor%dealloc(ndime,npoin,a%disp_cp, 'disp_cp','sld_turnof')
   endif

   if (a%kfl_timei==1) then
       call a%Memor%dealloc(ndime,npoin,2,a%veloc,'veloc','sld_turnof')
       call a%Memor%dealloc(ndime,npoin,2,a%accel,'accel','sld_turnof')
   endif

   call a%Memor%dealloc(npoin,a%kfl_fixrs,'kfl_fixrs','sld_turnof')
   call a%Memor%dealloc(nboun,a%kfl_bours,'kfl_bours','sld_turnof')

   if (a%kfl_sourcesPresent) then
       call a%Memor%dealloc(a%pwSourceDim,a%pwSourceSize,a%pwSource,'pwSource','sld_turnof')
       call a%Memor%dealloc(a%pwSourceSize,a%pwSourceId,'pwSourceId','sld_turnof')
   endif

   if(a%kfl_isSlave .eqv. .false.) then
      if (a%nptra/=0) then
         if (a%MPIrank == a%MPIroot) then
            call iofile(two,a%lun_trap,a%fil_trap,'SLD TRACKING OF POINTS')
         endif
      endif
   endif

   do ielem=1,nelem
       call a%Mesh%GetElemArraySize(ielem,pnode,pgaus)
       deallocate(a%extForce(ielem)%a)
       if(a%kfl_printStress) deallocate(a%stress_g(ielem)%a)
       if(a%kfl_printStrain) deallocate(a%strain_g(ielem)%a)
       deallocate(a%btraction(ielem)%a)
   end do

   call a%Memor%dealloc(nelem,a%extForce,'extForce','sld_turnof')
   if(a%kfl_printStress) call a%Memor%dealloc(nelem,a%stress_g,'stress_g','sld_turnof')
   if(a%kfl_printStrain) call a%Memor%dealloc(nelem,a%strain_g,'strain_g','sld_turnof')
   call a%Memor%dealloc(nelem,a%btraction,'btraction','sld_turnof')

   call a%SolidSpecificTurnof

end subroutine sld_turnofGeneral

subroutine sld_turnofSpecific(a)
   ! DESCRIPTION
   !    This routine closes the run for the elastic solid prob
   !-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_Solids
   implicit none
   class(SolidsProblem) :: a
   integer(ip)          :: ndime,ielem,nelem,npoin
   integer(ip)          :: sz,pgaus,pnode

   !Unknowns
   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNelem(nelem)

   sz = (ndime*(ndime+1))/2

   call a%Memor%dealloc(sz,sz,a%c_elastic,   'c_elastic','sld_turnof')

   if(a%kfl_printNodalSigma) call a%Memor%dealloc(sz   ,npoin,1,a%sigma,'sigma','sld_memall')
   if(a%kfl_printNodalPress) call a%Memor%dealloc(      npoin,1,a%press,'press','sld_memall')

   do ielem=1,nelem
       call a%Mesh%GetElemArraySize(ielem,pnode,pgaus)
       if (a%kfl_printGaussSigma) then
           deallocate(a%sigma_g(ielem)%a)
       endif

       if (a%kfl_printGaussPress) then
           deallocate(a%press_g(ielem)%a)
       endif

       !Memory error here, beware
       deallocate(a%pStrain(ielem)%a)

       if(a%kfl_printPrincipalStresses) then
           deallocate(a%printPStrain(ielem)%a)
       end if

   end do

   if(a%kfl_printPrincipalStresses) then
       call a%Memor%dealloc(nelem,a%printPStrain,'printPStrain','sld_turnof')
   end if
   call a%Memor%dealloc(nelem,a%pStrain,'pStrain','sld_turnof')
   if (a%kfl_printGaussSigma) call a%Memor%dealloc(nelem,a%sigma_g,'sigma_g','sld_turnof')
   if (a%kfl_printGaussPress) call a%Memor%dealloc(nelem,a%press_g,'press_g','sld_turnof')

end subroutine sld_turnofSpecific

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
! ###################################################################
! Copyright (c) 2013-2022, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:math.f90
!--------------------------------------------------------------------------
!
! MODULE: math
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief collection of mathematical/numerical routines that don't fit anywhere else
!
!> @date 10/13/98 MDG 1.0 original
!> @date 05/19/01 MDG 2.0 f90
!> @date 11/27/01 MDG 2.1 added kind support
!> @date 03/19/13 MDG 3.0 updated all routines
!> @date 11/13/13 MDG 4.0 added MatrixExponential routine
!> @date 11/23/15 MDG 4.1 moved several routines from other mods into this one
!> @date 10/24/17 MDG 4.2 added infty()/inftyd() functions to return the IEEE infinity value
!> @date 08/23/19 MDG 4.3 removed spaces around "kind" statements to facilitate f90wrap python wrapper generation
!> @date 10/04/19 MDG 4.4 adds vecnorm to replace non-standard NORM2 calls  (F2003 compliance)
!> @date 10/04/19 MDG 4.5 adds nan() function, returning a single or double precision IEEE NaN value
!> @date 11/01/19 MDG 4.6 adds Jaccard_Distance routine (moved from Indexingmod)
!--------------------------------------------------------------------------
! ###################################################################
!

module mod_math
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! collection of mathematical/numerical routines that don't fit anywhere else (no classes, just routines)

use mod_kinds
use mod_global

! public :: mInvert, cross3, infty, inftyd, nan, nan_d,

interface mInvert
        module procedure mInvert
        module procedure mInvert_d
end interface

interface cross3
        module procedure cross3
        module procedure cross3_d
end interface

interface vecnorm
        module procedure vecnorm
        module procedure vecnorm_d
        module procedure vecnorm2
        module procedure vecnorm2_d
end interface

contains

!--------------------------------------------------------------------------
recursive function vecnorm(vec) result(veclen)
!DEC$ ATTRIBUTES DLLEXPORT :: vecnorm
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the single precision length of a 1D vector

real(kind=sgl),INTENT(IN)        :: vec(:)
real(kind=sgl)                   :: veclen

integer(kind=irg)                :: sz(1)

sz = size(vec)

veclen = sqrt(sum(vec(1:sz(1))*vec(1:sz(1))))

end function vecnorm

!--------------------------------------------------------------------------
recursive function vecnorm_d(vec) result(veclen)
  !DEC$ ATTRIBUTES DLLEXPORT :: vecnorm_d
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the double precision length of a 1D vector

real(kind=dbl),INTENT(IN)       :: vec(:)
real(kind=dbl)                  :: veclen

integer(kind=irg)               :: sz(1)

sz = size(vec)

veclen = sqrt(sum(vec(1:sz(1))*vec(1:sz(1))))

end function vecnorm_d

!--------------------------------------------------------------------------
recursive function vecnorm2(vec) result(veclen)
!DEC$ ATTRIBUTES DLLEXPORT :: vecnorm2
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the single precision length of a 2D array

real(kind=sgl),INTENT(IN)        :: vec(:,:)
real(kind=sgl)                   :: veclen

integer(kind=irg)                :: sz(2)

sz = size(vec)

veclen = sqrt(sum(vec(1:sz(1),1:sz(2))*vec(1:sz(1),1:sz(2))))

end function vecnorm2

!--------------------------------------------------------------------------
recursive function vecnorm2_d(vec) result(veclen)
!DEC$ ATTRIBUTES DLLEXPORT :: vecnorm2_d
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the double precision length of a 2D array

real(kind=dbl),INTENT(IN)        :: vec(:,:)
real(kind=dbl)                   :: veclen

integer(kind=irg)                :: sz(2)

sz = size(vec)

veclen = sqrt(sum(vec(1:sz(1),1:sz(2))*vec(1:sz(1),1:sz(2))))

end function vecnorm2_d

!--------------------------------------------------------------------------
recursive function infty() result(infinity)
!DEC$ ATTRIBUTES DLLEXPORT :: infty
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the single precision IEEE value for infinity

real(kind=sgl)      :: infinity
real(kind=sgl)      :: big

big = HUGE(1.0)
infinity = big + HUGE(1.0)

end function infty

!--------------------------------------------------------------------------
recursive function inftyd() result(infinity)
!DEC$ ATTRIBUTES DLLEXPORT :: inftyd
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the double precision IEEE value for infinity

real(kind=dbl)      :: infinity
real(kind=dbl)      :: big

big = HUGE(1.D0)
infinity = big + HUGE(1.D0)

end function inftyd

!--------------------------------------------------------------------------
recursive function nan() result(x)
!DEC$ ATTRIBUTES DLLEXPORT :: nan
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the single precision IEEE value for nan

 use, intrinsic :: iso_fortran_env
 use, intrinsic :: ieee_arithmetic

 IMPLICIT NONE

real(kind=sgl)        :: x

x = ieee_value(x, ieee_quiet_nan)

end function nan

!--------------------------------------------------------------------------
recursive function nan_d() result(x)
!DEC$ ATTRIBUTES DLLEXPORT :: nan_d
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! return the double precision IEEE value for nan

 use, intrinsic :: iso_fortran_env
 use, intrinsic :: ieee_arithmetic

 IMPLICIT NONE

real(kind=dbl)        :: x

x = ieee_value(x, ieee_quiet_nan)

end function nan_d

!--------------------------------------------------------------------------
!
! FUNCTION: cross3
!
!> @author Saransh Singh, Carnegie Mellon University
!
!> @brief  cross product of two 3D vector in the order of input
!
!
!> @param u input vector 1
!> @param v input vector 2
!
!> @date 03/03/16   SS 1.0 original
!> @date 12/01/16  MDG 1.1 split in single and double precision versions
!--------------------------------------------------------------------------
recursive function cross3(u, v) result(res)
!DEC$ ATTRIBUTES DLLEXPORT :: cross3

IMPLICIT NONE

real(kind=sgl),INTENT(IN)      :: u(3)
real(kind=sgl),INTENT(IN)      :: v(3)
real(kind=sgl)                 :: res(3)

res(1) = u(2)*v(3) - u(3)*v(2)
res(2) = u(3)*v(1) - u(1)*v(3)
res(3) = u(1)*v(2) - u(2)*v(1)

end function cross3

!--------------------------------------------------------------------------
!
! FUNCTION: cross3_d
!
!> @author Saransh Singh, Carnegie Mellon University
!
!> @brief  cross product of two 3D vector in the order of input (double precision)
!
!
!> @param u input vector 1
!> @param v input vector 2
!
!> @date 03/03/16   SS 1.0 original
!--------------------------------------------------------------------------
recursive function cross3_d(u, v) result(res)
!DEC$ ATTRIBUTES DLLEXPORT :: cross3_d

IMPLICIT NONE

real(kind=dbl),INTENT(IN)      :: u(3)
real(kind=dbl),INTENT(IN)      :: v(3)
real(kind=dbl)                 :: res(3)

res(1) = u(2)*v(3) - u(3)*v(2)
res(2) = u(3)*v(1) - u(1)*v(3)
res(3) = u(1)*v(2) - u(2)*v(1)

end function cross3_d

!--------------------------------------------------------------------------
recursive subroutine mInvert_d(a,b,uni)
!DEC$ ATTRIBUTES DLLEXPORT :: mInvert_d
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! Invert a 3x3 matrix; if unitary, simply transpose

use mod_io

IMPLICIT NONE

real(kind=dbl),INTENT(IN)               :: a(3,3)
 !! input matrix
real(kind=dbl),INTENT(OUT)              :: b(3,3)
 !! output matrix
logical,INTENT(IN)                      :: uni
 !! unitary logical

type(IO_T)                              :: Message
real(kind=dbl)                          :: d
integer(kind=irg)                       :: i, j

! it is a regular (non-unitary) matrix
 if (.not.uni) then
  d = a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+ &
         a(1,3)*a(2,1)*a(3,2)-a(1,3)*a(2,2)*a(3,1)- &
         a(1,2)*a(2,1)*a(3,3)-a(1,1)*a(2,3)*a(3,2)
  if (d.ne.0.D0) then
   b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2)
   b(1,2)=a(1,3)*a(3,2)-a(1,2)*a(3,3)
   b(1,3)=a(1,2)*a(2,3)-a(1,3)*a(2,2)
   b(2,1)=a(2,3)*a(3,1)-a(2,1)*a(3,3)
   b(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1)
   b(2,3)=a(1,3)*a(2,1)-a(1,1)*a(2,3)
   b(3,1)=a(2,1)*a(3,2)-a(2,2)*a(3,1)
   b(3,2)=a(1,2)*a(3,1)-a(1,1)*a(3,2)
   b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1)
   b = b/d
  else
    do i=1,3
      write (*,*) (a(i,j),j=1,3)
    end do
!  call Message%printError('mInvert','matrix has zero determinant')
   call Message%printMessage('mInvert: matrix has zero determinant')
   b = a
  end if
 else
! it is a unitary matrix, so simply get the transpose
  b = transpose(a)
 endif

end subroutine mInvert_d

!--------------------------------------------------------------------------
recursive subroutine mInvert(a,b,uni)
!DEC$ ATTRIBUTES DLLEXPORT :: mInvert
  !! author: MDG
  !! version: 1.0
  !! date: 01/17/20
  !!
  !! Invert a single precision 3x3 matrix; if unitary, simply transpose

use mod_io

IMPLICIT NONE

real(kind=sgl),INTENT(IN)               :: a(3,3)
 !! input matrix
real(kind=sgl),INTENT(OUT)              :: b(3,3)
 !! output matrix
logical,INTENT(IN)                      :: uni
 !! unitary logical

type(IO_T)                              :: Message
real(kind=sgl)                          :: d                    !< auxiliary variable
integer(kind=irg)                       :: i, j
! it is a regular (non-unitary) matrix
 if (.not.uni) then
  d = a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+ &
         a(1,3)*a(2,1)*a(3,2)-a(1,3)*a(2,2)*a(3,1)- &
         a(1,2)*a(2,1)*a(3,3)-a(1,1)*a(2,3)*a(3,2)
  if (d.ne.0.0) then
   b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2)
   b(1,2)=a(1,3)*a(3,2)-a(1,2)*a(3,3)
   b(1,3)=a(1,2)*a(2,3)-a(1,3)*a(2,2)
   b(2,1)=a(2,3)*a(3,1)-a(2,1)*a(3,3)
   b(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1)
   b(2,3)=a(1,3)*a(2,1)-a(1,1)*a(2,3)
   b(3,1)=a(2,1)*a(3,2)-a(2,2)*a(3,1)
   b(3,2)=a(1,2)*a(3,1)-a(1,1)*a(3,2)
   b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1)
   b = b/d
  else
    do i=1,3
      write (*,*) (a(i,j),j=1,3)
    end do
!  call Message%printError('mInvert','matrix has zero determinant')
   call Message%printMessage('mInvert: matrix has zero determinant')
   b = a
  end if
 else
! it is a unitary matrix, so simply get the transpose
  b = transpose(a)
 endif

end subroutine mInvert

end module mod_math

max_stars_repo_name: marcdegraef/sampleRFZ

max_issues_repo_name: marcdegraef/sampleRFZ

max_forks_repo_name: marcdegraef/sampleRFZ


text:
!--------------------------------------------------------------------------------
!
!       GEOFRACTAL version 1.0
!
!--------------------------------------------------------------------------------
!
! GEOFRACTAL v1.0 computes the geometrical cross section of randomly orientated
! fractal dust aggregates based on a statistical distribution model of monomer
! particles developed in Tazaki (in prep.).
!
! If you wish to publish a paper that contains results obtained by this code,
! please acknowledge this code and cite relevant papers.
!
!
! Disclaimer:
! I reject all responsibility for the use of this code. Although it has been
! tested, the code may still contain a bug. I am not responsible for any damages
! caused by the use of the code. If you find a bug, please let me know.
!
!                                                   Ryo Tazaki (2021/Jan/1)
!                                                   email: r.tazaki -at- uva.nl
!
!--------------------------------------------------------------------------------
! INPUT PARAMETERS
!--------------------------------------------------------------------------------
!
! Df            : fractal dimension
! PN            : number of monomers
! k0            : fractal prefactor defined by PN = k0*(Rg/R0)^Df
!
! iqapp         : Select a method for solving radial and angular integration
!                 when computing the mean overlapping efficiency.
!                --------------------------------------------------------
!                          |    radial  (x)        |  angular (u)        |
!                ---------------------------------------------------------
!                iqapp = 1 |   numerical           |  numerical          |
!                iqapp = 3 |   numerical  ( D<=2 ) |  analytical         |
!                          |   analytical ( D> 2 ) |                     |
!                --------------------------------------------------------
!
! iqcon         : Switch to specify the numerical factor A.
!                 iqcon  = 1 : A = 1
!                 iqcon  = 2 : A is chosen to connect two regimes.
!
! iqcor         : Switch to specify a cut-off model of the correlation function
!                 iqcon  = 1 : The Gaussian cut-off function
!                 iqcon  = 2 : The exponential cut-off function
!                 iqcon  = 3 : The fractal dimension cut-off function
!
! Prefered set of options is iqcor=3, iqapp=3 and iqcon=2.
! Note that iqcor=1 is also possible, but avoid using iqcor=2.
!
!--------------------------------------------------------------------------------
! OUTPUT PARAMETER
!--------------------------------------------------------------------------------
!
! G  : The geocross section of the aggregate normalized by a factor PN*pi*R0^2.
!      Thus, "G" is non-dimensional quantity.
!
!--------------------------------------------------------------------------------

!module types
!implicit none
!integer, parameter      :: dp = selected_real_kind(P=15)
!end module types
!
!module const
!use types
!implicit none
!real(kind=dp),parameter :: pi = 3.1415926535897932384_dp
!end module const

subroutine geofrac(iqapp,iqcon,iqcor,PN,k0,df,G)
use types
implicit none
integer                 :: iqapp,iqcon,iqcor
real(kind=dp)           :: PN,k0,df,GS,G
real(kind=dp)           :: pnth,gth,sigmath,A,sigma
logical,parameter::debug=.false.

!--------------------------------------------------------------------------------
! some checks
!--------------------------------------------------------------------------------
if(PN .le. 0.9999_dp) then
        print *, 'number of monomers =',PN
        print *, 'error: The number of monomers should be larger than 1.'
        print *, 'stop'
        stop
endif
if(df .lt. 0.9999_dp .or. df .gt. 3.00001) then
        print *, 'fractal dimension df =',df
        print *, 'error: Fractal dimension df should be 1<=df<=3.'
        print *, 'stop'
        stop
endif
if(iqcon .ne. 1 .and. iqcon .ne. 2) then
        print *, 'error: iqcon should be either 1 or 2.'
        print *, 'stop'
        stop
endif
if(iqcor .ne. 1 .and. iqcor .ne. 2 .and. iqcor .ne. 3) then
        print *, 'error: iqcor should be either 1, 2, or 3.'
        print *, 'stop'
        stop
endif
if(iqapp .ne. 1 .and. iqapp .ne. 2 .and. iqapp .ne. 3) then
        print *, 'error: iqapp should be either 1 or 3.'
        print *, 'stop'
        stop
endif
if(iqapp .eq. 2 .and. (debug .eqv. .false.)) then
        print *, 'error: iqapp = 2 is only allowed for a debug mode.'
        print *, 'stop'
        stop
endif

!--------------------------------------------------------------------------------
! start calculations
!--------------------------------------------------------------------------------
A = 1.0_dp
if(iqcon .eq. 2) then
        PNTH = 11.0_dp*df-8.5_dp
        PNTH = min(PNTH,8.0_dp)
        if(PN .lt. PNTH) then
                G=GS(PN)
                return
        endif
        call mean_overlap_efficiency(iqapp,iqcor,k0,df,PNTH,sigmath)
        A = (1.0_dp+(PNTH-1.0_dp)*sigmath)*GS(PNTH)
endif
call mean_overlap_efficiency(iqapp,iqcor,k0,df,PN,sigma)
G = A / (1.0_dp + (PN-1.0_dp)*sigma)

return
end subroutine geofrac

!--------------------------------------------------------------------------------
!
! Minato et al. (2006)'s fitting formula for N < 16
!
!--------------------------------------------------------------------------------
function GS(PN)
use types
implicit none
real(kind=dp)::GS,PN
        GS = 12.5_dp*PN**(-0.315_dp)*exp(-2.53_dp/PN**0.0920_dp)
return
end function GS

!--------------------------------------------------------------------------------
!
! This subroutine calculations the mean overlapping efficiency: sigma
! (In Tazaki in prep., this quantity corresponds to \tilde{sigma} )
!
! For iqcor=1 [Gaussian cut-off model]
!
!               xmin       /ln(xmax)
! sigma  = --------------  |   dlnx frad(x)*S(rho),
!          16*Gamma(df/2)  /ln(xmin)
!
!       xmin = df*(k0/PN)^(2/df); a=(df-2)/2; rho = sqrt(x/xmin)
!
! For iqcor=2 [Exponential cut-off model]
!
!              xmin^2      /ln(xmax)
! sigma  = --------------  |   dlnx frad(x)*S(rho),
!           16*Gamma(df)   /ln(xmin)
!
!       xmin = 2*sqrt(df(df+1)/2)*(k0/PN)^(1.0/df); a=df-2; rho = x/xmin
!
! For iqcor=3 [Fractal dimension cut-off model]
!
!          xmin^{2/df} /ln(xmax)
! sigma  = ----------  |  dlnx frad(x)*S(rho),
!             16      /ln(xmin)
!
!       xmin = 2^(df-1)*k0/PN; a=(df-2)/df; rho = (x/xmin)^(1/df)
!
! where df is fractal dimension, R0 is the monomer radius, Rg is
! the radius of gyration, frad(x) is the radial integrand function,
! and S(rho) is the angular integral function.
!
! frad(x)dlnx = x^{a}*exp(-x) dlnx = x^{a-1}exp(-x) dx
!
!           16*rho^2  / (1/rho)
! S(rho)  = --------  |  du fang(rho,u),
!              pi     /0
!
! fang(rho,u) = [Arcsin{sqrt(1-rho^2u^2)} - rho*u*sqrt{1-rho^2u^2}]*u/sqrt(1-u^2)
!
! Note that the angular integral function S(rho) ≈1 when rho >> 1.
!
!--------------------------------------------------------------------------------
subroutine mean_overlap_efficiency(iqapp,iqcor,k0,df,PN,sigma)
use types; use const
implicit none
integer,parameter       :: nx = 1000
real(kind=dp),parameter :: xmax = 25.0_dp
integer                 :: i,iqapp,iqcor
real(kind=dp)           :: xmin,dlnx,stmp,k0,df,PN,sigma,frad
real(kind=dp)           :: factor,aicgm,gammq
real(kind=dp),allocatable::x(:),sang(:)

if(iqcor .eq. 1) then                           ! Gaussian
        aicgm  = 0.5_dp*(df-2.0_dp)
        xmin   = df * (k0/PN) ** (2.0/df)
        factor = xmin/16.0_dp/gamma(0.5_dp*df)
elseif(iqcor .eq. 2) then                       ! Exponential
        aicgm  = df - 2.0_dp
        xmin   = 2.0_dp * sqrt(0.5_dp*df*(df+1.0_dp))*(k0/PN)**(1.0/df)
        factor = xmin*xmin/16.0_dp/gamma(df)
elseif(iqcor .eq. 3) then                       ! Fractal dimension
        aicgm  = (df - 2.0_dp) / df
        xmin   = 2.0_dp**(df-1.0_dp) * k0 / PN
        factor = xmin ** (2.0_dp/df) / 16.0_dp
endif

if(df .gt. 2.0_dp .and. iqapp .eq. 3) then
        sigma  = factor*gamma(aicgm)*gammq(aicgm,xmin)
else
        dlnx = log(xmax/xmin)/real(nx-1,kind=dp)
        allocate(x(1:nx),sang(1:nx))
        do i=1,nx
                x(i) = exp(log(xmin)+real(i-1,kind=dp)*dlnx)
                if(iqapp .eq. 1) then
                        call angular_integration(iqcor,x(i),xmin,df,stmp)
                        sang(i) = stmp
                else
                        sang(i) = 1.0_dp
                endif
        enddo
        sigma = 0.0_dp
        do i=1,nx-1
        sigma = sigma + 0.5_dp * (frad(aicgm,x(i))*sang(i)+&
                frad(aicgm,x(i+1))*sang(i+1))
        enddo
        sigma = sigma * dlnx * factor
        deallocate(x,sang)
endif
return
end subroutine mean_overlap_efficiency

!--------------------------------------------------------------------------------
! The radial integrand function in ln(x) space
!--------------------------------------------------------------------------------
function frad(aicgm,x)
use types
implicit none
real(kind=dp)::aicgm,x,frad
        frad = x ** aicgm * exp(-x)
return
end function frad

!--------------------------------------------------------------------------------
! This subroutine computes the angular integral function S(rho)
!--------------------------------------------------------------------------------
subroutine angular_integration(iqcor,x,xmin,df,sang)
use types; use const
implicit none
integer            :: i,j,iqcor
integer,parameter  :: nmax_u = 500
real(kind=dp)      :: x,xmin,df,rho,umin,umax,du,fang,sang
real(kind=dp)      :: u(1:nmax_u)

if(iqcor .eq. 1) then
        rho = sqrt(x/xmin)
elseif(iqcor .eq. 2) then
        rho = x/xmin
elseif(iqcor .eq. 3) then
        rho = (x/xmin) ** (1.0_dp/df)
endif

sang = 0.0_dp
umin = 0.0_dp
umax = 1.0_dp / rho
du   = (umax - umin) / real(nmax_u-1,kind=dp)
do j=1,nmax_u
        u(j) = umin + du * real(j-1,kind=dp)
enddo
do j=1,nmax_u-1
        sang = sang + 0.5_dp*(fang(rho,u(j))+fang(rho,u(j+1)))*du
enddo
sang = 16.0_dp * rho * rho * sang / pi
return
end subroutine angular_integration

!--------------------------------------------------------------------------------
!
! This function returns the integrand of the angular integral.
!  fang (rho,u) =  [ Arcsin(sqrt(1-rho^2u^2)) - &
!               rho*u*sqrt(1-rho^2u^2)] * u / sqrt(1-u^2)
!
!--------------------------------------------------------------------------------
function fang(rho,u)
use types
implicit none
real(kind=dp)::rho,u,fang

if(abs(1.0_dp-rho*u) .le. 1.0e-10_dp) then
        fang = 0.0_dp
else
        fang = (asin(sqrt(1.0_dp-rho*rho*u*u))-&
                rho*u*sqrt(1.0_dp-rho*rho*u*u))*u/sqrt(1.0_dp-u*u)
endif

return
end function fang

!--------------------------------------------------------------------------------
!
! Following subroutines (gammq, gser, gcf, gammln) compute
! the incomplete Gamma function Q(a,x)
!
!             1     /infinity
! Q(a,x) = -------- | x^{a-1} exp(-x) dx
!          Gamma(a) /x
!
! where a>0. The definition of Q(a,x) obeys the definition in
!
! "Numerical Recipes in Fortran 77"
! Press, W. H., Teukolsky, S. A., Vetterling, W. T., et al.1992,
! Cambridge: University Press, |c1992, 2nd ed
!
! Following subroutines are taken from Numerical Recippes in Fortran 77
!
!--------------------------------------------------------------------------------
!
! History (R.T.)
!
! The original subroutine was slightly modified so that the old F77 grammers
! are rewritten by modern grammers.
!
! Convergence accuracy eps is changed from 3e-7 to 3e-14.
!
!--------------------------------------------------------------------------------
function gammq(a,x)
use types
real(kind=dp) a,gammq,x
real(kind=dp) gammcf,gamser,gln
if(x.lt.0..or.a.le.0.) then
      write(*,*) 'error'
      stop
endif
if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammq=1.-gamser
else
        call gcf(gammcf,a,x,gln)
        gammq=gammcf
endif
return
end function gammq

subroutine gser(gamser,a,x,gln)
use types
integer itmax
real(kind=dp) a,gamser,gln,x,eps
parameter (itmax=100,eps=3.e-14)
integer n
real(kind=dp) ap,del,sum,gammln
gln=gammln(a)
if(x.le.0.)then
        if(x.lt.0.) then
                write(*,*) "error"
                stop
        endif
        gamser=0.
        return
endif
ap=a
sum=1./a
del=sum
do n=1,itmax
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*eps) exit
        if(n .eq. itmax) then
                print *, 'a too large, itmax too small in gser'
                stop
        endif
enddo
gamser=sum*exp(-x+a*log(x)-gln)
return
end subroutine gser
subroutine gcf(gammcf,a,x,gln)
use types
integer itmax
real(kind=dp) a,gammcf,gln,x,eps,fpmin
parameter (itmax=100,eps=3.e-14,fpmin=1.e-30)
integer i
real(kind=dp) an,b,c,d,del,h,gammln
gln=gammln(a)
b=x+1.-a
c=1./fpmin
d=1./b
h=d
do i=1,itmax
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if(abs(d).lt.fpmin)d=fpmin
        c=b+an/c
        if(abs(c).lt.fpmin)c=fpmin
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.eps) exit
        if(n .eq. itmax) then
                print *, 'a too large, itmax too small in gcf'
                stop
        endif
enddo
gammcf=exp(-x+a*log(x)-gln)*h
return
end subroutine gcf
function gammln(xx)
use types
real(kind=dp) gammln,xx
integer j
real(kind=dp) ser,stp,tmp,x,y,cof(6)
save cof,stp
data cof,stp/76.18009172947146d0,-86.50532032941677d0,&
        24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,&
        -.5395239384953d-5,2.5066282746310005d0/
x=xx
y=x
tmp=x+5.5d0
tmp=(x+0.5d0)*log(tmp)-tmp
ser=1.000000000190015d0
do j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
enddo
gammln=tmp+log(stp*ser/x)
return
end function gammln

max_stars_repo_name: cdominik/optool


!
!  Copyright 2017 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
SUBROUTINE OUT_data_rt
use salmon_parallel, only: nproc_id_global, nproc_group_global, nproc_size_global
use salmon_communication, only: comm_is_root, comm_summation, comm_bcast
use scf_data
use writebox_rt_sub
use allocate_mat_sub
implicit none
integer       :: i1,i2,i3,jj,iob,is,it2,kk,iik
integer       :: ix,iy,iz
integer :: ibox
character(100) :: file_OUT_rt_data
integer :: ii,j1,j2,j3
integer :: nproc_id_global_datafiles
integer :: ista_Mxin_datafile(3)
integer :: iend_Mxin_datafile(3)
integer :: inum_Mxin_datafile(3)
integer :: nproc_xyz_datafile(3)
character(8) :: fileNumber_data
integer :: iob_myob
integer :: icorr_p

if(comm_is_root(nproc_id_global)) open(99,file=file_OUT_rt,form='unformatted')

if(comm_is_root(nproc_id_global)) write(99) itotNtime

if(num_datafiles_OUT>=2.and.num_datafiles_OUT<=nproc_size_global)then
  if(nproc_id_global<num_datafiles_OUT)then
    nproc_id_global_datafiles=nproc_id_global

    ibox=1
    nproc_xyz_datafile=1
    do ii=1,19
      do jj=3,1,-1
        if(ibox<num_datafiles_OUT)then
          nproc_xyz_datafile(jj)=nproc_xyz_datafile(jj)*2
          ibox=ibox*2
        end if
      end do
    end do

    do j3=0,nproc_xyz_datafile(3)-1
    do j2=0,nproc_xyz_datafile(2)-1
    do j1=0,nproc_xyz_datafile(1)-1
      ibox = j1 + nproc_xyz_datafile(1)*j2 + nproc_xyz_datafile(1)*nproc_xyz_datafile(2)*j3
      if(ibox==nproc_id_global_datafiles)then
        ista_Mxin_datafile(1)=j1*lg_num(1)/nproc_xyz_datafile(1)+lg_sta(1)
        iend_Mxin_datafile(1)=(j1+1)*lg_num(1)/nproc_xyz_datafile(1)+lg_sta(1)-1
        ista_Mxin_datafile(2)=j2*lg_num(2)/nproc_xyz_datafile(2)+lg_sta(2)
        iend_Mxin_datafile(2)=(j2+1)*lg_num(2)/nproc_xyz_datafile(2)+lg_sta(2)-1
        ista_Mxin_datafile(3)=j3*lg_num(3)/nproc_xyz_datafile(3)+lg_sta(3)
        iend_Mxin_datafile(3)=(j3+1)*lg_num(3)/nproc_xyz_datafile(3)+lg_sta(3)-1
      end if
    end do
    end do
    end do
    inum_Mxin_datafile(:)=iend_Mxin_datafile(:)-ista_Mxin_datafile(:)+1

    write(fileNumber_data, '(i6.6)') nproc_id_global_datafiles
    file_Out_rt_data = trim(adjustl(sysname))//"_rt_"//trim(adjustl(fileNumber_data))//".bin"
    open(89,file=file_Out_rt_data,form='unformatted')

  end if
end if

do iik=1,num_kpoints_rd
do iob=1,itotMST
  call calc_myob(iob,iob_myob)
  call check_corrkob(iob,iik,icorr_p)

  cmatbox_l2=0.d0
    if(mod(itotNtime,2)==1)then
      if(icorr_p==1)then
        do iz=mg_sta(3),mg_end(3)
        do iy=mg_sta(2),mg_end(2)
        do ix=mg_sta(1),mg_end(1)
          cmatbox_l2(ix,iy,iz)=zpsi_out(ix,iy,iz,iob_myob,iik)
        end do
        end do
        end do
      end if
    else
      if(icorr_p==1)then
        do iz=mg_sta(3),mg_end(3)
        do iy=mg_sta(2),mg_end(2)
        do ix=mg_sta(1),mg_end(1)
          cmatbox_l2(ix,iy,iz)=zpsi_in(ix,iy,iz,iob_myob,iik)
        end do
        end do
        end do
      end if
    end if

  call comm_summation(cmatbox_l2,cmatbox_l,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)

  if(num_datafiles_OUT==1.or.num_datafiles_OUT>nproc_size_global)then
    if(comm_is_root(nproc_id_global))then
      write(99) ((( cmatbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
  else
    if(nproc_id_global<num_datafiles_OUT)then
      write(89) ((( cmatbox_l(ix,iy,iz),ix=ista_Mxin_datafile(1),iend_Mxin_datafile(1)),   &
                                        iy=ista_Mxin_datafile(2),iend_Mxin_datafile(2)),   &
                                        iz=ista_Mxin_datafile(3),iend_Mxin_datafile(3))
    end if
  end if

end do
end do

if(num_datafiles_OUT>=2.and.num_datafiles_OUT<=nproc_size_global)then
  if(nproc_id_global<num_datafiles_OUT)then
    close(89)
  end if
end if

matbox_l=0.d0
do i3=ng_sta(3),ng_end(3)
do i2=ng_sta(2),ng_end(2)
do i1=ng_sta(1),ng_end(1)
   matbox_l(i1,i2,i3)=rho(i1,i2,i3)
end do
end do
end do

call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
if(comm_is_root(nproc_id_global))then
  write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if

if(ilsda == 1)then
  do is=1,2

    matbox_l=0.d0
    do i3=ng_sta(3),ng_end(3)
    do i2=ng_sta(2),ng_end(2)
    do i1=ng_sta(1),ng_end(1)
      matbox_l(i1,i2,i3)=rho_s(i1,i2,i3,is)
    end do
    end do
    end do

    call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
    if(comm_is_root(nproc_id_global))then
      write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
  end do
end if

matbox_l=0.d0
do i3=ng_sta(3),ng_end(3)
do i2=ng_sta(2),ng_end(2)
do i1=ng_sta(1),ng_end(1)
   matbox_l(i1,i2,i3)=Vh(i1,i2,i3)
end do
end do
end do
call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
if(comm_is_root(nproc_id_global))then
  write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if

if(ilsda == 0)then
  matbox_l=0.d0
  do i3=ng_sta(3),ng_end(3)
  do i2=ng_sta(2),ng_end(2)
  do i1=ng_sta(1),ng_end(1)
    matbox_l(i1,i2,i3)=Vxc(i1,i2,i3)
  end do
  end do
  end do
  call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
  if(comm_is_root(nproc_id_global))then
    write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
  end if
else if(ilsda == 1)then
  do is=1,2

    matbox_l=0.d0
    do i3=ng_sta(3),ng_end(3)
    do i2=ng_sta(2),ng_end(2)
    do i1=ng_sta(1),ng_end(1)
      matbox_l(i1,i2,i3)=Vxc_s(i1,i2,i3,is)
    end do
    end do
    end do

    call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
    if(comm_is_root(nproc_id_global))then
      write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
  end do
end if

matbox_l=0.d0
do i3=ng_sta(3),ng_end(3)
do i2=ng_sta(2),ng_end(2)
do i1=ng_sta(1),ng_end(1)
   matbox_l(i1,i2,i3)=Vh_stock1(i1,i2,i3)
end do
end do
end do
call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
if(comm_is_root(nproc_id_global))then
  write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if
matbox_l=0.d0
do i3=ng_sta(3),ng_end(3)
do i2=ng_sta(2),ng_end(2)
do i1=ng_sta(1),ng_end(1)
   matbox_l(i1,i2,i3)=Vh_stock2(i1,i2,i3)
end do
end do
end do
call comm_summation(matbox_l,matbox_l2,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
if(comm_is_root(nproc_id_global))then
  write(99) ((( matbox_l2(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if

if(comm_is_root(nproc_id_global)) then

  write(99) (vecDs(jj),jj=1,3)
  do it2=1,itotNtime
     write(99) (Dp(jj,it2),jj=1,3)
  end do
  do it2=1,itotNtime
     write(99) tene(it2)
  end do
  if(iflag_dip2==1)then
    write(99) ((vecDs2(jj,kk),jj=1,3),kk=1,num_dip2)
    do it2=1,itotNtime
      write(99) ((Dp2(jj,it2,kk),jj=1,3),kk=1,num_dip2)
    end do
  end if

  close(99)

end if

END SUBROUTINE OUT_data_rt

!---------------------------------------------------------------------------
SUBROUTINE IN_data_rt(Ntime)
use salmon_parallel, only: nproc_id_global, nproc_group_global, nproc_size_global
use salmon_communication, only: comm_is_root, comm_summation, comm_bcast
use scf_data
use new_world_sub
use readbox_rt_sub
use allocate_mat_sub
implicit none
integer       :: i1,i2,i3,jj,iob,is,it2,kk,iik
integer       :: ix,iy,iz
integer       :: Ntime
character(100) :: file_IN_rt_data
integer :: ibox
integer :: ii,j1,j2,j3
integer :: nproc_id_global_datafiles
integer :: ista_Mxin_datafile(3)
integer :: iend_Mxin_datafile(3)
integer :: inum_Mxin_datafile(3)
integer :: nproc_xyz_datafile(3)
character(8) :: fileNumber_data
integer :: iob_myob
integer :: icorr_p

if(comm_is_root(nproc_id_global))then

   open(98,file=file_IN_rt,form='unformatted')

end if

if(comm_is_root(nproc_id_global)) read(98) Miter_rt
call comm_bcast(Miter_rt,nproc_group_global)
itotNtime=Ntime+Miter_rt

allocate(rIe(0:itotNtime))
allocate(Dp(3,0:itotNtime))
allocate(tene(0:itotNtime))
if(iflag_dip2==1) then
    allocate(rIe2(0:itotNtime,1:num_dip2))
    allocate(Dp2(3,0:itotNtime,1:num_dip2))
    allocate(Qp2(3,3,0:itotNtime,1:num_dip2))
end if
allocate(Qp(3,3,0:itotNtime))

if(num_datafiles_IN>=2.and.num_datafiles_IN<=nproc_size_global)then
  if(nproc_id_global<num_datafiles_IN)then
    nproc_id_global_datafiles=nproc_id_global

    ibox=1
    nproc_xyz_datafile=1
    do ii=1,19
      do jj=3,1,-1
        if(ibox<num_datafiles_IN)then
          nproc_xyz_datafile(jj)=nproc_xyz_datafile(jj)*2
          ibox=ibox*2
        end if
      end do
    end do

    do j3=0,nproc_xyz_datafile(3)-1
    do j2=0,nproc_xyz_datafile(2)-1
    do j1=0,nproc_xyz_datafile(1)-1
      ibox = j1 + nproc_xyz_datafile(1)*j2 + nproc_xyz_datafile(1)*nproc_xyz_datafile(2)*j3
      if(ibox==nproc_id_global_datafiles)then
        ista_Mxin_datafile(1)=j1*lg_num(1)/nproc_xyz_datafile(1)+lg_sta(1)
        iend_Mxin_datafile(1)=(j1+1)*lg_num(1)/nproc_xyz_datafile(1)+lg_sta(1)-1
        ista_Mxin_datafile(2)=j2*lg_num(2)/nproc_xyz_datafile(2)+lg_sta(2)
        iend_Mxin_datafile(2)=(j2+1)*lg_num(2)/nproc_xyz_datafile(2)+lg_sta(2)-1
        ista_Mxin_datafile(3)=j3*lg_num(3)/nproc_xyz_datafile(3)+lg_sta(3)
        iend_Mxin_datafile(3)=(j3+1)*lg_num(3)/nproc_xyz_datafile(3)+lg_sta(3)-1
      end if
    end do
    end do
    end do
    inum_Mxin_datafile(:)=iend_Mxin_datafile(:)-ista_Mxin_datafile(:)+1

    write(fileNumber_data, '(i6.6)') nproc_id_global_datafiles
    file_IN_rt_data = trim(adjustl(sysname))//"_rt_"//trim(adjustl(fileNumber_data))//".bin"
    open(88,file=file_IN_rt_data,form='unformatted')

  end if
end if

do iik=k_sta,k_end
do iob=1,itotMST
  call calc_myob(iob,iob_myob)
  call check_corrkob(iob,iik,icorr_p)

  if(num_datafiles_IN==1.or.num_datafiles_IN>nproc_size_global)then
    if(comm_is_root(nproc_id_global))then
      read(98) ((( cmatbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
    call comm_bcast(cmatbox_l(lg_sta(1):lg_sta(1)+lg_num(1),lg_sta(2):lg_sta(2)+lg_num(2),lg_sta(3):lg_sta(3)+lg_num(3)),  &
                    nproc_group_global)
  else
    cmatbox_l2=0.d0
    if(nproc_id_global<num_datafiles_IN)then
      read(88) ((( cmatbox_l2(ix,iy,iz),ix=ista_Mxin_datafile(1),iend_Mxin_datafile(1)),   &
                                        iy=ista_Mxin_datafile(2),iend_Mxin_datafile(2)),   &
                                        iz=ista_Mxin_datafile(3),iend_Mxin_datafile(3))
    end if
    call comm_summation(cmatbox_l2,cmatbox_l,lg_num(1)*lg_num(2)*lg_num(3),nproc_group_global)
  end if
  if(mod(Miter_rt,2)==1)then
    if(icorr_p==1)then
      zpsi_out(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),  &
          mg_sta(3):mg_end(3),iob_myob,iik)=  &
      cmatbox_l(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),   &
             mg_sta(3):mg_end(3))
    end if
  else
    if(icorr_p==1)then
      zpsi_in(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),  &
          mg_sta(3):mg_end(3),iob_myob,iik)=  &
      cmatbox_l(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),   &
             mg_sta(3):mg_end(3))
    end if
  end if

end do
end do

if(num_datafiles_IN>=2.and.num_datafiles_IN<=nproc_size_global)then
  if(nproc_id_global<num_datafiles_IN)then
    close(88)
  end if
end if

if(comm_is_root(nproc_id_global))then
  read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if
call comm_bcast(matbox_l,nproc_group_global)
do i3=mg_sta(3),mg_end(3)
do i2=mg_sta(2),mg_end(2)
do i1=mg_sta(1),mg_end(1)
  rho(i1,i2,i3)=matbox_l(i1,i2,i3)
end do
end do
end do
if(ilsda == 1)then
  do is=1,2
    if(comm_is_root(nproc_id_global))then
      read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
    call comm_bcast(matbox_l,nproc_group_global)
    do i3=lg_sta(3),lg_end(3)
    do i2=lg_sta(2),lg_end(2)
    do i1=lg_sta(1),lg_end(1)
      rho_s(i1,i2,i3,is)=matbox_l(i1,i2,i3)
    end do
    end do
    end do
  end do
end if

if(comm_is_root(nproc_id_global))then
  read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if
call comm_bcast(matbox_l,nproc_group_global)
do iz=mg_sta(3),mg_end(3)
do iy=mg_sta(2),mg_end(2)
do ix=mg_sta(1),mg_end(1)
  Vh(ix,iy,iz)=matbox_l(ix,iy,iz)
end do
end do
end do

if(ilsda == 0)then
  if(comm_is_root(nproc_id_global))then
    read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
  end if
  call comm_bcast(matbox_l,nproc_group_global)
  do iz=mg_sta(3),mg_end(3)
  do iy=mg_sta(2),mg_end(2)
  do ix=mg_sta(1),mg_end(1)
    Vxc(ix,iy,iz)=matbox_l(ix,iy,iz)
  end do
  end do
  end do
else if(ilsda == 1)then
  do is=1,2
    if(comm_is_root(nproc_id_global))then
      read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
    end if
    call comm_bcast(matbox_l,nproc_group_global)
    do iz=mg_sta(3),mg_end(3)
    do iy=mg_sta(2),mg_end(2)
    do ix=mg_sta(1),mg_end(1)
      Vxc_s(ix,iy,iz,is)=matbox_l(ix,iy,iz)
    end do
    end do
    end do

  end do
end if

if(comm_is_root(nproc_id_global))then
  read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if
call comm_bcast(matbox_l,nproc_group_global)
do iz=mg_sta(3),mg_end(3)
do iy=mg_sta(2),mg_end(2)
do ix=mg_sta(1),mg_end(1)
  Vh_stock1(ix,iy,iz)=matbox_l(ix,iy,iz)
end do
end do
end do

if(comm_is_root(nproc_id_global))then
  read(98) ((( matbox_l(ix,iy,iz),ix=lg_sta(1),lg_end(1)),iy=lg_sta(2),lg_end(2)),iz=lg_sta(3),lg_end(3))
end if
call comm_bcast(matbox_l,nproc_group_global)
do iz=mg_sta(3),mg_end(3)
do iy=mg_sta(2),mg_end(2)
do ix=mg_sta(1),mg_end(1)
  Vh_stock2(ix,iy,iz)=matbox_l(ix,iy,iz)
end do
end do
end do

itt=Miter_rt
call allgatherv_vlocal

if(comm_is_root(nproc_id_global))then
  read(98) (vecDs(jj),jj=1,3)
  do it2=1,Miter_rt
    read(98) (Dp(jj,it2),jj=1,3)
  end do
  do it2=1,Miter_rt
    read(98) tene(it2)
  end do
end if

call comm_bcast(vecDs,nproc_group_global)

do jj=1,3
do it2=1,Miter_rt
   call comm_bcast(Dp(jj,it2),nproc_group_global)
end do
end do
do it2=1,Miter_rt
   call comm_bcast(tene(it2),nproc_group_global)
end do

if(iflag_dip2==1)then
  if(comm_is_root(nproc_id_global))then
    read(98) ((vecDs2(jj,kk),jj=1,3),kk=1,num_dip2)
    do it2=1,Miter_rt
      read(98) ((Dp2(jj,it2,kk),jj=1,3),kk=1,num_dip2)
    end do
  end if
  call comm_bcast(vecDs2,nproc_group_global)
  do kk=1,num_dip2
    call comm_bcast(Dp2(:,:,kk),nproc_group_global)
  end do
end if

if(comm_is_root(nproc_id_global)) close(98)

END SUBROUTINE IN_data_rt