text:
      subroutine asidgif
     & (
     &   res,
     &   y,        x,       umesh,
     &   shpif0,   shpif1,  shgif0,  shgif1,
     &   qwtif0,   qwtif1,
     &   ienif0,   ienif1,
     &   sum_vi_area, if_normals
     & )
          use hierarchic_m
          use local_m
          use e3if_m
          use e3if_geom_m
c
          implicit none
c
          real*8, dimension(nshg,nflow), intent(inout) :: res
          real*8, dimension(nshg,ndof),  intent(in)    :: y
          real*8, dimension(nshg,nsd),   intent(in)    :: x
          real*8, dimension(nshg, nsd), intent(inout) :: umesh
          real*8, dimension(nshl0,nqpt),intent(in)   :: shpif0
          real*8, dimension(nshl1,nqpt),intent(in)   :: shpif1
          real*8, dimension(nsd,nshl0,nqpt), intent(in)  :: shgif0
          real*8, dimension(nsd,nshl1,nqpt), intent(in)  :: shgif1
          real*8, dimension(nqpt), intent(in) :: qwtif0, qwtif1
          integer, dimension(:,:), pointer, intent(in)   :: ienif0, ienif1
          real*8, pointer, intent(inout) :: sum_vi_area(:,:)
          real*8, pointer, intent(in) :: if_normals(:,:)
c
          call malloc_e3if
          call malloc_e3if_geom
c
c.... create the matrix of mode signs for the hierarchic basis
c     functions.
c
          if (ipord .gt. 1) then
           call getsgn(ienif0,sgn0,nshl0,nenl0)
           call getsgn(ienif1,sgn1,nshl1,nenl1)
        endif
c
c... localize
c
        call localy(y, ycl0, ienif0, ndof, 'gather  ', nshg, nshl0, npro, ipord, gbytes)
        call localy(y, ycl1, ienif1, ndof, 'gather  ', nshg, nshl1, npro, ipord, gbytes)
c
        call localx(x, xl0,  ienif0, nsd, 'gather  ', nshg, nenl0, npro, gbytes)
        call localx(x, xl1,  ienif1, nsd, 'gather  ', nshg, nenl1, npro, gbytes)
c
        call localx(umesh, umeshl0,  ienif0, nsd, 'gather  ', nshg, nenl0, npro, gbytes)
        call localx(umesh, umeshl1,  ienif1, nsd, 'gather  ', nshg, nenl1, npro, gbytes)
c
        call localx(if_normals, if_normal_l0,  ienif0, nsd, 'gather  ', nshg, nenl0, npro, gbytes)
        call localx(if_normals, if_normal_l1,  ienif1, nsd, 'gather  ', nshg, nenl1, npro, gbytes)
c
       call e3if(shpif0,shpif1,shgif0,shgif1,qwtif0,qwtif1)
c
c.... assemble the local residual arrays
c
        call local (res, rl0, ienif0, nflow, 'scatter ', nshg,nshl0,npro,ipord,sbytes,flops)
        call local (res, rl1, ienif1, nflow, 'scatter ', nshg,nshl1,npro,ipord,sbytes,flops)
c
        call local (sum_vi_area, sum_vi_area_l0, ienif0, nsd+1, 'scatter ', nshg, nshl0,npro,ipord,sbytes,flops)
        call local (sum_vi_area, sum_vi_area_l1, ienif1, nsd+1, 'scatter ', nshg, nshl1,npro,ipord,sbytes,flops)
c
        call mfree_e3if
        call mfree_e3if_geom
c
      end subroutine asidgif

ext: f


text:
!Author: Huang Yihan
!This module is used to read some input parameters
!include: Radius, length of side, dancoff(This could be calculated)
!         Ressonance table file name
!The 69 group of energy structure of WIMS is given in this module
module input_mod
    implicit none
    real*8,parameter::WIMS_boundary(70)=(/&
    1.00E+07,6.0655e6,3.679e6,2.231e6,1.353e6,&
    8.21E+05,5.e5,3.025e5,1.83e5,1.11e5,6.734e4, &
    4.09E+04,2.478e4,1.503e4,9118.e0,5530.e0,&
    3.52E+03,2239.45e0,1425.1e0,906.898e0,367.262e0,&
    1.49E+02,75.501e0,48.052e0,27.7e0,15.968e0,9.88E+00,&
    4.e0,3.3e0,2.6e0,2.1e0,1.5e0,1.3e0,1.15E+00,1.123e0,&
    1.097e0,1.071e0,1.045e0,1.02e0,9.96E-01,.972e0,.95e0,&
    .91e0,.85e0,.78e0,6.25E-01,.5e0,.4e0,.35e0,.32e0,.3e0,&
    .28e0,2.50E-01,.22e0,.18e0,.14e0,.1e0,.08e0,6.70E-02,&
    .058e0,.05e0,.042e0,.035e0,.03e0,2.50E-02,.02e0,.015e0,&
     .01e0,.005e0,1.e-5/)    !the energy structure of WIMS 69 groups
    character*5,parameter    ::   input='input'!the filename of inputcard
    character*30             ::   libname      !the lib name
    real*8                   ::   radius       !the radius of pin
    real*8                   ::   length       !the length of side
    real*8                   ::   dancoff(2)   !Read dancoff
    real*8                   ::   density_H1   !the density of H1
    real*8                   ::   density_U8   !the density of U8
    integer*4,parameter      ::   N1=14        !the number of fast groups
    integer*4,parameter      ::   N2=13        !the number of res groups
    integer*4,parameter      ::   first_tag=15 !the first resonance group tag
    contains
        !===========================================
        !<Read inputs
        subroutine read_inputs
            call read_lib(libname)
            call read_radius(radius)
            call read_dancoff(dancoff(1),dancoff(2))
            call read_densityU8(density_U8)
            call read_densityU8(density_H1)
        end subroutine read_inputs
        !===========================================
        !<Read Filename subroutine
        subroutine read_lib(libname)
            character*30,intent(inout)   :: libname
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Lib:')
                if (k .gt. 0) then
                    read(words(k+4:80),*) libname
                    exit
                end if
            end do
            close(21)
        end subroutine read_lib
        !>===========================================
        !<Read Radius subroutine
        subroutine read_radius(radius)
            real*8,intent(inout)         :: radius
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Radius:')
                if (k .gt. 0) then
                    read(words(k+7:80),*) radius
                    exit
                end if
            end do
            close(21)
        end subroutine read_radius
        !>============================================
        !<Read Length of side subroutine
        subroutine read_length(length)
            real*8,intent(inout)         :: length
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Length:')
                if (k .gt. 0) then
                    read(words(k+7:80),*) length
                    exit
                end if
            end do
            close(21)
        end subroutine read_length
        !>============================================
        !<Read Temperature subroutine
        subroutine read_temperature(temperature)
            real*8,intent(inout)         :: temperature
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Temperature:')
                if (k .gt. 0) then
                    read(words(k+13:80),*) temperature
                    exit
                end if
            end do
            close(21)
        end subroutine read_temperature
        !>============================================
        !<Read Dancoff subroutine
        subroutine read_Dancoff(dancoff1,dancoff2)
            real*8,intent(inout)         :: dancoff1
            real*8,intent(inout)         :: dancoff2
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Dancoff:')
                if (k .gt. 0) then
                    read(words(k+8:80),*) dancoff1,dancoff2
                    exit
                end if
            end do
            close(21)
        end subroutine read_Dancoff
        !>============================================
        !<Read density_U8 subroutine
        subroutine read_densityU8(density_U8)
            real*8,intent(inout)         :: density_U8
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Density_U8:')
                if (k .gt. 0) then
                    read(words(k+11:80),*) density_U8
                    exit
                end if
            end do
            close(21)
        end subroutine read_densityU8
        !>============================================
        !<Read density_H1 subroutine
        subroutine read_densityH1(density_H1)
            real*8,intent(inout)         :: density_H1
            character*80                 :: words
            integer*4                    :: k
            open(21,file=input)
            do
                read(21,'(A80)') words
                k=index(words,'Density_U8:')
                if (k .gt. 0) then
                    read(words(k+11:80),*) density_H1
                    exit
                end if
            end do
            close(21)
        end subroutine read_densityH1
end module input_mod

max_stars_repo_name: Huang-Yihan/Recentyl-Physicals-homework

max_issues_repo_name: Huang-Yihan/Recentyl-Physicals-homework

max_forks_repo_name: Huang-Yihan/Recentyl-Physicals-homework


text:
***********************************************************************
c                    ORBEL_EHIE.F
***********************************************************************
*     PURPOSE:  Solves Kepler's eqn.   e is ecc.   m is mean anomaly.
*
*             Input:
*                           e ==> eccentricity anomaly. (real scalar)
*                           m ==> mean anomaly. (real scalar)
*             Returns:
*              orbel_ehybrid ==>  eccentric anomaly. (real scalar)
*
*     ALGORITHM: Use Danby's quartic for 3 iterations.
*                Eqn. is f(x) = x - e*sin(x+M). Note  that
*	         E = x + M. First guess is very good for e near 1.
*	         Need to first get M between 0. and PI and use
*		 symmetry to return right answer if M between PI and 2PI
*     REMARKS: Modifies M so that both E and M are in range (0,TWOPI)
*     AUTHOR: M. Duncan
*     DATE WRITTEN: May 25,1992.
*     REVISIONS:
***********************************************************************

      real*8 function orbel_ehie(e,m)

      include '../swift.inc'

c...  Inputs Only:
	real*8 e,m

c...  Internals:
      integer iflag,nper,niter,NMAX
      real*8 dx,x,sa,ca,esa,eca,f,fp

      parameter (NMAX = 3)

c----
c...  Executable code

c In this section, bring M into the range (0,TWOPI) and if
c the result is greater than PI, solve for (TWOPI - M).
	iflag = 0
	nper = m/TWOPI
	m = m - nper*TWOPI
	if (m .lt. 0.d0) m = m + TWOPI

	if (m.gt.PI) then
	   m = TWOPI - m
	   iflag = 1
	endif

c Make a first guess that works well for e near 1.
	x = (6.d0*m)**(1.d0/3.d0) - m
	niter =0

c Iteration loop
	do niter =1,NMAX
	    call orbel_scget(x + m,sa,ca)
	    esa = e*sa
	    eca = e*ca
	    f = x - esa
	    fp = 1.d0 -eca
	    dx = -f/fp
	    dx = -f/(fp + 0.5d0*dx*esa)
	    dx = -f/(fp + 0.5d0*dx*(esa+0.3333333333333333d0*eca*dx))
	    x = x + dx
	enddo

	orbel_ehie = m + x

	if (iflag.eq.1) then
	  orbel_ehie = TWOPI - orbel_ehie
	  m = TWOPI - m
	endif

	return
	end         !orbel_ehie
c------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_types_banddos
  USE m_juDFT
  USE m_types_fleurinput_base
  IMPLICIT NONE
  PRIVATE
  PUBLIC:: t_banddos
  TYPE,EXTENDS(t_fleurinput_base):: t_banddos
     LOGICAL :: dos =.FALSE.
     LOGICAL :: band =.FALSE.
     LOGICAL :: l_mcd =.FALSE.
     LOGICAL :: l_orb =.FALSE.
     LOGICAL :: l_jDOS = .FALSE.
     LOGICAL :: vacdos =.FALSE.
     INTEGER :: ndir =0
     INTEGER :: orbCompAtom=0
     INTEGER :: jDOSAtom=0
     REAL    :: e1_dos=0.5
     REAL    :: e2_dos=-0.5
     REAL    :: sig_dos=0.015
     INTEGER :: projdos !selects one atomtype and prints the projected dos if there are to many atoms
     REAL    :: e_mcd_lo =-10.0
     REAL    :: e_mcd_up= 0.0
     LOGICAL :: unfoldband =.FALSE.
     INTEGER :: s_cell_x=1
     INTEGER :: s_cell_y=1
     INTEGER :: s_cell_z=1
     REAL    :: alpha,beta,gamma !For orbital decomp. (was orbcomprot)
   CONTAINS
     PROCEDURE :: read_xml=>read_xml_banddos
     PROCEDURE :: mpi_bc=>mpi_bc_banddos
  END TYPE t_banddos
CONTAINS
  SUBROUTINE mpi_bc_banddos(this,mpi_comm,irank)
    USE m_mpi_bc_tool
    CLASS(t_banddos),INTENT(INOUT)::this
    integer,INTENT(IN):: mpi_comm
    INTEGER,INTENT(IN),OPTIONAL::irank
    INTEGER ::rank
    if (present(irank)) THEN
       rank=irank
    else
       rank=0
    end if
    CALL mpi_bc(this%dos ,rank,mpi_comm)
    CALL mpi_bc(this%band ,rank,mpi_comm)
    CALL mpi_bc(this%l_mcd ,rank,mpi_comm)
    CALL mpi_bc(this%l_orb ,rank,mpi_comm)
    CALL mpi_bc(this%l_jDOS,rank,mpi_comm)
    CALL mpi_bc(this%vacdos ,rank,mpi_comm)
    CALL mpi_bc(this%ndir ,rank,mpi_comm)
    CALL mpi_bc(this%orbCompAtom,rank,mpi_comm)
    CALL mpi_bc(this%jDOSAtom,rank,mpi_comm)
    CALL mpi_bc(this%e1_dos,rank,mpi_comm)
    CALL mpi_bc(this%e2_dos,rank,mpi_comm)
    CALL mpi_bc(this%sig_dos,rank,mpi_comm)
    CALL mpi_bc(this%e_mcd_lo ,rank,mpi_comm)
    CALL mpi_bc(this%e_mcd_up,rank,mpi_comm)
    CALL mpi_bc(this%unfoldband ,rank,mpi_comm)
    CALL mpi_bc(this%s_cell_x,rank,mpi_comm)
    CALL mpi_bc(this%s_cell_y,rank,mpi_comm)
    CALL mpi_bc(this%s_cell_z,rank,mpi_comm)
    CALL mpi_bc(this%alpha,rank,mpi_comm)
    CALL mpi_bc(this%beta,rank,mpi_comm)
    CALL mpi_bc(this%gamma,rank,mpi_comm)

  END SUBROUTINE mpi_bc_banddos
  SUBROUTINE read_xml_banddos(this,xml)
    USE m_types_xml
    CLASS(t_banddos),INTENT(INOUT)::this
    TYPE(t_xml),INTENT(INOUT)::xml

    CHARACTER(len=300) :: xPathA, xPathB
    INTEGER::numberNodes,iType,i,na
    LOGICAL::l_orbcomp,l_jDOS
    this%band = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@band'))
    this%dos = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@dos'))
    this%vacdos = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@vacdos'))
    this%l_mcd = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/@mcd'))

    numberNodes = xml%GetNumberOfNodes('/fleurInput/output/densityOfStates')

    IF ((this%dos).AND.(numberNodes.EQ.0)) THEN
       CALL juDFT_error("dos is true but densityOfStates parameters are not set!")
    END IF

    IF (numberNodes.EQ.1) THEN
       this%ndir = evaluateFirstIntOnly(xml%GetAttributeValue('/fleurInput/output/densityOfStates/@ndir'))
       this%e2_dos = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/densityOfStates/@minEnergy'))
       this%e1_dos = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/densityOfStates/@maxEnergy'))
       this%sig_dos = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/densityOfStates/@sigma'))
    END IF
    IF (this%band) THEN
       this%dos=.TRUE.
       this%ndir = -4
       WRITE(*,*) 'band="T" --> Overriding "dos" and "ndir"!'
    ENDIF

    ! Read in optional magnetic circular dichroism parameters
    numberNodes = xml%GetNumberOfNodes('/fleurInput/output/magneticCircularDichroism')

    IF ((this%l_mcd).AND.(numberNodes.EQ.0)) THEN
       CALL juDFT_error("mcd is true but magneticCircularDichroism parameters are not set!", calledby = "read_xml_banddos")
    END IF

    IF (numberNodes.EQ.1) THEN
       this%e_mcd_lo = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/magneticCircularDichroism/@energyLo'))
       this%e_mcd_up = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/magneticCircularDichroism/@energyUp'))
    END IF

    na = 0
    DO iType = 1, xml%GetNumberOfNodes('/fleurInput/atomGroups/atomGroup')
       WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']'
       DO i = 1, xml%GetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos')
          na = na + 1
          WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/relPos[',i,']'
          l_orbcomp = evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp'))
          IF(l_orbcomp) THEN
             IF(this%l_orb) THEN
                CALL juDFT_error("Multiple orbcomp flags set.", calledby = "read_xml_banddos")
             END IF
             this%l_orb = .TRUE.
             this%orbCompAtom = na
          ENDIF
          l_jDOS = evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@jDOS'))
          IF(l_jDOS) THEN
             IF(this%l_jDOS) THEN
                CALL juDFT_error("Multiple jDOS flags set.", calledby="read_xml_banddos")
             ENDIF
             this%l_jDOS = .TRUE.
             this%jDOSAtom = na
          ENDIF
       ENDDO
    ENDDO
    IF(this%l_orb.AND.this%l_jDOS) THEN
       CALL juDFT_error("Both jDOS and orbcomp flag set", calledby="read_xml_banddos")
    ENDIF

    ! Read in optional parameter for unfolding bandstructure of supercell
    numberNodes = xml%GetNumberOfNodes('/fleurInput/output/unfoldingBand')

    IF (numberNodes.EQ.1) THEN
       this%unfoldband = evaluateFirstBoolOnly(xml%GetAttributeValue('/fleurInput/output/unfoldingBand/@unfoldBand'))
       this%s_cell_x = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/unfoldingBand/@supercellX'))
       this%s_cell_y = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/unfoldingBand/@supercellY'))
       this%s_cell_z = evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/output/unfoldingBand/@supercellZ'))
    END IF
  END SUBROUTINE read_xml_banddos

END MODULE m_types_banddos

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
      SUBROUTINE CFACTR (A,LL,UL,SCR1,SCR2,SCR3,IOPT)
C
      INTEGER           FA        ,FL       ,FU       ,SR1      ,
     1                  SR2       ,SR3      ,UL       ,SCR1     ,
     2                  SCR2      ,SCR3     ,A        ,
     3                  MCB(7)    ,NAME(2)
      DOUBLE PRECISION  DET       ,MIND
      COMMON   /CDCMPX/ FA(7)     ,FL(7)    ,FU(7)    ,SR1      ,
     1                  SR2       ,SR3      ,DET(2)   ,POWR     ,
     2                  NX        ,MIND     ,IB       ,IBBAR
      COMMON   /SFACT / MFA(7)    ,MFL(7)   ,MFC(7)   ,M1FIL    ,
     1                  M2FIL     ,MXX      ,D(5)     ,M3FIL    ,
     2                  D1(2)     ,ICHOL
      COMMON   /SDCCSP/ JFA(7)    ,JFL(7)   ,JFC(7)   ,J1FIL    ,
     1                  J2FIL     ,JX
      COMMON   /ZZZZZZ/ IZ(1)
      DATA      NAME  / 4HCFAC,4HTR   /
C
C
      NZ = KORSZ(IZ)
      MCB(1) = A
      CALL RDTRL (MCB)
      IF (MCB(4) .NE. 6) GO TO 200
C
C     SYMMETRIC  COMPLEX
C
      DO 10 I = 1,7
      MFA(I) = MCB(I)
      MFL(I) = MCB(I)
      MFC(I) = MCB(I)
   10 CONTINUE
      MFL(1) = LL
      MFC(1) = UL
      MFL(4) = 4
      MFC(4) = 5
      M1FIL  = SCR1
      M2FIL  = SCR2
      MXX    = NZ
      M3FIL  = SCR3
      ICHOL  = 0
      CALL SDCOMP (*900,IZ,IZ,IZ)
      CALL WRTTRL (MFL)
      IOPT  = 2
      GO TO  60
C
C     UNSYMMETRIC  COMPLEX
C
  200 DO 210 I = 1,7
      FA(I) = MCB(I)
      FL(I) = MCB(I)
      FU(I) = MCB(I)
  210 CONTINUE
      FL(1) = LL
      FU(1) = UL
      FL(4) = 4
      FU(4) = 5
      SR1   = SCR1
      SR2   = SCR2
      SR3   = SCR3
      NX    = NZ
C     IB    = 0
C
C     IF IB IS SET TO ZERO HERE, T08021 PRINTS 27 MORE MESSAGES 3027
C     AND 3028 FROM GENVEC WHICH IS CALLED BY CFACTR, WHCIH IS CALLED BY
C     FRD2C, IN FRRD2 MODULE
C
CIBMI 6/93
      IBBAR = 0
      CALL CDCOMP (*900,IZ,IZ,IZ)
      CALL WRTTRL (FU)
      CALL WRTTRL (FL)
      IOPT  = 1
   60 RETURN
C
C     ERRORS
C
  900 CALL MESAGE (-5,A,NAME)
C
      RETURN
      END


text:
!-----------------------------------------------------------------------
! Fortran interface to CUDA utility Library
! written by Viktor K. Decyk, UCLA
      module gpulib2_h
      implicit none
!
      interface
         subroutine gpu_setgbsize(nblock)
         implicit none
         integer :: nblock
         end subroutine
      end interface
!
      interface
         integer function getmmcc()
         implicit none
         end function
      end interface
!
      interface
         subroutine gpu_fallocate(gp_f,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_iallocate(gp_i,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_i
         end subroutine
      end interface
!
      interface
         subroutine gpu_callocate(gp_c,nsize,irc)
         implicit none
         integer :: nsize, irc
         integer, dimension(2) :: gp_c
         end subroutine
      end interface
!
      interface
         subroutine gpu_deallocate(gp_d,irc)
         implicit none
         integer :: irc
         integer, dimension(2) :: gp_d
         end subroutine
      end interface
!
      interface
         subroutine gpu_fcopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         real, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_fcopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         real, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_icopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_icopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_ccopyin(f,gp_f,nsize)
         implicit none
         integer :: nsize
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_ccopyout(f,gp_f,nsize)
         implicit none
         integer :: nsize
         complex, dimension(*) :: f
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_zfmem(gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_zcmem(gp_f,nsize)
         implicit none
         integer :: nsize
         integer, dimension(2) :: gp_f
         end subroutine
      end interface
!
      interface
         subroutine gpu_set_cache_size(nscache)
         implicit none
         integer :: nscache
         end subroutine
      end interface
!
      interface
         subroutine emptykernel()
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine init_cu(dev,irc)
         implicit none
         integer :: dev, irc
         end subroutine
      end interface
!
      interface
         subroutine end_cu()
         implicit none
         end subroutine
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


module integrators

use, intrinsic :: iso_fortran_env

implicit none

! Fortran 2008
integer, parameter, private :: sp = REAL32
integer, parameter, private :: dp = REAL64
integer, parameter, private :: qp = REAL128

contains

      subroutine propagate(ptcls,tini,tfin,nsteps,indevar)
      use mpi
      use evalroutines, only : evalt
      implicit none
      real*8, dimension(:,:) :: ptcls !particle array
      real*8 :: tini,tfin
      integer :: nsteps,indevar
      integer :: mprocs,myrank,mpierr
      call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,mpierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,mpierr)
      if(indevar.eq.0)call adams11(ptcls,tini,tfin,nsteps,evalt)
      return
      end

      subroutine rk4(y,tini,tfin,nsteps,eval)
      use mpi
      use lienwiech, only : getfield_atallpts
      use diagnostics, only : getcentroid
      use lwprobcons, only : tauhistmin,tauhistmax,refptcltime,laststoredrefpart,testtoperformfieldcalc,istopafter1fieldcalc
      implicit none
      real*8, dimension(:,:) :: y
      real*8 :: tini,tfin
      integer :: nsteps,mlo,mhi,mlomin,mhimax,gblmlomin,gblmhimax
      real*8, allocatable, dimension(:,:) :: yt,f,a,b,c,d
      integer :: n1,np,n,irtn,k,i,m,lth,idofieldcalc
      real*8 :: h,t,tt,xdum,zdum
      integer :: mprocs,myrank,mpierr,ierr

      real*8, dimension(6) :: av  ! fix (6) later !!!!!!

      interface
          subroutine eval(t,y,f)
            real*8 :: t
            real*8, dimension(:,:) :: y,f
          end subroutine
      end interface

      call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,mpierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,mpierr)
      n1=ubound(y,1)
      np=ubound(y,2)
      if(myrank.eq.0)write(6,*)'(rk4) n1,np=',n1,np
      allocate(yt(n1,np),f(n1,np),a(n1,np),b(n1,np),c(n1,np),d(n1,np))
      h=(tfin-tini)/nsteps
      t=tini
!     call moments(y,n1,np,t,20,24,25,26,9,0,0,0,0,irtn)
      do n=1,nsteps
        if(laststoredrefpart.eq.-1)then
          if(t+h.ge.tauhistmin)then
            call afterstep(t,y) !store the initial values in the history arrays
          endif
        endif
        if(myrank.eq.0.and.mod(n,1000).eq.0)write(6,*)'STARTING STEP ',n
        call eval(t,y,f)
        a(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+0.5d0*a(1:n1,1:np)
        tt=t+0.5d0*h
        call eval(tt,yt,f)
        b(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+0.5d0*b(1:n1,1:np)
        tt=t+0.5d0*h
        call eval(tt,yt,f)
        c(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+c(1:n1,1:np)
        tt=t+h
        call eval(tt,yt,f)
        d(1:n1,1:np)=h*f(1:n1,1:np)
        y(1:n1,1:np)=y(1:n1,1:np)+(a(1:n1,1:np)+2.d0*b(1:n1,1:np)+2.d0*c(1:n1,1:np)+d(1:n1,1:np))/6.d0
        t=tini+n*h
!       call moments(y,n1,np,t,20,24,25,26,9,0,0,0,0,irtn)
        if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y) !store latest values in the history arrays

!       at this point the particle history has been stored up to time t.
!       getfieldsatpoint() can be called to compute the LW fields.

        call getcentroid(y,av,n1,np) !do this to compute the LW fields when
        call testtoperformfieldcalc(av,n,idofieldcalc)
        if(idofieldcalc.eq.1)then
          if(myrank.eq.0)write(6,*)'Step ',n,' complete; computing fields using history values'
          if(myrank.eq.0)write(6,*)'refptcltime(0)=',refptcltime(0)
          if(myrank.eq.0)write(6,*)'refptcltime(',laststoredrefpart,')=',refptcltime(laststoredrefpart)
          call getfield_atallpts(t,av,np,n)
          if(myrank.eq.0)write(6,*)'Done computing fields.'
          if(istopafter1fieldcalc.eq.1)exit !test to exit from loop
        endif
      enddo
      return
      end

      subroutine adams11(y,tini,tfin,nsteps,eval)
      use mpi
      use lienwiech, only : getfield_atallpts
      use putils, only : pwrite
      use diagnostics, only : getcentroid,moments
      use lwprobcons, only : tauhistmin,tauhistmax,refptcltime,laststoredrefpart,testtoperformfieldcalc,istopafter1fieldcalc
      implicit none
      real*8 :: tini,tfin
      integer :: nsteps
      real*8, dimension(:,:) :: y
      real*8, allocatable, dimension(:,:) :: yp,yc,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11
      real*8, dimension(10) :: am,bm
      real*8, parameter, dimension(10) :: a(1:10)= &
     &  (/57281.d0,-583435.d0,2687864.d0,-7394032.d0,13510082.d0,-17283646.d0,16002320.d0,-11271304.d0,9449717.d0,2082753.d0/)
      real*8, parameter, dimension(10) :: b(1:10)= &
     &  (/-2082753.d0,20884811.d0,-94307320.d0,252618224.d0,-444772162.d0,538363838.d0,-454661776.d0,265932680.d0,-104995189.d0,&
     &     30277247.d0/)
      integer :: n1,np,nsa,i,j,idofieldcalc
      real*8 :: h,t,tint,hdiv
      integer :: mprocs,myrank,mpierr,irtnerr

      real*8, dimension(6) :: av  ! fix (6) later !!!!!!

      interface
          subroutine eval(t,y,f)
            real*8 :: t
            real*8, dimension(:,:) :: y,f
          end subroutine
      end interface

      call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,mpierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,mpierr)

      n1=ubound(y,1)
      np=ubound(y,2)
      allocate(yp(n1,np),yc(n1,np))
      allocate(f1(n1,np),f2(n1,np),f3(n1,np),f4(n1,np),f5(n1,np),f6(n1,np),f7(n1,np),f8(n1,np),f9(n1,np),f10(n1,np),f11(n1,np))

      call getcentroid(y,av,n1,np)
      if(myrank.eq.0)write(6,*)'Entered Adams routine...'

      h=(tfin-tini)/nsteps
      t=tini
      if(laststoredrefpart.eq.-1)then
        if(t+h.ge.tauhistmin)then
          call afterstep(t,y) !store the initial values in the history arrays
        endif
      endif
      call eval(t,y,f1)
      call rk4cpy(y,tini,tini+h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f2)
      call rk4cpy(y,tini+h,tini+2.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f3)
      call rk4cpy(y,tini+2.d0*h,tini+3.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f4)
      call rk4cpy(y,tini+3.d0*h,tini+4.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f5)
      call rk4cpy(y,tini+4.d0*h,tini+5.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f6)
      call rk4cpy(y,tini+5.d0*h,tini+6.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f7)
      call rk4cpy(y,tini+6.d0*h,tini+7.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f8)
      call rk4cpy(y,tini+7.d0*h,tini+8.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f9)
      call rk4cpy(y,tini+8.d0*h,tini+9.d0*h,10,eval)
      t=t+h; if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y)
      call eval(t,y,f10)
      nsa=nsteps-9
      hdiv=h/7257600.0d+00

      am(1:10)=hdiv*a(1:10)
      bm(1:10)=hdiv*b(1:10)
      tint=t
      call moments(y,n1,np,tint,20,24,25,26,8,0,0,0,0,irtnerr)

      do 100 i=1,nsa
      yp(1:n1,1:np)=y(1:n1,1:np)+bm(1)*f1(1:n1,1:np)+bm(2)*f2(1:n1,1:np)+bm(3)*f3(1:n1,1:np) &
     & +bm(4)*f4(1:n1,1:np)+bm(5)*f5(1:n1,1:np)+bm(6)*f6(1:n1,1:np)+bm(7)*f7(1:n1,1:np) &
     & +bm(8)*f8(1:n1,1:np)+bm(9)*f9(1:n1,1:np)+bm(10)*f10(1:n1,1:np)
      call eval(t+h,yp,f11)
      yp(1:n1,1:np)=y(1:n1,1:np)+am(1)*f2(1:n1,1:np)+am(2)*f3(1:n1,1:np)+am(3)*f4(1:n1,1:np)+am(4)*f5(1:n1,1:np)+ &
     &                   am(5)*f6(1:n1,1:np)+am(6)*f7(1:n1,1:np)+am(7)*f8(1:n1,1:np)+am(8)*f9(1:n1,1:np)+am(9)*f10(1:n1,1:np)
      yc(1:n1,1:np)=yp(1:n1,1:np)+am(10)*f11(1:n1,1:np)
      call eval(t+h,yc,f11)
      y(1:n1,1:np)=yp(1:n1,1:np)+am(10)*f11(1:n1,1:np)
      f1(1:n1,1:np)=f2(1:n1,1:np)
      f2(1:n1,1:np)=f3(1:n1,1:np)
      f3(1:n1,1:np)=f4(1:n1,1:np)
      f4(1:n1,1:np)=f5(1:n1,1:np)
      f5(1:n1,1:np)=f6(1:n1,1:np)
      f6(1:n1,1:np)=f7(1:n1,1:np)
      f7(1:n1,1:np)=f8(1:n1,1:np)
      f8(1:n1,1:np)=f9(1:n1,1:np)
      f9(1:n1,1:np)=f10(1:n1,1:np)
      f10(1:n1,1:np)=f11(1:n1,1:np)
      t=tint+i*h
      call moments(y,n1,np,t,20,24,25,26,8,0,0,0,0,irtnerr)
      if(t.ge.tauhistmin .and. t.le.tauhistmax)call afterstep(t,y) !store latest values in the history arrays
      call getcentroid(y,av,n1,np)
      call testtoperformfieldcalc(av,i+9,idofieldcalc)
      if(idofieldcalc.eq.1)then
        if(myrank.eq.0)write(6,*)'Step ',i+9,' complete; computing fields at t=',t
        call getfield_atallpts(t,av,np,i+9)
        if(myrank.eq.0)write(6,*)'Done computing fields.'
        if(istopafter1fieldcalc.eq.1)exit !test to exit from loop
      endif
  100 continue
      return
      end

      subroutine afterstep(t,y)
      use mpi
      use lwprobcons, only : refptclhist,refptcltime,laststoredrefpart,nhistpoints
      use evalroutines, only : evalt1
      implicit none
      real*8 :: t
      real*8, dimension(:,:) :: y
      real*8, dimension(6) :: av
      integer :: n1,np,n,nstep
      integer :: myrank,mpierr

      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,mpierr)

      n1=ubound(y,1)
      np=ubound(y,2)
      nstep=laststoredrefpart+1 !assumes initilized to -1 in module  !!!!!!

      if(nstep.eq.nhistpoints+1)then
        if(myrank.eq.0)write(6,*)'shifting the arrays of stored quantities'
        do n=1,nhistpoints
          refptcltime(n-1)=refptcltime(n)
          refptclhist(1:6,n-1,1:np)=refptclhist(1:6,n,1:np)
        enddo
        nstep=nhistpoints
      endif

      refptcltime(nstep)=t

      do n=1,np
      refptclhist(1:6,nstep,n)=y(1:6,n)
      enddo
      if(myrank.eq.0)write(91,'(7(1pe16.9,1x))')refptcltime(nstep),refptclhist(1:6,nstep,1)
      laststoredrefpart=nstep
      return
      end

      subroutine rk4cpy(y,tini,tfin,nsteps,eval)
      use mpi
      implicit none
      real*8, dimension(:,:) :: y
      real*8 :: tini,tfin
      integer :: nsteps,mlo,mhi,mlomin,mhimax,gblmlomin,gblmhimax
      real*8, allocatable, dimension(:,:) :: yt,f,a,b,c,d
      integer :: n1,np,n,irtn,k,i,m,lth
      real*8 :: h,t,tt,xdum,zdum
      integer :: mprocs,myrank,mpierr,ierr

      real*8, dimension(6) :: av  ! fix (6) later !!!!!!

      interface
          subroutine eval(t,y,f)
            real*8 :: t
            real*8, dimension(:,:) :: y,f
          end subroutine
      end interface

      call MPI_COMM_SIZE(MPI_COMM_WORLD,mprocs,mpierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,mpierr)
      n1=ubound(y,1)
      np=ubound(y,2)
      allocate(yt(n1,np),f(n1,np),a(n1,np),b(n1,np),c(n1,np),d(n1,np))
      h=(tfin-tini)/nsteps
      t=tini
      do n=1,nsteps
        call eval(t,y,f)
        a(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+0.5d0*a(1:n1,1:np)
        tt=t+0.5d0*h
        call eval(tt,yt,f)
        b(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+0.5d0*b(1:n1,1:np)
        tt=t+0.5d0*h
        call eval(tt,yt,f)
        c(1:n1,1:np)=h*f(1:n1,1:np)
        yt(1:n1,1:np)=y(1:n1,1:np)+c(1:n1,1:np)
        tt=t+h
        call eval(tt,yt,f)
        d(1:n1,1:np)=h*f(1:n1,1:np)
        y(1:n1,1:np)=y(1:n1,1:np)+(a(1:n1,1:np)+2.d0*b(1:n1,1:np)+2.d0*c(1:n1,1:np)+d(1:n1,1:np))/6.d0
        t=tini+n*h
      enddo
      return
      end

      subroutine tdriftexact(tauini,taufin,y)
      implicit none
      real*8 :: tauini,taufin
      real*8, dimension(:,:) :: y
      real*8 :: gam2,gam
      integer :: n1,np,n
      real*8, parameter :: clite=299792458.d0
      n1=ubound(y,1)
      np=ubound(y,2)
      do n=1,np
        gam2=1.d0+y(2,n)**2+y(4,n)**2+y(6,n)**2
        gam=sqrt(gam2)
        y(1,n)=y(1,n)+(taufin-tauini)*y(2,n)*clite/gam
        y(3,n)=y(3,n)+(taufin-tauini)*y(4,n)*clite/gam
        y(5,n)=y(5,n)+(taufin-tauini)*y(6,n)*clite/gam
      enddo
      return
      end

end module


! { dg-do run }
!
! Tests the fix for PR68216
!
! Reported on clf: https://groups.google.com/forum/#!topic/comp.lang.fortran/eWQTKfqKLZc
!
PROGRAM hello
!
! This is based on the first testcase, from Francisco (Ayyy LMAO). Original
! lines are commented out. The second testcase from this thread is acalled
! at the end of the program.
!
    IMPLICIT NONE

    CHARACTER(LEN=:),DIMENSION(:),ALLOCATABLE :: array_lineas
    CHARACTER(LEN=:),DIMENSION(:),ALLOCATABLE :: array_copia
    character (3), dimension (2) :: array_fijo = ["abc","def"]
    character (100) :: buffer
    INTEGER :: largo , cant_lineas , i

    write (buffer, "(2a3)") array_fijo

!    WRITE(*,*) ' Escriba un numero para el largo de cada linea'
!    READ(*,*) largo
    largo = LEN (array_fijo)

!    WRITE(*,*) ' Escriba la cantidad de lineas'
!    READ(*,*) cant_lineas
    cant_lineas = size (array_fijo, 1)

    ALLOCATE(CHARACTER(LEN=largo) :: array_lineas(cant_lineas))

!    WRITE(*,*) 'Escriba el array', len(array_lineas), size(array_lineas)
    READ(buffer,"(2a3)") (array_lineas(i),i=1,cant_lineas)

!    WRITE(*,*) 'Array guardado: '
!    DO i=1,cant_lineas
!    WRITE(*,*) array_lineas(i)
!    ENDDO
     if (any (array_lineas .ne. array_fijo)) STOP 1

! The following are additional tests beyond that of the original.
!
! Check that allocation with source = another deferred length is OK
     allocate (array_copia, source = array_lineas)
     if (any (array_copia .ne. array_fijo)) STOP 2
     deallocate (array_lineas, array_copia)

! Check that allocation with source = a non-deferred length is OK
     allocate (array_lineas, source = array_fijo)
     if (any (array_lineas .ne. array_fijo)) STOP 3
     deallocate (array_lineas)

! Check that allocation with MOLD = a non-deferred length is OK
     allocate (array_copia, mold = [array_fijo(:)(1:2), array_fijo(:)(1:2)])
     if (size (array_copia, 1) .ne. 4) STOP 4
     if (LEN (array_copia, 1) .ne. 2) STOP 5

! Check that allocation with MOLD = another deferred length is OK
     allocate (array_lineas, mold = array_copia)
     if (size (array_copia, 1) .ne. 4) STOP 6
     if (LEN (array_copia, 1) .ne. 2) STOP 7
     deallocate (array_lineas, array_copia)

!    READ(*,*)
     call testdefchar
contains
     subroutine testdefchar
!
! This is the testcase in the above thread from Blokbuster
!
          implicit none
          character(:), allocatable :: test(:)

          allocate(character(3) :: test(2))
          test(1) = 'abc'
          test(2) = 'def'
          if (any (test .ne. ['abc', 'def'])) STOP 8

          test = ['aa','bb','cc']
          if (any (test .ne. ['aa', 'bb', 'cc'])) STOP 9

     end subroutine testdefchar

END PROGRAM


text:
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/yoj/arc/CCTM/src/vdiff/acm2/matrix.F,v 1.5 2011/10/21 16:11:45 yoj Exp $

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE SA_MATRIX1 ( KL, A, B, E, D, X )

C Rather than solving the ACM2 banded tridiagonal matrix using LU decomposition,
C it is much faster to split the solution into the ACM1 convective solver followed
C by the tridiagonal solver
C MATRIX1 is the ACM1 solver. When the PBL is convective, this solver is called
C followed by TRI. If not convective, only TRI is called.

C-- ACM1 Matrix is in this form (there is no subdiagonal:
C   B1 E2                     <- note E2 (flux from layer above), not E1
C   A2 B2 E3
C   A3    B3 E4
C   A4       B4 E5
C   A5          B5 E6
C   A6             B6

      USE VGRD_DEFN           ! vertical layer specifications
      USE SA_DEFN
!      USE CGRID_SPCS          ! CGRID mechanism species
      USE UTILIO_DEFN

      IMPLICIT NONE

!C Includes:
!      INTEGER, SAVE :: N_SPC_DIFF    ! global diffusion species

C Arguments:
      INTEGER, INTENT( IN )  :: KL         ! CBL sigma height
      REAL,    INTENT( IN )  :: A( : )     ! matrix column one
      REAL,    INTENT( IN )  :: B( : )     ! diagonal
      REAL,    INTENT( IN )  :: E( : )     ! superdiagonal
      REAL,    INTENT( IN )  :: D( :,: )   ! R.H.S
      REAL,    INTENT( OUT ) :: X( :,: )   ! returned solution

C Locals:
      REAL, ALLOCATABLE, SAVE :: BETA( : )
!     REAL :: BETA( N_SPCTAG )
      REAL  ALPHA, GAMA

      INTEGER L, V, IOS

      CHARACTER( 120 ) :: XMSG = ' '
      LOGICAL :: FIRSTIME = .TRUE.

C-----------------------------------------------------------------------

       IF ( FIRSTIME ) THEN
          FIRSTIME = .FALSE.
!         N_SPC_DIFF = N_GC_TRNS + N_AE_TRNS + N_NR_TRNS + N_TR_DIFF
          ALLOCATE ( BETA( N_SPCTAG ), STAT = IOS )
          IF ( IOS .NE. 0 ) THEN
             XMSG = 'Failure allocating BETA'
             CALL M3EXIT( 'SA_MATRIX', 0, 0, XMSG, XSTAT1 )
          END IF
       END IF   ! FIRSTIME

C-- ACM1 matrix solver

      DO V = 1, N_SPCTAG
         BETA( V ) = D( V,1 )
      END DO
      GAMA = B( 1 )
      ALPHA = 1.0

      DO L = 2, KL
         ALPHA = -ALPHA * E( L ) / B( L )
         DO V = 1, N_SPCTAG
            BETA( V ) = ALPHA * D( V,L ) + BETA( V )
         END DO
         GAMA = GAMA + ALPHA * A( L )
      END DO

      DO V = 1, N_SPCTAG
         X( V,1 )  = BETA( V ) / GAMA
         X( V,KL ) = ( D( V,KL ) - A( KL ) * X( V,1 ) ) / B( KL )
      END DO

C-- Back sub for Ux=y

      DO L = KL-1, 2, -1
         DO V = 1, N_SPCTAG
            X( V,L ) = ( D( V,L ) - A( L ) * X( V,1 ) - E( L+1 ) * X( V,L+1 ) ) / B( L )
         END DO
      END DO

      RETURN
      END

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
!!----
!!----
!!----
SubModule (Test_2) SubTest3
   Contains

   module subroutine set_m ( this, m )
     class(type_1), intent(inout) :: this
     integer,       intent(in) :: m

     this%n=m*3
  end subroutine set_m

End SubModule SubTest3

max_stars_repo_name: bdbaddog/scons-bugswork

max_issues_repo_name: bdbaddog/scons-bugswork

max_forks_repo_name: bdbaddog/scons-bugswork


text:
! { dg-do compile }
! we used to save the wrong components of a gfc_expr describing a
! substring of a constant string.  This yielded a segfault on
! translating the expressions read from the module.
module m
  character (*), parameter :: a = "AABBCC"(1:4)
end module m

use m
character(4) :: b
b = a
end

! { dg-final { cleanup-modules "m" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      subroutine qqb_WZjj(p,msq)
C-----Matrix squared for
C-----p(p1)+p(p2)-->W(l(p3)+a(p4))+Z(l(p5)+a(p6)+p(p7)+p(p8)
C-----Author: R.K.Ellis February 2013
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      include 'nwz.f'
      double precision p(mxpart,4),msq(-nf:nf,-nf:nf)
     & ,WZbbmsq,WZccmsq,WZggmsq,WZddidmsq,WZuuidmsq,
     & ms(-nf:nf,-nf:nf,-nf:nf,-nf:nf)

      integer i,j,k,l,ba,ac,sa,ua,da,gg,dq,uq,sq,cq,bq
      parameter(ba=-5,ac=-4,sa=-3,ua=-2,da=-1,
     &     gg=0,bq=+5,cq=+4,sq=+3,uq=+2,dq=+1)

      msq(:,:)=0d0
      ms(:,:,:,:)=0d0
c      open(unit=66,file='madpoints.f',status='unknown')
c      call writeformadgraph(p)
c      write (*,*)
c      write (*,*) " Phase space point:"
c      write (*,*)
c      write (*,*) "-----------------------------------------------------------------------------"
c      write (*,*)  "n        E             px             py              pz               m "
c      do i=1,8
c      if (i .lt. 3)
c     &    write (*,'(i2,1x,5e15.7)') i, -P(i,4),-P(i,1),-P(i,2),-P(i,3),
c     & p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2
c      if (i .ge. 3)
c     &    write (*,'(i2,1x,5e15.7)') i, P(i,4),P(i,1),P(i,2),P(i,3),
cc     . dsqrt(dabs(DOT(p(0,i),p(0,i))))
c     & p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2
c      enddo
c      write (*,*) "-----------------------------------------------------------------------------"
c      write (*,*) "-----------------------------------------------------------------------------"

      call spinoru(8,p,za,zb)

      if (nwz .eq. 1) then

C***********************************************************
C-----Identity tables antiquark-quark
c--- d~  d > d  u~ and s  s~ > s  c~
      ms(da,dq,dq,ua)=WZddidmsq(8,7,3,4,5,6,1,2)
      ms(sa,sq,sq,ac)=ms(da,dq,dq,ua)
c--- u~  u > d  u~ and c~  c > s  c~
      ms(ua,uq,dq,ua)=WZuuidmsq(2,7,3,4,5,6,1,8)
      ms(ac,cq,sq,ac)=ms(ua,uq,dq,ua)
c--- d~  u > d  d~ and s~  c > s  s~
      ms(da,uq,dq,da)=WZddidmsq(2,1,3,4,5,6,7,8)
      ms(sa,cq,sq,sa)=ms(da,uq,dq,da)
c--- s~  s > d  u~ and d~  d > s  c~ and b~  b > d  u~ and b~  b > s  c~
      ms(da,dq,sq,ac)=WZbbmsq(8,7,3,4,5,6,1,2)
      ms(sa,sq,dq,ua)=ms(da,dq,sq,ac)
      ms(ba,bq,dq,ua)=ms(da,dq,sq,ac)
      ms(ba,bq,sq,ac)=ms(da,dq,sq,ac)
c--- d~  u > s  s~ and s~  c > d  d~ and s~  c > b  b~ and d~  u > b  b~
      ms(da,uq,sq,sa)=WZbbmsq(2,1,3,4,5,6,7,8)
      ms(sa,cq,dq,da)=ms(da,uq,sq,sa)
      ms(sa,cq,bq,ba)=ms(da,uq,sq,sa)
      ms(da,uq,bq,ba)=ms(da,uq,sq,sa)
c--- d~  u > u  u~ and s~  c > c  c~
      ms(da,uq,uq,ua)=WZuuidmsq(2,1,3,4,5,6,7,8)
      ms(sa,cq,cq,ac)=ms(da,uq,uq,ua)
c--- s~  u  > d  s~
      ms(sa,uq,dq,sa)=WZbbmsq(2,7,3,4,5,6,1,8)
c--- b~  u   > d  b~
      ms(ba,uq,dq,ba)=ms(sa,uq,dq,sa)
c--- d~  c  > s  d~
      ms(da,cq,sq,da)=ms(sa,uq,dq,sa)
c--- s~  d > d  c~
      ms(sa,dq,dq,ac)=WZbbmsq(8,1,3,4,5,6,7,2)
c--- d~  s  > s  u~
      ms(da,sq,sq,ua)=ms(sa,dq,dq,ac)
c--- d~  b  > b  u~
      ms(da,bq,bq,ua)=ms(sa,dq,dq,ac)
c--- s~  b  > b  c~
      ms(sa,bq,bq,ac)=ms(sa,dq,dq,ac)
c--- u~  u > s  c~ and c~  c > d  u~
      ms(ua,uq,sq,ac)=WZccmsq(8,7,3,4,5,6,1,2)
      ms(ac,cq,dq,ua)=ms(ua,uq,sq,ac)
c--- d~  u  > c  c~
      ms(da,uq,cq,ac)=WZccmsq(2,1,3,4,5,6,7,8)
c--- s~  c  > u  u~
      ms(sa,cq,uq,ua)=ms(da,uq,cq,ac)
c--- c~  u  > d  c~
      ms(ac,uq,dq,ac)=WZccmsq(2,7,3,4,5,6,1,8)
c--- u~  c  > s  u~
      ms(ua,cq,sq,ua)=ms(ac,uq,dq,ac)
c--- s~ u  > u  c~
      ms(sa,uq,uq,ac)=WZccmsq(8,1,3,4,5,6,7,2)
c--- d~  c  > c  u~
      ms(da,cq,cq,ua)=ms(sa,uq,uq,ac)
c--- b~  c > s  b~
      ms(ba,cq,sq,ba)=ms(sa,uq,dq,sa)

C***********************************************************
C-----Identity tables quark-antiquark
c--- d  d~ > d  u~ and s  s~ > s  c~
      ms(dq,da,dq,ua)=WZddidmsq(8,7,3,4,5,6,2,1)
      ms(sq,sa,sq,ac)=ms(dq,da,dq,ua)
c--- u  d~ > d  d~ and c  s~ > s  s~
      ms(uq,da,dq,da)=WZddidmsq(1,2,3,4,5,6,7,8)
      ms(cq,sa,sq,sa)=ms(uq,da,dq,da)
c--- u  d~ > u  u~ and c  s~ > c  c~
      ms(uq,da,uq,ua)=WZuuidmsq(1,2,3,4,5,6,7,8)
      ms(cq,sa,cq,ac)=ms(uq,da,uq,ua)
c--- u  u~ > d  u~ and c  c~ > s  c~
      ms(uq,ua,dq,ua)=WZuuidmsq(1,7,3,4,5,6,2,8)
      ms(cq,ac,sq,ac)=ms(uq,ua,dq,ua)
c--- s  s~ > d  u~ and d  d~ > s  c~ and b  b~ > d  u~ and b  b~ > s  c~
      ms(dq,da,sq,ac)=WZbbmsq(8,7,3,4,5,6,2,1)
      ms(sq,sa,dq,ua)=ms(dq,da,sq,ac)
      ms(bq,ba,dq,ua)=ms(dq,da,sq,ac)
      ms(bq,ba,sq,ac)=ms(dq,da,sq,ac)
c--- u  d~ > s  s~ and c  s~ > d  d~ and c  s~ > b  b~ and u  d~ > b  b~
      ms(uq,da,sq,sa)=WZbbmsq(1,2,3,4,5,6,7,8)
      ms(cq,sa,dq,da)=ms(uq,da,sq,sa)
      ms(cq,sa,bq,ba)=ms(uq,da,sq,sa)
      ms(uq,da,bq,ba)=ms(uq,da,sq,sa)
c--- u  s~ > d  s~
      ms(uq,sa,dq,sa)=WZbbmsq(1,7,3,4,5,6,2,8)
c--- u  b~ > d  b~
      ms(uq,ba,dq,ba)=ms(uq,sa,dq,sa)
c--- c  d~ > s  d~
      ms(cq,da,sq,da)=ms(uq,sa,dq,sa)
c--- d  s~ > d  c~
      ms(dq,sa,dq,ac)=WZbbmsq(8,2,3,4,5,6,7,1)
c--- s  d~ > s  u~
      ms(sq,da,sq,ua)=ms(dq,sa,dq,ac)
c--- b  d~ > b  u~
      ms(bq,da,bq,ua)=ms(dq,sa,dq,ac)
c--- b  s~ > b  c~
      ms(bq,sa,bq,ac)=ms(dq,sa,dq,ac)
c--- u  u~ > s  c~ and c  c~ > d  u~
      ms(uq,ua,sq,ac)=WZccmsq(8,7,3,4,5,6,2,1)
      ms(cq,ac,dq,ua)=ms(uq,ua,sq,ac)
c--- u  d~ > c  c~
      ms(uq,da,cq,ac)=WZccmsq(1,2,3,4,5,6,7,8)
c--- c  s~ > u  u~
      ms(cq,sa,uq,ua)=ms(uq,da,cq,ac)
c--- u  c~ > d  c~
      ms(uq,ac,dq,ac)=WZccmsq(1,7,3,4,5,6,2,8)
c--- c  u~ > s  u~
      ms(cq,ua,sq,ua)=ms(uq,ac,dq,ac)
c--- u  s~ > u  c~
      ms(uq,sa,uq,ac)=WZccmsq(8,2,3,4,5,6,7,1)
c--- c  d~ > c  u~
      ms(cq,da,cq,ua)=ms(uq,sa,uq,ac)
c--- c  b~ > s  b~
      ms(cq,ba,sq,ba)=ms(uq,sa,dq,sa)
C***********************************************************
C-----Identity tables antiquark-antiquark
c--- u~ d~ > u~ u~
      ms(ua,da,ua,ua)=half*WZuuidmsq(7,2,3,4,5,6,1,8)
c--- c~ s~ > c~ c~
      ms(ac,sa,ac,ac)=ms(ua,da,ua,ua)
c--- d~ u~ > u~ u~
      ms(da,ua,ua,ua)=half*WZuuidmsq(7,1,3,4,5,6,2,8)
c--- s~ c~ > c~ c~
      ms(sa,ac,ac,ac)=ms(da,ua,ua,ua)
c--- d~ d~ > d~ u~
      ms(da,da,da,ua)=WZddidmsq(8,1,3,4,5,6,2,7)
c--- s~ s~ > s~ c~
      ms(sa,sa,sa,ac)=ms(da,da,da,ua)
c--- d~ c~ > u~ c~
      ms(da,ac,ua,ac)=WZccmsq(7,1,3,4,5,6,2,8)
c--- s~ u~ > c~ u~
      ms(sa,ua,ac,ua)=ms(da,ac,ua,ac)
c--- u~ s~ > c~  u~
      ms(ua,sa,ac,ua)=WZccmsq(7,2,3,4,5,6,1,8)
c--- c~ d~ > u~ c~
      ms(ac,da,ua,ac)=ms(ua,sa,ac,ua)
c--- s~ d~ u~ s~ and s~ d~ > c~ d~
c--- d~ s~ > u~ s~ and d~ s~ > c~ d~
      ms(sa,da,ua,sa)=WZbbmsq(7,2,3,4,5,6,1,8)
      ms(sa,da,ac,da)=WZbbmsq(7,1,3,4,5,6,2,8)
      ms(da,sa,ua,sa)=ms(sa,da,ac,da)
      ms(da,sa,ac,da)=ms(sa,da,ua,sa)
c--- d~ b~ > u~ b~
      ms(da,ba,ua,ba)=ms(da,sa,ua,sa)
c--- s~ b~ > c~ b~
      ms(sa,ba,ac,ba)=ms(da,sa,ua,sa)
c--- b~ d~ > u~ b~
      ms(ba,da,ua,ba)=ms(sa,da,ua,sa)
c--- b~ s~ > c~ b~
      ms(ba,sa,ac,ba)=ms(sa,da,ua,sa)
C***********************************************************
C-----Identity tables quark-quark
c--- u  d  > d  d
      ms(uq,dq,dq,dq)=half*WZddidmsq(1,7,3,4,5,6,8,2)
c--- c  s  > s  s
      ms(cq,sq,sq,sq)=ms(uq,dq,dq,dq)
c--- d  u  > d  d
      ms(dq,uq,dq,dq)=half*WZddidmsq(2,7,3,4,5,6,8,1)
c--- s  c  > s  s
      ms(sq,cq,sq,sq)=ms(dq,uq,dq,dq)
c--- u  u  > d  u
      ms(uq,uq,dq,uq)=WZuuidmsq(1,7,3,4,5,6,8,2)
c--- c  c  > s  c
      ms(cq,cq,sq,cq)=ms(uq,uq,dq,uq)
c--- u  b  > d  b
      ms(uq,bq,dq,bq)=WZbbmsq(1,7,3,4,5,6,8,2)
c--- u  s  > d  s
      ms(uq,sq,dq,sq)=ms(uq,bq,dq,bq)
c--- c  b  > s  b
      ms(cq,bq,sq,bq)=ms(uq,bq,dq,bq)
c--- c  d  > s  d
      ms(cq,dq,sq,dq)=ms(uq,bq,dq,bq)
c--- b  u  > d  b
      ms(bq,uq,dq,bq)=WZbbmsq(2,7,3,4,5,6,8,1)
c--- s  u  > d  s
      ms(sq,uq,dq,sq)=ms(bq,uq,dq,bq)
c--- b  c  > s  b
      ms(bq,cq,sq,bq)=ms(bq,uq,dq,bq)
c--- d  c  > s  d
      ms(dq,cq,sq,dq)=ms(bq,uq,dq,bq)
c--- u  c  > d  c
      ms(uq,cq,dq,cq)=WZccmsq(1,7,3,4,5,6,8,2)
c--- u  c  > u  s
      ms(uq,cq,uq,sq)=WZccmsq(2,8,3,4,5,6,7,1)
c--- c  u  > c  d
      ms(cq,uq,cq,dq)=ms(uq,cq,uq,sq)
c--- c  u  > s  u
      ms(cq,uq,sq,uq)=ms(uq,cq,dq,cq)

C***********************2q-2g-processes*********************************
c--- u   d~ > g  g
      ms(uq,da,gg,gg)=aveqq*half*WZggmsq(1,2,3,4,5,6,7,8)
c--- c   s~ > g  g
      ms(cq,sa,gg,gg)=ms(uq,da,gg,gg)
c--- d~  u  > g  g
      ms(da,uq,gg,gg)=aveqq*half*WZggmsq(2,1,3,4,5,6,7,8)
c--- s~  c  > g  g
      ms(sa,cq,gg,gg)=ms(da,uq,gg,gg)

c--- g  g  > u~  d
      ms(gg,gg,ua,dq)=avegg*WZggmsq(7,8,3,4,5,6,2,1)
c--- g  g  > c~  s
      ms(gg,gg,ac,sq)=ms(gg,gg,ua,dq)

c--- u  g  > d g
      ms(uq,gg,dq,gg)=aveqg*WZggmsq(1,7,3,4,5,6,2,8)
c--- c  g  > s g
      ms(cq,gg,sq,gg)=ms(uq,gg,dq,gg)
c--- g  u > d g
      ms(gg,uq,dq,gg)=aveqg*WZggmsq(2,7,3,4,5,6,1,8)
c--- g  c > s g
      ms(gg,cq,sq,gg)=ms(gg,uq,dq,gg)

c--- g  d~  > u~  g
      ms(gg,da,ua,gg)=aveqg*WZggmsq(7,2,3,4,5,6,1,8)
c--- g  s~  > c~  g
      ms(gg,sa,ac,gg)=ms(gg,da,ua,gg)
c--- d~ g  > u~  g
      ms(da,gg,ua,gg)=aveqg*WZggmsq(7,1,3,4,5,6,2,8)
c--- s~ g > c~  g
      ms(sa,gg,ac,gg)=ms(da,gg,ua,gg)

c      call madcheckp(ms)
c      write(6,*) 'qq'
c      include 'qq+'
c      write(6,*) 'aa'
c      include 'aa+'
c      write(6,*) 'qa'
c      include 'qa+'
c      write(6,*) 'aq'
c      include 'aq+'

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

      elseif (nwz .eq. -1) then

C***********************************************************
C-----Identity table antiquark-quark
c--- u~  u > u  d~ and c  c~ > c  s~
      ms(ua,uq,uq,da)=WZuuidmsq(8,7,3,4,5,6,1,2)
      ms(ac,cq,cq,sa)=ms(ua,uq,uq,da)
c--- u~  d > u  u~ and c~  s > c  c~
      ms(ua,dq,uq,ua)=WZuuidmsq(2,1,3,4,5,6,7,8)
      ms(ac,sq,cq,ac)=ms(ua,dq,uq,ua)
c--- d~  d > u  d~ and s~  s > c  s~
      ms(da,dq,uq,da)=WZddidmsq(2,7,3,4,5,6,1,8)
      ms(sa,sq,cq,sa)=ms(da,dq,uq,da)
c--- u~  d > d  d~ and c~  s > s  s~
      ms(ua,dq,dq,da)=WZddidmsq(2,1,3,4,5,6,7,8)
      ms(ac,sq,sq,sa)=ms(ua,dq,dq,da)

c--- d~  d > c  s~ and s~  s > u  d~
      ms(da,dq,cq,sa)=WZbbmsq(8,7,3,4,5,6,1,2)
      ms(sa,sq,uq,da)=ms(da,dq,cq,sa)
      ms(ba,bq,uq,da)=ms(da,dq,cq,sa)
      ms(ba,bq,cq,sa)=ms(da,dq,cq,sa)
c--- u~  d  > s  s~ and c~  s  > d  d~ and c~  s > b  b~ and u~  d > b  b~
      ms(ua,dq,sq,sa)=WZbbmsq(2,1,3,4,5,6,7,8)
      ms(ac,sq,dq,da)=ms(ua,dq,sq,sa)
      ms(ua,dq,bq,ba)=ms(ua,dq,sq,sa)
      ms(ac,sq,bq,ba)=ms(ua,dq,sq,sa)
c--- s~  d  > u  s~
      ms(sa,dq,uq,sa)=WZbbmsq(2,7,3,4,5,6,1,8)
c--- d~  s  > c  d~
      ms(da,sq,cq,da)=ms(sa,dq,uq,sa)
c--- b~  s > c  b~
      ms(ba,sq,cq,ba)=ms(sa,dq,uq,sa)
c--- b~  d   > u  b~
      ms(ba,dq,uq,ba)=ms(sa,dq,uq,sa)
c--- c~ d  > d  s~
      ms(ac,dq,dq,sa)=WZbbmsq(8,1,3,4,5,6,7,2)
c--- u~  s  > s  d~
      ms(ua,sq,sq,da)=ms(ac,dq,dq,sa)
c--- u~  b  > b  d~
      ms(ua,bq,bq,da)=ms(ua,sq,sq,da)
c--- c~  b  > b  s~
      ms(ac,bq,bq,sa)=ms(ua,sq,sq,da)
c--- c~  c > u  d~ and u~  u > c  s~ and b~  b > u  d~ and b~  b > c  s~
      ms(ua,uq,cq,sa)=WZccmsq(8,7,3,4,5,6,1,2)
      ms(ac,cq,uq,da)=ms(ua,uq,cq,sa)
c--- u~  d > c  c~ and c~  s > u  u~ and
      ms(ua,dq,cq,ac)=WZccmsq(2,1,3,4,5,6,7,8)
      ms(ac,sq,uq,ua)=ms(ua,dq,cq,ac)
      ms(ac,sq,bq,ba)=ms(ua,dq,cq,ac)
c--- c~  u > u  s~
      ms(ac,uq,uq,sa)=WZccmsq(8,1,3,4,5,6,7,2)
c--- u~  c  > c  d~
      ms(ua,cq,cq,da)=ms(ac,uq,uq,sa)
c--- c~  d  > u  c~
      ms(ac,dq,uq,ac)=WZccmsq(2,7,3,4,5,6,1,8)
c--- u~  s  > c  u~
      ms(ua,sq,cq,ua)=ms(ac,dq,uq,ac)

C***********************************************************
C-----Identity table antiquark-antiquark
c--- d~ u~ > d~ d~
      ms(da,ua,da,da)=half*WZddidmsq(7,2,3,4,5,6,1,8)
c--- s~ c~ > s~ s~
      ms(sa,ac,sa,sa)=ms(da,ua,da,da)
c--- u~ d~ > d~ d~
      ms(ua,da,da,da)=half*WZddidmsq(7,1,3,4,5,6,2,8)
c--- c~ s~ > s~ s~
      ms(ac,sa,sa,sa)=ms(ua,da,da,da)
c--- u~ u~ > u~ d~
      ms(ua,ua,ua,da)=WZuuidmsq(8,1,3,4,5,6,2,7)
c--- c~ c~ > c~ s~
      ms(ac,ac,ac,sa)=ms(ua,ua,ua,da)
c--- u~ s~ > d~ s~
      ms(ua,sa,da,sa)=WZbbmsq(7,1,3,4,5,6,2,8)
c--- c~ b~ > s~ b~
      ms(ac,ba,sa,ba)=ms(ua,sa,da,sa)
c--- u~ b~ > d~ b~
      ms(ua,ba,da,ba)=ms(ua,sa,da,sa)
c--- c~ d~ > s~ u~
      ms(ac,da,sa,da)=ms(ua,sa,da,sa)
c--- d~ c~ > s~  d~
      ms(da,ac,sa,da)=WZbbmsq(7,2,3,4,5,6,1,8)
c--- b~ u~ > d~ b~
      ms(ba,ua,da,ba)=ms(da,ac,sa,da)
c--- b~ c~ > s~ b~
      ms(ba,ac,sa,ba)=ms(da,ac,sa,da)
c--- s~ u~ > d~ s~
      ms(sa,ua,da,sa)=ms(da,ac,sa,da)
c--- c~ u~ > d~ c~ and c~ u~ > s~ u~
c--- u~ c~ > d~ c~ and u~ c~ > s~ u~
      ms(ac,ua,da,ac)=WZccmsq(7,2,3,4,5,6,1,8)
      ms(ac,ua,sa,ua)=WZccmsq(7,1,3,4,5,6,2,8)
      ms(ua,ac,da,ac)=ms(ac,ua,sa,ua)
      ms(ua,ac,sa,ua)=ms(ac,ua,da,ac)

C***********************************************************
C-----Identity tablec quark-quark
c--- d  u  > u  u
      ms(dq,uq,uq,uq)=half*WZuuidmsq(1,7,3,4,5,6,8,2)
c--- s  c  > c  c
      ms(sq,cq,cq,cq)=ms(dq,uq,uq,uq)
c--- u  d  > u  u
      ms(uq,dq,uq,uq)=half*WZuuidmsq(2,7,3,4,5,6,8,1)
c--- c  s  > c  c
      ms(cq,sq,cq,cq)=ms(uq,dq,uq,uq)
c--- d  d  > u  d
      ms(dq,dq,uq,dq)=WZddidmsq(1,7,3,4,5,6,8,2)
c--- s  s  > c  s
      ms(sq,sq,cq,sq)=ms(dq,dq,uq,dq)
c--- d  b  > u  b
      ms(dq,bq,uq,bq)=WZbbmsq(1,7,3,4,5,6,8,2)
c--- d  c  > u  c
      ms(dq,cq,uq,cq)=WZccmsq(1,7,3,4,5,6,8,2)
c--- s  b  > c  b
      ms(sq,bq,cq,bq)=ms(dq,bq,uq,bq)
c--- s  u  > c  u
      ms(sq,uq,cq,uq)=ms(dq,cq,uq,cq)
c--- b  d  > u  b
      ms(bq,dq,uq,bq)=WZbbmsq(2,7,3,4,5,6,8,1)
c--- b  s  > c  b
      ms(bq,sq,cq,bq)=ms(bq,dq,uq,bq)
c--- c  d  > u  c
      ms(cq,dq,uq,cq)=WZccmsq(2,7,3,4,5,6,8,1)
c--- u  s  > c  u
      ms(uq,sq,cq,uq)=ms(cq,dq,uq,cq)
c--- d  s  > u  s
      ms(dq,sq,uq,sq)=WZbbmsq(1,7,3,4,5,6,8,2)
c--- d  s  > d  s
      ms(dq,sq,dq,cq)=WZbbmsq(2,8,3,4,5,6,7,1)
c--- s  d  > s  u
      ms(sq,dq,sq,uq)=ms(dq,sq,dq,cq)
c--- s  d  > c  u
      ms(sq,dq,cq,dq)=ms(dq,sq,uq,sq)

C***********************************************************
C-----Identity tablec quark-antiquark
c--- u  u~ > u  d~ and c  c~ > c  s~
      ms(uq,ua,uq,da)=WZuuidmsq(8,7,3,4,5,6,2,1)
      ms(cq,ac,cq,sa)=ms(uq,ua,uq,da)
c--- d  d~ > u  d~ and s  s~ > c  s~
      ms(dq,da,uq,da)=WZddidmsq(1,7,3,4,5,6,2,8)
      ms(sq,sa,cq,sa)=ms(dq,da,uq,da)
c--- d  d~ > c  s~ and s  s~ > u  d~
      ms(dq,da,cq,sa)=WZbbmsq(8,7,3,4,5,6,2,1)
      ms(sq,sa,uq,da)=ms(dq,da,cq,sa)
      ms(bq,ba,uq,da)=ms(dq,da,cq,sa)
      ms(bq,ba,cq,sa)=ms(dq,da,cq,sa)
c--- d  u~ > u  u~ and s  c~ > c  c~
      ms(dq,ua,uq,ua)=WZuuidmsq(1,2,3,4,5,6,7,8)
      ms(sq,ac,cq,ac)=ms(dq,ua,uq,ua)
c--- d  u~ > c  c~ and s  c~ > u  u~ and s  c~ > b  b~ and d  u~ > b  b~
      ms(dq,ua,cq,ac)=WZccmsq(1,2,3,4,5,6,7,8)
      ms(sq,ac,uq,ua)=ms(dq,ua,cq,ac)
c--- d  u~ > d  d~ and s  c~ > s  s~
      ms(dq,ua,dq,da)=WZddidmsq(1,2,3,4,5,6,7,8)
      ms(sq,ac,sq,sa)=ms(dq,ua,dq,da)
c--- d  u~ > s  s~
      ms(dq,ua,sq,sa)=WZbbmsq(1,2,3,4,5,6,7,8)
c--- s  c~ > d  d~
      ms(sq,ac,dq,da)=ms(dq,ua,sq,sa)
      ms(sq,ac,bq,ba)=ms(dq,ua,sq,sa)
      ms(dq,ua,bq,ba)=ms(dq,ua,sq,sa)
c--- d  s~ > u  s~
      ms(dq,sa,uq,sa)=WZbbmsq(1,7,3,4,5,6,2,8)
c--- d  b~ > u  b~
      ms(dq,ba,uq,ba)=ms(dq,sa,uq,sa)
c--- s  d~ > c  d~
      ms(sq,da,cq,da)=ms(dq,sa,uq,sa)
c--- s  b~ > c  b~
      ms(sq,ba,cq,ba)=ms(dq,sa,uq,sa)
c--- d  c~ > d  s~
      ms(dq,ac,dq,sa)=WZbbmsq(8,2,3,4,5,6,7,1)
c--- s  u~ > s  d~
      ms(sq,ua,sq,da)=ms(dq,ac,dq,sa)
c--- b  u~ > b  d~
      ms(bq,ua,bq,da)=ms(dq,ac,dq,sa)
c--- b  c~ > b  s~
      ms(bq,ac,bq,sa)=ms(dq,ac,dq,sa)
c--- d  c~ > u  c~
      ms(dq,ac,uq,ac)=WZccmsq(1,7,3,4,5,6,2,8)
c--- s  u~ > c  u~
      ms(sq,ua,cq,ua)=ms(dq,ac,uq,ac)
c--- c  c~ > u  d~ and u  u~ > c  s~ and b  b~ > u  d~ and b  b~ > c  s~
      ms(uq,ua,cq,sa)=WZccmsq(8,7,3,4,5,6,2,1)
      ms(cq,ac,uq,da)=ms(uq,ua,cq,sa)
c--- u  c~ > u  s~
      ms(uq,ac,uq,sa)=WZccmsq(8,2,3,4,5,6,7,1)
c--- c  u~ > c  d~
      ms(cq,ua,cq,da)=ms(uq,ac,uq,sa)

C***********************2q-2g-processes*********************************
c--- d   u~ > g  g
      ms(dq,ua,gg,gg)=aveqq*half*WZggmsq(1,2,3,4,5,6,7,8)
c--- s   c~ > g  g
      ms(sq,ac,gg,gg)=ms(dq,ua,gg,gg)
c--- u~  d  > g  g
      ms(ua,dq,gg,gg)=aveqq*half*WZggmsq(2,1,3,4,5,6,7,8)
c--- c~  s  > g  g
      ms(ac,sq,gg,gg)=ms(ua,dq,gg,gg)

c--- g  g  > d~  u
      ms(gg,gg,da,uq)=avegg*WZggmsq(7,8,3,4,5,6,2,1)
c--- g  g  > s~  c
      ms(gg,gg,sa,cq)=ms(gg,gg,da,uq)

c--- d  g  > u g
      ms(dq,gg,uq,gg)=aveqg*WZggmsq(1,7,3,4,5,6,2,8)
c--- s  g  > c g
      ms(sq,gg,cq,gg)=ms(dq,gg,uq,gg)
c--- g  d > u g
      ms(gg,dq,uq,gg)=aveqg*WZggmsq(2,7,3,4,5,6,1,8)
c--- g  s > c g
      ms(gg,sq,cq,gg)=ms(gg,dq,uq,gg)

c--- g  u~  > d~  g
      ms(gg,ua,da,gg)=aveqg*WZggmsq(7,2,3,4,5,6,1,8)
c--- g  c~  > s~  g
      ms(gg,ac,sa,gg)=ms(gg,ua,da,gg)
c--- u~ g  > d~  g
      ms(ua,gg,da,gg)=aveqg*WZggmsq(7,1,3,4,5,6,2,8)
c--- c~ g > s~  g
      ms(ac,gg,sa,gg)=ms(ua,gg,da,gg)

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

c      write(6,*) 'qq'
c      include 'qq-'
c      write(6,*) 'aa'
c      include 'aa-'
c      write(6,*) 'qa'
c      include 'qa-'
c      write(6,*) 'aq'
c      include 'aq-'
c      call madcheckm(ms)

      endif

      msq(:,:)=0d0
      do i=-nf,nf
      do j=-nf,nf
      do k=-nf,nf
      do l=-nf,nf
      msq(i,j)=msq(i,j)+ms(i,j,k,l)
      enddo
      enddo
      enddo
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      DOUBLE PRECISION FUNCTION D1MACH(I)

C  DOUBLE-PRECISION MACHINE CONSTANTS (SEE R1MACH FOR DOCUMENTATION)

C  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST
C  TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE.

      INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4), SC
      DOUBLE PRECISION DMACH(5), EPS, EPSNEW, S

      EQUIVALENCE (DMACH(1),SMALL(1)), (DMACH(2),LARGE(1)),
     $            (DMACH(3),RIGHT(1)), (DMACH(4),DIVER(1)),
     $            (DMACH(5),LOG10(1))

      LOGICAL  PASS1
      SAVE     PASS1
      DATA     PASS1/.TRUE./

C IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T 3B SERIES AND
C MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T PC 7300),
C IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST.
cj HP, SUN and SIG are big-endian

      DATA (SMALL(N),N=1,2)/1048576,0/, (LARGE(N),N=1,2)/2146435071,-1/,
     $  (RIGHT(N),N=1,2)/1017118720,0/, (DIVER(N),N=1,2)/1018167296,0/,
     $  (LOG10(N),N=1,2)/1070810131,1352628735/, SC/987/

C IEEE ARITHMETIC MACHINES AND 8087-BASED MICROS, SUCH AS THE IBM PC
C AND AT&T 6300, IN WHICH THE LEAST SIGNIFICANT BYTE IS STORED FIRST.
cj ie. little-endian

cj      DATA (SMALL(N),N=1,2)/0,1048576/, (LARGE(N),N=1,2)/-1,2146435071/,
cj     $  (RIGHT(N),N=1,2)/0,1017118720/, (DIVER(N),N=1,2)/0,1018167296/,
cj     $  (LOG10(N),N=1,2)/1352628735,1070810131/, SC/987/

C AMDAHL MACHINES.

C      DATA (SMALL(N),N=1,2)/1048576,0/, (LARGE(N),N=1,2)/2147483647,-1/,
C     $ (RIGHT(N),N=1,2)/856686592,0/, (DIVER(N),N=1,2)/ 873463808,0/,
C     $ (LOG10(N),N=1,2)/1091781651,1352628735/, SC/987/

C BURROUGHS 1700 SYSTEM.

C      DATA (SMALL(N),N=1,2)/ZC00800000,Z000000000/,
C     $ (LARGE(N),N=1,2)/ZDFFFFFFFF,ZFFFFFFFFF/,
C     $ (RIGHT(N),N=1,2)/ZCC5800000,Z000000000/,
C     $ (DIVER(N),N=1,2)/ZCC6800000,Z000000000/,
C     $ (LOG10(N),N=1,2)/ZD00E730E7,ZC77800DC0/, SC/987/

C BURROUGHS 5700 SYSTEM.

C      DATA (SMALL(N),N=1,2)/O1771000000000000,O0000000000000000/,
C     $  (LARGE(N),N=1,2)/O0777777777777777,O0007777777777777/,
C     $  (RIGHT(N),N=1,2)/O1461000000000000,O0000000000000000/,
C     $  (DIVER(N),N=1,2)/O1451000000000000,O0000000000000000/,
C     $  (LOG10(N),N=1,2)/O1157163034761674,O0006677466732724/, SC/987/

C BURROUGHS 6700/7700 SYSTEMS.

C      DATA (SMALL(N),N=1,2)/O1771000000000000,O7770000000000000/,
C     $  (LARGE(N),N=1,2)/O0777777777777777,O7777777777777777/,
C     $  (RIGHT(N),N=1,2)/O1461000000000000,O0000000000000000/,
C     $  (DIVER(N),N=1,2)/O1451000000000000,O0000000000000000/,
C     $  (LOG10(N),N=1,2)/O1157163034761674,O0006677466732724/, SC/987/

C FTN4 ON THE CDC 6000/7000 SERIES.

C      DATA
C     $  (SMALL(N),N=1,2)/00564000000000000000B,00000000000000000000B/,
C     $  (LARGE(N),N=1,2)/37757777777777777777B,37157777777777777774B/,
C     $  (RIGHT(N),N=1,2)/15624000000000000000B,00000000000000000000B/,
C     $  (DIVER(N),N=1,2)/15634000000000000000B,00000000000000000000B/,
C     $  (LOG10(N),N=1,2)/17164642023241175717B,16367571421742254654B/,
C     $  SC/987/

C FTN5 ON THE CDC 6000/7000 SERIES.

C      DATA
C     $(SMALL(N),N=1,2)/O"00564000000000000000",O"00000000000000000000"/,
C     $(LARGE(N),N=1,2)/O"37757777777777777777",O"37157777777777777774"/,
C     $(RIGHT(N),N=1,2)/O"15624000000000000000",O"00000000000000000000"/,
C     $(DIVER(N),N=1,2)/O"15634000000000000000",O"00000000000000000000"/,
C     $(LOG10(N),N=1,2)/O"17164642023241175717",O"16367571421742254654"/,
C     $ SC/987/

C CONVEX C-1

C      DATA (SMALL(N),N=1,2)/'00100000'X,'00000000'X/,
C     $  (LARGE(N),N=1,2)/'7FFFFFFF'X,'FFFFFFFF'X/,
C     $  (RIGHT(N),N=1,2)/'3CC00000'X,'00000000'X/,
C     $  (DIVER(N),N=1,2)/'3CD00000'X,'00000000'X/,
C     $  (LOG10(N),N=1,2)/'3FF34413'X,'509F79FF'X/, SC/987/

C CRAY 1, XMP, 2, AND 3.

C      DATA
C     $ (SMALL(N),N=1,2)/201354000000000000000B,000000000000000000000B/,
C     $ (LARGE(N),N=1,2)/577767777777777777777B,000007777777777777776B/,
C     $ (RIGHT(N),N=1,2)/376434000000000000000B,000000000000000000000B/,
C     $ (DIVER(N),N=1,2)/376444000000000000000B,000000000000000000000B/,
C     $ (LOG10(N),N=1,2)/377774642023241175717B,000007571421742254654B/,
C     $ SC/987/

C DATA GENERAL ECLIPSE S/200
C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE -
C STATIC DMACH(5)

C      DATA SMALL/20K,3*0/, LARGE/77777K,3*177777K/,
C     $  RIGHT/31420K,3*0/, DIVER/32020K,3*0/,
C     $  LOG10/40423K,42023K,50237K,74776K/, SC/987/

C HARRIS SLASH 6 AND SLASH 7

C      DATA (SMALL(N),N=1,2)/'20000000,'00000201/,
C     $  (LARGE(N),N=1,2)/'37777777,'37777577/,
C     $  (RIGHT(N),N=1,2)/'20000000,'00000333/,
C     $  (DIVER(N),N=1,2)/'20000000,'00000334/,
C     $  (LOG10(N),N=1,2)/'23210115,'10237777/, SC/987/

C HONEYWELL DPS 8/70 SERIES.

C      DATA (SMALL(N),N=1,2)/O402400000000,O000000000000/,
C     $  (LARGE(N),N=1,2)/O376777777777,O777777777777/,
C     $  (RIGHT(N),N=1,2)/O604400000000,O000000000000/,
C     $  (DIVER(N),N=1,2)/O606400000000,O000000000000/,
C     $  (LOG10(N),N=1,2)/O776464202324,O117571775714/, SC/987/

C IBM 360/370 SERIES, XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86.

C      DATA (SMALL(N),N=1,2)/Z00100000,Z00000000/,
C     $  (LARGE(N),N=1,2)/Z7FFFFFFF,ZFFFFFFFF/,
C     $  (RIGHT(N),N=1,2)/Z33100000,Z00000000/,
C     $  (DIVER(N),N=1,2)/Z34100000,Z00000000/,
C     $  (LOG10(N),N=1,2)/Z41134413,Z509F79FF/, SC/987/

C INTERDATA 8/32 WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.

C      DATA (SMALL(N),N=1,2)/Z'00100000',Z'00000000'/,
C     $  (LARGE(N),N=1,2)/Z'7EFFFFFF',Z'FFFFFFFF'/,
C     $  (RIGHT(N),N=1,2)/Z'33100000',Z'00000000'/,
C     $  (DIVER(N),N=1,2)/Z'34100000',Z'00000000'/,
C     $  (LOG10(N),N=1,2)/Z'41134413',Z'509F79FF'/, SC/987/

C PDP-10 (KA PROCESSOR).

C      DATA (SMALL(N),N=1,2)/"033400000000,"000000000000/,
C     $  (LARGE(N),N=1,2)/"377777777777,"344777777777/,
C     $  (RIGHT(N),N=1,2)/"113400000000,"000000000000/,
C     $  (DIVER(N),N=1,2)/"114400000000,"000000000000/,
C     $  (LOG10(N),N=1,2)/"177464202324,"144117571776/, SC/987/

C PDP-10 (KI PROCESSOR).

C      DATA (SMALL(N),N=1,2)/"000400000000,"000000000000/,
C     $  (LARGE(N),N=1,2)/"377777777777,"377777777777/,
C     $  (RIGHT(N),N=1,2)/"103400000000,"000000000000/,
C     $  (DIVER(N),N=1,2)/"104400000000,"000000000000/,
C     $  (LOG10(N),N=1,2)/"177464202324,"047674776746/, SC/987/

C PDP-11 FORTRANS SUPPORTING 32-BIT INTEGERS
C (EXPRESSED IN INTEGER AND OCTAL).

C      DATA (SMALL(N),N=1,2)/8388608,0/, (LARGE(N),N=1,2)/2147483647,-1/,
C     $  (RIGHT(N),N=1,2)/612368384,0/, (DIVER(N),N=1,2)/620756992,0/,
C     $  (LOG10(N),N=1,2)/1067065498,-2063872008/, SC/987/

C      DATA (SMALL(N),N=1,2)/O00040000000,O00000000000/,
C     $  (LARGE(N),N=1,2)/O17777777777,O37777777777/,
C     $  (RIGHT(N),N=1,2)/O04440000000,O00000000000/,
C     $  (DIVER(N),N=1,2)/O04500000000,O00000000000/,
C     $  (LOG10(N),N=1,2)/O07746420232,O20476747770/, SC/987/

C PDP-11 FORTRANS SUPPORTING 16-BIT INTEGERS
C (EXPRESSED IN INTEGER AND OCTAL).

C      DATA SMALL/128,3*0/, LARGE/32767,3*-1/, RIGHT/9344,3*0/,
C     $  DIVER/9472,3*0/, LOG10/16282,8346,-31493,-12296/, SC/987/

C      DATA SMALL/O000200,3*O000000/, LARGE/O077777,3*O177777/,
C     $  RIGHT/O022200,3*O000000/, DIVER/O022400,3*O000000/,
C     $  LOG10/O037632,O020232,O102373,O147770/, SC/987/

C PRIME 50 SERIES SYSTEMS WITH 32-BIT INTEGERS AND 64V MODE
C INSTRUCTIONS, SUPPLIED BY IGOR BRAY.

C      DATA (SMALL(N),N=1,2)/:10000000000,:00000100001/,
C     $  (LARGE(N),N=1,2)/:17777777777,:37777677775/,
C     $  (RIGHT(N),N=1,2)/:10000000000,:00000000122/,
C     $  (DIVER(N),N=1,2)/:10000000000,:00000000123/,
C     $  (LOG10(N),N=1,2)/:11504046501,:07674600177/, SC/987/

C SEQUENT BALANCE 8000

C      DATA (SMALL(N),N=1,2)/$00000000, $00100000/,
C     $  (LARGE(N),N=1,2)/$FFFFFFFF, $7FEFFFFF/,
C     $  (RIGHT(N),N=1,2)/$00000000, $3CA00000/,
C     $  (DIVER(N),N=1,2)/$00000000, $3CB00000/,
C     $  (LOG10(N),N=1,2)/$509F79FF, $3FD34413/, SC/987/

C UNIVAC 1100 SERIES.

C      DATA (SMALL(N),N=1,2)/O000040000000,O000000000000/,
C     $  (LARGE(N),N=1,2)/O377777777777,O777777777777/,
C     $  (RIGHT(N),N=1,2)/O170540000000,O000000000000/,
C     $  (DIVER(N),N=1,2)/O170640000000,O000000000000/,
C     $  (LOG10(N),N=1,2)/O177746420232,O411757177572/, SC/987/

C VAX UNIX F77 COMPILER

C      DATA (SMALL(N),N=1,2)/128,0/, (LARGE(N),N=1,2)/-32769,-1/,
C     $  (RIGHT(N),N=1,2)/9344,0/, (DIVER(N),N=1,2)/9472,0/,
C     $  (LOG10(N),N=1,2)/546979738,-805796613/, SC/987/

C VAX-11 WITH FORTRAN IV-PLUS COMPILER

C      DATA (SMALL(N),N=1,2)/Z00000080,Z00000000/,
C     $  (LARGE(N),N=1,2)/ZFFFF7FFF,ZFFFFFFFF/,
C     $  (RIGHT(N),N=1,2)/Z00002480,Z00000000/,
C     $  (DIVER(N),N=1,2)/Z00002500,Z00000000/,
C     $  (LOG10(N),N=1,2)/Z209A3F9A,ZCFF884FB/, SC/987/

C VAX/VMS VERSION 2.2

C      DATA (SMALL(N),N=1,2)/'80'X,'0'X/,
C     $  (LARGE(N),N=1,2)/'FFFF7FFF'X,'FFFFFFFF'X/,
C     $  (RIGHT(N),N=1,2)/'2480'X,'0'X/, (DIVER(N),N=1,2)/'2500'X,'0'X/,
C     $  (LOG10(N),N=1,2)/'209A3F9A'X,'CFF884FB'X/, SC/987/

      IF( PASS1 )  THEN

         PASS1 = .FALSE.
         IF (SC.NE.987)
     $       CALL ERRMSG( 'D1MACH--NO DATA STATEMENTS ACTIVE',.TRUE.)

C                        ** CALCULATE MACHINE PRECISION
         EPSNEW = 0.01D0
   10    EPS = EPSNEW
            EPSNEW = EPSNEW / 1.1D0
C                               ** IMPORTANT TO STORE 'S' SINCE MAY BE
C                               ** KEPT IN HIGHER PRECISION IN REGISTERS
            S = 1.D0 + EPSNEW
            IF( S.GT.1.D0 ) GO TO 10
         IF( EPS/DMACH(4).LT.0.5D0 .OR. EPS/DMACH(4).GT.2.D0 )
     $       CALL ERRMSG( 'D1MACH--TABULATED PRECISION WRONG',.TRUE.)

      END IF

      IF (I.LT.1.OR.I.GT.5)
     $    CALL ERRMSG( 'D1MACH--ARGUMENT OUT OF BOUNDS',.TRUE.)
      D1MACH = DMACH(I)
      RETURN
      END


text:
module probdata_module

!     These determine the refinement criteria
      double precision, save ::   specerr, specgrad
      integer         , save ::  max_specerr_lev   ,max_specgrad_lev

!     Diffusion coefficient
      double precision, save :: diff_coeff

!     These help specify which specific problem
      integer         , save :: probtype

end module probdata_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
! +------------------+
! | Precision        |
! +------------------+
!
module set_precision

    implicit none
    integer, parameter :: i18=selected_int_kind(18);
    integer, parameter :: sp =selected_real_kind(6, 37)
    integer, parameter :: dp =selected_real_kind(15, 307)
end module set_precision

! +------------------+
! |  Timing          |
! +------------------+
module time

    implicit none
    integer, save :: t_m, t_n, t_o
    integer, save :: val0(8), val1(8), val2(8)
    real, save    :: etime1(2), etime2(2), tresult
    real, save    :: t0, t1, t2, t3 ! t4, tcpu1, tcpu2, twall1, twall2

    real, save    :: tloop_beg_dat, tloop_beg_mul, tloop_beg_snd, tloop_beg_cpy, &
                     tloop_beg_alp, tloop_beg_bet, tloop_beg_swp, tloop_beg_mv

    real, save    :: tloop_end_dat, tloop_end_mul, tloop_end_snd, tloop_end_cpy, &
                     tloop_end_alp, tloop_end_bet, tloop_end_swp, tloop_end_mv

    real, save    :: tloop_dat, tloop_mul, tloop_snd, tloop_cpy, &
                     tloop_cpu, tloop_alp, tloop_bet, tloop_swp, &
                     tloop_mv
    logical(kind=1), parameter :: print_more=.false.
  contains

  real function time_dt()
    integer :: val(8)
    call date_and_time(VALUES=val)
    time_dt=float(val(2)+val(3)*3600*24+val(5)*3600+val(6)*60+val(7))&
       + val(8)*.001
  end function

end module time

! +------------------+
! | MPI declarations |
! +------------------+
  module mpi_stuff

    use omp_lib

! ATTN: Using module conflicts in F90 bindings
    !use mpi
    implicit none
    INCLUDE 'mpif.h'

! Mpi variables
    integer mpierr, rank, procs, tag, sendreq, recvreq, req(2)
    integer sendst, recvst
    integer status(MPI_STATUS_SIZE)
! Openmp variables
    !integer OMP_GET_THREAD_NUM, OMP_GET_NUM_THREADS
    integer tid !, OMP_GET_THREAD_NUM
    integer ntd !, OMP_GET_NUM_PROCS, OMP_get_num_threads
    integer iter

  end module mpi_stuff

max_stars_repo_name: bthakur/benchmarks

max_issues_repo_name: bthakur/benchmarks

max_forks_repo_name: bthakur/benchmarks


text:
! simple code for random number sequences
! from numerical recipes..
function ranfb(ix)
  real :: y,ranfb
  integer :: ix,ix1,ix2,ia1,ia2,ic,j1,j2,j3

  ia1=273
  ia2=25301
  ic=226908345

  if (ix.lt.0) ix=-ix
  if (ix.gt.2**30) ix=MOD(ix,2**30)

  ix1=ix/(2**15)
  ix2=ix-ix1*(2**15)
  j1=MOD(ix1*ia2,2**15)*2**15
  j2=MOD(ix2*ia1,2**15)*2**15
  j3=ix2*ia2
  ix=MOD(j1+j2,2**30)
  ix=MOD(ix+j3,2**30)
  ix=MOD(ix+ic,2**30)
  y=FLOAT(ix)

  ranfb=y/2**30

end function

function ran2(ix)
 real ran2,ranfb
 integer ix
 ran2=ranfb(ix)
end function

FUNCTION gasdev(iseed)
 REAL :: gasdev,x(2),r2,fac,ran2
 INTEGER :: iseed
 INTEGER, SAVE :: iset=0
 REAL, SAVE :: saved
 IF(iset.EQ.0) THEN
   r2=2.
   DO WHILE(r2.GE.1.OR.r2.EQ.0)
     x(1)=ran2(iseed)
     x(2)=ran2(iseed)
     x=2*x-1
     r2=sum(x**2)
   ENDDO
   fac=sqrt(-2.*log(r2)/r2)
   saved=x(1)*fac
   gasdev=x(2)*fac
   iset=1
 ELSE
   iset=0
   gasdev=saved
 ENDIF
END FUNCTION

!
   SUBROUTINE makefractal(nstar,rstar,vstar,fdim,idum)
!
! makes fractal of nstar particles of dimension fdim, using ndiv
! subunits forcing the number of cells if force=.true.
! Simon Goodwin & Ant Whitworth (2004, A&A, 413, 929)
!
! fdim (real) = fractal dimension
! ndiv (int) = number of divisions of box (>=2)
! nstar (int) = number of stars in *final* distribution
! rstar (3,n) = positions to return
! vstar (3,n) = velocities to return
! force (logical) = make box dimension such that integer no. children
! vtype (logical) = .true. - inherit velocities
! vir (real) = virial ratio (virialised = 0.5)
! idum = random number seed
! outfile = output file name
   IMPLICIT NONE
   INTEGER :: ndiv,nstar,idum,isnext,iforce
   INTEGER :: is,nsmax,islast,flagabort,k
   INTEGER :: k1,k2,k3,i,isfirst,nswap,ndo,nsnow
   REAL :: fdim,dum,delta,dr(3),prob,rad,cov(3)
   REAL :: how,how2,ran2,com(3)
   REAL :: gasdev
   LOGICAL :: force,vtype
   INTEGER :: n,isurv,level
   REAL :: rstar(1:3,nstar),vstar(1:3,nstar)
   INTEGER, DIMENSION(:), ALLOCATABLE :: forced
   REAL, DIMENSION(:,:), ALLOCATABLE :: r,v
!
! set initial parameters
   ndiv=2
   force=.TRUE.
   vtype=.TRUE.
   ALLOCATE(r(1:3,1:nstar*50))
   ALLOCATE(v(1:3,1:nstar*50))
   ALLOCATE(forced(1:ndiv**3))
!
! ===================================================================
! probability of a child surviving
   prob=REAL(ndiv)**(fdim-3.)
!
! list characteristics
   WRITE(6,*) 'fractal dimension', fdim
   WRITE(6,*) 'ndiv', ndiv
   WRITE(6,*) 'number of stars', nstar
   WRITE(6,*) 'survival chance', prob
   WRITE(6,*) 'forcing?', force
!
! warning
   how=REAL(ndiv)**fdim
   WRITE(6,*) '# of children expected', how
   how2=how - INT(how)
   IF (how2>0.2 .AND. how2<0.8) THEN
     WRITE(6,*) 'WARNING:'
     IF (how2<0.5) THEN
       dum=LOG(REAL(INT(how)))/LOG(REAL(ndiv))
       WRITE(6,*) 'could look like fdim=', dum
     ELSE
       dum=LOG(REAL(INT(how)+1))/LOG(REAL(ndiv))
       WRITE(6,*) 'could look like fdim=', dum
     END IF
   END IF
!
   IF (force) THEN
     iforce=INT(how + 0.5)
     WRITE(6,*) 'forcing to survivors=', iforce
     dum=LOG(REAL(iforce))/LOG(REAL(ndiv))
     WRITE(6,*) 'which means fdim=', dum
   END IF
!
! build fractal cube
   r=0.
   v=0.
   nsmax=1
   islast=1
   flagabort=0
   delta=1./REAL(ndiv)
   level=0
   DO
     level=level + 1
     isfirst=nsmax + 1
     isnext=isfirst
! loop over all parents
     DO is=1,nsmax
!
! if forced decide on survivors now
       IF (force) THEN
         forced=0
         DO n=1,iforce
4          isurv=INT(ran2(idum)*REAL(ndiv**3)) + 1
           IF (forced(isurv)==1) GOTO 4
           forced(isurv)=1
         END DO
       END IF
! loop over all ndiv**3 possible children
       n=0
       DO k1=1,ndiv
         dr(1)=REAL(2*k1-1-ndiv)*delta
         DO k2=1,ndiv
           dr(2)=REAL(2*k2-1-ndiv)*delta
           DO k3=1,ndiv
             n=n+1
             dr(3)=REAL(2*k3-1-ndiv)*delta
             r(1:3,isnext)=r(1:3,is) + dr(1:3)
             IF (vtype) THEN
! keep parent's velocity
               v(1:3,isnext)=v(1:3,is)
             END IF
! do children survive?
! forced or not forced here
             IF (force) THEN
               IF (forced(n)==1) THEN
                 islast=isnext
                 isnext=isnext + 1
               END IF
             ELSE
               IF (ran2(idum)<prob) THEN
                 islast=isnext
                 isnext=isnext + 1
               END IF
             END IF
           END DO
         END DO
       END DO
     END DO
! are there children (we hope so)
     IF (islast/=nsmax) THEN
! delete parents
       DO is=isfirst,islast
         r(1:3,is-nsmax)=r(1:3,is)
         v(1:3,is-nsmax)=v(1:3,is)
       END DO
! new total number
       nsmax=islast - nsmax
       DO is=nsmax+1,islast
         r(1:3,is)=0.
         v(1:3,is)=0.
       END DO
! resize delta
       delta=delta/REAL(ndiv)
     ELSE
       flagabort=1
       EXIT
     END IF
     IF (flagabort==1) STOP 'no children'
! add noise to positions and velocity
     DO is=1,nsmax
       DO k=1,3
         r(k,is)=r(k,is) + (2.*ran2(idum) - 1.)*delta
!         v(k,is)=v(k,is) + 2.**(fdim-1.)*gasdev(idum)
         v(k,is)=v(k,is) + delta*gasdev(idum)
       END DO
     END DO
! exit if enough particles
     IF (nsmax>3.5*nstar) EXIT
!     WRITE(6,*) 'done', delta, nsmax
   END DO
   WRITE(6,*) 'Made initial cubic distribution'
!
! cut out sphere
   nsnow=nsmax
   DO is=1,nsmax
     rad=r(1,is)**2 + r(2,is)**2 + r(3,is)**2
     IF (rad>1.) THEN
       r(1,is)=69.
       nsnow=nsnow - 1
     END IF
   END DO
   nswap=nsmax
   ndo=1
   DO
8    IF (r(1,ndo)==69.) THEN
       r(1:3,ndo)=r(1:3,nswap)
       nswap=nswap - 1
       GOTO 8
     END IF
     ndo=ndo + 1
     IF (ndo==nsnow) EXIT
   END DO
   nsmax=ndo - 1
   IF (nsmax<nstar) THEN
     WRITE(6,*) 'nsmax<nstar - change random number seed'
     nstar=nsmax
  ENDIF
!
! we now have a fractal of nsmax stars
! need to randomly remove stars until we get nstar
   DO
! star to remove
     IF (nsmax==nstar) EXIT
     is=INT(ran2(idum)*REAL(nsmax)) + 1
     r(1:3,is)=r(1:3,nsmax)
     v(1:3,is)=v(1:3,nsmax)
     nsmax=nsmax - 1
   END DO
   WRITE(6,*) 'Cut to size', nsmax
!
! now put r in rstar array
   DO is=1,nstar
     rstar(1:3,is)=r(1:3,is)
     vstar(1:3,is)=v(1:3,is)
   END DO
!
   IF (.NOT.vtype) THEN
     DO is=1,nstar
       DO k=1,3
         vstar(k,is)=(2.*ran2(idum) - 1.)
       END DO
     END DO
   END IF
!
! set coms to zero
   com=0.
   cov=0.
   DO i=1,nstar
     com(1:3)=com(1:3) + rstar(1:3,i)
     cov(1:3)=cov(1:3) + vstar(1:3,i)
   END DO
   com=com/REAL(nstar)
   cov=cov/REAL(nstar)
   WRITE(6,*) 'COM', com(1), com(2), com(3)
   DO i=1,nstar
     rstar(1:3,i)=rstar(1:3,i) - com(1:3)
     vstar(1:3,i)=vstar(1:3,i) - cov(1:3)
   END DO
!
   END SUBROUTINE makefractal
!

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
module qlknn_all_nets
    use net_efeetg_gb
    use net_efeitg_gb_div_efiitg_gb
    use net_efetem_gb
    use net_efiitg_gb
    use net_efitem_gb_div_efetem_gb
    use net_pfeitg_gb_div_efiitg_gb
    use net_pfetem_gb_div_efetem_gb
    use net_dfeitg_gb_div_efiitg_gb
    use net_dfetem_gb_div_efetem_gb
    use net_vteitg_gb_div_efiitg_gb
    use net_vtetem_gb_div_efetem_gb
    use net_vceitg_gb_div_efiitg_gb
    use net_vcetem_gb_div_efetem_gb
    use net_dfiitg_gb_div_efiitg_gb
    use net_dfitem_gb_div_efetem_gb
    use net_vtiitg_gb_div_efiitg_gb
    use net_vtitem_gb_div_efetem_gb
    use net_vciitg_gb_div_efiitg_gb
    use net_vcitem_gb_div_efetem_gb
    use net_gam_leq_gb
    use net_efeitg_gb_div_efeitg_gb_rot0
    use net_efetem_gb_div_efetem_gb_rot0
    use net_efiitg_gb_div_efiitg_gb_rot0
    use net_efitem_gb_div_efitem_gb_rot0
    use net_pfeitg_gb_div_pfeitg_gb_rot0
    use net_pfetem_gb_div_pfetem_gb_rot0
    use net_dfeitg_gb_div_dfeitg_gb_rot0
    use net_dfetem_gb_div_dfetem_gb_rot0
    use net_vteitg_gb_div_vteitg_gb_rot0
    use net_vtetem_gb_div_vtetem_gb_rot0
    use net_vceitg_gb_div_vceitg_gb_rot0
    use net_vcetem_gb_div_vcetem_gb_rot0
    use net_dfiitg_gb_div_dfiitg_gb_rot0
    use net_dfitem_gb_div_dfitem_gb_rot0
    use net_vtiitg_gb_div_vtiitg_gb_rot0
    use net_vtitem_gb_div_vtitem_gb_rot0
    use net_vciitg_gb_div_vciitg_gb_rot0
    use net_vcitem_gb_div_vcitem_gb_rot0

    implicit none

    type(networktype), dimension(20), save :: nets
    type(networktype), dimension(19), save :: rotdiv_nets
contains
    subroutine load_nets()
        logical, save :: first_call = .TRUE.
        if (first_call) then
            nets( 1) = efeetg_gb()
            nets( 2) = efeitg_gb_div_efiitg_gb()
            nets( 3) = efetem_gb()
            nets( 4) = efiitg_gb()
            nets( 5) = efitem_gb_div_efetem_gb()
            nets( 6) = pfeitg_gb_div_efiitg_gb()
            nets( 7) = pfetem_gb_div_efetem_gb()
            nets( 8) = dfeitg_gb_div_efiitg_gb()
            nets( 9) = dfetem_gb_div_efetem_gb()
            nets(10) = vteitg_gb_div_efiitg_gb()
            nets(11) = vtetem_gb_div_efetem_gb()
            nets(12) = vceitg_gb_div_efiitg_gb()
            nets(13) = vcetem_gb_div_efetem_gb()
            nets(14) = dfiitg_gb_div_efiitg_gb()
            nets(15) = dfitem_gb_div_efetem_gb()
            nets(16) = vtiitg_gb_div_efiitg_gb()
            nets(17) = vtitem_gb_div_efetem_gb()
            nets(18) = vciitg_gb_div_efiitg_gb()
            nets(19) = vcitem_gb_div_efetem_gb()
            nets(20) = gam_leq_gb()

            rotdiv_nets(2) = efeitg_gb_div_efeitg_gb_rot0()
            rotdiv_nets(3) = efetem_gb_div_efetem_gb_rot0()
            rotdiv_nets(4) = efiitg_gb_div_efiitg_gb_rot0()
            rotdiv_nets(5) = efitem_gb_div_efitem_gb_rot0()
            rotdiv_nets(6) = pfeitg_gb_div_pfeitg_gb_rot0()
            rotdiv_nets(7) = pfetem_gb_div_pfetem_gb_rot0()
            rotdiv_nets(8) = dfeitg_gb_div_dfeitg_gb_rot0()
            rotdiv_nets(9) = dfetem_gb_div_dfetem_gb_rot0()
            rotdiv_nets(10) = vteitg_gb_div_vteitg_gb_rot0()
            rotdiv_nets(11) = vtetem_gb_div_vtetem_gb_rot0()
            rotdiv_nets(12) = vceitg_gb_div_vceitg_gb_rot0()
            rotdiv_nets(13) = vcetem_gb_div_vcetem_gb_rot0()
            rotdiv_nets(14) = dfiitg_gb_div_dfiitg_gb_rot0()
            rotdiv_nets(15) = dfitem_gb_div_dfitem_gb_rot0()
            rotdiv_nets(16) = vtiitg_gb_div_vtiitg_gb_rot0()
            rotdiv_nets(17) = vtitem_gb_div_vtitem_gb_rot0()
            rotdiv_nets(18) = vciitg_gb_div_vciitg_gb_rot0()
            rotdiv_nets(19) = vcitem_gb_div_vcitem_gb_rot0()
            first_call = .FALSE.
        end if
    end subroutine load_nets
end module qlknn_all_nets

max_stars_repo_name: QuaLiKiz-group/QLKNN-fortran

max_issues_repo_name: QuaLiKiz-group/QLKNN-fortran

max_forks_repo_name: QuaLiKiz-group/QLKNN-fortran